Hi there! I’m super new to Stan, so apologies in advance if I’ve missed something obvious here.
I’m fitting a model, and I want to specify varying upper and lower bounds for my parameters, and I figured out how to do that via section 21.4 in the User’s Guide. (using Q here instead of N which they used in the User’s Guide since I’m already using N for something else)
data{
.
int Q;
vector[Q] L; // lower bounds
vector[Q] U; // upper bounds
}
parameters{
vector<lower=0, upper=1>[Q] alpha_raw;
.
}
transformed parameters{
vector[Q] alpha = L + (U - L) .* alpha_raw; // this is the part that's supposed to give me my varying bounds
// stiff solver for ODE
real Z[N,2]
= integrate_ode_bdf(dZ_dt, Z_init, 0, ts, alpha, rep_array(0.0,0), rep_array(0,0));
}
However, I now want to pass my parameters to the stiff ODE solver, but it’s giving me an error as it requires the parameters argument to be defined as real[ ] but my parameters are now a vector.
The last line here is the problem,
real Z[N,2]
= integrate_ode_bdf(dZ_dt, Z_init, 0, ts, alpha, rep_array(0.0,0), rep_array(0,0));
The exact error is:
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Fifth argument to integrate_ode_bdf must have type real[ ]; found type = vector.
error in 'modelbb879b28_Updated_Fitting_Script' at line 44, column 85
-------------------------------------------------
42: // stiff solver for ODE
43: real Z[N,2]
44: = integrate_ode_bdf(dZ_dt, Z_init, 0, ts, alpha, rep_array(0.0,0), rep_array(0,0));
^
45: }
-------------------------------------------------
PARSER EXPECTED: ")"
Error in stanc(filename, allow_undefined = TRUE) :
failed to parse Stan model 'Updated_Fitting_Script' due to the above error.
Does anyone know how to either pass a vector to the stiff ODE solver, or tweak how I’ve defined my parameters/transformed parameters to meet the solvers real[ ] requirements?
Full Code:
functions{
real[] dZ_dt(
real t, // Time
real[] Z, // System state {Parasite, Host}
real[] alpha, // Parameters
real[] x_r, // Real data; below is integer data
int[] x_i){
real P = Z[1]; // System state coded as an array, such that Z = (P,H)
real H = Z[2];
real r = alpha[1]; // Parameters of the system, in order they appear
real O = alpha[2];
real h = alpha[3];
real b = alpha[4];
real c = alpha[5];
real u = alpha[6];
real dP_dt = P*r - H*(O*P/1 + O*P*h); // Mechanistic model
real dH_dt = b + H*(c*(O*P/1 + O*P*h)-u);
return{dP_dt,dH_dt}; // Return the system state
}
}
data{
int<lower=0>N; // Define N as non-negative integer
real ts[N]; // Assigns time points to N
real y_init[2]; // Initial conditions for ODE
real<lower=0>y[N,2]; // Define y as real and non-negative
int Q;
vector[Q] L; // lower bounds
vector[Q] U; // upper bounds
}
parameters{
vector<lower=0, upper=1>[Q] alpha_raw;
real<lower=0>Z_init[2]; // Initial population size non-neg
}
transformed parameters{
vector[Q] alpha = L + (U - L) .* alpha_raw;
// stiff solver for ODE
real Z[N,2]
= integrate_ode_bdf(dZ_dt, Z_init, 0, ts, alpha, rep_array(0.0,0), rep_array(0,0));
}
model{ // Ignore the means/variances for now...
alpha[{1}]~lognormal(0,1); // r
alpha[{2}]~lognormal(0,1); // O; bounded between 0 and 1
alpha[{3}]~lognormal(0,1); // h
alpha[{4}]~lognormal(0,1); // b
alpha[{5}]~lognormal(0,1); // c
alpha[{6}]~lognormal(0,1); // u; bounded between 0 and 1
Z_init~lognormal(log(10),1);
for (k in 1:2){
y_init[k]~poisson(Z_init[k]);
y[ ,k]~poisson(Z[ ,k]);
}
}
}
Thanks in advance for any assistance, apologies if I’m missing something really obvious! I’ve combed through other posts on here/the manual but can’t find anything that seems to answer my question.