原文鏈接：http://tecdat.cn/?p=19664?

原文出處：拓端數(shù)據(jù)部落

MCMC是從復雜概率模型中采樣的通用技術。

1. 蒙特卡洛

2. 馬爾可夫鏈

3. Metropolis-Hastings算法

問題

如果需要計算有復雜后驗pdf p（θ| y）的隨機變量θ的函數(shù)f（θ）的平均值或期望值。

您可能需要計算后驗概率分布p（θ）的最大值。

解決期望值的一種方法是從p（θ）繪制N個隨機樣本，當N足夠大時，我們可以通過以下公式逼近期望值或最大值

將相同的策略應用于通過從p（θ| y）采樣并取樣本集中的最大值來找到argmaxp（θ| y）。

解決方法

1.1直接模擬

1.2逆CDF

1.3拒絕/接受抽樣

如果我們不知道精確/標準化的pdf或非常復雜，則MCMC會派上用場。

馬爾可夫鏈

為了模擬馬爾可夫鏈，我們必須制定一個?過渡核T（xi，xj）。過渡核是從狀態(tài)xi遷移到狀態(tài)xj的概率。

?馬爾可夫鏈的收斂性意味著它具有平穩(wěn)分布π。馬爾可夫鏈的統(tǒng)計分布是平穩(wěn)的,那么它意味著分布不會隨著時間的推移而改變。

Metropolis算法

?對于一個Markov鏈是平穩(wěn)的。基本上表示

處于狀態(tài)x并轉換為狀態(tài)x'的概率必須等于處于狀態(tài)x'并轉換為狀態(tài)x的概率

或者

方法是將轉換分為兩個子步驟；候選和接受拒絕。

令q（x'| x）表示?候選密度，我們可以使用概率?α（x'| x）來調整q? 。

候選分布?Q（X'| X）是給定的候選X的狀態(tài)X'的條件概率，

和?接受分布?α（x'| x）的條件概率接受候選的狀態(tài)X'-X'。我們設計了接受概率函數(shù)，以滿足詳細的平衡。

該?轉移概率?可以寫成：

插入上一個方程式，我們有

Metropolis-Hastings算法?

A的選擇遵循以下邏輯。

在q下從x到x'的轉移太頻繁了。因此，我們應該選擇α（x | x'）=1。但是，為了滿足?細致平穩(wěn)，我們有

下一步是選擇滿足上述條件的接受。Metropolis-Hastings是一種常見的?選擇：

即，當接受度大于1時，我們總是接受，而當接受度小于1時，我們將相應地拒絕。因此，Metropolis-Hastings算法包含以下內容：

1. 初始化：隨機選擇一個初始狀態(tài)x；

2. 根據(jù)q（x'| x）隨機選擇一個新狀態(tài)x';

3.接受根據(jù)α（x'| x）的狀態(tài)。如果不接受，則不會進行轉移，因此無需更新任何內容。否則，轉移為x'；

4.轉移到2，直到生成T狀態(tài)；

5.保存狀態(tài)x，執(zhí)行2。

原則上，我們從分布P（x）提取保存的狀態(tài)，因為步驟4保證它們是不相關的。必須根據(jù)候選分布等不同因素來選擇T的值。 重要的是，尚不清楚應該使用哪種分布q（x'| x）；必須針對當前的特定問題進行調整。

屬性

Metropolis-Hastings算法的一個有趣特性是它?僅取決于比率

是候選樣本x'與先前樣本xt之間的概率，

是兩個方向（從xt到x'，反之亦然）的候選密度之比。如果候選密度對稱，則等于1。

馬爾可夫鏈從任意初始值x0開始，并且算法運行多次迭代，直到“初始狀態(tài)”被“忘記”為止。這些被丟棄的樣本稱為預燒（burn-in）。其余的x可接受值集代表分布P（x）中的樣本

Metropolis采樣

一個簡單的Metropolis-Hastings采樣

讓我們看看從?伽瑪分布?模擬任意形狀和比例參數(shù)，使用具有Metropolis-Hastings采樣算法。

下面給出了Metropolis-Hastings采樣器的函數(shù)。該鏈初始化為零，并在每個階段都建議使用N（a / b，a /（b * b））個候選對象。

基于正態(tài)分布且均值和方差相同gamma的Metropolis-Hastings獨立采樣

1. 從某種狀態(tài)開始xt。代碼中的x。

2. 在代碼中提出一個新的狀態(tài)x'候選

3. 計算“接受概率”

1. 從[0,1] 得出一些均勻分布的隨機數(shù)u；如果u <α接受該點，則設置xt + 1 = x'。否則，拒絕它并設置xt + 1 = xt。

MH可視化

1. set.seed(123)

2. 


3. for (i in 2:n) {

4. can <- rnorm(1, mu, sig)

5. aprob <- min(1, (dgamma(can, a, b)/dgamma(x,

6. a, b))/(dnorm(can, mu, sig)/dnorm(x,

7. mu, sig)))

8. u <- runif(1)

9. if (u < aprob)

10. x <- can

11. vec[i] <- x

畫圖

設置參數(shù)。

1. nrep<- 54000

2. burnin<- 4000

3. shape<- 2.5

4. rate<-2.6

修改圖，僅包含預燒期后的鏈

1. vec=vec[-(1:burnin)]

2. #vec=vec[burnin:length(vec)]

1. par(mfrow=c(2,1)) # 更改主框架，在一幀中有多少個圖形

2. plot(ts(vec), xlab="Chain", ylab="Draws")

3. abline(h = mean(vec), lwd="2", col="red" )

1. Min. ?1st Qu. ? Median ? ? Mean ?3rd Qu. ? ? Max.

2. 0.007013 0.435600 0.724800 0.843300 1.133000 3.149000

var(vec[-(1:burnin)])[1] 0.2976507

初始值

第一個樣本?vec?是我們鏈的初始/起始值。我們可以更改它，以查看收斂是否發(fā)生了變化。

1. x <- 3*a/b

2. vec[1] <- x

選擇方案

如果候選密度與目標分布P（x）的形狀匹配，即q（x'| xt）≈P（x'）q（x'|），則該算法效果最佳。 xt）≈P（x'）。如果使用正態(tài)候選密度q，則在預燒期間必須調整方差參數(shù)σ2。

通常，這是通過計算接受率來完成的，接受率是在最后N個樣本的窗口中接受的候選樣本的比例。

如果σ2太大，則接受率將非常低，因為候選可能落在概率密度低得多的區(qū)域中，因此a1將非常小，且鏈將收斂得非常慢。

示例2：回歸的貝葉斯估計

Metropolis-Hastings采樣用于貝葉斯估計回歸模型。

設定參數(shù)

DGP和圖

1. # 創(chuàng)建獨立的x值，大約為零

2. x <- (-(Size-1)/2):((Size-1)/2)

3. # 根據(jù)ax + b + N（0，sd）創(chuàng)建相關值

4. y <- ?trueA * x + trueB + rnorm(n=Size,mean=0,sd=trueSd)

正態(tài)分布擬然

1. 


2. pred = a*x + b

3. singlelikelihoods = dnorm(y, mean = pred, sd = sd, log = T)

4. sumll = sum(singlelikelihoods)

為什么使用對數(shù)

似然函數(shù)中概率的對數(shù)，這也是我求和所有數(shù)據(jù)點的概率（乘積的對數(shù)等于對數(shù)之和）的原因。

我們?yōu)槭裁匆鲞@個？強烈建議這樣做，因為許多小概率相乘的概率會變得很小。在某個階段，計算機程序會陷入數(shù)值四舍五入或下溢問題。

因此，?當您編寫概率時，請始終使用對數(shù)

示例：繪制斜率a的似然曲線

1. # 示例：繪制斜率a的似然曲線

2. plot (seq(3, 7, by=.05), slopelikelihoods , type="l")

先驗分布

這三個參數(shù)的均勻分布和正態(tài)分布。

1. # 先驗分布

2. 


3. # 更改優(yōu)先級，log為True，因此這些均為log

4. density/likelihood

5. aprior = dunif(a, min=0, max=10, log = T)

6. bprior = dnorm(b, sd = 2, log = T)

7. sdprior = dunif(sd, min=0, max=30, log = T)

后驗

先驗和概率的乘積是MCMC將要處理的實際量。此函數(shù)稱為后驗函數(shù)。同樣，這里我們使用和，因為我們使用對數(shù)。

1. posterior <- function(param){

2. return (likelihood(param) + prior(param))

3. }

Metropolis算法

該算法是從?后驗密度中采樣最常見的貝葉斯統(tǒng)計應用之一?。

上面定義的后驗。

1. 從隨機參數(shù)值開始

2. 根據(jù)某個候選函數(shù)的概率密度，選擇一個接近舊值的新參數(shù)值

3. 以概率p（new）/ p（old）跳到這個新點，其中p是目標函數(shù)，并且p> 1也意味著跳躍

4. 請注意，我們有一個?對稱的跳躍/候選分布?q（x'| x）。

標準差σ是固定的。

所以接受概率等于

1. ######## Metropolis 算法 ################

2. 


3. 


4. for (i in 1:iterations){

5. 


6. probab = exp(posterior(proposal) - posterior(chain[i,]))

7. if (runif(1) < probab){

8. chain[i+1,] = proposal

9. }else{

10. chain[i+1,] = chain[i,]

11. }

實施

（e）輸出接受的值，并解釋。

1. 


2. chain = metrMCMC(startvalue, 5500)

3. 


4. burnIn = 5000

5. accep = 1-mean(duplicated(chain[-(1:burnIn),]))

算法的第一步可能會因初始值而有偏差，因此通常會被丟棄來進行進一步分析（預燒期）。令人感興趣的輸出是接受率：候選多久被算法接受拒絕一次？候選函數(shù)會影響接受率：通常，候選越接近，接受率就越大。但是，非常高的接受率通常是無益的：這意味著算法在同一點上“停留”，這導致對參數(shù)空間（混合）的處理不夠理想。

我們還可以更改初始值，以查看其是否更改結果/是否收斂。

startvalue = c(4,0,10)

小結

1. V1 ? ? ? ? ? ? ?V2 ? ? ? ? ? ? ? ?V3

2. Min. ? :4.068 ? Min. ? :-6.7072 ? Min. ? : 6.787

3. 1st Qu.:4.913 ? 1st Qu.:-2.6973 ? 1st Qu.: 9.323

4. Median :5.052 ? Median :-1.7551 ? Median :10.178

5. Mean ? :5.052 ? Mean ? :-1.7377 ? Mean ? :10.385

6. 3rd Qu.:5.193 ? 3rd Qu.:-0.8134 ? 3rd Qu.:11.166

7. Max. ? :5.989 ? Max. ? : 4.8425 ? Max. ? :19.223

1. #比較:

2. summary(lm(y~x))

1. 
   

2. Call:

3. lm(formula = y ~ x)

4. 
   

5. Residuals:

6. Min ? ? ?1Q ?Median ? ? ?3Q ? ? Max

7. -22.259 ?-6.032 ?-1.718 ? 6.955 ?19.892

8. 
   

9. Coefficients:

10. Estimate Std. Error t value Pr(>|t|)

11. (Intercept) ?-3.1756 ? ? 1.7566 ?-1.808 ? ?0.081 .

12. x ? ? ? ? ? ? 5.0469 ? ? 0.1964 ?25.697 ? <2e-16 ***

13. ---

14. Signif. codes: ?0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

15. 
    

16. Residual standard error: 9.78 on 29 degrees of freedom

17. Multiple R-squared: ?0.9579, ? ?Adjusted R-squared: ?0.9565

18. F-statistic: 660.4 on 1 and 29 DF, ?p-value: < 2.2e-16

summary(lm(y~x))$sigma[1] 9.780494coefficients(lm(y~x))[1]

1. (Intercept)

2. -3.175555

coefficients(lm(y~x))[2]

1. x

2. 5.046873

總結：

1. ### 總結: #######################

2. 


3. par(mfrow = c(2,3))

4. hist(chain[-(1:burnIn),1],prob=TRUE,nclass=30,col="109"

5. abline(v = mean(chain[-(1:burnIn),1]), lwd="2")

最受歡迎的見解

1.用R語言模擬混合制排隊隨機服務排隊系統(tǒng)

2.R語言中使用排隊論預測等待時間

3.R語言中實現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型

4.R語言中的馬爾科夫機制轉換(Markov regime switching)模型

5.matlab貝葉斯隱馬爾可夫hmm模型

6.用R語言模擬混合制排隊隨機服務排隊系統(tǒng)

7.Python基于粒子群優(yōu)化的投資組合優(yōu)化

8.R語言馬爾可夫轉換模型研究交通傷亡人數(shù)事故預測

9.用機器學習識別不斷變化的股市狀況——隱馬爾可夫模型的應用