# ODE model with Boolean control variable

Hello! I am trying to reproduce a published ODE model of methane production ​by cows. The model is quite simple:

\frac{dy}{dt}=c\,u-d\,y

where u is a Boolean control variable that takes values of 1 when the cow is fed, and 0 otherwise. c and d are unknown parameters and y is the production of methane by a single cow.

I have tried to recreate this model with fake data but I have not been successful so far. I think I am missing something in the sense I don’t know how to tell Stan about the initial value for the control variable.

Here is my Stan code:

functions{
real[] ch4_model(real t, real[] y, real[] theta, real[] x_r, int[] x_i){
int u;
real dydt[1];
if(u==0){
dydt[1] = -theta[2]*y[1];
}
else {
dydt[1] = theta[1]*u - theta[2]*y[1];
}
return dydt;
}
}

data {
int<lower=1> T;
int<lower=0,upper=1> u[T,1];
real<lower=0> y[T,1];
real t0;
real ts[T];
real y0[1];
real mu_c;
real s_c;
real mu_d;
real s_d;
}

transformed data{
real x_r[0];
int x_i[0];
}

parameters {
real<lower=0> sigma;
real<lower=0> c;
real<lower=0> d;
}

transformed parameters{
real y_hat[T,1];
{
real theta[2];
theta[1] = c;
theta[2] = d;
y_hat = integrate_ode_bdf(ch4_model, y0, t0, ts, theta, x_r, x_i, 1.0E-6, 1.0E-6, 1.0E4 );
}
}

model {
sigma ~ cauchy(0,1);
c ~ normal(mu_c,s_c);
d ~ normal(mu_d,s_d);
for (t in 1:T)
y[t] ~ normal(y_hat[t], sigma);
}


And here is the R script for calling it.

library(deSolve)

#generating simulated data

parameters <- c(c=1.8,d=1.8) #"real" parameter values
state <- c(y=0.05) #initial value

times <- seq(0, 48, length.out=200)

set.seed(1985)
u<-rbinom(200, 1, 0.5) #boolean control variable

signal <- data.frame(times=times,u=u)

input <- approxfun(signal, rule = 2)
#input(seq(from = 0, to = 48, length.out=400))

methane <- function(t,state,parameters){
with(as.list(c(state,parameters)),{
u <- input(t)
dy <- c*u - d*y
list(c(dy))

})
}

out <- ode(y=state,times=times,func=methane,parms=parameters)

out[,2] <- out[,2]+abs(rnorm(1,0,0.015)) #adding a little random noise

library(rstan)

T <- 199 #number time points excluding initial value
y0 <- 0.05 #initial methane measurement
dim(y0) <- 1
t0 <- 0 #initial time
ts <- out[2:nrow(out),1] #time points
y <- as.matrix(out[2:nrow(out),2]) #methane
u <- as.matrix(u[-1]) #control variable
mu_c <- 1
s_c <- 1
mu_d <- 1
s_d <- 1

fit <- stan("methane.stan",
data=c("y0","t0","y","T","ts","u","mu_c","s_c","mu_d","s_d"),
#cores=min(2,parallel::detectCores()),
init = 0.5,
chains=1,iter=2000)

fit


1 Like

I don’t think the built-in ODE solver is capable of handling discrete dynamics of this type directly, but I’ll admit I have limited knowledge of the ODE solver, so I’d ask @wds15 to confirm if I’m right.

There are however indirect options! Since for constant u, the system is linear, it can be solved analytically for each segment with constant u (or simply for each time step). That might in fact be much more efficient than using the solver as the solution at each timestep would have just a single exp call. (dy/dt = c - d*y - Wolfram|Alpha).

Alternatively, you can interpolate u smoothly (non-smooth interpolations, like the one you use in R will cause trouble for Stan) and then you can use the built-in integrator.

Does that make sense?

Best of luck with your model!

2 Likes

putting discontinuities into the ODE RHS may work ok, but it is quite “brutal” to the ODE solver. The better approach is to do two solves in a sequence.

Looking at your code though tells me that there is no discontinuity in the ODE RHS. Can you please elaborate a bit what goes wrong? Maybe a plot helps or a print of differences you are not expecting to be there. Something so that I don’t. have. to run the code if possible. Thanks.

3 Likes

thank you for your answer. Now you mentioned it, indeed I don’t see discontinuity in the RHS. This is the first time I am trying to fit a model with such a discontinuous variable so I did not know where to start.

In summary, after running the model I get parameter estimates equal to zero and my the predictions are stuck as a single line with mean equal to the initial value of y (look at a part of the summary of this fit below).

Inference for Stan model: methane.
1 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=1000.

mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
sigma         0.45    0.00 0.02  0.41  0.43  0.44  0.46  0.50   848    1
c             0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.00   867    1
d             0.00    0.00 0.00  0.00  0.00  0.00  0.00  0.01   780    1
y_hat[1,1]    0.05    0.00 0.00  0.05  0.05  0.05  0.05  0.05   677    1
y_hat[2,1]    0.05    0.00 0.00  0.05  0.05  0.05  0.05  0.05   677    1
y_hat[3,1]    0.05    0.00 0.00  0.05  0.05  0.05  0.05  0.05   677    1
...



You mentioned that one approach would be to do two solves in a sequence, how is that done?

I don’t think that you need to deal with two solves in a sequence, but have a read here:

https://computing.llnl.gov/projects/sundials/usage-notes

illustrating what I mean in more detail.

functions{
real[] ch4_model(real t, real[] y, real[] theta, real[] x_r, int[] x_i){
int u;
real dydt[1];
if(u==0){
dydt[1] = -theta[2]*y[1];
}
else {
dydt[1] = theta[1]*u - theta[2]*y[1];
}
return dydt;
}
}


Variable u is left uninitialised in this context! Since this an integer, you have to use the x_i argument from the integrate ode interface. So it should read like

int u = x_i[1];

for example. The integrate ode call then must vary the x_i argument according to what is in the u array.

Thank you again for your help. I have made a little progress I think, but there is something still missing (or more likely I did not implement correctly your suggestion). So well, now it looks like this (I explain in the comments some things):

functions{
real[] ch4_model(real t, real[] y, real[] theta, real[] x_r, int[] x_i){
int u = x_i[1];
real dydt[1];
// I deleted the if else statement since it seems to have no difference in the output
dydt[1] = theta[1]*u - theta[2]*y[1];

return dydt;
}
}

data {
int<lower=1> T;
int<lower=0,upper=1> u[T+1] ; //plus one because the number of observations of y does not consider the first measure
real<lower=0> y[T,1];
real t0;
real ts[T];
real y0[1];
real mu_c;
real s_c;
real mu_d;
real s_d;
}

transformed data{
real x_r[0];
int x_i[1];
}

parameters {
real<lower=0> sigma;
real<lower=0> c;
real<lower=0> d;
}

transformed parameters{
real y_hat[T,1];
{
real theta[2];
theta[1] = c;
theta[2] = d;
// I put u in the integrator instead of x_i, with x_i I get a very weird posterior for one of the parameters
y_hat = integrate_ode_bdf(ch4_model, y0, t0, ts, theta, x_r, u, 1.0E-6, 1.0E-6, 1.0E4 );
}
}

model {
sigma ~ cauchy(0,1);
c ~ normal(mu_c,s_c);
d ~ normal(mu_d,s_d);
for (t in 1:T)
y[t] ~ normal(y_hat[t], sigma);
}


Running this model as it is above, using very tight priors around the real values (N\sim(1.8,0.01)) for c and d I get the following output:

Inference for Stan model: methane.
1 chains, each with iter=6000; warmup=3000; thin=1;
post-warmup draws per chain=3000, total post-warmup draws=3000.

mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
sigma        0.59    0.00 0.03 0.53 0.57 0.59 0.61  0.65  2635    1
c            1.78    0.00 0.01 1.76 1.78 1.78 1.79  1.80  2432    1
d            1.82    0.00 0.01 1.80 1.81 1.82 1.82  1.84  2594    1
y_hat[1,1]   0.38    0.00 0.00 0.38 0.38 0.38 0.38  0.38  2425    1
y_hat[2,1]   0.59    0.00 0.00 0.59 0.59 0.59 0.60  0.60  2424    1
y_hat[3,1]   0.73    0.00 0.00 0.72 0.73 0.73 0.73  0.74  2426    1
y_hat[4,1]   0.82    0.00 0.01 0.81 0.82 0.82 0.82  0.83  2427    1


Thereafter, at some point, the posterior predictions stabilize at mean 0.98 which does not match the original data that looks like this:

As I mentioned in the comments of the code, when I put x_i in the integrator, instead of u, I get these parameter posteriors:

Inference for Stan model: methane.
1 chains, each with iter=6000; warmup=3000; thin=1;
post-warmup draws per chain=3000, total post-warmup draws=3000.

mean se_mean   sd      2.5%       25%       50%       75%     97.5% n_eff Rhat
sigma             0.49    0.00 0.02      0.44      0.47      0.49      0.50      0.54  2162    1
c                 0.00    0.00 0.00      0.00      0.00      0.00      0.00      0.00  2495    1
d                 1.80    0.00 0.01      1.78      1.79      1.80      1.81      1.82  2401    1
y_hat[1,1]        0.03    0.00 0.00      0.03      0.03      0.03      0.03      0.03  2640    1
y_hat[2,1]        0.02    0.00 0.00      0.02      0.02      0.02      0.02      0.02  2607    1
y_hat[3,1]        0.01    0.00 0.00      0.01      0.01      0.01      0.01      0.01  2579    1
y_hat[4,1]        0.01    0.00 0.00      0.00      0.01      0.01      0.01      0.01  2557    1


Which astonishes me, because regardless the tight prior I gave to c, the posterior is simply 0.

I hope you have time to guide me in what I am doing wrong.

It looks to me as if you want to have u change along the time. If that is what you want, then you have to loop over the ts vector and do the integration one-by-one, each time passing in the currently “active” u for the given time interval.

1 Like

You are looking at ODEs with static events. Since you are using desolve, I suggest take a look at deSolve: Forcing functions and Events and its user manual regarding events. Then you’ll see, as @wds15 suggested, you’ll need use u to control feeding events and define the RHS accordingly. Specifically, you can either

• loop through time points and use time-specific u in ode solver call, or
• do a single ode solver call while in the RHS definition check which time interval t falls into and use u for that interval.

@martinmodrak is also right that you don’t have to use ode solver, but the two options above still apply.

3 Likes