0

I am a beginner with PyMC (https://github.com/pymc-devs/pymc) and am trying to construct a model with a dynamic component, essentially solving a small system of ordinary differential equations (ODEs) each time the model is called.

I have searched google and the (deprecated) PyMC mailing list to the best of my abilities and have come up with the model structure (pk_model.py) and caller (pk_fit.py) here:

https://gist.github.com/gyromagnetic/6097271

Running pk_fit seems to work at first, but then dumps many error messages related to the ODE solver. A standalone version of the ODEs and solver (not integrated with PyMC code) works fine.

Putting in various print statements, it appears that the code is at first working, but then at some point, the unknown parameter arguments (kcp, kpc, ke) change from scalars to numpy.ndarrays. This appears to be part of the problem.

Being a neophyte with PyMC, I expect that I am doing something obviously wrong.

I would greatly appreciate any help on this.

rcs
  • 67,191
  • 22
  • 172
  • 153
cytochrome
  • 549
  • 1
  • 4
  • 16

1 Answers1

2

You're correct: PyMC supplies the stochastic parameter values to your model as numpy.ndarrays, so all you need to do is convert back to a python scalar (get the first item) before computing the differentials.

For example, in your code:

dcc_dt = k1.item() * y_cp - (k2.item() + k3.item()) * y_cc
dcp_dt = k2.item() * y_cc - k1.item() * y_cp

Also, if you're dealing with large arrays of observed data, I suggest you look into converting your ODE model into a Python extension, using Cython, Boost.Python or similar. It will speed up the solution markedly, which is important if you need to sample many thousands of times in your MCMC. I used the latter, with the odeint library, for my ODE model and was able to get a speed up of 150x over the pure Python version.

Roger Stuckey
  • 31
  • 3
  • Would you mind adding the cython code here i.e. a complete example one can copy and paste?! That would be fantastic, thanks! – Cleb Mar 30 '18 at 14:21
  • 1
    @Cleb, there is a complete example [here](https://github.com/rstuckey/msd). There is a version that uses Numba in there too, which is another option. – Roger Stuckey Apr 12 '18 at 11:17