Bayesian parameter estimation and model comparison ==================================================== For the example scripts below to run, you will need to import the following modules: :: import numpy as np from numpy import sum import matplotlib.pyplot as plt from scipy.stats import distributions Is this coin fair? -------------------- Uniform prior ~~~~~~~~~~~~~~~ Imagine we had :math:`n=20, r=7`. We can plot the `posterior PDF `_ using :: n = 20 r = 7 Nsamp = 201 # no of points to sample at p = np.linspace(0,1,Nsamp) pdense = distributions.binom.pmf(r,n,p) # probability mass function plt.plot(p,pdense) plt.savefig('fig01.png') .. image:: fig01.png :scale: 50% Implementing this in Python and plotting :: pdense = distributions.binom.pmf(r,n,p) # probability mass function p_mean = sum(p*pdense)/sum(pdense) plt.axvline(p_mean,color='r') plt.savefig('fig02.png') plt.close() .. image:: fig02.png :scale: 50% It is instructive to repeat this example for :math:`n=20` for a range of values of :math:`r`. :: n = 20 Nsamp = 201 p = np.linspace(0,1,Nsamp) deltap = 1./(Nsamp-1) # step size between samples of p for r in range(0,20,2): # compute posterior pdense = distributions.binom.pmf(r,n,p) pdense = pdense / (deltap*sum(pdense)) # normalize posterior # compute mean p_mean = sum(p*pdense)/sum(pdense) # make the figure plt.figure() plt.plot(p,pdense) plt.axvline(p_mean,color='r') plt.xlabel('p') plt.ylabel('P(p|D,M)') plt.title('r={}'.format(r)) plt.savefig('fig03_r{}'.format(r)) plt.close() +--------------------------+------------------------------+-----------------------------+ | .. image:: fig03_r0.png | .. image:: fig03_r2.png | .. image:: fig03_r4.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig03_r6.png | .. image:: fig03_r8.png | .. image:: fig03_r10.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig03_r12.png | .. image:: fig03_r14.png | .. image:: fig03_r16.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig03_r18.png | | | | :width: 400px | | | +--------------------------+------------------------------+-----------------------------+ Beta prior ~~~~~~~~~~~~~~ Let's now redo the coin analysis with a `beta prior `_ with :math:`\alpha=\beta=10` :: alpha = 10 beta = 10 Nsamp = 201 # no of points to sample at p = np.linspace(0,1,Nsamp) pdense = distributions.beta.pdf(p,alpha,beta) # probability mass function plt.plot(p,pdense) plt.savefig('fig04.png') .. image:: fig04.png :scale: 50% :: alpha = 10 beta = 10 n = 20 Nsamp = 201 # no of points to sample at deltap = 1./(Nsamp-1) # step size between samples of p for r in range(0,20,2): # compute posterior pdense = distributions.beta.pdf(p,alpha+r,beta+n-r) pdense2 = distributions.binom.pmf(r,n,p) * distributions.beta.pdf(p,alpha,beta) pdense2 = pdense2 / (deltap*sum(pdense2)) # normalize posterior # make the figure plt.figure() plt.plot(p,pdense) plt.plot(p,pdense2,'--',lw=2) plt.xlabel('p') plt.ylabel('PDF') plt.title('r={}'.format(r)) plt.savefig('fig05_r{}'.format(r)) +--------------------------+------------------------------+-----------------------------+ | .. image:: fig05_r0.png | .. image:: fig05_r2.png | .. image:: fig05_r4.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig05_r6.png | .. image:: fig05_r8.png | .. image:: fig05_r10.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig05_r12.png | .. image:: fig05_r14.png | .. image:: fig05_r16.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig05_r18.png | | | | :width: 400px | | | +--------------------------+------------------------------+-----------------------------+ It is instructive to plot the likelihood, prior and posterior together: :: alpha = 10 beta = 10 n = 20 Nsamp = 201 # no of points to sample at deltap = 1./(Nsamp-1) # step size between samples of p prior = distributions.beta.pdf(p,alpha,beta) for r in range(0,20,2): like = distributions.binom.pmf(r,n,p) like = like/(deltap*sum(like)) # for plotting convenience only: see below post = distributions.beta.pdf(p,alpha+r,beta+n-r) # make the figure plt.figure() plt.plot(p,post,'k',label='posterior') plt.plot(p,like,'r',label='likelihood') plt.plot(p,prior,'b',label='prior') plt.xlabel('p') plt.ylabel('PDF') plt.legend(loc='best') plt.title('r={}'.format(r)) plt.savefig('fig06_r{}'.format(r)) +--------------------------+------------------------------+-----------------------------+ | .. image:: fig06_r0.png | .. image:: fig06_r2.png | .. image:: fig06_r4.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig06_r6.png | .. image:: fig06_r8.png | .. image:: fig06_r10.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig06_r12.png | .. image:: fig06_r14.png | .. image:: fig06_r16.png | | :width: 400px | :width: 400px | :width: 400px | +--------------------------+------------------------------+-----------------------------+ | .. image:: fig06_r18.png | | | | :width: 400px | | | +--------------------------+------------------------------+-----------------------------+ It is interesting to see how the likelihood and posterior evolve as we get more data. :: alpha = 10 beta = 10 n = 20 Nsamp = 201 # no of points to sample at deltap = 1./(Nsamp-1) # step size between samples of p prior = distributions.beta.pdf(p,alpha,beta) for i in range(1,9): r = 2**i n = (3.0/2.0)*r like = distributions.binom.pmf(r,n,p) like = like/(deltap*sum(like)) # for plotting convenience only post = distributions.beta.pdf(p,alpha+r,beta+n-r) # make the figure plt.figure() plt.plot(p,post,'k',label='posterior') plt.plot(p,like,'r',label='likelihood') plt.plot(p,prior,'b',label='prior') plt.xlabel('p') plt.ylabel('PDF') plt.legend(loc='best') plt.title('r/n={}/{:.0f}'.format(r,n)) plt.savefig('fig07_r{}'.format(r)) +----------------------------+------------------------------+-----------------------------+ | .. image:: fig07_r2.png | .. image:: fig07_r4.png | .. image:: fig07_r8.png | | :width: 400px | :width: 400px | :width: 400px | +----------------------------+------------------------------+-----------------------------+ | .. image:: fig07_r16.png | .. image:: fig07_r32.png | .. image:: fig07_r64.png | | :width: 400px | :width: 400px | :width: 400px | +----------------------------+------------------------------+-----------------------------+ | .. image:: fig07_r128.png | .. image:: fig07_r256.png | | | :width: 400px | :width: 400px | | +----------------------------+------------------------------+-----------------------------+