?原文：http://tecdat.cn/?p=3772

?

創(chuàng)建測試數(shù)據(jù)

作為第一步，我們創(chuàng)建一些測試數(shù)據(jù)，用于擬合我們的模型。讓我們假設(shè)預(yù)測變量和響應(yīng)變量之間存在線性關(guān)系，因此我們采用線性模型并添加一些噪聲。

我將x值平衡在零附近以“去相關(guān)”斜率和截距。?

?

1. trueA <- 5

2. 


3. trueB <- 0

4. 


5. trueSd <- 10

6. 


7. sampleSize <- 31

8. 


9. 


10. 


11. # 創(chuàng)建獨(dú)立的x值

12. 


13. x <- (-(sampleSize-1)/2):((sampleSize-1)/2)

14. 


15. # 根據(jù)ax + b + N（0，sd）創(chuàng)建因變量

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

17. 


18. 


19. 


20. plot(x,y, main="Test Data")

圖

定義統(tǒng)計(jì)模型

下一步是指定統(tǒng)計(jì)模型。我們已經(jīng)知道數(shù)據(jù)是用x和y之間的線性關(guān)系y = a * x + b和帶有標(biāo)準(zhǔn)差sd的正常誤差模型N（0，sd）創(chuàng)建的，所以讓我們使用相同的模型進(jìn)行擬合，看看如果我們可以檢索我們的原始參數(shù)值。

從模型中導(dǎo)出似然函數(shù)

為了估計(jì)貝葉斯分析中的參數(shù)，我們需要導(dǎo)出我們想要擬合的模型的似然函數(shù)?？赡苄允俏覀兤谕^察到的數(shù)據(jù)以我們所看到的模型的參數(shù)為條件發(fā)生的概率（密度）。因此，假設(shè)我們的線性模型y = b + a * x + N（0，sd）將參數(shù)（a，b，sd）作為輸入，我們必須返回在此模型下獲得上述測試數(shù)據(jù)的概率（這聽起來更復(fù)雜，正如你在代碼中看到的那樣，我們只是計(jì)算預(yù)測y = b + a * x和觀察到的y之間的差異，然后我們必須查找概率密度（使用dnorm）發(fā)生這種偏差。

1. likelihood <- function(param){

2. 


3. ????a = param[1]

4. 


5. ????b = param[2]

6. 


7. ????sd = param[3]

8. 


9. 


10. 


11. ????pred = a*x + b

12. 


13. ? ? ?sumll = sum(singlelikelihoods)

14. 


15. ? ? ?(sumll)??

16. 


17. }

18. 


19. 


20. 


21. ?slopevalues <- function(x){return(likelihood(c(x, trueB, trueSd)))}

?

?

斜率參數(shù)的對(duì)數(shù)似然曲線

作為說明，代碼的最后幾行繪制了斜率參數(shù)a的一系列參數(shù)值的似然性。?

為什么我們使用對(duì)數(shù)

您可能已經(jīng)注意到我返回似然函數(shù)中概率的對(duì)數(shù)，這也是我對(duì)所有數(shù)據(jù)點(diǎn)的概率求和的原因（乘積的對(duì)數(shù)等于對(duì)數(shù)之和）。我們?yōu)槭裁匆鲞@個(gè)？ 因?yàn)楹芏嘈「怕食艘缘目赡苄院芸炀蜁?huì)變得非常?。ū热?0 ^ -34）。在某些階段，計(jì)算機(jī)程序正在進(jìn)入數(shù)字舍入問題。?

定義先驗(yàn)

作為第二步，與貝葉斯統(tǒng)計(jì)中一樣，我們必須為每個(gè)參數(shù)指定先驗(yàn)分布。為了方便起見，我對(duì)所有三個(gè)參數(shù)使用了均勻分布和正態(tài)分布。?無信息先驗(yàn)通常是1 / sigma的比例（如果你想了解原因，請(qǐng)看這里）。?

1. #先驗(yàn)分布

2. 


3. prior <- function(param){

4. 


5. ????a = param[1]

6. 


7. ????b = param[2]

8. 


9. ????sd = param[3]

10. 


11. ????aprior =??(a, min=0, max=10, log = T)

12. 


13. ????bprior = dnorm(b, sd = 5, log = T)

14. 


15. ?}

?后驗(yàn)

先驗(yàn)和可能性的乘積是MCMC將要處理的實(shí)際數(shù)量。這個(gè)函數(shù)被稱為后驗(yàn)?。同樣，在這里我們使用總和，因?yàn)槲覀兪褂脤?duì)數(shù)。

1. posterior <- function(param){

2. 


3. ???return (?(param) + prior(param))

4. 


5. }

?

?MCMC

現(xiàn)在，實(shí)際上是Metropolis-Hastings算法。該算法最常見的應(yīng)用之一（如本例所示）是從貝葉斯統(tǒng)計(jì)中的后驗(yàn)密度中提取樣本。然而，原則上，該算法可用于從任何可積函數(shù)中進(jìn)行采樣。因此，該算法的目的是在參數(shù)空間中跳轉(zhuǎn)，但是以某種方式使得在某一點(diǎn)上的概率與我們采樣的函數(shù)成比例（這通常稱為目標(biāo)函數(shù)）。在我們的例子中，這是上面定義的后驗(yàn)。

這是通過

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

2. 根據(jù)稱為提議函數(shù)的某個(gè)概率密度，選擇接近舊值的新參數(shù)值

3. 以概率p（新）/ p（舊）跳到這個(gè)新點(diǎn)，其中p是目標(biāo)函數(shù)，p> 1表示跳躍

當(dāng)我們運(yùn)行這個(gè)算法時(shí)，它訪問的參數(shù)的分布會(huì)收斂到目標(biāo)分布p。那么，讓我們在R中得到 ：

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

2. 


3. 


4. 


5. proposalfunction <- function(param){

6. 


7. ????return(rnorm(3,mean = param, sd= c(0.1,0.5,0.3)))

8. 


9. }

10. 


11. 


12. 


13. run_metropolis_MCMC <- function(startvalue, iterations){

14. 


15. ? ? ??for (i in 1:iterations){

16. 


17. 


18. 


19. ? ? ? ? ?if (runif(1) < probab){

20. 


21. ????????????chain[i+1,] = proposal

22. 


23. ????????}else{

24. 


25. ????????????chain[i+1,] = chain[i,]

26. 


27. ????????}

28. 


29. ????}

30. 


31. ????return(chain)

32. 


33. }

34. 


35. 


36. 


37. ?chain = run_metropolis_MCMC(startvalue, 10000)

38. 


39. 


40. 


41. burnIn = 5000

42. 


43. acceptance = 1-mean(duplicated(chain[-(1:burnIn),]))

使用后驗(yàn)的對(duì)數(shù)可能在開始時(shí)有點(diǎn)混亂，特別是當(dāng)您查看計(jì)算接受概率的行時(shí)（probab = exp（后驗(yàn)（建議） - 后驗(yàn)（鏈[i，]）） ）。要理解我們?yōu)槭裁催@樣做，請(qǐng)注意p1 / p2 = exp [log（p1）-log（p2）]。

算法的第一步可能受初始值的偏差，因此通常被丟棄用于進(jìn)一步分析 。要看的一個(gè)有趣的輸出是接受率： 接受標(biāo)準(zhǔn)拒絕提案的頻率是多少？接受率可以受提議函數(shù)的影響：通常，提案越接近，接受率越大。然而，非常高的接受率通常是無益的：這意味著算法“停留”在同一點(diǎn) ?？梢宰C明，20％到30％的接受率對(duì)于典型應(yīng)用來說是最佳的 。

?

1. ###概要#######################

2. 


3. 


4. 


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

6. 


7. hist( [-(1:burnIn),1],nclass=30, , main="Posterior of a", xlab="True value = red line" )

8. 


9. abline(v = mean(chain[-(1:burnIn),1]))

10. 


11. 


12. 


13. 


14. 


15. #進(jìn)行比較：

16. 


17. summary(lm(y~x))

生成的圖應(yīng)該類似于下圖。您可以看到我們或多或少地檢索了用于創(chuàng)建數(shù)據(jù)的原始參數(shù)，并且您還看到我們獲得了一個(gè)圍繞最高后驗(yàn)值的特定區(qū)域，這些區(qū)域也有一些數(shù)據(jù)支持，這是貝葉斯相當(dāng)于置信度的間隔。

圖：上排顯示斜率（a）的后驗(yàn)估計(jì)，截距（b）和誤差的標(biāo)準(zhǔn)偏差（sd）。下一行顯示馬爾可夫鏈參數(shù)值。

?

還有問題嗎？請(qǐng)?jiān)谙旅媪粞裕?/strong>

最受歡迎的見解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)

3.R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真

4.R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型

7.R語言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析

8.R語言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)

標(biāo)簽：