clear all;
close all;
%Model 0: C 
N=10^6;

repeat=1000;
n1=floor((10*2*N)^(0.5));  
for j=1:10;
    ec=10+10*(j-1);
    n=floor(sqrt(2*ec*N));
    ns(j)=n;
    for i=1:repeat        
        idxrn=floor(rand(n,1).*N);
        C=getCollisions(idxrn);
        hatNsrn(j,i)=n^2/2/C;
        Cs(i,j)=C;
    end
    hatNrn=hatNsrn(j,:);
    vn_observed(j)=var(hatNrn);
         
    EC(j)=n^2/2/N;
    EC_obs(j)=mean(Cs(:,j))
    VC_obs(j)=var(Cs(:,j))
    
end
rsen_obs=vn_observed.^(0.5)./N;

f=figure;


subplot(2,2,1);
j=1;
%bin=floor(0.3*EC(j)):1:floor(2.5*EC(j));
[h,b]=hist(Cs(:,j));

bino1=binopdf(floor(b),floor(ns(j)^2/2),1/N).*repeat.*(b(2)-b(1));  
plot(b,h, 'o-', b,bino1,'-r', 'LineWidth',2);
legend('Hist', 'Bino'); legend('boxoff');
xlabel('Collisions C');
ylabel('Frequency');
title('(A)');
box off;


j=10;
[h,bin]=hist(Cs(:,j));
bino=binopdf(floor(bin),floor(ns(j)^2/2),1/N).*repeat.*(bin(2)-bin(1));  
subplot(2,2,2);
plot(bin, h, 'o-',bin,bino,'-r', 'LineWidth',2)
legend('Hist', 'Bino'); legend('boxoff');
xlabel('Collisions C');
ylabel('Frequency');
title('(B)');
box off;

subplot(2,2,3);
plot( ns, VC_obs,'-o',ns, EC_obs,'-r','LineWidth',2);
xlabel('Sample size n');
ylabel('var of C');
legend('Obs.','Est.',2);
legend('boxoff');
title('(C)');
box off;

subplot(2,2,4);
plot( ns, rsen_obs,'-o',ns, 1./EC_obs.^(0.5),'-r','LineWidth',2);
xlabel('Sample size n');
ylabel('RSE of N');
legend('Obs.', 'Est.');
legend('boxoff');
title('(D)');
box off;

saveas(f,'c:\size\size_simC0.eps', 'psc2');


