Open in CoCalc with one click!

PyMC Demo

First off: In case you are intersted in Bayesian Statistical Methods for Python, start with

--> Probabilistic Programming & Bayesian Methods for Hackers <--

The following is taken from this book and here shows how it works inside a Sage Worksheet.

The only caveat is, that this is installed from the master branch on github, which is version 3. To get the older 2.3, just install pymc locally.

General remark: Sage converts numbers to it's basic Types ("Integer", ...). To avoid this, 1r is a "raw" Python 1-int, or %python or even %default_mode python sets the cell or everything to pure Python mode.

import pymc as pm import numpy as np import matplotlib.pyplot as plt pm.__version__
'3.0'
ndims = 2 nobs = 20 xtrue = np.random.normal(scale=2., size=1) ytrue = np.random.normal(loc=np.exp(xtrue), scale=1, size=(ndims, 1)) zdata = np.random.normal(loc=xtrue + ytrue, scale=.75, size=(ndims, nobs))
with pm.Model() as model: x = pm.Normal('x', mu=0., sd=1) y = pm.Normal('y', mu=pm.exp(x), sd=2., shape=(ndims, 1)) # here, shape is telling us it's a vector rather than a scalar. z = pm.Normal('z', mu=x + y, sd=.75, observed=zdata) # shape is inferred from zdata
with model: start = pm.find_MAP() print("MAP found:") print("x:", start['x']) print("y:", start['y']) print("Compare with true values:") print("ytrue", ytrue) print("xtrue", xtrue)
MAP found: ('x:', array(0.42206175558831205)) ('y:', array([[ 2.00539329], [ 1.7133976 ]])) Compare with true values: ('ytrue', array([[ 1.88172061], [ 1.70898794]])) ('xtrue', array([ 0.45289684]))
with model: step = pm.NUTS()
with model: trace = pm.sample(3000, step, start)
[- 4% ] 132 of 3000 complete in 0.5 sec [--- 9% ] 280 of 3000 complete in 1.0 sec [----- 14% ] 441 of 3000 complete in 1.5 sec [------- 20% ] 623 of 3000 complete in 2.0 sec [--------- 26% ] 785 of 3000 complete in 2.5 sec [----------- 31% ] 932 of 3000 complete in 3.0 sec [------------- 36% ] 1084 of 3000 complete in 3.5 sec [--------------- 41% ] 1258 of 3000 complete in 4.0 sec [-----------------46% ] 1408 of 3000 complete in 4.5 sec [-----------------51% ] 1555 of 3000 complete in 5.0 sec [-----------------57%- ] 1713 of 3000 complete in 5.5 sec [-----------------62%--- ] 1874 of 3000 complete in 6.0 sec [-----------------67%----- ] 2031 of 3000 complete in 6.5 sec [-----------------73%------- ] 2202 of 3000 complete in 7.0 sec [-----------------78%--------- ] 2365 of 3000 complete in 7.5 sec [-----------------84%----------- ] 2524 of 3000 complete in 8.1 sec [-----------------89%------------- ] 2679 of 3000 complete in 8.6 sec [-----------------94%--------------- ] 2835 of 3000 complete in 9.1 sec [-----------------99%----------------- ] 2998 of 3000 complete in 9.6 sec [-----------------100%-----------------] 3000 of 3000 complete in 9.6 sec
pm.traceplot(trace)