2
2
# <nbformat>3.0</nbformat>
3
3
4
4
# <codecell>
5
- from pylab import *
6
5
import pandas as pd
7
- from pymc3 import *
8
- from pymc3 .distributions .timeseries import *
6
+ from pylab import *
7
+
8
+ from pymc3 import StudentT , Model , NUTS , Normal , find_MAP , trace , get_data_file
9
+ from pymc3 .distributions .timeseries import GaussianRandomWalk
9
10
10
11
# <markdowncell>
11
12
15
16
# <codecell>
16
17
17
18
data = pd .read_csv (get_data_file ('pymc3.examples' , 'data/pancreatitis.csv' ))
18
- countries = ['CYP' , 'DNK' , 'ESP' , 'FIN' ,'GBR' , 'ISL' ]
19
+ countries = ['CYP' , 'DNK' , 'ESP' , 'FIN' , 'GBR' , 'ISL' ]
19
20
data = data [data .area .isin (countries )]
20
21
21
22
age = data ['age' ] = np .array (data .age_start + data .age_end )/ 2
28
29
# <codecell>
29
30
30
31
for i , country in enumerate (countries ):
31
- subplot (2 ,3 , i + 1 )
32
+ subplot (2 , 3 , i + 1 )
32
33
title (country )
33
34
d = data [data .area == country ]
34
35
plot (d .age , d .value , '.' )
35
36
36
- ylim (0 ,rate .max ())
37
+ ylim (0 , rate .max ())
37
38
38
39
# <markdowncell>
39
40
43
44
# <codecell>
44
45
45
46
nknots = 10
46
- knots = np .linspace (data .age_start .min (),data .age_end .max (), nknots )
47
+ knots = np .linspace (data .age_start .min (), data .age_end .max (), nknots )
47
48
48
49
49
- def interpolate (x0 ,y0 , x , group ):
50
+ def interpolate (x0 , y0 , x , group ):
50
51
x = np .array (x )
51
52
group = np .array (group )
52
53
53
54
idx = np .searchsorted (x0 , x )
54
55
dl = np .array (x - x0 [idx - 1 ])
55
56
dr = np .array (x0 [idx ] - x )
56
- d = dl + dr
57
+ d = dl + dr
57
58
wl = dr / d
58
59
59
60
return wl * y0 [idx - 1 , group ] + (1 - wl )* y0 [idx , group ]
@@ -68,7 +69,7 @@ def interpolate(x0,y0, x, group):
68
69
69
70
sd = StudentT ('sd' , 10 , 2 , 5 ** - 2 )
70
71
71
- vals = Normal ('vals' , p , sd = sd , observed = rate )
72
+ vals = Normal ('vals' , p , sd = sd , observed = rate )
72
73
73
74
# <markdowncell>
74
75
@@ -80,31 +81,30 @@ def interpolate(x0,y0, x, group):
80
81
with model :
81
82
s = find_MAP (vars = [sd , y ])
82
83
83
- step = NUTS (scaling = s )
84
+ step = NUTS (scaling = s )
84
85
trace = sample (100 , step , s )
85
86
86
87
s = trace [- 1 ]
87
88
88
89
step = NUTS (scaling = s )
89
90
91
+
90
92
def run (n = 3000 ):
91
93
if n == "short" :
92
94
n = 150
93
95
with model :
94
96
trace = sample (n , step , s )
95
-
96
-
97
97
# <codecell>
98
98
99
99
for i , country in enumerate (countries ):
100
- subplot (2 ,3 , i + 1 )
100
+ subplot (2 , 3 , i + 1 )
101
101
title (country )
102
102
103
103
d = data [data .area == country ]
104
104
plot (d .age , d .value , '.' )
105
- plot (knots , trace [y ][::5 ,:, i ].T , color = 'r' , alpha = .01 );
105
+ plot (knots , trace [y ][::5 , :, i ].T , color = 'r' , alpha = .01 )
106
106
107
- ylim (0 ,rate .max ())
107
+ ylim (0 , rate .max ())
108
108
109
109
110
110
if __name__ == '__main__' :
0 commit comments