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.3 Multi-Dimensional Models

It is straightforward then to include more than one parameter and to do regression using this technique. For example, here is an example with some artificial data,

def linear(x,a,b):
    return a*x+b

model=MCMCModel_Regression(x,y,linear,
                a=Uniform(-10,10),
                b=Uniform(0,100),
                )
Three stacked chain plots for the straight-line model parameters a, b and the noise σ over 500 steps. Each set of traces starts widely scattered (b as far as 0-100, σ in the thousands) and collapses within about 100 steps into narrow horizontal bands near a≈0.3, b≈40 and small σ.
Figure 11.3. Chains for parameters a, b, and the noise σ.
model.run_mcmc(500)
model.plot_chains()
plot(x,y,'o')
model.plot_predictions(xfit,color='g')
Scatter of about 1,000 blue data points between roughly 37 and 43.5 with a dense green band of model predictions running horizontally near 40; the band's thinness against the scatter shows the fit's small uncertainty.
Figure 11.4. Data (blue) and predictions (green) for the model - the width of the predictions demonstrates the uncertainty.
model.plot_distributions()
Sampled posterior for the slope parameter a: a noisy blue histogram trace with a smooth green curve over it, titled a=0.249 (+0.210/-0.216), peaking near 0.25 and spanning about -0.1 to 0.6.
Figure 11.5. Distributions for parameters a and b (slope and intercept).

And we can look at best estimates, quartiles, and probability comparisons,

model.percentiles([5,50,95])
Show expected output
{'_sigma': array([ 0.97143798,  1.00744104,  1.0467333 ]),
 'a': array([ 0.07063144,  0.24939562,  0.42523751]),
 'b': array([ 39.88446461,  39.98633744,  40.09010139])}
model.P('a>0')
Show expected output
0.98936000000000002

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.