%Coupons
%parameters

for (n=[5:15])
%n=60;

%kmax=floor(4*n*log(n));
%kmax=45;

kmax=100;

pr=zeros(kmax);
maxsummand=zeros(kmax);

for (k=[n:kmax])
    coeff=1;
    bin=1;
    for (i=[0:n-1])
        summand=bin*(1-i/n)^k;
        pr(k)=pr(k)+coeff*summand;
        if (summand>maxsummand(k)) 
            maxsummand(k)=summand;
        end;
            
        coeff=-coeff;
        bin=bin*(n-i)/(i+1);
    end;
end;
 
for (k=[n:kmax])
     display(sprintf('%d %1.8f %10.1f', k, pr(k), maxsummand(k)));
end

%clf reset;
xlim([5,kmax]);
ylim([0,1]);
hold on;
plot([5:kmax], pr(5:kmax), 'Color', [n/60,0,1-n/60]);
%plot([1200:kmax], pr(1200:kmax), 'Color', [1,0,1]);
 
 
str=strcat('Prob. that all birthdays are represented with ', sprintf(' 5<=n<=%d', n));
title(str);
xlabel('k');
ylabel('prob.');
hold on;
end

hold off;

 
 
