How to specify ODE output when there are replicates

Modeling
1

Hi,
I have been confused for a while about the example 13.4 in the user guide. In the generated quantities section, the ODE solution is defined as below

vector[2] y_sim[T] = ode_rk45(sho, y0, t0, ts, theta);

I understand that the solution of ODE is a vector of 2 and the ODE is solved T times. My first question is that this kind of specification of y_sim is not considered as the optimal way (there will be a warning when the model is compiled). But it doesn’t work if it is array[T] vector[2] y_sim.

In this example, there is one solution based on one initial condition. What if I have a number of different initial conditions, then y_sim should be a array with [n_conditions, T, 2]. But the compiling won’t work. I wonder if anyone could give me a hint.

2

Can you please include the full Stan code generating that example and say which version of which Stan interface you’re using? My guess is that you’re using a version of RStan that’s using an old version of Stan that doesn’t support the new array syntax. The old array syntax has been deprecated and is scheduled for removal soon. I would suggest moving to CmdStanR or CmdStanPy, which are up to date with Stan and will support the array[T] syntax.

3

Hi Bob,

I am using cmdstanR 0.5.3 (CmdStan 2.30.1).

The full model is

functions {
vector sho(real t,
vector y,
real theta) {
vector[2] dydt;
dydt[1] = y[2];
dydt[2] = -y[1] - theta * y[2];
return dydt;
}
}
data {
int<lower=1> T;
array[T] vector[2] y;
real t0;
array[T] real ts;
}
parameters {
vector[2] y0;
vector<lower=0>[2] sigma;
real theta;
}
model {
vector[2] mu[T] = ode_rk45(sho, y0, t0, ts, theta);
sigma ~ normal(0, 2.5);
theta ~ std_normal();
y0 ~ std_normal();
for (t in 1:T) {
y[t] ~ normal(mu[t], sigma);
}
}

When I compile the example there is a warning:
Compiling Stan program…
Warning in ‘/var/folders/9x/86f88pr155s4lsxpdt623gx57n0ss1/T/Rtmpim4LXT/model-61687d701e5d.stan’, line 23, column 2: Declaration
of arrays by placing brackets after a variable name is deprecated and
will be removed in Stan 2.32.0. Instead use the array keyword before the
type. This can be changed automatically using the auto-format flag to
stanc

This is caused by the vector[2] mu[T] = ode_rk45(sho, y0, t0, ts, theta);

When I replace this with
array[T] vector [2] mu = ode_rk45(sho, y0, t0, ts, theta);
The warning is gone. This is good.

Now, if I have a number of replicates (n_rep), and I want to fit the ODE model to each of them. I attempted to do the model below:

functions {
vector sho(real t,
vector y,
real theta) {
vector[2] dydt;
dydt[1] = y[2];
dydt[2] = -y[1] - theta * y[2];
return dydt;
}
}
data {
int<lower=1> T;
array[T] vector[2] y;
real t0;
array[T] real ts;
int<lower=0> n_rep;
}
parameters {
vector[2] y0;
vector<lower=0>[2] sigma;
real theta;
}
model {
array[T,n_rep] vector[2] mu = ode_rk45(sho, y0, t0, ts, theta);
sigma ~ normal(0, 2.5);
theta ~ std_normal();
y0 ~ std_normal();
for (t in 1:T) {
y[t] ~ normal(mu[t], sigma);
}
}

But I got an error message
Semantic error in ‘/var/folders/9x/86f88pr155s4lsxpdt623gx57n0ss1/T/Rtmpim4LXT/model-61682954e7c1.stan’, line 24, column 2 to column 65:

22:  }
23:  model {
24:    array[T,n_rep] vector[2] mu = ode_rk45(sho, y0, t0, ts, theta);
       ^
25:    sigma ~ normal(0, 2.5);
26:    theta ~ std_normal();

Ill-typed arguments supplied to assignment operator =: lhs has type array[,] vector and rhs has type array vector

I couldn’t figure out myself and I hope you could provide me some hint.

4

I think I have found the solution
The key is to define mu as
array[n_rep, T] vector[2] mu;

Then
for (i in 1:n_rep) {
mu[i] = ode_rk45(sho, y0, t0, ts, theta);
}

mu[i] is an array (dimension = T) of vectors (length = 2).

Hopefully it is helpful for others who are interested