Best coding practices and examples for inference conditional on multidimensional ODE observations?

Modeling
1

Hi Everyone,

I’m presently looking to conduct inference from noisy observations of the states of an ODE system. I’ve previously worked on ODE parameter inferences in Stan conditional on a single vector of an observed “explicit variable” computed from the ODE states (so an observation element y_t would be compared with the model \hat{y}_t computed from multidimensional ODE states \hat{x}_t). I’m now looking to directly infer from multidimensional state observations, so y_t would be a vector of observations at time t, rather than a single element. And I would then compare y_t with a \hat{y}_t that is the ODE states \hat{x}_t plus observation noise. Would it be best to melt y_t and \hat{y}_t into single vectors at the likelihood evaluation stage, or would it be better to directly compare multidimensional arrays/matrix data types? Would folks like @wds15 and @charlesm93 happen to have examples of code where they’ve inferred from multidimensional state observations?

As an aside, after a long hiatus from using Stan, it’s real cool to see the additional updates and developments that have happened on the differential equation functionality front, such as the inclusion of the ckrk solver and the roll-out of adjoint sensitivity. @yizhang, @bbbales2, @wds15, @charlesm93, and others I’ve missed who contributed, much respect on the work you all have put in.

2

And from an old “best ODE practices” thread I dug up, I read the following:

@wds15, how do you switch solvers between the phases? I wasn’t previously aware you could do this.

3

You stop Stan and restart it. When restarting you give it the mass matrix manually…not very user friendly, I know.

4

Ah, I see. I think I may just try my luck with ode_ckrk for now then before doing any kind of manual tinkering.

Also, w.r.t. inference of multidimensional ODE states, would it be possible for me to do

y ~ normal(y_hat, sigma)

in the model block, where y and y_hat are multidimensional array[state_dim, N] real, or do I need to melt to a univariate structure? I recall from a couple of years ago that this vectorization was not possible with matrix types, but did not know if the operation was possible for a multidimensional array of reals.

Thanks again!

5

You can have an array of vectors.

6

Dug up some of your examples across various case studies and reports and put together the following Stan code:

functions {

  // Temperature function for ODE forcing.
  real temp_func(real t, real temp_ref, real temp_rise) {
    return temp_ref + (temp_rise * t) / (80 * 24 * 365) + 10 * sin((2 * pi() / 24) * t) + 10 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Exogenous SOC input function.
  real i_s_func(real t) {
    return 0.001 + 0.0005 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Exogenous DOC input function.
  real i_d_func(real t) {
    return 0.0001 + 0.00005 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Function for enforcing Arrhenius temperature dependency of ODE parameter.
  real arrhenius_temp(real input, real temp, real Ea, real temp_ref) {
    return input * exp(-Ea / 0.008314 * (1 / temp - 1 / temp_ref));
  }

  // Function for enforcing linear temperature dependency of ODE parameter.
  real linear_temp(real input, real temp, real Q, real temp_ref) {
    return input - Q * (temp - temp_ref);
  }

  /*
  AWB-ECA (equilibrium chemistry approximation) variant of AWB model.
  C[1] is soil organic carbon (SOC) density.
  C[2] is dissolved organic carbon (DOC) density.
  C[3] is microbial biomass carbon (MBC) density.
  C[4] is extracellular enzyme carbon (EEC) density.
  */

  vector AWB_ECA_ODE(real t, vector C,
                     real u_Q_ref, // Reference carbon use efficiency.
                     real Q, // Carbon use efficiency linear dependence negative slope.
                     real a_MSA, // AWB MBC-to-SOC transfer fraction.
                     real K_DE, // SOC decomposition K_m.
                     real K_UE, // DOC uptake K_m.
                     real V_DE_ref, // Reference SOC decomposition V_max.
                     real V_UE_ref, // Reference DOC uptake V_max.
                     real Ea_V_DE, // SOC V_max activation energy.
                     real Ea_V_UE, // DOC V_max activation energy.
                     real r_M, // MBC turnover rate.
                     real r_E, // Enzyme production rate.
                     real r_L, // Enzyme loss rate.
                     real temp_ref,
                     real temp_rise) {

    // Initiate exogenous input and forcing variables for future assignment.
    real temp;
    real i_s;
    real i_d;
    real u_Q; // Forced u_Q_ref.
    real V_DE; // Forced V_DE.
    real V_UE; // Forced V_UE.

    // Initiate dependent expressions.
    real F_S; // SOC decomposition
    real F_D; // DOC uptake

    // Assign input and forcing variables to appropriate value at time t.
    temp = temp_func(t, temp_ref, temp_rise); // x_r[1] is temp_ref 283.
    i_s = i_s_func(t);
    i_d = i_d_func(t);

    // Force temperature dependent parameters.
    u_Q = linear_temp(u_Q_ref, temp, Q, temp_ref);
    V_DE = arrhenius_temp(V_DE_ref, temp, Ea_V_DE, temp_ref);
    V_UE = arrhenius_temp(V_UE_ref, temp, Ea_V_UE, temp_ref);

    // Assign dependent expressions.
    F_S = V_DE * C[4] * C[1] / (K_DE + C[4] + C[1]);
    F_D = V_UE * C[3] * C[2] / (K_UE + C[3] + C[2]);

    // Compute derivatives.
    vector[4] dCdt; // Initiate vector for storing derivatives.
    dCdt[1] = i_s + a_MSA * r_M * C[3] - F_S;
    dCdt[2] = i_d + (1 - a_MSA) * r_M * C[3] + F_S + r_L * C[4] - F_D;
    dCdt[3] = u_Q * F_D * - (r_M + r_E) * C[3];
    dCdt[4] = r_E * C[3] - r_L * C[4];

    return dCdt;
  }

  // From https://discourse.mc-stan.org/t/rng-for-truncated-distributions/3122/12.
  real normal_lb_ub_rng(real mu, real sigma, real lb, real ub) {
      real p1 = normal_cdf(lb, mu, sigma);  // cdf with lower bound
      real p2 = normal_cdf(ub, mu, sigma);  // cdf with upper bound
      real u = uniform_rng(p1, p2);
      return (sigma * inv_Phi(u)) + mu;  // inverse cdf 
  }
}

data {
  int<lower=1> state_dim; // Number of state dimensions (4 for AWB).
  int<lower=1> N_t; // Number of observations.
  array[N_t] real<lower=0> ts; // Univariate array of observation time steps.
  array[N_t] vector<lower=0>[state_dim] y; // Multidimensional array of state observations bounded at 0.
  real<lower=0> temp_ref; // Reference temperature for temperature forcing.
  real<lower=0> temp_rise; // Assumed increase in temperature (°C/K) over next 80 years.
  real<lower=0> prior_scale_factor; // Factor multiplying parameter means to obtain prior standard deviations.
  real<lower=0> obs_error_scale; // Observation noise factor multiplying observations of model output x_hat.
  vector<lower=0>[state_dim] x_hat0; // Initial ODE conditions. [60.35767864, 5.16483124, 2.0068896, 0.99331202] for AWB-ECA instance of data.
  real<lower=0> u_Q_ref_mean;
  real<lower=0> Q_mean;
  real<lower=0> a_MSA_mean;
  real<lower=0> K_DE_mean;
  real<lower=0> K_UE_mean;
  real<lower=0> V_DE_ref_mean;
  real<lower=0> V_UE_ref_mean;
  real<lower=0> Ea_V_DE_mean;
  real<lower=0> Ea_V_UE_mean;
  real<lower=0> r_M_mean;
  real<lower=0> r_E_mean;
  real<lower=0> r_L_mean;
}

transformed data {
  real t0 = 0; // Initial time.
}

parameters {
  real<lower=0> u_Q_ref; // Reference carbon use efficiency.
  real<lower=0> Q; // Carbon use efficiency linear dependence negative slope.
  real<lower=0> a_MSA; // AWB MBC-to-SOC transfer fraction.
  real<lower=0> K_DE; // SOC decomposition K_m.
  real<lower=0> K_UE; // DOC uptake K_m.
  real<lower=0> V_DE_ref; // Reference SOC decomposition V_max.
  real<lower=0> V_UE_ref; // Reference DOC uptake V_max.
  real<lower=0> Ea_V_DE; // SOC V_max activation energy.
  real<lower=0> Ea_V_UE; // DOC V_max activation energy.
  real<lower=0> r_M; // MBC turnover rate.
  real<lower=0> r_E; // Enzyme production rate.
  real<lower=0> r_L; // Enzyme loss rate.
}

transformed parameters {
  array[N_t] vector<lower=0>[state_dim] x_hat = ode_ckrk(AWB_ECA_ODE, x_hat0, t0, ts, u_Q_ref, Q, a_MSA, K_DE, K_UE, V_DE_ref, V_UE_ref, Ea_V_DE, Ea_V_UE, r_M, r_E, r_L, temp_ref, temp_rise); 
}

model {
  u_Q_ref ~ normal(u_Q_ref_mean, u_Q_ref_mean * prior_scale_factor) T[0, 1];
  Q ~ normal(Q_mean, Q_mean * prior_scale_factor) T[0, 0.1];
  a_MSA ~ normal(a_MSA_mean, a_MSA_mean * prior_scale_factor) T[0, 1];
  K_DE ~ normal(K_DE_mean, K_DE_mean * prior_scale_factor) T[0, 5000];
  K_UE ~ normal(K_UE_mean, K_UE_mean * prior_scale_factor) T[0, 1];
  V_DE_ref ~ normal(V_DE_ref_mean, V_DE_ref_mean * prior_scale_factor) T[0, 1];
  V_UE_ref ~ normal(V_UE_ref_mean, V_UE_ref_mean * prior_scale_factor) T[0, 0.1];
  Ea_V_DE ~ normal(Ea_V_DE_mean, Ea_V_DE_mean * prior_scale_factor) T[5, 80];
  Ea_V_UE ~ normal(Ea_V_UE_mean, Ea_V_UE_mean * prior_scale_factor) T[5, 80];
  r_M ~ normal(r_M_mean, r_M_mean * prior_scale_factor) T[0, 0.1];
  r_E ~ normal(r_E_mean, r_E_mean * prior_scale_factor) T[0, 0.1];
  r_L ~ normal(r_L_mean, r_L_mean * prior_scale_factor) T[0, 0.1];

  // Likelihood evaluation.
  // y ~ normal(x_hat, obs_error_scale * x_hat);
  y ~ normal(x_hat, 1);
}

generated quantities {
  array[N_t] vector<lower=0>[state_dim] x_hat_post_pred;
  array[N_t] vector<lower=0>[state_dim] y_hat_post_pred;
  x_hat_post_pred = ode_ckrk(AWB_ECA_ODE, x_hat0, t0, ts, u_Q_ref, Q, a_MSA, K_DE, K_UE, V_DE_ref, V_UE_ref, Ea_V_DE, Ea_V_UE, r_M, r_E, r_L);
  y_hat_post_pred = normal_rng(x_hat_post_pred, obs_error_scale * x_hat_post_pred);
  real u_Q_ref_post = normal_lb_ub_rng(u_Q_ref_mean, u_Q_ref_mean * prior_scale_factor, 0, 1);
  real Q_post = normal_lb_ub_rng(Q_mean, Q_mean * prior_scale_factor, 0, 0.1);
  real a_MSA_post = normal_lb_ub_rng(a_MSA_mean, a_MSA_mean * prior_scale_factor, 0, 1);
  real K_DE_post = normal_lb_ub_rng(K_DE_mean, K_DE_mean * prior_scale_factor, 0, 5000);
  real K_UE_post = normal_lb_ub_rng(K_UE_mean, K_UE_mean * prior_scale_factor, 0, 1);
  real V_DE_ref_post = normal_lb_ub_rng(V_DE_ref_mean, V_DE_ref_mean * prior_scale_factor, 0, 1);
  real V_UE_ref_post = normal_lb_ub_rng(V_UE_ref_mean, V_UE_ref_mean * prior_scale_factor, 0, 0.1);
  real Ea_V_DE_post = normal_lb_ub_rng(Ea_V_DE_mean, Ea_V_DE_mean * prior_scale_factor, 5, 80);
  real Ea_V_UE_post = normal_lb_ub_rng(Ea_V_UE_mean, Ea_V_UE_mean * prior_scale_factor, 5, 80);
  real r_M_post = normal_lb_ub_rng(r_M_mean, r_M_mean * prior_scale_factor, 0, 0.1);
  real r_E_post = normal_lb_ub_rng(r_E_mean, r_E_mean * prior_scale_factor, 0, 0.1);
  real r_L_post = normal_lb_ub_rng(r_L_mean, r_L_mean * prior_scale_factor, 0, 0.1);
}

Eventually, I want to be able to do something like

 y ~ normal(x_hat, obs_error_scale * x_hat);

where the observation noise layered over the model output is proportional to the size of the observation. Of course, I’ve commented out for now because of a type issue. It looks like Stan doesn’t like having an array of vectors for the standard deviation. (It seems like I may need to do a loop for this?)

But then, trying something simpler with

y ~ normal(x_hat, 1);

I get the error when compiling the model,

Building: Semantic error:   -------------------------------------------------
   159:    // y ~ normal(x_hat, obs_error_scale * x_hat);
   160:    print(y);
   161:    y ~ normal(x_hat, 1);
           ^
   162:  }
   163:  
   -------------------------------------------------

Ill-typed arguments to '~' statement. No distribution 'normal' was found with the correct signature.

so it seems like normal isn’t taking arrays of vectors and I’m missing a step?

Thanks for helping me look into this, and I apologize if I made a glaring mistake that should have been obvious.

7

I thought you would want a multi normal cholesky data likelihood, no? Aren’t your observations like multi dimensional observations? Or is it not needed to model the correlation?

8

Yes, I made it so that the data y being fit was sampled from independent normal distributions with an ODE solution x serving as the mean and the observation noise being a constant proportion of x. So, each y_{i, t} \sim \mathcal{N}(x_{i,t}, 0.1 * x_{i,t}).

I’m trying to see if I can get back the true parameter values for the known data generating process with Stan (using the Stan results as a baseline to compare with other algorithms). So I’m looking for just a vectorised univariate normal for the likelihood and want to see if that were possible. I suppose I could use the multivariate normal Cholesky likelihood with L = 0.1 * \hat{x}_{i, t} with diag_matrix(vector x)? That said, based on the documentation, does multi_normal_cholesky take an array of vectors, or do I need to do further type conversions, since the documentation suggests an input of vectors y for mu.

Would using multi_normal_cholesky with a diagonal matrix be the best approach for what I’m trying to do? In PyTorch, using the vectorised normal distribution proved to be a smoother process for testing out a similar inference goal, but I’m not sure the Stan normal distribution can do the same kind of vectorisation?

9

So multi normal cholesky does take an array of vectors, sure.

However, as you go anyway univariate, you should create one vector per output which contains all the time points. Then you do one call to normal loss per output (with all time points each time).

10

I’m having a bit of trouble solidifying this, my apologies. To clarify, this sounds like I should melt all the observations from the [state_dim, N_t] shape to [1, state_dim * N_t]? Would it be possible to expound on what you mentioned in pseudocode, particularly how to handle the transformation from the array[N_t] vector[state_dim] data structure that the model output comes in right now to the single vector? And does it make sense to have the ODE function instead return a different data type that is not array[N_t] vector[state_dim] to facilitate the transformation to 1D vector?

Thank you for your patient help, Sebastian!

11

The ode output cannot be changed. You get an array of length t with a vector of length as you have states. You should repack that into a reversed ordering…so an array of length as many states you have and each entry is a vector of length t (one for each state)…let me know if pseudo code is (still) needed…on my ipad this does not work too well…

12

Thanks for getting back to me so quickly. I’ll try to go without pseudo-code for now. So, basically, I need to reshape from array[N_t] vector[state_dim] to array[state_dim] vector[N_t]. What function(s) in Stan will be able to handle the transformation for the model output? (For the data, I can just read that in as array[state_dim] vector[N_t] directly, so transforming the data won’t be an issue.) I’m hoping there is a way which would not require a loop. I’m familiar with going from array → matrix, for example, and Stan has easy functions for that, but reshaping arrays of vectors I’m not familiar with.

And then y ~ normal(x_hat, obs_error_scale * x_hat) will work as is without needing to rework that statement if y and x_hat are in the array[state_dim] vector[N_t] format? Looking at this example from @charlesm93, I suppose a for loop for the likelihood evaluation could also do, where

for (i in 1:state_dim) {
  y[i, ] ~ normal(x_hat[i, ], obs_error_scale * x_hat[i, ]);
}

if the likelihood evaluation vectorisation only extends to univariate vectors.

13

I have some code that initially compiles in the following:

functions {

  // Temperature function for ODE forcing.
  real temp_func(real t, real temp_ref, real temp_rise) {
    return temp_ref + (temp_rise * t) / (80 * 24 * 365) + 10 * sin((2 * pi() / 24) * t) + 10 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Exogenous SOC input function.
  real i_s_func(real t) {
    return 0.001 + 0.0005 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Exogenous DOC input function.
  real i_d_func(real t) {
    return 0.0001 + 0.00005 * sin((2 * pi() / (24 * 365)) * t);
  }

  // Function for enforcing Arrhenius temperature dependency of ODE parameter.
  real arrhenius_temp(real input, real temp, real Ea, real temp_ref) {
    return input * exp(-Ea / 0.008314 * (1 / temp - 1 / temp_ref));
  }

  // Function for enforcing linear temperature dependency of ODE parameter.
  real linear_temp(real input, real temp, real Q, real temp_ref) {
    return input - Q * (temp - temp_ref);
  }

  /*
  AWB-ECA (equilibrium chemistry approximation) variant of AWB model.
  C[1] is soil organic carbon (SOC) density.
  C[2] is dissolved organic carbon (DOC) density.
  C[3] is microbial biomass carbon (MBC) density.
  C[4] is extracellular enzyme carbon (EEC) density.
  */

  vector AWB_ECA_ODE(real t, vector C,
                     real u_Q_ref, // Reference carbon use efficiency.
                     real Q, // Carbon use efficiency linear dependence negative slope.
                     real a_MSA, // AWB MBC-to-SOC transfer fraction.
                     real K_DE, // SOC decomposition K_m.
                     real K_UE, // DOC uptake K_m.
                     real V_DE_ref, // Reference SOC decomposition V_max.
                     real V_UE_ref, // Reference DOC uptake V_max.
                     real Ea_V_DE, // SOC V_max activation energy.
                     real Ea_V_UE, // DOC V_max activation energy.
                     real r_M, // MBC turnover rate.
                     real r_E, // Enzyme production rate.
                     real r_L, // Enzyme loss rate.
                     real temp_ref,
                     real temp_rise) {

    // Initiate exogenous input and forcing variables for future assignment.
    real temp;
    real i_s;
    real i_d;
    real u_Q; // Forced u_Q_ref.
    real V_DE; // Forced V_DE.
    real V_UE; // Forced V_UE.

    // Initiate dependent expressions.
    real F_S; // SOC decomposition
    real F_D; // DOC uptake

    // Assign input and forcing variables to appropriate value at time t.
    temp = temp_func(t, temp_ref, temp_rise); // x_r[1] is temp_ref 283.
    i_s = i_s_func(t);
    i_d = i_d_func(t);

    // Force temperature dependent parameters.
    u_Q = linear_temp(u_Q_ref, temp, Q, temp_ref);
    V_DE = arrhenius_temp(V_DE_ref, temp, Ea_V_DE, temp_ref);
    V_UE = arrhenius_temp(V_UE_ref, temp, Ea_V_UE, temp_ref);

    // Assign dependent expressions.
    F_S = V_DE * C[4] * C[1] / (K_DE + C[4] + C[1]);
    F_D = V_UE * C[3] * C[2] / (K_UE + C[3] + C[2]);

    // Initiate vector for storing derivatives.
    vector[4] dCdt;

    // Compute derivatives.
    dCdt[1] = i_s + a_MSA * r_M * C[3] - F_S;
    dCdt[2] = i_d + (1 - a_MSA) * r_M * C[3] + F_S + r_L * C[4] - F_D;
    dCdt[3] = u_Q * F_D * - (r_M + r_E) * C[3];
    dCdt[4] = r_E * C[3] - r_L * C[4];

    return dCdt;
  }

  // From https://discourse.mc-stan.org/t/rng-for-truncated-distributions/3122/12.
  real normal_lb_ub_rng(real mu, real sigma, real lb, real ub) {
      real p1 = normal_cdf(lb, mu, sigma);  // cdf with lower bound
      real p2 = normal_cdf(ub, mu, sigma);  // cdf with upper bound
      real u = uniform_rng(p1, p2);
      return (sigma * inv_Phi(u)) + mu;  // inverse cdf 
  }
}

data {
  int<lower=1> state_dim; // Number of state dimensions (4 for AWB).
  int<lower=1> N_t; // Number of observations.
  array[N_t] real<lower=0> ts; // Univariate array of observation time steps.
  array[state_dim] vector<lower=0>[N_t] y; // Multidimensional array of state observations bounded at 0. y in [state_dim, N_t] shape to facilitate likelihood sampling.
  real<lower=0> temp_ref; // Reference temperature for temperature forcing.
  real<lower=0> temp_rise; // Assumed increase in temperature (°C/K) over next 80 years.
  real<lower=0> prior_scale_factor; // Factor multiplying parameter means to obtain prior standard deviations.
  real<lower=0> obs_error_scale; // Observation noise factor multiplying observations of model output x_hat.
  vector<lower=0>[state_dim] x_hat0; // Initial ODE conditions.
  real<lower=0> u_Q_ref_mean;
  real<lower=0> Q_mean;
  real<lower=0> a_MSA_mean;
  real<lower=0> K_DE_mean;
  real<lower=0> K_UE_mean;
  real<lower=0> V_DE_ref_mean;
  real<lower=0> V_UE_ref_mean;
  real<lower=0> Ea_V_DE_mean;
  real<lower=0> Ea_V_UE_mean;
  real<lower=0> r_M_mean;
  real<lower=0> r_E_mean;
  real<lower=0> r_L_mean;
}

transformed data {
  real t0 = 0; // Initial time.
}

parameters {
  real<lower=0> u_Q_ref; // Reference carbon use efficiency.
  real<lower=0> Q; // Carbon use efficiency linear dependence negative slope.
  real<lower=0> a_MSA; // AWB MBC-to-SOC transfer fraction.
  real<lower=0> K_DE; // SOC decomposition K_m.
  real<lower=0> K_UE; // DOC uptake K_m.
  real<lower=0> V_DE_ref; // Reference SOC decomposition V_max.
  real<lower=0> V_UE_ref; // Reference DOC uptake V_max.
  real<lower=0> Ea_V_DE; // SOC V_max activation energy.
  real<lower=0> Ea_V_UE; // DOC V_max activation energy.
  real<lower=0> r_M; // MBC turnover rate.
  real<lower=0> r_E; // Enzyme production rate.
  real<lower=0> r_L; // Enzyme loss rate.
}

transformed parameters {
  array[N_t] vector<lower=0>[state_dim] x_hat_intmd = ode_ckrk(AWB_ECA_ODE, x_hat0, t0, ts, u_Q_ref, Q, a_MSA, K_DE, K_UE, V_DE_ref, V_UE_ref, Ea_V_DE, Ea_V_UE, r_M, r_E, r_L, temp_ref, temp_rise);
  // Transform model output to match observations y in shape, [state_dim, N_t].
  array[state_dim] vector<lower=0>[N_t] x_hat;
  for (i in 1:N_t) {
    for (j in 1:state_dim) {
      x_hat[j, i] = x_hat_intmd[i, j];
    }
  }
}

model {
  u_Q_ref ~ normal(u_Q_ref_mean, u_Q_ref_mean * prior_scale_factor) T[0, 1];
  Q ~ normal(Q_mean, Q_mean * prior_scale_factor) T[0, 0.1];
  a_MSA ~ normal(a_MSA_mean, a_MSA_mean * prior_scale_factor) T[0, 1];
  K_DE ~ normal(K_DE_mean, K_DE_mean * prior_scale_factor) T[0, 5000];
  K_UE ~ normal(K_UE_mean, K_UE_mean * prior_scale_factor) T[0, 1];
  V_DE_ref ~ normal(V_DE_ref_mean, V_DE_ref_mean * prior_scale_factor) T[0, 1];
  V_UE_ref ~ normal(V_UE_ref_mean, V_UE_ref_mean * prior_scale_factor) T[0, 0.1];
  Ea_V_DE ~ normal(Ea_V_DE_mean, Ea_V_DE_mean * prior_scale_factor) T[5, 80];
  Ea_V_UE ~ normal(Ea_V_UE_mean, Ea_V_UE_mean * prior_scale_factor) T[5, 80];
  r_M ~ normal(r_M_mean, r_M_mean * prior_scale_factor) T[0, 0.1];
  r_E ~ normal(r_E_mean, r_E_mean * prior_scale_factor) T[0, 0.1];
  r_L ~ normal(r_L_mean, r_L_mean * prior_scale_factor) T[0, 0.1];

  // Likelihood evaluation.
  for (i in 1:state_dim) {
    y[i,] ~ normal(x_hat[i,], obs_error_scale * x_hat[i,]);
  }
}

generated quantities {
  array[N_t] vector<lower=0>[state_dim] x_hat_post_pred;
  array[state_dim] vector<lower=0>[N_t] x_hat_post_pred_intmd;  
  array[N_t, state_dim] real<lower=0> y_hat_post_pred;
  x_hat_post_pred_intmd = ode_ckrk(AWB_ECA_ODE, x_hat0, t0, ts, u_Q_ref, Q, a_MSA, K_DE, K_UE, V_DE_ref, V_UE_ref, Ea_V_DE, Ea_V_UE, r_M, r_E, r_L, temp_ref, temp_rise);
  // Transform posterior predictive model output to match observations y in dimensions, [state_dim, N_t].
  for (i in 1:N_t) {
    for (j in 1:state_dim) {
      x_hat_post_pred[j, i] = x_hat_post_pred_intmd[i, j];
    }
  }
  for (i in 1:state_dim) {
    y_hat_post_pred[i,] = normal_rng(x_hat_post_pred[i,], obs_error_scale * x_hat_post_pred[i,]);
  }
  real u_Q_ref_post = normal_lb_ub_rng(u_Q_ref_mean, u_Q_ref_mean * prior_scale_factor, 0, 1);
  real Q_post = normal_lb_ub_rng(Q_mean, Q_mean * prior_scale_factor, 0, 0.1);
  real a_MSA_post = normal_lb_ub_rng(a_MSA_mean, a_MSA_mean * prior_scale_factor, 0, 1);
  real K_DE_post = normal_lb_ub_rng(K_DE_mean, K_DE_mean * prior_scale_factor, 0, 5000);
  real K_UE_post = normal_lb_ub_rng(K_UE_mean, K_UE_mean * prior_scale_factor, 0, 1);
  real V_DE_ref_post = normal_lb_ub_rng(V_DE_ref_mean, V_DE_ref_mean * prior_scale_factor, 0, 1);
  real V_UE_ref_post = normal_lb_ub_rng(V_UE_ref_mean, V_UE_ref_mean * prior_scale_factor, 0, 0.1);
  real Ea_V_DE_post = normal_lb_ub_rng(Ea_V_DE_mean, Ea_V_DE_mean * prior_scale_factor, 5, 80);
  real Ea_V_UE_post = normal_lb_ub_rng(Ea_V_UE_mean, Ea_V_UE_mean * prior_scale_factor, 5, 80);
  real r_M_post = normal_lb_ub_rng(r_M_mean, r_M_mean * prior_scale_factor, 0, 0.1);
  real r_E_post = normal_lb_ub_rng(r_E_mean, r_E_mean * prior_scale_factor, 0, 0.1);
  real r_L_post = normal_lb_ub_rng(r_L_mean, r_L_mean * prior_scale_factor, 0, 0.1);
}

Nonetheless, @wds15, if you know if a more efficient way of reshaping the arrays of vectors than

  for (i in 1:N_t) {
    for (j in 1:state_dim) {
      x_hat[j, i] = x_hat_intmd[i, j];
    }
  }

and evaluating the likelihood than

  for (i in 1:state_dim) {
    y[i,] ~ normal(x_hat[i,], obs_error_scale * x_hat[i,]);
  }

let me know. If not, there’s some other issues I have with the Stan model post-compiling, but they aren’t related to multidimensional vectorisation or reshaping (that I know of), so I’ll be creating a separate thread for those issues.