Bayesian Parameter Estimation for SimBiology

Bayesian parameter estimation to fit SimBiology models to data using a constant (Gaussian) error model.

Bayesian Parameter Estimation

Copyright 2023 The MathWorks, Inc.

Bayesian parameter estimation to fit SimBiology models to data using a constant (Gaussian) error model. Create a mcmc.MetropolisHastingsSampler with a SimBiology model, data, prior distribution of parameters to estimate, as well as a proposal distribution. Then generate samples and compute the maximum aposterior parameter estimates using the sample method.

Getting Started

Read and run examples driver_PK.mlx and driver_PCSK9.mlx.

Help

After installing mcmc.MetropolisHastingsSampler as MATLAB® Add-On, you can also type

>> help mcmc.MetropolisHastingsSampler

>> help mcmc.MetropolisHastingsSampler.sample

>> help mcmc.MetropolisHastingsSampler.plot

in the MATLAB® Command Window to get instructions for of how to use the sampler.

Software Requirements

This feature is supported in MATLAB® versions R2022a and later. Required products are SimBiology® and Statistics And Machine Learning Toolbox™.

Limitation

• mcmc.MetropolisHastingsSampler assumes that response and data names do not contain =.

Markov-Chain Monte Carlo Sampling

The goal is to sample from the posterior distribution of parameters to estimate given a SimBiology model and data. We assume only a prior distribution and the likelihood function is known. According to Bayes' theorem, the posterior distribution is then proportional to the product of the prior and the likelihood:

posterior(sample) ~ likelihood(sample | data) * prior(sample)

Metropolis Algorithm

We use a Markov-Chain Monte Carlo (MCMC) algorithm, proposed by N. Metropolis [1], to draw samples from the posterior distribution. The Metropolis algorithm to draw N samples can be summarized as follows:

Start at initial sample S with posterior(S) > 0.

1. Draw new sample N from proposal distribution Proposal.
2. Draw value U from uniform distribution U[0,1].
3. If min(1, posterior(N)/posterior(S)) > U % Accept sample: add N to set of samples. S = N;end
4. Samples{end+1} = S;

Repeat until a specified number of samples are generated.

The proposal is a joint probability distribution used to "propose" new samples. The distribution need not be symmetric. The proposal is symmetric if Proposal(N | S) == Proposal(S | N), i.e. the probability of proposing the sample N when the current sample of the Markov chain is S is equal to the probability of proposing the sample S if the current sample of the Markov chain was N.

Illustration

The central figure shows how the Metropolis algorithm draws samples from an artificial joint probability distribution of two parameters. The yellow star is the current sample, the green/red star indicates the new sample. The color green inidcates that the new sample is accepted, red indicates that the new sample is rejected.

The two Markov chains of the individual parameters are shown to the right and bottom of the central plot. One sample is added per iteration of the Metropolis algorithm. The outermost plots show the marginal distributions of the two parameters.

References

[1] Metropolis, N., Rosenbluth, A.W., Rosenbluth, M., Teller, A., Teller, E. (1953). "Equation of State Calculations by Fast Computing Machines". Journal of Chemical Physics. 21 (6): 1087–1092.