Metropolis-Hasting Algorithm

One efficient way to get samples from complex distribution is by using Markov Chain. Metropolis-Hasting Algorithm is a famous one for this purpose.

The candidate point for the next sample is accepted as the next state of the chain with the probability given by

\[\alpha (X_t,Y) = \min \left\{ {1,\frac{\pi (Y)q(X_t|Y)}{\pi (X_t)q(Y|X_t)} \right\}\]
If the point $Y$ is not accepted, the chain does not move and  $X_{t+1}=X_t$

$q(.|X_t)$ is a proposal distribution. It can be multivariate normal with mean $X_t$ and fixed covariance matrix.

For example, generate random variables from a standard Cauchy distribution (target distribution) given by

\[f(x) = \frac{1}{\pi} \left[\frac{1}{(1 + {x^2})}\right]; \;\:\:-\infty < x < \infty \]
or
\[f(x) \propto \frac{1}{{(1 + {x^2})}}\]

We can use a normal distribution for the proposal distribution.

The MATLAB code for Metropolis-Hasting algorithm is as follow.

n = 10000;

sig = 2;

x = zeros(1,n);

x(1)=normrnd(0,sig);

cauchy = @(x) 1./(1+x.^2);

norm = @(x,mu,sig) 1/sig * exp(-0.5*((x-mu)/sig).^2);

for i=2:n  

    y = normrnd(x(i-1),sig);

    alpha = min([1, cauchy(y)*norm(x(i-1),y,sig)/...

                   (cauchy(x(i-1))*norm(y,x(i-1),sig))]);

    if rand() < alpha

        x(i)=y;

    else

        x(i)=x(i-1);

    end

end

%discard the first 500 points

cv = x(501:end);

% Plot the distribution and compare to the true cauchy distribution

[n,x]=hist(cv,50);

width=x(2)-x(1);

bar(x,n/(sum(n)*width),1); hold on;

lx = linspace(min(x),max(x),100);

plot(x,cauchy(x)/pi,'r','Linewidth',2); hold off;

axis tight;