Login
📚 Statistical Inference for Everyone
Chapters ▾
1 Introduction to Probability2 Applications of Probability3 Random Sequences and Visualization4 Introduction to Model Comparison5 Applications of Model Comparison6 Introduction to Parameter Estimation7 Priors, Likelihoods, and Posteriors8 Common Statistical Significance Tests9 Applications of Parameter Estimation and Inference10 Multi-parameter Models11 Introduction to MCMC12 Concluding Thoughts13 Appendix A: Computational Analysis14 Appendix B: Notation and Standards15 Appendix C: Common Distributions and Their Properties16 Appendix D: Tables

11.2 One-Dimensional Models

The model we first look at is the coin flip model: given 17 heads in 25 flips, what is the probability distribution of the the measure of the coin's bent-ness, θ. We know the solution is of the form of a beta distribution, but we perform the same analysis with the MCMC technique.

h,N=data=17,25
def P_data(data,theta):
    h,N=data
    distribution=Bernoulli(h,N)
    return distribution(theta)
model=MCMCModel(data,P_data,
                theta=Uniform(0,1))
About a dozen overlaid MCMC chains for θ plotted against time step, 0 to 500. The traces start spread across the whole 0-1 range, then within roughly 25 steps settle into a dense band between about 0.4 and 0.9 centered near 0.7.
Figure 11.1. So-called MCMC “chains” for parameter θ versus time. Observe that the values of θ start spread evenly from 0 to 1 at the beginning and then thin down to a range of about 0.5-0.8 with the middle around 0.7 (17/25 = 0.68).
model.run_mcmc(500)
model.plot_chains()

Reading the Output

We can now plot the distributions of the parameters, just θ in this case, yielding best-fits, uncertainties, etc...

model.plot_distributions()

(see Figure 11.2)

Jagged sampled histogram of the posterior for θ (blue dots) with two smooth fitted curves over it, titled with the estimate θ=0.669 (+0.157/-0.187 for the 2.5-97.5 percentile range). The hump peaks near θ=0.67.
Figure 11.2. Distribution of θ, and the 95% credible interval.

We can further perform some simple calculations on the probabilities for the parameters, such as

model.P('theta>0.5')
Show expected output
0.96173333333333333
model.P('(0.2<theta) & (theta<.5)')
Show expected output
0.038266666666666664

Adapted from Statistical Inference for Everyone, by Brian Blais (Bryant University), licensed under CC BY-SA 4.0 (dual-licensed under the GNU FDL 1.2 or later; this adaptation uses the CC BY-SA grant). Changes were made; this adaptation is distributed under the same license. License: CC-BY-SA-4.0.