关于数据挖掘:拓端tecdat视频R语言中的Stan概率编程MCMC采样的贝叶斯模型

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

原文出处：拓端数据部落公众号

概率编程使咱们可能实现统计模型，而不用放心技术细节。这对于基于 MCMC 采样的贝叶斯模型特地有用。

R 语言中 RStan 贝叶斯层次模型剖析示例

stan 简介

Stan 是用于贝叶斯推理的 C ++ 库。它基于 No-U-Turn 采样器（NUTS），该采样器用于依据用户指定的模型和数据预计后验散布。应用 Stan 执行剖析波及以下步骤：

2. 应用 Stan 建模语言指定统计模型。通过专用的_.stan_  文件实现此操作。

4. 筹备要提供给模型的数据。

6. 应用该 stan 函数从后验散布中采样。

8. 剖析后果。

在本文中，我将通过两个层次模型展现 Stan 的用法。我将应用第一个模型探讨 Stan 的基本功能，并应用第二个示例演示更高级的利用。

 学校数据集

咱们要应用的第一个数据集是  学校的数据集。该数据集掂量了教练打算对大学入学考试（在美国应用的学业能力测验（SAT））的影响。数据集如下所示：

正如咱们所看到的：对于八所学校中的大多数，短期教练打算确实进步了 SAT 分数。对于此数据集，咱们有趣味估算与每所学校相干的实在教练打算成果大小。咱们思考两种代替办法。首先，咱们能够假如所有学校彼此独立。然而，这将难以解释，因为学校的后验区间因为高标准差而在很大水平上重叠。第二，假如所有学校的实在成果都雷同，则能够汇总所有学校的数据。然而，这也是不合理的，因为该打算有针对学校的不同成果（例如，不同的老师和学生应该有不同的打算）。

因而，须要另一个模型。分层模型的长处是能够合并来自所有八所学校的信息，而无需假设它们具备独特的实在成果。咱们能够通过以下形式指定档次贝叶斯模型：

依据该模型，教练的成果遵循正态分布，其均值是实在成果 θj，其标准偏差为 σj（从数据中得悉）。真正的影响 θj 遵循参数 μ 和 τ 的正态分布。

定义 Stan 模型文件

在指定了要应用的模型之后，咱们当初能够探讨如何在 Stan 中指定此模型。在为上述模型定义 Stan 程序之前，让咱们看一下 Stan 建模语言的构造。

变量

在 Stan 中，能够通过以下形式定义变量：

int<lower=0> n; # 下界是 0

int<upper=5> n; # 下限是 5

int<lower=0,upper=5> n; # n 的范畴是 \[0,5\]

留神，如果先验已知变量，则应指定变量的高低边界。

多维数据能够通过方括号指定：

vector\[n\] numbers; // 长度为 n 的向量

real\[n\] numbers;  // 长度为 n 的浮点数组

matrix\[n,n\] matrix; // n 乘 n 矩阵 

程序 

Stan 中应用以下程序：

• _data_：用于指定以贝叶斯规定为条件的数据

• _转换后的数据_：用于预处理数据

• 参数 （必填）：用于指定模型的参数

• _转换后的参数_：用于计算后验之前的参数解决

• 模型 （必填）：用于指定模型

• _生成数量_：用于对后果进行后处理

对于  模型  程序块，能够两种等效形式指定散布。第一个，应用以下统计符号：

y ~ normal(mu, sigma); # y 遵从正态分布 

第二种办法应用基于对数概率密度函数（lpdf）的程序化表示法：

target += normal_lpdf(y | mu, sigma); # 减少正态对数密度 

Stan 反对大量的概率分布。通过 Stan 指定模型时，该  lookup 函数会派上用场：它提供从 R 函数到 Stan 函数的映射。思考以下示例：

library(rstan) # 加载 stan 包

lookup(rnorm)

##     StanFunction             Arguments ReturnType Page

## 355   normal_rng (real mu, real sigma)       real  494

在这里，咱们看到 R 中的 rnorm 等价于 Stan 的 normal_rng。

模型

当初，咱们理解了 Stan 建模语言的基础知识，咱们能够定义模型，并将其存储在一个名为的文件中  schools.stan：

留神，θ 永远不会呈现在参数中。这是因为咱们没有显式地对 θ 进行建模，而是对 η（各个学校的标准化成果）进行了建模。而后，依据 μ，τ 和 η 在_变换后的参数_局部结构 θ。此参数化使采样器更高效。

筹备数据进行建模

在拟合模型之前，咱们须要将输出数据编码为一个列表，其参数应与 Stan 模型的数据局部绝对应。对于学校数据，数据如下：

schools.data <- list(

  n = 8,

  y = c(28,  8, -3,  7, -1,  1, 18, 12),

  sigma = c(15, 10, 16, 11,  9, 11, 10, 18)

)

从后验散布抽样

咱们能够应用 stan 函数从后验散布中采样，函数执行以下三个步骤：

2. 它将模型标准转换为 C ++ 代码。

4. 它将 C ++ 代码编译为共享对象。

6. 它依据指定的模型，数据和设置从后验散布中采样。

如果  rstan_options(auto_write = TRUE)，则雷同模型的后续调用将比第一次调用快得多，因为该  stan 函数随后跳过了前两个步骤（转换和编译模型）。此外，咱们将设置要应用的内核数：

options(mc.cores = parallel::detectCores()) # 并行化

rstan\_options(auto\_write = TRUE)  # 存储编译的 stan 模型 

当初，咱们能够从后验中编译模型和样本。

模型解释

咱们将首先对模型进行根本解释，而后钻研 MCMC 程序。

根本模型解释

要应用拟合模型执行推断，咱们能够应用  print 函数。

print(fit1) # 可选参数：pars，probs

## Inference for Stan model: schools.

## 4 chains, each with iter=2000; warmup=1000; thin=1; 

## post-warmup draws per chain=1000, total post-warmup draws=4000.

## 

##            mean se\_mean   sd   2.5%    25%    50%    75%  97.5% n\_eff Rhat

## mu         7.67    0.15 5.14  -2.69   4.42   7.83  10.93  17.87  1185    1

## tau        6.54    0.16 5.40   0.31   2.52   5.28   9.05  20.30  1157    1

## eta\[1\]     0.42    0.01 0.92  -1.47  -0.18   0.44   1.03   2.18  4000    1

## eta\[2\]     0.03    0.01 0.87  -1.74  -0.54   0.03   0.58   1.72  4000    1

## eta\[3\]    -0.18    0.02 0.92  -1.95  -0.81  -0.20   0.45   1.65  3690    1

## eta\[4\]    -0.03    0.01 0.92  -1.85  -0.64  -0.02   0.57   1.81  4000    1

## eta\[5\]    -0.33    0.01 0.86  -2.05  -0.89  -0.34   0.22   1.43  3318    1

## eta\[6\]    -0.20    0.01 0.87  -1.91  -0.80  -0.21   0.36   1.51  4000    1

## eta\[7\]     0.37    0.02 0.87  -1.37  -0.23   0.37   0.96   2.02  3017    1

## eta\[8\]     0.05    0.01 0.92  -1.77  -0.55   0.05   0.69   1.88  4000    1

## theta\[1\]  11.39    0.15 8.09  -2.21   6.14  10.30  15.56  30.22  2759    1

## theta\[2\]   7.92    0.10 6.25  -4.75   4.04   8.03  11.83  20.05  4000    1

## theta\[3\]   6.22    0.14 7.83 -11.41   2.03   6.64  10.80  20.97  3043    1

## theta\[4\]   7.58    0.10 6.54  -5.93   3.54   7.60  11.66  20.90  4000    1

## theta\[5\]   5.14    0.10 6.30  -8.68   1.40   5.63   9.50  16.12  4000    1

## theta\[6\]   6.08    0.10 6.62  -8.06   2.21   6.45  10.35  18.53  4000    1

## theta\[7\]  10.60    0.11 6.70  -0.94   6.15  10.01  14.48  25.75  4000    1

## theta\[8\]   8.19    0.14 8.18  -8.13   3.59   8.01  12.48  25.84  3361    1

## lp__     -39.47    0.07 2.58 -45.21 -41.01 -39.28 -37.70 -34.99  1251    1

## 

## Samples were drawn using NUTS(diag_e) at Thu Nov 29 11:17:50 2018.

## For each parameter, n_eff is a crude measure of effective sample size,

## and Rhat is the potential scale reduction factor on split chains (at 

## convergence, Rhat=1).

在此，行名称示意预计的参数：mu 是后验散布的平均值，而 tau 是其标准偏差。eta 和 theta 的条目别离示意矢量 η 和 θ 的估计值。这些列示意计算值。百分比示意置信区间。例如，教练打算的总体成果的 95%可信区间 μ 为 [-1.27,18.26]。因为咱们不确定平均值，因而 θj 的 95%置信区间也很宽。例如，对于第一所学校，95%置信区间为 [−2.19,32.33]。

咱们能够应用以下 plot 函数来可视化预计中的不确定性：

黑线示意 95%的距离，而红线示意 80%的距离。圆圈示意平均值的预计。

咱们能够应用以下 extract 函数获取生成的样本：

# 获取样本

samples <- extract(fit1, permuted = TRUE) # 每个参数 1000 个样本 

MCMC 诊断

 通过绘制采样过程的轨迹图，咱们能够确定采样期间是否出了问题。例如，链条在一个地位停留的工夫过长或在一个方向上走了太多步，就会有问题。咱们能够应用 traceplot 函数绘制模型中应用的四个链的轨迹：

# 诊断:

要从各个马尔可夫链中获取样本，咱们能够 extract 再次应用函数：

##          parameters

## chains           mu       tau     eta\[1\]     eta\[2\]     eta\[3\]     eta\[4\]

##   chain:1  1.111120  2.729124 -0.1581242 -0.8498898  0.5025965 -1.9874554

##   chain:2  3.633421  2.588945  1.2058772 -1.1173221  1.4830778  0.4838649

##   chain:3 13.793056  3.144159  0.6023924 -1.1188243 -1.2393491 -0.6118482

##   chain:4  3.673380 13.889267 -0.0869434  1.1900236 -0.0378830 -0.2687284

##          parameters

## chains        eta\[5\]     eta\[6\]     eta\[7\]      eta\[8\]   theta\[1\]

##   chain:1  0.3367602 -1.1940843  0.5834020 -0.08371249  0.6795797

##   chain:2 -1.8057252  0.7429594  0.9517675  0.55907356  6.7553706

##   chain:3 -1.5867789  0.6334288 -0.4613463 -1.44533007 15.6870727

##   chain:4  0.1028605  0.3481214  0.9264762  0.45331024  2.4657999

##          parameters

## chains     theta\[2\] theta\[3\]    theta\[4\]  theta\[5\]  theta\[6\]  theta\[7\]

##   chain:1 -1.208335 2.482769 -4.31289292  2.030181 -2.147684  2.703297

##   chain:2  0.740736 7.473028  4.88612054 -1.041502  5.556902  6.097494

##   chain:3 10.275294 9.896345 11.86930758  8.803971 15.784656 12.342510

##   chain:4 20.201935 3.147213 -0.05906019  5.102037  8.508530 16.541455

##          parameters

## chains     theta\[8\]      lp__

##   chain:1 0.8826584 -41.21499

##   chain:2 5.0808317 -41.17178

##   chain:3 9.2487083 -40.35351

##   chain:4 9.9695268 -36.34043

为了对采样过程进行更高级的剖析，咱们能够应用该  shinystan 软件包。应用该软件包，能够通过以下形式启动 Shiny 应用程序来剖析拟合模型：

library(shinystan)

launch_shinystan(fit1)

档次回归

当初，咱们对 Stan 有了根本的理解，咱们能够深入研究更高级的应用程序：让咱们尝试一下档次回归。在惯例回归中，咱们对以下模式的关系进行建模

此示意假如所有样本都具备雷同的散布。如果只存在一组样本，那么咱们就会遇到问题，因为将疏忽组内和组之间的潜在差别。

另一种抉择是为每个组建设一个回归模型。然而，在这种状况下，预计单个模型时，小样本量会带来问题。

档次回归是两个极其之间的折衷。该模型假如组是类似的，但存在差别。

假如每个样本都属于 K 组之一。而后，档次回归指定如下：

其中 Yk 是第 k 组的后果，αk 是截距，Xk 是特色，β（k）示意权重。层次模型不同于其中 Yk 别离拟合每个组的模型，因为假设参数 αk 和 β（k）源自独特的散布。

 数据集

分层回归的经典示例是 老鼠数据集。该数据集蕴含 5 周内测得的 鼠体重。让咱们加载数据：

##   day8 day15 day22 day29 day36

## 1  151   199   246   283   320

## 2  145   199   249   293   354

## 3  147   214   263   312   328

## 4  155   200   237   272   297

## 5  135   188   230   280   323

## 6  159   210   252   298   331

让咱们考察数据：

library(ggplot2)

ggplot(ddf, aes(x = variable, y = value, group = Group)) + geom\_line() + geom\_point()

 数据显示线性增长趋势对于不同的大鼠十分类似。然而，咱们还看到，大鼠的初始体重不同，须要不同的截距，并且成长速度也须要不同的斜率。因而，分层模型仿佛是适当的。

档次回归模型的标准

该模型能够指定如下：

第 i 个大鼠的截距由 αi 示意，斜率由 βi 示意。留神，测量工夫的核心是 x = 22，它是工夫序列数据的中值测量值（第 22 天）。

当初，咱们能够指定模型并将其存储在名为 rats.stan 的文件中：

请留神，模型代码估算的是方差（sigmasq  变量）而不是标准差。

材料筹备

为了筹备模型数据，咱们首先将测量点提取为数值，而后将所有内容编码为列表构造：

data <- list(N = nrow(df), T = ncol(df), x = days,

                 y = df, xbar = median(days))

拟合回归模型

当初，咱们能够为老鼠体重数据集拟合贝叶斯档次回归模型：

# 模型蕴含截距（alpha）和斜率（beta）的预计 

档次回归模型的预测

在确定了每只大鼠的 α 和 β 之后，咱们当初能够预计任意工夫点单个大鼠的体重。在这里，咱们寻找从第 0 天到第 100 天的大鼠体重。

ggplot(pred.df\[pred.df$Rat %in% sel.rats, \], 

       aes(x = Day, y = Weight, group = Rat, 



    geom_line()  +

与原始数据相比，该模型的预计是平滑的，因为每条曲线都遵循线性模型。钻研最初一个图中所示的置信区间，咱们能够看到方差预计是正当的。咱们对采样时（第 8 至 36 天）的老鼠体重充满信心，然而随着来到采样区域，不确定性会减少。

最受欢迎的见解

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 模型实现