贝叶斯推断简介

作者: 虚胖一场 | 来源:发表于2021-09-14 17:04 被阅读0次

贝叶斯推断简介
python实现贝叶斯推断——垃圾邮件分类
小强学AI之 - 2你患癌的概率（朴素贝叶斯）
【转】贝叶斯推断及其互联网应用---阮一峰
贝叶斯分类器（1）贝叶斯决策论概述、贝叶斯和频率、概率和似然
贝叶斯分类器(10X单细胞和10X空间转录组的基础算法)
【机器学习实战（四）】sklearn包中朴素贝叶斯库的使用
使用pymc3进行贝叶斯推断的简单示例
频率学派与贝叶斯学派
Python数据分析与机器学习26-贝叶斯算法

你是一名程序员，今天你实现了一个非常复杂的算法。你不知道这段代码有没有 bug。根据你以往的经验，写出 bug 是在所难免的，毕竟你不是 Jeff Dean。你先拿一些简单的测试用例 run 了一下，通过了。你又试了一些复杂的测试用例，也通过了。接下来，你连续试了一大堆用例，都通过了。于是你松了一口气，开始觉得这段程序可能没有 bug 了。

如果你这样思考问题，那么恭喜你，你是一个贝叶斯主义者。所谓贝叶斯推断，说白了，就是在观察到新的证据之后，更新你的信念。贝叶斯主义者往往不会对某件事确信无疑，但他可以对某件事很有信心。拿上面的例子来说，我们无法保证代码 100% 没有 bug（除非穷尽所有可能的情况）。但当代码通过了大量的测试用例之后，我们对它更有信心了。

我们经常用到“概率”这个词，但对什么是“概率”，其实是有不同看法的。

在贝叶斯学派看来，一个事件的概率是人们根据已有的知识和经验对该事件发生的可能性所给出的个人信念^[1]，这种信念用区间 $[0,1]$ 中的一个数来表示，可能性大的对应较大的数。
与之相对的，在频率学派看来，一个事件在独立重复实验中发生的频率会趋于一个极限，这个极限就是该事件的概率。

乍一看，频率学派的观点好像更严谨。举例来说，你想知道抛一枚硬币国徽朝上的概率，你就不停地抛这枚硬币，统计国徽朝上的次数的占比，随着你抛的次数的增多，你会发现这个占比渐渐趋于稳定，这个极限值就是国徽朝上的概率。问题是，有些情形下，你是无法做独立重复实验的。比如老板问你目前负责的某个项目有多大的概率能拿到业务成果。这个项目只能做一次，要么成功，要么失败，怎么计算成功的频率呢？对此，频率学派可能会告诉你，假设有很多个平行世界，这个项目在一些平行世界里成功了，而在另一些平行世界里失败了。这个项目成功的概率就是在平行世界里成功的频率的极限。但作为贝叶斯主义者，我们不需要诉诸虚无缥缈的平行世界。基于以往的工作经验和对这个项目现状的分析，我们依然能够给出我们的信念，“老板，我有 7 成把握，这个项目能成。”

我们把对事件 $A$ 发生的可能性的信念记为 $P(A)$ ，称作事件 $A$ 的先验概率（prior probability）。当观察到证据 $E$ 之后，我们更新了对事件 $A$ 的信念。更新后的信念记为 $P(A|E)$ ，称作事件 $A$ 的后验概率（posterior probability）^[2]。那么，在观察到证据 $E$ 之后，我们怎么才能更新先验概率来得到后验概率呢？很简单，我们使用贝叶斯定理^[3]：
$\begin{aligned} P(A|E) &= \frac{P(E|A)P(A)}{P(E)} \\ &\propto P(E|A)P(A) \end{aligned}$
$P(E|A)$ 是当事件 $A$ 发生时，我们观察到证据 $E$ 的概率。它是关于证据 $E$ 的函数，称作证据 $E$ 的似然函数（likelihood）。

例 1：2020 年 12 月 23 日，你看完勇士 VS 篮网的比赛，库里三分球 10 投 2 中。你想预测一下库里在本赛季的三分球命中率。因为这是库里本赛季的第一场比赛，直接用 $2/10=20\%$ 显然是不合理的。那么该怎么做呢？

查阅资料，我发现库里在最近的两个赛季三分球一共是 859 投 366 中。用 $\Theta$ 表示库里本赛季三分球的命中率。根据历史数据，我认为 $\Theta$ 服从参数为 $\alpha=366$ 和 $\beta=859-366=493$ 的 Beta 分布。即我对库里本赛季三分球命中率 $\Theta=\theta$ 的信念为
$P(\Theta=\theta)=\mathrm{Beta}(\theta;\alpha=366, \beta=493)$
Beta 分布的期望为
$\mathbb E(\Theta)=\frac{\alpha}{\alpha+\beta}$
所以，在开始的时候，我觉得库里本赛季的三分球命中率会在 $366/(366+493)\approx42.61\%$ 附近波动。

已知库里的三分命中率为 $\theta$ ，那么他投 $n$ 次三分球命中的次数 $X$ 应该服从参数为 $n$ 和 $\theta$ 的二项分布，即似然为
$P(X=x|\Theta=\theta)=\mathrm{Binomial}(x;n,\theta)$
在看到本场比赛库里 10 投 2 中的结果后，我就可以更新我的信念了。我们有
$\mathrm{Beta}(\theta;\alpha,\beta) = \frac{1}{\mathrm{B}(\alpha, \beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}$
其中
$\mathrm B(\alpha, \beta) = \int_0^1\theta^{\alpha-1}(1-\theta)^{\beta-1}\mathrm d\theta$
又有
$\mathrm{Binomial}(x;n,\theta) = \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x}$
因此
$\begin{aligned} P(\Theta=\theta|X=x) &= \frac{P(X=x|\Theta=\theta)P(\Theta=\theta)}{P(X=x)} \\ &= \frac{P(X=x|\Theta=\theta)P(\Theta=\theta)}{\int_0^1 P(X=x|\Theta=\theta)P(\Theta=\theta)\mathrm d\theta} \\ &= \frac{\mathrm{Binomial}(x;n,\theta)\mathrm{Beta}(\theta;\alpha,\beta)}{\int_0^1 \mathrm{Binomial}(x;n,\theta)\mathrm{Beta}(\theta;\alpha,\beta)\mathrm d\theta} \\ &=\frac{\frac{1}{\mathrm B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\cdot\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x}}{\int_0^1 \frac{1}{\mathrm B(\alpha,\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\cdot\frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \theta^x(1-\theta)^{n-x} \mathrm d\theta }\\ &=\frac{\theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}}{\int_0^1 \theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}\mathrm d\theta} \\ &=\frac{1}{\mathrm B(\alpha+x, \beta+n-x)}\theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}\\ &\equiv \mathrm{Beta}(\theta;\alpha+x,\beta+n -x) \end{aligned}$
于是，我的信念更新为
$\begin{aligned} P(\Theta=\theta|X=2) &= \mathrm{Beta}(\theta;\alpha=366+2,\beta=493+10-2)\\ &= \mathrm{Beta}(\theta;\alpha=368,\beta=501) \end{aligned}$
也就是说，在观察到库里首场比赛三分球 10 投 2 中之后，我更新了我的信念，认为他本赛季的三分球命中率将会在 $368/(368+501)\approx42.35\%$ 附近波动。

看到这些计算过程，有些读者可能已经崩溃了。然而这已经算是最简单的情形之一了。但是别着急，我们并不一定要把后验分布计算出来！我们所关心的无非是后验分布的期望、方差、偏度、峰度之类的统计量，如果我们能够设法从后验分布中采样，就可以利用这些样本来估计这些统计量。那么，有什么方法可以从一个复杂的概率分布中采样呢？非常幸运的是，马尔可夫链蒙特卡洛（Markov chain Monte Carlo, MCMC）方法可以做到这一点。关于这个方法，后续我会专门介绍。这里我们只讲一下实际操作。

在 Python 下，我们可以使用 PyMC3 轻松实现 MCMC 采样。拿上面的例子来说，

import pymc3 as pm
import seaborn as sns

with pm.Model() as model:
    # theta 的先验分布为 Beta 分布
    theta = pm.Beta(
        name='theta', 
        alpha=366, 
        beta=493
    )
    # 似然函数为二项分布
    likelihood = pm.Binomial(
        name='likelihood', 
        p=theta,
        n=10,
        observed=[2]
    )
    # 利用 MCMC 方法从后验分布中采样
    posterior = pm.sample()

# 利用后验分布的样本估计后验分布的期望
print(posterior.theta.mean())
# 0.4234174176263356

# 样本的直方图
sns.distplot(posterior.theta)

采样结果

例 2： 下图展示了某商品在一段时间内每天的销量。请问在这段时间内销量的平均水平是否发生过变化？

商品的销量

用 $S_t$ 表示第 $t$ 天的销量，这种计数类型的随机变量适合用 Poisson 分布来刻画，即似然为
$S_t\sim\mathrm{Poisson}(\lambda)$
假设在第 $\tau$ 天，销量的水平发生了变化，则
$\lambda=\begin{cases} \lambda_1,&t<\tau\\ \lambda_2,&t\geq\tau \end{cases}$
我们只需要给出 $\lambda_1$ 、 $\lambda_2$ 和 $\tau$ 的先验分布，然后利用 MCMC 方法从他们的后验分布中采样即可。
关于 $\lambda_1$ 和 $\lambda_2$ ，我们只知道他们一定是正数，且和这段时间的平均销量应该相差不多。因此，我们假设他们的先验分布为指数分布
$\begin{aligned} \lambda_1&\sim\mathrm{Exp}(\alpha)\\ \lambda_2&\sim\mathrm{Exp}(\alpha) \end{aligned}$
其中 $\frac{1}{\alpha}$ 为这段时间的平均销量。
关于 $\tau$ ，我们所知甚少，只能假设它的先验分布是离散均匀分布
$\tau\sim\mathrm{DiscreteUniform}(0, 89)$

同样地，我们使用 PyMC3 从后验分布中采样

with pm.Model() as model:
    alpha = 1.0 / sales.mean()
    lambda_1 = pm.Exponential('lambda_1', alpha)
    lambda_2 = pm.Exponential('lambda_2', alpha)
    
    tau = pm.DiscreteUniform('tau', lower=0, upper=len(sales)-1)
    
    t = np.arange(len(sales))
    lambda_ = pm.math.switch(t < tau, lambda_1, lambda_2)
    
    likelihood = pm.Poisson('likelihood', lambda_, observed=sales)
    
    posterior = pm.sample(20000, tune=8000, step=pm.Metropolis())

结果如下：

采样结果
根据采样的结果，

\lambda_1

、

\lambda_2

和

\tau

的均值分别是

\begin{aligned} \lambda_1 &\approx 17.96\\ \lambda_2 &\approx 24.00\\ \tau &\approx 50.40 \end{aligned}

销量的变化

这几乎就是正确答案了，因为所谓的“销量”实际上是我用下面的代码随机生成的

from scipy.stats import poisson
import numpy as np

sales = []

for _ in range(52):
    sales.append(poisson.rvs(mu=18))
for _ in range(38):
    sales.append(poisson.rvs(mu=23))

sales = np.array(sales)

参考资料

个人信念？那岂不是每个人都不同？的确如此，因为每个人拥有的知识和经验都不同，对同一件事的判断肯定不一样。举例来说，股(jiu)民(cai)和知道内部消息的大佬对某只股票涨跌的可能性就有不同的信念，正所谓“时间的朋友不如领导的朋友”。 ↩
注意到贝叶斯推断可以迭代使用：后验概率可以当作新的先验概率，从而在观察到更新的证据之后，再次更新我们的信念。 ↩
需要说明的是，这个公式并不是贝叶斯推断所独有的。学习过概率论的朋友一眼就能看出来，这个就是课本上的“逆概公式”。 ↩

贝叶斯推断简介
你是一名程序员，今天你实现了一个非常复杂的算法。你不知道这段代码有没有 bug。根据你以往的经验，写出 bug 是...
python实现贝叶斯推断——垃圾邮件分类
理论理论强推阮一峰大神的个人网站1.贝叶斯推断及其互联网应用（一）：定理简介2.贝叶斯推断及其互联网应用（二）：...
小强学AI之 - 2你患癌的概率（朴素贝叶斯）
前提请先认真阅读下面两篇博文：阮一峰博客（贝叶斯推断及其互联网应用（一）：定理简介）阮一峰博客（贝叶斯推断及其互...
【转】贝叶斯推断及其互联网应用---阮一峰
一、什么是贝叶斯推断贝叶斯推断（Bayesian inference）是一种统计学方法，用来估计统计量的某种性质...
贝叶斯分类器（1）贝叶斯决策论概述、贝叶斯和频率、概率和似然
贝叶斯分类器，即是以贝叶斯决策理论为基础的分类器，什么是贝叶斯决策理论呢？贝叶斯决策论 1 统计推断中的贝叶斯学...
贝叶斯分类器(10X单细胞和10X空间转录组的基础算法)
贝叶斯分类器，即是以贝叶斯决策理论为基础的分类器，什么是贝叶斯决策理论呢？贝叶斯决策论 1 统计推断中的贝叶斯学...
【机器学习实战（四）】sklearn包中朴素贝叶斯库的使用
目录朴素贝叶斯相关知识点回顾1.1. 什么是朴素贝叶斯分类器1.2. 朴素贝叶斯推断1.3. 朴素贝叶斯学习sk...
使用pymc3进行贝叶斯推断的简单示例
本文是学习《贝叶斯方法概率编程与贝叶斯推断》重点参考了这个教程[https://docs.pymc.io/en/...
频率学派与贝叶斯学派
序本次记录频率学派的主要思想与贝叶斯的区别，并给出贝叶斯公式频率学派与贝叶斯学派的差别频率派把需要推断的参数...
Python数据分析与机器学习26-贝叶斯算法
一. 贝叶斯算法介绍 1.1 贝叶斯简介贝叶斯(约1701-1761) Thomas Bayes，英国数学家贝...