【R】基于R实现贝叶斯分析（一）

贝叶斯简介

贝叶斯学派和频率学派的核心区别：用信念度而不是长期频率来衡量某件事发生的概率，这个信念度可能会随着新信息的加入而发生变化。贝叶斯的三个核心概念：

• 先验概率（prior）：在没有进行任何尝试时一个人对某件事的信念，例如人群中X病毒携带率为3%；
• 似然（likelihood）：观测到的数据，例如A检测结果为阳性；
• 后验（posterior）：通过观测到的数据对某件事的信念进行修正，例如A携带X病毒的概率为75%。

贝叶斯法则：根据数据更新你关于模型的概率的公式
$$P(MODEL\mid DATA)\propto P(MODEL)\times P(DATA\mid MODEL)$$

Why R

众多贝叶斯计算引擎都可以在R中直接调用。

贝叶斯计算引擎：
贝叶斯引擎是用于贝叶斯统计建模和推断的工具，这些引擎能帮助用户构建复杂的概率模型，并利用马尔可夫链蒙特卡罗（MCMC）等方法进行推断。

• BUGS
• WinBUGS/OpenBUGS
• JAGS（后续会更新相关文章专门介绍JAGS）
• Stan

用于贝叶斯分析/计算的R包：

1. MCMCpack: 包含针对常用模型的 MCMC 算法，通常在社会和行为科学中使用
2. coda: 包含一套用于收敛诊断和输出分析的函数，可用于总结、绘图和诊断 MCMC 样本的收敛性
3. rjags: 实现在 R 中完全交互式地使用JAGS的 R 包
4. r2jags: 基于 rjags，但提供了更为简便的接口，简化了 JAGS 的调用流程
5. runjags: 一个包含许多 JAGS 增强功能的 R 包，包括自定义 JAGS 模块
6. rstan: Stan 的 R 语言接口。它允许用户方便地从 R 中拟合 Stan 模型并访问输出结果
7. Nimble: 一个提供通用 MCMC 系统的 R 包，允许对用 BUGS/JAGS 模型语言编写的模型进行自定义 MCMC。用户可以选择采样器并编写新的采样器。模型和采样器通过生成的C++ 代码自动编译以提高速度
8. LaplacesDemon: 一个旨在提供完整贝叶斯环境的 R 包，包括众多 MCMC算法、具有多个优化算法的拉普拉斯近似、大量示例、数十种额外的概率分布以及众多 MCMC 诊断工具等

理论基础

单参数贝叶斯模型：仅包含一个待估计参数的概率模型的R语言建模实现，一个直观的例子：总体比例p的贝叶斯推断，希望从样本数据中得到p的点估计和区间估计（不确定性）。本文将从总体比例p估计问题出发，探讨单参数贝叶斯建模过程。

贝叶斯分析的三步走：

• 指定先验分布
• 写出似然函数
• 计算参数的后验分布

上面三个步骤中，我们需要研究的问题：

1. 先验分布：如何合理表达先验信念

• 选择什么先验分布方便推导？（无信息先验、层级先验、共轭先验等）

2. 似然函数：如何用分布描述数据生成过程

• 使用什么分布建模？（正态、二项、泊松等）

3. 后验分布：如何计算后验分布和后验预测

• 是否能够解析计算？（有共轭形式时可以解析；否则需数值方法如 MCMC）

对于总体比例估计问题，我们的似然函数就是二项分布，即独立重复试验中成功次数。若我们观察到y 次成功，在 n 次试验中：$y\sim Binomial(n,p)$，其似然函数为：
$$L(p)=P(y\mid p)=(n,y)p^y(1-p{)}^{n-y}$$

贝叶斯模型：
$$p(𝜃,y)=𝜋(𝜃)p(y|𝜃)$$

• $𝜋(𝜃)$是先验分布
• $p(y|𝜃)$是似然函数
• 贝叶斯模型即𝜃和y的联合分布$p(𝜃,y)$

我们关心的是参数的后验分布$𝜋(𝜃\mid y)$，它可以由贝叶斯公式计算出：
$$𝜋(𝜃\mid y)=p(𝜃,y)/p(y)=𝜋(𝜃)p(y\mid 𝜃)/p(y)$$
其中𝑝(𝑦) 不依赖于𝜃，当给定观测数据时可视作常数，从而可以看作是一个归一化常数，用于将𝜋(𝜃 ∣ 𝑦)化成有效的概率密度函数。

一、三种先验分布和对应后验的计算

构建先验分布的常见方法是假设先验分布属于特定的参数密度函数族。然后选择先验分布的参数（称为超参数 hyper parameters），使先验分布尽可能地表示先验信念。

在研究总体概率p时，有三种典型先验：1. 离散先验；2. Beta先验（共轭先验）；3. 直方图先验。本文下面将详细讲解这三种先验以及对应后验的计算。

1. 离散先验

假设𝑝只能取有限个值（如 0.1, 0.2, …, 0.9），我们给这些值一个先验概率分布：$P(p=p_i)=w_i$， $w_i$求和为1。

离散先验分布：

p = seq(0.05, 0.95, by = 0.1)
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior = prior/sum(prior)
plot(p, prior, type = "h", ylab="Prior Probability")

• p为概率值向量
• prior为p对应的先验概率向量，归一化处理
• type = "h"为直方图选项，用于绘制概率质量函数（pmf）

输出：

已知观测数据后计算后验分布：

library(LearnBayes)
data = c(11, 16)
post = pdisc(p, prior, data)
round(cbind(p, prior, post),2)

• data为观测数据：11次成功，16次失败
• pdisc用于计算离散先验下的后验概率，三个参数分别是概率点向量、先验概率、观测数据

输出：

对比先验和后验：

library(lattice)
PRIOR=data.frame("prior",p,prior)
POST=data.frame("posterior",p,post)
names(PRIOR)=c("Type","P","Probability")
names(POST)=c("Type","P","Probability")
data=rbind(PRIOR,POST)
xyplot(Probability~P|Type,data=data, layout=c(1,2), type="h", lwd=3,col="black")

• data.frame用于创建一个三列的数据框，names用于给数据框命名
  rbind表示将两个数据框按行合并
  xyplot是lattice包的核心函数，用于绘制散点图或线图，这里Probability~P|Type表示纵轴为Probability,横轴为P ，按照Type分面板绘制

输出：

2.Beta先验（共轭先验）

从上图介绍我们可以得到为什么要选择Beta先验：它和似然函数（二项分布）有相同形式的密度函数，这样似然*先验得到的后验还是这个形式。

问题：如何确定Beta分布的两个参数？

• 策略1: 画出来好多条密度函数，找一条满足期望的
• 策略2: 根据想要的均值/样本量确定
  - $Beta(\alpha ,\beta )$分布的均值为$\alpha /(\alpha +\beta )$
  - $Beta(\alpha ,\beta )$分布的等效先验样本量为$\alpha +\beta -2$（其中$\alpha -1$次成功，$\beta -1$次失败）

在R中可以用LearnBayes包中的beta.select()函数找到与先验知识（对分位数的认知）匹配的两个超参数。

quantile1=list(p=.5,x=.3) # 中位数
quantile2=list(p=.9,x=.5) # 90%分位数
beta.select(quantile1,quantile2)

输出：
[1] 3.26 7.19

• 先验信息：研究者认为比例 𝑝 的中位数为 0.3，90% 分位数为 0.5
• 选出的$\alpha =3.26$，$\beta =7.19$

假设我们选出了一个合适的先验$Beta(\alpha ,\beta )$，观测$n$次中有$y$次成功，就可以用贝叶斯公式求出后验：
$$\pi \propto p^y(1-p{)}^{n-y}{p}^{\alpha -1}(1-p{)}^{\beta -1}=p^{\alpha +y-1}(1-p{)}^{\beta +n-y-1}$$

这是$Beta(\alpha +y,\beta +n-y)$的核，这就是共轭性的含义：先验分布和后验分布来自同一种分布类型。

3. 直方图先验

就是每个区间内指定一个先验概率：

midpt = seq(0.05, 0.95, by = 0.1) # 按照0.1的间隔划分0.05到0.95的区间
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0) # 每个区间的先验概率
prior = prior/sum(prior) # 归一化
curve(histprior(x,midpt,prior),
      from=0, to=1,
      ylab="Prior density",
      ylim=c(0,.3))

输出：

计算后验分布一样，乘上似然函数的等价beta分布表示：

curve(histprior(x,midpt,prior) * dbeta(x,s+1,f+1),
      from=0, to=1,
      ylab="Posterior density")

输出：

二. 后验抽样

后验抽样是指在获取参数后验分布后，从后验分布中采样参数，进而用于后续贝叶斯推断。

1. 网格点采样法

网格点采样法，也称为网格近似法 (Grid Approximation)，是一种通过在参数空间中定义一个离散的网格，然后计算每个网格点上的未归一化后验概率（即似然函数值乘以先验概率密度值）来近似后验分布的方法。一旦计算出这些概率，就可以将它们归一化，使其总和为 1，从而得到一个离散化的后验概率分布，从这个离散分布中抽取样本。

p = seq(0, 1, length=500) # 网格点，500个0到1之间均匀的数值点
post = histprior(p, midpt, prior)*dbeta(p, s+1, f+1) # 以直方图先验为例，计算每个网格点的后验概率
post = post/sum(post) # 后验概率归一化
ps = sample(p, 5000, replace = TRUE,
            prob = post) # 从post分布冲抽取5000个样本
hist(ps, breaks= 30, xlab="p",main="")

输出：

2. 其他方法

后续会专门写一篇文章讲各种MCMC方法。

三、贝叶斯推断

通过确定先验分布后得到的后验分布包含了关于未知参数的所有当前信息，所有贝叶斯推断都基于后验分布，包括：

1. 参数估计
   • 点估计：后验均值 (posterior mean)、众数 (posterior mode / MAP)、中位数 (posterior median)
   • 离散程度：后验方差
   • 可信区间(credible interval)：后验分位数
2. 假设检验：后验概率
3. 预测：估计潜在可观测但当前未观测的量
4. 模型选择/变量选择：Bayes factor,DIC

1. 参数估计

(1) 后验均值

后验分布的均值常用作$p$点估计，$Beta(\alpha +y,\beta +n-y)$的均值为$\alpha +y/\alpha +\beta +n$

(2) 后验方差

后验方差是后验分布离散程度的一种度量，后验方差越大，我们对参数的不确定性就越大。

$Beta(\alpha ,\beta )$的方差为：$\alpha \beta /(\alpha +\beta {)}^2(\alpha +\beta +1)$

对于 $Beta(\alpha ,\beta )$先验和二项似然（n 次试验中有 y 次成功），后验分布为 $Beta(\alpha +y,\beta +n-y)$，从而后验方差为$(\alpha +y)(\beta +n-y)/(\alpha +\beta +n{)}^2(\alpha +\beta +n+1)$

注意：对于均匀先验和二项似然，后验方差（几乎）总是小于先验方差

(3) 后验区间

两种类型的可信区间：

1. 等尾可信区间(equal-tail credible interval):

100(1 − α)% 的等尾可信区间是从后验分布的 α/2 分位数到 (1 − α/2) 分位数的区间，若要构建 95% 的等尾可信区间，我们需要计算后验分布的 0.025 和 0.975 分位数。

例如，求后验均值的95%置信区间有两种方法，一是直接根据beta分布来计算：

a = 3.26; b = 7.19 # 选出的先验参数
s = 11; f = 16 # 观测数据（似然）
qbeta(c(0.05, 0.95), a + s, b + f)

输出：[1] 0.2555267 0.5133608

二是用模拟的方法：

ps = rbeta(1000, a + s, b + f)
quantile(ps, c(0.05, 0.95))

输出：5% 95%
0.2490099 0.5152768

• beta表示从beta分布中生成指定个随机数，
• 这里生成了1000个样本，然后计算样本的0.05和0.95分位数

对于贝叶斯学派，95% 可信区间的解释是：观察到观测数据后，真实参数 p 落在该区间内的概率为 0.95

2. 最高后验密度区间(highest posterior density, HPD):

区间内任意点的密度大于区间外任意点的密度，在包含指定概率的所有区间中长度最短（有时称为最小长度可信区间）。

当后验分布高度偏斜或多峰时，优于等尾可信区间

计算可以使用 TeachingDemos 包中的 R 函数hpd()或emp.hpd()

2. 假设检验

假设我们要检验关于参数 p 的以下假设：
$$H0:p\le 0.10v.sH1:p>0.10$$
我们只需计算后验分布在这两个区间上的概率，使用R函数 pbeta即可。

3. 预测

(1) 先验预测

以离散先验为例：

library(LearnBayes)
p = seq(0.05, 0.95, by=.1)
prior = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.8, 0.3, 0.1, 0.05)
prior = prior/sum(prior)
ns = 20
ys = 0:20
pred = pdiscp(p, prior, ns, ys)
round(cbind(0:20, pred), 3)

• 前面都和第一部分的离散先验一样
• ns表示未来样本量，ys表示未来可能的成功个数（从0～20）
• pdiscp用于计算离散先验下的后验预测分布，四个参数分别为概率点向量、先验概率、未来样本量、未来成功次数

输出：

(2) 后验预测

以Beta先验为例，有解析计算和模拟两种方法：

1. 解析计算：本质就是更新beta分布的两个参数

library(LearnBayes)
ab = c(3.26, 7.19)
ns = 20; ys = 0:20
pred = pbetap(ab, ns, ys)
round(cbind(0:20, pred), 3)

输出：

• pbetap函数计算 Beta 先验下的后验预测分布（Beta-Binomial 分布），三个参数分别是Beta先验的参数向量、未来试验次数、需要预测的成功次数

2.基于模拟：模拟的方法适合任何先验分布，下例展示了beta先验的后验模拟

m = 1000
a = 3.26; b = 7.19 # 先验Beta分布的参数
s = 11; f = 16
ns = 20
p = rbeta(m, a+s, b+f)
y = rbinom(m, ns, p)
freq = table(y) # 统计各成功次数y出现的频次
ys = as.integer(names(freq)) # 提取成功次数（去重）
predprob = freq/sum(freq) # y的概率分布
# 成功概率分布的直方图绘制
plot(ys, predprob, type="h", xlab="y",
ylab="Predictive Probability",
cex.lab=1, cex.axis=1.0)
# 解析解验证
library(LearnBayes)
TruePred1 = pbetap(c(a+s,b+f), ns, ys) # beta-二项后验分布的概率
points(ys, TruePred1, type="p", col='red')

• rbeta从后验分布中生成1000个p的随机样本
• rbinom对每个采样的 p，生成 ns=20 次试验的成功次数 y

输出：

如果后验是混合分布，那么后验预测分布是一个连续混合分布或复合分布，下例中展示了用模拟的方法从混合正态分布0.3N(0, 1) + 0.7N(3, 1/2)中随机抽样：

n <- 10000; alpha <- 0.3
k <- sample(1:2, size = n, replace = T, prob=c(alpha, 1-alpha)) # 根据概率分到1组或2组
mu <- numeric(length(k)) 
sd <- numeric(length(k)) 
mu <- c(0,3)[k] # 根据k选择均值（0或3）
sd <- c(1, 1/2)[k]# 根据k选择方差（1或1/2）
x <- rnorm(n, mu, sd) # 生成混合样本
hist(x, breaks='FD', prob=T, xlim=c(-4,6), ylim=c(0,0.8), main='Histogram of a mixture normal')
# 混合密度函数
p <- function(x,alpha){
alpha*dnorm(x, 0, 1)+(1-alpha)*dnorm(x, 3, 1/2)
}
curve(dnorm(x,0,1), col=4, lty=2, lwd=2, add=T) # N（0，1）
curve(dnorm(x,3,1/2), col=4, lty=4, lwd=2, add=T) # N（3，1/2）
curve(p(x,0.3), col=2, lty=1, lwd= 2, add=T) # 混合密度
legend("topleft", legend = c("N(0,1)","N(1,1/2)", "0.3N(0,1)+0.7N(1,1/2)"), cex=0.6, lty = c(2,4,1), col=c(4, 4, 2), lwd =c(2,2,2))

输出：