After using the code below to create a series of .png images, from the command line, use, e.g., ImageMagick convert to create an animated gif from the .pngs
convert -delay 20 -loop 0 img*.png bsim.gif
from scipy.misc import comb
r=exp(log(2)/10)
rest=map(int,32*cumprod(r*ones(60))+.5) #log-spaced steps from 32 to 2048
trials=10000
f=2 # for p=1/2
#f=6 #for p=1/6
p=1./f
sigma=p*(1-p)
for k,n in enumerate(range(1,33)+rest):
figure(figsize=(6,4))
hist([sum(np.random.rand(n) < p) for i in xrange(trials)],bins=range(0,n+2),
linewidth=0,align='left')
#simulation done entirely in the above line, rest just adds text annotations
if n<=32:
m=arange(n+1)
ms=trials*max(p**m*(1-p)**(n-m)*comb(n,m))
else:
ms = trials/sqrt(2*pi*sigma*n)
#ms is used to set the upper ylim, so that it decreases smoothly with n
text(0.7, 0.5,'n={}'.format(n), horizontalalignment='left',
verticalalignment='center',fontsize=24,
transform=gca().transAxes) #write the n=...
text(0.7, 0.3,'$\sigma/\mu={}$'.format(round(sqrt((1-p)/(p*n)),2)), horizontalalignment='left',
verticalalignment='center',fontsize=18,
transform=gca().transAxes) #write the sigma/mu=...
annotate('n/{}'.format(f),(p*n,0),(p*n,-ms/6.4),arrowprops=dict(arrowstyle="->"),ha='center',fontsize=20,color='b')
xlim(-.5,n+.5*(1-p)/p)
ylim(0,1.1*ms)
savefig('pngs/img{:03d}.png'.format(k)) #save 92 images in pngs subdir
close()