2525def glm_hierarchical_model (random_seed = 123 ):
2626 """Sample glm hierarchical model to use in benchmarks"""
2727 np .random .seed (random_seed )
28- data = pd .read_csv (pm .get_data (' radon.csv' ))
29- data [' log_radon' ] = data [' log_radon' ].astype (theano .config .floatX )
28+ data = pd .read_csv (pm .get_data (" radon.csv" ))
29+ data [" log_radon" ] = data [" log_radon" ].astype (theano .config .floatX )
3030 county_idx = data .county_code .values
3131
3232 n_counties = len (data .county .unique ())
3333 with pm .Model () as model :
34- mu_a = pm .Normal (' mu_a' , mu = 0. , sd = 100 ** 2 )
35- sigma_a = pm .HalfCauchy (' sigma_a' , 5 )
36- mu_b = pm .Normal (' mu_b' , mu = 0. , sd = 100 ** 2 )
37- sigma_b = pm .HalfCauchy (' sigma_b' , 5 )
38- a = pm .Normal ('a' , mu = 0 , sd = 1 , shape = n_counties )
39- b = pm .Normal ('b' , mu = 0 , sd = 1 , shape = n_counties )
34+ mu_a = pm .Normal (" mu_a" , mu = 0.0 , sd = 100 ** 2 )
35+ sigma_a = pm .HalfCauchy (" sigma_a" , 5 )
36+ mu_b = pm .Normal (" mu_b" , mu = 0.0 , sd = 100 ** 2 )
37+ sigma_b = pm .HalfCauchy (" sigma_b" , 5 )
38+ a = pm .Normal ("a" , mu = 0 , sd = 1 , shape = n_counties )
39+ b = pm .Normal ("b" , mu = 0 , sd = 1 , shape = n_counties )
4040 a = mu_a + sigma_a * a
4141 b = mu_b + sigma_b * b
42- eps = pm .HalfCauchy (' eps' , 5 )
42+ eps = pm .HalfCauchy (" eps" , 5 )
4343 radon_est = a [county_idx ] + b [county_idx ] * data .floor .values
44- pm .Normal (' radon_like' , mu = radon_est , sd = eps , observed = data .log_radon )
44+ pm .Normal (" radon_like" , mu = radon_est , sd = eps , observed = data .log_radon )
4545 return model
4646
4747
@@ -50,24 +50,27 @@ def mixture_model(random_seed=1234):
5050 np .random .seed (1234 )
5151 size = 1000
5252 w_true = np .array ([0.35 , 0.4 , 0.25 ])
53- mu_true = np .array ([0. , 2. , 5. ])
54- sigma = np .array ([0.5 , 0.5 , 1. ])
53+ mu_true = np .array ([0.0 , 2.0 , 5.0 ])
54+ sigma = np .array ([0.5 , 0.5 , 1.0 ])
5555 component = np .random .choice (mu_true .size , size = size , p = w_true )
5656 x = np .random .normal (mu_true [component ], sigma [component ], size = size )
5757
5858 with pm .Model () as model :
59- w = pm .Dirichlet ('w' , a = np .ones_like (w_true ))
60- mu = pm .Normal ('mu' , mu = 0. , sd = 10. , shape = w_true .shape )
61- enforce_order = pm .Potential ('enforce_order' , tt .switch (mu [0 ] - mu [1 ] <= 0 , 0. , - np .inf ) +
62- tt .switch (mu [1 ] - mu [2 ] <= 0 , 0. , - np .inf ))
63- tau = pm .Gamma ('tau' , alpha = 1. , beta = 1. , shape = w_true .shape )
64- pm .NormalMixture ('x_obs' , w = w , mu = mu , tau = tau , observed = x )
59+ w = pm .Dirichlet ("w" , a = np .ones_like (w_true ))
60+ mu = pm .Normal ("mu" , mu = 0.0 , sd = 10.0 , shape = w_true .shape )
61+ enforce_order = pm .Potential (
62+ "enforce_order" ,
63+ tt .switch (mu [0 ] - mu [1 ] <= 0 , 0.0 , - np .inf )
64+ + tt .switch (mu [1 ] - mu [2 ] <= 0 , 0.0 , - np .inf ),
65+ )
66+ tau = pm .Gamma ("tau" , alpha = 1.0 , beta = 1.0 , shape = w_true .shape )
67+ pm .NormalMixture ("x_obs" , w = w , mu = mu , tau = tau , observed = x )
6568
6669 # Initialization can be poorly specified, this is a hack to make it work
6770 start = {
68- 'mu' : mu_true .copy (),
69- ' tau_log__' : np .log (1. / sigma ** 2 ),
70- ' w_stickbreaking__' : np .array ([- 0.03 , 0.44 ])
71+ "mu" : mu_true .copy (),
72+ " tau_log__" : np .log (1.0 / sigma ** 2 ),
73+ " w_stickbreaking__" : np .array ([- 0.03 , 0.44 ]),
7174 }
7275 return model , start
7376
@@ -77,26 +80,34 @@ class OverheadSuite:
7780 Just tests how long sampling from a normal distribution takes for various
7881 samplers
7982 """
83+
8084 params = [pm .NUTS , pm .HamiltonianMC , pm .Metropolis , pm .Slice ]
8185 timer = timeit .default_timer
8286
8387 def setup (self , step ):
8488 self .n_steps = 10000
8589 with pm .Model () as self .model :
86- pm .Normal ('x' , mu = 0 , sd = 1 )
90+ pm .Normal ("x" , mu = 0 , sd = 1 )
8791
8892 def time_overhead_sample (self , step ):
8993 with self .model :
90- pm .sample (self .n_steps , step = step (), random_seed = 1 ,
91- progressbar = False , compute_convergence_checks = False )
94+ pm .sample (
95+ self .n_steps ,
96+ step = step (),
97+ random_seed = 1 ,
98+ progressbar = False ,
99+ compute_convergence_checks = False ,
100+ )
92101
93102
94103class ExampleSuite :
95104 """Implements examples to keep up with benchmarking them."""
105+
96106 timeout = 360.0 # give it a few minutes
97107 timer = timeit .default_timer
98108
99109 def time_drug_evaluation (self ):
110+ # fmt: off
100111 drug = np .array ([101 , 100 , 102 , 104 , 102 , 97 , 105 , 105 , 98 , 101 ,
101112 100 , 123 , 105 , 103 , 100 , 95 , 102 , 106 , 109 , 102 , 82 ,
102113 102 , 100 , 102 , 102 , 101 , 102 , 102 , 103 , 103 , 97 , 97 ,
@@ -107,46 +118,52 @@ def time_drug_evaluation(self):
107118 100 , 101 , 102 , 103 , 97 , 101 , 101 , 100 , 101 , 99 ,
108119 101 , 100 , 100 , 101 , 100 , 99 , 101 , 100 , 102 , 99 ,
109120 100 , 99 ])
110-
111- y = pd .DataFrame ({
112- 'value' : np .r_ [drug , placebo ],
113- 'group' : np .r_ [['drug' ]* len (drug ), ['placebo' ]* len (placebo )]
114- })
121+ # fmt: on
122+
123+ y = pd .DataFrame (
124+ {
125+ "value" : np .r_ [drug , placebo ],
126+ "group" : np .r_ [["drug" ] * len (drug ), ["placebo" ] * len (placebo )],
127+ }
128+ )
115129 y_mean = y .value .mean ()
116130 y_std = y .value .std () * 2
117131
118132 sigma_low = 1
119133 sigma_high = 10
120134 with pm .Model ():
121- group1_mean = pm .Normal (' group1_mean' , y_mean , sd = y_std )
122- group2_mean = pm .Normal (' group2_mean' , y_mean , sd = y_std )
123- group1_std = pm .Uniform (' group1_std' , lower = sigma_low , upper = sigma_high )
124- group2_std = pm .Uniform (' group2_std' , lower = sigma_low , upper = sigma_high )
125- lambda_1 = group1_std ** - 2
126- lambda_2 = group2_std ** - 2
127-
128- nu = pm .Exponential (' ν_minus_one' , 1 / 29. ) + 1
129-
130- pm .StudentT (' drug' , nu = nu , mu = group1_mean , lam = lambda_1 , observed = drug )
131- pm .StudentT (' placebo' , nu = nu , mu = group2_mean , lam = lambda_2 , observed = placebo )
132- diff_of_means = pm .Deterministic (' difference of means' , group1_mean - group2_mean )
133- pm .Deterministic (' difference of stds' , group1_std - group2_std )
135+ group1_mean = pm .Normal (" group1_mean" , y_mean , sd = y_std )
136+ group2_mean = pm .Normal (" group2_mean" , y_mean , sd = y_std )
137+ group1_std = pm .Uniform (" group1_std" , lower = sigma_low , upper = sigma_high )
138+ group2_std = pm .Uniform (" group2_std" , lower = sigma_low , upper = sigma_high )
139+ lambda_1 = group1_std ** - 2
140+ lambda_2 = group2_std ** - 2
141+
142+ nu = pm .Exponential (" ν_minus_one" , 1 / 29.0 ) + 1
143+
144+ pm .StudentT (" drug" , nu = nu , mu = group1_mean , lam = lambda_1 , observed = drug )
145+ pm .StudentT (" placebo" , nu = nu , mu = group2_mean , lam = lambda_2 , observed = placebo )
146+ diff_of_means = pm .Deterministic (" difference of means" , group1_mean - group2_mean )
147+ pm .Deterministic (" difference of stds" , group1_std - group2_std )
134148 pm .Deterministic (
135- 'effect size' , diff_of_means / np .sqrt ((group1_std ** 2 + group2_std ** 2 ) / 2 ))
136- pm .sample (draws = 20000 , cores = 4 , chains = 4 ,
137- progressbar = False , compute_convergence_checks = False )
149+ "effect size" , diff_of_means / np .sqrt ((group1_std ** 2 + group2_std ** 2 ) / 2 )
150+ )
151+ pm .sample (
152+ draws = 20000 , cores = 4 , chains = 4 , progressbar = False , compute_convergence_checks = False
153+ )
138154
139155 def time_glm_hierarchical (self ):
140156 with glm_hierarchical_model ():
141- pm .sample (draws = 20000 , cores = 4 , chains = 4 ,
142- progressbar = False , compute_convergence_checks = False )
157+ pm .sample (
158+ draws = 20000 , cores = 4 , chains = 4 , progressbar = False , compute_convergence_checks = False
159+ )
143160
144161
145162class NUTSInitSuite :
146- """Tests initializations for NUTS sampler on models
147- """
163+ """Tests initializations for NUTS sampler on models"""
164+
148165 timeout = 360.0
149- params = (' adapt_diag' , ' jitter+adapt_diag' , ' jitter+adapt_full' , ' adapt_full' )
166+ params = (" adapt_diag" , " jitter+adapt_diag" , " jitter+adapt_full" , " adapt_full" )
150167 number = 1
151168 repeat = 1
152169 draws = 10000
@@ -159,32 +176,49 @@ def time_glm_hierarchical_init(self, init):
159176
160177 def track_glm_hierarchical_ess (self , init ):
161178 with glm_hierarchical_model ():
162- start , step = pm .init_nuts (init = init , chains = self .chains , progressbar = False , random_seed = 123 )
179+ start , step = pm .init_nuts (
180+ init = init , chains = self .chains , progressbar = False , random_seed = 123
181+ )
163182 t0 = time .time ()
164- trace = pm .sample (draws = self .draws , step = step , cores = 4 , chains = self .chains ,
165- start = start , random_seed = 100 , progressbar = False ,
166- compute_convergence_checks = False )
183+ trace = pm .sample (
184+ draws = self .draws ,
185+ step = step ,
186+ cores = 4 ,
187+ chains = self .chains ,
188+ start = start ,
189+ random_seed = 100 ,
190+ progressbar = False ,
191+ compute_convergence_checks = False ,
192+ )
167193 tot = time .time () - t0
168- ess = float (pm .ess (trace , var_names = [' mu_a' ])[' mu_a' ].values )
194+ ess = float (pm .ess (trace , var_names = [" mu_a" ])[" mu_a" ].values )
169195 return ess / tot
170196
171197 def track_marginal_mixture_model_ess (self , init ):
172198 model , start = mixture_model ()
173199 with model :
174- _ , step = pm .init_nuts (init = init , chains = self .chains ,
175- progressbar = False , random_seed = 123 )
200+ _ , step = pm .init_nuts (
201+ init = init , chains = self .chains , progressbar = False , random_seed = 123
202+ )
176203 start = [{k : v for k , v in start .items ()} for _ in range (self .chains )]
177204 t0 = time .time ()
178- trace = pm .sample (draws = self .draws , step = step , cores = 4 , chains = self .chains ,
179- start = start , random_seed = 100 , progressbar = False ,
180- compute_convergence_checks = False )
205+ trace = pm .sample (
206+ draws = self .draws ,
207+ step = step ,
208+ cores = 4 ,
209+ chains = self .chains ,
210+ start = start ,
211+ random_seed = 100 ,
212+ progressbar = False ,
213+ compute_convergence_checks = False ,
214+ )
181215 tot = time .time () - t0
182- ess = pm .ess (trace , var_names = ['mu' ])['mu' ].values .min () # worst case
216+ ess = pm .ess (trace , var_names = ["mu" ])["mu" ].values .min () # worst case
183217 return ess / tot
184218
185219
186- NUTSInitSuite .track_glm_hierarchical_ess .unit = ' Effective samples per second'
187- NUTSInitSuite .track_marginal_mixture_model_ess .unit = ' Effective samples per second'
220+ NUTSInitSuite .track_glm_hierarchical_ess .unit = " Effective samples per second"
221+ NUTSInitSuite .track_marginal_mixture_model_ess .unit = " Effective samples per second"
188222
189223
190224class CompareMetropolisNUTSSuite :
@@ -200,15 +234,21 @@ def track_glm_hierarchical_ess(self, step):
200234 if step is not None :
201235 step = step ()
202236 t0 = time .time ()
203- trace = pm .sample (draws = self .draws , step = step , cores = 4 , chains = 4 ,
204- random_seed = 100 , progressbar = False ,
205- compute_convergence_checks = False )
237+ trace = pm .sample (
238+ draws = self .draws ,
239+ step = step ,
240+ cores = 4 ,
241+ chains = 4 ,
242+ random_seed = 100 ,
243+ progressbar = False ,
244+ compute_convergence_checks = False ,
245+ )
206246 tot = time .time () - t0
207- ess = float (pm .ess (trace , var_names = [' mu_a' ])[' mu_a' ].values )
247+ ess = float (pm .ess (trace , var_names = [" mu_a" ])[" mu_a" ].values )
208248 return ess / tot
209249
210250
211- CompareMetropolisNUTSSuite .track_glm_hierarchical_ess .unit = ' Effective samples per second'
251+ CompareMetropolisNUTSSuite .track_glm_hierarchical_ess .unit = " Effective samples per second"
212252
213253
214254class DifferentialEquationSuite :
@@ -223,30 +263,34 @@ def freefall(y, t, p):
223263
224264 # Times for observation
225265 times = np .arange (0 , 10 , 0.5 )
226- y = np .array ([
227- - 2.01 ,
228- 9.49 ,
229- 15.58 ,
230- 16.57 ,
231- 27.58 ,
232- 32.26 ,
233- 35.13 ,
234- 38.07 ,
235- 37.36 ,
236- 38.83 ,
237- 44.86 ,
238- 43.58 ,
239- 44.59 ,
240- 42.75 ,
241- 46.9 ,
242- 49.32 ,
243- 44.06 ,
244- 49.86 ,
245- 46.48 ,
246- 48.18
247- ]).reshape (- 1 , 1 )
248-
249- ode_model = pm .ode .DifferentialEquation (func = freefall , times = times , n_states = 1 , n_theta = 2 , t0 = 0 )
266+ y = np .array (
267+ [
268+ - 2.01 ,
269+ 9.49 ,
270+ 15.58 ,
271+ 16.57 ,
272+ 27.58 ,
273+ 32.26 ,
274+ 35.13 ,
275+ 38.07 ,
276+ 37.36 ,
277+ 38.83 ,
278+ 44.86 ,
279+ 43.58 ,
280+ 44.59 ,
281+ 42.75 ,
282+ 46.9 ,
283+ 49.32 ,
284+ 44.06 ,
285+ 49.86 ,
286+ 46.48 ,
287+ 48.18 ,
288+ ]
289+ ).reshape (- 1 , 1 )
290+
291+ ode_model = pm .ode .DifferentialEquation (
292+ func = freefall , times = times , n_states = 1 , n_theta = 2 , t0 = 0
293+ )
250294 with pm .Model () as model :
251295 # Specify prior distributions for some of our model parameters
252296 sigma = pm .HalfCauchy ("sigma" , 1 )
@@ -263,4 +307,4 @@ def freefall(y, t, p):
263307 return np .mean ([ess .sigma , ess .gamma ]) / tot
264308
265309
266- DifferentialEquationSuite .track_1var_2par_ode_ess .unit = ' Effective samples per second'
310+ DifferentialEquationSuite .track_1var_2par_ode_ess .unit = " Effective samples per second"
0 commit comments