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

Stan是一種用于指定統計模型的概率編程語言。Stan通過馬爾可夫鏈蒙特卡羅方法（例如No-U-Turn采樣器，一種漢密爾頓蒙特卡洛采樣的自適應形式）為連續變量模型提供了完整的貝葉斯推斷。

可以通過R使用rstan?包來調用Stan，也可以?通過Python使用?pystan?包。這兩個接口都支持基于采樣和基于優化的推斷，并帶有診斷和后驗分析。

在本文中，簡要展示了Stan的主要特性。還顯示了兩個示例：第一個示例與簡單的伯努利模型相關，第二個示例與基于常微分方程的Lotka-Volterra模型有關。

什么是Stan？

• Stan是命令式概率編程語言。
• Stan程序定義了概率模型。
• 它聲明數據和（受約束的）參數變量。
• 它定義了對數后驗。
• Stan推理：使模型擬合數據并做出預測。
• 它可以使用馬爾可夫鏈蒙特卡羅（MCMC）進行完整的貝葉斯推斷。
• 使用變分貝葉斯（VB）進行近似貝葉斯推斷。
• 最大似然估計（MLE）用于懲罰最大似然估計。

Stan計算什么？

• 得出后驗分布?。
• MCMC采樣。
• 繪制，其中每個繪制都按后驗概率的邊緣分布。
• 使用直方圖，核密度估計等進行繪圖

安裝?rstan

要在R中運行Stan，必須安裝?rstan?C ++編譯器。在Windows上，?Rtools?是必需的。

最后，安裝?rstan：

install.packages(rstan)

Stan中的基本語法

定義模型

Stan模型由六個程序塊定義?：

• 數據（必填）。
• 轉換后的數據。
• 參數（必填）。
• 轉換后的參數。
• 模型（必填）。
• 生成的數量。

數據塊讀出的外部信息。

data {int N;int x[N];int offset; }

變換后的數據?塊允許數據的預處理。

transformed data {int y[N];for (n in 1:N)y[n] = x[n] - offset; }

?參數?塊定義了采樣的空間。

parameters {real<lower=0> lambda1;real<lower=0> lambda2; }

變換參數?塊定義計算后驗之前的參數處理。

transformed parameters {real<lower=0> lambda;lambda = lambda1 + lambda2; }

在?模型?塊中，我們定義后驗分布。

model {y ~ poisson(lambda);lambda1 ~ cauchy(0, 2.5);lambda2 ~ cauchy(0, 2.5); }

最后，?生成的數量?塊允許進行后處理。

generated quantities {int x_predict;x_predict = poisson_rng(lambda) + offset; }

類型

Stan有兩種原始數據類型，?并且兩者都是有界的。

• int?是整數類型。
• real?是浮點類型。

int<lower=1> N;real<upper=5> alpha; real<lower=-1,upper=1> beta;real gamma; real<upper=gamma> zeta;

實數擴展到線性代數類型。

vector[10] a; // 列向量 matrix[10, 1] b;row_vector[10] c; // 行向量 matrix[1, 10] d;

整數，實數，向量和矩陣的數組均可用。

real a[10];vector[10] b;matrix[10, 10] c;

Stan還實現了各種?約束?類型。

simplex[5] theta; // sum(theta) = 1ordered[5] o; // o[1] < ... < o[5] positive_ordered[5] p;corr_matrix[5] C; // 對稱和 cov_matrix[5] Sigma; // 正定的

關于Stan的更多信息

所有典型的判斷和循環語句也都可用。

if/then/elsefor (i in 1:I)while (i < I)

有兩種修改?后驗的方法。

y ~ normal(0, 1);target += normal_lpdf(y | 0, 1);# 新版本的Stan中已棄用： increment_log_posterior(log_normal(y, 0, 1))

而且許多采樣語句都是?矢量化的。

parameters {real mu[N];real<lower=0> sigma[N]; }model {// for (n in 1:N)// y[n] ~ normal(mu[n], sigma[n]);y ~ normal(mu, sigma); // 向量化版本 }

貝葉斯方法

概率是?認知的。例如，?約翰·斯圖亞特·米爾?（John Stuart Mill）說：

事件的概率不是事件本身，而是我們或其他人期望發生的情況的程度。每個事件本身都是確定的，不是可能的；如果我們全部了解，我們應該或者肯定地知道它會發生，或者它不會。

對我們來說，概率表示對它發生的期望程度。

概率可以量化不確定性。

Stan的貝葉斯示例：重復試驗模型

我們解決一個小例子，其中的目標是給定從伯努利分布中抽取的隨機樣本，以估計缺失參數的后驗分布??（成功的機會）。

步驟1：問題定義

在此示例中，我們將考慮以下結構：

• 數據：

  • ，試用次數。
  • ，即試驗n的結果??（已知的建模數據）。

• 參數：

• 先驗分布

• 概率

• 后驗分布

步驟2：Stan

我們創建Stan程序，我們將從R中調用它。

data {int<lower=0> N; // 試驗次數int<lower=0, upper=1> y[N]; // 試驗成功 }model {theta ~ uniform(0, 1); // 先驗y ~ bernoulli(theta); // 似然 }

步驟3：數據

在這種情況下，我們將使用示例隨機模擬一個隨機樣本，而不是使用給定的數據集。

# 生成數據y = rbinom(N, 1, 0.3) y ## [1] 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1

?根據數據計算?MLE作為樣本均值：

## [1] 0.25

步驟4：rstan使用貝葉斯后驗估計?

最后一步是使用R中的Stan獲得我們的估算值。

## ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 1). ## Chain 1: ## Chain 1: Gradient evaluation took 7e-06 seconds ## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds. ## Chain 1: Adjust your expectations accordingly! ## Chain 1: ## Chain 1: ## Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 1: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 1: ## Chain 1: Elapsed Time: 0.012914 seconds (Warm-up) ## Chain 1: 0.013376 seconds (Sampling) ## Chain 1: 0.02629 seconds (Total) ## Chain 1: ...## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 4). ## Chain 4: ## Chain 4: Gradient evaluation took 3e-06 seconds ## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds. ## Chain 4: Adjust your expectations accordingly! ## Chain 4: ## Chain 4: ## Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup) ## Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup) ## Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup) ## Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup) ## Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup) ## Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup) ## Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling) ## Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling) ## Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling) ## Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling) ## Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling) ## Chain 4: Iteration: 5000 / 5000 [100%] (Sampling) ## Chain 4: ## Chain 4: Elapsed Time: 0.012823 seconds (Warm-up) ## Chain 4: 0.014169 seconds (Sampling) ## Chain 4: 0.026992 seconds (Total) ## Chain 4: ## Inference for Stan model: 6dcfbccbf2f063595ccc9b93f383e221. ## 4 chains, each with iter=5000; warmup=2500; thin=1; ## post-warmup draws per chain=2500, total post-warmup draws=10000. ## ## mean se_mean sd 10% 90% n_eff Rhat ## theta 0.27 0.00 0.09 0.16 0.39 3821 1 ## lp__ -13.40 0.01 0.73 -14.25 -12.90 3998 1 ## # 提取后驗抽樣 # 計算后均值（估計） mean(theta_draws) ## [1] 0.2715866 # 計算后驗區間 ## 10% 90% ## 0.1569165 0.3934832 ggplot(theta_draws_df, aes(x=theta)) +geom_histogram(bins=20, color="gray")

RStan：MAP，MLE

Stan的估算優化；兩種觀點：

• 最大后驗估計(MAP)。
• 最大似然估計（MLE）。

optimizing(model, data=c("N", "y")) ## $par ## theta ## 0.4 ## ## $value ## [1] -3.4 ## ## $return_code ## [1] 0

種群競爭模型 ---Lotka-Volterra模型

• 洛特卡（Lotka，1925）和沃爾泰拉（Volterra，1926）制定了參數化微分方程，描述了食肉動物和獵物的競爭種群。
• 完整的貝葉斯推斷可用于估計未來（或過去）的種群數量。
• Stan用于對統計模型進行編碼并執行完整的貝葉斯推理，以解決從噪聲數據中推斷參數的逆問題。

在此示例中，我們希望根據公司每年收集的毛皮數量，將模型擬合到1900年至1920年之間各自種群的加拿大貓科食肉動物和野兔獵物。

數學模型

我們表示U（t）和V（t）作為獵物和捕食者種群數量?分別。與它們相關的微分方程為：

這里：

• α：獵物增長速度。
• β：捕食引起的獵物減少速度。
• γ：自然的捕食者減少速度。
• δ：捕食者從捕食中增長速度。

stan中的Lotka-Volterra

real[] dz_dt(data real t, // 時間real[] z, // 系統狀態real[] theta, // 參數data real[] x_r, // 數值數據data int[] x_i) // 整數數據 {real u = z[1]; // 提取狀態real v = z[2];

觀察到已知變量：

• ：表示在時間?的物種數量

必須推斷未知變量）：

• 初始狀態：?：k的初始物種數量。
• 后續狀態：在時間t的物種數量k。
• 參量?。

假設誤差是成比例的（而不是相加的）：

等效：

與

建立模型

已知常數和觀測數據的變量。

data {int<lower = 0> N; // 數量測量real ts[N]; // 測量次數>0real y0[2]; // 初始數量real<lower=0> y[N,2]; // 后續數量 }

未知參數的變量。

parameters {real<lower=0> theta[4]; // alpha, beta, gamma, deltareal<lower=0> z0[2]; // 原始種群real<lower=0> sigma[2]; // 預測誤差 }

先驗分布和概率。

model {// 先驗sigma ~ lognormal(0, 0.5);theta[{1, 3}] ~ normal(1, 0.5);// 似然（對數正態）for (k in 1:2) {y0[k] ~ lognormal(log(z0[k]), sigma[k]);

我們必須為預測的總體定義變量?：

• 初始種群（z0）。
• 初始時間（0.0），時間（ts）。
• 參數（theta）。
• 最大迭代次數（1e3）。

Lotka-Volterra參數估計

print(fit, c("theta", "sigma"), probs=c(0.1, 0.5, 0.9))

獲得結果：

mean se_mean sd 10% 50% 90% n_eff Rhat ## theta[1] 0.55 0 0.07 0.46 0.54 0.64 1168 1 ## theta[2] 0.03 0 0.00 0.02 0.03 0.03 1305 1 ## theta[3] 0.80 0 0.10 0.68 0.80 0.94 1117 1 ## theta[4] 0.02 0 0.00 0.02 0.02 0.03 1230 1 ## sigma[1] 0.29 0 0.05 0.23 0.28 0.36 2673 1 ## sigma[2] 0.29 0 0.06 0.23 0.29 0.37 2821 1

分析所得結果：

• Rhat接近1表示收斂；n_eff是有效樣本大小。
• 10％，后驗分位數；例如。
• 后驗均值是貝葉斯點估計：α=0.55。
• 后驗平均估計的標準誤為0。
• α的后驗標準偏差為0.07。

---

總結

以上是生活随笔為你收集整理的R语言RStan贝叶斯示例：重复试验模型和种群竞争模型Lotka Volterra的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。