[stat.simulation] Hasting-Metropolis Algorithm

问题背景：我们有一些观测数据X，这些数据假设是取值为1,...,m；我们还知道每个数据观测到的频数为: $b(j), j=1,...,m. B=\sum_{j=1}^mb(j)$

　　但是我们现在无法计算B的大小。（这是一个假设，毕竟计算一串数字的和不是难事）

问题：我们需要通过仿真产生一串随机变量，并且它们的概率分布函数为：

$\pi(j)=b(j)/B, j=1,...,m$

分析：如果B是可以计算的，那么 $\pi$ (j)自然也是可以计算的。然后自然很容易随机生成服从这个概率分布的一串随机数。但是B不能计算。。。我们可以采用一个曲线救国的方案。

　　直观上来考虑这个问题，产生随机数时是一个一个地产生随机数，每个随机数取值为1,...,m中某一个。我们可以：

1）把一个随机数看成一个状态；

2）一个随机数的产生取决于前一个随机数，那么每个状态取决于前一个状态；

问题等效：对于状态{1,...,m}，有某个Markov Chain状态转移矩阵P， $\text{P}(i,j)=\Pr(i\rightarrow j)$ ,随机选择一个初始状态 $X[1]\in{1,...,m}$ ，根据状态转移矩阵P,，依次产生一个状态序列 $\{X[1], X[2], X[3], ...\}$ 。但是，重点是我们需要设计这个转移矩阵P，使得最终的产生的状态序列的概率分布为 $\pi$ ，也就是说这个Markov Chain的最终稳定状态时各个状态的概率分布为 $\pi$ 。

Hasting-Metropolis算法描述

　　算法引入了一个随机选择的Markov Chain转移矩阵Q， $\text{Q}(i,j)=\Pr(i\rightarrow j)$ 。注意Q不是上面讲的P，而是用来构造P的一个辅助状态转移矩阵。构造的过程：

1）根据Q，从当前状态，设为i，转移到状态下一个状态，设为j；

2）引入另一个概率函数 $\alpha(i, j)$ ，这个概率决定转移到下一个状态j或者停留在现在的状态i；

结合上述1,2）可以得到一个新的转移矩阵，并且可以经过精心构造 $\alpha(i, j)$ ，使得这个矩阵作为解决“等效问题”中的P矩阵，下面来看如何“精心构造”：

根据1,2）

$\text{P}(i,j)=q(i,j) \alpha(i,j), \text{when } i\ne j;$

$\text{P}(i,i)=q(i,i)+q(i,j) (1-\alpha(i,j))$

根据平稳Markov Chain的“时序可逆”性质：

$\pi(i)P(i,j)=\pi(j)P(j,i) \text{ for }i\ne j$

$\Leftarrow \pi(i)q(i,j)\alpha(i,j)=\pi(j)q(j,i)\alpha(j,i)\\ \Leftarrow \left\{ \begin{array}{c} \alpha(i,j)=1, \\ \alpha(j,i)=\frac{\pi(i)q(i,j)}{\pi(j)q(j,i)} \\ \end{array} \right$

$\left\{ \begin{array}{c} \alpha(j,i)=1, \\ \alpha(i,j)=\frac{\pi(j)q(j,i)}{\pi(i)q(i,j)} \\ \end{array} \right$

因为 $\alpha(i, j)$ 是一个概率值，必须小于1，所以最终 $\alpha(i, j)$ 的构造为：

$\alpha(i,j)=\min\{1,\frac{\pi(j)q(j,i)}{\pi(i)q(i,j)}\}$

但是 $\pi$ 是未知的（因为B无法求值），好在上式中是 $\pi$ (i)/ $\pi$ (j)=b(i)/b(j)，那么：

$\alpha(i,j)=\min\{1,\frac{b(j)q(j,i)}{b(i)q(i,j)}\}$

Hasting-Metropolis算法流程

1）选择一个不可约Markov Chain概率转移矩阵 $\text{Q}(i,j)=\Pr(i\rightarrow j)$ ；随机选择初始状态 $k\in \{1,...,m\}$ ;

2) let n=1, X[n]=k;

3) 生成随机数 $X \text{, such that } \text{P}\{X=j\}=q(X[n],j)$ ，生成随机数 U∈（0, 1）

4) 如果 $U<$b(X)q(X,X[n])$/$b(X[n]q(X[n],X)$$ ，则选择NS=X；否则选择NS=X[n];

5) n=n+1, X[n]=NS;

6) go to 3)

说明，以上在讨论函数 $\alpha(i, j)$ 和算法流程中，都没有特别考虑i = j的情况。原因是i=j时“不失一般性”。此时 $\alpha(i, j)$ =1，算法步骤4）一定会选择NS=X,而此时X==X[n]。所以，算法步骤3）中如果产生的X==X[n]，那么4）定然会保持原来状态从而X[n+1]<-X[n]。

仿真例子

通过R语言实现用H-M算法采样服从混合高斯分布密度函数。

主程序：

source("intgMatrix.R")
source("sampleByPr.R")

FF = 10000# sample points

br<-rnorm(FF/2, 0, 10)
br2<-rnorm( FF/2, 60, 10)
br <- c(br, br2)
br<-round(br)

tbl <- table(br)

b<-unname(tbl)   # a sequence of number

m <- length(b)
#generate Q
r<-c()
for(i in 1:m)
{
  r_ <- runif(m, 0, 1)
  r_ <- r_ / sum(r_)  # nomalize
  r<- c(r, r_)
}
Q<- t(matrix(r, m, m))
Qint <- intgMatrix(Q)     # row-based integral

# start to sample sn data
sn = FF  # #sample
X <- c()  # sample data (index of b)

# initial step
k<- round(runif(1, 1, m)) 
X[1] <- k

# interation steps
for(n in 1:(sn-1))
{
  PrX_j <- Qint[X[n], ]  # last sample's index as row-index in Q
  while(TRUE){
    X_tmp <- X[n]
    X_tmp <- sampleByPr(PrX_j) # the index of X[n]
    if(X_tmp != X[n]){
      break
    }
  }
  
  U <- runif(1, 0, 1) 
  # for alph(i, j)
  al_ <- (b[X_tmp] * Q[X_tmp, X[n]]) / (b[X[n]] * Q[X[n], X_tmp])
  if(al_ > 1){al_ <- 1}
  NS <- X[n]
  if(U < al_) {
    NS <- X_tmp
  }
  
  X[n + 1] <- NS
  
}

# plot


bnn<-as.numeric(b)
plot(1:length(bnn),bnn,col="red")
Xnn <- as.numeric(table(X))
lines(1:length(Xnn),Xnn/rt,col="green")

概率积分矩阵，Qi[i,k]=∑_j=1,..,kQ[i,j]

intgMatrix <- function(Q)
{
  Qi<-matrix(nrow=dim(Q)[1], ncol = dim(Q)[2])
  Qi[, 1] = Q[,1]
  for(i in 1:dim(Q)[1]){
    for(j in 2:dim(Q)[2]){
      Qi[i,j] <- Qi[i, j- 1] + Q[i, j]
    }
  }
  return(Qi)
}

根据概率积分向量生成采样点。

概率积分向量pi[i]=∑_j=1,..,ip[i],其中p[i]=Pr{X=i}为原始概率分布

sampleByPr <- function(Pr_vector)
{
  ru <- runif(1,0,1)
  
  for(i in 1:length(Pr_vector)){
    if(Pr_vector[i] > ru){
      return(i)
    }
  }
}