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

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

由Kaizong Ye,Liao Bao撰写

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

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

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

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


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

其实这个模型可以分两层,也可以分很多很多层。如果仅考虑两层的话,类似于此: [公式]

三阶段的话,则会像是: [公式]

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

总归得先看公式的.

如果还没忘记可交换性的概念的话(忘记了的请举手,并自觉往上看),从一个群体中独立采样的结果是可以看作是可交换的。这样的话群体的属性就可以用一个固定的未知参数 [公式] 来描述,即: [公式]

从而当考虑分层数据 [公式] ,其中[公式] 时,有 [公式]

而对于 [公式] 同样而可以假定: [公式]

这样就有了三个模型

组内: [公式]

组间: [公式]

先验分布: [公式]


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

组内模型: [公式]

组间模型: [公式]


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

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

为了计算 [公式] ,取一列服从 [公式] 上均匀分布的随机变量 [公式] ,则 [公式] 且 [公式]

从而 [公式] ,而由大数定理 [公式]

而当结果为 [公式] 时,对 [公式]

MCMC入门

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

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

  • 初始化马尔科夫链初始状态 [公式]

  • 对 [公式] 循环以下过程

  • 第 [公式] 个时刻马氏链状态为 [公式] ,采样 [公式]

  • 从均匀分布采样 [公式]

  • 如果 [公式] 则接受转移 [公式] ,即 [公式]

  • 否则不接受转移,即 [公式]


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

通过构造转移矩阵Q [公式]

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

  • 随机初始化 [公式]

  • 对 [公式] 循环采样

  • [公式]

  • [公式]

  • [公式]

  • [公式]

  • [公式]

  • [公式]

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



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

学校数据集

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

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

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


视频

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

探索见解

根据该模型,教练的效果遵循正态分布,其均值是真实效果θj,其标准偏差为σj(从数据中得知)。真正的影响θj遵循参数μ和τ的正态分布。

定义Stan模型文件

在指定了要使用的模型之后,我们现在可以讨论如何在Stan中指定此模型。在为上述模型定义Stan程序之前,让我们看一下Stan建模语言的结构。

变量

在Stan中,可以通过以下方式定义变量:

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

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


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

阅读文章


程序 

Stan中使用以下程序 :

  • data:用于指定以贝叶斯规则为条件的数据
  • 转换后的数据:用于预处理数据
  • 参数  (必填):用于指定模型的参数
  • 转换后的参数:用于计算后验之前的参数处理
  • 模型  (必填):用于指定模型
  • 生成数量:用于对结果进行后处理


随时关注您喜欢的主题


对于  模型  程序块,可以两种等效方式指定分布。第一个,使用以下统计符号:

第二种方法使用基于对数概率密度函数(lpdf)的程序化表示法:

Stan支持大量的概率分布。通过Stan指定模型时,该  lookup 函数会派上用场:它提供从R函数到Stan函数的映射。考虑以下示例:

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

模型

现在,我们了解了Stan建模语言的基础知识,我们可以定义模型,并将其存储在一个名为的文件中  schools.stan

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

准备数据进行建模

在拟合模型之前,我们需要将输入数据编码为一个列表,其参数应与Stan模型的数据部分相对应。对于学校数据,数据如下:

从后验分布抽样

我们可以使用stan 函数从后验分布中采样,函数执行以下三个步骤:

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

如果  rstan_options(auto_write = TRUE),则相同模型的后续调用将比第一次调用快得多,因为该  stan 函数随后跳过了前两个步骤(转换和编译模型)。

此外,我们将设置要使用的内核数:

现在,我们可以从后验中编译模型和样本。

模型解释

我们将首先对模型进行基本解释,然后研究MCMC程序。

基本模型解释

要使用拟合模型执行推断,我们可以使用  print 函数。

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

我们可以使用以下plot 函数来可视化估计中的不确定性  :

黑线表示95%的间隔,而红线表示80%的间隔。圆圈表示平均值的估计。

我们可以使用以下extract 函数获取生成的样本  :

MCMC诊断

 通过绘制采样过程的轨迹图,我们可以确定采样期间是否出了问题。例如,链条在一个位置停留的时间过长或在一个方向上走了太多步,就会有问题。我们可以使用traceplot 函数绘制模型中使用的四个链的轨迹  :

要从各个马尔可夫链中获取样本,我们可以extract 再次使用函数:

为了对采样过程进行更高级的分析,我们可以使用该  shinystan 软件包 。使用该软件包,可以通过以下方式启动Shiny应用程序来分析拟合模型:

层次回归

现在,我们对Stan有了基本的了解,我们可以深入研究更高级的应用程序:让我们尝试一下层次回归。在常规回归中,我们对以下形式的关系进行建模

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

另一种选择是为每个组建立一个回归模型。但是,在这种情况下,估计单个模型时,小样本量会带来问题。

层次回归是两个极端之间的折衷。该模型假设组是相似的,但存在差异。

假设每个样本都属于K组之一。然后,层次回归指定如下:

其中Yk是第k组的结果,αk是截距,Xk是特征,β(k)表示权重。层次模型不同于其中Yk分别拟合每个组的模型,因为假定参数αk和β(k)源自共同的分布。

数据集

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

让我们调查数据:

 数据显示线性增长趋势对于不同的大鼠非常相似。但是,我们还看到,大鼠的初始体重不同,需要不同的截距,并且生长速度也需要不同的斜率。因此,分层模型似乎是适当的。

层次回归模型的规范

该模型可以指定如下:

第i个大鼠的截距由αi表示,斜率由βi表示。注意,测量时间的中心是x = 22,它是时间序列数据的中值测量值(第22天)。

现在,我们可以指定模型并将其存储在名为 rats.stan的文件中 :

请注意,模型代码估算的是方差(  sigmasq  变量)而不是标准差。 

资料准备

为了准备模型数据,我们首先将测量点提取为数值,然后将所有内容编码为列表结构:

拟合回归模型

现在,我们可以为老鼠体重数据集拟合贝叶斯层次回归模型:

层次回归模型的预测

在确定了每只大鼠的α和β之后,我们现在可以估计任意时间点单个大鼠的体重。在这里,我们寻找从第0天到第100天的大鼠体重。

与原始数据相比,该模型的估计是平滑的,因为每条曲线都遵循线性模型。研究最后一个图中所示的置信区间,我们可以看到方差估计是合理的。我们对采样时(第8至36天)的老鼠体重充满信心,但是随着离开采样区域,不确定性会增加。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

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

​非常感谢您阅读本文,如需帮助请联系我们!

 
QQ在线咨询
售前咨询热线
15121130882
售后咨询热线
0571-63341498