Chapter 5 MCMC techniques

5.1 Markov chains

A discrete time Markov chain \(\{X_k\}_k\) is a sequence of rvs such that the distribution of \(X_k\) depens only on \(X_{k-1}\) and not on the previous rvs.

The Markov chain is homogeneous if the conditional distribution does not depend on time point. In such a case, the probabilities \(P_{i,j}=P(X_k=j|X_{k-1}=i)\) are called transition probabilities.

A pmf \(\pi\) is stationary distribution of the Markov chain if for every state \(\pi(j)=\sum P_{i,j}\pi(i)\).

Ergodic Markov chains have a unique stationary distribution and it will be reached as a limit distribution from any initial distribution.

MCMC methods consist in generating ergodic Markov chains whose limit distribution is a goal distribution. They do NOT produce sequences of independent observations!

5.2 Metropolis-Hastings

Our goal is to simulate from dmf \(f\) and we use an instrumental conditional dmf \(g(\cdot|\cdot)\) from which it is simple to simulate. The support of \(g\) must contain that of \(f\).

Set \(X_0\) such that \(f(X_0)>0\) and \(k=0\).
Set \(k=k+1\). Generate \(Y_k\sim g(\cdot|X_{k-1})\) and \(U\) random number in \((0,1)\) independent.
If \(U\leq\alpha(X_{k-1},Y_k)\), set \(X_k=Y_k\), otherwise \(X_k=X_{k-1}\).
If \(k<n\) go to Step \(2.\), otherwise return \(X_0,\ldots,X_n\).

\[\alpha(x,y)=\min\left\{\frac{f(y)g(x|y)}{f(x)g(y|x)},1\right\}\]

Example (Poisson)

In order to simulate a Poisson rv, we take as instrumental probability mass function (given \(x\neq 0\)), the one that assesses probability \(1/2\) to \(x-1\) and probability \(1/2\) to \(x+1\). If \(x=0\), then the probability of \(0\) is \(1/2\) and the probability of \(1\) is \(1/2\)

rinst <- function(x) {
  if(x==0) return(sample(1,x=c(0,1)))
  else return(sample(1,x=c(x-1,x+1)))
}

poisson.metro=function(lamb,ini,n){ 
 x <- c(ini)
 y <- rinst(ini)
 for(k in 2:n) {
  u <- runif(1)
  alpha<-lamb^y/factorial(y)/(lamb^x[k-1]/factorial(x[k-1]))
  if(u<=alpha) x <- c(x,y)
  else x <- c(x,x[k-1])
  y <- rinst(x[k])
 } 
return(x)
}

lmbd <- 1
plot(cumsum(poisson.metro(lamb=lmbd,ini=2,n=1000))/(1:1000),
     type="l")
abline(h=lmbd,col="red")

5.3 Gibbs sampling

Consider a random vector \({\bf X}=(X^{(1)},\ldots,X^{(d)})\) with joint dmf \(f\) difficult to simulate from, but whose conditional densities of each component given the remaining componets \(f_i(\cdot|x_k^{(1)},\ldots,x_k^{(i-1)},x_k^{(i+1)},\ldots,x_k^{(d)})\) are easy to simulate from.

The Gibbs sampler is Metropolis-Hastings algorithm with \(g(\cdot|{\bf X})\) having marginal densities as above. In such a case, \(f({\bf y})g({\bf x}|{\bf y})=f({\bf x})f({\bf y})\), so all candidates will be accepted.

Set \({\bf X}_0\) such that \(f({\bf X}_0)>0\) and \(k=0\).
Set \(k=k+1\). For \(i=1,\ldots,d\), generate \(X_k^{(i)}\) from \(f_i(\cdot|X_k^{(1)},\ldots,X_k^{(i-1)},X_{k-1}^{(i+1)},\ldots,X_{k-1}^{(d)})\).
If \(k<n\) go to Step \(2.\), otherwise return \({\bf X}_0,\ldots,{\bf X}_n\).

Example

If the random vector \((X,Y)^t\) follows a bivariate normal distribution with mean vector \((0,0)^t\) and covariance matrix having 1s on the maing diagonal and \(\rho\)s elsewhere (\(X\) and \(Y\) are standard normal random variables with correlation \(\rho\))

\(X\sim{\rm N}(\rho Y,\sqrt{1-rho^2})\) and \(Y\sim{\rm N}(\rho X,\sqrt{1-\rho^2})\)

set.seed(1) ; n <- 10^4 ; x <- c(0) ; y <- c(0)
rho <- 0.8 ; sigma <- sqrt(1 - rho^2)
for (i in 2:n) {                                
  x <- c(x,rnorm(1,mean=rho*y[i-1],sd=sigma))
  y <- c(y,rnorm(1,mean=rho*x[i],sd=sigma))
}
c(mean(x),mean(y))

## [1] -0.01986475 -0.01208011

c(sd(x),sd(y))

## [1] 0.9947670 0.9899634

cor(x,y)

## [1] 0.7984191

par(mfrow=c(2,2))
plot(x[1:100],type="l") ; acf(x)
s.norm <- rnorm(100) ; plot(s.norm,type="l") ; acf(s.norm)