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.