Pseudo-marginal Metropolis–Hastings algorithm

Given a current state ${\displaystyle \theta _{n}}$ the Metropolis–Hastings algorithm proposes a new state according to some density ${\displaystyle \theta '\sim Q(\cdot \mid \theta _{n})}$. The algorithm then sets ${\displaystyle \theta _{n+1}=\theta '}$ with probability

${\displaystyle a(\theta _{n},\theta ')=\min \left(1,{\frac {\pi (\theta ')}{\pi (\theta _{n})}}{\frac {Q(\theta _{n}\mid \theta ')}{Q(\theta '\mid \theta _{n})}}\right)}$

otherwise the old state is kept, that is, ${\displaystyle \theta _{n+1}=\theta _{n}}$.

If the density ${\displaystyle \pi }$ is not available analytically the above algorithm cannot be employed. The pseudo-marginal Metropolis–Hastings algorithm in contrast only assumes the existence of an unbiased estimator ${\displaystyle {\hat {\pi }}_{\theta }}$, i.e. the estimator must satisfy the equation ${\displaystyle \mathbb {E} [{\hat {\pi }}_{\theta }]=\pi (\theta ).}$ Now, given ${\displaystyle \theta _{n}}$ and the respective estimate ${\displaystyle {\hat {\pi }}_{\theta _{n}}}$ the algorithm proposes a new state according to some density ${\displaystyle \theta '\sim Q(\cdot \mid \theta _{n})}$. Next, compute an estimate ${\displaystyle {\hat {\pi }}_{\theta '}}$ and set ${\displaystyle \theta _{n+1}=\theta '}$ with probability

${\displaystyle a(\theta _{n},\theta ')=\min \left(1,{\frac {{\hat {\pi }}_{\theta '}}{{\hat {\pi }}_{\theta _{n}}}}{\frac {Q(\theta _{n}\mid \theta ')}{Q(\theta '\mid \theta _{n})}}\right)}$

otherwise the old state is kept, that is, ${\displaystyle \theta _{n+1}=\theta _{n}}$.

Consider a model consisting of i.i.d. latent real-valued random variables ${\displaystyle Z_{1},\ldots ,Z_{n}}$ with ${\displaystyle Z_{i}\sim f_{\theta }(\cdot )}$ and suppose one can only observe these variables through some additional noise ${\displaystyle Y_{i}\mid Z_{i}=z\sim g_{\theta }(\cdot \mid z)}$ for some conditional density ${\displaystyle g}$. (This could be due to measurement error, for instance.) We are interested in Bayesian analysis of this model based on some observed data ${\displaystyle y_{1},\ldots ,y_{n}}$. Therefore, we introduce some prior distribution ${\displaystyle p(\theta )}$ on the parameter. In order to compute the posterior distribution

${\displaystyle p(\theta \mid y_{1},\ldots ,y_{n})\propto p_{\theta }(y_{1},\ldots ,y_{n})p(\theta )}$

we need to find the likelihood function ${\displaystyle p_{\theta }(y_{1},\ldots ,y_{n})}$. The likelihood contribution of any observed data point ${\displaystyle y}$ is then

${\displaystyle p_{\theta }(y)=\int g_{\theta }(y\mid z)f_{\theta }(z)\,dz}$

and the joint likelihood of the observed data ${\displaystyle y_{1},\ldots ,y_{n}}$ is

${\displaystyle p_{\theta }(y_{1},\ldots ,y_{n})=\prod _{i=1}^{n}p_{\theta }(y_{i})=\prod _{i=1}^{n}\int g_{\theta }(y_{i}\mid z_{i})f_{\theta }(z_{i})\,dz_{i}.}$

If the integral on the right-hand side is not analytically available, importance sampling can be used to estimate the likelihood. Introduce an auxiliary distribution ${\displaystyle q}$ such that ${\displaystyle g_{\theta }(y\mid z)f_{\theta }(z)>0\Rightarrow q(z)>0}$ for all ${\displaystyle z}$ then

${\displaystyle {\hat {p}}_{\theta }(y_{i})={\frac {1}{N}}\sum _{k=1}^{N}{\frac {g_{\theta }(y_{i}\mid Z_{k})f_{\theta }(Z_{k})}{q(Z_{k})}},\qquad Z_{k}{\overset {i.i.d.}{\sim }}q(\cdot )}$

is an unbiased estimator of ${\displaystyle p_{\theta }(y_{i})}$ and the joint likelihood can be estimated unbiasedly by

${\displaystyle {\hat {p}}_{\theta }(y_{1},\ldots ,y_{n})=\prod _{i=1}^{n}{\hat {p}}_{\theta }(y_{i})=\prod _{i=1}^{n}{\frac {1}{N}}\sum _{k=1}^{N}{\frac {g_{\theta }(y_{i}\mid Z_{i,k})f_{\theta }(Z_{i,k})}{q(Z_{i,k})}},\qquad Z_{i,k}{\overset {i.i.d.}{\sim }}q(\cdot ).}$