Skip to content

Latest commit

 

History

History
110 lines (68 loc) · 4.13 KB

README.md

File metadata and controls

110 lines (68 loc) · 4.13 KB
Raw

Multi-Parameter Bayesian Inference Using Markov Chain Monte Carlo (MCMC) Sampling and the Metropolis-Hastings Algorithm

Features

  • The code has been written in plain vanilla C++ and tested using g++ 8.1.0 (MinGW-W64).

  • Likelihood (pdf) can be defined as an arbitrary function with any number of independent parameters.

  • Prior functions are defined using an array of structures, and can be any pdf from function prior_dist in file Metropolis.cpp (other priors can be easily added).

  • Jumps in the Metropolis-Hastings algorithm are proposed using a normal distribution of the parameters.

  • Function random_number in file Metropolis.cpp can be used to generate random numbers from any arbitrary pdf.

  • Results are passed back using a structure.

  • Usage: test.exe example.

Main Parameters

example Name of the example to run.

likelihood Name of the likelihood function.

par Array with the parameters of the likelihood function.

n_data Number of data to be sampled from the likelihood function.

data Array with the data sampled from the likelihood function.

a0, b0 Support interval for the likelihood function.

priors Array of structures with the priors. Each prior is assigned to one of the likelihood parameter following the same order as in par.

samples Number of jumps to perform in the Metropolis-Hastings algorithm.

par_init Initial value for the parameters in the Metropolis-Hastings algorithm.

width_prop Standard deviation of the normal distribution used to search the neighboroud of a parameter. A good value should give about 50% of accepted jumps.

i0 Index specifying the burn-in/warm-up amount.

posterior Array containing the jumps (accepted or rejected) of all parameters.

jumps Number of jumps actually accepted.

res Structure with the results returned by the Metropolis-Hastings algorithm.

save_res Name of the file where to save the results. If empty, results are not saved.

Examples

There are three examples: Coin, Normal, and Random (see test.cpp for the specific parameters and results). A brief description is given below. Plots similar to the ones in the Python version (see the Readme) can be easily generated using the results saved in the file specified by save_res.

Coin:

One parameter (theta), Bernoulli distribution as likelihood, beta distribution as prior, admit an analytical solution.

// Numerical results:
// - accepted jumps = 53.9%
// - <theta> mean and std = 0.332, 0.057

// Analytical results:
// - <theta> mean and std = 0.333, 0.058

Normal:

Two parameters (mean and standard deviation), normal distribution as likelihood, normal distribution as prior for the mean, gamma distribution as prior for the standard deviation.

// Numerical results:
// - accepted jumps = 52.7%
// - <mu> mean and std = -1.207, 0.170
// - <sigma> mean and std = 0.728, 0.150

// Numerical (pymc3) results:
// - <mu> mean and std = -1.210, 0.171
// - <sigma> mean and std = 0.730, 0.156

The pymc3 results have been obtained loading the same random data generated by the C++ version (array data) in the Python version.

Random:

A generic pdf (a piece-wise function in the example) is correctly approximated using 50000 randomly generated numbers.

// Piece-wise function:
//
// pdf =  0.0                        in [-inf, 0]
// pdf =  0.3 * x                    in (0, 1]
// pdf = -0.2 * x + 0.5              in (1, 2]
// pdf =  0.1                        in (2, 3]
// pdf =  0.1 * x - 0.2              in (3, 4]
// pdf =  0.2                        in (4, 5]
// pdf = -0.1 * x + 0.7              in (5, 7]
// pdf =  0.0                        in (7, +inf]
//
// The integral in [-inf, +inf] is 1.

References