Hastings-Metropolis演算法(離散型)

01-25

1. 問題

給定 $m$ 個正數 $b(j),j=1,2,cdots,m$ ,( $m$ 可以取無窮大）令 $B=sum^m_{j=1}b(j)$ 。一般來說 $m$ 很大，這樣 $B$ 就很難計算。如何模擬一個隨機變數其分布率為 $pi(j)=b(j)/B,quad j=1,2,cdots,m$ ？

2.Hastings-Metropolis演算法

下面將我們構造一個可逆的馬爾科夫鏈，其極限分布為 $pi(j)=b(j)/B,quad j=1,2,cdots,m$ 。假設 $Q=(q_{ij})_{m imes m}$ 為狀態空間 ${1,2,cdots,m}$ 上的不可約隨機矩陣，即(1) $q_{ij}geq0$ ；(2) $sum^m_{j=1}q_{ij}=1,j=1,2,cdots,m$ 。

生成一個可逆的馬爾科夫鏈 $X={X_n:n=0,1,2,cdots}$ 方法為：

(1) 假設 $X_n=i,$ 那麼根據分布 $q(i,.)={q_{ij},j=1,2,cdots,m}$ 生成一個 $j$ ,

(2) 然後再以 $alpha_{i,j}$ 的概率令 $X_{n+1}=j$ ,以 $1-alpha_{i,j}$ 的概率令 $X_{n+1}=i$ ，其中 $alpha_{i,j}=minleft{frac{pi(j)q(j,i)}{pi(i)q(i,j)},1 ight}=minleft{frac{b(j)q(j,i)}{b(i)q(i,j)},1 ight}$ 。

該方法稱為Hastings-Metropolis演算法。

3. 命題

$X={X_n:n=0,1,2,cdots}$ 是可逆的馬爾科夫鏈，其轉移矩陣為 $P=(p_{i,j})$ ,其中 $p_{ij}=q_{i,j}alpha_{i,j}$ , $pi$ 為 $X$ 的平穩分布。

4.Hastings-Metropolis演算法步驟

(1). 選擇一個狀態空間為 ${1,2,cdots,m}$ 的不可約Markov轉移概率矩陣 $Q=(q_{ij})_{m imes m}$ ，從 ${1,2,cdots,m}$ 中任選一個數 $k$ ；

(2). 令 $n=0,X_0=k$ ;

(3). 根據分布 $q(X_n,.)$ 生成一個隨機數 $X$ ，使得 $P{X=j}=q(X_n,j)$ ,並且生成一個隨機數 $U$ ；

(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). 返回到(3).

5.例子

假設 ${x_1,x_2,cdots,x_n}$ 為 ${1,2,cdots,n}$ 的一個排列， $l$ 為滿足 $sum^n_{j=1}jx_j>a$ 的排列的個數，其中 $a>0$ 為給定的常數，如何估計 $l$ 的個數？下面我們利用Hastings-Metropolis演算法進行估計。

(1) 如何在排列集合上定義一個轉移矩陣 $Q=(q(s,t)),s,t$ 為任意兩個排列?

首先我們先定義鄰域的概念。假設 $s$ 是一個排列，若 $t$ 為 $s$ 上任意兩個位置的數交換後得到的排列，那麼稱 $t$ 為 $s$ 的鄰居。例如 $(1,2,4,3)$ 就是 $(1,2,3,4)$ 的一個鄰居。令 $N(s)={t:t$ 為 $s$ 的鄰居 $}$ , 這樣我們定義 $q(s,t)=frac{1}{|N(s)|}$ , $tin N(s)$ ,即從 $s$ 出發，等可能地達到它任意一個鄰居。

(2) 定義 $alpha$

因為極限分布為 $pi={pi_j=1/l:j=1,2,cdots,l}$ ,這樣 $forall s,t,pi(s)=pi(t)$ .於是

$alpha(s,t)=min(|N(s)/N(t)|,1)$ .如果 $N(s)geq N(t)$ ,則 $alpha(s,t)=1$ .

(3) 演算法步驟

步驟1 任意選擇一個排列 $s$ ；
步驟2 令 $n=0,X_n=s$ ;
步驟3 從 $s$ 的鄰域中隨機選取一個排列 $t$ ;
步驟4 如果 $N(s)geq N(t)$ ，則 $X_{n+1}=s$ ;如果 $N(s)< N(t)$ ，則生成一個隨機數 $U$ ,若 $U<|N(s)|/N(t)$ ,則 $X_{n+1}=t$ ,否則 $X_{n+1}=s$ .
步驟5 $n=n+1$
步驟6 返回到步驟3

(4) 模擬的馬爾科夫鏈 $X$ 的極限分布為 $pi={pi_j=1/l:j=1,2,cdots,l}$ 。

6.作業

（1）給出上面例子的R程序；

（2）假設 $pi(0)=frac{1}{3},pi(x)=frac{2^{-x+1}}{3},x=1,2,cdots$ ,請利用Hastings-Metropolis演算法模擬分布 $pi$ 的方差。

解： $pi$ , $xin Z$ ,首先利用Hastings-Metropolis演算法模擬分布 $pi$ .令馬爾科夫鏈為 $X_1,X_2,cdots,X_n$ 其極限分布為 $pi$ ，那麼 $|X_1|,|X_2|,cdots,|X_n|$ 的極限分布為 $pi$ .