Hi,
I’m having trouble running a multi-level gamma distributed regression model with the SMC step function. I keep getting a “Domian Error in Arguments” error, which I believe comes from the fact that y_hat can be evaluated as zero in some cases, so when I calculate beta it gives less that or equal to zero. Is there anyway to fix this? The code and error is below
with pm.Model() as full_model:
int_mu = pm.Gamma('int_mu', 1, 3)
int_sd = pm.HalfCauchy('int_sd',10)
lvCT_mu = pm.Normal('lvCT_mu', 1, 1)
lvCT_sd = pm.HalfCauchy('lvCT_sd',1)
hvCT_mu = pm.Normal('hvCT_mu', 1, 1)
hvCT_sd = pm.HalfCauchy('hvCT_sd',5)
intercept = pm.Normal('intercept', mu=int_mu, sd=int_sd, shape=subjs)
lvCT = pm.Normal('lvCT', mu=lvCT_mu, sd=lvCT_sd, shape=subjs)
hvCT = pm.Normal('hvCT', mu=hvCT_mu, sd=hvCT_sd, shape=subjs)
# Model error
epsilon = pm.HalfCauchy('epsilon',1e2)
# Expected value
y_hat = intercept[rt_regression_data.subjs_idx]+\
lvCT[rt_regression_data.subjs_idx] * rt_regression_data.lvCT/10+ \
hvCT[rt_regression_data.subjs_idx] * rt_regression_data.hvCT/10
alpha_est=(y_hat**2/epsilon)
beta_est=(epsilon/y_hat)
y_like = pm.Gamma('y_like\n', alpha=alpha_est, beta=beta_est, observed=rt_regression_data.rt.values)
full_model_trace=pm.sample(500, chains=4, cores=10, n_init=1000, tune=3000, target_accept=.95,
random_seed=99, step=pm.SMC())
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
~/.local/lib/python3.6/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point,
givens, size)
575 try:
--> 576 return dist_tmp.random(point=point, size=size)
577 except (ValueError, TypeError):
~/.local/lib/python3.6/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
2534 dist_shape=self.shape,
-> 2535 size=size)
2536
~/.local/lib/python3.6/site-packages/pymc3/distributions/distribution.py in
generate_samples(generator, *args, **kwargs)
734 if dist_bcast_shape[:len(size_tup)] == size_tup:
--> 735 samples = generator(size=dist_bcast_shape, *args, **kwargs)
736 else:
~/.local/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in rvs(self, *args, **kwds)
961 if not np.all(cond):
--> 962 raise ValueError("Domain error in arguments.")
963
ValueError: Domain error in arguments.
During handling of the above exception, another exception occurred:
ValueError Traceback (most recent call last)
<ipython-input-251-15dba46ea8fc> in <module>
30
31 full_model_trace=pm.sample(500, chains=4, cores=10, n_init=1000, tune=3000,
target_accept=.95,
---> 32 random_seed=99, step=pm.SMC())
~/.local/lib/python3.6/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)
338 progressbar=progressbar,
339 model=model,
--> 340 random_seed=random_seed)
341 else:
342 if 'njobs' in kwargs:
~/.local/lib/python3.6/site-packages/pymc3/step_methods/smc.py in sample_smc(draws, step,
cores, progressbar, model, random_seed)
153
154 pm._log.info("Sample initial stage: ...")
--> 155 posterior, var_info = _initial_population(draws, model, variables)
156
157 while beta < 1:
~/.local/lib/python3.6/site-packages/pymc3/step_methods/smc_utils.py in _initial_population(draws,
model, variables)
17 var_info = {}
18 start = model.test_point
---> 19 init_rnd = pm.sample_prior_predictive(draws, model=model)
20 for v in variables:
21 var_info[v.name] = (start[v.name].shape, start[v.name].size)
~/.local/lib/python3.6/site-packages/pymc3/sampling.py in sample_prior_predictive(samples, model,
vars, var_names, random_seed)
1320 names = get_default_varnames(model.named_vars, include_transformed=False)
1321 # draw_values fails with auto-transformed variables. transform them later!
-> 1322 values = draw_values([model[name] for name in names], size=samples)
1323
1324 data = {k: v for k, v in zip(names, values)}
~/.local/lib/python3.6/site-packages/pymc3/distributions/distribution.py in draw_values(params,
point, size)
393 point=point,
394 givens=temp_givens,
--> 395 size=size)
396 givens[next_.name] = (next_, value)
397 drawn[(next_, size)] = value
~/.local/lib/python3.6/site-packages/pymc3/distributions/distribution.py in _draw_value(param, point,
givens, size)
583 with _DrawValuesContextBlocker():
584 val = np.atleast_1d(dist_tmp.random(point=point,
--> 585 size=None))
586 # Sometimes point may change the size of val but not the
587 # distribution's shape
~/.local/lib/python3.6/site-packages/pymc3/distributions/continuous.py in random(self, point, size)
2533 return generate_samples(stats.gamma.rvs, alpha, scale=1. / beta,
2534 dist_shape=self.shape,
-> 2535 size=size)
2536
2537 def logp(self, value):
~/.local/lib/python3.6/site-packages/pymc3/distributions/distribution.py in
generate_samples(generator, *args, **kwargs)
733 )
734 if dist_bcast_shape[:len(size_tup)] == size_tup:
--> 735 samples = generator(size=dist_bcast_shape, *args, **kwargs)
736 else:
737 samples = generator(size=size_tup + dist_bcast_shape, *args, **kwargs)
~/.local/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in rvs(self, *args, **kwds)
960 cond = logical_and(self._argcheck(*args), (scale >= 0))
961 if not np.all(cond):
--> 962 raise ValueError("Domain error in arguments.")
963
964 if np.all(scale == 0):
ValueError: Domain error in arguments.
Thanks in advance!