5
5
import theano
6
6
import pymc3 as pm
7
7
from tqdm import tqdm
8
+ import warnings
8
9
9
10
from .arraystep import metrop_select
10
11
from .metropolis import MultivariateNormalProposal
@@ -57,7 +58,6 @@ class SMC():
57
58
816-832. `link <http://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9399
58
59
%282007%29133:7%28816%29>`__
59
60
"""
60
- default_blocked = True
61
61
62
62
def __init__ (self , n_steps = 5 , scaling = 1. , p_acc_rate = 0.01 , tune = True ,
63
63
proposal_name = 'MultivariateNormal' , threshold = 0.5 ):
@@ -88,14 +88,15 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
88
88
random_seed : int
89
89
random seed
90
90
"""
91
+ warnings .warn ("Warning: SMC is experimental, hopefully it will be ready for PyMC 3.6" )
91
92
model = modelcontext (model )
92
93
93
94
if random_seed != - 1 :
94
95
np .random .seed (random_seed )
95
96
96
97
beta = 0
97
98
stage = 0
98
- acc_rate = 0
99
+ acc_rate = 1
99
100
model .marginal_likelihood = 1
100
101
variables = model .vars
101
102
discrete = np .concatenate ([[v .dtype in pm .discrete_types ] * (v .dsize or 1 ) for v in variables ])
@@ -113,7 +114,6 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
113
114
likelihoods = np .array ([likelihood_logp (sample ) for sample in posterior ]).squeeze ()
114
115
beta , old_beta , weights , sj = _calc_beta (beta , likelihoods , step .threshold )
115
116
model .marginal_likelihood *= sj
116
- pm ._log .info ('Beta: {:f} Stage: {:d}' .format (beta , stage ))
117
117
# resample based on plausibility weights (selection)
118
118
resampling_indexes = np .random .choice (np .arange (draws ), size = draws , p = weights )
119
119
posterior = posterior [resampling_indexes ]
@@ -129,8 +129,10 @@ def sample_smc(draws=5000, step=None, progressbar=False, model=None, random_seed
129
129
if acc_rate == 0 :
130
130
acc_rate = 1. / step .n_steps
131
131
step .scaling = _tune (acc_rate )
132
- step .n_steps = 1 + (np .ceil ( np . log (step .p_acc_rate ) / np .log (1 - acc_rate )). astype ( int ))
132
+ step .n_steps = 1 + int (np .log (step .p_acc_rate ) / np .log (1 - acc_rate ))
133
133
134
+ pm ._log .info ('Stage: {:d} Beta: {:f} Steps: {:d} Acc: {:f}' .format (stage , beta ,
135
+ step .n_steps , acc_rate ))
134
136
# Apply Metropolis kernel (mutation)
135
137
proposed = 0.
136
138
accepted = 0.
0 commit comments