全文連結 http://tecdat.cn/?p=23255
最近我們被客戶要求撰寫關於rstan的研究報告，包括一些圖形和統計輸出。
本文將談論Stan以及如何在R中使用rstan創建Stan模型
儘管Stan提供了使用其程式語言的文檔和帶有例子的用戶指南，但對於初學者來說，這可能是很難理解的。
Stan
Stan是一種用於指定統計模型的程式語言。它最常被用作貝葉斯分析的MCMC採樣器。馬爾科夫鏈蒙特卡洛（MCMC）是一種抽樣方法，允許你在不知道分布的所有數學屬性的情況下估計一個機率分布。它在貝葉斯推斷中特別有用，因為後驗分布往往不能寫成表達式。要使用Stan，用戶要寫一個Stan程序，代表他們的統計模型。這個程序指定了模型中的參數和目標後驗密度。Stan代碼被編譯並與數據一起運行，輸出一組參數的後驗模擬。Stan與最流行的數據分析語言，如R、Python、shell、MATLAB、Julia和Stata的接口。我們將專注於在R中使用Stan。
rstan
rstan允許R用戶實現貝葉斯模型。你可以使用熟悉的公式和data.frame語法（如lm()）來擬合模型。通過為常用的模型類型提供預編譯的stan代碼來實現這種更簡單的語法。它使用起來很方便，但只限於特定的 "常用 "模型類型。如果你需要擬合不同的模型類型，那麼你需要自己用rstan編碼。
模型擬合函數以前綴stan_開始，以模型類型結束。建模函數有兩個必要的參數。
此外，還有一個可選的先驗參數，它允許你改變默認的先驗分布。
stan()函數讀取和編譯你的stan代碼，並在你的數據集上擬合模型。
stan()函數有兩個必要參數。
例子
作為一個簡單的例子來演示如何在這些包中指定一個模型，我們將使用汽車數據來擬合一個線性回歸模型。我們的因變量是mpg，所有其他變量是自變量。
mtcars %>%
head()
首先，我們將擬合模型。對於線性回歸，我們使用stan函數。
點擊標題查閱往期內容
R語言RStan貝葉斯示例：重複試驗模型和種群競爭模型Lotka Volterra
左右滑動查看更多
01
02
03
04
summary(fit)
輸出顯示參數摘要，包括平均值、標準差和量值。此外，它還顯示了MCMC的診斷統計Rhat和有效樣本量。這些統計數據對於評估MCMC算法是否收斂非常重要。
接下來，我們將用rstan來擬合同一個模型。下面是我們模型的stan代碼，保存在一個名為stan的文件中（你可以在RStudio中創建一個.stan文件，或者使用任何文本編輯器，並保存擴展名為.stan的文件）。
數據
int
int
matrix[N, K] X; // 預測矩陣
參數
real alpha; // 截距
模型
y ~ normal(alpha + X * beta, sigma); // 目標密度
Stan代碼在 "程序塊 "中結構化。每個Stan模型都需要三個程序塊，即數據、參數和模型。
數據塊是用來聲明作為數據讀入的變量的。在我們的例子中，我們有結果向量（y）和預測矩陣（X）。當把矩陣或向量聲明為一個變量時，你需要同時指定對象的維度。因此，我們還將讀出觀測值的數量（N）和預測器的數量（K）。
在參數塊中聲明的變量是將被Stan採樣的變量。在線性回歸的情況下，感興趣的參數是截距項（alpha）和預測因子的係數（beta）。此外，還有誤差項，sigma。
模型區塊是定義變量機率聲明的地方。在這裡，我們指定目標變量具有正態分布，其平均值為α+X*β，標準差為sigma。在這個塊中，你還可以指定參數的先驗分布。默認情況下，參數被賦予平坦的（非信息性）先驗。
此外，還有一些可選的程序塊：函數、轉換的數據、轉換的參數和生成的數量。
接下來，我們需要以Stan程序所期望的方式來格式化我們的數據。stan()函數要求將數據作為一個命名的列表傳入，其中的元素是你在數據塊中定義的變量。對於這個程序，我們創建一個元素為N、K、X和Y的列表。
list(
N = 32,
K = 10,
X = predictors,
y = mpg
)
現在我們已經準備好了我們的代碼和數據，我們把它們傳給函數來擬合模型。
fit_rstan
輸出類似的匯總統計數據，包括每個參數的平均值、標準偏差和量值。這些結果可能相似但不完全相同。它們之所以不同，是因為統計數據是根據後驗的隨機抽樣來計算的。
評估收斂性
當使用MCMC擬合一個模型時，檢查鏈是否收斂是很重要的。我們推薦可視化來直觀地檢查MCMC的診斷結果。我們將創建軌跡圖，Rhat值圖。
首先，讓我們創建軌跡圖。軌跡圖顯示了MCMC疊代過程中參數的採樣值。如果模型已經收斂，那麼軌跡圖應該看起來像一個圍繞平均值的隨機散點。如果鏈在參數空間中蜿蜒，或者鏈收斂到不同的值，那就證明有問題了。我們來演示。
mcmctrace()
這些軌跡圖表明，兩個模型都已經收斂了。對於所有的參數，四條鏈都是混合的，沒有明顯的趨勢。
接下來，我們將檢查Rhat值。Rhat是一種收斂診斷方法，它比較了各條鏈的參數估計值。如果鏈已經收斂並且混合良好，那麼Rhat值應該接近1。如果鏈沒有收斂到相同的值，那麼Rhat值將大於1。Rhat值為1.05或更高，表明存在收斂問題。rhat()函數需要一個Rhat值的向量作為輸入，所以我們首先提取Rhat值。
rhat() +
yaxis_text()
所有的Rhat值都低於1.05，說明沒有收斂問題。
Stan是一個建立貝葉斯模型的強大工具，這些包使R用戶可以很容易地使用Stan。
點擊文末 「閱讀原文」
獲取全文完整資料。
本文選自《R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數據和可視化診斷》。
點擊標題查閱往期內容
【視頻】馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現|數據分享
R語言實現MCMC中的Metropolis–Hastings算法與吉布斯採樣
R語言貝葉斯METROPOLIS-HASTINGS GIBBS 吉布斯採樣器估計變點指數分布分析泊松過程車站等待時間
R語言馬爾可夫MCMC中的METROPOLIS HASTINGS，MH算法抽樣（採樣）法可視化實例
python貝葉斯隨機過程：馬爾可夫鏈Markov-Chain，MC和Metropolis-Hastings，MH採樣算法可視化
Python貝葉斯推斷Metropolis-Hastings（M-H）MCMC採樣算法的實現
Metropolis Hastings採樣和貝葉斯泊松回歸Poisson模型
Matlab用BUGS馬爾可夫區制轉換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H採樣分析時間序列R語言RSTAN MCMC：NUTS採樣算法用LASSO 構建貝葉斯線性回歸模型分析職業聲望數據
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉換隨機波動率SV模型、粒子濾波、Metropolis Hasting採樣時間序列分析
R語言Metropolis Hastings採樣和貝葉斯泊松回歸Poisson模型
R語言貝葉斯MCMC：用rstan建立線性回歸模型分析汽車數據和可視化診斷
R語言貝葉斯MCMC：GLM邏輯回歸、Rstan線性回歸、Metropolis Hastings與Gibbs採樣算法實例
R語言貝葉斯Poisson泊松-正態分布模型分析職業足球比賽進球數
R語言用Rcpp加速Metropolis-Hastings抽樣估計貝葉斯邏輯回歸模型的參數
R語言邏輯回歸、Naive Bayes貝葉斯、決策樹、隨機森林算法預測心臟病
R語言中貝葉斯網絡（BN）、動態貝葉斯網絡、線性模型分析錯頜畸形數據
R語言中的block Gibbs吉布斯採樣貝葉斯多元線性回歸
Python貝葉斯回歸分析住房負擔能力數據集
R語言實現貝葉斯分位數回歸、lasso和自適應lasso貝葉斯分位數回歸分析
Python用PyMC3實現貝葉斯線性回歸模型
R語言用WinBUGS 軟體對學術能力測驗建立層次（分層）貝葉斯模型
R語言Gibbs抽樣的貝葉斯簡單線性回歸仿真分析
R語言和STAN,JAGS：用RSTAN,RJAG建立貝葉斯多元線性回歸預測選舉數據
R語言基於copula的貝葉斯分層混合模型的診斷準確性研究
R語言貝葉斯線性回歸和多元線性回歸構建工資預測模型
R語言貝葉斯推斷與MCMC：實現Metropolis-Hastings 採樣算法示例
R語言stan進行基於貝葉斯推斷的回歸模型
R語言中RStan貝葉斯層次模型分析示例
R語言使用Metropolis-Hastings採樣算法自適應貝葉斯估計與可視化
R語言隨機搜索變量選擇SSVS估計貝葉斯向量自回歸（BVAR）模型
WinBUGS對多元隨機波動率模型：貝葉斯估計與模型比較
R語言實現MCMC中的Metropolis–Hastings算法與吉布斯採樣
R語言貝葉斯推斷與MCMC：實現Metropolis-Hastings 採樣算法示例
R語言使用Metropolis-Hastings採樣算法自適應貝葉斯估計與可視化
視頻：R語言中的Stan機率編程MCMC採樣的貝葉斯模型
R語言MCMC:Metropolis-Hastings採樣用於回歸的貝葉斯估計