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

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

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

可以通過(guò)R使用rstan?包來(lái)調(diào)用Stan，也可以?通過(guò)Python使用?pystan?包。這兩個(gè)接口都支持基于采樣和基于優(yōu)化的推斷，并帶有診斷和后驗(yàn)分析。

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

什么是Stan？

• Stan是命令式概率編程語(yǔ)言。

• Stan程序定義了概率模型。

• 它聲明數(shù)據(jù)和（受約束的）參數(shù)變量。

• 它定義了對(duì)數(shù)后驗(yàn)。

• Stan推理：使模型擬合數(shù)據(jù)并做出預(yù)測(cè)。

• 它可以使用馬爾可夫鏈蒙特卡羅（MCMC）進(jìn)行完整的貝葉斯推斷。

• 使用變分貝葉斯（VB）進(jìn)行近似貝葉斯推斷。

• 最大似然估計(jì)（MLE）用于懲罰最大似然估計(jì)。

Stan計(jì)算什么？

• 得出后驗(yàn)分布?。

• MCMC采樣。

• 繪制

• ，其中每個(gè)繪制

• 都按后驗(yàn)概率

• 的邊緣分布。

• 使用直方圖，核密度估計(jì)等進(jìn)行繪圖

安裝?rstan

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

最后，安裝?rstan：

install.packages(rstan)

Stan中的基本語(yǔ)法

定義模型

Stan模型由六個(gè)程序塊定義?：

• 數(shù)據(jù)（必填）。

• 轉(zhuǎn)換后的數(shù)據(jù)。

• 參數(shù)（必填）。

• 轉(zhuǎn)換后的參數(shù)。

• 模型（必填）。

• 生成的數(shù)量。

數(shù)據(jù)塊讀出的外部信息。

1. data {

2. int N;

3. int x[N];

4. int offset;

5. }

變換后的數(shù)據(jù)?塊允許數(shù)據(jù)的預(yù)處理。

1. transformed data {

2. int y[N];

3. for (n in 1:N)

4. y[n] = x[n] - offset;

5. }

?參數(shù)?塊定義了采樣的空間。

1. parameters {

2. real<lower=0> lambda1;

3. real<lower=0> lambda2;

4. }

變換參數(shù)?塊定義計(jì)算后驗(yàn)之前的參數(shù)處理。

1. transformed parameters {

2. real<lower=0> lambda;

3. lambda = lambda1 + lambda2;

4. }

在?模型?塊中，我們定義后驗(yàn)分布。

1. model {

2. y ~ poisson(lambda);

3. lambda1 ~ cauchy(0, 2.5);

4. lambda2 ~ cauchy(0, 2.5);

5. }

最后，?生成的數(shù)量?塊允許進(jìn)行后處理。

1. generated quantities {

2. int x_predict;

3. x_predict = poisson_rng(lambda) + offset;

4. }

類(lèi)型

Stan有兩種原始數(shù)據(jù)類(lèi)型，?并且兩者都是有界的。

• int?是整數(shù)類(lèi)型。

• real?是浮點(diǎn)類(lèi)型。

1. int<lower=1> N;

2. 


3. real<upper=5> alpha;

4. real<lower=-1,upper=1> beta;

5. 


6. real gamma;

7. real<upper=gamma> zeta;

實(shí)數(shù)擴(kuò)展到線性代數(shù)類(lèi)型。

1. vector[10] a; ? ? // 列向量

2. matrix[10, 1] b;

3. 


4. row_vector[10] c; // 行向量

5. matrix[1, 10] d;

整數(shù)，實(shí)數(shù)，向量和矩陣的數(shù)組均可用。

1. real a[10];

2. 


3. vector[10] b;

4. 


5. matrix[10, 10] c;

Stan還實(shí)現(xiàn)了各種?約束?類(lèi)型。

1. simplex[5] theta; ? ? ? ?// sum(theta) = 1

2. 


3. ordered[5] o; ? ? ? ? ? ?// o[1] < ... < o[5]

4. positive_ordered[5] p;

5. 


6. corr_matrix[5] C; ? ? ? ?// 對(duì)稱(chēng)和

7. cov_matrix[5] Sigma; ? ? // 正定的

關(guān)于Stan的更多信息

所有典型的判斷和循環(huán)語(yǔ)句也都可用。

1. if/then/else

2. 


3. for (i in 1:I)

4. 


5. while (i < I)

有兩種修改?后驗(yàn)的方法。

1. y ~ normal(0, 1);

2. 


3. target += normal_lpdf(y | 0, 1);

4. 


5. # 新版本的Stan中已棄用：

6. increment_log_posterior(log_normal(y, 0, 1))

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

1. parameters {

2. real mu[N];

3. real<lower=0> sigma[N];

4. }

5. 


6. model {

7. // for (n in 1:N)

8. // y[n] ~ normal(mu[n], sigma[n]);

9. 


10. y ~ normal(mu, sigma); ?// 向量化版本

11. }

貝葉斯方法

概率是?認(rèn)知的。例如，?約翰·斯圖亞特·米爾?（John Stuart Mill）說(shuō)：

事件的概率不是事件本身，而是我們或其他人期望發(fā)生的情況的程度。每個(gè)事件本身都是確定的，不是可能的；如果我們?nèi)苛私?，我們?yīng)該或者肯定地知道它會(huì)發(fā)生，或者它不會(huì)。

對(duì)我們來(lái)說(shuō)，概率表示對(duì)它發(fā)生的期望程度。

概率可以量化不確定性。

Stan的貝葉斯示例：重復(fù)試驗(yàn)?zāi)Ｐ?/h1>

我們解決一個(gè)小例子，其中的目標(biāo)是給定從伯努利分布中抽取的隨機(jī)樣本，以估計(jì)缺失參數(shù)的后驗(yàn)分布?

?（成功的機(jī)會(huì)）。

步驟1：?jiǎn)栴}定義

在此示例中，我們將考慮以下結(jié)構(gòu)：

• 數(shù)據(jù)：

• ，試用次數(shù)。

• ，即試驗(yàn)n的結(jié)果??（已知的建模數(shù)據(jù)）。

• 參數(shù)：

• 先驗(yàn)分布

• 概率

• 后驗(yàn)分布

步驟2：Stan

我們創(chuàng)建Stan程序，我們將從R中調(diào)用它。

1. 


2. data {

3. int<lower=0> N; ? ? ? ? ? ? ? // 試驗(yàn)次數(shù)

4. int<lower=0, upper=1> y[N]; ? // 試驗(yàn)成功

5. }

6. 


7. 


8. model {

9. theta ~ uniform(0, 1); ? ? ? ?// 先驗(yàn)

10. y ~ bernoulli(theta); ? ? ? ? // 似然

11. }

步驟3：數(shù)據(jù)

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

1. # 生成數(shù)據(jù)

2. 


3. y = rbinom(N, 1, 0.3)

4. y

## ?[1] 0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1

?根據(jù)數(shù)據(jù)計(jì)算?MLE作為樣本均值：

## [1] 0.25

步驟4：rstan使用貝葉斯后驗(yàn)估計(jì)?

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

1. ##

2. ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 1).

3. ## Chain 1:

4. ## Chain 1: Gradient evaluation took 7e-06 seconds

5. ## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.

6. ## Chain 1: Adjust your expectations accordingly!

7. ## Chain 1:

8. ## Chain 1:

9. ## Chain 1: Iteration: ? ?1 / 5000 [ ?0%] ?(Warmup)

10. ## Chain 1: Iteration: ?500 / 5000 [ 10%] ?(Warmup)

11. ## Chain 1: Iteration: 1000 / 5000 [ 20%] ?(Warmup)

12. ## Chain 1: Iteration: 1500 / 5000 [ 30%] ?(Warmup)

13. ## Chain 1: Iteration: 2000 / 5000 [ 40%] ?(Warmup)

14. ## Chain 1: Iteration: 2500 / 5000 [ 50%] ?(Warmup)

15. ## Chain 1: Iteration: 2501 / 5000 [ 50%] ?(Sampling)

16. ## Chain 1: Iteration: 3000 / 5000 [ 60%] ?(Sampling)

17. ## Chain 1: Iteration: 3500 / 5000 [ 70%] ?(Sampling)

18. ## Chain 1: Iteration: 4000 / 5000 [ 80%] ?(Sampling)

19. ## Chain 1: Iteration: 4500 / 5000 [ 90%] ?(Sampling)

20. ## Chain 1: Iteration: 5000 / 5000 [100%] ?(Sampling)

21. ## Chain 1:

22. ## Chain 1: ?Elapsed Time: 0.012914 seconds (Warm-up)

23. ## Chain 1: ? ? ? ? ? ? ? ?0.013376 seconds (Sampling)

24. ## Chain 1: ? ? ? ? ? ? ? ?0.02629 seconds (Total)

25. ## Chain 1:

26. ...

28. ## SAMPLING FOR MODEL '6dcfbccbf2f063595ccc9b93f383e221' NOW (CHAIN 4).

29. ## Chain 4:

30. ## Chain 4: Gradient evaluation took 3e-06 seconds

31. ## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.03 seconds.

32. ## Chain 4: Adjust your expectations accordingly!

33. ## Chain 4:

34. ## Chain 4:

35. ## Chain 4: Iteration: ? ?1 / 5000 [ ?0%] ?(Warmup)

36. ## Chain 4: Iteration: ?500 / 5000 [ 10%] ?(Warmup)

37. ## Chain 4: Iteration: 1000 / 5000 [ 20%] ?(Warmup)

38. ## Chain 4: Iteration: 1500 / 5000 [ 30%] ?(Warmup)

39. ## Chain 4: Iteration: 2000 / 5000 [ 40%] ?(Warmup)

40. ## Chain 4: Iteration: 2500 / 5000 [ 50%] ?(Warmup)

41. ## Chain 4: Iteration: 2501 / 5000 [ 50%] ?(Sampling)

42. ## Chain 4: Iteration: 3000 / 5000 [ 60%] ?(Sampling)

43. ## Chain 4: Iteration: 3500 / 5000 [ 70%] ?(Sampling)

44. ## Chain 4: Iteration: 4000 / 5000 [ 80%] ?(Sampling)

45. ## Chain 4: Iteration: 4500 / 5000 [ 90%] ?(Sampling)

46. ## Chain 4: Iteration: 5000 / 5000 [100%] ?(Sampling)

47. ## Chain 4:

48. ## Chain 4: ?Elapsed Time: 0.012823 seconds (Warm-up)

49. ## Chain 4: ? ? ? ? ? ? ? ?0.014169 seconds (Sampling)

50. ## Chain 4: ? ? ? ? ? ? ? ?0.026992 seconds (Total)

51. ## Chain 4:

1. ## Inference for Stan model: 6dcfbccbf2f063595ccc9b93f383e221.

2. ## 4 chains, each with iter=5000; warmup=2500; thin=1;

3. ## post-warmup draws per chain=2500, total post-warmup draws=10000.

4. ##

5. ## ? ? ? ? mean se_mean ? sd ? ?10% ? ?90% n_eff Rhat

6. ## theta ? 0.27 ? ?0.00 0.09 ? 0.16 ? 0.39 ?3821 ? ?1

7. ## lp__ ?-13.40 ? ?0.01 0.73 -14.25 -12.90 ?3998 ? ?1

8. ##

1. # 提取后驗(yàn)抽樣

2. # 計(jì)算后均值（估計(jì)）

3. mean(theta_draws)

## [1] 0.2715866# 計(jì)算后驗(yàn)區(qū)間

1. ## ? ? ? 10% ? ? ? 90%

2. ## 0.1569165 0.3934832

1. ggplot(theta_draws_df, aes(x=theta)) +

2. geom_histogram(bins=20, color="gray")

RStan：MAP，MLE

Stan的估算優(yōu)化；兩種觀點(diǎn)：

• 最大后驗(yàn)估計(jì)(MAP)。

• 最大似然估計(jì)（MLE）。

optimizing(model, data=c("N", "y"))

1. ## $par

2. ## theta

3. ## ? 0.4

4. ##

5. ## $value

6. ## [1] -3.4

7. ##

8. ## $return_code

9. ## [1] 0

種群競(jìng)爭(zhēng)模型 ---Lotka-Volterra模型

• 洛特卡（Lotka，1925）和沃爾泰拉（Volterra，1926）制定了參數(shù)化微分方程，描述了食肉動(dòng)物和獵物的競(jìng)爭(zhēng)種群。

• 完整的貝葉斯推斷可用于估計(jì)未來(lái)（或過(guò)去）的種群數(shù)量。

• Stan用于對(duì)統(tǒng)計(jì)模型進(jìn)行編碼并執(zhí)行完整的貝葉斯推理，以解決從噪聲數(shù)據(jù)中推斷參數(shù)的逆問(wèn)題。

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

數(shù)學(xué)模型

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

這里：

• α：獵物增長(zhǎng)速度。

• β：捕食引起的獵物減少速度。

• γ：自然的捕食者減少速度。

• δ：捕食者從捕食中增長(zhǎng)速度。

stan中的Lotka-Volterra

1. real[] dz_dt(data real t, ? ? ? // 時(shí)間

2. real[] z, ? ? ? ? ? ? ? ? ? ? // 系統(tǒng)狀態(tài)

3. real[] theta, ? ? ? ? ? ? ? ? // 參數(shù)

4. data real[] x_r, ? ? ? ? ? ? ?// 數(shù)值數(shù)據(jù)

5. data int[] x_i) ? ? ? ? ? ? ? // 整數(shù)數(shù)據(jù)

6. {

7. real u = z[1]; ? ? ? ? ? ? ? ?// 提取狀態(tài)

8. real v = z[2];

觀察到已知變量：

• ：表示在時(shí)間?

• 的

• 物種數(shù)量

必須推斷未知變量）：

• 初始狀態(tài)：?

• ：k的初始物種數(shù)量。

• 后續(xù)狀態(tài)

• ：在時(shí)間t的物種數(shù)量k。

• 參量?

• 。

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

等效：

與

建立模型

已知常數(shù)和觀測(cè)數(shù)據(jù)的變量。

1. data {

2. int<lower = 0> N; ? ? ? // 數(shù)量測(cè)量

3. real ts[N]; ? ? ? ? ? ? // 測(cè)量次數(shù)>0

4. real y0[2]; ? ? ? ? ? ? // 初始數(shù)量

5. real<lower=0> y[N,2]; ? // 后續(xù)數(shù)量

6. }

未知參數(shù)的變量。

1. parameters {

2. real<lower=0> theta[4]; ? ?// alpha, beta, gamma, delta

3. real<lower=0> z0[2]; ? ? ? // 原始種群

4. real<lower=0> sigma[2]; ? ?// 預(yù)測(cè)誤差

5. }

先驗(yàn)分布和概率。

1. model {

2. // 先驗(yàn)

3. sigma ~ lognormal(0, 0.5);

4. theta[{1, 3}] ~ normal(1, 0.5);

5. 


6. // 似然（對(duì)數(shù)正態(tài)）

7. for (k in 1:2) {

8. y0[k] ~ lognormal(log(z0[k]), sigma[k]);

我們必須為預(yù)測(cè)的總體定義變量?：

• 初始種群（z0）。

• 初始時(shí)間（0.0），時(shí)間（ts）。

• 參數(shù)（theta）。

• 最大迭代次數(shù)（1e3）。

Lotka-Volterra參數(shù)估計(jì)

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

獲得結(jié)果：

1. mean ?se_mean ? sd ?10% ?50% ?90% ?n_eff ?Rhat

2. ## theta[1] ? ?0.55 ? ?0 ? ? 0.07 0.46 0.54 0.64 ? 1168 ? ? 1

3. ## theta[2] ? ?0.03 ? ?0 ? ? 0.00 0.02 0.03 0.03 ? 1305 ? ? 1

4. ## theta[3] ? ?0.80 ? ?0 ? ? 0.10 0.68 0.80 0.94 ? 1117 ? ? 1

5. ## theta[4] ? ?0.02 ? ?0 ? ? 0.00 0.02 0.02 0.03 ? 1230 ? ? 1

6. ## sigma[1] ? ?0.29 ? ?0 ? ? 0.05 0.23 0.28 0.36 ? 2673 ? ? 1

7. ## sigma[2] ? ?0.29 ? ?0 ? ? 0.06 0.23 0.29 0.37 ? 2821 ? ? 1

分析所得結(jié)果：

• Rhat接近1表示收斂；n_eff是有效樣本大小。

• 10％，后驗(yàn)分位數(shù)；例如

• 。

• 后驗(yàn)均值是貝葉斯點(diǎn)估計(jì)：α=0.55。

• 后驗(yàn)平均估計(jì)的標(biāo)準(zhǔn)誤為0。

• α的后驗(yàn)標(biāo)準(zhǔn)偏差為0.07。

最受歡迎的見(jiàn)解

1.matlab使用貝葉斯優(yōu)化的深度學(xué)習(xí)

2.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)

3.R語(yǔ)言Gibbs抽樣的貝葉斯簡(jiǎn)單線性回歸仿真

4.R語(yǔ)言中的block Gibbs吉布斯采樣貝葉斯多元線性回歸

5.R語(yǔ)言中的Stan概率編程MCMC采樣的貝葉斯模型

6.Python用PyMC3實(shí)現(xiàn)貝葉斯線性回歸模型

7.R語(yǔ)言使用貝葉斯 層次模型進(jìn)行空間數(shù)據(jù)分析

8.R語(yǔ)言隨機(jī)搜索變量選擇SSVS估計(jì)貝葉斯向量自回歸（BVAR）模型

9.matlab貝葉斯隱馬爾可夫hmm模型實(shí)現(xiàn)