Metropolis-Hasting算法

Step 1: 選擇一個不可約Markov鏈轉移概率q(i,j),i,j$\in$S.再從S={1，2，… ,m}中選擇某個整數.
Step 2: 令n=0以及$x_0=k$.
Step 3: 產生一個隨機變量X使得$P\{X=j\}=q(X_n,j)$, 再產生一個隨機數U.
Step 4: 如果$U<\frac{\pi (X)q(X,X_n)}{\pi (X_n)q(X_N,X)}$, 則MH = X； 否則MH = $X_n$.
Step 5：$n\leftarrow n+1,X_n=MH$.
Step 6: 返回Step 3.

建議分布的選擇方法：

（1）Metropolis選擇；（2）獨立抽樣；（3）單元素MH算法

舉例說明：

設觀測變量$Y=(Y_1,Y_2,...,Y_n)$具有如下分布：
$$\pi (y_i|k,\theta ,\lambda )=\frac{{\theta }^{y_i}{e}^{?\theta }}{y_i!},i=1,2,...,k,$$
$$\pi (y_i|k,\theta ,\lambda )=\frac{{\lambda }^{y_i}{e}^{?\lambda }}{y_i!},i=k+1,k+2,...n.$$
假定各參數的先驗分布為：
$$\begin{gathered} \pi (\theta |b_1)=Gamma(0.5,b_1), \\ \pi (\lambda |b_2)=Gamma(0.5,b_2), \\ \pi (b_1)=IG(0,1), \\ \pi (b_2)=IG(0,1), \\ \pi (k)=Uniform(1,2,...,n), \end{gathered}$$
其中，$k,\theta ,\lambda$條件獨立，$b_1,b_2$獨立，以及\
$$\begin{gathered} Gamma(\alpha ,\beta )=\frac{1}{\Gamma (\alpha ){\beta }^{\alpha }}{x}^{\alpha ?1}{e}^{?\frac{x}{\beta }}, \\ IG(\alpha ,\beta )=\frac{e^{?\frac{1}{\beta x}}}{\Gamma (\alpha ){\beta }^{\alpha }{x}^{\alpha +1}}, \\ Uniform(1,2,...,n):P\{X=j\}=1/n,j=1,2,...,n \end{gathered}$$
其中$\alpha$為形狀參數，$\beta$為刻度參數。
解答：
首先求解后驗分布：
$$\begin{gathered} \pi (k,\theta ,\lambda ,b_1,b_2|Y)={\Pi }_{i=1}^k\pi (Y_i|k,\theta ,\lambda ){\Pi }_{i=k+1}^n\pi (Y_i|k,\theta ,\lambda ) \\ \times \pi (\theta |b_1)\pi (\lambda |b_2)\pi (b_1)\pi (b_2)\pi (k) \\ ={\Pi }_{i=1}^k\frac{{\theta }^{Y_i}{e}^{?\theta }}{Y_i!}{\Pi }_{i=k+1}^n\frac{{\lambda }^{Y_i}{e}^{?\lambda }}{Y_i!}\times \frac{1}{\Gamma (0.5)b_1}{\theta }^{?0.5}{e}^{?\frac{\theta }{b_1}} \\ \times \frac{1}{\Gamma (0.5)b_2^{0.5}}{\lambda }^{?0.5}{e}^{?\frac{\lambda }{b_2}}\times \frac{e^{?1/b_1}{e}^{?1/b_2}}{b_1{b}_2}\frac{1}{n}. \end{gathered}$$
各參數的滿條件分布為
$$\begin{gathered} \pi (\theta |k,\lambda ,b_1,b_2,Y)\propto Gamma(\sum _{i=1}^k{Y}_i+0.5,\frac{b_1}{k{b}_1+1}), \\ \pi (\lambda |k,\theta ,b_1,b_2,Y)\propto Gamma(\sum _{i=k+1}^n{Y}_i+0.5,\frac{b_2}{(n?k)b_2+1}), \\ \pi (k|\theta ,\lambda ,b_1,b_2,Y)\propto {\theta }^{{\sum }_{i=1}^k{Y}_i}{\lambda }^{{\sum }_{i=k+1}^n{Y}_i}{e}^{?k\theta ?(n?k)\lambda }, \\ \pi (b_1|k,\theta ,\lambda ,b_2,Y)\propto IG(0.5,\frac{1}{1+\theta }), \\ \pi (b_1|k,\theta ,\lambda ,b_1,Y)\propto IG(0.5,\frac{1}{1+\lambda }). \end{gathered}$$
<注:>由于編輯公式麻煩，只給出推導最后結果，中間步驟省略.
###Y的觀測數據

---

3 5 9 3 4 5 5 5 5 13 18 27 8 4 10 8 3 12 10 10 3 9 8 5 9 4 6 1 5 14 7 9 10 8 13 8 11 11 10 11 13 10 3 8 5

---

由于k滿條件分布不是標準分布，因此，用MH算法抽取$k^{i+1}$,而其他分部則使用Gibbs抽樣，其k的建議分布設為$q(k,k^{\prime })=\frac{1}{m?1}$.

使用R語言實現抽樣

mhsampler<-function(NUMIT=10000,dat=Y){n<-length(dat)mchain<-matrix(NA,nr=5,nc=NUMIT)kinit<-floor(n/2)mchain[,1]<-c(1,1,kinit,1,1)for(i in 2:NUMIT){currtheta<-mchain[1,i-1]currlambda<-mchain[2,i-1]currk<-mchain[3,i-1]currb1<-mchain[4,i-1]currb2<-mchain[5,i-1]##sample from full conditional distribution of theta(Gibbs update)currtheta<-rgamma(1,shape=sum(Y[1:currk])+.5,scale=currb1/(currk*currb1+1))##sample from full conditional distribution of lambda(Gibbs update)currlambda<-rgamma(1,shape=sum(Y[(currk+1):n])+.5,scale=currb2/((n-currk)*currb2+1))##sample from full conditional distribution of k(MH update)propk<-sample(x=seq(2,n-1),size=1) #draw one sample at random from uniform(2,..(n-1))##Matropolis accept-reject step(in log scale)logMHratio=sum(Y[1:propk])*log(currtheta)+sum(Y[(propk+1):n])*log(currlambda)-propk*currtheta-(n-propk)*currlambda-(sum(Y[1:currk])*log(currtheta)+sum(Y[(currk+1):n])*log(currlambda)-currk*currtheta-(n-currk)*currlambda)logalpha<-min(0,logMHratio) #alpha=min(1,MHratio)if (log(runif(1))<logalpha){currk<-propk}currk=currk #if we do not sample k (k fixed)## sample from full conditional distribution of b1 (Gibbs update: draw from inverse Gamma)currb1<-1/rgamma(1,shape=.5,scale=1/(currtheta+1))## sample from full conditional distribution of b2 (Gibbs update: draw from inverse Gamma)currb2<-1/rgamma(1,shape=.5,scale=1/(currlambda+1))## update chain with new valuesmchain[,i]=c(currtheta,currlambda,currk,currb1,currb2)}return(mchain) }

最后可以運行apply(mhsampler(),1,mean)獲得$\theta ,\lambda ,k,b_1,b_2$的后驗均值估計。

總結

以上是生活随笔為你收集整理的R语言 MCMC算法及其实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。