R语言RStan贝叶斯示例:重复试验模型和种群竞争模型Lotka Volterra

Stan是一种用于指定统计模型的概率编程语言。Stan通过马尔可夫链蒙特卡罗方法(例如No-U-Turn采样器,一种汉密尔顿蒙特卡洛采样的自适应形式)为连续变量模型提供了完整的贝叶斯推断。

由Kaizong Ye,Liao Bao撰写

可以通过R使用rstan 包来调用Stan,也可以 通过Python使用 pystan 包。这两个接口都支持基于采样和基于优化的推断,并带有诊断和后验分析


在本文中,简要展示了Stan的主要特性。还显示了两个示例:第一个示例与简单的伯努利模型相关,第二个示例与基于常微分方程的Lotka-Volterra模型有关。


什么是Stan?

  • Stan是命令式概率编程语言。
  • Stan程序定义了概率模型。
  • 它声明数据和(受约束的)参数变量。
  • 它定义了对数后验。
  • Stan推理:使模型拟合数据并做出预测。
  • 它可以使用马尔可夫链蒙特卡罗(MCMC)进行完整的贝叶斯推断。
  • 使用变分贝叶斯(VB)进行近似贝叶斯推断。
  • 最大似然估计(MLE)用于惩罚最大似然估计。

 

×

先验概率与后验概率

1、先验概率 – 定义:直观理解,所谓“先”,就是在事情之前,即在事情发生之前事情发生的概率。是根据以往经验和分析得到的概率。 – 例子:比如抛硬币,我们都认为正面朝上的概率是0.5,这就是一种先验概率,在抛硬币前,我们只有常识。这个时候事情还没发生,我们进行概率判断。所谓的先验概率是对事情发生可能性猜测的数学表示

2、后验概率 – 定义:事情已经发生了,事情发生可能有很多原因,判断事情发生时由哪个原因引起的概率。 – 例子:比如今天你没去学校,原因有两个,可能是生病了,也可能是自行车坏了。然后上课时老师发现你没来。(这里是一个结果,就是你没来学校这件事情已经发生了)老师叫同学计算一下概率,分别是因为生病了没来学校的概率和自行车坏了没来学校的概率。

贝叶斯定理

开始谈贝叶斯定理之前,必须首先从条件概率 (conditional probability) 说起。所谓条件概率,就是在一个事件发生的情况下,去判断另一个相关联的事件发生的概率,或者简单说,就是指在事件 B 发生的情况下,事件 A 发生的概率。通常记为 P(A|B) 。接下来对贝叶斯公式做一个简单的推导,根据概率知识,我们可以求得 P(A|B) 为:

[公式]

同理可得:

[公式]

所以可以得到:

[公式]

最终可以得到大名鼎鼎的贝叶斯公式:

[公式]


实例、偷窃与狗吠

住在一座别墅里的一家人,在过去的 20 年中,发生了 2 次盗窃,这家人养了一只狗,这只狗平均每周晚上叫 3 次,而且,当发生盗窃时,这只狗会叫的概率是 0.9,那么问题是:在这只狗叫的情况下,发生盗窃的概率是多少?

我们假设 A 事件是狗晚上叫,则 P(A) = 3/7 ,假设 B 事件是发生盗窃,则P(B) = 2/365*20 。我们还知道,当 B 事件发生的条件下,A 事件发生的概率 P(A|B) = 0.9 。因此,根据贝叶斯公式,我们推断得到,狗叫时,发生盗窃的概率,即:

[公式]



Stan计算什么?

  • 得出后验分布 。
  • MCMC采样。
  • 绘制顺序,其中每个绘制都按后验概率的边缘分布。
  • 使用直方图,核密度估计等进行绘图

安装 rstan

要在R中运行Stan,必须安装 rstan C ++编译器。在Windows上, Rtools 是必需的。

最后,安装 rstan

Stan中的基本语法


视频

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

定义模型

Stan模型由六个程序块定义 :

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

数据块读出的外部信息。

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

 

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

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

 

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

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

类型

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

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

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

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


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

阅读文章


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

关于Stan的更多信息


随时关注您喜欢的主题


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

有两种修改 后验的方法。

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

贝叶斯方法

概率是 认知的。例如, 约翰·斯图亚特·米尔 (John Stuart Mill)说:

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

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

概率可以量化不确定性。

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

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

步骤1:问题定义

在此示例中,我们将考虑以下结构:

  • 数据:
    • ,试用次数。
    • ,即试验n的结果  (已知的建模数据)。
  • 参数:
  • 先验分布
  • 概率
  • 后验分布

步骤2:Stan

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

步骤3:数据

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

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

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

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

RStan:MAP,MLE

Stan的估算优化;两种观点:

  • 最大后验估计(MAP)
  • 最大似然估计(MLE)。

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

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

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

数学模型

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

这里:

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

stan中的Lotka-Volterra

观察到已知变量:

  • :表示在时间 物种数量

必须推断未知变量):

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

假设误差是成比例的(而不是相加的):

等效:

建立模型

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

未知参数的变量。

先验分布和概率。

我们必须为预测的总体定义变量 :

  • 初始种群(z0)。
  • 初始时间(0.0),时间(ts)。
  • 参数(theta)。
  • 最大迭代次数(1e3)。

Lotka-Volterra参数估计

获得结果:

分析所得结果:

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


可下载资源

关于作者

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

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

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

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