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);
}