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

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

?

現(xiàn)在有了對(duì)貝葉斯方法的概念理解，我們將實(shí)際研究使用它的回歸模型。為了簡(jiǎn)單起見(jiàn)，我們從回歸的標(biāo)準(zhǔn)線性模型開(kāi)始。然后添加對(duì)采樣分布或先驗(yàn)的更改。我們將通過(guò) R 和相關(guān)的 R 包 rstan 使用編程語(yǔ)言 Stan。

示例：線性回歸模型

在下文中，我們將設(shè)置一些初始數(shù)據(jù)，并使用標(biāo)準(zhǔn) lm 函數(shù)運(yùn)行模型比較。

設(shè)置

首先，我們需要?jiǎng)?chuàng)建我們將在此處使用的數(shù)據(jù)以及本文檔中的大多數(shù)其他示例。

1. 


2. # 設(shè)置可復(fù)制種子

3. set.seed(85309)

4. 


5. 


6. # 運(yùn)行 lm 以供稍后比較； 但如果需要，請(qǐng)立即檢查

7. modlm = lm(y~., data=data.frame)

此時(shí)我們有三個(gè)協(xié)變量和一個(gè) y，它是正態(tài)分布線性函數(shù)，標(biāo)準(zhǔn)差等于 2。系數(shù)的總體值包括截距分別為 5、0.2、-1.5 和 0.9，盡管添加了噪聲，但樣本的實(shí)際估計(jì)值略有不同。現(xiàn)在我們準(zhǔn)備好為輸入到 Stan 的數(shù)據(jù)設(shè)置一個(gè) R 列表對(duì)象，以及對(duì)這些數(shù)據(jù)進(jìn)行建模的相應(yīng) Stan 代碼。

我將展示在 R 中通過(guò)單個(gè)字符串實(shí)現(xiàn)的所有 Stan 代碼，然后提供每個(gè)相應(yīng)模型塊的一些細(xì)節(jié)。但是，這里的目標(biāo)不是專注于工具，而是專注于概念。

Stan 的數(shù)據(jù)列表應(yīng)包括 Stan 代碼中可能使用的任何矩陣、向量或值。例如，與數(shù)據(jù)一起，可以包括樣本大小、組指標(biāo)（例如混合模型）等。在這里，我們可以只使用樣本大小 (N)、模型矩陣中的列數(shù) (K)、目標(biāo)變量 (y) 和模型矩陣 (X)。

1. # 為stan輸入創(chuàng)建數(shù)據(jù)列表對(duì)象

2. dat = list

接下來(lái)是 Stan 代碼。在 R2OpenBugs 或 rjags 中，可以使用代碼調(diào)用單獨(dú)的文本文件，并且可以對(duì) rstan 執(zhí)行相同操作，但出于我們的目的，我們?cè)?R 代碼中顯示它。首先要注意的是模型代碼。接下來(lái)，Stan 有必須按順序調(diào)用的編程塊。我將在代碼中列出所有塊來(lái)記錄它們的順序并依次討論每個(gè)塊，即使我們不會(huì)全部使用它們。// 或 # 之后或 /* */ 之間的任何內(nèi)容都是與代碼相關(guān)的注釋。而分布用 ～～ 指定，例如，?y ~ normal(0, 1)?表示 y 正態(tài)分布，均值為 0，標(biāo)準(zhǔn)差為 1。

此外，要安裝 rstan，需要通過(guò) CRAN 或 GitHub。它不需要單獨(dú)安裝 Stan 本身，但它確實(shí)需要幾個(gè)步驟并且需要 C++ 編譯器。一旦你安裝了 rstan，它就會(huì)像任何其他 R 包一樣被調(diào)用。

1. # 使用語(yǔ)法創(chuàng)建模型對(duì)象

2. stanmodelcode = "

3. data { // 數(shù)據(jù)塊

4. int<lower=1> N; // 樣本大小

5. int<lower=1> K; // 模型矩陣的尺寸

7. /*

8. 轉(zhuǎn)換后的數(shù)據(jù) { // 轉(zhuǎn)換后的數(shù)據(jù)塊。目前不使用。

9. }

10. */

12. 參數(shù)// 參數(shù)塊

13. real<lower=0> sigma; // 誤差比例

14. }

16. 模型 // 模型塊

17. mu = X * beta; // 創(chuàng)建線性預(yù)測(cè)器

18. 


19. // 先驗(yàn)指標(biāo)

20. sigma ~ cauchy(0, 5); // 由于sigma被限定在0

21. 


22. // 似然

23. y ~ normal(mu, sigma);

25. /*

26. 生成的數(shù)量 { // 生成的數(shù)量塊

Stan代碼

第一部分是數(shù)據(jù)塊，我們告訴 Stan 它應(yīng)該從數(shù)據(jù)列表中獲得的數(shù)據(jù)。設(shè)置邊界作為對(duì)數(shù)據(jù)輸入的檢查，這就是?< >?。聲明的前兩個(gè)變量是?N?和?K，都是整數(shù)。接下來(lái)代碼分別聲明模型矩陣和目標(biāo)向量。我們聲明變量的類型和維度，然后聲明其名稱。在 Stan 中，在一個(gè)塊中聲明的所有內(nèi)容都可用于后續(xù)塊，但在一個(gè)塊中聲明的內(nèi)容可能不會(huì)在更早的塊中使用。即使在一個(gè)塊中，任何聲明的東西，例如?N?和?K, 然后可以隨后使用，就像我們指定模型矩陣的維度一樣?X。

作為參考，以下內(nèi)容來(lái)注明了感興趣的變量以及將在其中聲明它們的相關(guān)塊。

種類聲明塊建模的，未建模的數(shù)據(jù)數(shù)據(jù)，轉(zhuǎn)換后的數(shù)據(jù)建模參數(shù)，缺失數(shù)據(jù)參數(shù)，轉(zhuǎn)換參數(shù)未建模參數(shù)數(shù)據(jù)，轉(zhuǎn)換后的數(shù)據(jù)產(chǎn)生量轉(zhuǎn)換數(shù)據(jù)、轉(zhuǎn)換參數(shù)、生成量循環(huán)索引循環(huán)語(yǔ)句

轉(zhuǎn)換后的數(shù)據(jù)塊是您可以執(zhí)行對(duì)數(shù)或中心化等類似操作的地方，即您可以根據(jù)輸入數(shù)據(jù)或一般情況下創(chuàng)建新數(shù)據(jù)。但是，如果您使用的是 R，那么首先在 R 中執(zhí)行這些操作并將它們包含在數(shù)據(jù)列表中。您還可以在此處聲明任何未建模的參數(shù)，例如您希望固定為某個(gè)值的參數(shù)。

要估計(jì)的主要感興趣的參數(shù)位于參數(shù)塊中。與數(shù)據(jù)塊一樣，您只能聲明這些變量，不能進(jìn)行任何賦值。在這里，我們注意到要估計(jì)的 β?和 σ，后者的下限為零。在實(shí)踐中，如果截距或其他系數(shù)在顯著不同的尺度上，您可能更愿意將它們分開(kāi)建模。

轉(zhuǎn)換后的參數(shù)塊是可能包含可選參數(shù)的地方。為了提高效率，您通常只想放置依賴于參數(shù)塊的特定興趣的東西。

模型塊是指定您的先驗(yàn)和可能性以及任何必要變量的聲明的地方。例如，此處包含線性預(yù)測(cè)器，因?yàn)樗鼘②呄蛴谒迫? 請(qǐng)注意，我們可以將線性預(yù)測(cè)器放在轉(zhuǎn)換后的參數(shù)部分，但這會(huì)減慢過(guò)程，而且我們對(duì)這些特定值不太感興趣。

我對(duì)系數(shù)使用的是正態(tài)先驗(yàn)，平均值為零，標(biāo)準(zhǔn)差很大。對(duì)于σ的估計(jì)，我使用的是Cauchy?分布。許多使用BUGS的回歸例子都會(huì)使用反伽馬先驗(yàn)，這對(duì)這個(gè)模型來(lái)說(shuō)是完全可以的，盡管它對(duì)其他方差參數(shù)的效果并不理想。如果我們沒(méi)有為參數(shù)的先驗(yàn)分布指定任何東西，均勻分布將是默認(rèn)的。

最后，你想計(jì)算的任何東西都可以放在這里--對(duì)新數(shù)據(jù)的預(yù)測(cè)、參數(shù)的比率、參數(shù)大于x的次數(shù)、為報(bào)告目的的參數(shù)轉(zhuǎn)換等等。我們稍后將對(duì)此進(jìn)行演示。

運(yùn)行模型

現(xiàn)在我們對(duì)代碼的作用有了一個(gè)概念。貝葉斯估計(jì)，像最大似然法一樣，以初始猜測(cè)為起點(diǎn)，然后以迭代的方式運(yùn)行，每一步都從后驗(yàn)分布中產(chǎn)生模擬抽樣，然后糾正這些抽樣，直到最后達(dá)到某個(gè)目標(biāo)，或平穩(wěn)分布。這一部分是關(guān)鍵，與經(jīng)典的統(tǒng)計(jì)學(xué)不同。我們的目標(biāo)是一個(gè)分布，而不是一個(gè)點(diǎn)估計(jì)。

這個(gè)模擬過(guò)程被稱為馬爾科夫鏈蒙特卡洛，簡(jiǎn)稱MCMC。這個(gè)過(guò)程的具體細(xì)節(jié)使許多貝葉斯編程語(yǔ)言/方法與眾不同。在MCMC中，所有來(lái)自后驗(yàn)的模擬抽樣都是基于以前的抽樣并與之相關(guān)的，因?yàn)檫@個(gè)過(guò)程是沿著走向平穩(wěn)分布的道路前進(jìn)的。我們通常會(huì)讓這個(gè)過(guò)程預(yù)燒，或者說(shuō)從最初的起點(diǎn)開(kāi)始有點(diǎn)穩(wěn)定下來(lái)，這可能會(huì)有很大的偏差，因此在最初的幾次迭代中，后續(xù)的估計(jì)值也會(huì)有很大的偏差。然而，作為進(jìn)一步的檢查，我們將多次運(yùn)行整個(gè)過(guò)程，也就是說(shuō)，有一個(gè)以上的鏈。由于這些鏈將從不同的地方開(kāi)始，如果多個(gè)鏈最后都到達(dá)了同一個(gè)結(jié)果，我們就可以對(duì)結(jié)果更有信心。

你會(huì)注意到Stan將其代碼編譯為C++的時(shí)間可能比運(yùn)行模型的時(shí)間要長(zhǎng)，而在我的電腦上，每條鏈只需要一秒鐘多一點(diǎn)的時(shí)間。然而，貝葉斯方法曾經(jīng)需要很長(zhǎng)的時(shí)間，即使是像這樣的標(biāo)準(zhǔn)回歸，這也許是貝葉斯分析在過(guò)去幾十年里才流行起來(lái)的主要原因；我們根本沒(méi)有機(jī)器來(lái)有效地做這件事。即使是現(xiàn)在，對(duì)于高度復(fù)雜的模型和大數(shù)據(jù)集來(lái)說(shuō)，它仍然需要很長(zhǎng)的時(shí)間來(lái)運(yùn)行，盡管通常不會(huì)太長(zhǎng)。

在下面的代碼中，我們注意到包含stan模型代碼的對(duì)象，數(shù)據(jù)列表，我們想要多少次迭代（5000），我們想要這個(gè)過(guò)程在開(kāi)始保留任何估計(jì)值之前運(yùn)行多長(zhǎng)時(shí)間（warmup=2500），我們想要保留多少次后驗(yàn)的抽?。╰hin=10意味著每十次抽?。?，以及鏈的數(shù)量（chains=4）。最后，我們將有四條鏈，從參數(shù)的后驗(yàn)分布中抽取1000次。即使在verbose = FALSE的情況下，Stan也會(huì)向R控制臺(tái)輸出大量的輸出，我在這里省略，但是你會(huì)看到一些關(guān)于編譯過(guò)程的初始信息，當(dāng)每條鏈通過(guò)iter參數(shù)中指定的10%的迭代時(shí)的更新，以及最后對(duì)耗時(shí)的估計(jì)。你可能還會(huì)看到一些信息，除非它們是高度重復(fù)的，否則不應(yīng)該被視為錯(cuò)誤。

1. library(rstan)

2. 


3. ### 運(yùn)行模型并檢查結(jié)果

4. fit = stan(

5. 


6. warmup = 2500,

7. chains = 4)

隨著模型的運(yùn)行，我們現(xiàn)在可以檢查結(jié)果。在下文中，我們指定要顯示的數(shù)字精度，我們想要哪些參數(shù)，以及我們想要哪些后驗(yàn)抽樣的量級(jí)，在本例中是中位數(shù)和那些會(huì)產(chǎn)生95%區(qū)間估計(jì)的參數(shù)。

1. # 摘要

2. print(fit

到目前為止還不錯(cuò)。平均估計(jì)值反映了感興趣的參數(shù)的后驗(yàn)結(jié)果的平均值，是標(biāo)準(zhǔn)回歸分析中報(bào)告的典型系數(shù)。值得注意的是95%的概率或置信區(qū)間，因?yàn)樗鼈儾皇悄闼赖闹眯艆^(qū)間。這里沒(méi)有重復(fù)抽樣的解釋。概率區(qū)間是更直觀的。它的意思很簡(jiǎn)單，根據(jù)這個(gè)模型的結(jié)果，真實(shí)值有95%的可能性會(huì)落在這兩點(diǎn)之間。

將這些結(jié)果與R的lm函數(shù)的結(jié)果相比較，我們可以看到我們得到了類似的估計(jì)值，因?yàn)樗鼈冊(cè)谛?shù)點(diǎn)后兩位是相同的。事實(shí)上，如果我們使用統(tǒng)一先驗(yàn)，模型與用標(biāo)準(zhǔn)最大似然估計(jì)所做的基本相同。在這里，對(duì)于一個(gè)并不復(fù)雜的模型來(lái)說(shuō)，我們有相當(dāng)多的數(shù)據(jù)，所以我們會(huì)期望可能性明顯超過(guò)先驗(yàn)。
?

summary

但是我們?cè)趺粗牢覀兊哪Ｐ褪欠襁\(yùn)作良好呢？有幾個(gè)標(biāo)準(zhǔn)的診斷方法，但讓我們看一下目前的一些情況。在摘要中，se_mean是蒙特卡洛誤差，是對(duì)只有有限數(shù)量的后驗(yàn)抽樣所帶來(lái)的不確定性的估計(jì)。n_eff是給定所有鏈的有效樣本量，基本上占了鏈中的自相關(guān)，即當(dāng)我們從一次抽樣到下一次抽樣時(shí)估計(jì)的相關(guān)性。它實(shí)際上不需要很大，但如果它相對(duì)于所需的總抽樣數(shù)來(lái)說(shuō)很小，那就可能引起關(guān)注了。 Rhat是衡量鏈的混合程度的指標(biāo)，當(dāng)鏈被允許運(yùn)行無(wú)限次抽簽時(shí)，它就會(huì)變成1。在這種情況下，n_eff和Rhat表明我們有很好的收斂性，但我們也可以用跟蹤圖來(lái)直觀地檢查。

1. # 視覺(jué)化

2. srace(fit'))

我們可以看到跟蹤圖中，每條鏈的估計(jì)值都能很快地從起點(diǎn)找到一個(gè)或多或少的穩(wěn)定狀態(tài)（灰色的初始預(yù)燒迭代）。此外，所有三條鏈（每條鏈都用不同的顏色表示）都混合得很好，并在同一結(jié)論附近跳動(dòng)。

stan開(kāi)發(fā)團(tuán)隊(duì)通過(guò)shinystan包使交互式探索診斷變得很容易。此外，coda包中還有其他診斷方法，Stan模型的結(jié)果可以很容易地轉(zhuǎn)換為與之配合。下面的代碼演示了如何開(kāi)始。

1. 


2. bets = extract$beta

所以，你已經(jīng)有了。除了制作數(shù)據(jù)列表和產(chǎn)生特定語(yǔ)言的模型代碼的初始設(shè)置之外，相對(duì)于標(biāo)準(zhǔn)模型，運(yùn)行貝葉斯回歸模型并不一定需要太多的時(shí)間。最主要的也許是采用不同的思維方式，并從這個(gè)新的角度來(lái)解釋結(jié)果。對(duì)于你所熟悉的標(biāo)準(zhǔn)模型，可能不需要太長(zhǎng)的時(shí)間就能像使用那些模型一樣自如，現(xiàn)在你將有一個(gè)靈活的工具，帶你進(jìn)入新的領(lǐng)域，加深理解。

最受歡迎的見(jiàn)解

1.使用R語(yǔ)言進(jìn)行METROPLIS-IN-GIBBS采樣和MCMC運(yùn)行

2.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型

3.R語(yǔ)言實(shí)現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣

4.R語(yǔ)言BUGS JAGS貝葉斯分析 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣

5.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

6.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真分析

7.R語(yǔ)言用Rcpp加速M(fèi)etropolis-Hastings抽樣估計(jì)貝葉斯邏輯回歸模型的參數(shù)

8.R語(yǔ)言使用Metropolis- Hasting抽樣算法進(jìn)行邏輯回歸

9.R語(yǔ)言中基于混合數(shù)據(jù)抽樣(MIDAS)回歸的HAR-RV模型預(yù)測(cè)GDP增長(zhǎng)