Application to system with multiple parameters, inputs and outputs

[Zhe-Jessica]

Brief of the question: For a system with multiple parameters, inputs and outputs, what is the format of the argument values in the functions SimulatorForwardProblem () and SimulatorInferenceProblem ()?

More details of the question below:
I tried the code in the "examples" file and had a general idea of how to use the package. The "linearode" model is pretty straightforward, but I could not run the examples of "Tsurf_inv" and "ddm". I think the primary reason is I can't load the data. When I tried the SimulationBasedInference.jl, I used the "linearode" as a prototype. There are some errors when I tried the package to a simple lotka volterra model (with 1 parameters and two outputs).

1. In the SimulatorForwardProblem(odeprob, observable)
   Can I give a function which solve the odeprob as an argument to the function?

If I have multiple outputs, how to use SimulatorObservable () function for the argument "observable"? Get the argument for each output individually or together in oneSimulatorObservable ()function? What format I should follow to set the arguement of all outputs?

2. In the function SimulatorInferenceProblem (forward_prob, ode_solver, model_prior, lik)
   Similar questions of how to set "model_prior" and "lik" for multiple parameters and likelihood.

3. In my case, I did not get any error before I run "enis_sol = solve(inference_prob, EnIS(), ensemble_size=128, rng=rng);" and get error and stuck when I run the solve () function.
   The error shows:

I think the error point out that I probably defined wrong argument types in functions for inference_prob. It is the main reason why I asked previous two questions.

It might be difficult to check more details behind the problem, so I attached the code if you would like to have a look.
https://github.com/Zhe-Jessica/LV_SBI/blob/main/LV_SBI.jl

I do understand you are probably busy and do not have time to check all the code. It would be very helpful if there are some other examples of the application of the SimulationBasedInference.jl package to the systems with multiple parameters and multiple outputs. I probably can have more ideas from this kind of examples.

Thank you very much for your time. Look forward to your response.

[bgroenks96]

Hi, thanks for raising the issue.

1. Can I give a function which solve the odeprob as an argument to the function?

You can either provide a function which takes the parameters (i.e. f(theta)) and returns the relevant state/outputs, or you can provide an ODEProblem or similar SciMLProblem and provide the solver at solve-time in the typical fashion for the SciML interface, e.g:

forward_prob = SimulatorForwardProblem(odeprob, observable)
sol = solve(forward_prob, Tsit5())

When constructing a SimulatorInferenceProblem the solver must also be provided.

2. If I have multiple outputs, how to use SimulatorObservable () function for the argument "observable"? Get the argument for each output individually or together in oneSimulatorObservable ()function? What format I should follow to set the arguement of all outputs?

This depends somewhat on your use case. For the Lotka-Volterra example, I would probably just make one observable that returns both outputs. You could then add more observables that compute relevant derived quantities (if applicable). Multiple parameter priors can be provided simply as arguments to the prior method, e.g:

model_prior = prior(a=Beta(1,1), b=Beta(1,1), ...)

Note that the order must match that of the parameter vector in the ODEProblem.

3. In my case, I did not get any error before I run "enis_sol = solve(inference_prob, EnIS(), ensemble_size=128, rng=rng);" and get error and stuck when I run the solve () function.

The error you showed here arises because the SimulatorForwardProblem assumes that your parameters are a vector, not a scalar. So this can be resolved by changing α = p to α = p[1].

There are some other issues with your configuration of the observables, however. Since you want a time-sampled observable, you need to provide the initial time and save times as such:

y1 = SimulatorObservable(:u1, state -> state.u[1,:], tspan[1], tsave, (1,))
y2 = SimulatorObservable(:u2, state -> state.u[2,:], tspan[1], tsave, (1,))

After fixing that, the code should work.

I will leave this issue open for now in case you have further questions. There are also clearly some improvements that need to be made to the documentation to make these things more clear.

If you can spare a bit of time, I would also appreciate if you could open a pull request adding this LV ODE example to the examples/ folder.

[Zhe-Jessica]

Thank you very much for your help. You pointed out the problem of the code. The parameter was not given as a vector, which leads to the problem.

Besides, the arguments of SimulatorObservable () as you mentioned would lead to errors.

y1 = SimulatorObservable(:u1, state -> state.u[1,:], tspan[1], tsave, (1,))
y2 = SimulatorObservable(:u2, state -> state.u[2,:], tspan[1], tsave, (1,))

I managed to get a result by using SimulatorObservable(:u2, state -> state.u[2,:], (Nt,)), where Nt is the length of the time step.

I will add the LV ODE example to the folder. It is a simple time-related example to show the application to the system with multiple parameters and outputs. Hope it can help people apply to similar systems.

The last question I have for the package is about the noise_scale_prior (for example in linearode.jl example). Could you please explain more details of the distribution of this prior? Does this prior describle the noise range (amplitude), or the noise ratio or the standard deviation (such as if the noise follow normal distribution)?

[bgroenks96]

Could you please explain more details of the distribution of this prior? Does this prior describle the noise range (amplitude), or the noise ratio or the standard deviation (such as if the noise follow normal distribution)?

It is the prior on the noise scale parameter for the likelihood/predictive distribution. For a Gaussian likelihood, this corresponds to the standard deviation. Note that for some algorithms, this prior is (partially) ignored. For example, the ensemble Kalman methods (EKS and ESMDA) necessarily assume this noise term to be constant, so the expected value of the prior is taken as the constant.

I managed to get a result by using SimulatorObservable(:u2, state -> state.u[2,:], (Nt,)), where Nt is the length of the time step.

This is very surprising. I am curious to see the working code, as I don't think that should work unless you were to construct the SimulatorForwardProblem using a function (which solves the ODEProblem) instead of supplying the ODEProblem directly.

[Zhe-Jessica]

Thanks for your explain. I used a function to solve the ODEProblem in the working code.

I am also thinking about another question. I know PosteriorStats.summarize shows the mean, std, etc. Meanwhile, I also want to check the estimation of paramter from each particle/ensemble member, so that I can see more detail of the estimation. It is possible to check the estimation of each particle using enis_sol.result.ens in ensemble importance sampling method. However, I couldn't find similar result for ensemble smoother and ensemble Kalman sampler method. Could you please give some suggestions on how to check them?

Besides, I pulled a request to add the example last week, but it seems didn't work. I will also add a link in case you would like to check or add to the repository.

[bgroenks96]

The function you want is get_transformed_ensemble; see the degree-day model example here.

The importance sampling algorithm is actually a bit special in this case. The ensemble stored in enis_sol.result.ens is actually the prior ensemble, i.e. the samples from the prior. The output of importance sampling is a set of weights (determined by the likelihood) which can then be used in statistical calculations like mean and std. You can find an example of this here.

In general, I think all of the APIs related to handling the results need further improvement and documentation for user friendliness... but it's a tricky problem since each class of algorithm has very different outputs. Feel free to create more issues if you encounter further problems.

[bgroenks96]

Actually, PosteriorStats.summarize for EnIS is broken... it doesn't account for the weights. I updated it to throw an error for now.

You can find the updated example script here which uses the recommended pattern for the observables (i.e. treating them as time sampled).