Hi all!
I’m trying to build a fitting model with an ODE solver, wherein the data I’m feeding the model is a stochastic time series that includes oscillations, while the process model is a deterministic time series that begins with dampening oscillations before reaching an equilibrium value. For clarity, the stochastic data was generated using the same process model included in my fitting model, but was made stochastic via the inclusion of demographic stochasticity.
Because of the mismatch between the stochastic data and the deterministic model, I’ve had a lot of trouble with model-fitting. It was then suggested to me that I may be able resolve my fitting issue by going for one-step-ahead prediction, rather than a global model fit. Essentially, I want estimation at a given time step to be dependent on the previous time step, rather than on the model’s estimation history.
My initial thought is that this could be done using a loop in the transformed parameters block, but I’m really not sure.
Conceptually, this all makes a lot of sense to me, but coding it? Another story entirely.
I’ve included my model below, and would desperately appreciate any advice people may have! :)
(P.S. I’m using rstan)
functions{
real[] dZ_dt(
real t, // Time
real[] Z, // System state {Parasite, Host}
real[] theta, // Parms
real[] x_r, // Real data; below is integer data
int[] x_i){
real P = Z[1]; // System state coded as an array, such that Z = (P,H)
real H = Z[2];
real r = theta[1]; // Parameters of the system, in order they appear
real O = theta[2];
real h = theta[3];
real b = theta[4];
real c = theta[5];
real u = theta[6];
real dP_dt = P*r - H*(O*P/(1 + O*P*h)); // Mechanistic model
real dH_dt = b + H*(c*(O*P/(1 + O*P*h))-u);
return({dP_dt,dH_dt}); // Return the system state
}
}
data{
int<lower=0>N; // Define N as non-negative integer
real ts[N]; // Assigns time points to N
real y0[2]; // Initial conditions for ODE
real t0; // Initial time point
int<lower=0>y[N,2]; // Define y as real and non-negative
}
parameters{
real<lower=0>r;
real<lower=0,upper=1>O;
real<lower=0>h;
real<lower=0>b;
real<lower=0>c;
real<lower=0,upper=1>u;
}
transformed parameters{
// Stiff solver for ODE
real Z[N,2]
= integrate_ode_bdf(dZ_dt,y0,t0,ts,{r, O, h, b, c, u},rep_array(0.0,2),rep_array(0,2),1e-10,1e-10,2e4);
for (i in 1:N) {
Z[i,1] = Z[i,1] + 1e-6;
Z[i,2] = Z[i,2] + 1e-6;
}
}
model{ // Ignore the means/variances for now...
r~normal(2.5,1); // r
O~beta(2,2); // O; bounded between 0 and 1
h~normal(0.06,1); // h
b~normal(35,1); // b
c~normal(0.2,1); // c
u~beta(2,2); // u; bounded between 0 and 1
for (k in 1:2){
y[ ,k]~poisson(Z[ ,k]);
}
}
generated quantities {
int y_rep[N, 2];
for (k in 1:2) {
for (n in 1:N)
y_rep[n,k] = poisson_rng(Z[n,k]);
}
}",
"Stan_Model_TypeII.stan")
stanc("Stan_Model_TypeII.stan") # To check that we wrote a file
Stan_Model_TypeII <- stan_model("Stan_Model_TypeII.stan")
# Squeezing the data into a form that Stan gets
Stoch_Data_TypeII = read.csv('/Users/mjarviscross/Desktop/GitHub/FcnalRespProj/data/Type II/FcnalRespProj_Data/Stoch_Data_TypeII.csv')
N <- length(Stoch_Data_TypeII$t)-1 # df is 2923; N is 2922
ts <- 1:N
y_init <- c(10,350) # Initial states, P = 10; H = 350
y0 <- array(y_init) # For the ODE solver
t0 <- 0.1
y <- as.matrix(Stoch_Data_TypeII[2:(N+1),3:4])
y <- cbind(y[ ,1],y[ ,2]); # This worked, sick; where y[,1] is P, and y[,2] is H
Stan_StochData_TypeII <- list(N=N,ts=ts,y0=y0,t0=t0,y=y)
# Fitting the data to the model
Stan_Model_TypeII <- stan_model("Stan_Model_TypeII.stan")
fit9 <- sampling(Stan_Model_TypeII,
data = Stan_StochData_TypeII,
warmup = 200,
iter = 400,
chains = 2,
cores = 2,
thin = 1,
control = list(max_treedepth=15,
adapt_delta=0.99),
seed = 123,
check_data = TRUE,
diagnostic_file = "/Users/mjarviscross/Desktop/GitHub/FcnalRespProj/Fit9.csv",
show_messages = TRUE,
verbose = TRUE)
