R語言中實現馬爾可夫鏈蒙特卡羅MCMC模型

原文連接：http://tecdat.cn/?p=2687

什麼是MCMC，何時使用它？

MCMC只是一個從分佈抽樣的算法。

這只是衆多算法之一。這個術語表明「馬爾可夫鏈蒙特卡洛」，由於它是一種使用「馬爾可夫鏈」（咱們將在後面討論）的「蒙特卡羅」（即隨機）方法。MCMC只是蒙特卡洛方法的一種，儘管能夠將許多其餘經常使用方法看做是MCMC的簡單特例。

正如上面的段落所示，這個話題有一個引導問題，咱們會慢慢解決。

我爲何要從分配中抽樣？

你可能沒有意識到你想（實際上，你可能並不想）。可是，從分佈中抽取樣本是解決一些問題的最簡單的方法。

可能MCMC最經常使用的方法是從貝葉斯推理中的某個模型的後驗機率分佈中抽取樣本。經過這些樣本，你能夠問一些問題：「參數的平均值和可信度是多少？」。

若是這些樣本是來自分佈的獨立樣本，則 估計均值將會收斂在真實均值上。

假設咱們的目標分佈是一個具備均值m和標準差的正態分佈s。顯然，這種分佈的意思是m，但咱們試圖經過從分佈中抽取樣原本展現。

做爲一個例子，考慮用均值m和標準誤差s來估計正態分佈的均值（在這裏，我將使用對應於標準正態分佈的參數）：

咱們能夠很容易地使用這個rnorm 函數從這個分佈中抽樣

 seasamples<-rn 000,m,s)

樣本的平均值很是接近真實平均值（零）：

mean(sa es) ## [1] -0. 537

事實上，在這種狀況下，$ n $樣本估計的預期方差是$ 1 / n $，因此咱們預計大部分值在$ \\ pm 2 \\，/ \\ sqrt {n} = 0.02 $ 10000分的真實意思。

summary(re 0,mean(rnorm(10000,m,s)))) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -0.03250 -0.00580 0.00046 0.00042 0.00673 0.03550

這個函數計算累積平均值（即元素$ k $，元素$ 1,2，\\ ldots，k $除以$ k $）之和。

cummean<-fun msum(x)/seq_along(x) plot(cummaaSample",ylab="Cumulative mean",panel.aabline(h=0,col="red"),las=1)

將x軸轉換爲對數座標並顯示另外30個隨機方法：

能夠從您的一系列採樣點中抽取樣本分位數。

這是分析計算的點，其機率密度的2.5％低於：

 p<-0.025a.true<-qnorm(p,m,s)a.true1## [1] -1.96

咱們能夠經過在這種狀況下的直接整合來估計這個（使用上面的論點）

aion(x)dnorm(x,m,s)g<-function(a)integrate(f,-Inf,a)$valuea.int<-uniroot(function(x)g(a10,0))$roota.int1## [1] -1.96

並用Monte Carlo積分估計點：

1 2a.mc<-unnasamples,p))a.mc1## [1] -2.0231a.true-a.mc1## [1] 0.06329

可是，在樣本量趨於無窮大的極限內，這將會收斂。此外，有可能就錯誤的性質做出陳述; 若是咱們重複採樣過程100次，那麼咱們獲得一系列與均值附近的偏差相同幅度的偏差的估計：

a.mc<-replicate(anorm(10000,m,s),p))summary(a.true-a.mc) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -0.05840 -0.01640 -0.00572 -0.00024 0.01400 0.07880

這種事情真的很常見。在大多數貝葉斯推理中，後驗分佈是一些（可能很大的）參數向量的函數，您想對這些參數的子集進行推理。在一個等級模型中，你可能會有大量的隨機效應項被擬合，可是你最想對一個參數作出推論。在貝葉斯框架中，您能夠計算您感興趣的參數在全部其餘參數上的邊際分佈（這是咱們上面要作的）。

爲了說明這個問題，

考慮在邊長爲$ 2r $的方格內半徑爲$ r $的圓; 空間的「有趣」區域是一個隨機選擇的點位於圓圈內的一個很好的機會。

對於半徑爲$ 2r $ 的立方體中半徑爲$ r $ 的球體，球體的體積爲$ 4 /（3 \\ pi r ^ 3）$，立方體的體積爲$（2d）^ 3 $  ，做爲問題的維數，$ d $增長（使用超立方體中的超球面）

 d<-2:10plot(d,pi^(d/2)/(d*2^(d-1)*gamma(d/2)),logahere")

因此咱們不須要增長不少維度來主要對潛在空間的一小部分感興趣。

a<-function(d,n)mean(repa00)) ## dimension p.interesting## 1 1 0.5219## 2 2 0.2173## 3 3 0.0739## 4 4 0.0218## 5 5 0.0070## 6 6 0.0025## 7 7 0.0006## 8 8 0.0000## 9 9 0.0000## 10 10 0.0000

即便只看4-5個維度，若是咱們試圖對參數空間進行完全整合

爲何「正常統計」不使用蒙特卡洛方法？

對於傳統教學統計中的許多問題，而不是從分佈中抽樣，可使函數最大化或最大化。因此咱們須要一些函數來描述可能性並使其最大化（最大似然推理），或者一些計算平方和並使其最小化的函數。

然而，蒙特卡羅方法在貝葉斯統計中的做用與頻率統計中的優化程序相同，這只是執行推理的算法。因此，一旦你基本知道MCMC正在作什麼，你能夠像大多數人把他們的優化程序看成黑匣子同樣對待它，像一個黑匣子。

馬爾可夫鏈蒙特卡羅

假設咱們想要抽取一些目標分佈，可是咱們不能像從前那樣抽取獨立樣本。有一個使用馬爾科夫鏈蒙特卡洛（MCMC）來作這個的解決方案。首先，咱們必須定義一些事情，以便下一句話是有道理的：咱們要作的是試圖構造一個馬爾科夫鏈，它的難以抽樣的目標分佈做爲它的平穩分佈。

定義

讓$ X\_t $表示在時間$ t $時的一些隨機變量的值。馬爾可夫鏈從某個點$ X\_0 $開始，生成一系列樣本$ \[X0，X1，X2，\\ ldots，Xt\] $，而後遵循一系列隨機步驟。

 假設咱們有一個三態馬爾科夫過程。讓咱們P爲鏈中的轉移機率矩陣：

P<-rbind(a(.2,.1,.7),c(.25,.25,.5))P ## [,1] [,2] [,3]## [1,] 0.50 0.25 0.25## [2,] 0.20 0.10 0.70## [3,] 0.25 0.25 0.50 rowSums(P) ## [1] 1 1 1

條目P[i,j]給出了從狀態i到狀態的機率j（因此這是上面提到的$ P（i \\到j）$。

請注意，與行不一樣，列不必定總和爲1：

colSums(P) ## [1] 0.95 0.60 1.45

這個函數採用一個狀態向量x（其中x[i]是處於狀態的機率i），並經過將其與轉移矩陣相乘來迭代它P，使系統前進到n步驟。

iterate.P<-function(x,P,n){res<-matrix(NA,n+1,lena<-xfor(iinseq_len(n))res[i+1,]<-x<-x%*%P res}

從處於狀態1的系統開始（x向量$ \[1,0,0\] $也是如此，表示處於狀態1的機率爲100％，而且不處於任何其餘狀態）

一樣，對於另外兩種可能的起始狀態：

y2<-iterate.P(c(0,1,0),P,n)y3<-iterate.P(c(0,0,1),P,n)

這代表了平穩分佈的收斂性。

ma=1,xlab="Step",ylab="y",las=1)matlines(0:n,y2,lty=2)matlines(0:n,y3,lty=3)

咱們可使用R的eigen函數來提取系統的主要特徵向量（t()這裏轉置矩陣以便獲得左特徵向量）。

v<-eigen(t(P)ars[,1]v<-v/sum(v)# normalise eigenvector

而後在以前的數字上加上一點，代表咱們有多接近收斂：

matplot(0:n,y1a3,lty=3)points(rep(10,3),v,col=1:3)

上面的過程迭代了不一樣狀態的整體機率; 而不是經過系統的實際轉換。因此，讓咱們迭代系統，而不是機率向量。

 run<-function(i,P,n){res<-integer(n)for(a(n))res[[t]]<-i<-sample(nrow(P),1,pr=P[i,]) res}

這鏈條運行了100個步驟：

 samples<-run(1,P,100)ploaes,type="s",xlab="Step",ylab="State",las=1)

繪製咱們在每一個狀態隨時間變化的時間分數，而不是繪製狀態：

 plot(cummean(samplesa2)lines(cummean(samples==3),col=3)

再運行一下（5000步）

n<-5000set.seed(1)samples<-run(1,P,n)plot(cummeanasamples==2),col=2)lines(cummean(samples==3),col=3)abline(h=v,lty=2,col=1:3)

因此這裏的關鍵是：馬爾可夫鏈是整潔和理解的東西，有一些不錯的屬性。馬爾可夫鏈有固定的分佈，若是咱們運行它們足夠長的時間，咱們能夠看看鏈條在哪裏花費時間，並對該平穩分佈進行合理的估計。

Metropolis算法

這是最簡單的MCMC算法。本節不打算展現如何設計高效的MCMC採樣器，而只是爲了看到他們確實在工做。

算法進行以下。

從$ x\_t $開始。

建議一個新的狀態$ x ^ \\ prime $

計算「接受機率」

從$ \[0,1\] $中抽取一些均勻分佈的隨機數$ u $; 若是$ u <\\ alpha $接受該點，則設置$ x {t + 1} = x ^ \\ prime $。不然拒絕它並設置$ x {t + 1} = x\_t $。

請注意，在上面的步驟3中，未知歸一化常數由於而退出

這將產生一系列樣本$ {x 0，x 1，\\ ldots} $。請注意，若是建議的樣本被拒絕，相同的值將出如今連續的樣本中。

還要注意，這些不是來自目標分佈的獨立樣本; 他們是依賴樣本 ; 也就是說，示例$ x\_t $取決於$ x\_ {t-1} $等等。然而，因爲鏈條接近平穩分佈，只要咱們抽取足夠的點數，這種依賴性就不會有問題。

MCMC採樣1d（單參數）問題

這是一個目標分配樣本。這是兩個正態分佈的加權和。這種分佈至關簡單，能夠從MCMC中抽取樣本。

至關隨意的，這裏是一些參數和目標密度的定義。

 p<-0.4ma1,2)sd<-c(.5,2)f<-function(x)p*dnora],sd[1])+(1-p)*dnorm(x,mu[2],sd[2])

機率密度繪製

咱們來定義一個很是簡單的提議算法，該算法從以當前點爲中心的標準誤差爲4的正態分佈中抽樣

而這隻須要運行MCMC的幾個步驟。它將從點x返回一個矩陣，其nsteps行數和列數與x元素的列數相同。若是在標量上運行， x它將返回一個向量。

 run<-funagth(x))for(iinseq_len(nsteps))res[i,]<-x<-step(x,f,q)drop(res)}

咱們選擇一個地方開始（如何-10，只是選擇一個很是糟糕的一點）

這裏是馬爾可夫鏈的前1000步，目標密度在右邊：

 layout(matrix(ca,type="s",xpd=NA,ylab="Parameter",xlab="Sample",las=1)usr<-par("usr")xx<-seq(usr[a4],length=301)plot(f(xx),xx,type="l",yaxs="i",axes=FALSE,xlab="")

即便只有一千個（非獨立的）樣本，咱們也開始至關相似於目標分佈。

hist(res,5aALSE,main="",ylim=c(0,.4),las=1,xlab="x",ylab="Probability density")z<-integrate(f,-Inf,Inf)$valuecurve(f(x)/z,add=TRUE,col="red",n=200)

運行更長時間，事情開始看起來更好：

res.long<-run(-10,f,q,50000)hist(res.long,100,freq=FALSE,main="",ylim=c(0,.4),las=1,xlab="x",ylabaadd=TRUE,col="red",n=200)

如今，運行不一樣的提案機制 - 一個標準差很大（33個單位），另外一個標準差很小（3個單位）。

res.fast<-run(-10action(x)rnorm(1,x,33),1000)res.slow<-run(-10,f,functanorm(1,x,.3),1000)

這裏是與上面相同的情節 - 注意三條軌跡正在移動的不一樣方式。

相反，紅色的痕跡（大的提案）正在提示可能性空間中的可怕空間，並拒絕其中的大部分空間。這意味着它每每會一次空間留下來。

藍色的蹤影提出了傾向於被接受的小動做，可是它隨着大部分的軌跡隨機行走。它須要數百次迭代才能達到機率密度的大部分。

您能夠在隨後的參數中看到不一樣提議步驟在自相關中的效果 - 這些圖顯示了不一樣延遲步驟之間自相關係數的衰減，藍線表示統計獨立性。

 par(mfrow=c(1,3ain="Intermediate")acf(res.fast,las=1,main="Large steps")

由此能夠計算獨立樣本的有效數量：

1coda::effectiveSize(res)1 2## var1 ## 1871coda::effectiveSize(res.fast)1 2## var1 ## 33.191coda::effectiveSize(res.slow)1 2## var1 ## 5.378

這更清楚地顯示了鏈條運行時間更長的狀況：

 naun(-10,f,q,n))xlim<-range(sapply(saa100)hh<-lapply(samples,function(x)hist(x,br,plot=FALSE))ylim<-c(0,max(f(xx)))

顯示100,1,000,10,000和100,000步：

 par(mfrow=c(2,2),mar=rep(.5,4),oma=c(4,4,0,0))for(hinhh){plot(h,main="",freq=a=300)}

MCMC在兩個維度

給出了一個多元正態密度，給定一個均值向量（分佈的中心）和方差 - 協方差矩陣。

 make.mvn<-function(mean,vcv){logdet<-as.numeric(detea+logdetvcv.i<-solve(vcv)function(x){dx<-x-meanexp(-(tmp+rowSums((dx%*%vcv.i)*dx))/2)}}

如上所述，將目標密度定義爲兩個mvns的總和（此次未加權）：

 mu1<-c(-1,1)mu2<-c(2,-2)vcv1<-ma5,.25,1.5),2,2)vcv2<-matrix(c(2,-.5,-.5,2aunctioax)+f2(x)x<-seq(-5,6,length=71)y<-seq(-7,6,lena-expand.grid(x=x,y=y)z<-matrix(aaTRUE)

從多元正態分佈取樣也至關簡單，但咱們將使用MCMC從中抽取樣本。

這裏有一些不一樣的策略 - 咱們能夠同時在兩個維度上提出動做，或者咱們能夠獨立地沿着每一個軸進行採樣。這兩種策略都能奏效，雖然它們的混合速度會有所不一樣。

假設咱們實際上並不知道如何從mvn中抽樣 ，讓咱們提出一個在兩個維度上一致的提案分佈，從每邊的寬度爲「d」的正方形取樣。

比較抽樣分佈與已知分佈：

例如，參數1 的邊際分佈是多少？

hisales[,1],freq=FALSa",xlab="x",ylab="Probability density")

咱們須要整合第一個參數的第二個參數的全部可能值。那麼，由於目標函數自己並非標準化的，因此咱們必須將其分解爲第一維積分值 。

 m<-function(x1){g<-Vectorize(function(x2)f(c(x1,ae(g,-Inf,Inf)$value}xx<-seq(mina]),max(sales[,1]),length=201)yy<-sapply(xx,m)z<-integrate(splinefun(xx,yy),min(xx),max(xx))$valuehist(samples[,1],freq=FALSE,ma,0.25))lines(xx,yy/z,col="red")

相關文章

R語言用Backfitting _MCMC_抽樣算法進行貝葉斯推理案例

 R語言中實現馬爾可夫鏈蒙特卡羅_MCMC_模型 

R語言中的Stan機率編程_MCMC_採樣的貝葉斯模型

R語言實現_MCMC_中的Metropolis–Hastings算法與吉布斯...

matlab實現_MCMC_的馬爾可夫切換ARMA - GARCH模型估計