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

最近我們被客戶要求撰寫關于rstan的研究報告，包括一些圖形和統(tǒng)計輸出。

本文將談論Stan以及如何在R中使用rstan創(chuàng)建Stan模型

盡管Stan提供了使用其編程語言的文檔和帶有例子的用戶指南，但對于初學者來說，這可能是很難理解的。

Stan

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

rstan

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

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

• 公式。一個指定因變量和自變量的公式（y ~ x1 + x2）。

• data。一個包含公式中變量的數(shù)據(jù)框。

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

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

stan()函數(shù)有兩個必要參數(shù)。

• 文件。包含你的Stan程序的.stan文件的路徑。

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

例子

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

mtcars?%>% ??head()

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

點擊標題查閱往期內容

R語言RStan貝葉斯示例：重復試驗模型和種群競爭模型Lotka Volterra

左右滑動查看更多

01

02

03

04

summary(fit)

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

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

數(shù)據(jù) ??int<lower=0>?N;???//?觀測值的數(shù)量??int<lower=0>?K;???//?預測的數(shù)量??matrix[N,?K]?X;???//?預測矩陣... 參數(shù) ??real?alpha;???????????//?截距... 模型 ??y?~?normal(alpha?+?X?*?beta,?sigma);??//?目標密度

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

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

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

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

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

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

list( ????N?=?32, ????K?=?10, ????X?=?predictors, ????y?=?mpg ??)

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

fit_rstan

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

評估收斂性

當使用MCMC擬合一個模型時，檢查鏈是否收斂是很重要的。我們推薦可視化來直觀地檢查MCMC的診斷結果。我們將創(chuàng)建軌跡圖，Rhat值圖。

首先，讓我們創(chuàng)建軌跡圖。軌跡圖顯示了MCMC迭代過程中參數(shù)的采樣值。如果模型已經(jīng)收斂，那么軌跡圖應該看起來像一個圍繞平均值的隨機散點。如果鏈在參數(shù)空間中蜿蜒，或者鏈收斂到不同的值，那就證明有問題了。我們來演示。

??mcmctrace()

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

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

?rhat()??+ ??yaxis_text()

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

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

點擊文末?“閱讀原文”

獲取全文完整資料。

本文選自《R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷》。

點擊標題查閱往期內容

【視頻】馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現(xiàn)|數(shù)據(jù)分享
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯METROPOLIS-HASTINGS GIBBS 吉布斯采樣器估計變點指數(shù)分布分析泊松過程車站等待時間
R語言馬爾可夫MCMC中的METROPOLIS HASTINGS，MH算法抽樣（采樣）法可視化實例
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH采樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC采樣算法的實現(xiàn)
Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區(qū)制轉換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列R語言RSTAN MCMC：NUTS采樣算法用LASSO 構建貝葉斯線性回歸模型分析職業(yè)聲望數(shù)據(jù)
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉換隨機波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
R語言Metropolis Hastings采樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數(shù)據(jù)和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs采樣算法實例
R語言貝葉斯Poisson泊松-正態(tài)分布模型分析職業(yè)足球比賽進球數(shù)
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數(shù)
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言中貝葉斯網(wǎng)絡（BN）、動態(tài)貝葉斯網(wǎng)絡、線性模型分析錯頜畸形數(shù)據(jù)
R語言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數(shù)據(jù)集
R語言實現(xiàn)貝葉斯分位數(shù)回歸、lasso和自適應lasso貝葉斯分位數(shù)回歸分析
Python用PyMC3實現(xiàn)貝葉斯線性回歸模型
R語言用WinBUGS 軟件對學術能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數(shù)據(jù)
R語言基于copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言stan進行基于貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現(xiàn)MCMC中的Metropolis–Hastings算法與吉布斯采樣
R語言貝葉斯推斷與MCMC：實現(xiàn)Metropolis-Hastings 采樣算法示例
R語言使用Metropolis-Hastings采樣算法自適應貝葉斯估計與可視化
視頻：R語言中的Stan概率編程MCMC采樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings采樣用于回歸的貝葉斯估計