趁著周末，學習了此算法。一個重要的作用就是用來模擬目標分布的樣本。下面看看具體情況。

---

1.名詞解釋
MCMC方法就是*構造合適的馬爾科夫鏈進行抽樣而使用蒙特卡洛方法進行積分計算,既然馬爾科夫鏈可以收斂到平穩分布。我們可以建立一個以π為平穩分布的馬爾科夫鏈，對這個鏈運行足夠長時間之后，可以達到平穩狀態。此時馬爾科夫鏈的值就相當于在分布π(x)中抽取樣本。利用馬爾科夫鏈進行隨機模擬的方法就是MCMC。

第一個MC: Monte Carlo(蒙特卡洛)。這個簡單來說是讓我們使用隨機數（隨機抽樣）來解決計算問題。在MCMC中意味著：后驗分布作為一個隨機樣本生成器，我們利用它來生成樣本(simulation)，然后通過這些樣本對一些感興趣的計算問題（特征數，預測）進行估計。

第二個MC：Markov Chain(馬爾科夫鏈)。第二個MC是這個方法的關鍵，因為我們在第一個MC中看到，我們需要利用后驗分布生成隨機樣本，但后驗分布太復雜，當這些樣本獨立時，利用大數定律樣本均值會收斂到期望值。如果得到的樣本是不獨立的，那么就要借助于馬爾科夫鏈進行抽樣，利用Markov Chain的平穩分布這個概念實現對復雜后驗分布的抽樣。

2.馬爾科夫鏈的定義如下：
設θ(t)是一個隨機過程，如果它滿足下面的性質：
p(θ(t+h)=xt+h|θ(s)=xs,s≤t)=p(θ(t+h)=xt+h|θ(t)=xt), 任意的h>0。

這個定義又稱為馬爾科夫性質，對一個馬爾科夫鏈來說，未來狀態只與當前t時刻有關，而與t時刻之前的歷史狀態無關（條件獨立）。

馬爾科夫鏈的一個很重要的性質是平穩分布。簡單的說，主要統計性質不隨時間而變的馬爾科夫鏈就可以認為是平穩的。數學上有馬氏鏈收斂定理，當步長n足夠大時，一個非周期且任意狀態聯通的馬氏鏈可以收斂到一個平穩分布π(x)。這個定理就是所有的MCMC方法的理論基礎。

結論：一個Markov鏈可以由它的初始狀態以及轉移概率矩陣P完全確定。

3.什么是平穩分布？它和求極限概率分布有什么關系呢？

定義：Markov鏈有轉移概率矩陣P，如果有一個概率分布{πi}滿足,則稱為這個Markov鏈的平穩分布。這個定義用矩陣形式寫出來就是π*P=π.

這個定義的含義：如果一個過程的初始狀態X0有平穩分布π,我們可以知道對所有n，Xn有相同的分布π。再根據Markov性質可以得到，對任何k，有Xn,Xn+1,...,Xn+k的聯合分布不依賴于n，顯然這個過程是嚴格平穩的，平穩分布也由此得名！！

4.拒絕抽樣

基本思想是，我們需要對一個分布f(x)進行采樣，但是卻很難直接進行采樣，所以我們想通過另外一個容易采樣的分布g(x)的樣本， 用某種機制去除掉一些樣本，從而使得剩下的樣本就是來自與所求分布f(x)的樣本。

具體的采樣過程如下：

對于g(x)進行采樣得到一個樣本xi, xi ~ g(x);

對于均勻分布采樣 ui ~ U(a,b);

如果ui<= f(x)/[M*g(x)], 那么認為xi是有效的樣本；否則舍棄該樣本； （# 這個步驟充分體現了這個方法的名字：接受-拒絕）

反復重復步驟1~3，直到所需樣本達到要求為止。

示例：產生服從beta(2,7)的隨機數。提議分布g取為均勻分布,常數M取為beta(2,7)的密度函數的最大值。

a <- 2 b <- 7 xmax <- (a-1)/(a+b-2) dmax <- xmax^(a-1)*(1-xmax)^(b-1)*gamma(a+b)/(gamma(a)*gamma(b)) y <- runif(1000) x <- na.omit(ifelse(runif(1000) <= dbeta(y,a,b)/dmax,y,NA))

查看ks檢驗:

z <- x[1:323] ks.test(z,"pbeta",2,7) #接受域圖 plot(dbeta(x,2,7)~x,type = "l")

5.M-H抽樣

可逆馬氏鏈的可逆性經常表示為(細致平衡方程,detailed balance equations) ,從而如果一個目標分布滿足此細致平衡方程，則容易驗證

根據 平穩分布的定義。

MH算法：
下面按如下方式定義一個馬氏鏈：
1.從時刻t的狀態i轉移到下個時刻的狀態，由轉移核生成一個候選的狀態j；
2.以概率min{1,pj/pi}接受下一時刻的狀態為Xt+1=j,否則Xt+1=i
這里用到了馬爾科夫鏈的另一個性質，如果具有轉移矩陣P和分布π(x)的馬氏鏈對所有的狀態i,j滿足下面的等式：π(i)p(i,j)=π(j)p(j,i)
這個等式稱為細致平衡方程。滿足細致平衡方程的分布π(x)是平穩的。 所以我們希望抽樣的馬爾科夫鏈是平穩的，可以把細致平衡方程作為出發點。
示例：
使用MH抽樣法,從Rayleigh分布中抽樣,Rayleigh分布的密度為：
，取自由度為Xt的卡方分布為提議分布，則使用MH算法如下：
1). 令g(.|x)為X2(df=x)
2).從X2(1)中產生X0，并存在X[1]中。
3). 對i=2,….,N，重復,
(a) 產生備選樣本,從X2(df=Xt)=X2(df=X[i-1])中產生Y
(b) 產生U~U(0,1)
(c)在Xt=X[i-1]，計算
若U <= r(Xt,Y ),則接受Y,令Xt+1 =Y,否則令Xt+1=Xt

例如：
利用M-H抽樣方法從Rayleigh分布中抽樣，此分布的密度函數為：

# 計算Rayleigh密度函數在某點的值 f <- function(x,s) {if(x < 0) return (0)stopifnot(s > 0) # if not all true,stop is calledreturn((x/s^2)*exp(-x^2/(2*s^2))) }N <- 10000 s <- 4 x <- numeric(N) x[1] <- rchisq(1,df=1) #初始化提議分布 k <- 0 u <- runif(N) for(i in 2:N) {y <- rchisq(1,df = x[i-1]) #候選點num <- f(y,s)*dchisq(x[i-1],df = y) den <- f(x[i-1],s)*dchisq(y,df = x[i-1])if (u[i] <= num/den) x[i] <- yelse {x[i] <- x[i-1]k <- k+1 #y is rejected} } print(k)# 做樣本路徑圖(trace plot) index <- 500:1000 y1 <- x[index] plot(index,y1,type = "l",main = "",ylab = "x")#在候選點被拒絕的時間點上鏈沒有移動,因此圖中有很多短的水平平移

stopifnot()對函數參數進行檢驗，可看幫助文檔

接下來比較Rayleigh分布的分位數和MH算法下得到樣本分位數

b <- 2001 # discard the burnin(2000個) sample y <- x[b:N] a <- ppoints(100) #產生概率點,用來計算分位點 QR <- s*sqrt(-2*log(1-a)) # quantiles of Rayleigh Q <- quantile(x,a) #quantitles of MH qqplot(QR,Q,main = "",xlab = "Rayleigh Quantiles",ylab = "Sample Quuantiles") hist(y,breaks = "scott",main = "",xlab = "",freq = FALSE) lines(QR,f(QR,4))

6.隨機游動的MH算法

使用提議分布N(Xt,s^2)和隨機游動Metropolis算法產生自由度為V的t分布隨機數

rw.Metropolis <- function(n,sigma,x0,N){# n: degree of freedom of t distribution# sigma: standard variance of proposal distribution N(xt,sigma)# x0: initial value# N: size of random numbers required.x <- numeric(N)x[1] <- x0u <- runif(N)k <- 0for(i in 2:N){y <- rnorm(1, x[i-1], sigma)if(u[i] <= dt(y,n)/dt(x[i-1],n))x[i] <- yelse{x[i] <- x[i-1]k <- k+1}}return (list(x=x,k=k)) }n <- 4 N <- 2000 sigma <- c(.05,.5,2,16) x0 <- 25 rw1 <- rw.Metropolis(n,sigma[1],x0,N) rw2 <- rw.Metropolis(n,sigma[2],x0,N) rw3 <- rw.Metropolis(n,sigma[3],x0,N) rw4 <- rw.Metropolis(n,sigma[4],x0,N)# rate of candidate points rejected print(c(rw1$k,rw2$k,rw3$k,rw4$k)/N)

拒絕率

[1] 0.0040 0.1220 0.4645 0.9085

只有第二個鏈的拒絕率在區間[0.15,0.5]
在不同的提議分布方差下,檢驗所得鏈的收斂性

par(mfrow=c(2,2)) refline <- qt(c(0.25,0.975),df=n) rw <- cbind(rw1$x,rw2$x,rw3$x,rw4$x) for(j in 1:4){plot(rw[,j],type = "l", #trace plotxlab = bquote(sigma==.(round(sigma[j],3))),ylab = "x",ylim = range(rw[,j]))abline(h = refline) } par(mfrow=c(1,1)) # reset to default

路徑圖

可以看出,sigma^2=0.05時,增量太小，幾乎每個候選點都被接受了，鏈在2000次迭代后還沒有收斂。
Sigma^2=0.5,鏈的收斂較慢;sigma^2=2時，鏈很快收斂；而sigma^2=16時，接受的概率太小，使得大部分候選點都被拒絕了(圖形放大看,有很多小區間-)。

bquote(expr) #引用表達式,其中封裝在.()中的要被計算出來

7.Gibbs抽樣

可以看做MH算法當alpha=1的一個特例,用于目標分布為多元分布的情況。

假設在多元分布中所有的一元條件分布都是可以確定的,記m維隨機向量X=(X1,X2,…,Xm)`
X-i表示X中去掉分量Xi后剩余的m-1維向量,那么一元條件分布就是f(xi|x-i)

Gibbs抽樣就是在這m個條件分布中迭代產生樣本，算法：
1)給出初值X(0);
2)對t=1,…,T進行迭代

• (a)令x1=X1(t?1)；
• (b)依次更新每一個分量，即對i=1,…,m,
  (i)從f(Xi∣X?i(t?1))中產生抽樣Xi(t)；
  (ii)更新xi(t)=Xi(t)；
• (c)令X(t)=(X1(T),X2(T),…,Xm(T))′(每個抽取的樣本都被接受了)；
• (d)更新t。
  在這個算法里，對每一個狀態t,X(t)的分量是依次更新的。這個分量更新的過程是在一元分布f(Xi∣X?i)中進行的，所以抽樣是比較容易的。

示例：使用Gibbs抽樣抽取二元正態分布N(μ1,μ2,σ12,σ22,ρ)的隨機數 在二元正態分布的條件下，兩個分量的一元條件分布依然是正態分布：
f(x1∣x2)～N(μ1+ρσ1σ2(x2?μ2),(1?ρ2)σ12)
f(x2∣x1)～N(μ2+ρσ2σ1(x1?μ1),(1?ρ2)σ22)

p_ygivenx <- function(x,m1,m2,s1,s2) {return(rnorm(1,m2+rho*s2/s1*(x-m1),sqrt(1-rho^2)*s2)) } p_xgiveny <- function(y,m1,m2,s1,s2) {return (rnorm(1,m1+rho*s1/s2*(y-m2),sqrt(1-rho^2)*s1)) }N = 5000 K = 20 # iteration in each sampling x_res = vector(length = N) y_res = vector(length = N) m1 = 10; m2 = -5;s1 = 5;s2 = 2 rho = 0.5 y = m2for(i in 1:N) #采樣N次 { for(j in 1:k) # 每次采樣迭代20次{x = p_xgiveny(y,m1,m2,s1,s2)y = p_ygivenx(x,m1,m2,s1,s2)}#print(x)x_res[i] = x;y_res[i] = y; }hist(x_res,freq = F) dev.new()library(MASS) valid_range = seq(from = N/2,to = N,by = 1) MVN.kdensity <- kde2d(x_res[valid_range],y_res[valid_range],h = 10) plot(x_res[valid_range],y_res[valid_range],col = "blue",xlab = "x",ylab = "y") contour(MVN.kdensity,add = TRUE)#二元正態分布等高線圖

8.關于鏈的收斂有這樣一些檢驗方法。

（1）圖形方法 這是簡單直觀的方法。我們可以利用這樣一些圖形：
（a）跡圖（trace plot）：將所產生的樣本對迭代次數作圖，生成馬氏鏈的一條樣本路徑。如果當t足夠大時，路徑表現出穩定性沒有明顯的周期和趨勢，就可以認為是收斂了。
(b)自相關圖（Autocorrelation plot）：如果產生的樣本序列自相關程度很高，用跡圖檢驗的效果會比較差。一般自相關隨迭代步長的增加而減小，如果沒有表現出這種現象，說明鏈的收斂性有問題。
（c）遍歷均值圖（ergodic mean plot）：MCMC的理論基礎是馬爾科夫鏈的遍歷定理。因此可以用累積均值對迭代步驟作圖，觀察遍歷均值是否收斂。
其它還有 （2）蒙特卡洛誤差 （3）Gelman-Rubin方法

參考1
參考2

總結

以上是生活随笔為你收集整理的马尔科夫蒙特卡洛算法(MCMC)的全部內容，希望文章能夠幫你解決所遇到的問題。

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