Error when solving ODE in transformed parameters block but not in model block

Modeling
,
1

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 for dy_dt (from the parameters block)
  • kappa is a real number that is also a parameter for dy_dt but is entered in the data 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.

2

Looks like a bug. Do you know how I can reproduce it in a minimalist model?

3

Hi @yizhang, sorry for the delay, I have been able to make a fairly simple model that still returns the same error.

The stan file I used is as below:

functions {
  vector dy_dt(real t, vector y, real[] theta){
    vector[1] dy;
    dy[1] = theta[1] * (1 - y[1]);
    
    return dy;
  }
}

data {
  int n; // number of data points
  real t[n]; // time of data points
  matrix[n, 1] y; // data at times t
  
  matrix[1, 1] y0; // Initial conditions
  
  int Nt; // number of time points at which we solve ODE
  real ts[Nt]; // time points at which we solve ODE
  int t_idx[n]; // index of ts for each entry of t
  
  real sigma_lik; // standard deviation
  }

parameters {
  real <lower=0.45, upper=0.55> theta[1];
}

transformed parameters {
  
  vector[Nt] y_solution[1];
  y_solution = ode_bdf_tol(dy_dt, to_vector(y0), 0, ts,
                           1e-6, 1e-6, 1000000, theta);
}

model {
  real sigma[1];
  real y_hat[Nt, 1];
  
  theta[1] ~ uniform(0.45, 0.55);
  
  for (i in 1:Nt){
    y_hat[i,] = to_array_1d(y_solution[i,]);
  }
  
  for (j in 1:n){
    sigma[1] = y_hat[t_idx[j], 1] * sigma_lik + 1e-6;
    y[j] ~ normal(y_hat[t_idx[j], 1], sigma[1]);
  }
}

This just tries to fit the parameter \theta in the ode \frac{\textrm{d}y}{\textrm{d}t} = \theta (1 - y). I have used slightly odd types considering we only have one parameter and one equation but this is just to emulate my full model as closely as possible.

The R script I used to run the model is as follows, I just generate some data from the solution and add a bit of noise:

n <- 64

theta <- 0.5
sigma_lik <- 0.05

t <- rep(seq(0.5, 16, by=0.5), each=n/32)
mu <- 1 + exp(-theta * t)
set.seed(25222)
y <- matrix(rnorm(n=64, mean = mu, sd = sigma_lik * mu), nrow=n)

y0 <- matrix(c(2), nrow=1)

Nt <- 32
ts <- seq(0.5, 16, by=0.5)
t_idx <- match(t, ts)

error_data <- list(n = n, t = t, y = y, y0 = y0, 
                   Nt = Nt, ts = ts, t_idx = t_idx,
                   sigma_lik = sigma_lik)

mod_file <- file.path(cmdstan_path(), "Error_example", "Error_model.stan")
err_mod <- cmdstan_model(mod_file)

vi_sample <- err_mod$variational(
  data = error_data,
  iter=1000,
  tol_rel_obj = 0.001,
  seed = 25222)

HMC_sample <- err_mod$sample(data = error_data, 
                                    seed = 25222,
                                    iter_sampling=100, iter_warmup=50)

Running either the VI or HMC model will give an error. If you cut out lines 30 and 31 from the stan file and move them to the appropriate position in the model block (immediately after theta[1] ~ uniform(0.45, 0.55)) the error disappears.

Thanks again for looking at this, I really appreciate it.

1 Like
4

This is the cause. For a 1D ODE with Nt outputs it should be vector[1] y_solution[Nt] or preferably with the array syntax array[Nt] vector[1] y_solution. The latest Stan 2.29 should emit the correct error msg but previous versions may not (hence the runtime error, but I don’t remember when it was fixed).

5

Thank you so much that’s sorted it! I’ve gone through and updated R, rstan and cmdstanr as well and it does give the correct error as well. Thank you again.