︠11b5544e-e209-4a0b-927b-1bea51a863b3si︠ %md # PyMC Demo First off: In case you are intersted in Bayesian Statistical Methods for Python, start with --> **[Probabilistic Programming & Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/)** <-- The following is taken from this book and [here](http://nbviewer.ipython.org/github/pymc-devs/pymc/blob/master/pymc/examples/tutorial.ipynb) 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. ︡34394a4f-13fd-4199-8aeb-e8bf901e888f︡{"md":"# PyMC Demo\n\nFirst off: In case you are intersted in Bayesian Statistical Methods for Python, start with\n\n--> **[Probabilistic Programming & Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/)** <--\n\nThe following is taken from this book and [here](http://nbviewer.ipython.org/github/pymc-devs/pymc/blob/master/pymc/examples/tutorial.ipynb)\nshows how it works inside a Sage Worksheet.\n\nThe only caveat is, that this is installed from the master branch on github, which is version 3.\nTo get the older 2.3, just install pymc locally.\n\nGeneral remark: Sage converts numbers to it's basic Types (\"Integer\", ...).\nTo 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."}︡ ︠d7564ea9-7140-4b89-bc98-b627f0e4ad06s︠ import pymc as pm import numpy as np import matplotlib.pyplot as plt pm.__version__ ︡926f5b21-0cbe-4db6-8595-282d3f297981︡{"stdout":"'3.0'\n"}︡ ︠3cac1d0b-960b-43a2-8372-aa2840864d9fs︠ 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)) ︡efed6316-e1c2-46b6-bcc0-f6cb6175f204︡ ︠efc670be-9236-469f-af4c-344b6c60149bs︠ 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 ︡85d0934b-166d-4270-ac05-effe080dddb9︡ ︠40f0ae1a-5965-401a-bc65-632628d73e47s︠ 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) ︡28ffa512-9488-4f10-a901-f58a092135e6︡{"stdout":"MAP found:\n"}︡{"stdout":"('x:', array(0.42206175558831205))\n"}︡{"stdout":"('y:', array([[ 2.00539329],\n [ 1.7133976 ]]))\n"}︡{"stdout":"Compare with true values:\n"}︡{"stdout":"('ytrue', array([[ 1.88172061],\n [ 1.70898794]]))\n"}︡{"stdout":"('xtrue', array([ 0.45289684]))\n"}︡ ︠2d3803d9-09cf-4508-be25-3b23bc162ab0s︠ with model: step = pm.NUTS() ︡7ba91e7e-b5f3-47ac-9186-591913eafa92︡ ︠a207ccad-e5f0-4bce-959a-cb93babc27c7s︠ with model: trace = pm.sample(3000, step, start) ︡01e1177b-cb09-45ac-b9c7-865a1c1b57eb︡{"stdout":"[- 4% ] 132 of 3000 complete in 0.5 sec"}︡{"stdout":"\n[--- 9% ] 280 of 3000 complete in 1.0 sec"}︡{"stdout":"\n[----- 14% ] 441 of 3000 complete in 1.5 sec"}︡{"stdout":"\n[------- 20% ] 623 of 3000 complete in 2.0 sec"}︡{"stdout":"\n[--------- 26% ] 785 of 3000 complete in 2.5 sec"}︡{"stdout":"\n[----------- 31% ] 932 of 3000 complete in 3.0 sec"}︡{"stdout":"\n[------------- 36% ] 1084 of 3000 complete in 3.5 sec"}︡{"stdout":"\n[--------------- 41% ] 1258 of 3000 complete in 4.0 sec"}︡{"stdout":"\n[-----------------46% ] 1408 of 3000 complete in 4.5 sec"}︡{"stdout":"\n[-----------------51% ] 1555 of 3000 complete in 5.0 sec"}︡{"stdout":"\n[-----------------57%- ] 1713 of 3000 complete in 5.5 sec"}︡{"stdout":"\n[-----------------62%--- ] 1874 of 3000 complete in 6.0 sec"}︡{"stdout":"\n[-----------------67%----- ] 2031 of 3000 complete in 6.5 sec"}︡{"stdout":"\n[-----------------73%------- ] 2202 of 3000 complete in 7.0 sec"}︡{"stdout":"\n[-----------------78%--------- ] 2365 of 3000 complete in 7.5 sec"}︡{"stdout":"\n[-----------------84%----------- ] 2524 of 3000 complete in 8.1 sec"}︡{"stdout":"\n[-----------------89%------------- ] 2679 of 3000 complete in 8.6 sec"}︡{"stdout":"\n[-----------------94%--------------- ] 2835 of 3000 complete in 9.1 sec"}︡{"stdout":"\n[-----------------99%----------------- ] 2998 of 3000 complete in 9.6 sec"}︡{"stdout":"\n[-----------------100%-----------------] 3000 of 3000 complete in 9.6 sec\n"}︡ ︠5aac1396-69d4-45ff-acf0-e0e4d1d30008s︠ pm.traceplot(trace) ︡562e8e3b-233f-4d00-8127-28ad0291d914︡{"once":false,"file":{"show":true,"uuid":"9d69c81a-f99b-418f-bcad-1419902bcacf","filename":"/projects/20e4a191-73ea-4921-80e9-0a5d792fc511/.sage/temp/compute7dc2/15416/tmp_qDKNsT.svg"}}︡ ︠e6558934-75ae-4e07-8407-7f4ff9971be2︠