Hello everyone,

I have what ought to be a simple question. I’m trying to include multiple binary indicator variables (dummy variables) within an ODE system: `B1`

and `B2`

either take the value `0`

or `1`

, enabling two different estimates for a parameter being fit to two groups of data; e.g., `(a1*B1 + a2*B2)`

.

I include `B1`

and `B2`

as real data `x_r`

in the `functions {}`

block, as well as the `data {}`

and `transformed data {}`

blocks (Stan code below).

The script compiles, but something is clearly going awry:

```
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: integrate_ode_rk45: continuous data[1] is nan, but must be finite! (in . . . at line 56)
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."
```

I’ve also tried only including a single indicator variable `B1`

that equals a single value of `1`

, such that part of my ODE is essentially: `a1 * 1`

, but this produces the same error. So I’m clearly doing something wrong in enabling `x_r`

, but I’ve tried all possible different ways of introducing and defining `x_r`

that I’m aware of.

Below is a simplified version of the model that produces the same error: I first simulate data in R, and then I display the Stan model. Note that in this simplified example, I don’t actually expect the parameter `a1`

and `a2`

to differ at all . . . this is just a simple example being used to diagnose the issue.

Finally, in case it’s at all useful, I also paste the full, actual Stan model I’m running, since the ODE structure and manner of calculation is a bit different (but again, even the most simplified ODE kicks up the same error).

Thank you – I very much appreciate the assistance!!

Simplified R code to simulate ODE system:

```
library(rstan)
library(ggplot2)
library(reshape2)
## Simulate ODE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## params
params <- c(a=3,
h=.013)
## create a list of initial conditions
R0 <- seq(50, 300, length.out=10)
## specify ODE function
resourceLoss <- function(t, R0, params) {
with(as.list(c(R0, params)),{
dR <- -((a * R0)/(1 + h * a * R0))
list(c(dR))
})
}
## set t
length.tlist <- 2
t.list <- seq(1, 2, length.out=length.tlist)
## run ODE
out <- ode(y=c(R0), times=t.list, func=resourceLoss, parms=params)
## save as data frame
loss <- as.data.frame(out)
## save list without Times for later use
dRdt.list <- round(out[,-1],3)
## create and add prefix for col names
loss.dat <- as.data.frame(dRdt.list)
colnames(loss.dat) <- paste("T", colnames(loss.dat), sep="")
## plot all dRdt
all_dRdt <- melt(loss, id.vars="time")
tplot <- ggplot(all_dRdt, aes(time, value, col=variable,group=variable)) +
geom_point() + geom_line(color="black") +
xlab("Time") + ylab("Resource density") + ggtitle("Kelp remaining = R0 - dRdT")
print(tplot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## add noise ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## reparameterized Gamma
my.rgamma<-function(n, mu, sigma){
rgamma(n, shape = (mu^2)/(sigma^2), scale = (sigma^2)/mu)}
## set dims for matrix and set sigma
rws=length(t.list) ## number of rows ie data points to simulate
cls=ncol(dRdt.list) ## number of columns ie separate FR runs
sig=9.50
## apply rgamma to deterministic data
mult_ODE <- function(n, dat, noise){
noise.out <- my.rgamma(n=length(dRdt.list), mu=dat, sigma=noise)
return(matrix(noise.out, nrow=rws, ncol=cls, byrow=FALSE))
}
noisy.dat <- round(mult_ODE(n, dat=dRdt.list, noise=sig),digits=3)
## replace noisy row 1 with deterministic R0
noisy.dat[1,] <- dRdt.list[1,]
## create data frame for plotting, add time as a column
noise.dat.df <- as.data.frame(noisy.dat)
noise.dat.df <- cbind(loss$time, noise.dat.df)
names(noise.dat.df)[1]<-"time"
## plot
melt_noise <- melt(noise.dat.df, id.vars="time")
noise_plot <- ggplot(melt_noise, aes(time, value, col=variable, group=variable)) +
geom_point() + geom_line(data=all_dRdt, aes(time, value, group=variable), col="black") +
xlab("Time") + ylab("Resource density") + ggtitle("Resources remaining = R0 - dRdT + noise")
print(noise_plot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## configure data list for STAN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## extract deterministic list of initial conditions
R0.list <- as.numeric(loss.dat[1,])
## number of columns ie ODE runs
dim.cols <- ncol(noisy.dat)
## create binary indicator variables
ind.length = 5
B1 <- c(rep(1, length.out=ind.length), rep(0, length.out=ind.length))
B2 <- c(rep(0, length.out=ind.length), rep(1, length.out=ind.length))
## list for Stan
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
t0 = t.list[1],
th = array(t.list[-1], dim=c(1)),
nLev = ncol(noisy.dat),
B1 = array(B1, dim=c(ncol(noisy.dat))),
B2 = array(B1, dim=c(ncol(noisy.dat))),
R0_obs = R0.list,
R_obs = array(noisy.dat[-1,], dim=c(1,dim.cols)))
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
```

and here’s a simplified Stan model:

```
functions {
real[] resourceLoss(
real t, // time
real[] R, // state
real[] theta, // params
real[] x_r,
int[] x_i) {
real dRdt[x_i[1]];
real a1 = theta[1];
real a2 = theta[2];
real h = theta[3];
real B1 = x_r[1]; // indicator variable 1
real B2 = x_r[2]; // indicator variable 2
for (i in 1:x_i[1]) {
dRdt[i] = -(((a1 * B1) + (a2 * B2)) * R[i])/(1 + ((a1 * B1) + (a2 * B2)) * h * R[i]);
}
return dRdt;
}
}
data {
int <lower=1> Tm; // Tm = Tmax = T-1 (T-initial pt) = 10
int <lower=1> nLev; // number of diff R0 i.e. # ODEs =
real R0_obs[nLev]; // Intial conditions R0 list
real <lower=1> t0; // starting time t = 0
real <lower=t0> th[Tm]; // array of length Tm i.e. 1
real R_obs[Tm, nLev]; // simulated drift remaining
real B1[nLev];
real B2[nLev];
}
transformed data{
real x_r[2];
int x_i[1];
//real x_r[1] = B1;
//real x_r[2] = B2;
x_i[1] = nLev;
}
parameters {
real <lower=0> a1;
real <lower=0> a2;
real <lower=0> h;
real <lower=0> sigma;
}
transformed parameters {
real theta[3];
theta[1] = a1;
theta[2] = a2;
theta[3] = h;
}
model {
real R[Tm, nLev];
real alpha[Tm, nLev];
real beta[Tm, nLev];
R = integrate_ode_rk45(resourceLoss, R0_obs, t0, th, theta, x_r, x_i);
for (t in 1:Tm) {
for (i in 1:nLev) {
alpha[t,i] = pow (R[t,i], 2) / pow (sigma, 2);
}
}
for (t in 1:Tm) {
for (i in 1:nLev) {
beta[t,i] = (1 / (pow (sigma, 2) / R[t,i]));
}
}
a1 ~ exponential(0.1);
a2 ~ exponential(0.1);
h ~ exponential(0.1);
sigma ~ exponential(0.1);
for (t in 1:Tm) {
for (i in 1:nLev) {
R_obs[t,i] ~ gamma(alpha[t,i], beta[t,i]);
}
}
}
```

In the interest of completeness, here’s the “real” system I’m wanting to fit in Stan.

dSdt = - a*S[t]*\frac{v - F[t]}{v}*q

dAdt = - a*A[t]*\frac{v - F[t]}{v}*(1 - q)

dFdt = (a*S[t]*\frac{v - F[t]}{v}*q) + (a* A[t]*\frac{v - F[t]}{v}*(1 - q)) - \frac{a*F[t]}{1 + a*b*F[t]}

And the full Stan model is below, though note I haven’t incorporated the third variable dAdt into the Stan model. But otherwise, this is the actual model I’m trying to run (I want to figure out the indicator variable before adding additional complexity).

```
functions {
real[] resourceLoss(
real t, // time
real[] Y, // state 2
real[] theta, // params
real[] x_r,
int[] x_i) {
real dR_dt;
real dF_dt;
real R = Y[1];
real F = Y[2];
real a = theta[1];
real v1 = theta[2];
real v2 = theta[3];
real v3 = theta[4];
real B1 = x_r[1];
real B2 = x_r[2];
real B3 = x_r[3];
dR_dt = (-(a * R) * (((v1*B1 + v2*B2 + v3*B3) - F) / (v1*B1 + v2*B2 + v3*B3)));
dF_dt = ((a * R) * (((v1*B1 + v2*B2 + v3*B3) - F) / (v1*B1 + v2*B2 + v3*B3)));
return {dR_dt, dF_dt};
}
}
data {
int <lower=1> Tm; // Tm = Tmax
int <lower=1> nLev; // number of diff R0 i.e. # ODEs = 35
real Y0_obs[nLev, 2]; // Intial conditions for R and F
real <lower=1> t0; // starting time t = 0
real <lower=t0> th[Tm]; // array of length Tm
real R_obs[Tm, nLev]; // experimental observations
real B1[nLev]; // indicator variable 1: 48hrs
real B2[nLev]; // indicator variable 2: 96hrs
real B3[nLev]; // indicator variable 3: 144hrs
}
transformed data{
real x_r[3];
int x_i[0];
//real x_r[1] = B1;
//real x_r[2] = B2;
//real x_r[3] = B3;
}
parameters {
real <lower=0> a;
real <lower=0> v1;
real <lower=0> v2;
real <lower=0> v3;
real <lower=0> sigma;
}
transformed parameters {
real theta[4];
theta[1] = a;
theta[2] = v1;
theta[3] = v2;
theta[4] = v3;
}
model {
real R[Tm, nLev, 2];
real alpha[Tm, nLev];
real beta[Tm, nLev];
for (n in 1:nLev) {
R[ , n, ] = integrate_ode_rk45(resourceLoss, Y0_obs[n], t0, th, theta, x_r, x_i);
}
for (t in 1:Tm) {
for (n in 1:nLev) {
alpha[t,n] = pow (R[t,n,1], 2) / pow (sigma, 2);
}
}
for (t in 1:Tm) {
for (n in 1:nLev) {
beta[t,n] = (1 / (pow (sigma, 2) / R[t,n,1]));
}
}
a ~ exponential(0.1);
v1 ~ normal(50, 15);
v2 ~ normal(40, 25);
v3 ~ normal(30, 25);
sigma ~ exponential(0.1);
for (t in 1:Tm) {
for (n in 1:nLev) {
R_obs[t,n] ~ gamma(alpha[t,n], beta[t,n]);
}
}
}
```

Thanks again for any insight you all may have!!