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

由Kaizong Ye，Liao Bao撰写

Stan是用于贝叶斯推理的C ++库。它基于No-U-Turn采样器（NUTS），该采样器用于根据用户指定的模型和数据估计后验分布。

使用Stan执行分析涉及以下步骤：

使用Stan建模语言指定统计模型。通过专用的.stan 文件完成此操作。
准备要提供给模型的数据。
使用该stan 函数从后验分布中采样。
分析结果。

分层贝叶斯模型 (Bayesian hieratical model) 首先是基于可交换性的基础的。不记得的自觉看前文。

在理解了可交换性这个基础之后，一起来看到这里的模型的外在。

其实这个模型可以分两层，也可以分很多很多层。如果仅考虑两层的话，类似于此： $P(\theta,\phi|Y)=\frac{P(Y|\theta,\phi)P(\theta,\phi)}{P(Y)}=\frac{P(Y|\theta)P(\theta|\phi)P(\phi)}{P(Y)}\\$

三阶段的话，则会像是： $P(\theta,\phi,X|Y)=\frac{P(Y|\theta)P(\theta|\phi)P(\phi|X)P(X)}{P(Y)}\\$

但这个式子怎么用来做预测呢？

总归得先看公式的.

如果还没忘记可交换性的概念的话(忘记了的请举手，并自觉往上看)，从一个群体中独立采样的结果是可以看作是可交换的。这样的话群体的属性就可以用一个固定的未知参数 $\phi$ 来描述，即： $\begin{align} \phi &\sim p(\phi)\\ \left\{ Y_1,Y_2,\cdots,Y_n|\phi \right\}&\sim^{i.i.d.} p(y|\phi)\\ \end{align}\\$

从而当考虑分层数据 $\left\{ Y_1,Y_2,\cdots,Y_n\right\}$ ，其中 $Y_j= \left\{ Y_{1,j},Y_{2.j},\cdots,Y_{n_j,j} \right\}$ 时，有 $\left\{ Y_{1,j},\cdots,Y_{n_j,j}\right\}\sim^{i.i.d.}p\left(y|\phi_j\right)\\$

而对于 $\phi_j$ 同样而可以假定： $\left\{\phi_1,\phi_2,\cdots,\phi_n|\phi\right\}\sim^{i.i.d.}p\left(\phi|\psi\right)\\$

这样就有了三个模型

组内： $\left\{ Y_{1,j},\cdots,Y_{n_j,j}\right\}\sim^{i.i.d.}p\left(y|\phi_j\right)$

组间： $\left\{\phi_1,\phi_2,\cdots,\phi_n|\phi\right\}\sim^{i.i.d.}p\left(\phi|\psi\right)$

先验分布： $\psi \sim P(\psi)$

而为了表现各组的差异，回归手段就需要被派上用场了

组内模型： $\phi_j=N(\theta_j,\sigma^2 ),\qquad P(y|\phi_j )=normal(\theta_j,\sigma^2)$

组间模型： $\psi=N(\mu,\tau^2 ),\qquad P(\theta_j |\psi)=normal(\mu,\tau^2 )$

而要估计这些参数，方法还是挺多的，常用的有蒙特卡洛方法 (Monte Carlo Method, MC) 马尔科夫链蒙特卡洛方法 (Markov Chain Monte Carlo Method, MCMC)

蒙特卡洛方法即使用随机数解决实际问题，比如求个积分。

为了计算 $J=\int_a^b{f(x)dx}$ ，取一列服从 $[a,b]$ 上均匀分布的随机变量 $\{X_n\}$ ，则 $f(X_n) i.i.d.$ 且 $E\left[f(X_i )\right]=\frac{1}{b-a} \int_a^b{f(x)dx}=\frac{1}{b-a}J\\$

从而 $J=(b-a)E[f(X_i)]$ ，而由大数定理 $\frac{1}{n} \sum_{k=1}^n\xrightarrow{~~p~~}{f(X_k)}E[f(X_i)]\\$

而当结果为 $a$ 时，对 $\forall \varepsilon>0″ style=”vertical-align: middle; margin-right: 3px; margin-left: 3px; display: inline-block;”> ，有 <img decoding=$

MCMC入门

而MCMC则是基于马尔科夫链的蒙特卡洛方法，(完全不了解或已经忘了马尔科夫链的话戳这里 )

在引入Gibbs采样方法之前，一直使用的是一个 M-H (Metropolis-Hastings) 采样方法，一共分两步走

初始化马尔科夫链初始状态 $X_0=x_0$
对 $t=0,1,2,⋯$ 循环以下过程
第 $t$ 个时刻马氏链状态为 $X_t=x_t$ ，采样 $y\sim q(x|x_t )$
从均匀分布采样 $u\sim Uniform[0,1]$
如果 $u<\alpha(x_t,y)=\min{\left\{\frac{p(y)p(x_t|y)}{p(x_t)p(y|x_t)} ,1\right\} }\\$ 则接受转移 $x_t\rightarrow y$ ，即 $X_{t+1}\rightarrow y$
否则不接受转移，即 $X_{t+1}=x_t$

收敛之后就得到了所需的样本。然而接受率 $\alpha$ 通常小于1，使 M-H 算法效率不够高，为了提高效率有了 Gibbs 方法。

通过构造转移矩阵Q $\begin{align} &Q(A\rightarrow B)=p(y_b|x_1)\qquad &若x_A=x_B=x_1\\ &Q(A\rightarrow C)=p(x_c|y_1) &若y_A=y_C=y_1\\ &Q(A\rightarrow D)=0 &其他 \end{align}\\$

从而该转移矩阵对平面内任意两点均满足细致平稳条件（并可推广到多维）

随机初始化 $\left\{x_i:i=1,2,⋯,n\right\}$
对 $t=1,2,\cdots$ 循环采样
$x_1^{(t+1)}\sim p\left(x_1 |x_2^{(t)},x_3^{(t)},⋯,x_n^{(t)}\right)$
$x_2^{(t+1)}\sim p\left(x_2 |x_1^{(t+1)},x_3^{(t)},⋯,x_n^{(t)}\right)$
$\cdots$
$x_j^{(t+1)}\sim p\left(x_j |x_1^{(t+1)},\cdots,x_{j-1}^{(t+1)},x_{j+1}^{(t)}⋯,x_n^{(t)}\right)$
$\cdots$
$x_ n^{(t+1)}\sim p\left(x_n |x_1^{(t+1)},x_2^{(t+1)},⋯,x_{n-1}^{(t+1)}\right)$

收敛后就可以得到样本了。

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

学校数据集

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

可下载资源

完整程序、数据和文档（word）

最受欢迎的见解

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

正如我们所看到的：对于八所学校中的大多数，短期教练计划的确提高了SAT分数。对于此数据集，我们有兴趣估算与每所学校相关的真实教练计划效果大小。我们考虑两种替代方法。首先，我们可以假设所有学校彼此独立。但是，这将难以解释，因为学校的后验区间由于高标准差而在很大程度上重叠。第二，假设所有学校的真实效果都相同，则可以汇总所有学校的数据。但是，这也是不合理的，因为该计划有针对学校的不同效果（例如，不同的老师和学生应该有不同的计划）。

因此，需要另一个模型。分层模型的优点是可以合并来自所有八所学校的信息，而无需假定它们具有共同的真实效果。我们可以通过以下方式指定层次贝叶斯模型：

视频

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

根据该模型，教练的效果遵循正态分布，其均值是真实效果θ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]

注意，如果先验已知变量，则应指定变量的上下边界。

多维数据可以通过方括号指定：

R语言stan进行贝叶斯推理分析

阅读文章 ➜


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 函数从后验分布中采样，函数执行以下三个步骤：

它将模型规范转换为C ++代码。
它将C ++代码编译为共享对象。
它根据指定的模型，数据和设置从后验分布中采样。

如果 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天）的老鼠体重充满信心，但是随着离开采样区域，不确定性会增加。

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。在此对他对本文所作的贡献表示诚挚感谢，他在上海财经大学完成了统计学专业的硕士学位，专注人工智能领域。擅长Python.Matlab仿真、视觉处理、神经网络、数据分析。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

R语言Rstan概率编程规划MCMC采样的贝叶斯模型简介

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

学校数据集

定义Stan模型文件

变量

R语言stan进行贝叶斯推理分析

程序

随时关注您喜欢的主题

模型

注意，θ 永远不会出现在参数中。这是因为我们没有显式地对θ进行建模，而是对η（各个学校的标准化效果）进行了建模。然后，根据μ，τ和η在变换后的参数部分构造θ 。此参数化使采样器更高效。

准备数据进行建模

从后验分布抽样

模型解释

基本模型解释

MCMC诊断

层次回归

此表示假设所有样本都具有相同的分布。如果只存在一组样本，那么我们就会遇到问题，因为将忽略组内和组之间的潜在差异。

数据集

层次回归模型的规范

资料准备

拟合回归模型

层次回归模型的预测

相关文章

R语言Rstan概率编程规划MCMC采样的贝叶斯模型简介

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

学校数据集

定义Stan模型文件

变量

R语言stan进行贝叶斯推理分析

程序

随时关注您喜欢的主题

模型

注意，θ 永远不会出现在参数中。这是因为我们没有显式地对θ进行建模，而是对η（各个学校的标准化效果）进行了建模。然后， 根据μ，τ和η在变换后的参数部分构造θ 。此参数化使采样器更高效。

准备数据进行建模

从后验分布抽样

模型解释

基本模型解释

MCMC诊断

层次回归

此表示假设所有样本都具有相同的分布。如果只存在一组样本，那么我们就会遇到问题，因为将忽略组内和组之间的潜在差异。

数据集

层次回归模型的规范

资料准备

拟合回归模型

层次回归模型的预测

相关文章

关注有关新文章的微信公众号

注意，θ 永远不会出现在参数中。这是因为我们没有显式地对θ进行建模，而是对η（各个学校的标准化效果）进行了建模。然后，根据μ，τ和η在变换后的参数部分构造θ 。此参数化使采样器更高效。