from scipy.misc import factorial,comb
factorial(10)/(factorial(6)*factorial(4)),comb(10,4,exact=True)

plot(range(1,71),1-cumprod(arange(365,365-70,-1)/365.),'r+')
plot(range(2,71),1-cumprod(ones(69)*364./365.),'g+');

figure(figsize=(8,6))
plot(range(1,71),1-cumprod(arange(365,295,-1)/365.),'r+')
plot(range(2,71),1-cumprod(ones(69)*364./365.),'g+')
yticks(arange(0,1.1,.1));
xlim(0,70)
xticks(range(0,71,5))
grid('on')
legend([r'$1-\frac{365!}{(365-n)!\ 365^n}$',
        r'$1-(\frac{364}{365})^{n-1}$'],
       numpoints=3,loc='upper left',fontsize=18);

def bdayplot(n,inc=25,sym='',loc=None,):
  plot(range(1,n+1),1-cumprod(arange(365,365-n,-1)/365.),'r'+sym)
  plot(range(2,n+1),1-cumprod(ones(n-1)*364./365.),'g'+sym)
  xticks(range(0,n+1,inc))
  yticks(arange(0,1.1,.1))
  xlim(0,n)
  grid('on')
  legend([r'$1-\frac{365!}{(365-n)!\ 365^n}$',
        r'$1-(\frac{364}{365})^{n-1}$'],
       numpoints=3,loc=loc,fontsize=17);
    
figure(figsize=(16,6))
subplot(1,2,1)
bdayplot(70,5,'+','upper left')
subplot(1,2,2)
bdayplot(365)
plot([70,70],[0,1],c='gray',linestyle='--');

figure(figsize=(8,6))
plot(range(1,71),1-cumprod(arange(365,365-70,-1)/365.),'r+')
plot(range(2,71),1-cumprod(ones(69)*364./365.),'g+')
xlim(0,70)
ylim(0,1)
grid('on')
plot([n*(n-1)/(2.*365) for n in range(16)],'b')
plot([1-exp(-(n-1)/365.) for n in range(366)],'y')
xticks(range(0,75,5))
yticks(arange(0,1.1,.1));
legend([r'$1-\frac{365!}{(365-n)!\ 365^n}$',
        r'$1-(\frac{364}{365})^{n-1}$',
        r'$\frac{n(n-1)}{2}\frac{1}{365}$',
        r'$1-{\rm e}^{-(n-1)/365}$'],
        numpoints=3,loc='upper left',fontsize=14);