Domain Error in Arguments for Gamma MLM

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!

1 Like