Note
Click here to download the full example code
Bayesian inference
This example shows how to perform Bayesian inference given the uniform prior
\(p(\mathbf{x})=p(x_1,x_2)=0.01\)
where \(x_i \in [-5,5], i=1,2\), and likelihood
\(L(\mathbf{x}|\mathbf{d})=\exp \left(-\frac{1}{100}\left(x_1-1\right)^2-\left(x_1^2-x_2\right)^2\right)\).
import numpy as np
ndim = 2
bounds = np.array([[-5,5],[-5,5]])
def prior(x):
return 0.01
def likelihood(x):
return np.exp(-(x[0]-1)**2/100 - (x[0]**2-x[1])**2)
To estimate the posterior using grid estimation, we need to import the
GridEstimation
class, create an instance, and call the
GridEstimation.run()
method.
from psimpy.inference.bayes_inference import GridEstimation
grid_estimator = GridEstimation(ndim, bounds, prior, likelihood)
posterior, x_ndim = grid_estimator.run(nbins=[50,40])
The following figure plots the estimated posterior.
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1,figsize=(6,4))
# mask insignificant values
posterior = np.where(posterior < 1e-10, None, posterior)
contour = ax.contour(x_ndim[0], x_ndim[1], np.transpose(posterior), levels=10)
plt.colorbar(contour, ax=ax)
ax.set_xlim(bounds[0,0], bounds[0,1])
ax.set_ylim(bounds[1,0], bounds[1,1])
ax.set_title('Grid estimation')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
plt.tight_layout()

To estimate the posterior using Metropolis Hastings estimation, we need to import
the MetropolisHastingsEstimation
class, create an instance, and call the
MetropolisHastingsEstimation.run()
method. The
MetropolisHastingsEstimation.run()
method has a parameter, mh_sampler
,
which takes an instance of MetropolisHastings
as argument.
from psimpy.inference.bayes_inference import MetropolisHastingsEstimation
mh_estimator = MetropolisHastingsEstimation(ndim, bounds, prior, likelihood)
# create a mh_sampler
from psimpy.sampler.metropolis_hastings import MetropolisHastings
from scipy.stats import multivariate_normal
init_state = np.array([-4,-4])
f_sample = multivariate_normal.rvs
nburn = 100
nthin = 10
seed = 1
kwgs_f_sample = {'random_state': np.random.default_rng(seed)}
mh_sampler = MetropolisHastings(ndim=ndim, init_state=init_state,
f_sample=f_sample, bounds=bounds, nburn=nburn, nthin=nthin, seed=seed,
kwgs_f_sample=kwgs_f_sample)
nsamples = 5000
mh_samples, mh_accept = mh_estimator.run(nsamples, mh_sampler)
The following figure plots the samples drawn from the unnormalized posterior, which can be used to estimate the posterior and its poperties.
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,1,figsize=(5,4))
ax.scatter(mh_samples[:,0], mh_samples[:,1], s=10, c='r', marker='o', alpha=0.1)
ax.set_xlim(bounds[0,0], bounds[0,1])
ax.set_ylim(bounds[1,0], bounds[1,1])
ax.set_title('MH estimation')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
plt.tight_layout()

Note
Besides prior
and likelihood
, one can also instantiate
the MetropolisHastingsEstimation
class with
ln_prior
andln_likelihood
: Natural logarithm ofprior
andlikelihood
.
ln_pxl
: Natural logarithm of the product ofprior
andlikelihood
.
Total running time of the script: ( 0 minutes 9.739 seconds)