clear all;
close all;
%%%% vaiance of size estimation
%%%
f=figure;
for data=1:2
    N=10^4;
    mud=10;
    p=((1:N)+0)'.^(-0.8-0.10*data);
    
    zeta=sum(p);
    A=N*mud/zeta;
    d=p.*A;
    d=ceil(d);
    [1./mud, var(1./d)];
    true_mud=mean(d);
    mud_recp=sum(1./d)/N;
    gammasqr=mean(d.^2)./mean(d)^2-1;
    Gamma=gammasqr+1;
    
    for j=1:10;
        n=500*j;
        ns(j)=n;
        gammasqrs(j)=gammasqr;
        repeat=50;
        [idxs,dxs]=getPPSSample(d, repeat*n);
        for i=1:repeat
            dx=dxs((i-1)*n+1:i*n);
            idx=idxs((i-1)*n+1:i*n);
            C(j,i)=getCollisions(idx');
            hatN(j,i)=n^2/(2*C(j,i))*Gamma;
        end
        hatNs=hatN(j,:);
        Cs=C(j,:);
        rse_observed(j)=(var(hatNs)/N^2)^0.5;
        rse_est(j)=1./mean(Cs)^(0.5);  
        
    end
    %f1=figure;
    %plot(Ns,gammasqr,'-o', 'LineWidth',2);
    %xlabel('Data size N','fontSize',12);
    %ylabel('\gamma^2','fontSize',12);
    %saveas(f1,'c:\community\n_gamma_.eps', 'psc2');
   
    t={'(A)', '(B)', '(C)', '(D)', '(E)', '(F)'};
    
    subplot(2,3,(data-1)*3+1);
    %loglog(flipud(sortrows(d)),'o');
    bin=1:max(d);
    h=hist(d, bin);
    loglog(bin, h,'o');
    xlabel('degree');
    ylabel('Frequency');
    title(t((data-1)*3+1));
    box off;
    
    subplot(2,3,(data-1)*3+2);
    loglog(flipud(sortrows(d)),'o');
    xlabel('degree');
    ylabel('Frequency');
    title(t((data-1)*3+2));
    box off;
    
    subplot(2,3,(data-1)*3+3);
    plot(ns,rse_observed,'-o',ns,rse_est, '+-','LineWidth',2);
    xlabel('n');
    ylabel('RSE');
    legend('True value', 'Estimated');
    legend('boxoff')
    title(t((data-1)*3+3));
    box off;
    
end

 subplot(1,3,1);
    %loglog(flipud(sortrows(d)),'o');
    bin=1:max(d);
    h=hist(d, bin);
    loglog(bin, h,'o');
    xlabel('Degree');
    ylabel('Frequency');
    title(t((data-1)*3+1));
    box off;
f=figure
   
 subplot(1,2,1);
    loglog(flipud(sortrows(d)),'o');
    xlabel('Rank');
    ylabel('Degree');
    title('(A)');
    box off;
    
    rse_est2=(2.*N./(ns.^2.*Gamma).*(1+ns./(0.8.*mud))).^(0.5);
    subplot(1,2,2);
    plot(ns,rse_observed,'-o',ns,rse_est2, '+-',ns, rse0,'r-','LineWidth',2);
    xlabel('n');
    ylabel('RSE');
    legend('RE Observed', 'RE Estimated', 'RN');
    legend('boxoff')
    title('(B)');
    box off;
   
saveas(f,'c:\community\size_sim_ns.eps', 'psc2');

