Skip to content

Instantly share code, notes, and snippets.

@fonnesbeck

fonnesbeck/Cox Model (PyMC 3).ipynb

Created July 13, 2015 15:14
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save fonnesbeck/7b99108ba23703bbe48b to your computer and use it in GitHub Desktop.
Save fonnesbeck/7b99108ba23703bbe48b to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
Cox Model (PyMC 3).ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/fonnescj/GitHub/notebook/src/ipython/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.\n",
" \"`IPython.html.widgets` has moved to `ipywidgets`.\", ShimWarning)\n",
"/Users/fonnescj/GitHub/notebook/src/ipython/IPython/utils/traitlets.py:5: UserWarning: IPython.utils.traitlets has moved to a top-level traitlets package.\n",
" warn(\"IPython.utils.traitlets has moved to a top-level traitlets package.\")\n"
]
}
],
"source": [
"%matplotlib inline\n",
"from pymc3 import Gamma, Poisson, Normal, Model, sample, forestplot, NUTS, Metropolis, find_MAP, starting, Slice\n",
"import theano.tensor as tt\n",
"from theano import function as fn\n",
"from theano import printing\n",
"import numpy as np\n",
"import scipy as sp"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is the original model, implemented in BUGS:\n",
"\n",
" model\n",
" { \n",
" # Set up data\n",
" for(i in 1:Nsubj) {\n",
" for(j in 1:T) {\n",
" # risk set = 1 if obs.t >= t\n",
" Y[i,j] <- step(obs.t[i] - t[j] + eps)\n",
" # counting process jump = 1 if obs.t in [ t[j], t[j+1] )\n",
" # i.e. if t[j] <= obs.t < t[j+1]\n",
" dN[i, j] <- Y[i, j] * step(t[j + 1] - obs.t[i] - eps) * FAIL[i]\n",
" }\n",
" }\n",
"\n",
" # Model \n",
" for(j in 1:T) {\n",
" for(i in 1:Nsubj) {\n",
" dN[i, j] ~ dpois(Idt[i, j]) # Likelihood\n",
" Idt[i, j] <- Y[i, j] * exp(beta[1]*pscenter[i] + beta[2]*\n",
" hhcenter[i] + beta[3]*ncomact[i] + beta[4]*rleader[i] + beta[5]*dleader[i] + beta[6]*inter1[i] + beta[7]*inter2[i]) * dL0[j] # Intensity\n",
" } \n",
" dL0[j] ~ dgamma(mu[j], c)\n",
" mu[j] <- dL0.star[j] * c # prior mean hazard \n",
" }\n",
"\n",
"\n",
" c ~ dgamma(0.0001, 0.00001)\n",
" r ~ dgamma(0.001, 0.0001)\n",
"\n",
"\n",
" for (j in 1 : T) { dL0.star[j] <- r * (t[j + 1] - t[j]) } \n",
" # next line indicates number of covariates and is for the corresponding betas\n",
" for(i in 1:7) {beta[i] ~ dnorm(0.0,0.00001)} \n",
"\n",
"\n",
" }"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dta = dict(T=73, Nsubj=430, eps=0.0, t=[1, 21, 85, 128, 129, 148, 178, 204,\n",
" 206, 210, 211, 212, 225, 238, 241,\n",
" 248, 259, 273, 275, 281, 286, 289,\n",
" 301, 302, 303, 304, 313, 317, 323,\n",
" 344, 345, 349, 350, 351, 355, 356,\n",
" 359, 364, 385, 386, 389, 390, 391,\n",
" 392, 394, 395, 396, 397, 398, 399,\n",
" 400, 406, 415, 416, 426, 427, 434,\n",
" 435, 437, 441, 447, 448, 449, 450,\n",
" 451, 453, 455, 456, 458, 459, 460,\n",
" 461, 462, 463],\n",
"obs_t = [460, 313, 435, 350, 435, 350, 350, 460, 460, 448, 225, 225, 396, 435, 396, 396, 453, 396, 456, 397, 397, 396, 395, 275, 449, 395, 395, 462, 302, 302, 458, 461, 396, 241, 389, 458, 304, 304, 395, 395, 364, 460, 415, 463, 396, 459, 441, 435, 396, 458, 437, 396, 356, 356, 396, 455, 396, 462, 399, 400, 350, 350, 395, 395, 441, 355, 85, 458, 128, 396, 386, 386, 386, 462, 458, 390, 390, 396, 396, 396, 427, 458, 395, 275, 275, 395, 359, 395, 395, 441, 395, 463, 178, 275, 463, 396, 396, 259, 396, 396, 458, 441, 396, 463, 396, 463, 435, 396, 437, 396, 398, 463, 460, 462, 460, 460, 210, 396, 435, 458, 385, 323, 323, 359, 396, 396, 460, 238, 441, 450, 392, 458, 396, 458, 396, 396, 462, 435, 396, 394, 396, 435, 458, 1, 395, 395, 451, 462, 458, 462, 396, 286, 396, 349, 449, 462, 455, 21, 463, 461, 461, 456, 435, 396, 460, 462, 462, 435, 435, 460, 386, 396, 458, 386, 461, 441, 435, 435, 463, 456, 396, 275, 460, 406, 460, 406, 317, 406, 461, 396, 359, 458, 463, 435, 462, 458, 396, 396, 273, 396, 435, 281, 275, 396, 447, 225, 447, 396, 435, 416, 396, 248, 396, 435, 435, 396, 461, 385, 396, 458, 458, 396, 461, 396, 448, 396, 396, 460, 455, 456, 463, 462, 458, 463, 396, 462, 395, 456, 396, 463, 396, 435, 459, 396, 396, 396, 395, 435, 455, 395, 461, 344, 396, 395, 396, 317, 396, 395, 426, 461, 396, 289, 441, 395, 396, 458, 396, 396, 435, 396, 395, 396, 441, 345, 396, 359, 435, 435, 396, 396, 395, 458, 461, 458, 212, 301, 458, 456, 395, 396, 395, 435, 396, 396, 303, 458, 460, 400, 396, 462, 359, 458, 396, 206, 441, 396, 458, 396, 462, 396, 396, 275, 396, 395, 435, 435, 462, 225, 458, 462, 396, 396, 289, 396, 303, 455, 400, 400, 359, 461, 396, 462, 460, 463, 463, 463, 204, 435, 435, 396, 396, 396, 463, 458, 396, 455, 435, 396, 396, 463, 396, 461, 463, 460, 441, 460, 435, 435, 460, 455, 460, 395, 460, 460, 460, 435, 449, 463, 462, 129, 391, 396, 391, 391, 434, 356, 462, 396, 349, 225, 396, 435, 461, 391, 391, 351, 211, 461, 212, 434, 148, 356, 458, 456, 455, 435, 463, 463, 462, 435, 463, 437, 460, 396, 406, 451, 460, 435, 396, 460, 455, 396, 398, 456, 458, 396, 456, 449, 396, 128, 396, 462, 463, 396, 396, 396, 435, 460, 396, 458],\n",
"FAIL= [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],\n",
"pscenter= [\n",
" -.01434325, -.01460965, .01322687, .00971885, -.03223412, -.01113493, -.01359567, -.03357866, -.0387039, -.0553269, -.03238896, -.07464545, -.07325128, -.07062459, -.07464545, -.07032613, -.0703005, .00965232, -.01408955, .00577483, -.00219072, -.00084567, .01643198, .06509522, .06824313, .07300876, .07300876, .01394272, .06824313, .02063087, .00383186, -.02573045, -.02410864, -.02272752, .05120398, -.00997729, -.00550709, -.02062663, -.03077685, -.01688493, .01035959, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .01149963, .0034338, .0376236, .00733331, .01520069, .03832785, .03832785, -.02622275, -.02622275, -.02622275, -.01492678, -.02897806, -.02897806, -.02897806, -.02847666, -.031893, -.03919478, -.04224754, -.04743705, -.0510477, -.031893, -.01129093, .01706207, .00193999, -.01503116, .003101, -.00083466, .02395027, -.07952866, -.08559135, -.07251801, -.06586029, -.08432532, -.0613939, -.081205, -.07540084, -.08488011, -.08488011, -.08488011, -.07492433, -.08907269, -.09451609, -.05301854, -.08980743, -.0771635, -.0771635, -.08650947, -.07856082, -.0771635, -.08204606, -.08178245, -.05263504, -.05355574, -.05109092, -.04696729, -.04696729, -.04696729, -.05257489, -.05303248, -.05348096, -.04983674, -.04699414, .00584956, -.00792241, -.01719816, -.02138029, -.01576016, -.04274812, -.04014061, .0471441, .0471441, .0471441, .0471441, .0471441, .0471441, .0471441, .04233112, .0471441, .04233112, .050568, .07388823, .0493324, .04512087, .03205975, .02913185, .06010427, .05324252, .06973204, .05579907, .01212243, .07962459, .05054695, .06672142, .14026688, .01734403, .06078221, .06543709, .06438115, .20126908, -.03138622, -.02180659, .01637333, -.02415774, .01828684, .03106104, .04268495, .01897239, .01591935, -.02367065, -.0619156, -.06403028, -.06851645, -.04821694, -.03889525, -.05023452, -.05013452, -.01557191, -.01171948, -.01362136, -.01174715, -.02707938, -.02634164, -.02634164, -.02634164, -.00692153, -.02381614, -.00890537, -.00611669, -.00894752, -.03551984, -.0252678, -.01513384, -.01016569, -.03551984, -.03773227, -.01978032, .06803483, .06706496, .10551275, .15091534, .03092981, .06556855, .10781559, .12671031, .0936299, .09362991, .09362991, .08294538, .09362991, .09362991, .09362991, .01177025, .02610553, .03546937, .03546937, .03546937, .034415, -.00305626, .04973665, .05103208, .07546701, .05306436, .00824125, .01961115, .01202359, -.02919447, -.01016712, .01756074, -.04035511, -.04753104, -.04463152, -.04845615, -.05010044, .00031411, -.07911871, -.08799869, -.07980882, -.09393142, -.08000018, -.07666632, -.07817401, -.07444922, -.07226554, -.08216553, -.0777643, -.07752042, -.05767992, -.04727952, -.03774814, -.06870384, -.05999847, -.05947695, .02989959, .04627543, .02772475, .02883079, .03642944, .02871235, .04148949, .04240279, .07747082, .07626323, .04268012, .03225577, .06468724, -.05140995, -.05399637, -.05351515, .07302427, .02432223, .0490674, .0490674, .0490674, .0490674, .09013112, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .10476315, .07008056, .08666077, .01546215, .01667466, .03417671, .05253941, .04293926, .01496588, .02692172, -.03827151, .04809769, .08742411, .04533176, .01455173, .01831875, .02710811, .09834951, .09952456, .06993483, .02945534, .038731, .1181948, .04435538, .04435538, -.02357505, .05824019, .05820741, -.02357505, .09324722, .15534712, .07207468, .04692869, -.03490683, -.04404809, -.05054474, -.05325826, -.0474724, -.04905931, .01068221, .02879751, .00852646, .02693032, .01835589, .02989959, .02989959, .02989959, .04976377, .04439012, .03397319, .02989959, .02989959, .05468828, .04463226, .05886378, .06311052, .02989959, .04595331, .04203459, .01231324, -.01399783, .04595331, .00145386, .04601278, .06459354, -.0007196, .00012216, -.07614055, -.08435525, -.07957162, -.10299519, -.08156988, -.08225659, -.07449063, -.00210284, -.00797183, -.025355, -.01258251, -.04372031, -.03985972, -.03545086, -.03384566, -.04025533, -.07523724, -.05947702, -.061286, -.07666647, -.07663169, -.05902354, -.07652324, -.07645561, -.06258684, -.09604834, -.08813326, -.03292062, -.07848112, -.08239502, -.08316891, -.07244316, -.075417, -.07652324, -.07922532, -.08755959, -.08583414, -.07450142, -.08066016, -.06057205, -.07652324, -.06249051, -.08781742, -.086076, -.07652324, -.07696518, -.0618688, -.06073988, -.06524737, -.04419825, -.04489509, -.04390368, -.04358438, -.04489509, -.04520512, -.04187583, -.03653955, -.03973426, -.03753508, -.03569439, -.06789339, .06689456, .05526327, .05139003, .02641841, .04891529, .07078697, .06862645, .06832582, .04104258, -.00120631, .01947345, .04891779, .04891779, .03561932, .02576244, .03158225, .03608047, .08685057, .04632537, .06841581, -.02899643],\n",
"hhcenter= [ -.78348798, -.63418788, -.91218799, -.98388809, -.23518796, .11481193, -1.415588, -1.2535881, -.55738801, -.88128799, -1.109488, .05721192, -1.045788, -.30888793, .29651192, -.36688802, -.50058788, .02271203, -.59088796, -.04198809, .50561196, -.07418796, .98481184, .78921205, .09431199, -.06488796, 2.1662121, .08891205, 1.4004121, 1.316112, 1.9362121, 2.0107121, 1.150712, .31951192, -.23918791, -.1562881, -.9575879, -.07728811, .29641202, 1.2273121, 1.7717118, 1.5764117, .14181189, .72131211, 1.279212, .68241197, -.72808808, -.00488802, -.23938794, -1.000788, .55081207, -.52348799, 1.780612, -.35888812, .36481193, 1.5480118, -.03078791, 1.389112, .30211189, .70901209, -.16668792, 1.435812, .47001198, 2.0838118, 1.1673121, .18461208, -.30608794, 1.4470119, .23301201, -.58458799, .44011191, -.61948794, -.41388795, .263212, .66171199, .92451197, .78081208, .90991193, 1.6920118, 1.334012, 1.2101121, .41591194, -.48498794, -.73278803, -1.093588, .09911207, -.93418807, -.46908805, .0205119, .0535119, -.14228792, -.55708808, -.45498797, -.54008788, -.30998799, -.10958811, -.0960879, -.01338812, -.88168806, -.51788801, .36801198, .46621206, .13271193, -.11208793, -.76768798, -.54508799, -1.2773881, .16641192, .95871216, -.48238799, 1.6281118, -.18848796, -.49718806, -.41348812, -.31628796, -.59528798, -.11718794, -.57058805, -.59488791, -.21248789, -.65658802, -.56298798, -.52698797, -.65758795, -.04988809, .55341202, -.76328796, .254612, 1.3500118, -.54958791, 1.665812, .14671211, 1.963912, .29161194, -.56838793, 1.9371119, .90991193, -.39558789, .39521196, -.55208796, -.05268808, -.77368802, -.45428798, .05841212, -.45308802, -.12458798, .01431207, -.28228804, .79281193, -.26358792, -.54738802, -.38158795, -.54118794, -.72828788, -.58128804, .355912, -.24078794, -1.0384881, -.75038809, -.41018793, -.43538806, -1.566388, -.53388804, -.28388807, -1.2348881, -.69028801, -1.620088, -.78128809, -.54648799, -.92738789, .11871199, .26851204, .61571199, .82891208, 1.1985121, 1.012012, 1.0602121, -.02988811, .79301196, .67731196, .43991187, .9404121, .5254119, 1.0365119, 1.6220121, .61671191, -.50318807, 2.6073117, .02361206, -.60438794, -.79278797, -.18108793, -.48178813, -.44038793, -.22628804, -.07398792, .519512, .40211204, .582012, 1.830512, .80441195, .58801204, -.56368798, -1.5451881, .45991209, -.23448797, -.36918804, 1.3247118, .19541197, -.20818801, 1.163012, -.78228801, -.6048879, -.575288, 1.3241119, .0147119, -.76518792, -.37478802, -.35508797, -.90038794, -1.250888, -.46608803, -.98488802, -1.5185881, -.90908808, -1.048188, -.90138787, -.77278799, -1.248988, -.34448811, -.61628789, .38531187, -.51728791, -.00878807, -.60078806, -.45358798, .46301201, -.22048803, -.71518797, -.76478809, -.75028795, -.4952881, .01731209, -.83718795, .57951194, .54291207, .45341209, .16941194, 1.054112, .61721212, 2.2717118, 1.1593118, 2.0280118, .92281204, 1.0100121, -.1866879, 2.6503119, 2.3914118, -.19948788, -.36418793, -.9259879, -.71058792, -.1104879, .16971211, 1.474812, 1.9360118, 2.5344119, 2.0171118, 1.9387121, .55071193, -.03918811, .20681195, .40421203, -.75518793, -.45678803, -1.0271881, .77211195, 1.146812, -1.147788, -1.565588, -.34888789, 1.303812, 1.952312, 1.639112, .07731203, .25901201, -.45608804, -.5028879, .03641204, -.03808804, .38571194, .31831196, -.17648788, -.44528791, -.55918807, -.53108805, .39721206, -.06328794, -.34038803, -.05988808, -.89548796, -.03518792, .045512, -.1859879, -.039288, -.82568806, .01431207, .40091208, -.2531881, .030412, -.31918809, -.54958791, -.79078788, .36691192, -.324388, -1.0082881, -1.232188, -.53248805, -.23678799, -.89188808, .25111201, -.6766879, -.3565881, -.61228794, -.21078797, -1.0343881, -.58358806, -.15588804, -.39238808, -.67818803, -.19498797, 1.099412, 1.2767119, -.64068788, -.50678796, -.64058799, -.86918801, 1.4048119, -.59648794, .23331194, .68371207, .11251191, -.17128797, .17081194, -.44218799, -.48708794, .09591202, .20131211, -.20108791, -.02158805, -.48188803, -.3012881, -.55008787, -1.146188, -.82128805, -.87638801, -.54488796, -.60288805, -1.003088, -.25078794, -.14818807, -.14738794, -.80938786, -.85988802, -.90188807, -.94998807, -.75718802, -.37418792, -.66708797, 1.0981121, 1.1441121, .47381189, -.12958808, -.34358808, -.84328789, -.33498809, -.98088807, -.6903879, -1.284988, -.80838794, -.91838807, -.81848806, -.34488794, -.83438796, .12971191, .99381214, -.91608804, -.31808802, -.01018806, .98171192, -.91638798, -1.043988, -1.0103881, 1.451612, -.01528808, .02441196, -.41458794, .25691202, .18601207, -.815988, -.02908798, -.59088796, -.35608789, .79691201, 1.8123121, -.98588794, 1.548912, 2.3653121, -.09238812, .96741205, .05891208, -.15618797, -.5660879, -.28338811, -.10088798, 1.1663117, .21981196, .07151202, -.009088, -.49578807, .15441208, -.44488809, -.2677879, -.54388803, -.25468799, .68631202, -.88128799, -.84628791, -1.2549881, -.36198804],\n",
"ncomact= [ 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1],\n",
"rleader= [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
"dleader= [ 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
"inter1= [ -.01434325, -.01460965, 0, 0, 0, -.01113493, 0, 0, 0, -.0553269, -.03238896, 0, 0, -.07062459, -.07464545, -.07032613, 0, 0, -.01408955, 0, -.00219072, 0, 0, 0, 0, 0, .07300876, .01394272, 0, 0, 0, 0, 0, 0, .05120398, 0, -.00550709, -.02062663, -.03077685, -.01688493, 0, .01149963, 0, .01149963, .01149963, 0, 0, 0, 0, 0, 0, 0, 0, 0, .01149963, .0034338, .0376236, .00733331, 0, .03832785, .03832785, -.02622275, -.02622275, -.02622275, -.01492678, 0, 0, -.02897806, -.02847666, 0, 0, -.04224754, -.04743705, -.0510477, -.031893, 0, 0, 0, -.01503116, .003101, -.00083466, .02395027, -.07952866, 0, 0, -.06586029, 0, -.0613939, -.081205, -.07540084, -.08488011, -.08488011, 0, -.07492433, -.08907269, -.09451609, 0, -.08980743, 0, -.0771635, 0, 0, -.0771635, -.08204606, 0, -.05263504, 0, -.05109092, -.04696729, 0, -.04696729, 0, -.05303248, -.05348096, 0, 0, .00584956, -.00792241, -.01719816, 0, -.01576016, 0, -.04014061, 0, 0, 0, 0, 0, .0471441, 0, .04233112, 0, .04233112, 0, 0, .0493324, .04512087, .03205975, .02913185, 0, .05324252, 0, 0, 0, 0, .05054695, 0, .14026688, .01734403, .06078221, 0, 0, 0, -.03138622, 0, .01637333, 0, 0, 0, 0, .01897239, .01591935, 0, -.0619156, 0, -.06851645, 0, -.03889525, -.05023452, -.05013452, 0, 0, -.01362136, 0, 0, -.02634164, 0, 0, 0, 0, -.00890537, -.00611669, 0, 0, 0, -.01513384, 0, -.03551984, 0, -.01978032, 0, .06706496, .10551275, 0, .03092981, .06556855, 0, 0, 0, .09362991, 0, 0, 0, 0, 0, 0, .02610553, .03546937, 0, 0, .034415, 0, 0, 0, .07546701, 0, 0, 0, 0, -.02919447, -.01016712, 0, 0, 0, 0, -.04845615, -.05010044, 0, 0, 0, 0, 0, 0, -.07666632, 0, 0, -.07226554, -.08216553, -.0777643, 0, 0, -.04727952, 0, -.06870384, -.05999847, 0, 0, 0, .02772475, .02883079, .03642944, 0, .04148949, 0, 0, 0, .04268012, .03225577, 0, -.05140995, -.05399637, 0, 0, .02432223, 0, .0490674, .0490674, .0490674, 0, 0, 0, 0, 0, 0, 0, 0, .10476315, 0, 0, 0, 0, 0, .07008056, 0, 0, .01667466, 0, .05253941, .04293926, 0, .02692172, 0, 0, .08742411, .04533176, 0, .01831875, 0, .09834951, .09952456, 0, .02945534, .038731, 0, .04435538, 0, -.02357505, 0, 0, -.02357505, .09324722, 0, 0, 0, -.03490683, 0, -.05054474, 0, -.0474724, -.04905931, 0, .02879751, 0, 0, 0, 0, 0, 0, 0, .04439012, 0, .02989959, .02989959, .05468828, .04463226, 0, 0, 0, 0, 0, .01231324, -.01399783, .04595331, .00145386, 0, .06459354, -.0007196, 0, -.07614055, -.08435525, 0, -.10299519, 0, 0, 0, -.00210284, -.00797183, 0, 0, 0, 0, -.03545086, 0, 0, 0, 0, -.061286, -.07666647, 0, -.05902354, -.07652324, -.07645561, 0, 0, 0, -.03292062, 0, 0, 0, 0, -.075417, 0, -.07922532, 0, -.08583414, -.07450142, -.08066016, 0, 0, -.06249051, 0, 0, 0, 0, -.0618688, 0, -.06524737, -.04419825, -.04489509, 0, 0, 0, -.04520512, -.04187583, 0, 0, -.03753508, 0, 0, 0, 0, 0, 0, 0, 0, .06862645, 0, 0, -.00120631, .01947345, 0, 0, .03561932, 0, .03158225, .03608047, 0, 0, 0, -.02899643],\n",
"\n",
"inter2= [-.78348798, -.63418788, 0, 0, 0, .11481193, 0, 0, 0, -.88128799, -1.109488, 0, 0, -.30888793, .29651192, -.36688802, 0, 0, -.59088796, 0, .50561196, 0, 0, 0, 0, 0, 2.1662121, .08891205, 0, 0, 0, 0, 0, 0, -.23918791, 0, -.9575879, -.07728811, .29641202, 1.2273121, 0, 1.5764117, 0, .72131211, 1.279212, 0, 0, 0, 0, 0, 0, 0, 0, 0, .36481193, 1.5480118, -.03078791, 1.389112, 0, .70901209, -.16668792, 1.435812, .47001198, 2.0838118, 1.1673121, 0, 0, 1.4470119, .23301201, 0, 0, -.61948794, -.41388795, .263212, .66171199, 0, 0, 0, 1.6920118, 1.334012, 1.2101121, .41591194, -.48498794, 0, 0, .09911207, 0, -.46908805, .0205119, .0535119, -.14228792, -.55708808, 0, -.54008788, -.30998799, -.10958811, 0, -.01338812, 0, -.51788801, 0, 0, .13271193, -.11208793, 0, -.54508799, 0, .16641192, .95871216, 0, 1.6281118, 0, -.49718806, -.41348812, 0, 0, -.11718794, -.57058805, -.59488791, 0, -.65658802, 0, -.52698797, 0, 0, 0, 0, 0, 1.3500118, 0, 1.665812, 0, 1.963912, 0, 0, 1.9371119, .90991193, -.39558789, .39521196, 0, -.05268808, 0, 0, 0, 0, -.12458798, 0, -.28228804, .79281193, -.26358792, 0, 0, 0, -.72828788, 0, .355912, 0, 0, 0, 0, -.43538806, -1.566388, 0, -.28388807, 0, -.69028801, 0, -.78128809, -.54648799, -.92738789, 0, 0, .61571199, 0, 0, 1.012012, 0, 0, 0, 0, .43991187, .9404121, 0, 0, 0, .61671191, 0, 2.6073117, 0, -.60438794, 0, -.18108793, -.48178813, 0, -.22628804, -.07398792, 0, 0, 0, 1.830512, 0, 0, 0, 0, 0, 0, -.36918804, 1.3247118, 0, 0, 1.163012, 0, 0, 0, 1.3241119, 0, 0, 0, 0, -.90038794, -1.250888, 0, 0, 0, 0, -1.048188, -.90138787, 0, 0, 0, 0, 0, 0, -.00878807, 0, 0, .46301201, -.22048803, -.71518797, 0, 0, -.4952881, 0, -.83718795, .57951194, 0, 0, 0, 1.054112, .61721212, 2.2717118, 0, 2.0280118, 0, 0, 0, 2.6503119, 2.3914118, 0, -.36418793, -.9259879, 0, 0, .16971211, 0, 1.9360118, 2.5344119, 2.0171118, 0, 0, 0, 0, 0, 0, 0, 0, .77211195, 0, 0, 0, 0, 0, 1.952312, 0, 0, .25901201, 0, -.5028879, .03641204, 0, .38571194, 0, 0, -.44528791, -.55918807, 0, .39721206, 0, -.34038803, -.05988808, 0, -.03518792, .045512, 0, -.039288, 0, .01431207, 0, 0, .030412, -.31918809, 0, 0, 0, -.324388, 0, -1.232188, 0, -.23678799, -.89188808, 0, -.6766879, 0, 0, 0, 0, 0, 0, 0, -.67818803, 0, 1.099412, 1.2767119, -.64068788, -.50678796, 0, 0, 0, 0, 0, .68371207, .11251191, -.17128797, .17081194, 0, -.48708794, .09591202, 0, -.20108791, -.02158805, 0, -.3012881, 0, 0, 0, -.87638801, -.54488796, 0, 0, 0, 0, -.14738794, 0, 0, 0, 0, -.75718802, -.37418792, 0, 1.0981121, 1.1441121, .47381189, 0, 0, 0, -.33498809, 0, 0, 0, 0, -.91838807, 0, -.34488794, 0, .12971191, .99381214, -.91608804, 0, 0, .98171192, 0, 0, 0, 0, -.01528808, 0, -.41458794, .25691202, .18601207, 0, 0, 0, -.35608789, .79691201, 0, 0, 1.548912, 0, 0, 0, 0, 0, 0, 0, 0, 1.1663117, 0, 0, -.009088, -.49578807, 0, 0, -.2677879, 0, -.25468799, .68631202, 0, 0, 0, -.36198804])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def load_data_cox(dta):\n",
" array = lambda x : np.array(dta[x], dtype=float)\n",
" t = array('t')\n",
" obs_t = array('obs_t')\n",
" pscenter = array('pscenter')\n",
" hhcenter = array('hhcenter')\n",
" ncomact = array('ncomact')\n",
" rleader = array('rleader')\n",
" dleader = array('dleader')\n",
" inter1 = array('inter1')\n",
" inter2 = array('inter2')\n",
" fail = array('FAIL')\n",
" return (t, obs_t, pscenter, hhcenter, ncomact,\n",
" rleader, dleader, inter1, inter2, fail)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"(t, obs_t, pscenter, hhcenter, ncomact, rleader,\n",
" dleader, inter1, inter2, fail) = load_data_cox(dta)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"X = np.array([pscenter, hhcenter, ncomact, rleader, dleader, inter1, inter2])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(7, 430)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X.shape"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"with Model() as model:\n",
" \n",
" T = len(t) - 1\n",
" nsubj = len(obs_t)\n",
"\n",
" # risk set equals one if obs_t >= t\n",
" Y = np.array([[int(obs >= time) for time in t] for obs in obs_t])\n",
" # counting process. jump = 1 if obs_t \\in [t[j], t[j+1])\n",
" dN = np.array([[Y[i,j]*int(t[j+1] >= obs_t[i])*fail[i] for j in range(T)] for i in\n",
" range(nsubj)])\n",
"\n",
" c = Gamma('c', .0001, .00001)\n",
" r = Gamma('r', .001, .0001)\n",
" \n",
" dL0_star = r*np.diff(t)\n",
" \n",
" mu = dL0_star * c # prior mean hazard\n",
" \n",
" dL0 = Gamma('dL0', mu, c, shape=T)\n",
"\n",
" beta = Normal('beta', np.zeros(7),\n",
" np.ones(7)*.00001, shape=7)\n",
"\n",
" linear_model = tt.exp(tt.dot(X.T, beta))\n",
" idt = Y[:, :-1] * tt.outer(linear_model, dL0)\n",
"\n",
" dn_like = Poisson('dn_like', idt, observed=dN)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" [-----------------100%-----------------] 2000 of 2000 complete in 22.2 sec"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/Cellar/python3/3.4.2_1/Frameworks/Python.framework/Versions/3.4/lib/python3.4/importlib/_bootstrap.py:321: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility\n",
" return f(*args, **kwds)\n"
]
}
],
"source": [
"with model:\n",
" start = find_MAP()\n",
" step = NUTS(scaling=start)\n",
" trace = sample(2000, step)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.gridspec.GridSpec at 0x10b598be0>"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEgCAYAAADhUed1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF1BJREFUeJzt3X2QZXV95/H3R9BZd0d0emLQRXRhjalMa5ZdVxNLCNEu\nH2YgoGSjpqR0Z9woxiiLUTIgOG4mizii8Wmjru6MD5RGSREHXTTFzoyArlHXUrBux0RlRvAxFj0Y\nUDEYvvvHvTO5XLqb7p6eOefcfr+quqbPOb97zrdvX+7nnnO+/SNVhSRJbXO/pguQJGk2BpQkqZUM\nKElSKxlQkqRWMqAkSa1kQEmSWsmA0iFLcr8kr0+yN8n3kuxI8sCRMS9MckeSm5PcMvh37WBbknxw\nsP6zSR42WP/mJK9cZC0PTfK+JN8fHKOX5MLl+2nnPO4Hk7x2jm3vSnLGLOv3JnnqPPt8wOA5+W6S\nu5Mctv9ek1wzdJybk7x9ifvZkuQDy13fckuyc/Dc3jXf70DNOrrpAjQWzgWeDDwO+CnwXuAtwEtG\nxn20qjbN8vhnAMdV1fFJXgecm+T9wBTwHxZZyyeAG4ATq+onSR4EPHKR+1hWVXXOEh/3j8DxSR4F\n3LS8Vd3rWE8bOs6jasz/QLKqzgRIclifVx0az6C0HH4TuKKq7qiqfwLOA56f5BdGxmWOx/8cOGpw\nhvCAwfJbgVcP9rcgSTYAjwDOqaqfAFTV7VXVGxm3Y3DG9/bB2cK3k0wObU+SzUm+PjjLeWeSVSP7\neO5g+82DM4bVs9TzpcGn9DuTzBbMAI9L8ukkP0hyfZJfWujPO3ScTUn+Jsm3knwkycRi9zG6y1mO\nsSfJqwfP3b4k30xy6siYL9H/3T976Cz5X4+MOT3JVwbbPpXkkSPbHzU4izslyZeTfCfJe4e2vzvJ\n20Ye8+EkW4eW1yT508HZ8/eS3JDk5IX+rEP72ZTkq4Nav5Xk1XON1WFSVX75dUhfwEXA9cDawfIT\ngf3AU4bG/O5g3V7gWuDUkX38D+CbwCeBFwMfW0IdrwX+cgHjdgA/AH5vsPygke3nAV8FHkb/DWw7\n8I6h7b9C/0zx1MHyk4BbgdfOcbw9wKZZ1u8d/LwPGSy/CfjSLOMeBfwTcL9Ztj0b+A7wmMHyfwM+\nscTf43zH2QNMA/9+sPzHwHWzjNsCfGCO/T8e+Afg1wfLm4AbZqnhbuBTwLGjvx/gPwJ/Dxw9WD4G\nuAN45NCYVcDTgAw9J705atoLPHWW9Y8e/I5/ZWjd6iPx35Nf//zlGZSWwyXALuC6JH8HnAP8DPjF\nAwOq6sP030RPoP+G8bEkJwxtf1lV/VvgWfQvGZ63hDqOA25f4Ngrquo9g2OPPubFwNaq+n7135le\nA7xoaPvvAB+vqmsHj/8ccPUS6gV4Y1XdNvj+tcC/G1xqW6gXA2+tqr8bLG8FfjPJw5dYz3wuq6ov\nD76/jsVfOv0vwOVV9dcAVbUdWJXk10fGFfD8qvrBYNzB309V/T/6gXz6YNVzgc9U1c1DY35WVdcM\nfncAfwE8ZpG1/gi4E3j6gbPAqrpjkfvQITKgdMiq6u6qel1VTVbVY+i/ER1D/yxleNwPB//upv8G\nd/q9dgavAq4EzknyuSQfTfKQBZbyDwyF4n2YmWfb8cCbktw0uEfxf4GfDL3pP4z+J+9lVVU/pn+W\n+dBFPOx44LyhWr9O/5P/4b7vdheLf/84HvidA7Um2QtMDNaP2j/Pft4L/OfB9y8A3jO8cXCJ9rwk\n1yW5nv7Z+f0W02QyeK0+kf4Z3Z7B5dqphT5ey8OA0uFwBv1Pn389z5j7D8YclOQRwPOBzwNPqKon\n0b+s9F8XeNzPA09I8i8XXfE9fRN4YVWdOPg6oarWVtX3Btu/DZww8pijlnisg41K6Xc1rgFunnv4\nrLW+ZqTWh1bV55dYz+F0E/DOkVp/saquWOR+Lqd/lvgk+pfido5s/wPgt4H/VFWn0A+xRauqr1fV\nK6vql+mf9X989F6kDi8DSocsyb9K8uDB9ycC2+jfj7lzaMwjD3yCTfI04FeBj43s6o303whGGyMe\nyMLsBL4LvDvJMYNjrU0y25nafP4UeEuSXx6q/8FD2z8CbEjya4Ntz6EfykvxhiTHJjma/vP2qar6\n+1nGzXUz/y3AliRPHKp1oWecsx1jzqaBBdpPv/Hj/kmOGrnU+GfAS5I84+ABZ6913hqq6kfAVcCH\ngfdX1c9HhhwPfB/44aBR542D9fdfzA8yfAmafvPOz7j3a1OHkQGl5fBY4IYk3wb+CnhbVb1jZMwf\nAt9Osg84H3j6gUt+AIMuq4dV1UcG+/h+kluAZwJvYwGq6i76rel3A18fPP56+jfn7zH0PvbzPvph\n8eeD7q2bgP8+tP2bwMvp30e7GTiV+e9BzXW8Aj5EP6h/AJxI/57SXGP3JfnKSK17Bo9566Db7Cbg\nffP9fLNJcg3w2aHjjP4d1ELbzj8E/AS4hX6jybOGav0a8FvAHw2e1730z0oesIRjvYd+EP2vWba9\nif5l0u8C/xv4n/SD5bhZxhbw4cFzd7CDcvAB54oD3YjAy4BnzhKGOowOdLlIktQqnkFJklrJgJIk\ntZIBJUlqJQNKktRKK2ay2F27dtkNIkktNDU1NeufFqyYLr79+/evjB9U9zAxMcHMzHyTRkhq0po1\na+b8uzcv8UmSWsmAkiS1kgElSWolA0pj7fzzz2+6BElLZJOEJKkxNklIkjrHgJIktZIBJUlqJQNK\nktRKBpTG2qWXXtp0CZKWyC4+jTWnOpLazS4+SVLn3GdAJfniYnea5PeWVs699rMlyY1Jzhpad3KS\nLyTZNjJ2a5K9STYsx7G1eFVFr9ej1+uxUs7MJR0+CzmDWso7zYuX8Ji5bK6qK4eWVwGXjA6qqouB\nHct4XC1CVbFx4yWsX7+P9ev3sWnT6w0pSYdkIf8/qAcleS/wOOCTVfU6gCTPALYAdwNvPhAiSd4P\nPCbJbmB3Vf3JYP1G4FnAvwHeW1VvX0rBVbUryalzbJ7zWqaWx8TEmnm2Xnbwu507z2bt2vn3NTOz\nf3mKkjSWFhJQ/wI4F/gJcH2SdwPfB94APBm4E9iT5BNV9Y9V9cIkX6iqp47s54NVtSPJKuALwJIC\nSks3f7gcectRz32FnHPxSd21kID6YVX9GCDJl4FHAD8f/Ptx+mctDwaOA/YOHjPbmcxvJDkd+DHw\nwEOsW0twOM9YDlzi2737JACmpm5g+/YLSJo9qd28eXOjx5e0dAsJqOOSrAFuBx4PvLaq9ieZBs6s\nqttnecz9k6TueRPibVX12CSPBJ576KV7Oa9NkrBjx4VMT08DsG7dhsbDSVK3LSSg9gNvBh4LfKCq\nDnwMvwD4eJICvlNVZw895hrgk0n2VtVLB+s+k+SzwFeAxfxhyj3e5ZL8EbAeODbJMVX1kkXsS4dR\nEiYnJ5suQ9KYaPUf6iY5F3ge8MaRTr7Zxm4FTgNeWVWfHt3uH+pKUvvM94e6rQ6o5WRASVL7OJOE\nVizn4pO6yzMojTXn4pPazTMoSVLnGFCSpFYyoCRJrWRASZJayYDSWHMuPqm77OKTJDXGLj5JUucY\nUJKkVjKgJEmtZEBJklrJgNJYcy4+qbvs4tNYcy4+qd3s4pMkdY4BJUlqJQNKktRKBpQkqZUMKI01\n5+KTussuPklSY+zikyR1jgElSWolA0qS1EoGlCSplQwojTXn4pO6yy4+jTXn4pPardNdfEm2JLkx\nyVmD5Xcl2ZPk00lOGBq3NcneJBuaq1ZdVVX0ej16vR4r5UOb1HZHN13AAm2uqqsBquocgCRPAc4H\nXjpYf3GSu5orUV1VVWzceAm7d58EwNTUTrZvv4Bkzg92ko6ArgTUbG4HfjayzncULdr09DRXXXXZ\nweVdu/rrJicnG6xKUusv8c3jRcC7mi5C4+eOO87mlFNObroMacXrZEAlOR3426r6WtO1qN0WMhff\nunXrOOOMV7F69eWsXn05Z575am691cYKqWmt7+JLsgX44oF7UEkeD/xuVb3qvsYOs4tP86kqpqen\ngX5gef9JOjLm6+Lryj2o4R/gCuCWJHuAG6vq3IZq0hhJ4j0nqWW6EFC3ARclWVVVV1bVibMNSrIV\nOA249ohWJ0k6LFp/iW+5eIlPktqn03+oK0lamQwojTXn4pO6y0t8GmvOxSe1m5f4JEmdY0BJklrJ\ngJIktZIBJUlqJQNKY20hc/FJaie7+CRJjbGLT5LUOQaUJKmVDChJUisZUJKkVjKgNNaci0/qLrv4\nNNaci09qN7v4JEmdY0BJklrJgJIktZIBJUlqJQNKY825+KTusotPktQYu/gkSZ1jQEmSWsmAkiS1\nkgElSWolA0pjzbn4pO5qfUAl2ZLkxiRnDZb/JMnuJNckOWFo3NYke5NsaK7ae6sqer0evV6PldIx\n2Sbbtm1rugRJS3R00wUs0Oaquhqgqi4CSPJkYDPwksH6i5Pc1VyJ91ZVbNx4Cbt3nwTA1NROtm+/\ngGTOrkpJ0kBXAmo2vwb8zci6Vr3zT09Pc9VVlx1c3rnzbNauhZmZ/Q1WJUnd0MmASnIt8HDg5KZr\nGTUxsWZoafby7jnGwJKk2XQyoKrq1CRPAD4APLPpeoYNh829L/Hd4CU+SVqgTgbUwA+AVncdJGHH\njguZnp4GYN26DYbTEeZcfFJ3dSWgDr6rJ/kI8AvAT4E/aKyiBUrC5ORk02WsWJs3b266BElL1IWA\nug24KMmqqrqyqp4726AkW4HTgGuPaHWSpMPC2cwlSY1xNnNJUucYUJKkVjKgNNaci0/qLu9BaaxN\nTEwwMzPTdBmS5uA9KElS5xhQkqRWMqAkSa1kQEmSWsmA0lhzLj6pu+zikyQ1xi4+SVLnGFCSpFYy\noCRJrWRASZJayYDSWHMuPqm77OLTWHMuPqnd7OKTJHWOASVJaiUDSpLUSgaUJKmVDCiNNefik7rL\nLj5JUmPs4pMkdY4BJUlqJQNKktRKBpQkqZVaH1BJtiS5MclZQ+sekGRfkt8fWrc1yd4kG5qptF2q\nil6vR6/XY6U0wszGufik7mp9QA1srqorh5bPAb40PKCqLgZ2HNGqWqqq2LjxEtav38f69fvYtOn1\nKzaktm3b1nQJkpbo6KYLWKwkDwSeDnwUWD26+chX1A4TE2tG1lx28LudO89m7dp7bp2Z2X/4i5Kk\nQ9C5gAJeAbwDOLbpQo6Ue4fPkdmnISapSZ0KqCTHAKdU1RuSvJAVcsa02KA4cIlv9+6TAJiauoHt\n2y8gWRFPl6Qx0amAAk4GViX5EHAicFSSPVU13XBdrZKEHTsuZHq6/7SsW7fBcJLUOV0JqABU1dXA\n1QBJXgCsNpxml4TJycmmy2icc/FJ3dX6ufiSnAs8D3jjSCff6LitwGnAK6vq06PbnYtPktpnvrn4\nWh9Qy8WAkqT2cbJYSVLnGFCSpFYyoCRJrWRAaaw5F5/UXTZJaKxNTEwwMzPTdBmS5mCThCSpcwwo\nSVIrGVCSpFYyoCRJrWRAaaw5F5/UXXbxSZIaYxefJKlzDChJUisZUJKkVjKgJEmtZEBprDkXn9Rd\ndvFprDkXn9RudvFJkjrHgJIktZIBJUlqJQNKktRKBpTGmnPxSd1lF58kqTF28UmSOseAkiS1kgEl\nSWql1gdUki1Jbkxy1mB5R5LPJdmd5AVD47Ym2ZtkQ3PV/rOqotfr0ev1WCn3+SRpObU+oAY2V9WV\ng+8LeE5VPbWqPnBgQFVdDOxopLoRVcXGjZewfv0+1q/fx6ZNrzekGuJcfFJ3tb6LL8kW4ItVdfVg\neQfwuqr61n2NHXYkuvgmJtYseOzMzP7DWIkOcC4+qd3m6+I7+kgWskxuBz6U5FbgvKr65pE68GIC\n6FD3ZYBJWuk6F1BV9QqAJCcBlwHPPlLHXmhoHLjEt3v3SQBMTd3A9u0XkMz5QUGSNKJzATXkTuCu\npouYTRJ27LiQ6elpANat22A4SdIidSWgDr67J/lz4OH0L/W9rLGK7kMSJicnmy5DkjqrCwF1G3BR\nklVVdWVVPW+2QUm2AqcB1x7R6tRqzsUndVfru/iWi3PxSVL7OBefJKlzDChJUisZUJKkVjKgJEmt\nZEBprDkXn9RddvFprDkXn9RudvFJkjrHgJIktZIBJUlqJQNKktRKBpTGmnPxSd1lF58kqTF28UmS\nOseAkiS1kgElSWolA0qS1EoGlMaac/FJ3WUXn8aac/FJ7WYXnySpcwwoSVIrGVCSpFYyoCRJrWRA\naaw5F5/UXXbxSZIaYxefJKlzDChJUiu1PqCSbElyY5KzBsvHJdmd5LokbxoatzXJ3iQbmqtWkg6P\nqqLX69Hr9Vgpt2aObrqABdpcVVcPvr8MeE1VfW54QFVdnOSuI1+aJB1eVcXGjZewe/dJAExN7WT7\n9gtI5rx9Mxa6ElAAJLkf8OjRcBoeciTrUftdeumlbN68uekypHuZmFizyEdcdvC7nTvPZu3axT16\nZmb/Io/XvNZ38SXZAnyxqq5Ocizwf4BvAMcA76iqv5xt7Oh+7OJbmZyLT4di8SGyMh1K+M3Xxdep\nMyjgVuA24Lfp1/7ZJJ+qqp82W5akcdSWs457X+K7wUt8bVNVP09yC/DwqvpOkjubrkmSDrck7Nhx\nIdPT0wCsW7dh7MMJuhNQw7+JzcB7khwDXOHZk6SVIAmTk5NNl3FEdSGgbgMuSrKqqq6sqpuBe7WS\nJ9kKnAZce6QLlCQtv9Y3SSwXmyRWJrv4pHabr0nCgJIkNca5+CRJnWNASZJayYCSJLWSASVJaiUD\nSmPt0ksvbboESUtkF5/GmnPxSe1mF58kqXMMKElSKxlQkqRWMqAkSa1kQGmsnX/++U2XIGmJ7OKT\nJDXGLj5JUucYUJKkVjKgJEmtZEBJklrJgNJYcy4+qbvs4tNYcy4+qd3s4pMkdY4BJUlqJQNKktRK\nBpQkqZUMKI015+KTussuPklSY+zikyR1TusDKsmWJDcmOSvJMUn2JNk9+Pe2oXFbk+xNsqHJerU0\nVUWv16PX67FSzuolze/opgtYoM1VdfXg+6cAJPlV4OUHBlTVxUnuaqI4HZqqYuPGS9i9+yQApqZ2\nsn37BSRznvlLWgG6ElCzeTnw9pF1vqN10PT0NFddddnB5Z07z2Z6+jNMTk42WJWkprX+Et9skkwA\nx1fVjU3XokN3yiknz7puYmLNIe/bufik7upkQAEvBt7TdBFaHrfeOsMZZ7yK1asvZ/XqyznzzFdz\n660zzMzsP+R9b9u2bRkqlNSEzl3iS3IU8FvAKU3XouWRhB07LmR6ehqAdes2eP9JUmcCavjd6tnA\nx6vq7qaK0fJL4j0nSffQhYC6DbgoyaqqurKq/mK2QUm2AqcB1x7R6iRJh4UzSWis+f+DktrNmSS0\nYjkXn9RdnkFJkhrjGZQkqXMMKElSKxlQkqRWMqAkSa1kQGmsORef1F128Wms+XdQUrvZxSdJ6hwD\nSpLUSgaUJKmVDChJUiutmCaJXbt2rYwfVJI6ZmpqatZGiRUTUJKkbvESnySplQwoSVIrGVCSpFYy\noCRJrWRAdUSSdyXZk+TTSU5oup4uSDKV5Pok1yV5atP1dIWvtUOX5AFJ9iX5/aZr6TK7+DomyVOA\n51TVS5uupc2SBPgMMAUE+Kuq+o1mq+oWX2tLl+QVwKnArqr6s6br6SrPoLrnduBnTRfRAb8E/G1V\n3VlVPwW+keTRTRfVMb7WliDJA4GnAzubrqXrjm66AN1TkqcB5wNF/5N/AX9YVV8dDHkR8NaGyuuS\ntcCPkryZ/vP4o8G6bzRaVbf4WluaVwDvAI5tupCuM6BapqquAa6ZbVuS0+mfFXztyFbVSbcCDwFe\nSj+g3jlYpwXwtbY0SY4BTqmqNyR5If3XnpbIgOqIJI8HfrOqXtV0LR3xDfqX+aD/JvHoqvLsaQF8\nrR2Sk4FVST4EnAgclWRPVU03XFcn2STREUluAm4B7gZurKpzGy6p9QaXS7fQv0z6x4OzU90HX2vL\nI8kLgNU2SSydASVJaiW7+CRJrWRASZJayYCSJLWSASVJaiUDSpLUSgaUJKmVDChJUisZUJKkVvr/\n/7lwZR/e+f0AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10b598940>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"forestplot(trace, vars=['beta'])"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.4.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment