[alpha,beta]=meshgrid(0:0.1:1,0:0.1:1);
% Alpha: weight given to merit over celebrity
% Beta: noise level

N=1000; % Total population
K=5; % Number of people to adore
N0=10; % Initial population

C=alpha*0; %Array for correlations

ntrials=10; % Number of trials to average

for i=1:size(alpha,1)
    for j=1:size(alpha,2)
        cccc=0; % store sum of correlations
        for trial=1:ntrials

            merit=1./rand(1,N); % real merit, a=2 power law
            connection=ceil(rand(N,K)*N0); % Initial admiration, beyond N0 will be overwritten

            for n=N0:N
                cc=connection(1:(n-1),:);
                [fame,dummy]=hist(cc(:),1:(n-1)); % Calc fame
                noisymerit=(1-beta(i,j))*merit+beta(i,j)./rand(1,N);  % Combine merit with noise into apparent merit
                score = alpha(i,j)*(noisymerit(1:(n-1))/sum(merit(1:(n-1)))) + (1-alpha(i,j))*fame/sum(fame);
                cumscore=cumsum(score); cumscore=cumscore/sum(score);
                for k=1:K
                    r=rand; % Select target with probaility proportional to score
                    c=1;
                    while(cumscore(c)<r)
                        c=c+1;
                    end
                    connection(n,k)=c;
                end
            end
            [fame,dummy]=hist(connection(:),1:N);
            ccc=corrcoef(log(1+merit),log(1+fame));
            cccc=cccc+ccc(1,2);
        end
        C(i,j)=cccc/ntrials;

        loglog(merit,fame,'.')
        xlabel('Merit');
        ylabel('Fame')
        drawnow
    end
end

surf(alpha,beta,C)
xlabel('Weight to merit')
ylabel('Randomness')
zlabel('Correlation')