I’m trying to run a simple model that uses pymc.ode.DifferentialEquation, but I’m getting an error that I don’t understand. I’ve modeled my code after the examples in this tutorial. Here’s my code:
import numpy as np
import pymc as pm
from pymc.ode import DifferentialEquation
# Simulate data.
t = np.arange(0, 2000, 50)
y_ = 30 - 20 * np.exp(-t/1000)
y_obs = np.random.normal(y_, 0.5)
def f(y, t, p):
return (1 / p[0]) * (30 - y[0])
ode_model = DifferentialEquation(func=f, times=t, n_states=1, n_theta=1, t0=0)
with pm.Model() as model:
tau = pm.HalfNormal("tau", 1000)
y0 = pm.Normal("y0", 10, 3)
sigma = pm.HalfNormal("sigma", 4)
ode_solution = ode_model(y0=[y0], theta=[tau])
y = pm.Normal("y", mu=ode_solution, sigma=sigma, observed=y_obs)
trace = pm.sample(2000, tune=1000, cores=1)
And here’s the traceback that I get when I run the model:
Traceback
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[7], line 7
4 sigma = pm.HalfNormal("sigma", 4)
6 ode_solution = ode_model(y0=[y0], theta=[tau])
----> 7 y = pm.Normal("y", mu=ode_solution, sigma=sigma, observed=y_obs)
9 trace = pm.sample(2000, tune=1000, cores=1)
File ~\Anaconda3\envs\pymc\Lib\site-packages\pymc\distributions\distribution.py:308, in Distribution.__new__(cls, name, rng, dims, initval, observed, total_size, transform, *args, **kwargs)
305 elif observed is not None:
306 kwargs["shape"] = tuple(observed.shape)
--> 308 rv_out = cls.dist(*args, **kwargs)
310 rv_out = model.register_rv(
311 rv_out,
312 name,
(...)
317 initval=initval,
318 )
320 # add in pretty-printing support
File ~\Anaconda3\envs\pymc\Lib\site-packages\pymc\distributions\continuous.py:521, in Normal.dist(cls, mu, sigma, tau, **kwargs)
515 sigma = pt.as_tensor_variable(sigma)
517 # tau = pt.as_tensor_variable(tau)
518 # mean = median = mode = mu = pt.as_tensor_variable(floatX(mu))
519 # variance = 1.0 / self.tau
--> 521 return super().dist([mu, sigma], **kwargs)
File ~\Anaconda3\envs\pymc\Lib\site-packages\pymc\distributions\distribution.py:387, in Distribution.dist(cls, dist_params, shape, **kwargs)
385 ndim_supp = cls.rv_op(*dist_params, **kwargs).owner.op.ndim_supp
386 create_size = find_size(shape=shape, size=size, ndim_supp=ndim_supp)
--> 387 rv_out = cls.rv_op(*dist_params, size=create_size, **kwargs)
389 rv_out.logp = _make_nice_attr_error("rv.logp(x)", "pm.logp(rv, x)")
390 rv_out.logcdf = _make_nice_attr_error("rv.logcdf(x)", "pm.logcdf(rv, x)")
File ~\Anaconda3\envs\pymc\Lib\site-packages\pytensor\tensor\random\basic.py:268, in NormalRV.__call__(self, loc, scale, size, **kwargs)
247 def __call__(self, loc=0.0, scale=1.0, size=None, **kwargs):
248 r"""Draw samples from a normal distribution.
249
250 Signature
(...)
266
267 """
--> 268 return super().__call__(loc, scale, size=size, **kwargs)
File ~\Anaconda3\envs\pymc\Lib\site-packages\pytensor\tensor\random\op.py:289, in RandomVariable.__call__(self, size, name, rng, dtype, *args, **kwargs)
288 def __call__(self, *args, size=None, name=None, rng=None, dtype=None, **kwargs):
--> 289 res = super().__call__(rng, size, dtype, *args, **kwargs)
291 if name is not None:
292 res.name = name
File ~\Anaconda3\envs\pymc\Lib\site-packages\pytensor\graph\op.py:295, in Op.__call__(self, *inputs, **kwargs)
253 r"""Construct an `Apply` node using :meth:`Op.make_node` and return its outputs.
254
255 This method is just a wrapper around :meth:`Op.make_node`.
(...)
292
293 """
294 return_list = kwargs.pop("return_list", False)
--> 295 node = self.make_node(*inputs, **kwargs)
297 if config.compute_test_value != "off":
298 compute_test_value(node)
File ~\Anaconda3\envs\pymc\Lib\site-packages\pytensor\tensor\random\op.py:334, in RandomVariable.make_node(self, rng, size, dtype, *dist_params)
329 elif not isinstance(rng.type, RandomType):
330 raise TypeError(
331 "The type of rng should be an instance of either RandomGeneratorType or RandomStateType"
332 )
--> 334 shape = self._infer_shape(size, dist_params)
335 _, static_shape = infer_static_shape(shape)
336 dtype = self.dtype or dtype
File ~\Anaconda3\envs\pymc\Lib\site-packages\pytensor\tensor\random\op.py:201, in RandomVariable._infer_shape(self, size, dist_params, param_shapes)
199 param_batched_dims = getattr(param, "ndim", 0) - param_ndim_supp
200 if param_batched_dims > size_len:
--> 201 raise ValueError(
202 f"Size length is incompatible with batched dimensions of parameter {i} {param}:\n"
203 f"len(size) = {size_len}, len(batched dims {param}) = {param_batched_dims}. "
204 f"Size length must be 0 or >= {param_batched_dims}"
205 )
207 if self.ndim_supp == 0:
208 return size
ValueError: Size length is incompatible with batched dimensions of parameter 0 DifferentialEquation{func=<function f at 0x000002772D8BCC20>, times=(0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650, 1700, 1750, 1800, 1850, 1900, 1950), n_states=1, n_theta=1, t0=0}.0:
len(size) = 1, len(batched dims DifferentialEquation{func=<function f at 0x000002772D8BCC20>, times=(0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250, 1300, 1350, 1400, 1450, 1500, 1550, 1600, 1650, 1700, 1750, 1800, 1850, 1900, 1950), n_states=1, n_theta=1, t0=0}.0) = 2. Size length must be 0 or >= 2
Can anyone tell me what I’m doing wrong?