Python3 code to show the convergence of a function to the Gaussian distribution using matplotlib animations

gaussian_distribution_approximation_graph

# 
# Animated graph of two functions whose difference tends to 0 as n -> inf
# Using matplotlib.animation
#
# author: Julien Harbulot
# article url: http://julienharbulot.com/data-science/bernoulli-urn/
# licence: MIT Licence
#

%matplotlib inline

import math
import numpy as np 
import scipy.special
from scipy.stats import norm
from scipy.special import gamma
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

n = 5.
f = .5
def real(x, f = f, n = n):
    r = n * f
    return (gamma(n+2) / (gamma(r+1) * gamma(n - r + 1))) * x**r * (1-x)**(n-r)

def nreal(x, f = f, n = n):
    return real(x, f, n) / real(f, f, n)

def approx(x, f = f, n = n):
    sigma2 = f * (1-f) / float(n)
    mean   = f
    return 1 / math.sqrt(2 * math.pi * sigma2) * math.exp(- (x - mean)**2 / (2 * sigma2))

def napprox(x, f = f, n = n):
    return approx(x, f, n) / real(f, f, n)

# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(0, 1.1))
line1, = ax.plot([], [], lw=2)
line2, = ax.plot([], [], lw=2)
text = ax.text(0.46, 0.5, '', transform=ax.transAxes)
site = ax.text(0.01, 0.95, 'julienharbulot.com', transform=ax.transAxes)


# initialization function: plot the background of each frame
def init():
    line1.set_data([], [])
    line2.set_data([], [])
    text.set_text('')
    return [line1, line2]

# animation function.  This is called sequentially
steps = 100
factor = 4
real_frames = factor * steps
start = 10

x = np.linspace(0, 1, 100)
def animate(i):
    if i <= start:
        i = 1.0
    else:
        i = float(i - start) / float(factor) + 1
    text.set_text("n = {}".format(math.floor(i)))
    line1.set_data(x, [nreal(xi, f, i) for xi in x])
    line2.set_data(x, [napprox(xi, f, i) for xi in x])
    return [line1, line2]

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=real_frames - steps, interval=25, blit=True)
anim.save('line.gif', dpi=80, writer='imagemagick')