Hi! I’m very new to this community and Bayesian statistics in general but I’m very excited about PYMC! I’m having a use case which is perhaps a bit uncommon. Indeed, I want to estimate parameters for an ODE model which has a time varying external input.
I was playing around with the free fall example from here: https://docs.pymc.io/en/v3/pymc-examples/examples/ode_models/ODE_API_introduction.html
Let’s say I want to add some pre-recorded wind pattern which influences the free fall. In that case, I need to access the time parameter of the ode function freefall(y, t, p) (which is ignored in the original GSoC example). However the ‘t’ parameter takes strange values. Here are the values it takes for each iteration:
There are two strange things happening in the end: first, time gets out of range. It is defined as: times = np.arange(0, 10, 0.5), so it should stop at 9.5. Even more strangely, on the last iteration, the parameter receives a non-numeric value ‘t’.
If I look at the types of these values, I see the following:
It is very strange that on the last iteration, the function receives an Aesara tensor instead of a float.
I don’t know if this is a bug or if I’m just not understanding how to use this class.
Here is the actual code for this example:
So the reason why you see different time values printed than what the elements of your time vector are, is that the ODE solvers (scipy.integrate.odeint in this case) have a variable stepsize. They adapt dynamically to how the state variable in your ODE system change and take smaller steps if there’s a lot of curvature, or larger steps when things are boring. In the end you get an interpolation back.
With time-varying external input, do you mean a parameter (elements of theta) that changes its value over time? How hard that’s going to be depends on the nature of this change: Is it a step function, or smooth? Does it depend on the state variables or just the time?
Thank you very much for the insightful reply! The external input of the system is a time series which does not have any analytical solution. Indeed, it corresponds to user input of a device. I want to do system identification where I fit model parameters given the system response to user input, the system being a set of differential equations.
The use case here is pharmacology where I want to model individual patient response by measuring amount of drug given and effect. I do have population parameter distribution which serve as priors and I want to use patient specific response to obtain more accurate individual parameter estimations.
I saw that scipy’s odeint is completely numerical so I didn’t anticipate that pymc needs analytical gradients. Do you think there is a way I could get this working with arbitrary input?
I couldn’t find any example online showing how to use the ‘t’ parameter in this context.
I understand now that analytical gradients are necessary. I reformulated my input as a symbolic SymPy function. However, I still don’t understand how to use the function inside the differential equation.
I modified the example code above code like so:
from sympy.abc import t
from sympy.printing.aesaracode import aesara_function
expr = t**2
input = aesara_function([t], [expr])
def freefall(y, t, p):
return 2.0 * p[1] - p[0] * y[0] + input(t)
This gives me the following errors:
TypeError: Bad input argument with name “t” to aesara function with name “/home/jaj/.conda/envs/pymc/lib/python3.10/site-packages/sympy/printing/aesaracode.py:509” at index 0 (0-based).
[…]
Expected an array-like object, but found a Variable: maybe you are trying to call a function on a (possibly shared) variable instead of a numeric array?
I also tried to formulate the model with sunode without success.
Is there some minimal example available making use of the ‘t’ parameter in an ODE model to get me started?
I seem to be making some progress but I still don’t have a working example. After reading some code, I found I could use the sympy function “aesara_code” in order to generate Aesara variables from sympy expressions:
I now actually provide a TensorVariable to PyMC when it expects one but I now get the error:
aesara.graph.utils.MissingInputError: Input 0 (t) of the graph (indices start from 0), used to compute Elemwise{pow,no_inplace}(t, TensorConstant{2}), was not provided and not given a value.
This happens in the codepath with the “aesara_code” call. The “return aef(time)” path works without issues.
It seems like the TensorVariable needs some additional initialization.
Can you post a small working example with some dummy data that includes the PyMC model? I have a suspicion of what is wrong – you’re mixing compiled aesara functions (aef) and symbolic tensors (generated by aesara_code). You can inject pre-defined PyMC variables into sympy code using a cache dictionary and the aesara_code function, but it’s not totally clear to me how t should be defined so I can’t give a concrete suggestion on how to implement it.
