metropolis#

View page source

Metropolis MCMC Algorithm#

Syntax#

( a , c ) = dismod_at.metropolis ( log_f , m , x0 , s )

log_f#

Given a numpy float vector of length n , the syntax

d = log_f ( x )

sets the float d to the log of the un-normalized density corresponding to the positive function \(f(x)\) mapping \(\B{R}^n\) to the non-negative real values. If \(f(x)\) is zero, the corresponding log-density value should equal - float("inf") .

m#

is the int number of vectors in the MCMC chain.

x0#

is a numpy float vector of length n that specifies the initial vector in the chain (denoted by \(x^0\) below).

s#

is a float or, a numpy float vector of length n , that specifies the scaling for each of the components of \(x\).

Vector Case#

If s is a vector, for \(i = 1 , \ldots, m-1\), and \(j = 0 , \ldots, n-1\), the j-th component of the i-th proposal vector \(y^i\) is given by

\[y_j^i = x_j^{i-1} + w_j^{i-1} s_j\]

where \(w_j^i \sim \B{N}(0, 1)\) are all independent.

Float Case#

If s is a float ,

\[y_j^i = x_j^{i-1} + w_j^{i-1} s\]

a#

is the int acceptance count; i.e. the number of indices \(i\) such that \(x^i = y^i\) (for the other indices \(x^i = x^{i-1}\)).

c#

is an \(m \times n\) numpy float array that contains the components of the Markov Chain. We use the notation \(x_j^i\) for c [ i , j ] . For any smooth function \(g : \B{R}^n \rightarrow \B{R}\), the Metropolis algorithm provides the following approximation as \(m \rightarrow \infty\),

\[\frac{1}{m} \sum_{i=0}^{m-1} g( x^i ) \rightarrow \frac{ \int g( x ) f ( x ) \B{d} x }{ \int f( x ) \B{d} x }\]

Example#

The file user_metropolis.py contains an example and test of this routine.