I am trying to recreate this example of bayesian PK/PD modelling using pymc3.....
The video shows the WinBUGS code and I am trying to convert to pymc3
https://www.youtube.com/watch?v=AQDXRoBan6Y
model here....
https://i.stack.imgur.com/3zfOv.jpg
WinBUGS code is here ....
https://i.stack.imgur.com/6LafO.jpg
My code is ....
from pymc3 import Model, Normal, Lognormal, Uniform
import numpy as np
import pandas as pd
data = pd.read_csv('/Users/Home/Documents/pymc3/fxa.data.csv' )
cobs = np.array(data['cobs'])
fxa = np.array(data['fxa.inh.obs'])
pkpd_model = Model()
with pkpd_model:
# Priors for unknown model parameters
emax = Uniform ('emax', lower =0, upper =100)
ec50 = Lognormal('ec50', mu=0, tau = 100000)
gamma = Uniform('gamma', lower=0, upper =10)
sigma = Uniform('sigma', lower = 0, upper = 1000 )
# Expected value of outcome
fxaMean = emax*(np.power(cobs, gamma)) / (np.power(ec50, gamma) + np.power(cobs, gamma))
# Likelihood (sampling distribution) of observations
fxa = Normal('fxa', mu=fxaMean, sd=sigma, observed=fxa )
But when I run the code I get the following error, which seems to relate to the way theano is interpreting the np.power function.
I am not sure how to proceed as I am a noob to pymc3 and theano and PK/PD modelling too!
Thanks in advance
Applied interval-transform to emax and added transformed emax_interval to model.
Applied log-transform to ec50 and added transformed ec50_log to model.
Applied interval-transform to gamma and added transformed gamma_interval to model.
Applied interval-transform to sigma and added transformed sigma_interval to model.
---------------------------------------------------------------------------
AsTensorError Traceback (most recent call last)
<ipython-input-28-1fa311a15ed0> in <module>()
14
15 # Likelihood (sampling distribution) of observations
---> 16 fxa = Normal('fxa', mu=fxaMean, sd=sigma, observed=fxa )
//anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
23 data = kwargs.pop('observed', None)
24 dist = cls.dist(*args, **kwargs)
---> 25 return model.Var(name, dist, data)
26 elif name is None:
27 return object.__new__(cls) # for pickle
//anaconda/lib/python2.7/site-packages/pymc3/model.pyc in Var(self, name, dist, data)
282 self.named_vars[v.name] = v
283 else:
--> 284 var = ObservedRV(name=name, data=data, distribution=dist, model=self)
285 self.observed_RVs.append(var)
286 if var.missing_values:
//anaconda/lib/python2.7/site-packages/pymc3/model.pyc in __init__(self, type, owner, index, name, data, distribution, model)
556 self.missing_values = data.missing_values
557
--> 558 self.logp_elemwiset = distribution.logp(data)
559 self.model = model
560 self.distribution = distribution
//anaconda/lib/python2.7/site-packages/pymc3/distributions/continuous.pyc in logp(self, value)
191 sd = self.sd
192 mu = self.mu
--> 193 return bound((-tau * (value - mu)**2 + T.log(tau / np.pi / 2.)) / 2.,
194 tau > 0, sd > 0)
195
//anaconda/lib/python2.7/site-packages/theano/tensor/var.pyc in __radd__(self, other)
232 # ARITHMETIC - RIGHT-OPERAND
233 def __radd__(self, other):
--> 234 return theano.tensor.basic.add(other, self)
235
236 def __rsub__(self, other):
//anaconda/lib/python2.7/site-packages/theano/gof/op.pyc in __call__(self, *inputs, **kwargs)
609 """
610 return_list = kwargs.pop('return_list', False)
--> 611 node = self.make_node(*inputs, **kwargs)
612
613 if config.compute_test_value != 'off':
//anaconda/lib/python2.7/site-packages/theano/tensor/elemwise.pyc in make_node(self, *inputs)
541 using DimShuffle.
542 """
--> 543 inputs = list(map(as_tensor_variable, inputs))
544 shadow = self.scalar_op.make_node(
545 *[get_scalar_type(dtype=i.type.dtype).make_variable()
//anaconda/lib/python2.7/site-packages/theano/tensor/basic.pyc in as_tensor_variable(x, name, ndim)
206 except Exception:
207 str_x = repr(x)
--> 208 raise AsTensorError("Cannot convert %s to TensorType" % str_x, type(x))
209
210 # this has a different name, because _as_tensor_variable is the
AsTensorError: ('Cannot convert [Elemwise{mul,no_inplace}.0 Elemwise{mul,no_inplace}.0\n Elemwise{mul,no_inplace}.0 ..., Elemwise{mul,no_inplace}.0\n Elemwise{mul,no_inplace}.0 Elemwise{mul,no_inplace}.0] to TensorType', <type 'numpy.ndarray'>)
