Fit ode parameters using multiple experiments

Modeling
1

Hi,

Thanks for developing this language and making it available!

I am rather new to Stan and using it for my masters thesis.

My project aims to estimate parameters of a 1d model for transient heat in fluid flow in pipes. To fit the model I have data from a number of experiments where the inputs and states (flow rates, fluid type, temperature and more) to the system are varied between experiments.

My question regards how i should go about inferring parameters of the model from multiple experiments.

My initial thought has been to randomly choose an experiment at each iteration in Stan. I believe this might explain it:

p(\theta | y_{1:N}, u_{1:N}) \propto p(\theta) p(y_{1:N}| \theta, u_{1:N}).

Below I’ve added a non-working model with some pseudo-code just to try to illustrate what im getting at.

The number of rows in each dataset vary quite alot, producing a ragged data structure. To get it all into Stan I stack the data vertically and add a vector of indexes to access each experiment in Stan.

functions {
real ode_fun(real t, vector y, theta, ... )
  real dydt = f(t, y, theta);
  return dydt;
}

data {
  int N            // Number of experiments
  int K            // Sum of timesteps in all experiments
  vector[N] N_T    // Number of timesteps for each experiment
  vector[N] N_i    // Indexes to each experiment (cumsum(N_T))
  vector[K] T       // Timesteps for all experiments stacked vertically
  vector[K] Z    // All experimental data stacked vertically
}

parameters {
 real<lower = 0> theta;
}

model {
 theta ~ normal(1, 0.5);
 
 choose n uniformly from <1 : N>
 // Extract time and experimental data
 int n_T       = N_T[n];
 int n_i       = N_i[n];
 vector[n_T] t = T[n_i - n_T + 1 : n_i];
 vector[n_T] z = Z[n_i - n_T + 1 : n_i];
 
 vector[1] y[n_T] = ode_bdf_tol(ode_fun, ... , theta, ... )

 z ~ normal(y, 1);
}

If this is a viable approach, how do i go about randomly selecting an experiment inside the stan model?

My current and loosely defined plan can be summed up like this:

  • Infer posterior for all experiments (bag optimize)
  • Infer posterior for groups of experiments based on similar system states
  • Infer posterior for each individual experiment
  • Use cross validation to validate the MAP estimate on groups of experiments.

I’m happy for any input that can be offered. I struggled a bit presenting this in a concise way while also explaining why i want to do what i do. So if anything is unclear, please point it out.

2

Why do you want to select just one experiment?

Switching experiments during the stan run unlikely to make sense, as this will change the posterior.

3

Why do you want to select just one experiment?

I thought that to get a posterior representative of all the experiments i’d have to do something like what i suggested, or possibly do some kind of prior update based on the posterior of a previously sampled experiment… But this question made me realise that i might as well just sample all the experiments at each iteration, which is probably what you are getting at? I cant believe i didnt think of that earlier, thank you so much!

Switching experiments during the stan run unlikely to make sense, as this will change the posterior.

Thanks for pointing that out, makes sense.

4

I think this is possible in principle, but in practice I don’t know how to (cheaply) “save” the posterior and update it using only the new data.

If all experiments concern the same system with the same ODE/PDE parameters, just different initial/boundary conditions or forcing functions, then the right approach should be to fit all experiments at once.

Of course this may be computationally expensive or difficult to fit. How do you discretize the problem?

Do report how everything works out, I’m curious.

5

Of course this may be computationally expensive or difficult to fit. How do you discretize the problem?

The problem is separated into segments of pipe with equal flow rate, meaning a segment starts and ends at an inflow location. Each segment is discretized using FEM on a 1d mesh. The theoretical foundation and implementation of this is done by others, whereas I’m trying to improve the model using data.

Computationally a single forward run of ~200 meters of pipe is in seconds to tens of seconds. So worst case runtime for 1000 warmups and 1000 samples is in days.

Good to see what I’m working on is of interest. I’m happy to elaborate and to give an update once i get some results.

6

Hi,

Are your experiments already available, or are you updating your parameter inference as experiments come in?

If it’s the former, what I do is fit the ode model to each experiment. That is in your model section, you fit the same parameter(s) to several observation vectors.

If it’s the latter, i don’t know. I would really like to learn how to use the posterior of after inferring from one experiment as the prior for a next. Especially if there is some parameters covary.

Good luck

7

Keep in mind that each hmc iteration consists of several solves, up to 1023 I believe with default settings. Also, runtime may vary considerably depending on the parameter values. In addition, the ODE solve with sensitivities generally takes longer than a standard solve without sensitivities.

8

They are available, but varying in length, so not really suitable to store as vectors in a matrix.

I think i can overcome this either by zero-padding everything to equal length or appending the simulated results to a single vector and fit against the concatenated measurements of equal length, like this:

model {
  theta ~ normal(mu, sigma); 
  vector[length of all experiments] sim_results;
  for (exp in experiments) {
    sim_results[exp.start : exp.end]  = simulate(..., theta);
  }
  experiment_data ~ normal(sim_results, 1);
}

Which i think should be the same as if I had matrices of equal length measurements and experiments to fit.

The data is preprocessed so that experiments start and end at 0 so this wont introduce any “jumps” in the data, though as long as the “jumps” dont introduce any new error I dont think this should matter anyway.

Thanks!

9

Thanks, this is good to know!

10

Reading the set-up of your experiments, it sounds to me similar to a model published about microbial dynamics in Stan. In that paper they also had time series with varying length. Check the supplementary material and you will find how it is implemented.

11

Hi,

You don’t need to use a matrix, you van use a list of vectors of varying lengths. What you wrote in your last post is similar to what I’ve done, which works. I see no reason for padding. Another advantage of the list of vectors approach, is that observations do not need to de equally spaced in different experiments.

Best,

Sergio

12

This is great stuff, thanks!

13

I think i see what you’re getting at now! I just didnt know it was possible to have a lis of vectors of different lengths in Stan, will try to implement this as it appears to be the approach that makes the most sense.

Thanks!