Skip to content

Instantly share code, notes, and snippets.

@jseabold

jseabold/pymc_spatial_surival_debug.ipynb

Created June 25, 2014 16:47
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 jseabold/b8e39e7751138675fe58 to your computer and use it in GitHub Desktop.
Save jseabold/b8e39e7751138675fe58 to your computer and use it in GitHub Desktop.
Download ZIP
Replicates Table 4 columns 3 and 5. Fails on Table 3 column 3 of Darmofal's "Bayesian Spatial Survival Models for Political Event Processes"
Display the source blob
Display the rendered blob
Raw
pymc_spatial_surival_debug.ipynb
{
"metadata": {
"name": "",
"signature": "sha256:524c06573310f090451fce1ad4e525e4332365e2c618124152bf7f94503d7bd1"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Background"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Model employed in Darmofal's \"Bayesian Spatial Survival Models for Political Event Processes\"\n",
"\n",
"Adapted from the WinBUGS code found [here](http://people.cas.sc.edu/darmofal/AJPSCodeandData.zip)\n",
"\n",
"CAR prior for PyMC is introduced [here](http://glau.ca/?p=340)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"import pandas as pd\n",
"from pymc import (Weibull, Gamma, Normal, Lambda, \n",
" deterministic, MCMC, stochastic, normal_like, Poisson,\n",
" AdaptiveMetropolis, MAP, Matplot)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### All Data"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"indiv_dta = dict(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, 396, 232, 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",
" pscenter = [ -.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, .14392118, .14392118, .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, .95411211, .65311199, -.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, 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, 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, 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, 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",
" 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, 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],\n",
" # adjacency matrix\n",
" adj = [[2, 7, 83, 229], [1, 3, 7, 83, 84, 107], [2, 4, 6, 7, 107, 108, 112], [3, 5, 6, 7, 108, 112, 114, 225, 227], [4, 114, 225, 360, 361], [3, 4, 7], [1, 2, 3, 4, 6, 227, 229], [9, 10, 11, 222, 223, 225, 226, 365, 366], [8, 10, 11], [8, 9, 11, 222, 318, 319], [8, 9, 10, 169, 170, 226, 319, 367], [13, 15, 17], [12, 14, 15, 16, 17, 69], [13, 15, 17, 57, 61, 69, 266, 397, 399], [12, 13, 14, 17], [13, 17, 263], [12, 13, 14, 15, 16, 72, 263, 264, 399], [19, 20, 23, 24, 324, 326], [18, 20, 21, 266, 324], [18, 19, 21, 22, 24, 28], [19, 20, 28, 35, 36, 57, 266], [20, 28], [18, 24], [18, 20, 23, 26, 27, 28], [29], [24, 27, 30], [24, 26, 28, 30, 33, 35], [20, 21, 22, 24, 27, 35, 48], [25, 31], [26, 27, 31, 32, 33], [29, 30, 32], [30, 31, 33, 34], [27, 30, 32, 34, 35], [32, 33, 35, 37, 39], [21, 27, 28, 33, 34, 36, 37], [21, 35, 37, 38, 57], [34, 35, 36, 38, 39], [36, 37, 39, 40, 42, 57], [34, 37, 38, 40], [38, 39, 41, 42], [40, 42, 43, 46], [38, 40, 41, 43, 44, 45, 57], [41, 42, 44, 46], [42, 43, 45, 46, 47, 48], [42, 44, 48, 51, 57, 58], [41, 43, 44, 47, 49, 53], [44, 46, 48, 49, 50, 62], [28, 44, 45, 47, 50, 51], [46, 47, 50, 52, 53], [47, 48, 49, 51, 52, 54, 55], [45, 48, 50, 55, 56, 58], [49, 50, 53, 54], [46, 49, 52, 54, 55], [50, 52, 53, 55], [50, 51, 53, 54, 56, 62], [51, 55, 58, 62, 63, 64], [14, 21, 36, 38, 42, 45, 58, 59, 60, 61, 266], [45, 51, 56, 57, 59, 60, 64], [57, 58, 60, 61], [57, 58, 59, 61, 64, 65], [14, 57, 59, 60, 65, 69], [47, 55, 56, 63, 64], [56, 62, 64], [56, 58, 60, 62, 63, 65], [60, 61, 64, 68, 69], [67, 68, 69], [66, 69], [65, 66, 69], [13, 14, 61, 65, 66, 67, 68], [71, 73, 75], [70, 72, 73, 75], [17, 71, 73, 74, 75, 264, 399, 432], [70, 71, 72, 74, 75, 156, 246, 264, 322, 432], [72, 73, 75], [70, 71, 72, 73, 74], [77, 80, 81], [76, 78, 80, 174, 349, 350], [77, 79, 80], [78, 80, 284, 285], [76, 77, 78, 79, 81, 285], [76, 80, 173, 174, 285, 288], [183, 249, 250, 334, 343], [1, 2, 84], [2, 83, 85, 86, 87, 88, 106, 107, 113], [84, 86, 87, 88, 89, 90], [84, 85, 88, 89, 106], [84, 85, 88, 91, 94], [84, 85, 86, 87, 90, 94, 97], [85, 86, 90, 97], [85, 88, 89, 97], [87, 92, 93, 94], [91, 93], [91, 92, 94, 95], [87, 88, 91, 93, 95, 96, 97, 98], [93, 94, 96], [94, 95, 98, 102], [88, 89, 90, 94, 98], [94, 96, 97, 101, 102, 104, 105], [100, 102, 103, 104, 105], [99, 102, 103, 104, 105], [98, 102, 105], [96, 98, 99, 100, 101, 102, 103, 105], [99, 100, 102], [98, 99, 100, 105], [98, 99, 100, 101, 102, 104], [84, 86, 113, 116, 352], [2, 3, 84, 108, 113], [3, 4, 107, 110, 112, 113, 116], [110, 111, 115, 116], [108, 109, 111, 112, 116], [109, 110, 112, 114, 115], [3, 4, 108, 110, 111, 114], [84, 106, 107, 108, 116], [4, 5, 111, 112, 115, 241, 353, 360], [109, 111, 114, 116, 353], [106, 108, 109, 110, 113, 115, 352, 353], [118], [117], [120, 121, 142], [119, 121, 123, 141, 142, 208, 422], [119, 120, 122, 123, 142, 221, 224], [121, 123, 221, 244, 245], [120, 121, 122, 208, 209, 244, 357], [125, 266, 230, 324], [124, 230, 266, 397, 432], [127, 128, 129, 132], [126, 128, 136, 138], [126, 127, 130, 131, 136, 138], [126, 130, 131, 132], [128, 129, 131, 132, 134], [128, 129, 130, 132, 133, 134, 135, 138, 139], [126, 129, 130, 131], [131, 135, 139, 141, 420], [130, 131, 135], [131, 133, 134, 420], [127, 128, 138, 139, 140, 142, 143, 146, 150], [144, 145, 160, 216, 217, 218, 223], [127, 128, 131, 136, 139], [131, 133, 136, 138, 141, 142], [136, 143, 144, 150, 152], [120, 133, 139, 142, 420, 421, 422], [119, 120, 121, 136, 139, 141, 143, 145, 224], [136, 140, 142, 144, 145], [137, 140, 143, 145, 152, 153, 160, 223], [137, 142, 143, 144, 217, 224], [136, 148, 150], [149, 150, 151, 154, 305], [146, 149, 150, 197], [147, 148, 150, 197, 198, 302, 305], [136, 140, 146, 147, 148, 149, 151, 152], [147, 150, 152, 154, 155], [140, 144, 150, 151, 153, 154, 155], [144, 152, 154, 160, 161], [147, 151, 152, 153, 161, 162, 163, 164, 298, 299, 305], [151, 152], [73, 157, 159, 244, 246, 322], [156, 158, 159, 219, 221, 222, 244, 318], [157, 219, 220, 221], [156, 157, 318, 321, 322], [137, 144, 153, 161, 164, 223, 361, 362, 363, 364, 365], [153, 154, 160, 162, 163, 164, 165], [154, 161, 163], [154, 161, 162, 164, 165, 298, 299, 303, 431], [154, 160, 161, 163, 165, 361, 408, 431], [161, 163, 164], [167, 168, 169, 228, 229], [166, 168], [166, 167, 169, 172], [11, 166, 168, 169, 170, 171, 172, 226, 228], [11, 169, 171, 367, 368], [169, 170, 172, 368], [168, 169, 171, 368, 375], [81, 174, 175, 177, 248, 288, 411], [77, 81, 173, 175, 349], [173, 174, 176, 177, 349], [175, 177, 179, 180, 181, 182], [173, 175, 176, 177, 178, 179, 248], [177, 179, 247, 248], [176, 177, 178, 180], [176, 179, 181], [176, 180, 182], [176, 181], [82, 184, 185, 187, 343, 400], [183, 185, 188, 189, 343, 346], [183, 184, 187, 188, 189], [187, 190, 407], [183, 185, 186, 188, 190, 400, 407, 410], [184, 185, 187, 189, 190, 336, 339, 346, 409, 429, 430], [184, 185, 188], [186, 187, 188, 407, 409], [192, 247], [191, 247, 248], [194, 195, 196, 426, 427], [193, 195, 197, 198], [193, 194, 196, 198, 199], [193, 195, 199, 200, 201], [148, 149, 194, 198], [149, 194, 195, 197, 199, 204, 207, 302, 306], [195, 196, 198, 200, 202, 204], [196, 199, 201, 202, 203], [196, 200, 203, 205], [199, 200, 203, 204, 205, 207], [200, 201, 202, 205], [198, 199, 202, 207], [201, 202, 203, 206, 207], [205, 207], [198, 202, 204, 205, 206, 306], [120, 123, 209, 210, 422], [123, 208, 210, 213, 214, 215, 357], [208, 209, 211, 212, 213, 422], [210, 212, 213], [210, 211, 213], [209, 210, 211, 212, 215, 422], [209, 215, 243, 357], [209, 213, 214, 422, 426], [137, 217, 218], [137, 145, 216, 218, 224], [137, 216, 217, 223, 224], [157, 158, 220, 221, 222, 223, 224], [158, 219, 221], [121, 122, 157, 158, 219, 220, 224, 244], [8, 10, 157, 219, 223, 318], [8, 137, 144, 160, 218, 219, 222, 224, 365], [121, 142, 145, 217, 218, 219, 221, 223], [4, 5, 8, 226, 227, 361, 362, 364, 366], [8, 11, 169, 225, 227, 228], [4, 7, 225, 226, 228, 229], [166, 169, 226, 227, 229], [1, 7, 166, 227, 228], [124, 125, 243, 357, 432], [232, 233, 235, 237, 403, 404], [231, 233, 234, 235, 237, 238, 242], [231, 232, 237, 401, 403], [232, 236, 238, 242], [231, 232, 236, 240, 242, 358, 404, 408], [234, 235, 238, 240], [231, 232, 233, 238, 351, 355], [232, 234, 236, 237, 239, 242, 355], [238, 240, 241, 242, 355], [235, 236, 239, 241, 242, 358], [114, 239, 240, 353, 354, 355, 358, 359, 360], [232, 234, 235, 238, 239, 240], [214, 230, 357], [122, 123, 156, 157, 221, 245, 246, 357], [122, 244], [73, 156, 244, 357, 432], [178, 191, 192, 248], [173, 177, 178, 192, 247, 411], [82, 250, 251, 328, 330, 334], [82, 249, 251], [249, 250, 252, 330, 335], [251, 254, 260, 335], [256, 257, 259, 260, 283, 284, 286, 335, 337, 342], [252, 255, 260, 261, 279], [254, 255, 256, 258, 259, 260, 261], [253, 255, 257, 258, 259, 261], [253, 256, 261, 274, 281, 283], [255, 256, 261], [253, 255, 256, 260], [252, 253, 254, 255, 259, 335], [254, 255, 256, 257, 258, 274, 279], [263, 264], [16, 17, 262, 264, 382, 385, 389], [17, 72, 73, 262, 263, 322, 385, 399], [266], [14, 19, 21, 57, 124, 125, 265, 324, 397], [268, 271], [267, 269, 271], [268, 270, 271], [269, 271, 272, 275], [267, 268, 269, 270, 272, 284], [270, 271, 275, 284], [275, 278, 280, 281, 282, 283, 284], [257, 261, 275, 276, 278, 279, 280, 281], [270, 272, 273, 274, 276, 277, 278, 284], [274, 275, 277, 278], [275, 276], [273, 274, 275, 276, 280, 284], [254, 261, 274], [273, 274, 278, 281], [257, 273, 274, 280, 282, 283], [273, 281, 283], [253, 257, 273, 281, 282, 284], [79, 253, 271, 272, 273, 275, 278, 283, 285, 286], [79, 80, 81, 284, 286, 288, 292], [253, 284, 285, 292, 337], [288, 289, 290], [81, 173, 285, 287, 289, 290, 292, 411], [287, 288, 290, 291, 292], [287, 288, 289, 291, 411], [289, 290, 292, 293, 297], [285, 286, 288, 289, 291, 297, 337], [291, 294, 295, 296, 297], [293, 295], [293, 294, 296], [293, 295, 297], [291, 292, 293, 296, 332, 337, 348], [154, 163, 299], [154, 163, 298, 300, 303, 304, 305], [299, 304, 305], [302, 304, 305, 309, 313, 315], [149, 198, 301, 305, 306, 310, 313], [163, 299, 304, 315, 429, 430, 431], [299, 300, 301, 303, 305, 309, 312, 315], [147, 149, 154, 299, 300, 301, 302, 304], [198, 207, 302], [308, 310, 316], [307, 316], [301, 304, 312, 315], [302, 307, 311, 313, 314, 316], [310, 313], [304, 309], [301, 302, 310, 311, 314, 315], [310, 313, 315, 316, 331, 348, 429], [301, 303, 304, 309, 313, 314, 429], [307, 308, 310, 314, 348], [318, 321], [10, 157, 159, 222, 317, 319, 321], [10, 11, 318, 320, 321, 367, 370, 379], [319, 321, 322, 379], [159, 317, 318, 319, 320, 322], [73, 156, 159, 264, 320, 321, 379, 385], [325, 327, 414], [18, 19, 124, 266, 325, 326, 327, 414, 415], [323, 324, 327, 414], [18, 324, 327], [323, 325, 326, 324], [249, 329, 330, 334, 340], [328, 334, 340], [249, 251, 328, 335, 340], [314, 339, 341, 345, 347, 348, 429], [297, 333, 336, 337, 339, 348], [332, 336, 337, 338, 340, 342, 343, 344], [82, 249, 328, 329, 340, 343], [251, 252, 253, 260, 330, 340, 342], [188, 332, 333, 339, 344, 346], [253, 286, 292, 297, 332, 333, 338, 342], [333, 337, 342], [188, 331, 332, 336, 347, 348, 429], [328, 329, 330, 333, 334, 335, 342, 343], [331, 345, 347], [253, 333, 335, 337, 338, 340], [82, 183, 184, 333, 334, 340, 344, 346], [333, 336, 343, 346], [331, 341, 347], [184, 188, 336, 343, 344], [331, 339, 341, 345, 429], [297, 314, 316, 331, 332, 339], [77, 174, 175, 350], [77, 349], [237, 355, 356], [106, 116, 353, 355, 356], [114, 115, 116, 241, 352, 354, 355], [241, 353, 355], [237, 238, 239, 241, 351, 352, 353, 354, 356], [351, 352, 355], [123, 209, 214, 230, 243, 244, 246, 432], [235, 240, 241, 359, 361, 408], [241, 358, 360, 361], [5, 114, 241, 359, 361], [5, 160, 164, 225, 358, 359, 360, 363, 408], [160, 225, 363, 364], [160, 361, 362, 364], [160, 225, 362, 363, 365, 366], [8, 160, 223, 364, 366], [8, 225, 364, 365], [11, 170, 319, 368, 370], [170, 171, 172, 367, 370, 371, 374, 375], [370, 371, 392, 396], [319, 367, 368, 369, 371, 379, 392], [368, 369, 370, 374, 377, 380, 390, 396], [378, 383, 390, 392, 396], [374, 380, 384, 388, 395], [368, 371, 373, 375, 380, 384], [172, 368, 374, 380, 384, 388, 391], [380, 387], [371, 378, 380, 383, 387, 390], [372, 377, 383, 390], [319, 320, 322, 370, 383, 385, 392], [371, 373, 374, 375, 376, 377, 381, 387, 388, 394], [380, 393, 394], [263, 389], [372, 377, 378, 379, 385, 387, 392], [373, 374, 375, 388, 391, 395], [263, 264, 322, 379, 383, 387, 389], [387, 389, 394], [376, 377, 380, 383, 385, 386, 389, 394], [373, 375, 380, 384, 391], [263, 382, 385, 386, 387, 394], [371, 372, 377, 378, 392, 396], [375, 384, 388, 395], [369, 370, 372, 379, 383, 390, 396], [381], [380, 381, 386, 387, 389], [373, 384, 391], [369, 371, 372, 390, 392], [14, 125, 266, 398, 399, 432], [397, 399], [14, 17, 72, 264, 397, 398, 432], [183, 187, 402, 406, 409, 410], [233, 403], [400, 403, 406], [231, 233, 401, 402, 404, 406], [231, 235, 403, 405, 406, 408], [404, 406, 408, 409, 430, 431], [400, 402, 403, 404, 405, 409], [186, 187, 190, 409, 410], [164, 235, 358, 361, 404, 405, 431], [188, 190, 400, 405, 406, 407, 410, 430], [187, 400, 407, 409], [173, 248, 288, 290], [413, 416, 417, 418], [412, 415, 418], [323, 324, 325, 415, 416, 418, 419], [324, 413, 414, 418], [412, 414, 417, 419], [412, 416, 418, 419], [412, 413, 414, 415, 417, 419], [414, 416, 417, 418], [133, 135, 141, 421, 423, 428], [141, 420, 422, 425, 428], [120, 141, 208, 210, 213, 215, 421, 425, 426], [420, 424, 428], [423, 428], [421, 422, 426, 427, 428], [193, 215, 422, 425, 427], [193, 425, 426], [420, 421, 423, 424, 425], [188, 303, 314, 315, 331, 339, 347, 430], [188, 303, 405, 409, 429, 431], [163, 164, 303, 405, 408, 430], [72, 73, 125, 230, 246, 357, 397, 399]])"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"shared_dta = dict(T=73, Nsubj=430, eps=0.0, t=[1, 21, 85, 128, 129, 148, 178, 204, 206, 210, 211, 212, 225, 238, 241, 248, 259, 273, 275, 281, 286, 289, 301, 302, 303, 304, 313, 317, 323, 344, 345, 349, 350, 351, 355, 356, 359, 364, 385, 386, 389, 390, 391, 392, 394, 395, 396, 397, 398, 399, 400, 406, 415, 416, 426, 427, 434, 435, 437, 441, 447, 448, 449, 450, 451, 453, 455, 456, 458, 459, 460, 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= [ -.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",
" 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],\n",
" adj = [[11, 8], [8, 6, 5, 4], [12], [10, 5, 2], [15, 13, 10, 6, 4, 2], [30, 23, 15, 8, 5, 2], [48, 26, 13, 10], [23, 22, 11, 6, 2, 1], [16, 14, 12], [13, 7, 5, 4], [24, 22, 8, 1], [14, 9, 3], [34, 26, 15, 10, 7, 5], [19, 18, 16, 12, 9], [34, 32, 30, 13, 6, 5], [20, 18, 17, 14, 9], [29, 28, 27, 25, 20, 16], [19, 16, 14], [18, 14], [27, 17, 16], [48, 31, 26, 25], [35, 24, 23, 11, 8], [40, 35, 30, 22, 8, 6], [35, 22, 11], [48, 31, 28, 21, 17], [34, 31, 21, 13, 7], [29, 20, 17], [33, 31, 29, 25, 17], [33, 28, 27, 17], [40, 36, 35, 32, 23, 15, 6], [38, 34, 33, 28, 26, 25, 21], [36, 34, 30, 15], [38, 37, 31, 29, 28], [45, 38, 36, 32, 31, 26, 15, 13], [40, 30, 24, 23, 22], [45, 40, 39, 34, 32, 30], [44, 43, 38, 33], [45, 43, 42, 41, 37, 34, 33, 31], [46, 45, 40, 36], [39, 36, 35, 30, 23], [47, 43, 42, 38], [46, 45, 41, 38], [47, 44, 41, 38, 37], [43, 37], [46, 42, 39, 38, 36, 34], [45, 42, 39], [43, 41], [25, 21, 7]]\n",
" )"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def adjust_adjacency(adj):\n",
"\n",
" adj2 = []\n",
" for row in adj:\n",
" adj2.append((np.array(row) - 1).tolist())\n",
" \n",
" adj = adj2\n",
" del adj2\n",
" return adj"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def load_data(individ=True):\n",
" if individ:\n",
" array = lambda x : np.array(indiv_dta[x], dtype=float)\n",
" adj = indiv_dta['adj']\n",
" else:\n",
" array = lambda x : np.array(shared_dta[x], dtype=float)\n",
" adj = shared_dta['adj']\n",
"\n",
" # make zero-based indexing\n",
" adj = adjust_adjacency(adj)\n",
" \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",
" return (obs_t, pscenter, hhcenter, ncomact, \n",
" rleader, dleader, inter1, inter2,\n",
" adj)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def load_data_cox():\n",
" array = lambda x : np.array(shared_dta[x], dtype=float)\n",
" adj = shared_dta['adj']\n",
"\n",
" # make zero-based indexing\n",
" adj = adjust_adjacency(adj)\n",
" \n",
" t = array('t')\n",
" fail = array('FAIL')\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",
" return (obs_t, t, pscenter, hhcenter, ncomact, \n",
" rleader, dleader, inter1, inter2, fail,\n",
" adj)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def weibull_spatialindivid_model():\n",
" (obs_t, pscenter, hhcenter, ncomact, rleader,\n",
" dleader, inter1, inter2, adj) = load_data()\n",
" \n",
" nsubj = 432\n",
" beta0 = Normal('beta0', 0, .001, value=1)\n",
" beta = Normal('beta', np.zeros(7), np.ones(7)*.00001, \n",
" value=np.array([-.36, -.26, -.29, -.22, -.61, -9.73, -.23]))\n",
" rho = Gamma('rho', alpha=.01, beta=100, value=.1)\n",
" tau = Gamma('tau', alpha=.01, beta=.01, value=.001)\n",
" sigma = 1./(tau**.5)\n",
" theta = 1./tau\n",
" #beta = (1./mu)**(1./rho) to agree with BUGS\n",
" \n",
" # neighbor weights for each district\n",
" ww = np.zeros((432, 432))\n",
" for i,row in enumerate(adj):\n",
" ww[i,row] = 1.\n",
" # make row-stochastic / row sums to 1\n",
" ww /= ww.sum(1)[:,None]\n",
" \n",
" @stochastic\n",
" def car(t=tau, value=np.zeros(nsubj)):\n",
" # the conditional mean is the average of the neighbors random effects\n",
" mu_ = np.dot(ww, value)\n",
" # scale precision for the number of neighbors\n",
" taux = t*(ww > 0).sum(1)\n",
" return normal_like(value, mu_, taux)\n",
" \n",
" # zero-centered W\n",
" \n",
" W = Lambda('W', lambda mu=car : mu - np.mean(mu))\n",
" X = np.column_stack((pscenter, hhcenter, ncomact, rleader, dleader, inter1, inter2))\n",
" @deterministic\n",
" def fit_beta(b0=beta0, b1=beta, r=rho, c=car):\n",
" mu = np.exp(b0 + np.dot(X, b1) + W.value)\n",
" return (1./mu)**(1./r)\n",
" \n",
" y = Weibull('obs_t', rho, beta=fit_beta, value=obs_t, observed=True) \n",
" \n",
" return locals()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"variables = weibull_spatialindivid_model()\n",
"m1 = MCMC(variables)\n",
"m1.use_step_method(AdaptiveMetropolis, m1.beta)\n",
"m1.sample(100000, burn=25000)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plt.plot(m1.beta.trace())"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"m1.beta.summary()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"m1.rho.summary()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\theta$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"np.mean(1./m1.tau.trace().mean())"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def weibull_spatialshared_model():\n",
" (obs_t, pscenter, hhcenter, ncomact, rleader,\n",
" dleader, inter1, inter2, adj) = load_data(False)\n",
" \n",
" # the state each observation belongs to\n",
" stfips = (np.array(\n",
" [41, 41, 41, 41, 41, 41, 41, 45, 45, 45, 45, 35, 35, 35, 35, \n",
" 35, 35, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 30, 30, 30, 30, 30, 30, \n",
" 18, 18, 18, 18, 18, 18, 27, 47, 47, 47, 47, 47, 47, 47, 47, \n",
" 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, \n",
" 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 13, 13, 13, 13, \n",
" 13, 8, 8, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26,\n",
" 26, 26, 26, 26, 26, 26, 26, 21, 21, 21, 21, 21, 21, 21, 21, \n",
" 21, 21, 32, 32, 32, 32, 31, 31, 31, 31, 31, 31, 46, 46, 46, \n",
" 46, 46, 46, 46, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 29, \n",
" 29, 29, 29, 29, 29, 29, 29, 3, 3 , 48, 48, 48, 48, 48, 48, \n",
" 48, 48, 48, 48, 48, 48, 48, 48, 48, 10, 10, 10, 10, 10, 10, \n",
" 10, 10, 34, 34, 34, 34, 34, 34, 34, 34, 34, 42, 42, 42, 42, \n",
" 42, 2, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 4, 15, \n",
" 15, 15, 12, 12, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, \n",
" 20, 20, 40, 40, 40, 22, 22, 16, 16, 16, 16, 16, 16, 16, 16, \n",
" 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, \n",
" 16, 16, 16, 16, 16, 16, 16, 16, 25, 25, 25, 25, 25, 25, 25, \n",
" 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 36, 36, 36, \n",
" 36, 36, 36, 11, 11, 11, 11, 11, 17, 17, 17, 17, 17, 17, 17, \n",
" 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 19, \n",
" 19, 44, 44, 44, 44, 44, 44, 5, 38, 38, 38, 38, 38, 38, 38, \n",
" 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, \n",
" 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, \n",
" 39, 39, 23, 23, 23, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, \n",
" 33, 9, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, \n",
" 28, 28, 28, 6]) - 1).tolist() # zero-based indexing\n",
"\n",
" # neighbor weights for each state\n",
" ww = np.zeros((48, 48))\n",
" for i, row in enumerate(adj):\n",
" ww[i, row] = 1.\n",
" # make row-stochastic / row sums to 1\n",
" ww /= ww.sum(1)[:,None]\n",
"\n",
" \n",
" nsubj = 430\n",
" beta0 = Normal('beta0', 0, .001, value=1)\n",
" beta = Normal('beta', np.zeros(7), np.ones(7)*.00001, \n",
" value=np.array([-.36, -.26, -.29, -.22, -.61, -9.73, -.23]))\n",
" rho = Gamma('rho', alpha=.01, beta=100, value=.1)\n",
" tau = Gamma('tau', alpha=.01, beta=.01, value=.001)\n",
" sigma = 1./(tau**.5)\n",
" theta = 1./tau\n",
" #beta = (1./mu)**(1./rho) to agree with BUGS\n",
"\n",
" @stochastic\n",
" def car(t=tau, value=np.zeros(48)):\n",
" # the conditional mean is the average of the neighbors random effects\n",
" mu_ = np.dot(ww, value)\n",
" # scale precision for the number of neighbors\n",
" taux = t*(ww > 0).sum(1)\n",
" return normal_like(value, mu_, taux)\n",
" \n",
" # zero-centered W\n",
" \n",
" W = Lambda('W', lambda mu=car : mu - np.mean(mu))\n",
" X = np.column_stack((pscenter, hhcenter, ncomact, rleader, dleader, inter1, inter2))\n",
" @deterministic\n",
" def fit_beta(b0=beta0, b1=beta, r=rho, c=car):\n",
" mu = np.exp(b0 + np.dot(X, b1) + W.value[stfips])\n",
" return (1./mu)**(1./r)\n",
" \n",
" y = Weibull('obs_t', rho, beta=fit_beta, value=obs_t, observed=True) \n",
" \n",
" return locals()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"variables = weibull_spatialshared_model()\n",
"m2 = MCMC(variables)\n",
"m2.use_step_method(AdaptiveMetropolis, m2.beta)\n",
"m2.sample(100000, burn=25000)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"m2.beta.summary()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"m1.rho.summary()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\theta$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"(1/m2.tau.trace().mean())"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def cox_spatialshared_model():\n",
" (obs_t, t, pscenter, hhcenter, ncomact, rleader,\n",
" dleader, inter1, inter2, fail, adj) = load_data_cox()\n",
" \n",
" T = len(t) - 1\n",
" nsubj = len(obs_t)\n",
"\n",
" stfips = (np.array(\n",
" [41, 41, 41, 41, 41, 41, 41, 45, 45, 45, 45, 35, 35, 35, 35, \n",
" 35, 35, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, \n",
" 24, 24, 24, 24, 24, 24, 24, 24, 24, 30, 30, 30, 30, 30, 30, \n",
" 18, 18, 18, 18, 18, 18, 27, 47, 47, 47, 47, 47, 47, 47, 47, \n",
" 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, 47, \n",
" 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 13, 13, 13, 13, \n",
" 13, 8, 8, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26,\n",
" 26, 26, 26, 26, 26, 26, 26, 21, 21, 21, 21, 21, 21, 21, 21, \n",
" 21, 21, 32, 32, 32, 32, 31, 31, 31, 31, 31, 31, 46, 46, 46, \n",
" 46, 46, 46, 46, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 29, \n",
" 29, 29, 29, 29, 29, 29, 29, 3, 3 , 48, 48, 48, 48, 48, 48, \n",
" 48, 48, 48, 48, 48, 48, 48, 48, 48, 10, 10, 10, 10, 10, 10, \n",
" 10, 10, 34, 34, 34, 34, 34, 34, 34, 34, 34, 42, 42, 42, 42, \n",
" 42, 2, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 37, 4, 15, \n",
" 15, 15, 12, 12, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, 20, \n",
" 20, 20, 40, 40, 40, 22, 22, 16, 16, 16, 16, 16, 16, 16, 16, \n",
" 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, \n",
" 16, 16, 16, 16, 16, 16, 16, 16, 25, 25, 25, 25, 25, 25, 25, \n",
" 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 36, 36, 36, \n",
" 36, 36, 36, 11, 11, 11, 11, 11, 17, 17, 17, 17, 17, 17, 17, \n",
" 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 19, \n",
" 19, 44, 44, 44, 44, 44, 44, 5, 38, 38, 38, 38, 38, 38, 38, \n",
" 38, 38, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, \n",
" 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, 39, \n",
" 39, 39, 23, 23, 23, 33, 33, 33, 33, 33, 33, 33, 33, 33, 33, \n",
" 33, 9, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, \n",
" 28, 28, 28, 6]) - 1).tolist() # zero-based indexing\n",
" \n",
" # neighbor weights for each state\n",
" ww = np.zeros((48, 48))\n",
" for i, row in enumerate(adj):\n",
" ww[i, row] = 1.\n",
" # make row-stochastic / row sums to 1\n",
" ww /= ww.sum(1)[:,None]\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]*(t[j + 1] >= obs_t[i]) * fail[i] \n",
" for i in range(nsubj)] for j in range(T)])\n",
" \n",
" c = Gamma('c', .0001, .00001, value=.01)\n",
" r = Gamma('r', .001, .0001, value=.01)\n",
" dL0_star = r*np.array([t[j+1] - t[j] for j in range(len(t) - 1)])\n",
" mu = dL0_star * c # prior mean hazard\n",
" dL0 = Gamma('dL0', mu, c, value=np.ones(len(t)-1))\n",
" \n",
" beta = Normal('beta', np.zeros(7), np.ones(7)*.00001, \n",
" value=np.array([-.36, -.26, -.29, -.22, -.61, -9.73, -.23]))\n",
" \n",
" tau = Gamma('tau', alpha=.01, beta=.01, value=.001)\n",
" sigma = 1./(tau**.5)\n",
" theta = 1./tau\n",
" \n",
" @stochastic\n",
" def car(t=tau, value=np.zeros(48)):\n",
" # the conditional mean is the average of the neighbors random effects\n",
" mu_ = np.dot(ww, value)\n",
" # scale precision for the number of neighbors\n",
" taux = t*(ww > 0).sum(1)\n",
" return normal_like(value, mu_, taux)\n",
" \n",
" # zero-centered W\n",
" W = Lambda('W', lambda mu=car : mu - np.mean(mu))\n",
"\n",
" X = np.column_stack((pscenter, hhcenter, ncomact, rleader, dleader, inter1, inter2))\n",
" @deterministic\n",
" def idt(b1=beta, dl0=dL0):\n",
" print beta.value\n",
" fitted = np.exp(np.dot(X, b1) + W.value[stfips])\n",
" yy = (Y[:,:T] * np.outer(fitted, dl0))\n",
" return np.transpose(yy)\n",
" \n",
" dn_like = Poisson('dn_like', idt, value=dN, observed=True)\n",
" \n",
" return locals()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"variables = cox_spatialshared_model()\n",
"m3 = MCMC(variables, db='pickle', dbname='junk.pkl')\n",
"m3.use_step_method(AdaptiveMetropolis, m3.beta)\n",
"m3.sample(10)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"m3.beta.trace()"
],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment