I have been trying to fit a model using CmdStanR that involves solving an ODE in order to find a posterior for its parameters. My model currently solves this ODE in the model
block and runs fine. I would like it to solve the ODE in the transformed parameters
block instead so that I can output the solution to the ODE in the generated quantities
block (I understand this would not be possible currently as the variables in model
are locally defined).
The lines on which I solve the ODE are:
transformed parameters{
vector[Nt] y_solution[16];
y_solution = ode_bdf_tol(dy_dt, Y0s_solve, 0, time, 1e-6, 1e-6, 1000000, theta, kappa);
}
To get an error I move the above lines of code from the model
block to transformed parameters
and run the below in my R script:
mod_exp13_fixed_VI <- AML_mod$variational(
data = stan_data,
iter=10000,
tol_rel_obj = 0.001,
seed = 4543)
Which gives the following error:
Chain 1 Assertion failed!
Chain 1
Chain 1 Program: C:\Users\Taylo\OneDrive\Documents\.cmdstanr\cmdstan-2.28.2\AML\AML_Simdata.exe
Chain 1 File: stan/lib/stan_math/lib/eigen_3.3.9/Eigen/src/Core/DenseCoeffsBase.h, Line 408
Chain 1
Chain 1 Expression: index >= 0 && index < size()
Warning: Fitting finished unexpectedly! Use the $output() method for more information.
This error only appears at the end of the variational
output - it performs gradient ascent without issue and then only fails when it attempts to sample from the approximated posterior. I am not sure what this means but I thought it may be a helpful clue. This error does not appear if I integrate the ODE in the model
block.
-
Nt
is the integer 32 (the number of time points at which the ODE is solved) -
Y0s_solve
is a vector of initial conditions of length 16 -
time
is an array containing the numbers from 0.5 to 16 in steps of 0.5 (it is at these time points I wish to solve the ODE) -
theta
is an array of length 30 containing parameters fordy_dt
(from theparameters
block) -
kappa
is a real number that is also a parameter fordy_dt
but is entered in thedata
block.
The function dy_dt
is as below:
functions {
vector dy_dt(real t, vector y, real[] theta, real kappa) {
//this function defines the system of ODEs
vector[16] dy;
real z = y[1]+y[2]+y[3]+y[4]+y[5] + 2 *(y[6] + y[7] + y[8] + y[9]); // 'size' of cells in niche
//print("z = ", z);
dy[1] = theta[1]*y[1]*(kappa-z) - theta[2]*y[1];
dy[2] = theta[2]*y[1] + theta[3]*y[2]*(kappa-z) - (theta[4] + theta[5] + theta[6])*y[2];
dy[3] = theta[4]*y[2] - (theta[7] + theta[8])*y[3];
dy[4] = theta[5]*y[2] - (theta[9] + theta[10] + theta[11])*y[4];
dy[5] = theta[9]*y[4] - (theta[12] + theta[13])*y[5];
dy[6] = theta[16]*y[6]*(kappa-z) - theta[17]*y[6];
dy[7] = theta[17]*y[6] + theta[18]*y[7]*(kappa-z) - (theta[19] + theta[20] + theta[21])*y[7];
dy[8] = theta[19]*y[7] - (theta[22] + theta[23])*y[8];
dy[9] = theta[20]*y[7] - (theta[24] + theta[25])*y[9];
dy[10] = theta[6]*y[2] - theta[14]*y[10];
dy[11] = theta[7]*y[3] - theta[15]*y[11];
dy[12] = theta[10]*y[4] - theta[28]*y[12];
dy[13] = theta[13]*y[5] - theta[29]*y[13];
dy[14] = theta[21]*y[7] - theta[26]*y[14];
dy[15] = theta[22]*y[8] - theta[27]*y[15];
dy[16] = theta[25]*y[9] - theta[30]*y[16];
return dy;
}
}
I have noticed that the error is similar to the error in this thread though I cannot see where any kind of indexing error could have crept in with the two lines of code that would not exist using the same lines in the model
block.
Any help would be greatly appreciated and please let me know if I should provide more/different information as I am relatively new to this forum (have been using Stan for around 3 months). Thanks in advance.