yunyin/main.m

%1-C(365,n)n!/(365)^n
init=365;
k=0;
for n=4:init
    pk=1-nchoosek(365,n)*factorial(n)/power(365,n);
    %pk=1-factorial(364)/(power(365,n)*factorial(n-2));
    if pk>=0.5
     k=n;
     break;
    end
end
k