40 - MCMC方法用于后验分布抽样

2022-03-18 本文已影响0人 Hello育种

此内容整理来自：澳大利亚学校Ruidong Xiang and Mike Goddard的讲授课程。

1 介绍MCMC

Monto Carlo and Markov chain(MCMC) 中文为蒙特卡洛-马尔科夫链。

它允许贝叶斯在难以计算时近似后验分布的参数

它需要一些计算，但但在计算机上非常容易编码。

1.1 蒙特卡洛

解释：从一个分布中抽取随机样本来估计这个分布的性质
核心是：从某个方程中抽取大量随机数；然后计算这些数字的平均值

例子-计算蓝色区域的面积

第二个计算方法是我们不知道求圆面积的公式，也是蒙特卡洛的方法

image.png

1.2 马尔科夫链

生成随机样本的过程
每个随机样本都依赖于它之前的随机样本，这就是为什么它是一个“链”

1.3 实现步骤

最简明的说明：

1. 对分布的值有一个（好的）初始猜测
2.在这个初始值上添加一个随机噪声，为新样本做出一个“提议”
1. 决定如何处理提案：“接受”或“拒绝”；如果被拒绝，我们保留旧样本
1. 从第 2 步开始重复

有许多算法可以确定：

如何向旧样本添加随机噪声
如何决定接受或拒绝新样品的提案

Gibbs 采样器是这些算法之一，它属于 Metropolis-Hastings 类别

2 MCMC简单的例子（Metropolis-Hastings)

问题：我们想知道人类身高的平均值，我们可以使用 MCMC 找出人类身高的平均值吗？
假设：我们知道人类身高的标准差是6cm，并且观察到一名身高为 170 厘米的人 (N=1)

解决步骤（虚假代码）：

假设身高180cm为起始值
为了生成“提案”，我们将 N(0,5) 的正态分布中的随机数添加到起始值 180
我们比较正态分布（密度高度）中新值的概率（Pnew）和正态分布中旧值的概率（Pold）；正态分布~N(170,6)，因为我们用一个170cm人类身高的样本观察，所以mean = 170/1
如果 Pnew > Pold；然后接受新值
如果 Pnew <= Pold; 然后计算新值概率与旧值概率之比（Pnew/Pold）；然后给定 Pnew/Pold，随机决定是接受还是拒绝该提案
从第 2 步开始重复，重复 500 次

3 Gibbs simpler例子

Gibbs simple可以处理需要估计多个参数的问题，即多元分布。
其可以直接从参数的条件分布或给定另一个参数的特定值的参数的概率分布中为每个参数抽取样本
当来自后验条件的信息都是是已知的，对每个参数进行抽样

问题：我们有关于疾病测试的敏感性和特异性的历史数据。其试验是做了一些测试，发现了阳性测试的总数，我们想找出这种疾病的患病率，即 (a + c) / 1000？但是我们懒得去实验室确认谁真的生病或健康（a、b、c、d 的值），是否可以使用gibbs simpler估计出在人群中这种病的患病率呢？

已有数据如下：

image.png

解决步骤（伪代码）

我们使用 Gibbs 采样器从 beta 分布中获取流行率的点和区间估计值：𝐵𝑒𝑡𝑎（N 个病人，N 个健康人）
为此，我们无需去实验室，就将患病（或健康）的总人数分为 2 个部分：1）X 患病（或健康）人数呈阳性，2）Y 患病（或健康）人数人们测试为阴性：
𝑃𝑟𝑒𝑣𝑎𝑙𝑒𝑛𝑐𝑒|𝑋，𝑌〜𝐵𝑒𝑡𝑎（𝑎𝑙𝑝ℎ𝑎。+𝑋+𝑌，𝑏𝑒𝑡𝑎。+𝑁.𝑡𝑜𝑡𝑡𝑜𝑡𝑋）， X + Y = 患病总人数； 𝑁.𝑡𝑜𝑡−𝑋−𝑌 = 健康人总数 𝑎𝑙𝑝ℎ𝑎.𝑝𝑟𝑖𝑜𝑟 = 𝑏𝑒𝑡𝑎.𝑝𝑟𝑖𝑜𝑟 = 1
𝑋|𝑟，𝑃𝑟𝑒𝑣𝑎𝑙𝑒𝑛𝑐𝑒〜𝐵𝑖𝑜𝑛𝑜𝑚（𝑟，𝑝），其中r =样品数量的阳性（220）和p =𝑃𝑟𝑒𝑣𝑎𝑙𝑒𝑛𝑐𝑒* sensitivity / p.pos.test（p.pos.test = 220/1000）
𝑌|𝑟,𝑃𝑟𝑒𝑣𝑎𝑙𝑒𝑛𝑐𝑒~𝐵𝑖𝑜𝑛𝑜𝑚(𝑁.𝑡𝑜𝑡−𝑟,𝑞)，其中𝑁.𝑡𝑜𝑡 - q=𝑃𝑟𝑒𝑣𝑎𝑙𝑒𝑛𝑐𝑒* (1 –灵敏度)/ (1 - P.pos.test )

注意些真正代码时。首先需要写成X，Y的方程，再写患病率的方程。

总结

MCMC 是一种强大的分析方法，可在难以计算时从后验分布估计参数

Gibbs 采样器是一种有效的多变量分布 MCMC 算法

非常危险：估计没有足够数据信息的未知变量，
选择初始值，去除值，迭代数，需要测试收敛性，
或者运行不同的链（初始值不同），比较，如果结果类似，说明可行