后验几率分布的抽样

【欢迎转发分享，转载请注明出处】

不少介绍贝叶斯后验几率的书中都经常使用医学化验为例子，好比一种方法的检出阳性率是多少之类的。 假设有一种检验方法有95%的灵敏度可以识别出一我的是否是吸血鬼vampirism,即Pr(+|vampire) = 0.95，也就是针对一个吸血鬼可以95%的可能性识别出来。固然，它也有误诊的时候，Pr(+|mortal) = 0.01，即对一个正常人，也有1%的可能被诊断为吸血鬼。此外，在人群中，咱们还知道，吸血鬼占的比例是很小很小的，只有0.1%，即Pr(vampire) = 0.001。在这种状况下，若是一个“人”被诊断为阳性，那么他是正常人仍是吸血鬼？几率是多少？若是咱们使用贝叶斯的方法，不难作出以下判断：

Pr(vampire|+)=Pr(+|vampire)Pr(vampire)/Pr(+)

其中，

Pr(+)=Pr(+|vampire)Pr(vampire)+Pr(+|mortal)(1−Pr(vampire))

在R语言中，咱们能够经过下面的过程来计算：

PrPV <- 0.95
PrPM <- 0.01
PrV <- 0.001
PrP <- PrPV*PrV + PrPM*(1-PrV)
(PrVP <- PrPV*PrV / PrP)

## [1] 0.08683729

从上面的结果能够看出，若是一个“人”被诊断为吸血鬼，那么他只有8.7%的可能性是一个真正的吸血鬼。这和咱们的感受有些相悖，这么高灵敏度(95%)的方法，只有8.7%的可能性。放在现实生活中，不少这样的检验，好比HIV检验等。这种检验的一个特色是，群体中阳性个体不多，因此即使是检验结果是阳性，也不能保证该个体就是患者。 不过若是我咱们换一个角度来看上面的例子，就很容易理解了，咱们不用几率问题，而是用具体数值。好比 （1）一个100,000的群体中，有100人是吸血鬼 （2）每100个吸血鬼可以检验出95个阳性 （3）在剩下99,900人中，有999我的的检验也是阳性 这时，若是对全部100,000我的都进行检验，那么有多少阳性结果是真正吸血鬼？

Pr(vampire|+)=95/(95+999) =0.087

这个时候就容易理解多了。与一些抽象的几率问题相比，人们对具体计数的问题更易于理解。 本章咱们继续研究几率分布的问题，不过咱们会在后验几率中随机抽样，经过具体计数的方式来理解几率分布。 为何要用抽样呢？首先，使用抽样能够避免复杂的微积分等数学问题，虽然不少科研研究人员有必定的数学基础，可是对于微积分依然敬而远之。而经过抽样，就能够把微积分问题为题转化为频率问题，接下来对频率问题的处理对不少人就驾轻就熟了。其次，不少计算后验几率的方法并无具体的公式，好比使用MCMC获得的后验几率就是经过随机抽样获得的，因此没法使用微积分来描述后验几率分布。

在网格估计的后验几率分布中抽样

以前介绍过网格估计法。下面是经过R语言，对掷地球实验的后验几率的网格估计，采用1000个网格点。

p_grid <- seq(from=0 , to=1 , length.out=1000)
prior <- rep(1 , 1000 )
likelihood <- dbinom( 6 , size=9 , prob=p_grid)
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)

下面咱们在上述模拟的后验几率分布中进行10000次随机抽样。

samples <- sample(p_grid, size = 10000, prob = posterior, replace = T)
plot(samples)

上图是这10000次随机抽样的分布状况。能够看出在几率为0.6的时候抽样点密度很大，在<0.25的时候，抽样点不多。咱们换一种形式来看这个图，改为密度分布图，以下：

library(rethinking)

## rethinking (Version 1.59)

dens(samples, xlab = "Proportion water(p)")

这时这张经过抽样估计的后验几率分布图就和咱们真实获得的后验几率分布图很类似了。并且随着抽样样本的增大，抽样对后验几率的估计会也来越精确。

抽样结果的应用

一旦你的模型获得了后验几率分布，那么模型的任务就算完成了。可是你的任务才算刚刚开始。你须要进一步对后验几率分布进行归纳总结，以及解释。具体要根据你的研究目的来决定，通常常见的问题有：

在某一界值内的后验几率是多少？ 在两个界值参数之间的后验几率是多少？ 选取什么界值会使获得后验几率分布的5%？ 等等。。。

1. 界值区间

在掷地球的例子中，若是我要问海洋覆盖率<0.5的后验几率是多少？那么，经过网格估计法获得的后验几率中，把小于0.5的部分的后验几率加起来就能够了。

sum(posterior[p_grid<0.5])

## [1] 0.1718746

也就是说，海洋覆盖率小于50%的后验几率是17%左右。看起来是很简单的，可是这是咱们用网格估计法获得的，而网格估计的缺点就是当参数特别多的时候，运算量迅速增长，因此应用范围很窄。下面咱们使用抽样的方法从新回答上面的问题，即只要把左右抽样样本中后验几率<0.5的样本相加，而后除以总抽样次数(10000次)，就能够估算出来咱们想要的结果。

sum(samples < 0.5)/10000

## [1] 0.1704

这就和咱们网格估计法的结果很接近了。(下图左)

一样地，若是是海洋覆盖率在50%-75%的后验几率，（下图右）用抽样估计以下：

sum(samples > 0.5 & samples < 0.75)/10000

## [1] 0.5981

2. 置信区间

好比咱们要想知道海洋覆盖率是多少的时候，后验几率的累积达到80%，这个时候用网格估计法就不太方便了，而使用抽样则很简单。下图左

quantile(samples, 0.8)

##       80% 
## 0.7647648

一样的，若是咱们想求中间80%的后验几率分布（下图右），能够

quantile(samples, c(0.1,0.9))

##       10%       90% 
## 0.4504505 0.8178178

上述后验几率都是近似正态分布的，因此经过百分区间就可以反映出整个分布的状况。若是后验几率分布不是正太分布，那么就有问题了。好比，咱们掷地球，前3次都掷到了海洋，这是，咱们使用网格进行后验几率分布，而后再估计中间50%的置信区间：

p_grid <- seq( from=0 , to=1 , length.out=1000 )
prior <- rep(1,1000)
likelihood <- dbinom( 3, size=3 , prob=p_grid )
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)
samples <- sample( p_grid , size=1e4 , replace=TRUE , prob=posterior )
quantile(samples,c(0.25,0.75))

##       25%       75% 
## 0.7087087 0.9309309

中间50%的百分区间在0.70-0.93之间，可是这个区间内并无包含最高的后验几率分布。在本例中，最高的后验几率分布在p=1处。以下图左。若是经过此区间来描述后验几率分布就存在严重的误导。

这个时候咱们能够用另外一个参数来描述，最高后验几率密度区间(HPDI)，它是包含必定累积几率的最窄区间(如上图右）。经过rethinking包中的HPDI函数能够来求解：

HPDI(samples, prob = 0.5)

##      |0.5      0.5| 
## 0.8418418 0.9989990

这是这个区间就会包含后验几率最大的值，此时的区间大小是最窄的，只有0.998-0.841=0.157。

百分区间和HPDI对于近似正态分布的后验几率分布来说，结果是很类似的，差异就是非正态分布的状况。HPDI在估计的时候计算量比较大，并且对抽样大小很敏感。

总之，在百分区间和HPDI估计差异很大的时候，那么你最好可以把整个后验分布几率的状况都画出来，以避免产生误导。

为何是95%？
咱们平时最经常使用到的就是95%的置信区间，也就是咱们的参数有5%的可能性不在这个区间内。一般统计学检验的水平也是5%，即p<0.05。除了习惯外，恐怕没有更科学的解释了。这个锅可能还要Ronald Fisher来背：“P=0.05的标准误差是1.96，大概接近2；用这个数做为判断标准误差是否有统计学差别是很方便的。”尽管不少人并不认同由于方便就使用0.05，但它彷佛成了约定俗成的一个判断标准。在Fisher的晚年，他也很积极的反对使用0.05做为界值。
可是，若是不用0.05，那么要用什么数来做为判断标准呢？目前彷佛尚未一致的结论。不过，你要知道0.05并非惟一的选择。

置信区间意味着什么？
咱们一般认为95%的置信区间就是真实的参数有95%的可能性落在这个区间范围以内。对于严格非贝叶斯者来讲，这是不对的，由于他们不容许经过几率来描述参数的不肯定性。他们认为，若是你大量重复这个实验，那么你获得区间的95%范围会包含这个真实的参数。也就是参数是固定的。关于95%置信区间的解释，不少人也没有弄清楚，有些用了非贝叶斯观点来解释95%置信区间，有些则用了贝叶斯几率来解释。
无论你是贝叶斯者仍是非贝叶斯者，95%的置信区间就意味着真的是95%吗？实际上是咱们过于自信了。还记得以前提到的“小世界和大世界”吗？这个95%仅仅是对你的“小世界”，即模型，来讲的；对于真实的“大世界”，未必是你推测的95%。

3. 点估计

使用点估计来描述贝叶斯后验几率分布并非一个很好的选择，一个点很难来描述整个分布状况。若是你就是想用一个点来描述呢？思考咱们以前提到的例子，在掷地球的实验中，咱们前三次连续观察到了海洋，你能够选择的点有最大数、平均数、中位数。

p_grid[which.max(posterior)]

## [1] 1

mean(samples)

## [1] 0.8010685

median(samples)

## [1] 0.8433433

这三个数差异很大，你应该选哪个呢？该怎么选择？

咱们能够选择一个损失函数loss function，而后基于整个分布和咱们选择的估计点，用该损失函数来估计选择该估计点形成的损失大小。统计学者和博弈论者对损失函数一直很感兴趣，不过在科学界用的并很少。那么什么是损失函数呢 举个例子，我要和你打赌，赌海洋表面覆盖率具体是多少，你告诉我一个你认为最正确的地球表面海洋覆盖率，若是你回答的正确，我就给你100美圆。若是你回答的不对，那么我会根据你的回答和正确答案之间的差异大小，扣除相应的钱。也就是越接近正确答案，你获得钱越多，反之越少。你由于没有答对而形成损失大小和你的答案距离真实值的大小成正比。 若是你有后验几率分布，你确定要根据后验几率分布来选择一个损失钱最少的答案。假如你选择了海洋覆盖率p=0.5做为你的答案，那么咱们根据后验几率分布计算一下此时损失多少：

sum(posterior*abs(0.5-p_grid))

## [1] 0.3128752

也就是，若是你选择0.5做为答案，你可能会损失3成的钱。那么选多少做为答案会损失最少呢？咱们能够把因此得可能都计算一遍。

loss <- sapply(p_grid, function(d) sum(posterior*abs(d - p_grid)))
p_grid[which.min(loss)]

## [1] 0.8408408

整个损失函数以下图：

最小值出如今0.84附近，也就是若是你回答0.84，那么你损失的钱会是最少的。

因此为了找一个相对准确的点估计值，你必须有一个损失函数来评估不一样选择形成的损失大小。可是不一样的损失函数可能获得不一样的结果。上述例子中，用差值的绝对值做为损失函数，获得的点估计值很接近中位数，若是选择差值的平方做为损失函数，可能获得的点估计值就很接近均数了。固然，若是后验几率分布为近似正态分布，那么这两个损失函数最后估计的结果可能差异不大。 在科学领域，人们彷佛还不太习惯使用损失函数，通常做点估计的时候一般习惯性的用均值或者中位数等。而更科学合理的方法是经过损失函数来选择一个最合适的点。

经过抽样进行模拟预测

1. 哑数据

在以前提到的掷地球的实验中，咱们经过实验，不只可以推断出后验几率分布，还可以模拟出模型的预测。似然函数是能够双向工做的，既能够经过观测来推断每个观测发生的可能性；也能够在给定参数下定义全部可能的观测分布，以此进行抽样模拟或预测。

咱们把模拟获得的数据称之为哑数据dummy data，即咱们真实观测数据的“替代品”。对于掷地球实验的哑数据，能够经过下面的二项分布似然函数来模拟。

其中w是观测到海洋的次数，n是总的投掷数。当n=2，即咱们投掷两次，咱们可能有3种结果：0次海洋、1次海洋和2次海洋。若是地球表面海洋覆盖率是0.7，那么咱们观测到这3种结果的可能性为：

dbinom(0:2,size = 2, prob = 0.7)

## [1] 0.09 0.42 0.49

也就是咱们有9%的几率观测到0次海洋，42%的可能观测到1次海洋，49%的可能观测到2次海洋。 上面是咱们经过似然函数模拟获得观众观测的几率，若是咱们想获得具体的模拟观测数据，即哑数据，则可使用rbinom函数来实现：

rbinom(10, size = 2, prob = 0.7)

##  [1] 1 2 1 2 2 1 2 1 2 1

即模拟十次实验，咱们获得0次0观测海洋，5次1观测到海洋，5次2观测海洋。 若是进行10000次模拟，每次模型投掷9次，则

dummy_w <- rbinom(10000, size = 9, prob = 0.7)
simplehist(dummy_w, xlab = "dummy water count")

经过使用似然函数进行模拟观测，能够看到咱们的实际观测结果很类似。

2. 模型核查

在上面的例子中，若是模拟观测结果和咱们实际观测结果差异很大，那么就要注意了，可能在某些地方会存在错误。此外，还应该核查一个模型是否充足。也就是咱们还要检查那些没可以被咱们模型所描述的数据，以便咱们可以对模型进行修正和改进。

没有模型是十全十美的，总会有一些使人不满意的地方，咱们所要判断的就是这些地方是不是很重要的地方。没有人但愿本身建造的模型仅仅可以描述现有数据，而不能进行数据预测。模型预测出现的不完美在所不免。有不完美的地方不可怕，只要咱们知道哪儿不完美就能够了。

再回到咱们的例子中，在上面模拟抽样的过程当中，咱们用了抽样几率0.7，咱们实际上使用对后验几率进行了点估计，以此来进行抽样的。可是，实际上，后验几率是一个在0-1之间的分布，而不是一个点。若是用点估计，那么咱们就会丢失不少信息，形成“过分自信”。

实际上在抽样过程当中，有两个不肯定性的来源，一个是观测的不肯定性，除非海洋覆盖率为0或者为1，不然咱们没法预测某一次投掷的结果。另外一个不肯定性来自参数不肯定性，也就是海洋覆盖率自己，在0-1之间，咱们没法肯定它真实的是哪个值。

若是咱们在一个参数不肯定的后验几率分布中进行抽样，那么咱们须要考虑参数的全部可能值，对于每个可能值，咱们获得获得一个抽样分布，最后把全部可能值得抽样分布进行平均，获得后验几率预测分布。以下图经过9个参数值的估计获得的后验几率预测分布：

咱们获得的后验几率预测分布是用来预测的，它考虑了参数不肯定性对预测结果的影响，因此预测会更准确一些。经过这个分布图，和上面咱们用0.7获得的分布图比较能够看出，这儿获得的分布图更加矮胖，方差更大。而若是你要是用0.7估计获得的高瘦的分布图来预测，就会出现“过分自信”的结果。

那么咱们怎样把参数的不肯定性考虑进去，来进行抽样呢？很简单，把上面的0.7换成咱们的后验几率分布就能够了：

w <- rbinom(10000, size = 9, prob = samples)
simplehist(w, xlab = "number water count")

其中的“sample”是咱们在后验几率分布中抽样的结果，表明了后验几率的分布状况。经过上面的结果，能够预测咱们9次抽样最有可能看到6次结果为“海洋”。和咱们实际观测的异质性较高。

到这儿咱们的对模型的核查就完了吗？尚未，还要考虑其余可能影响模型的缘由。好比，万一各个实验之间不是相互独立的怎么办？前一个实验的结果会影响到后一次实验。为了衡量独立性，咱们考虑两个方面：一个是实验中连续出现“海洋”或者“陆地”的最长长度；另外一个是相邻实验间，“陆地”和“海洋”之间的转换次数。 在咱们的实际实验中，咱们的观测结果是依次是“W L W W W LW L W”（W-海洋，L-陆地），因此上述两个指标在本实验中分别是3和6。一样咱们经过模拟，能够获得在相互独立这一预期状况下，这两个参数的分布状况。以下如：

【蓝色表示实际观测值】

在相互独立的假设下，咱们观测的最长长度为3和模拟的结果很吻合；而相互转化数咱们观测值为6，而预期最可能的值为4，这说明咱们的相邻实验结果之间可能互斥性。在这种状况下，各次实验之间存在的互斥性，可能使得每次实验结果带来的信息量大大减少，固然，咱们最终可能也能获得比较准确地海洋覆盖率的估计值，可是可能须要更多的试验次数。

总结

本部份内容介绍了关于后验几率分布的一些操做。咱们的目的是从后验几率分布中进行抽样，这样咱们就避免了使用微积分等复杂的方式来描述后验几率分布了。这些在后验几率中抽得的样本能够进一步用来区间估计、点估计、后验几率预测核查以及模拟。在点估计的时候，咱们还提到了一个重要概念，即损失函数。

后验几率预测核查综合考虑了有后验几率分布致使的参数不肯定性和似然函数定义的抽样自己的不肯定性，以此来核实结果以及判断模型是否充分。

随着模型自己变得愈来愈复杂，后验几率预测模拟的使用愈来愈普遍。

---