I am running the following model (first time Bayesian modelling/user of pymc). It estimates housing prices as a function of their attributes, spatially nested in ADA which are then spatially nested in a set of clusters. I’m using pymc 3.6 and theano 1.0.3.
with pm.Model() as contextual_effect:
# Priors
sigma_aC = pm.HalfNormal('sigma_aC', 16)
sigma_aA = pm.HalfNormal('sigma_aA', 277)
# Cluster branch model for slope
gammaC = pm.Normal('gammaC', mu=0., sigma=1, shape=9)
# Branch model for intercept
mu_aC = pm.Deterministic('mu_aC', gammaC[0] + gammaC[1]*size2 + gammaC[2]*size3 + gammaC[3]*size4 + gammaC[4]*size5 +gammaC[5]*size6
+ gammaC[6]*size7 + gammaC[7]*size8 + gammaC[8]*xbarC[clusterID])
# Cluster variation not explained by branch model
eps_aC = pm.Normal('eps_aC', mu=0, sigma=sigma_aC, shape=nclusters)
aC = pm.Deterministic('aC', mu_aC + eps_aC[clusterID])
# # ADA entropy model for slope
gammaA = pm.Normal('gammaA', mu=0., sigma=1, shape=10)
# Entropy model for intercept
mu_aA = pm.Deterministic('mu_aA', gammaA[0] + gammaA[1]*BArea + gammaA[2]*CArea + gammaA[3]*DArea + gammaA[4]*EArea + gammaA[5]*FArea + gammaA[6]*GArea + gammaA[7]*HArea + gammaA[8]*IArea + gammaA[9]*xbarA[ADA])
# # ADA variation not explained by entropy model
eps_aA = pm.Normal('eps_aA', mu=0, sigma=sigma_aA, shape=nADA)
aA = pm.Deterministic('aA', mu_aA + eps_aA[ADA])
# # Common slopes
b_cluster = pm.Normal('b_cluster', mu=0., sigma=10)
b_ADA = pm.Normal('b_ADA', mu=0., sigma=10)
# Model error
sigma_y = pm.Uniform('sigma_y', lower=0, upper=100)
# Individual level slopes
b_area = pm.Normal('b_area', mu=0.004, sigma=0.000018)
b_year = pm.Normal('b_year', mu=1.651504, sigma=0.03)
b_house2 = pm.Normal('b_house2', mu=-2.378035, sigma=0.1)
b_house3 = pm.Normal('b_house3', mu=-1.557363, sigma=0.1)
b_house4 = pm.Normal('b_house4', mu=-1, sigma=0.1)
# Expected value
y_hat = aC + aA + b_cluster * branch_measure + b_ADA * entropy_measure + b_area * area \
+ b_year * year + b_house2 * house2 + b_house3 * house3 + b_house4 * house4
# Data likelihood
y_like = pm.Normal('y_like', mu=y_hat, sigma=sigma_y, observed=price)
pm.model_to_graphviz(contextual_effect)`
I was having some issues with estimation of the model using NUTS (speed/convergence) so tried getting an initial estimate using ADVI. This seems to work (i.e. it converges and proceeds to using NUTS), but gives an error with the following
traceback:
RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py", line 110, in run
self._start_loop()
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py", line 160, in _start_loop
point, stats = self._compute_point()
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py", line 191, in _compute_point
point, stats = self._step_method.step(self._point)
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\step_methods\arraystep.py", line 247, in step
apoint, stats = self.astep(array)
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\step_methods\hmc\base_hmc.py", line 156, in astep
self.potential.update(hmc_step.end.q, hmc_step.end.q_grad, self.tune)
File "C:\Users\ITSLab\Anaconda3\lib\site-packages\pymc3\step_methods\hmc\quadpotential.py", line 278, in update
super().update(sample, grad)
TypeError: update() missing 1 required positional argument: 'tune'
"""
The above exception was the direct cause of the following exception:
TypeError Traceback (most recent call last)
TypeError: update() missing 1 required positional argument: 'tune'
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
<ipython-input-10-e3d5baf01b19> in <module>
1 with contextual_effect:
----> 2 contextual_effect_trace = pm.sample(1000, tune=1000, max_treedepth=20, target_accept=0.85, init='advi+adapt_diag_grad')
3 pm.summary(contextual_effect_trace, varnames=['gammaC', 'b_cluster', 'b_area', 'b_year', 'b_house2', 'b_house3', 'b_house4'])
~\Anaconda3\lib\site-packages\pymc3\sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, **kwargs)
435 _print_step_hierarchy(step)
436 try:
--> 437 trace = _mp_sample(**sample_args)
438 except pickle.PickleError:
439 _log.warning("Could not pickle model, sampling singlethreaded.")
~\Anaconda3\lib\site-packages\pymc3\sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, **kwargs)
967 try:
968 with sampler:
--> 969 for draw in sampler:
970 trace = traces[draw.chain - chain]
971 if (trace.supports_sampler_stats
~\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py in __iter__(self)
391
392 while self._active:
--> 393 draw = ProcessAdapter.recv_draw(self._active)
394 proc, is_last, draw, tuning, stats, warns = draw
395 if self._progress is not None:
~\Anaconda3\lib\site-packages\pymc3\parallel_sampling.py in recv_draw(processes, timeout)
295 else:
296 error = RuntimeError("Chain %s failed." % proc.chain)
--> 297 raise error from old_error
298 elif msg[0] == "writing_done":
299 proc._readable = True
RuntimeError: Chain 1 failed.
I’m running with init=‘auto’ but this is very slow (estimated time 10 hours). I initially estimated the model with this method, but it took a couple days and did not converge. I switched halfCauchy() to halfNormal(), and added max_treedepth=20, target_accept=0.85 based on suggestions I found online. The mu/sigma are also updated from initial values of mu=0 and sigma=1e5 (which I used based on a pymc tutorial with a similar specification).