where the normalization constant Z = f dxp*(x) is not known but p*(x) can be evaluated at any particular x . Suppose we want to estimate the expectation of a function f(x) under the distribution p ( x ) denoted as

$$\langle f(\mathbf{x})\rangle_{p(\mathbf{x})}=\int d x f(\mathbf{x})p(\mathbf{x})$$

e.g., the mean of x under p(x) is given by (x). The intractable integration can be approximated by an average if we can find N points x (i), i = 1 . . . N from p (x)

$$\langle f(\mathbf{x})\rangle_{p(\mathbf{x})}\approx{\frac{1}{N}}\sum_{i\,=\,1}^{N}f(\mathbf{x}^{(i)})$$
$$\left(\,14\,\right)$$

When x (i) are generated by independently sampling from p (x), it can be shown that as N
approaches infinity, the approximation becomes exact.

However, generating independent samples from p ( x ) is a difficult task in high dimensions but it is usually easier to generate dependent samples, that is we generate x (i+1) by making use of x (i) . It is somewhat surprising, that even if x (i) and x (i+1) are correlated (and provided ergodicity conditions are satisfied), Eq. 14 remains still valid and estimated quantities converge to their true values when number of samples N goes to infinity.

A sequence of dependent samples x (i) is generated by using a Markov chain that has the stationary distribution p (x). The chain is defined by a collection of transition probabilities, i.e., a transition kernel T(x(i+l)|x(i) ). The definition of the kernel is implicit, in the sense that one defines a procedure to generate the x (i+1) given x (i). The Metropolis algorithm (Metropolis & Ulam, 1949; Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953) provides a simple way of defining an ergodic kernel that has the desired stationary distribution p (x). Suppose we have a sample x (i). A candidate x ′ is generated by sampling from a symmetric proposal distribution q(x'|x(i) (for example a Gaussian centered at x ( i ). The candidate x ′ is accepted as the next sample x ( i + 1 ) if p(x ′ ) > p(x ( i ). If x has a lower probability, it can be still accepted, but only with probability p(x')/p(x'').

The algorithm is initialized by generating the first sample x (0) according to an (arbitrary) proposal distribution.

However for a given transition kernel T, it is hard to assess the time required to converge to the stationary distribution so in practice one has to run the simulation until a very large nmber of samples have been obtained, (see e.g., Roberts & Rosenthal, 1998). The choice of the proposal distribution q is also very critical. A poor choice may lead to the rejection of many candidates x' hence resulting in a very slow convergence to the stationary distribution.

For a large class of probability models, where the full posterior p ( x ) is intractable, one can still efficiently compute marginals of form p ( x k | x − k ), x − k = x 1 . . . x k − 1 , x k + 1 , . . x K
exactly. In this case one can apply a more specialized Markov chain Monte Carlo (MCMC)
algorithm, the Gibbs sampler given below.

 ## 1. Initialize $x_{1:K}^{[0\:]}$ by sa. 
1. Initialize x (0)

 $\text{2.For}i\,=0\,.\,.\,.\,N-\,1$ . 
$$\mathrm{g~from~a~proposal~}q(x_{1:K})$$