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

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

本文將談?wù)揝tan以及如何在R中使用rstan創(chuàng)建Stan模型。盡管Stan提供了使用其編程語言的文檔和帶有例子的用戶指南，但對于初學(xué)者來說，這可能是很難理解的。

Stan

Stan是一種用于指定統(tǒng)計模型的編程語言。它最常被用作貝葉斯分析的MCMC采樣器。馬爾科夫鏈蒙特卡洛（MCMC）是一種抽樣方法，允許你在不知道分布的所有數(shù)學(xué)屬性的情況下估計一個概率分布。它在貝葉斯推斷中特別有用，因為后驗分布往往不能寫成表達(dá)式。要使用Stan，用戶要寫一個Stan程序，代表他們的統(tǒng)計模型。這個程序指定了模型中的參數(shù)和目標(biāo)后驗密度。Stan代碼被編譯并與數(shù)據(jù)一起運(yùn)行，輸出一組參數(shù)的后驗?zāi)M。Stan與最流行的數(shù)據(jù)分析語言，如R、Python、shell、MATLAB、Julia和Stata的接口。我們將專注于在R中使用Stan。

rstan

rstan允許R用戶實(shí)現(xiàn)貝葉斯模型。你可以使用熟悉的公式和data.frame語法（如lm()）來擬合模型。通過為常用的模型類型提供預(yù)編譯的stan代碼來實(shí)現(xiàn)這種更簡單的語法。它使用起來很方便，但只限于特定的 "常用 "模型類型。如果你需要擬合不同的模型類型，那么你需要自己用rstan編碼。

模型擬合函數(shù)以前綴stan_開始，以模型類型結(jié)束。建模函數(shù)有兩個必要的參數(shù)。

- 公式。一個指定因變量和自變量的公式（y ~ x1 + x2）。
- data。一個包含公式中變量的數(shù)據(jù)框。

此外，還有一個可選的先驗參數(shù)，它允許你改變默認(rèn)的先驗分布。

stan()函數(shù)讀取和編譯你的stan代碼，并在你的數(shù)據(jù)集上擬合模型。

stan()函數(shù)有兩個必要參數(shù)。
- 文件。包含你的Stan程序的.stan文件的路徑。

- data。一個命名的列表，提供模型的數(shù)據(jù)。
?

例子

作為一個簡單的例子來演示如何在這些包中指定一個模型，我們將使用汽車數(shù)據(jù)來擬合一個線性回歸模型。我們的因變量是mpg，所有其他變量是自變量。

1. 


2. mtcars %>%

3. head()

首先，我們將擬合模型。對于線性回歸，我們使用stan函數(shù)。?

summary(fit)

輸出顯示參數(shù)摘要，包括平均值、標(biāo)準(zhǔn)差和量值。此外，它還顯示了MCMC的診斷統(tǒng)計Rhat和有效樣本量。這些統(tǒng)計數(shù)據(jù)對于評估MCMC算法是否收斂非常重要。

接下來，我們將用rstan來擬合同一個模型。下面是我們模型的stan代碼，保存在一個名為stan的文件中（你可以在RStudio中創(chuàng)建一個.stan文件，或者使用任何文本編輯器，并保存擴(kuò)展名為.stan的文件）。

1. 


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

3. int<lower=0> N; ? // 觀測值的數(shù)量

4. int<lower=0> K; ? // 預(yù)測的數(shù)量

5. matrix[N, K] X; ? // 預(yù)測矩陣

6. ...

7. 參數(shù)

8. real alpha; ? ? ? ? ? // 截距

9. ...

10. 模型

11. y ~ normal(alpha + X * beta, sigma); ?// 目標(biāo)密度

12. 


Stan代碼在 "程序塊 "中結(jié)構(gòu)化。每個Stan模型都需要三個程序塊，即數(shù)據(jù)、參數(shù)和模型。

數(shù)據(jù)塊是用來聲明作為數(shù)據(jù)讀入的變量的。在我們的例子中，我們有結(jié)果向量（y）和預(yù)測矩陣（X）。當(dāng)把矩陣或向量聲明為一個變量時，你需要同時指定對象的維度。因此，我們還將讀出觀測值的數(shù)量（N）和預(yù)測器的數(shù)量（K）。

在參數(shù)塊中聲明的變量是將被Stan采樣的變量。在線性回歸的情況下，感興趣的參數(shù)是截距項（alpha）和預(yù)測因子的系數(shù)（beta）。此外，還有誤差項，sigma。

模型區(qū)塊是定義變量概率聲明的地方。在這里，我們指定目標(biāo)變量具有正態(tài)分布，其平均值為α+X*β，標(biāo)準(zhǔn)差為sigma。在這個塊中，你還可以指定參數(shù)的先驗分布。默認(rèn)情況下，參數(shù)被賦予平坦的（非信息性）先驗。

此外，還有一些可選的程序塊：函數(shù)、轉(zhuǎn)換的數(shù)據(jù)、轉(zhuǎn)換的參數(shù)和生成的數(shù)量。

接下來，我們需要以Stan程序所期望的方式來格式化我們的數(shù)據(jù)。stan()函數(shù)要求將數(shù)據(jù)作為一個命名的列表傳入，其中的元素是你在數(shù)據(jù)塊中定義的變量。對于這個程序，我們創(chuàng)建一個元素為N、K、X和Y的列表。

1. list(

2. N = 32,

3. K = 10,

4. X = predictors,

5. y = mpg

6. )

現(xiàn)在我們已經(jīng)準(zhǔn)備好了我們的代碼和數(shù)據(jù)，我們把它們傳給函數(shù)來擬合模型。

fit_rstan

輸出類似的匯總統(tǒng)計數(shù)據(jù)，包括每個參數(shù)的平均值、標(biāo)準(zhǔn)偏差和量值。這些結(jié)果可能相似但不完全相同。它們之所以不同，是因為統(tǒng)計數(shù)據(jù)是根據(jù)后驗的隨機(jī)抽樣來計算的。

評估收斂性

當(dāng)使用MCMC擬合一個模型時，檢查鏈?zhǔn)欠袷諗渴呛苤匾?。我們推薦可視化來直觀地檢查MCMC的診斷結(jié)果。我們將創(chuàng)建軌跡圖，Rhat值圖。

首先，讓我們創(chuàng)建軌跡圖。軌跡圖顯示了MCMC迭代過程中參數(shù)的采樣值。如果模型已經(jīng)收斂，那么軌跡圖應(yīng)該看起來像一個圍繞平均值的隨機(jī)散點(diǎn)。如果鏈在參數(shù)空間中蜿蜒，或者鏈?zhǔn)諗康讲煌闹担蔷妥C明有問題了。我們來演示。?

1. 


2. mcmctrace()

這些軌跡圖表明，兩個模型都已經(jīng)收斂了。對于所有的參數(shù)，四條鏈都是混合的，沒有明顯的趨勢。

接下來，我們將檢查Rhat值。Rhat是一種收斂診斷方法，它比較了各條鏈的參數(shù)估計值。如果鏈已經(jīng)收斂并且混合良好，那么Rhat值應(yīng)該接近1。如果鏈沒有收斂到相同的值，那么Rhat值將大于1。Rhat值為1.05或更高，表明存在收斂問題。rhat()函數(shù)需要一個Rhat值的向量作為輸入，所以我們首先提取Rhat值。

1. rhat() ?+

2. yaxis_text()

所有的Rhat值都低于1.05，說明沒有收斂問題。

Stan是一個建立貝葉斯模型的強(qiáng)大工具，這些包使R用戶可以很容易地使用Stan。

最受歡迎的見解

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

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

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

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

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

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

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

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

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