I’m simulating a system of differential equations (ODEs). Initially I used integrate_ode_rk45
as tutorials suggested. I then save obtained posterior and do a simulation for new conditions.
According to the deprecation error I changed ODE function from integrate_ode_rk45
to ode_rk45
. However, I can’t run expose_functions(..)
with the new ODE interface, which is necesary to simulate with new conditions.
An example to reproduce the error is below. (ODEs are taken from here):
library(brms)
library(tidyverse)
library(here)
library(cmdstanr)
library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())
# simulated data for exponential decay: c=A0*exp(-k1*time)
set.seed(1)
time <- c(0,1,2,3,4,5,7,9,11)
y_exp <- 100*exp(-0.2*time) # c0 =100, kr=0.2, exact simulated data
df_sim <- data.frame(time, y_exp)
df_sim <- df_sim %>% mutate(y_obs=y_exp+rnorm(length(time),0,5)) %>% dplyr::select(-y_exp)# random normal error added to simulated data with mean = 0 and sd = 5
N_t <- length(df_sim$time) # number of datapoints
fo_model2 <- "
vector fo2(real t, //time
vector y, // the state
real k1){ // the parameter next to A0
vector[1] dydt; // dimension of ODE
dydt[1] = - k1 * y[1]; // ODE
return dydt; // returns a 1-dimensional vector
}
// function called from brms: integration of ODEs
vector fo_ode2(data vector time, vector vA0, vector vk1, int N) {
vector[1] y0 = rep_vector(vA0[1], 1); //one initial value
//int N = size(time);
vector[1] y_hat_ode[N] = ode_rk45(fo2, y0 ,0.0,to_array_1d(time),vk1[1]);
vector[N] y_hat;
for(i in 1:N) y_hat[i] = y_hat_ode[i,1];
return(y_hat);
}
"
df_sim_ode <- df_sim[-1,] #remove the first datapoint at t0
fo_formula <- bf(y_obs~fo_ode2(time,A0,k1, 8),
A0~1,
k1~1,
nl=TRUE, loop=FALSE)
fo_priors <- c(prior(normal(100,10),nlpar=A0),
prior(normal(0.2,0.1), nlpar=k1),
prior(cauchy(0,10), class=sigma)
)
fo_result2 <- brm(
data=df_sim_ode,
family = gaussian,
formula=fo_formula,
prior=fo_priors,
inits=0,
iter=4000,
chains=4,
cores=4,
stanvars = stanvar(scode=fo_model2, block="functions"),
backend="cmdstanr",
file="fo_result2")
# Everything above runs correctly
# This function causes the error
expose_functions(fo_result2, vectorize = FALSE)
Produces this error:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Variable "fo2" does not exist.
error in 'modelf4cf9c8e4d57_User_defined_functions' at line 18, column 45
-------------------------------------------------
16: vector[1] y0 = rep_vector(vA0[1], 1); //one initial value
17: //int N = size(time);
18: vector[1] y_hat_ode[N] = ode_rk45(fo2, y0 ,0.0,to_array_1d(time),vk1[1]);
^
19: vector[N] y_hat;
-------------------------------------------------
Error in stanc(model_code = mc, model_name = "User-defined functions", :
failed to parse Stan model 'User-defined functions' due to the above error.
I’m using:
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
brms_2.18.0
cmdstanr_0.5.2
I’m wondering if anyone knows if it’s a bug or am I doing something wrong?