 # Using ODE solver with varying upper and lower bounds

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; // System state coded as an array, such that Z = (P,H)
real H = Z;

real r = alpha; // Parameters of the system, in order they appear
real O = alpha;
real h = alpha;
real b = alpha;
real c = alpha;
real u = alpha;

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; // 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; // 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.

1 Like

To array 1d is your friend here.

Figured it out, thank you!

I’ve seen these equations before… @MadelineJC? =D

1 Like

Same model :)

Hi! Yes, we’re working on this project together, but aren’t in the same place, so are often trying to solve the same problem at the same time, blindly, lol.

1 Like

Good luck to you both!

1 Like