MH算法

MH算法的两个例子，第一个目标分布为双指数分布的MH算法实例，proposed分布是sd=1的正态分布。另一个是使用硬币去估计骰子的分布，目标分布也即每类的概率都是$\frac{1}{6}$，给出了状态转移矩阵作为propose。

1

MH-Algorithm for from the (standard) double exponential distribution.

\[\pi(x)=\frac{1}{2}e^{-|x|}\]

Input: target $\pi(x)$ and iteration T

Initialize: $t=0$ and $x=x_0=0$

Repeat until t=T:

My proposal distribution $g(x'|x_t)$ is $N(\mu,1)$, and $x'$ is generated from $N(\mu,1)$.
$A(x'|x_t)=min(1,\frac{\pi(x')g(x'|x_t)}{\pi(x_t)g(x_t|x')})$
Generate $U\sim unif([0,1])$, if $U<=A(x'|x_t)$, accept $x_{t+1}=x'$, if $U>A(x'|x_t)$, accept $x_{t+1}=x_t$. Above is the acceptance probability
t->t+1

pi_target=function(x){
  result=0.5*exp(-abs(x))
  return(result)
}


MH_algorithm=function(x){
  out_mh=rep(0,x)
  out_mh[1]=0
  for (i in 2:x){
    current=out_mh[i-1]
    proposed=rnorm(1,mean=current,sd=1)
    a=(pi_target(proposed)*rnorm(1,mean=current,sd=1))/(pi_target(current)*rnorm(1,mean=proposed,sd=1))
    a=min(1,a)
    if (a>runif(1)){
      out_mh[i]=proposed
    } 
    else{
      out_mh[i]=current
    }
  }
  return(out_mh)
}

I plot the hist of $x_t$ and the real curve.

a=MH_algorithm(1000000)
b=hist(a,probability = TRUE)
xx = seq(-30,30,length=1000000)
lines(xx,pi_target(xx),col="red")

Then we check the convergence of the sampling

error=rep(0,length(seq(from=1000,to=100000,by=1000)))
for (i in seq(from=1000,to=100000,by=1000)){
  a=MH_algorithm(i)
  a=hist(a, plot = FALSE)
  a_i=a$density
  xx = a$mids
  b=sort(pi_target(xx))
  error[i/1000]=sum((a_i-b)^2)/i
}

plot(x=seq(1000,100000,by=1000),y=error,main = "Plot of Convergence", xlab = "Iteration Number",xlim = c(1000,100000))

As we can see the covergence above as the iteration bigger, the error tend to be 0.

2

MH-Algorithm for 6-sided dice.

If we neither have a dice to roll or a computer to simulate the rolls. We only can flip fair coins. The rule we flip the coin is: go up if it is a head and go down if it is a tail. For example, if X = 1, we always propose X = 2 for the next state of the dice, while if X = 2, we have 50% chance to propose X = 1 and $50%$ chance to propose X = 3, and if X = 6, we always propose X = 5. This random experiment is essentially a Metropolis-Hastings algorithm to modify the simple random walk as flipping a fair coin to have the right probabilities for a 6-sided dice.

We can denote the transition matrix as $$P=

\[\begin{pmatrix} 0& 1& 0& 0& 0& 0\\ 0.5& 0& 0.5& 0& 0& 0\\ 0& 0.5& 0& 0.5& 0& 0\\ 0& 0& 0.5& 0& 0.5& 0\\ 0& 0& 0& 0.5& 0& 0.5\\ 0& 0& 0& 0& 1&0 \end{pmatrix}\]

Input: target $\pi(x)$ and iteration T

Initialize: $t=0$ and $x_0=1$

Repeat until t=T:

My proposal is transition matrix P, and $x'$ is generated from it.
$A(x'|x_t)=min(1,\frac{\pi(x')g(x'|x_t)}{\pi(x_t)g(x_t|x')})$
Generate $U\sim unif([0,1])$, if $U<=A(x'|x_t)$, accept $x_{t+1}=x'$, if $U>A(x'|x_t)$, accept $x_{t+1}=x_t$. Above is the acceptance probability
t->t+1

Trans_matrix=function(x){
  if (x==1){
    return(2)
  }
  if (x==6){
    return(5)
  }
  a=runif(1)
  if (x==2){
    if (a>0.5){
      return(1)
    } else{
      return(3)
    }
  }
  if (x==3){
    if (a>0.5){
      return(2)
    } else{
      return(4)
    }
  }
  if (x==4){
    if (a>0.5){
      return(3)
    } else{
      return(5)
    }
  }
  if (x==5){
    if (a>0.5){
      return(4)
    } else{
      return(6)
    }
  }
}

G_func=function(x,y){
  if (abs(x-y)>1){
    return(0)
  }else{
    
    if (x==1 && y==2){
      return(1)
    } else if (x==6 && y==5){
      return(1)
    }else{
      return(0.5)
    }
  }
  
}


MH_dice=function(x){
  out_mh=rep(0,x)
  out_mh[1]=1
  for (i in 2:x){
    current=out_mh[i-1]
    proposed=Trans_matrix(current)
    a=G_func(proposed,current)/G_func(current,proposed)
    a=min(1,a)
    if (a>runif(1)){
      out_mh[i]=proposed
    } 
    else{
      out_mh[i]=current
    }
  }
  return(out_mh)
}

Above is the MH algorithm code for a dice only flip fair coins

And then we put iteration as 1000000, summary the output, we can find that all of the result values are approximatly the same.

a=MH_dice(1000000)
hist(a,probability = TRUE)

table(a)

## a
##      1      2      3      4      5      6 
## 166373 166513 166717 166983 166946 166468

What’s more, I check the convergence of the algorithm.

error=rep(0,length(seq(from=1000,to=100000,by=1000)))
for (i in seq(from=1000,to=100000,by=1000)){
  a=MH_dice(i)
  b=table(a)/i
  error[i/1000]=sum((b-1/6)^2)
}

plot(x=seq(1000,100000,by=1000),y=error,main = "Plot of Convergence", xlab = "Iteration Number")

MH_Algorithm

MH算法

1

2