Hi all,
Is there a way to fit the parameters of an ODE that not only depends on t and y, but also on some exogenous features X(t), using PyMC3? I’m using this as an example, but I’d like to move from y'=f(t,y,p) to y'=f(t,y,X,p), where X is a matrix of exogenous features that vary with time.
You can try writing the function f(y, t, p) and inside the function index to X using t.
I’m not sure I understand what you want to do. Can you give a more concrete example?
You can add time dependent terms to the ode using p. eg
y' = -p_0 t + p_1 \sin(t)
You might want to have a look at sunode, which is much faster than then the internal pymc3 implementation:
https://sunode.readthedocs.io/en/stable/quickstart_pymc3.html
You can install it with conda install -c conda-forge sunode by now, no need to build the package yourself.
Overall, I’d like to have something like this:
def f(y, t, p, X):
return -p[0]*t + p[1]*sin(t) + p[2]*X[t,0] + p[3]*X[t,1]
assuming X is a matrix of size (T,2), where T should be equal to the number of observations of y. I think @junpenglao was describing something like this, but wondering if I got the idea (and the pseudocode) correctly.
t is a continuous variable, not an integer. You could use something like p[2] * X[floor(t), 0], but you’d get a discontinuous right hand side that way. Unless you take care of that on the side of the solver it will give you incorrect results.
This has nothing to do with X vs p however. X[0, 0] and so on can just be entries of p.
Can you share what this is about? Then maybe we can suggest a better solution.
I thought with step_size = 1 it will just turn into a discrete ODE?
This example may clarify the problem. Basically, I’ve data from some input or control features X(t) that interact with a system. The system produces an outcome y(t). For simplicity, I’m assuming that y'=X(t)\beta, but it could be a non-linear relationship like eqs 2 and 3 of this reference. Given that I have access to the observed data y_{obs}(t) (red line in plot below) and the input data X(t), I’d like to use PyMC3 to recover the vector of coefficients \beta = [2,3]. To emphasize, I’m modeling the rate of change in y, not y itself.
# For reproducibility
np.random.seed(20394)
def f(y, t, p):
return np.dot(p[0][int(t),:], np.array([2,3]))
# Times for observation
times = np.arange(0,20,1)
# Input matrix
X = np.random.randn(20,2)
y = odeint(f, t=times, y0=0, args=tuple([[X]]))
yobs = np.random.normal(y,1)
fig, ax = plt.subplots(dpi=120)
plt.plot(times, yobs, label='observed data', linestyle='dashed', marker='o', color='red')
plt.plot(times, y, label='true function', color='k', alpha=0.5)
plt.legend()
plt.xlabel('Time (Seconds)')
plt.ylabel(r'$y(t)$');
plt.show()
Neither sunode or pymc3.ode have support for fixed step-size solvers, so you can’t set step_size=1 as things are. We could implement that if we think it is really helpful, but maybe a different strategy for the gradients would be better then. (Maybe use the autodiff of the approximated ode, instead of approximating the derivatives of the true ode solution)
So your ode looks like this:
Given some values x_i we define a function X(t):
The ODE is then
The right-hand-side X(t)\beta is not continuous though, which breaks assumptions of most ode solvers. Usually, you should solve the ode for each interval t\in (0, 1), t\in (1, 2) and so on seperately, or otherwise restart the solver at those discontinuities. The sundials doc has a note about this: https://computing.llnl.gov/projects/sundials/usage-notes at A2) Discontinuities in the RHS function
This is already a problem for your odeint call, not just inside a pymc model.
I think nothing stops you from doing this exact thing from your code with pymc3.ode, in sunode you will probably run into some (I think solvable) problems because np.floor(0.5) is still a double and can’t be used to index an array.
But maybe you can find a continuous version of X? I don’t see anything in the paper you linked that looks like a discontinuous right-hand-side (haven’t read it in detail…)
Thanks @aseyboldt. This makes sense. I’ll try to reformulate the problem so it doesn’t have the discontinuities.