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

像任何統(tǒng)計(jì)建模一樣，貝葉斯建?？赡苄枰獮槟愕难芯繂栴}設(shè)計(jì)合適的模型，然后開發(fā)該模型，使其符合你的數(shù)據(jù)假設(shè)并運(yùn)行（點(diǎn)擊文末“閱讀原文”獲取完整代碼數(shù)據(jù)）。

1. 了解?Stan

統(tǒng)計(jì)模型可以在R或其他統(tǒng)計(jì)語(yǔ)言的各種包中進(jìn)行擬合。但有時(shí)你在概念上可以設(shè)計(jì)的完美模型，在限制了你可以使用的分布和復(fù)雜性的軟件包或程序中很難或不可能實(shí)現(xiàn)。這時(shí)你可能想轉(zhuǎn)而使用統(tǒng)計(jì)編程語(yǔ)言，如Stan。

Stan是一種新式的語(yǔ)言，它提供了一種更全面的學(xué)習(xí)和實(shí)現(xiàn)貝葉斯模型的方法，可以適應(yīng)復(fù)雜的數(shù)據(jù)結(jié)構(gòu)。Stan開發(fā)團(tuán)隊(duì)的一個(gè)目標(biāo)是通過清晰的語(yǔ)法、更好的采樣器（這里的采樣是指從貝葉斯后驗(yàn)分布中抽取樣本）以及與許多平臺(tái)（包括R、RStudio、ggplot2和Shiny）的集成，使貝葉斯建模更易于使用。

在這個(gè)入門教程中，我們將從一個(gè)線性模型開始，經(jīng)歷模型建立的迭代過程。在我們的高級(jí)stan教程中，我們將探索更復(fù)雜的模型結(jié)構(gòu)。

首先，在建立模型之前，你需要定義你的問題并了解你的數(shù)據(jù)。探索它們，繪制它們，計(jì)算一些匯總統(tǒng)計(jì)。

一旦你對(duì)你的數(shù)據(jù)和你想用統(tǒng)計(jì)模型回答的問題有了了解，你就可以開始建立貝葉斯模型的迭代過程。

• 設(shè)計(jì)你的模型。

• 選擇先驗(yàn)

• 對(duì)后驗(yàn)分布進(jìn)行采樣。

• 檢查模型收斂（traceplots、rhats ）

• 使用后驗(yàn)預(yù)測(cè)批判性地評(píng)估模型并檢查它們與您的數(shù)據(jù)的比較情況

• 重復(fù)…

模擬數(shù)據(jù)也是很好的做法，以確保你的模型正確，作為測(cè)試你的模型的另一種方式。

2. 數(shù)據(jù)

首先，讓我們找到一個(gè)可以擬合簡(jiǎn)單線性模型的數(shù)據(jù)集。?氣候變化對(duì)地球最顯著的影響之一是北半球每年海冰范圍的減少。讓我們使用 Stan 的線性模型探索海冰范圍如何隨時(shí)間變化。

通過運(yùn)行setwd("your-file-path")?包含您自己的文件路徑的代碼?，將您的工作目錄設(shè)置為您保存數(shù)據(jù)的文件夾??，F(xiàn)在，讓我們加載數(shù)據(jù)：

#?添加stringsAsFactors?=?F意味著數(shù)字變量將不會(huì)被 #?作為因子/分類變量讀入 ece?<-?red.cv("sv",?stinsAsFators?=?F)

我們來看一下數(shù)據(jù)：

我們可以用這些數(shù)據(jù)提出什么研究問題？以下情況如何：

研究問題：?北半球的海冰范圍是否會(huì)隨著時(shí)間的推移而減少？

為了探索這個(gè)問題的答案，首先我們可以做一個(gè)數(shù)字。

plot(?th?~?yr,?data)

圖 1. 北半球海冰范圍隨時(shí)間的變化。

現(xiàn)在，讓我們使用?lm().

l1?<-?lm(exnoh?~?yer,?data?=?sie) summary(l1)

我們可以將該模型添加到我們的繪圖中：

ablne(m1,?l?=?2,?ty?=?2,?w?=?3)

圖 2. 北半球海冰范圍隨時(shí)間的變化（加上線性模型擬合）。

記住線性模型的方程：

y = α + β ? x + 誤差

在?Stan?你需要指定你想模型。? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

也許我們已經(jīng)找到了問題的答案，但本教程的重點(diǎn)是探索使用編程語(yǔ)言?Stan，所以現(xiàn)在讓我們嘗試在 Stan 中編寫相同的模型。

準(zhǔn)備數(shù)據(jù)

讓我們重命名變量并將年份從 1 索引到 39。關(guān)于貝葉斯模型的一個(gè)關(guān)鍵是您必須使用信息分布來描述數(shù)據(jù)中的變化。因此，您希望確保您的數(shù)據(jù)符合這些分布，并且它們將適用于您的模型。在這種情況下，我們真的想知道從數(shù)據(jù)集的開始到數(shù)據(jù)集結(jié)束的海冰是否發(fā)生了變化，而不是 1979 年到 2017 年。我們不需要我們的模型估計(jì) 500 年或 600 年的海冰是什么樣的，就在我們的數(shù)據(jù)集的持續(xù)時(shí)間內(nèi)。因此，我們將年份數(shù)據(jù)設(shè)置為索引 1 到 30 年。

x?<-?I(year?-?1978)

我們可以使用新數(shù)據(jù)重新運(yùn)行該線性模型。

summary(lm1)

我們還可以從我們的簡(jiǎn)單模型中提取一些關(guān)鍵的匯總統(tǒng)計(jì)數(shù)據(jù)，以便我們Stan?稍后可以將它們與模型的輸出進(jìn)行比較?。

coeff\[1\]?#?截距值 coeff\[2\]?#?斜率 sigma(lm1)?#?殘差

現(xiàn)在讓我們將其轉(zhuǎn)換為用于輸入Stan?模型的數(shù)據(jù)框?。傳遞給 Stan 的數(shù)據(jù)需要是命名對(duì)象列表。此處給出的名稱需要與模型中使用的變量名稱相匹配。

庫(kù)
-

請(qǐng)確保安裝了以下庫(kù)（這些是本Stan?教程和下一個(gè)教程的庫(kù)?）。?rstan?是最重要的，如果您沒有 C++ 編譯器，則需要一些額外的東西。

3.?Stan?程序

我們將首先用語(yǔ)言編寫一個(gè)線性模型?Stan。這可以寫在你的 R 腳本中，或者單獨(dú)保存為一個(gè)?.stan?文件并調(diào)用到?R.

一個(gè)?Stan?程序具有三個(gè)必需的“塊”：

• “數(shù)據(jù)”?塊：您可以在其中聲明數(shù)據(jù)類型、它們的維度、任何限制（即 upper = 或 lower = ，用作檢查?Stan）及其名稱。

• “參數(shù)”?塊：您可以在此處指明要建模的參數(shù)和名稱。對(duì)于線性回歸，我們希望對(duì)回歸線周圍的誤差的截距、任何斜率和標(biāo)準(zhǔn)偏差進(jìn)行建模。

• “模型”?塊：這是包含任何抽樣語(yǔ)句的地方，包括正在使用的模型。模型塊是指明要為參數(shù)包含的任何先驗(yàn)分布的地方。如果未定義Stan?先驗(yàn)，則?使用默認(rèn)先驗(yàn)?uniform(-infinity, +infinity)。您可以在聲明參數(shù)時(shí)使用上限或下限來限制先驗(yàn)（即?lower = 0\> 確保參數(shù)為正）。

采樣由?~?符號(hào)表示，并且?Stan?已經(jīng)包含許多常見的分布作為矢量化函數(shù)。

還有四個(gè)可選塊：

• “功能”

• "轉(zhuǎn)化的數(shù)據(jù)"

• "轉(zhuǎn)換后的參數(shù)

• "生成的數(shù)量"

注釋//?在 Stan中用 表示?。該write("model code", "file_name")?允許我們?cè)?R 腳本中編寫 Stan 模型并將文件輸出到工作目錄（或者您可以設(shè)置不同的文件路徑）。

write("//簡(jiǎn)單線性回歸的模型 數(shù)據(jù) ?int?<?lower?=?1?>?N;?//?樣本大小 vector\[N\]?x;//?預(yù)測(cè) vecor\[N\]?y;//?結(jié)果 參數(shù)? ?real?alha;?//?截距 ?real?beta;?//?斜率(回歸系數(shù)) ?real?<?loer?=?0?>?sima;?//?誤差SD 模型? ?y?~?nrmal(alpha?+?x?*?beta?,?siga); 產(chǎn)生的數(shù)量? //?后驗(yàn)預(yù)測(cè)分布"?。 "md1.stan"

首先，我們應(yīng)該檢查我們的?Stan?模型以確保我們編寫了一個(gè)文件。

現(xiàn)在讓我們保存該文件路徑。

?"model1.stan"

在這里，我們隱式地將?uniform(-infinity, +infinity)?先驗(yàn)用于我們的參數(shù)。

4. 運(yùn)行Stan?模型

Stan?程序C++?在被使用之前被遵守?。這意味著在 R 可以使用模型之前需要運(yùn)行 C++ 代碼。為此，您必須?C++?安裝編譯器。編譯后，您可以在每個(gè)會(huì)話中多次使用模型，但在開始新R?會(huì)話時(shí)必須重新編譯?。有許多?C++?編譯器，而且它們?cè)诓煌到y(tǒng)中通常是不同的。如果您的模型一堆錯(cuò)誤，請(qǐng)不要擔(dān)心。只要模型可以與該stan()?函數(shù)一起使用?，它就可以正確編譯。如果我們想使用以前編寫的?.stan?文件，我們?cè)?code>file?函數(shù)中使用?參數(shù)?stan_model()?。

我們通過使用stan()?函數(shù)擬合我們的模型?，并為它提供模型、數(shù)據(jù)，并指示預(yù)熱的迭代次數(shù)（這些迭代稍后不會(huì)用于后驗(yàn)分布，因?yàn)樗鼈冎皇悄Ｐ汀邦A(yù)熱” ”），總迭代次數(shù)，我們要運(yùn)行的鏈數(shù)，我們要使用的內(nèi)核數(shù)（Stan?為并行化而設(shè)置），它表示同時(shí)運(yùn)行的鏈數(shù)（即，如果您的計(jì)算機(jī)有四個(gè)內(nèi)核） ，您可以在每個(gè)鏈上運(yùn)行一個(gè)鏈，同時(shí)創(chuàng)建四個(gè)鏈）和細(xì)化，這是我們想要存儲(chǔ)我們的預(yù)熱后迭代的頻率?！皌hin = 1”將保留每次迭代，“thin = 2”將保留每一秒，依此類推……

Stan?如果warmup =?未指定參數(shù)，則自動(dòng)使用一半的迭代作為預(yù)熱?。

stan(fie?=moel1,?data?=?data,?wrmup?=?500,?ier?=?1000,?chais?=?4,?cres?=?2,?thn?=?1)

訪問fit?對(duì)象的內(nèi)容?

結(jié)果?stan()?保存為?stanfit?對(duì)象（S4 類）。

我們可以通過執(zhí)行對(duì)象的名稱來獲取參數(shù)估計(jì)和采樣器診斷的匯總統(tǒng)計(jì)信息：

fit

模型輸出展示了什么？你怎么知道你的模型已經(jīng)收斂了？您能看到指示您的 C++ 編譯器已運(yùn)行的文本嗎？

從這個(gè)輸出中，我們可以通過查看Rhat?每個(gè)參數(shù)的值來快速評(píng)估模型收斂性?。當(dāng)這些值等于或接近 1 時(shí)，鏈已經(jīng)收斂。還有許多其他診斷方法，但這對(duì) Stan 來說很重要。

我們還可以通過從模型對(duì)象中提取參數(shù)來查看參數(shù)的完整后驗(yàn)。有很多方法可以查看后驗(yàn)。

poteir?<-?exrat(fit)

extract()?將每個(gè)參數(shù)的后驗(yàn)估計(jì)放入一個(gè)列表中。

讓我們與我們之前使用“l(fā)m”的估計(jì)進(jìn)行比較：

plot(y?~?x)

圖 3. 北半球海冰范圍隨時(shí)間的變化（比較?Stan?線性模型擬合和一般?lm?擬合）。

結(jié)果與lm?輸出相同?。這是因?yàn)槲覀兪褂昧艘粋€(gè)簡(jiǎn)單的模型，并且在我們的參數(shù)上放置了非信息先驗(yàn)。

將回歸線估計(jì)中的可變性可視化的一種方法是繪制來自后驗(yàn)的多個(gè)估計(jì)。

plot(y?~?x,?pch?=?20)

圖 4. 北半球海冰范圍隨時(shí)間的變化（Stan?線性模型擬合）。

5. 改變先驗(yàn)

我們?cè)僭囈淮?，但現(xiàn)在對(duì)海冰和時(shí)間之間的關(guān)系采用更有信息量的預(yù)設(shè)。我們將使用小標(biāo)準(zhǔn)差的正態(tài)先驗(yàn)。如果我們使用標(biāo)準(zhǔn)差非常大的正態(tài)先驗(yàn)（比如1000，或者10000），它們的作用與均勻先驗(yàn)非常相似。

參數(shù)? ?real?alha;?//?截距 ?real?bta;?//?斜率(回歸系數(shù)) ?real?<?lowr?=?0?>?sima;?//?誤差SD 模型 ?beta?~?nomal(1,?0.1); ?y?~?normal(apha?+?x?*?beta?,?siga);

我們將擬合該模型并將其與使用均勻先驗(yàn)的均值估計(jì)進(jìn)行比較。

?stan(stn_oel) plot(y?~?x)

圖 5. 北半球海冰范圍隨時(shí)間的變化（Stan?線性模型擬合）。

后驗(yàn)預(yù)測(cè)發(fā)生了什么變化？模型是否更好地?cái)M合數(shù)據(jù)？為什么模型擬合發(fā)生了變化？通過制作非常窄的先驗(yàn)分布，我們的模型改變了什么？

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

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

轉(zhuǎn)存失敗重新上傳取消

左右滑動(dòng)查看更多

轉(zhuǎn)存失敗重新上傳取消

01

02

03

04

?

嘗試自己將先驗(yàn)更改為一些不同的數(shù)字，看看會(huì)發(fā)生什么，這是貝葉斯建模中的一個(gè)常見問題，如果您的先驗(yàn)分布非常窄，但不符合您對(duì)系統(tǒng)或數(shù)據(jù)分布的理解，您可以運(yùn)行無法有意義地解釋數(shù)據(jù)變化的模型。然而，這并不是說您不應(yīng)該選擇一些信息豐富的先驗(yàn)，您確實(shí)希望使用先前的分析和對(duì)您的研究系統(tǒng)的理解來告知您的模型先驗(yàn)和設(shè)計(jì)。您只需要仔細(xì)考慮您做出的每個(gè)建模決策！

6. 收斂診斷

在繼續(xù)之前，我們應(yīng)該再次檢查模型參數(shù)的?Rhat?值、有效樣本大小 (?n_eff) 和跟蹤圖，以確保模型已收斂且可靠。

n_f?是有效樣本大小的粗略度量。通常只需要擔(dān)心這個(gè)數(shù)字小于迭代次數(shù)的 1/100 或 1/1000。

對(duì)于跟蹤圖，我們可以直接從后驗(yàn)查看它們：

plot(alpha,?type?=?"l") plot(beta,?type?=?"l") plot(sigma,?type?=?"l")

圖 6. alpha 的跡線圖，截距。

對(duì)于更簡(jiǎn)單的模型，收斂通常不是問題，除非您的代碼中有錯(cuò)誤，或者運(yùn)行采樣器的迭代次數(shù)太少。

收斂性差

嘗試僅運(yùn)行 50 次迭代的模型并檢查跟蹤圖。

fid?<-?stan(?iter?=?50)

這在預(yù)熱后也有一些“轉(zhuǎn)換”，表明模型指定錯(cuò)誤，或者采樣器未能完全采樣后驗(yàn)。

plot(alpha,?type?=?"l") plot(beta,?type?=?"l") plot(sigma,?type?=?"l")

圖 7. alpha 的錯(cuò)誤軌跡圖，截距。

參數(shù)匯總

我們也可以直接通過后驗(yàn)得到參數(shù)的匯總。讓我們還繪制非貝葉斯線性模型值，以確保我們的模型運(yùn)行正確

par(mfrow?=?c(1,3)) plot(dnsty(alpha)

圖 8.Stan?模型擬合的密度圖分布?與一般lm?擬合的估計(jì)值比較?。

從后驗(yàn)我們可以直接計(jì)算任何參數(shù)超過或低于某個(gè)感興趣值的概率。

beta > 0 的概率：

sum(beta>0)/lenth(beta) #?0

beta > 0.2 的概率：

sum(bea>0.2)/legth(beta) #?0

診斷圖?rstan

雖然我們可以直接使用后驗(yàn)，但?rstan?內(nèi)置了許多有用的功能。

plot(fit)

圖 9.Stan?模型不同鏈的跟蹤圖?。

我們還可以查看后驗(yàn)密度和直方圖。

dens(it) hist(ft)

圖 10.Stan?模型中截距、斜率和殘差方差的后驗(yàn)密度圖和直方圖?。

我們可以生成指示平均參數(shù)估計(jì)值和我們可能感興趣的任何置信區(qū)間的圖。請(qǐng)注意，beta?和?sigma?參數(shù)的 95% 置信區(qū)間?非常小，因此您只能看到點(diǎn)。根據(jù)您自己數(shù)據(jù)的差異，當(dāng)您進(jìn)行自己的分析時(shí)，您可能會(huì)看到更小或更大的置信區(qū)間。

plot(fit)

圖 11.Stan?模型的參數(shù)估計(jì)?。

后驗(yàn)預(yù)測(cè)檢查

對(duì)于預(yù)測(cè)和作為模型診斷的另一種形式，?Stan?可以使用隨機(jī)數(shù)生成器在每次迭代中為每個(gè)數(shù)據(jù)點(diǎn)生成預(yù)測(cè)值。通過這種方式，我們可以生成預(yù)測(cè)，這些預(yù)測(cè)也代表了我們模型和數(shù)據(jù)生成過程中的不確定性。可用于獲取我們想要的關(guān)于后驗(yàn)的任何其他信息，或?qū)π聰?shù)據(jù)進(jìn)行預(yù)測(cè)。

參數(shù) ?real?alpha;?//?截距 ?real?beta;?//?斜率(回歸系數(shù)) ?y?~?nomal(x?*?eta?+?alpa,?sgma); 產(chǎn)生的數(shù)量? ?for?(n?in?1:N)?{ ?yp\[n\]?= normlrg(x\[n\]?* bta + apha, sima)?。 ?}

請(qǐng)注意，GQ（生成量）塊不支持矢量化，因此我們必須將其放入循環(huán)中。但是由于它被編譯為?C++，循環(huán)實(shí)際上非?？?，并且 Stan 每次迭代只評(píng)估一次 GQ 塊，因此它不會(huì)為您的采樣增加太多時(shí)間。通常，數(shù)據(jù)生成函數(shù)將是您在模型塊中使用的分布，但帶有?_rng?后綴。

?stan(modl2GQ,?data?,?ier?=?1000,?hans?=?4,?cres?=?2,?tin?=?1)

y_rep?從后驗(yàn)中提取?值。

處理y_rep?值有很多選擇?。

每一行都是模型的一次迭代（單一后驗(yàn)估計(jì)）。

我們可以制作一些更漂亮的圖。這個(gè)包是ggplot2。

在200次后驗(yàn)抽樣中，比較y的密度和y的密度。

poy(y,?yrep\[1:200,?\])

圖 12. 比較隨機(jī)后驗(yàn)抽取的估計(jì)值。

在這里，我們看到數(shù)據(jù)（深藍(lán)色）與我們的后驗(yàn)預(yù)測(cè)非常吻合。

我們還可以使用它來比較匯總統(tǒng)計(jì)的估計(jì)值。

pp(y?=?y,?yep?=?yrep,?tat?=?"mean")

圖 13. 比較匯總統(tǒng)計(jì)的估計(jì)值。

我們可以更改傳遞給?stat?函數(shù)的函數(shù)，甚至可以自己編寫！

我們可以調(diào)查每個(gè)數(shù)據(jù)點(diǎn)的平均后驗(yàn)預(yù)測(cè)與每個(gè)數(shù)據(jù)點(diǎn)的觀察值（默認(rèn)線為 1:1）

scttrg(y?=?y,?yrp?=?yrep)

圖 14. 每個(gè)數(shù)據(jù)點(diǎn)的平均后驗(yàn)預(yù)測(cè)與每個(gè)數(shù)據(jù)點(diǎn)的觀測(cè)值。

所以現(xiàn)在您已經(jīng)學(xué)習(xí)了如何運(yùn)行線性模型?Stan?并檢查模型收斂性。

如有任何問題，請(qǐng)聯(lián)系我們！

點(diǎn)擊文末“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語(yǔ)言STAN貝葉斯線性回歸模型分析氣候變化影響北半球海冰范圍和可視化檢查模型收斂性》。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

R語(yǔ)言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷

R語(yǔ)言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實(shí)例

R語(yǔ)言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進(jìn)球數(shù)

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

R語(yǔ)言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機(jī)森林算法預(yù)測(cè)心臟病

R語(yǔ)言中貝葉斯網(wǎng)絡(luò)（BN）、動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)、線性模型分析錯(cuò)頜畸形數(shù)據(jù)

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

Python貝葉斯回歸分析住房負(fù)擔(dān)能力數(shù)據(jù)集

R語(yǔ)言實(shí)現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應(yīng)lasso貝葉斯分位數(shù)回歸分析

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

R語(yǔ)言用WinBUGS 軟件對(duì)學(xué)術(shù)能力測(cè)驗(yàn)建立層次（分層）貝葉斯模型

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

R語(yǔ)言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預(yù)測(cè)選舉數(shù)據(jù)

R語(yǔ)言基于copula的貝葉斯分層混合模型的診斷準(zhǔn)確性研究

R語(yǔ)言貝葉斯線性回歸和多元線性回歸構(gòu)建工資預(yù)測(cè)模型

R語(yǔ)言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例

R語(yǔ)言stan進(jìn)行基于貝葉斯推斷的回歸模型

R語(yǔ)言中RStan貝葉斯層次模型分析示例

R語(yǔ)言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化

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

WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較

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

R語(yǔ)言貝葉斯推斷與MCMC：實(shí)現(xiàn)Metropolis-Hastings 采樣算法示例

R語(yǔ)言使用Metropolis-Hastings采樣算法自適應(yīng)貝葉斯估計(jì)與可視化

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

R語(yǔ)言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計(jì)