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
I appreciate any help you may give me :)
