# 使用語法創建模型對象

stanmodelcode = "

data { // 數據塊

int N; // 樣本大小

int K; // 模型矩陣的尺寸

/*

轉換後的數據 { // 轉換後的數據塊。目前不使用。

}

*/

參數// 參數塊

real sigma; // 誤差比例

}

模型 // 模型塊

mu = X * beta; // 創建線性預測器

// 先驗指標

sigma ~ cauchy(0, 5); // 由於sigma被限定在0

// 似然

y ~ normal(mu, sigma);

/*

生成量 { // 生成量塊

Stan代碼

第一部分是數據塊，我們告訴 Stan 它應該從數據列表中獲得的數據。設置邊界作為對數據輸入的檢查，這就是 < > 。聲明的前兩個變量是 N 和 K，都是整數。接下來代碼分別聲明模型矩陣和目標向量。我們聲明變量的類型和維度，然後聲明其名稱。在 Stan 中，在一個塊中聲明的所有內容都可用於後續塊，但在一個塊中聲明的內容不會在更早的塊中使用，例如聲明 N 和 K, 然後可以隨後使用，就像我們指定模型矩陣的維度一樣 X。

作為參考，以下內容來註明了感興趣的變量以及將在其中聲明它們的相關塊。

種類聲明塊建模的，未建模的數據數據，轉換後的數據建模參數，缺失數據參數，轉換參數未建模參數數據，轉換後的數據產生量轉換數據、轉換參數、生成量循環索引循環語句轉換後的數據塊是您可以執行對數或中心化等類似操作的地方，即您可以根據輸入數據或一般情況下創建新數據。但是，如果您使用的是 R，那麼首先在 R 中執行這些操作並將它們包含在數據列表中。您還可以在此處聲明任何未建模的參數，例如您希望固定為某個值的參數。

點擊標題查閱往期內容

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

左右滑動查看更多

01

02

03

04

要估計的主要感興趣的參數位於參數塊中。與數據塊一樣，您只能聲明這些變量，不能進行任何賦值。在這裡，我們注意到要估計的 β 和 σ，後者的下限為零。在實踐中，如果截距或其他係數在顯著不同的尺度上，您可能更願意將它們分開建模。

轉換後的參數塊是可能包含可選參數的地方。為了提高效率，您通常只想放置依賴於參數塊的特定興趣的東西。

模型塊是指定您的先驗和可能性以及任何必要變量的聲明的地方。例如，此處包含線性預測器，因為它將趨向於似然. 請注意，我們可以將線性預測器放在轉換後的參數部分，但這會減慢過程，而且我們對這些特定值不太感興趣。

我對係數使用的是正態先驗，平均值為零，標準差很大。對於σ的估計，我使用的是Cauchy 分布。許多使用BUGS的回歸例子都會使用反伽馬先驗，這對這個模型來說是完全可以的，儘管它對其他方差參數的效果並不理想。如果我們沒有為參數的先驗分布指定任何東西，均勻分布是默認的。

最後，你想計算的任何東西都可以放在這裡--對新數據的預測、參數的比率、參數大於x的次數、為報告目的的參數轉換等等。我們稍後將對此進行演示。

運行模型

現在我們對代碼的作用有了一個概念。貝葉斯估計，像最大似然法一樣，以初始猜測為起點，然後以疊代的方式運行，每一步都從後驗分布中產生模擬抽樣，然後糾正這些抽樣，直到最後達到某個目標，或平穩分布。這一部分是關鍵，與經典的統計學不同。我們的目標是一個分布，而不是一個點估計。

這個模擬過程被稱為馬爾科夫鏈蒙特卡洛，簡稱MCMC。這個過程的具體細節使許多貝葉斯程式語言/方法與眾不同。在MCMC中，所有來自後驗的模擬抽樣都是基於以前的抽樣並與之相關的，因為這個過程是沿著走向平穩分布的道路前進的。我們通常會讓這個過程預燒，或者說從最初的起點開始有點穩定下來，這可能會有很大的偏差，因此在最初的幾次疊代中，後續的估計值也會有很大的偏差。然而，作為進一步的檢查，我們將多次運行整個過程，也就是說，有一個以上的鏈。由於這些鏈將從不同的地方開始，如果多個鏈最後都到達了同一個結果，我們就可以對結果更有信心。

你會注意到Stan將其代碼編譯為C++的時間可能比運行模型的時間要長，而在我的電腦上，每條鏈只需要一秒鐘多一點的時間。然而，貝葉斯方法曾經需要很長的時間，即使是像這樣的標準回歸，這也許是貝葉斯分析在過去幾十年里才流行起來的主要原因；我們根本沒有機器來有效地做這件事。即使是現在，對於高度複雜的模型和大數據集來說，它仍然需要很長的時間來運行，儘管通常不會太長。

在下面的代碼中，我們注意到包含stan模型代碼的對象，數據列表，我們想要多少次疊代（5000），我們想要這個過程在開始保留任何估計值之前運行多長時間（warmup=2500），我們想要保留多少次後驗的抽取（thin=10意味著每十次抽取），以及鏈的數量（chains=4）。最後，我們將有四條鏈，從參數的後驗分布中抽取1000次。即使在verbose = FALSE的情況下，Stan也會向R控制台輸出大量的輸出，我在這裡省略，但是你會看到一些關於編譯過程的初始信息，當每條鏈通過iter參數中指定的10%的疊代時的更新，以及最後對耗時的估計。

library(rstan)

### 運行模型並檢查結果

fit (

warmup = 2500,

chains = 4)

隨著模型的運行，我們現在可以檢查結果。在下文中，我們指定要顯示的數字精度，我們想要哪些參數，以及我們想要哪些後驗抽樣的量級，在本例中是中位數和那些會產生95%區間估計的參數。

# 摘要

print(fit

到目前為止還不錯。平均估計值反映了感興趣的參數的後驗結果的平均值，是標準回歸分析中報告的典型係數。值得注意的是95%的機率或置信區間，因為它們不是你所知道的置信區間。這裡沒有重複抽樣的解釋。機率區間是更直觀的。它的意思很簡單，根據這個模型的結果，真實值有95%的可能性會落在這兩點之間。

將這些結果與R的lm函數的結果相比較，我們可以看到我們得到了類似的估計值，因為它們在小數點後兩位是相同的。事實上，如果我們使用統一先驗，模型與用標準最大似然估計所做的基本相同。

summary

但是我們怎麼知道我們的模型是否運作良好呢？有幾個標準的診斷方法，但讓我們看一下目前的一些情況。在摘要中，se_mean是蒙特卡洛誤差，是對只有有限數量的後驗抽樣所帶來的不確定性的估計。n_eff是給定所有鏈的有效樣本量，基本上占了鏈中的自相關，即當我們從一次抽樣到下一次抽樣時估計的相關性。它實際上不需要很大，但如果它相對於所需的總抽樣數來說很小，那就可能引起關注了。Rhat是衡量鏈的混合程度的指標，當鏈被允許運行無限次抽樣時，它就會變成1。在這種情況下，n_eff和Rhat表明我們有很好的收斂性，但我們也可以用跟蹤圖來直觀地檢查。

# 可視化

srace(fit'))

我們可以看到跟蹤圖中，每條鏈的估計值都能很快地從起點找到一個或多或少的穩定狀態（灰色的初始預燒疊代）。此外，所有三條鏈（每條鏈都用不同的顏色表示）都混合得很好，並在同一結論附近跳動。

stan開發團隊通過shinystan包使交互式探索診斷變得很容易。此外，coda包中還有其他診斷方法，Stan模型的結果可以很容易地轉換為與之配合。下面的代碼演示了如何開始。

bets = extract$beta

除了製作數據列表和產生特定語言的模型代碼的初始設置之外，相對於標準模型，運行貝葉斯回歸模型並不一定需要太多的時間。最主要的也許是採用不同的思維方式，並從這個新的角度來解釋結果。對於你所熟悉的標準模型，可能不需要太長的時間就能像使用那些模型一樣自如。

本文摘選 《 R語言MCMC的rstan貝葉斯回歸模型和標準線性回歸模型比較 》 ，點擊「閱讀原文」獲取全文完整資料。

點擊標題查閱往期內容

R語言RSTAN MCMC：NUTS採樣算法用LASSO 構建貝葉斯線性回歸模型分析職業聲望數據

R語言STAN貝葉斯線性回歸模型分析氣候變化影響北半球海冰範圍和可視化檢查模型收斂性

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採樣用於回歸的貝葉斯估計