@aseyboldt Thank you for your code. I am running into an error when using your integration class.
Granted my integration is significantly more complicated than your example.
I am trying to recreate the Graphical Model for a Stop-Signal Reaction Time analysis found in the paper –http://dora.erbe-matzke.com/papers/Matzkeetal2012.pdf – There is some old code I found here -https://github.com/twiecki/stopsignal/blob/master/src/stop_likelihoods.pyx – that code has most of the liklihood functions in the older pymc and uses some of @twiecki’s cython code. I believe I have most of the setup correct except the integration. When running the pm.sample() I get a “TypeError: not all arguments converted during string formatting” error. I am unsure what that error means and where to begin debugging. I am sure it has to do with my integration code. If i need to make this its own post I will. But I am not sure if the issue is the integration code or my own.
def ExGaussCDF(value, mu, sigma, nu):
p = tt.switch(tt.gt(nu,0.05 * sigma),
tt.log(1-(std_cdf((value-mu)/sigma)
-std_cdf((value - mu - ((sigma**2)/nu))/sigma)
*tt.exp((((mu+((sigma**2)/nu))**2)-(mu**2)-2*value *((sigma**2)/nu))/(2*(sigma**2))))),
tt.log((1-(std_cdf(((value-mu)/sigma))))))
return p
def ExGaussPDF(value, mu, sigma, nu):
lp = tt.switch(tt.gt(nu, 0.05 * sigma),
- tt.log(nu) + (mu - value) / nu + 0.5 * (sigma / nu)**2
+ logpow(std_cdf((value - mu) / sigma - sigma / nu), 1.),
- tt.log(sigma * tt.sqrt(2 * np.pi))
- 0.5 * ((value - mu) / sigma)**2)
return lp
basic_model = pm.Model()
with basic_model:
# declare all variables for integration
x = tt.dscalar('x_')
x.tag.test_value = np.zeros(())
mu_go_ = tt.dscalar('mu_go_')
mu_go_.tag.test_value = np.ones(())
sigma_go_ = tt.dscalar('sigma_go_')
sigma_go_.tag.test_value = np.ones(())
nu_go_ = tt.dscalar('nu_go_')
nu_go_.tag.test_value = np.ones(())
mu_stop_ = tt.dscalar('mu_stop_')
mu_stop_.tag.test_value = np.ones(())
sigma_stop_ = tt.dscalar('sigma_stop_')
sigma_stop_.tag.test_value = np.ones(())
nu_stop_ = tt.dscalar('nu_stop_')
nu_stop_.tag.test_value = np.ones(())
ssd_ = tt.dscalar('ssd_')
ssd_.tag.test_value = np.ones(())
start = tt.dscalar('start')
start.tag.test_value = 0.001 # 1 ms
stop = tt.dscalar('stop')
stop.tag.test_value = 5. # 5 seconds
# Make function e^(ExGaussCDF(x,go)) * e^(ExGaussPDF(x-ssd,stop))
func = tt.exp(tt.switch(tt.gt(nu_go_,0.05 * sigma_go_),
tt.log(1-(std_cdf((x-mu_go_)/sigma_go_)
-std_cdf((x - mu_go_ - ((sigma_go_**2)/nu_go_))/sigma_go_)
*tt.exp((((mu_go_+((sigma_go_**2)/nu_go_))**2)-(mu_go_**2)-2*x *((sigma_go_**2)/nu_go_))/(2*(sigma_go_**2))))),
tt.log((1-(std_cdf(((x-mu_go_)/sigma_go_))))))) * \
tt.exp(tt.switch(tt.gt(nu_stop_, 0.05 * sigma_stop_),
- tt.log(nu_stop_) + (mu_stop_ - (x-ssd_)) / nu_stop_ + 0.5 * (sigma_stop_ / nu_stop_)**2
+ logpow(std_cdf(((x-ssd_) - mu_stop_) / sigma_stop_ - sigma_stop_ / nu_stop_), 1.),
- tt.log(sigma_stop_ * tt.sqrt(2 * np.pi))
- 0.5 * (((x-ssd_) - mu_stop_) / sigma_stop_)**2))
integrate = Integrate(func, x, mu_go_,sigma_go_,nu_go_,mu_stop_,sigma_stop_,nu_stop_,ssd_)
# Priors for unknown model parameters
mu_go = pm.Uniform('mu_go',.001,1)
sigma_go = pm.Uniform('sigma_go',.001,.3)
nu_go = pm.Uniform('nu_go',.001,.3)
# Priors for unknown stop model parameters
mu_stop = pm.Uniform('mu_stop',.001,.6)
sigma_stop = pm.Uniform('sigma_stop',.001,.25)
nu_stop = pm.Uniform('nu_stop',.001,.25)
def logp_SSRT(value,ssd_rts):
p = ExGaussPDF(value, mu_go, sigma_go, nu_go) + ExGaussCDF(value-ssd_rts, mu_stop, sigma_stop, nu_stop)
return tt.sum(p)
def logp_Inhibition(ssd_is):
p = 0
#sum over inhibit trials
for issd in ssd_is:
# print('SSD value: %.3f' % issd.tag.test_value)
p += tt.log(integrate(start,stop,mu_go,sigma_go,nu_go,mu_stop,sigma_stop,nu_stop,issd))
return p
goRTg = pm.ExGaussian('goRTg',mu=mu_go,sigma=sigma_go,nu=nu_go,observed = RTgo)
SSRT = pm.DensityDist('SSRT',logp_SSRT, observed={'value': np.array(SSRT_stop),'ssd_rts': ssd_rt})
Inhibition = pm.DensityDist('Inhib',logp_Inhibition, observed= {'ssd_is': ssd_inhib} )
trace = pm.sample(5000, burn = 1000)
Here is the total error output:
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-38-5b2c3f9ef829> in <module>()
87 Inhibition = pm.DensityDist('Inhib',logp_Inhibition, observed= {'ssd_is': ssd_inhib} )
88
---> 89 trace = pm.sample(5000, burn = 1000)
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, njobs, tune, nuts_kwargs, step_kwargs, progressbar, model, random_seed, live_plot, discard_tuned_samples, live_plot_kwargs, compute_convergence_checks, **kwargs)
371 start_, step = init_nuts(init=init, chains=chains, n_init=n_init,
372 model=model, random_seed=random_seed,
--> 373 progressbar=progressbar, **args)
374 if start is None:
375 start = start_
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/sampling.py in init_nuts(init, chains, n_init, model, random_seed, progressbar, **kwargs)
1352 'Unknown initializer: {}.'.format(init))
1353
-> 1354 step = pm.NUTS(potential=potential, model=model, **kwargs)
1355
1356 return start, step
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py in __init__(self, vars, max_treedepth, early_max_treedepth, **kwargs)
150 `pm.sample` to the desired number of tuning steps.
151 """
--> 152 super(NUTS, self).__init__(vars, **kwargs)
153
154 self.max_treedepth = max_treedepth
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/step_methods/hmc/base_hmc.py in __init__(self, vars, scaling, step_scale, is_cov, model, blocked, potential, integrator, dtype, Emax, target_accept, gamma, k, t0, adapt_step_size, step_rand, **theano_kwargs)
61
62 super(BaseHMC, self).__init__(vars, blocked=blocked, model=model,
---> 63 dtype=dtype, **theano_kwargs)
64
65 self.adapt_step_size = adapt_step_size
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py in __init__(self, vars, model, blocked, dtype, **theano_kwargs)
213
214 self._logp_dlogp_func = model.logp_dlogp_function(
--> 215 vars, dtype=dtype, **theano_kwargs)
216
217 def step(self, point):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/model.py in logp_dlogp_function(self, grad_vars, **kwargs)
656 varnames = [var.name for var in grad_vars]
657 extra_vars = [var for var in self.free_RVs if var.name not in varnames]
--> 658 return ValueGradFunction(self.logpt, grad_vars, extra_vars, **kwargs)
659
660 @property
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/pymc3/model.py in __init__(self, cost, grad_vars, extra_vars, dtype, casting, **kwargs)
397 self._cost, grad_vars, self._ordering.vmap)
398
--> 399 grad = tt.grad(self._cost_joined, self._vars_joined)
400 grad.name = '__grad'
401
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in grad(cost, wrt, consider_constant, disconnected_inputs, add_names, known_grads, return_disconnected, null_gradients)
603
604 rval = _populate_grad_dict(var_to_app_to_idx,
--> 605 grad_dict, wrt, cost_name)
606
607 for i in xrange(len(rval)):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in _populate_grad_dict(var_to_app_to_idx, grad_dict, wrt, cost_name)
1369 return grad_dict[var]
1370
-> 1371 rval = [access_grad_cache(elem) for elem in wrt]
1372
1373 return rval
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1369 return grad_dict[var]
1370
-> 1371 rval = [access_grad_cache(elem) for elem in wrt]
1372
1373 return rval
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in <listcomp>(.0)
1019 inputs = node.inputs
1020
-> 1021 output_grads = [access_grad_cache(var) for var in node.outputs]
1022
1023 # list of bools indicating if each output is connected to the cost
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_grad_cache(var)
1324 for idx in node_to_idx[node]:
1325
-> 1326 term = access_term_cache(node)[idx]
1327
1328 if not isinstance(term, gof.Variable):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gradient.py in access_term_cache(node)
1160
1161 input_grads = node.op.L_op(inputs, node.outputs,
-> 1162 new_output_grads)
1163
1164 if input_grads is None:
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/gof/op.py in L_op(self, inputs, outputs, output_grads)
709
710 def L_op(self, inputs, outputs, output_grads):
--> 711 return self.grad(inputs, output_grads)
712
713 def R_op(self, inputs, eval_points):
<ipython-input-13-44f8ebb860ef> in grad(self, inputs, grads)
39 dargs = []
40 for grad in grads:
---> 41 integrate = Integrate(grad, self._var, self._extra_vars)
42 darg = out * integrate(start, stop, *args)
43 dargs.append(darg)
<ipython-input-13-44f8ebb860ef> in __init__(self, expr, var, *extra_vars)
8 [var] + list(extra_vars),
9 self._expr,
---> 10 on_unused_input='ignore')
11
12 def make_node(self, start, stop, *extra_vars):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/compile/function.py in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
300 fn = orig_function(inputs, outputs,
301 mode=mode,
--> 302 accept_inplace=accept_inplace, name=name)
303 else:
304 # note: pfunc will also call orig_function -- orig_function is
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/compile/function_module.py in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys)
1816 mode = theano.compile.mode.get_mode(mode)
1817
-> 1818 inputs = list(map(convert_function_input, inputs))
1819 if outputs is not None:
1820 if isinstance(outputs, (list, tuple)):
~/Google_Drive/pymc3-dev/lib/python3.6/site-packages/theano/compile/function_module.py in convert_function_input(input)
1899 else:
1900 raise TypeError("Invalid input syntax: %s (check "
-> 1901 "documentation or use an In instance)" % orig)
1902 elif isinstance(input[0], SymbolicInput):
1903 if len(input) == 1:
TypeError: not all arguments converted during string formatting