    % C does nt follow binomial distribution
    % its variance is larger than binomial dist EC, smaller than EC+n^3p3
    %N=10^6, beta=50 EC=44, VC_obs=90, VC_p3=100
    %N=10^4, EC=33, VC_obs=44, VC_p3=72
    %                38   59  88
clear all;
close all;

repeat=100;
for j=1:10;
    N=10^4*j;
    Ns(j)=N;
    mud=10;
    p=(50+(1:N))'.^(-1);
    zeta=sum(p);
    A=N*mud/zeta;
    d=p.*A;
    d=floor(d+0.5);
    true_mud=mean(d);
    mud_recp=sum(1./d)/N;
    gammasqr=mean(d.^2)./mean(d)^2-1
    gammasqrs(j)=gammasqr;
    G=gammasqr+1;
    n10rn=floor((10*2*N)^(0.5));  
      [idxs,dxs]=getPPSSample(d, repeat*n10rn);
    [n10re, hatN]=nForRSE(dxs',idxs',N,G,0.2,repeat);
    
    ns_rn(j)=n10rn;
    ns_re(j)=n10re;
    for i=1:repeat
        %        dx=dxs((i-1)*n+1:i*n);
        idx=idxs((i-1)*n10re+1:i*n10re);
        C=getCollisions(idx');
        hatNs(j,i)=n10re^2*G/2/C;
       
        Cs(i,j)=C;
    end
    hatN=hatNs(j,:);
    %hatNrn=hatNsrn(j,:);
   % vn_observed(j)=var(hatNrn);
    ve_observed(j)=var(hatN);
    
    n=n10re;
    tau=sum(d);
    pi=d./tau;
    p2i=pi.^2;
    p3i=pi.^3;
    p2=sum(p2i);
    p3=sum(p3i);
    EC=n^2/2*p2
    EC_obs(j)=mean(Cs(:,j))
    VC_obs(j)=var(Cs(:,j))
    VC_p3(j)=EC+n^3*p3
    
      
end

f=figure;

subplot(2,3,1);
loglog(flipud(sortrows(d)),'o');
xlabel('Rank');
ylabel('Frequency');
title('(A)');
box off;

subplot(2,3,2);
hist(Cs(1,:));
%plot(Ns,gammasqrs+1,'-o','LineWidth',2);
xlabel('N');
ylabel('\Gamma');
title('(B)');
box off;

 bin=1:5:2*EC
    bino=binopdf(1:5:2*EC,floor(n^2/2),p2);  
    h=hist(Cs(:,j), bin);
   
subplot(2,3,3);
   plot(bin,h./repeat./5,'o',bin,bino,'-r')
 xlabel('Est N');
ylabel('Frequency');
title('(C)');
box off;
ns=1:10;
subplot(2,3,3);
plot(ns, EC_obs, 'o-',ns,VC_obs,'x', ns,VC_p3, '-');
legend('EC','VC_obs', 'VC_p3');
aveas(f,'c:\community\size_simC.eps', 'psc2');


