I am dealing with an ODE model that depends on a forcing variable F(t) measured at each time step so that the target variable y[t] has the same dimension of F and the model is:
\frac{dy(t)}{dt} = -k \cdot y(t) \cdot F(t)
I am struggling defining the right function signature to generate some fake data to be used as Principled Bayesian Workflow:
functions {
vector swp (real t, vector y, real k, real F) {
vector[1] dydt;
dydt[1] = - k * y[1] * F;
return dydt;
}
}
data {
int<lower = 1> T;
vector[1] y0;
real t0;
array[T] real ts;
real k;
array[T] real F;
}
generated quantities {
array[T] vector[1] y_sim = ode_rk45(swp, y0, t0, ts, k, F);
for (t in 1:T) {
y_sim[t, 1] += normal_rng(0, 0.1);
}
}
compiling I get the following error :
Available signatures:
(<F1>, vector, real, array[] real, real, real) => vector
The 6th argument must be real but got array[] real
This won’t work as you need a smooth function on the ODE RHS rather than a random thing. I’d guess that stochastic ODEs are what you’d need to look into instead.
The error message you are getting is because in your definition of the swp function you have specified that the F argument is of real type but the F variable you pass to ode_rk45, which will be passed in as that argument, is an array of real values of length T. I think you could fix the immediate problem by changing the real F argument specifier in the swp definition to array[] real F, however, this would then I think result in an error when trying to evaluate the line
dydt[1] = - k * y[1] * F;
as F would now be an array of reals and so not compatible for assignment to a single array element.
More generally there is a mismatch between your mathematical model and how you are currently setting up the Stan model / specifying the problem. In your definition of the ODE system, F is a function of time. You say that F is “measured at each time step” - its not clear what you mean by “each timestep” here, but assuming you mean each time at which you will have observations of the system state y then this means you will still only have values for F(t) at a discrete subset of possible time t values. The ODE solver will need to evaluate the right-hand-side function swp and so F at arbitrary times not corresponding to the observation times.
You therefore need to map from the discrete F values you have to a function of time, and getting back to the issues @wds15 alluded to, this function will need to be smooth (at least Lipschitz continuous) for the ODE solver to be able to compute a solution (and for there to be a guarantee that a unique solution exists at all). If the F values you access to do correspond to direct noiseless measurements of the function F then one option would be to approximate F as an interpolant of these discrete values. If you used piecewise linear interpolation then I think that the ODE system has a tractable analytic solution and so you could skip the use of the ODE solver altogether.
If the F values are noisy / indirect measurements of F then you would need to define a suitable model for F to jointly infer - for example placing a Gaussian process prior on F.
Potentially you may approximate F piece-wise linear or constant. But then your will need index each F entry with time so in the RHS function you may condition F value on t. Alternatively you can do a loop with each F value solved by on ode call, just keep tag of the initial condition.
@wds15 Does Stan support stochastic ODE or should I use a different tool? @matt-graham I think that placing a gaussian process prior on F could become numerically infeasible having month of mesures of F taken every 20 minutes. @yizhang do you think splines would be a good idea? Should I build the spline inside the swp function or is there a way to know the value of t and then build only the spline tract used by the ode_rk function (I will have more or less 10k points)?
Glancing over the respective Wikipedia entry for SDEs it looks to me as if SDEs are not really supported in Stan. That is, you need to map the SDEs to things like partial ODE systems is what I understand from the Wiki page. However, wikipedia also points out that MC integration is a possible technique to solve the matter. This would suggest that a cleverly chosen scheme in Stan will do just that, but it’s not obvious to me instantly how this would work (so it may work).
A main deciding point of what is the appropriate modelling approach to use would be how smooth you believe F to be apriori. If you specify F as the path of a univariate SDE
where W is a Wiener noise process, a and B are the drift and diffusion coefficient terms respectively and \theta represents any free parameters in the system, then this will lead to rough F paths that are no-where differentiable. Smoother paths can be simulated using an augmented SDE system with degenerate diffusion coefficient. For example the SDE system
would correspond to having F be the (smoother) temporal integral of an auxiliary process V which itself is goverened by a univariate SDE as above. This is again related to the works on state space formulation of temporal Gaussian processes linked above, with the process for approximating smooth Gaussian processes with infinitely differentiable kernels like the squared exponential effectively involving recursively applying augmentations like above to compute increasingly smooth approximations.
In case you are interested in a simpler approach of using a piecewise linear interpolation for F, then the initial value problem over one time interval t \in [t_0, t_1] for which F(t_0) = F_0 and F(t_1) = F_1 and we linearly interpolate for intermediate t can be written as