拓端数据tecdat|R语言RStan贝叶斯示例：重复试验模型和种群竞争模型Lotka Volterra

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

 

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模型由六个程序块定义 ：

• 数据（必填）。
• 转换后的数据。
• 参数（必填）。
• 转换后的参数。
• 模型（必填）。
• 生成的数量。

数据块读出的外部信息。

1.  
   data {
2.  
   int N;
3.  
   int x[N];
4.  
   int offset;
5.  
   }

变换后的数据 块允许数据的预处理。

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

 参数 块定义了采样的空间。

1.  
   parameters {
2.  
   real<lower=0> lambda1;
3.  
   real<lower=0> lambda2;
4.  
   }

变换参数 块定义计算后验之前的参数处理。

1.  
   transformed parameters {
2.  
   real<lower=0> lambda;
3.  
   lambda = lambda1 + lambda2;
4.  
   }

在 模型 块中，我们定义后验分布。

1.  
   model {
2.  
   y ~ poisson(lambda);
3.  
   lambda1 ~ cauchy(0, 2.5);
4.  
   lambda2 ~ cauchy(0, 2.5);
5.  
   }

最后， 生成的数量 块允许进行后处理。

1.  
   generated quantities {
2.  
   int x_predict;
3.  
   x_predict = poisson_rng(lambda) + offset;
4.  
   }

类型

Stan有两种原始数据类型， 并且两者都是有界的。

• int 是整数类型。
• real 是浮点类型。

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;

实数扩展到线性代数类型。

1.  
   vector[10] a; // 列向量
2.  
   matrix[10, 1] b;
3.  
    
4.  
   row_vector[10] c; // 行向量
5.  
   matrix[1, 10] d;

整数，实数，向量和矩阵的数组均可用。

1.  
   real a[10];
2.  
    
3.  
   vector[10] b;
4.  
    
5.  
   matrix[10, 10] c;

Stan还实现了各种 约束 类型。

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; // 对称和
7.  
   cov_matrix[5] Sigma; // 正定的

关于Stan的更多信息

所有典型的判断和循环语句也都可用。

1.  
   if/then/else
2.  
    
3.  
   for (i in 1:I)
4.  
    
5.  
   while (i < I)

有两种修改 后验的方法。

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))

而且许多采样语句都是 矢量化的。

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.  
    }

贝叶斯方法

概率是 认知的。例如， 约翰·斯图亚特·米尔 （John Stuart Mill）说：

事件的概率不是事件本身，而是我们或其他人期望发生的情况的程度。每个事件本身都是确定的，不是可能的；如果我们全部了解，我们应该或者肯定地知道它会发生，或者它不会。

对我们来说，概率表示对它发生的期望程度。

概率可以量化不确定性。

Stan的贝叶斯示例：重复试验模型

我们解决一个小例子，其中的目标是给定从伯努利分布中抽取的随机样本，以估计缺失参数的后验分布  （成功的机会）。

步骤1：问题定义

在此示例中，我们将考虑以下结构：

• 数据：

  • ，试用次数。
  • ，即试验n的结果  （已知的建模数据）。

• 参数：

• 先验分布

• 概率

• 后验分布

步骤2：Stan

我们创建Stan程序，我们将从R中调用它。

1.  
    
2.  
   data {
3.  
   int<lower=0> N; // 试验次数
4.  
   int<lower=0, upper=1> y[N]; // 试验成功
5.  
   }
6.  
    
7.  
    
8.  
   model {
9.  
   theta ~ uniform(0, 1); // 先验
10.  
    y ~ bernoulli(theta); // 似然
11.  
    }

步骤3：数据

在这种情况下，我们将使用示例随机模拟一个随机样本，而不是使用给定的数据集。

1.  
   # 生成数据
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

 根据数据计算 MLE作为样本均值：

## [1] 0.25

步骤4：rstan使用贝叶斯后验估计 

最后一步是使用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.  
    ...
27.  
     
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.  
   # 提取后验抽样
2.  
   # 计算后均值（估计）
3.  
   mean(theta_draws)

## [1] 0.2715866

# 计算后验区间

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的估算优化；两种观点：

• 最大后验估计(MAP)。
• 最大似然估计（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

种群竞争模型 ---Lotka-Volterra模型

• 洛特卡（Lotka，1925）和沃尔泰拉（Volterra，1926）制定了参数化微分方程，描述了食肉动物和猎物的竞争种群。
• 完整的贝叶斯推断可用于估计未来（或过去）的种群数量。
• Stan用于对统计模型进行编码并执行完整的贝叶斯推理，以解决从噪声数据中推断参数的逆问题。

 

在此示例中，我们希望根据公司每年收集的毛皮数量，将模型拟合到1900年至1920年之间各自种群的加拿大猫科食肉动物和野兔猎物。

数学模型

我们表示U（t）和V（t）作为猎物和捕食者种群数量 分别。与它们相关的微分方程为：

这里：

• α：猎物增长速度。
• β：捕食引起的猎物减少速度。
• γ：自然的捕食者减少速度。
• δ：捕食者从捕食中增长速度。

stan中的Lotka-Volterra

1.  
   real[] dz_dt(data real t, // 时间
2.  
   real[] z, // 系统状态
3.  
   real[] theta, // 参数
4.  
   data real[] x_r, // 数值数据
5.  
   data int[] x_i) // 整数数据
6.  
   {
7.  
   real u = z[1]; // 提取状态
8.  
   real v = z[2];

观察到已知变量：

• ：表示在时间 的物种数量

必须推断未知变量）：

• 初始状态： ：k的初始物种数量。
• 后续状态：在时间t的物种数量k。
• 参量 。

假设误差是成比例的（而不是相加的）：

等效：

与

建立模型

已知常数和观测数据的变量。

1.  
   data {
2.  
   int<lower = 0> N; // 数量测量
3.  
   real ts[N]; // 测量次数>0
4.  
   real y0[2]; // 初始数量
5.  
   real<lower=0> y[N,2]; // 后续数量
6.  
   }

未知参数的变量。

1.  
   parameters {
2.  
   real<lower=0> theta[4]; // alpha, beta, gamma, delta
3.  
   real<lower=0> z0[2]; // 原始种群
4.  
   real<lower=0> sigma[2]; // 预测误差
5.  
   }

先验分布和概率。

1.  
   model {
2.  
   // 先验
3.  
   sigma ~ lognormal(0, 0.5);
4.  
   theta[{1, 3}] ~ normal(1, 0.5);
5.  
    
6.  
   // 似然（对数正态）
7.  
   for (k in 1:2) {
8.  
   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))

获得结果：

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

分析所得结果：

• Rhat接近1表示收敛；n_eff是有效样本大小。
• 10％，后验分位数；例如。
• 后验均值是贝叶斯点估计：α=0.55。
• 后验平均估计的标准误为0。
• α的后验标准偏差为0.07。

---

最受欢迎的见解

1.matlab使用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简单线性回归仿真

4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言使用贝叶斯 层次模型进行空间数据分析

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.matlab贝叶斯隐马尔可夫hmm模型实现