# Vectorized algebraic differential equations

Hi everyone,

I’m relatively new to Stan and I’m trying to implement a pretty complicated differential algebraic equation model in cmdstanR.

I’m sharing the code below. I have a few questions that I hope you can help me with.

For efficiency reasons, I’m trying to vectorize the righthand side (rhs) function of the ODE. I first estimate the algebraic equations (auxiliaries) using previously determined parameters given in ‘aux_mat’. So in the rhs function I’m for estimating these, so that I can use them to estimate the derivatives of the ‘stocks’.

So my questions: (Thanks very much for any insight you can give me!)

• Currently, I’m getting an error using append_row in the rhs ODE function, " Ill-typed arguments supplied to function ‘append_row’ " but I haven’t been able to figure out how to fix it?
• Is this the right way to loop over ODE solutions of multiple individuals and assigning them to a single matrix that has the same shape as the data?
• Is my assignment of the measurement model good like this, or do I need to utilize another for-loop?
• Why do people add a “generated quantities” block to obtain a log likelihood? Is that necessary to run the model?
• Finally, Do you have any further suggestions to make the below code as efficient as possible? I will work with a model that contains about 10 stock variables (based on around 60 uncertain parameters), 15 auxiliaries and several constants. Furthermore, I will loop over 650 individuals, so efficiency is key.
``````functions {
real[] rhs(real time,               // Current timepoint
vector[] z,                // Initial state
vector[] constants,        // Constants values
matrix[] K,                // Matrix for stock parameters
matrix[] aux_mat,          // Matrix for auxiliary parameters
vector[] intercept_stocks, // Intercept of stocks
vector[] intercept_aux,    // Intercept of auxiliaries
int n_auxiliaries,         // Number of auxiliaries
int n_stocks) {            // Number of stocks
// First estimate the auxiliaries based on the stocks and the constants
vector[n_auxiliaries] auxiliaries = intercept_aux + aux_mat * append_row(z, constants);

// Next estimate the stocks
vector[n_stocks] dzdt = intercept_stocks + K * append_row(auxiliaries,
append_row(z, constants));
return dzdt;
}
}
data {
int<lower = 1> n_ind;                // Number of individuals
int<lower = 1> n_obs;                // Number of timesteps
int<lower = 1> n_stocks;             // Number of differential equations
int<lower = 1> n_auxiliaries;        // Number of auxiliary variables
int<lower = 1> n_constants;          // Number of constants
int<lower = 1> n_params;             // Number of parameters in equations
int<lower = 1> n_variables;          // Number of variables
real aux_mat[n_auxiliaries, n_stocks + n_constants]; // Auxiliary matrix
real y0[n_ind, n_stocks];            // The baseline data
real y[n_ind*(n_obs), n_stocks];     // The full data
real constants[n_ind, n_constants];  // The data of the constants
real t0;                             // The first time-point
real ts[n_obs];                      // The time-steps to be returned
}

parameters {
real intercept[n_stocks];
real theta[n_params];
real<lower = 0> sigma[n_stocks];
}

transformed parameters{
matrix[n_ind * n_obs, n_stocks] mu; // Output from the ODE solver
matrix[n_stocks, n_variables] K; // Matrix for stock parameters

K[1, 1] = theta[1];
K[1, 3] = -theta[2];
K[1, 5] = theta[3];
K[2, 2] = -theta[4];
K[2, 1] = -theta[5];
K[2, 4] = -theta[6];
K[2, 5] = -theta[7];

for (i in 1:n_ind) {
mu[i:n_obs*i,:] = ode_rk45(rhs, x0, t0, ts, K, aux_mat,
intercept_stocks, intercept_aux,
n_auxiliaries, n_stocks);
}
}

model {
// Priors
intercept[:] ~ normal(0, 1.0);
theta ~ normal(0, 0.5);
sigma ~ lognormal(0, 1);

// Measurement model
y ~ normal(mu, sigma);
}
``````

The ode_rk45 function is only available in cmdstan, not in rstan… as you need at least Stan version 2.24 (or so)…so try again with cmdstanr and use the latest and greatest 2.28.2 ?

Thanks for your reply! Sorry i said it wrong, I’m actually using cmdstanR. And I have cmdstan version 2.28.1. So my questions remain!

Any further comments on the above code?

For the `append_row` issue, it looks like what you want is to use `vector` for the function arguments but accidentally used `vector[]` which is for array of vectors.