Bayes Inference

Bayesian inference is a statistical method for making probabilistic inference about unknown parameters based on observed data and prior knowledge. The update from prior knowledge to posterior in light of observed data is based on Bayes’ theorem:

(1)\[p(\mathbf{x} \mid \mathbf{d}) = \frac{L(\mathbf{x} \mid \mathbf{d}) p(\mathbf{x})} {\int L(\mathbf{x} \mid \mathbf{d}) p(\mathbf{x}) d \mathbf{x}}\]

where \(\mathbf{x}\) represents the collection of unknown parameters and \(\mathbf{d}\) represents the collection of observed data. The prior probability distribution, \(p(\mathbf{x})\), represents the degree of belief in the parameters before any data is observed. The likelihood function, \(L(\mathbf{x} \mid \mathbf{d})\), represents the probability of observing the data given the parameters. The posterior distribution, \(p(\mathbf{x} \mid \mathbf{d})\), is obtained by multiplying the prior probability distribution by the likelihood function and then normalizing the result. Bayesian inference allows for incorporating subjective prior beliefs, which can be updated as new data becomes available.

For many real world problems, it is hardly possible to analytically compute the posterior due to the complexity of the denominator in equation (1), namely the nomalizing constant. In this module, two numerical approximations are implemented: grid estimation and Metropolis Hastings estimation.

In grid estimation, the denominator in equation (1) is approximated by numerical integration on a regular grid and the posterior value at each grid point is computed, as shown in equation (2). The number of grid points increases dramatically with the increase of the number of unknown parameters. Grid estimation is therefore limited to low-dimensional problems.

(2)\[p(\mathbf{x} \mid \mathbf{d}) \approx \frac{L(\mathbf{x} \mid \mathbf{d}) p(\mathbf{x})} {\sum_{i=1}^N L\left(\mathbf{x}_i \mid \mathbf{d}\right) p\left(\mathbf{x}_i\right) \Delta \mathbf{x}_i}\]

Metropolis Hastings estimation directly draw samples from the unnormalized posterior distribution, namely the numerator of equation (1). The samples are then used to estimate properties of the posterior distribution, like the mean and variance, or to estimate the posterior distribution.

GridEstimation Class

The GridEstimation class is imported by:

from psimpy.inference.bayes_inference import GridEstimation

Methods

class GridEstimation(ndim, bounds=None, prior=None, likelihood=None, ln_prior=None, ln_likelihood=None, ln_pxl=None, args_prior=None, kwgs_prior=None, args_likelihood=None, kwgs_likelihood=None, args_ln_pxl=None, kwgs_ln_pxl=None)[source]

Set up basic input for Bayesian inference.

Parameters
  • ndim (int) – Dimension of parameter x.

  • bounds (numpy array) – Upper and lower boundaries of each parameter. Shape (ndim, 2). bounds[:, 0] corresponds to lower boundaries of each parameter and bounds[:, 1] to upper boundaries of each parameter.

  • prior (Callable) – Prior probability density function. Call with prior(x, *args_prior, **kwgs_prior). Return the prior probability density value at x, where x is a one dimension numpy.ndarray of shape (ndim,).

  • likelihood (Callable) – Likelihood function. Call with likelihood(x, *args_likelihood, **kwgs_likelihood). Return the likelihood value at x.

  • ln_prior (Callable) – Natural logarithm of prior probability density function. Call with ln_prior(x, *args_prior, **kwgs_prior). Return the value of natural logarithm of prior probability density at x.

  • ln_likelihood (Callable) – Natural logarithm of likelihood function. Call with ln_likelihood(x, *args_likelihood, **kwgs_likelihood). Return the value of natural logarithm of likelihood at x.

  • ln_pxl (Callable) – Natural logarithm of the product of prior and likelihood. Call with ln_pxl(x, *args_ln_pxl, **kwgs_ln_pxl). Return the value of natural logarithm of the product of prior and likelihood at x.

  • args_prior (list, optional) – Positional arguments for prior or ln_prior.

  • kwgs_prior (dict, optional) – Keyword arguments for prior or ln_prior.

  • args_likelihood (list, optional) – Positional arguments for likelihood or ln_likelihood.

  • kwgs_likelihood (dict, optional) – Keyword arguments for likelihood or ln_likelihood.

  • args_ln_pxl (list, optional) – Positional arguments for ln_pxl.

  • kwgs_ln_pxl (dict, optional) – Keyword arguments for ln_pxl.

run(nbins)[source]

Use Grid approximation to estimate the posterior.

Parameters

nbins (int or list of ints) – Number of bins for each parameter. If int, the same value is used for each parameter. If list of int, it should be of length ndim.

Return type

tuple[ndarray, list[ndarray]]

Returns

  • posterior (numpy array) – Estimated posterior probability density values at grid points. Shape of (nbins[0], nbins[1], ..., nbins[ndim]).

  • x_ndim (list of numpy array) – Contain ndim 1d numpy.ndarray x1, x2, … Each xi is a 1d array of length nbins[i], representing the 1d coordinate array along the i-th axis.

MetropolisHastingsEstimation Class

The MetropolisHastingsEstimation class is imported by:

from psimpy.inference.bayes_inference import MetropolisHastingsEstimation

Methods

class MetropolisHastingsEstimation(ndim, bounds=None, prior=None, likelihood=None, ln_prior=None, ln_likelihood=None, ln_pxl=None, args_prior=None, kwgs_prior=None, args_likelihood=None, kwgs_likelihood=None, args_ln_pxl=None, kwgs_ln_pxl=None)[source]

Set up basic input for Bayesian inference.

Parameters
  • ndim (int) – Dimension of parameter x.

  • bounds (numpy array) – Upper and lower boundaries of each parameter. Shape (ndim, 2). bounds[:, 0] corresponds to lower boundaries of each parameter and bounds[:, 1] to upper boundaries of each parameter.

  • prior (Callable) – Prior probability density function. Call with prior(x, *args_prior, **kwgs_prior). Return the prior probability density value at x, where x is a one dimension numpy.ndarray of shape (ndim,).

  • likelihood (Callable) – Likelihood function. Call with likelihood(x, *args_likelihood, **kwgs_likelihood). Return the likelihood value at x.

  • ln_prior (Callable) – Natural logarithm of prior probability density function. Call with ln_prior(x, *args_prior, **kwgs_prior). Return the value of natural logarithm of prior probability density at x.

  • ln_likelihood (Callable) – Natural logarithm of likelihood function. Call with ln_likelihood(x, *args_likelihood, **kwgs_likelihood). Return the value of natural logarithm of likelihood at x.

  • ln_pxl (Callable) – Natural logarithm of the product of prior and likelihood. Call with ln_pxl(x, *args_ln_pxl, **kwgs_ln_pxl). Return the value of natural logarithm of the product of prior and likelihood at x.

  • args_prior (list, optional) – Positional arguments for prior or ln_prior.

  • kwgs_prior (dict, optional) – Keyword arguments for prior or ln_prior.

  • args_likelihood (list, optional) – Positional arguments for likelihood or ln_likelihood.

  • kwgs_likelihood (dict, optional) – Keyword arguments for likelihood or ln_likelihood.

  • args_ln_pxl (list, optional) – Positional arguments for ln_pxl.

  • kwgs_ln_pxl (dict, optional) – Keyword arguments for ln_pxl.

run(nsamples, mh_sampler)[source]

Use Metropolis Hastings sampling to draw samples from the posterior.

Parameters
  • nsamples (int) – Number of samples to be drawn.

  • mh_sampler (instance of MetropolisHastings) – Metropolis Hastings sampler.

Return type

tuple[ndarray, ndarray]

Returns

  • mh_samples (numpy array) – Samples drawn from the posterior. Shape of (nsamples, ndim).

  • mh_accept (numpy array) – Shape of (nsamples,). Each element indicates whether the corresponding sample is the proposed new state (value 1) or the old state (value 0). np.mean(mh_accept) thus gives the overall acceptance ratio.