Using ODE solver with varying upper and lower bounds

Modeling
, ,
1

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.

1 Like
2

To array 1d is your friend here.

3

Figured it out, thank you!

4

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

1 Like
5

Same model :)

6

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
7

Good luck to you both!

1 Like