This code is taken from https://matplotlib.org/examples/animation/index.html and extended for the normally distributed case.
%matplotlib inline
from scipy.stats import beta
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
class UpdateDistbeta(object):
def __init__(self, ax, prob=0.5):
self.success = 0
self.prob = prob
self.line, = ax.plot([], [], 'k-',color='blue')
self.data_point = ax.text(0.02, 0.95, '', transform=ax.transAxes)
self.success_rate = ax.text(0.02, 0.90, '', transform=ax.transAxes)
self.bayes_estimator = ax.text(0.02, 0.85, '', transform=ax.transAxes)
self.x = np.linspace(0, 1, 200)
self.ax = ax
# Set up plot parameters
self.ax.set_xlim(0, 1)
self.ax.set_ylim(0, 15)
self.ax.grid(True)
# This vertical line represents the theoretical value, to
# which the plotted distribution should converge.
self.ax.axvline(prob, linestyle='--', color='red')
def init(self):
self.success = 0
self.line.set_data([], [])
return self.line,
def __call__(self, i):
# This way the plot can continuously run and we just keep
# watching new realizations of the process
if i == 0:
return self.init()
# Choose success based on exceed a threshold with a uniform pick
data = 0
if np.random.rand(1,) < self.prob:
self.success += 1
data = 1
y = beta.pdf(self.x, self.success + 1, (i - self.success) + 1)
self.line.set_data(self.x, y)
rate = self.success/float(i)
bayes_estimator = (self.success+1)/(i+2)
self.data_point.set_text('Data point = %i' % data)
self.success_rate.set_text('MLE for success rate = %.5f' % rate)
self.bayes_estimator.set_text('Bayes estimator = %.5f' % bayes_estimator)
return self.line,
# Fixing random state for reproducibility
np.random.seed(19680801)
fig, ax = plt.subplots()
ud = UpdateDistbeta(ax, prob=0.3)
animbeta = FuncAnimation(fig, ud, frames=np.arange(100), init_func=ud.init,
interval=100, blit=True)
animbeta.save('moviebeta.mp4')
from scipy.stats import norm
import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
class UpdateDistnormal(object):
def __init__(self, ax, mu0, sigma0, sigma):
self.mu0 = mu0
self.sigma0 = sigma0
self.sigma = sigma
self.data = np.random.normal(self.mu0,self.sigma0,1)
self.cumdata = 0
self.line, = ax.plot([], [], 'k-',color='blue')
self.data_point = ax.text(0.02, 0.95, '', transform=ax.transAxes)
self.cumdata_point = ax.text(0.02, 0.90, '', transform=ax.transAxes)
self.true_value = ax.text(0.02, 0.85, '', transform=ax.transAxes)
self.bayes_estimator = ax.text(0.02,0.80, '',transform=ax.transAxes)
self.x = np.linspace(-10, 10, 1000)
self.ax = ax
# Set up plot parameters
self.ax.set_xlim(-10, 10)
self.ax.set_ylim(0, 5)
self.ax.grid(True)
self.ax.axvline(self.data, linestyle='--', color='red')
# This vertical line represents the theoretical value, to
# which the plotted distribution should converge.
def init(self):
self.line.set_data([], [])
return self.line,
def __call__(self, i):
# This way the plot can continuously run and we just keep
# watching new realizations of the process
if i == 0:
return self.init()
# Choose success based on exceed a threshold with a uniform pick
data = np.random.normal(self.data,self.sigma,1)
self.cumdata = self.cumdata + data
sigmaupdate = (self.sigma0**(-2)+self.sigma**(-2))**(-0.5)
self.mu0 = (data/self.sigma**2+self.mu0/self.sigma0**2)*sigmaupdate**2
self.sigma0 = sigmaupdate
y = norm.pdf(self.x, self.mu0, self.sigma0)
self.line.set_data(self.x, y)
average = self.cumdata/i
self.data_point.set_text('Data point = %.5f' % data)
self.cumdata_point.set_text('MLE = %.5f' % average)
self.true_value.set_text('True value = %.5f' % self.data)
self.bayes_estimator.set_text('Bayes estimator = %.5f' % self.mu0)
return self.line,
# Fixing random state for reproducibility
np.random.seed(19680802)
fig, ax = plt.subplots()
ud = UpdateDistnormal(ax, mu0 = 5.0, sigma0 = 1.0, sigma = 5.0)
animnormal= FuncAnimation(fig, ud, frames=np.arange(100), init_func=ud.init,
interval=100, blit=True)
animnormal.save('movienormal.mp4')