以 Python 为实例，介绍贝叶斯理论，python贝叶斯,未经许可，禁止转载！英文

本文由 编橙之家 - 巽离 翻译，蒋生武 校稿。未经许可，禁止转载！
英文出处：nbviewer.ipython.org。欢迎加入翻译组。

什么是贝叶斯理论？

贝叶斯理论就是允许我们根据抽样分布（或称似然度）和先验分布推导出后验分布。

什么是抽样分布？

抽样分布指的是在给定参数θ的条件下，X发生的概率。记做p(X|θ).

举例说明，我们记录了1000次抛硬币的数据。如果1表示硬币正面。可以用python这么表示：

输入 [2]:

Python

import numpy as np
data_coin_flips = np.random.randint(2, size=1000)
np.mean(data_coin_flips)

输出[2]:

0.46800000000000003

抽样分布允许我们指定数据的生成方式。对抛硬币来说，我们可以假设总体数据服从以 p 为参数的伯努利分布。这个参数 p 是抛出1的概率（抛硬币得到正面）。返回结果为1的概率p以及结果为0的概率（1-p）。

可以看出这种分布完美的切合抛硬币事件。如果硬币是均匀的情况下，我们知道p值为0.5，因为抛得1（正面）或0（反面）的机会是平均的。在python中我们可以通过下面两条语句创建这样的分布：

输入[3]:

Python

bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000)
np.mean(bernoulli_flips)

输出[3]:

0.46500000000000002

既然已经定义了我们所认为的数据生成方式，就可以根据给定参数来计算数据的似然概率P(X|θ)。 因为我们已经选定了伯努利分布，只需要一个参数p就可以了。

我们可以使用伯努利分布的概率质量函数（PMF）来计算单次抛硬币的期望概率。PMF根据一个数据值和给定的参数（本例中为p），返回相应的概率。对伯努利分布来说，这个pmf函数非常简单： 如果数据值为1， 那么它的概率返回p， 如果数据值为0， 那么返回概率为(1-p)。 我们可以定义如下函数：

In [4]:

Python

def bern_pmf(x, p):    
    if (x == 1):       
        return p    
    elif (x == 0):        
        return 1 - p   
    else:        
        return "Value Not in Support of Distribution"

现在我们使用这个函数，根据输入输参数来获取每个数据点的概率。如果参数p的值为0.5，那么返回值也总是0.5.

In [5]:

Python

print(bern_pmf(1, .5))
print(bern_pmf(0, .5))

0.5

0.5

这是一种非常简单的 PMF，不过其他分布的 PMF 会复杂得多。所以呢，我们应该了解这个模块 Scipy，它将这些都作为内建函数，直接调用就可以了。像下面这样描述PMF：

In [6]:

Python

import scipy.stats as st
print(st.bernoulli.pmf(1, .5))
print(st.bernoulli.pmf(0, .5))

0.5

0.5

这个模块很给力，但是我们实际上想知道的是这1000个样本数据的总体概率。怎么才能知道呢？这里有个小技巧，假设我们的样本数据是独立同分布的。那么我们就可以简单的认为所有数据的概率就是每个个体的独立概率的乘积：

p(x1,…,xn|β)=p(x1|β)∗…∗p(xn|β)。Python中可以这样实现：

In [7]:

Python

np.product(st.bernoulli.pmf(data_coin_flips, .5))

Out[7]:

9.3326361850321888e-302

得出这样的数据能做什么用呢？ 说实话，单个数据本身其实并没有什么卵用。我们现在需要做的是给我们的抽样模型生成更多的分布。现在，我们仅仅使用概率0.5来测试我们的模型，如果概率为0.8呢？ 或者0.2呢？我们数据的概率又会怎样？想要验证这些情况， 可以定义一个概率p 的网格。下面我将构造100个0到1之间的数据组成的网格（因为概率值都在0和1之间），之后呢，我会在这些值的基础上来计算我们抽样数据的概率。

In [8]:

Python

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inlines
ns.set(style='ticks', palette='Set2')params = np.linspace(0, 1, 100)
p_x = [np.product(st.bernoulli.pmf(data_coin_flips, p)) for p in params]
plt.plot(params, p_x)
sns.despine()

现在可以看出一点端倪了。 从图中很明显看出，当概率p=0.5的时候，数据分布的概率达到峰值，概率0.4到0.6之间也还算靠谱。Nice。假设我们的数据服从伯努利分布的话，我们可以很清楚的看出数据服从的参数 p 的值了。

先验分布

贝叶斯理论认为我们不仅要考虑抽样分布还得考虑先验分布。什么是先验分布呢？就是p(θ)，或者说为每个参数θ指定一个明确的概率。在我们的抽样分布中，我们给参数p定义了100个0到1之间的数据。现在我们需要对每个值都设定一个先验概率。这个概率是得到新数据前我们假设的概率。大多数情况下， 我们假设硬币是均匀的。看起来就像上图那样的分布。我们可以这样做：

In [9]:

Python

fair_flips = bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000)
p_fair = np.array([np.product(st.bernoulli.pmf(fair_flips, p)) for p in params])
p_fair = p_fair / np.sum(p_fair)
plt.plot(params, p_fair)
sns.despine()

基本上，我们就像之前一样创建1000次抛硬币数据并生成抽样分布（不过我们将数值除以抽样分布的总和使得所有值加起来等于1）.这样我们就得到了一个“均匀分布的硬币”的先验概率。这也就意味着在得到任何实际数据之前，我们主观上认定硬币是均匀的。我们也可以从上图中直接看出先验概率，0.5时达到顶峰，并且绝大多数处于0.4和0.6区间之内。

我知道你在想什么——你觉得这相当无聊对吧？抽样分布和先验分布看起来也没差啊。所以呢， 我们将它们组合在一起，我们仍然假设先验概率是均匀的，而将抽样数据换成非均匀的硬币：

In [27]:

Python

unfair_flips = bernoulli_flips = np.random.binomial(n=1, p=.8, size=1000)
p_unfair = np.array([np.product(st.bernoulli.pmf(unfair_flips, p)) for p in params])
fig, axes = plt.subplots(2, 1, sharex=True)
axes[0].plot(params, p_unfair)
axes[0].set_title("Sampling Distribution")
axes[1].plot(params, p_fair)
axes[1].set_title("Prior Distribution")
sns.despine()plt.tight_layout()

哎呀，事情变得有趣了呢。我们的数据证明硬币不是均匀分布的（因为我们根据非均匀的概率0.8来生成的数据）， 但是我们先验概率由坚称硬币是均匀的。现在我们该怎么办呢？

Bayes Theorem (Posterior Distribution)

贝叶斯理论允许我们从抽样和先验概率来推导出后验概率。后验概率 P(θ|x)，通俗一点地说，就是在给定数据情况下，我们要计算的该假设的概率。仔细想一想就会知道，我们一直想要的就是这个后验概率。这个给定的数据可能来自于问卷调查，也可能来自于互联网，在这些给定数据下，我们需要找出什么样的假设参数最有可能生成我们的样本数据。怎样得到后验分布呢？首先我们需要一些数学理论的说明（别担心，很简单）。

根据定义，我们可知（如果你信不过我的数学水平， 这个页面可以帮你复习）：

• P(A|B)=P(A,B)/P(B) 。通俗一点说，给定数据B的前提下A发生的概率等于他们的联合概率除以B的概率。
• P(B|A)=P(A,B)/P(A) 。与之类似的， 给定数据A的前提下B发生的概率等于他们的联合概率除以A的概率。

可能你也注意到了， 这俩公式有着相同的因子，所以呢：

• P(A,B)=P(A|B)∗P(B)
• P(A,B)=P(B|A)∗P(A)

因此：

P(A|B)∗P(B)=P(B|A)∗P(A)

即：

P(A|B)=P(B|A)∗P(A)/P(B)

这里就很明显啦，将θ和X带入公式，替换A、B，可知：

P(θ|X)=P(X|θ)∗P(θ)/P(X)

完美！现在们就可以带入我们知道的术语啦：

后验概率 = 似然度 * 先验概率P(X)

不过问题是：P(X)是个啥？ 换句话说，我们数据的概率是多少？有点诡异……

好吧，我们回到数学，继续使用B和A的公式：

我们知道P(B)=∑AP(A,B)  (如果忘了，这个页面可以帮你复习 )

从刚才的讨论，我们知道：

P(A,B)=P(A|B)∗P(A)

因此：

P(B)=∑AP(A|B)∗P(A)

带入 θ 和 X ：

P(X)=∑θP(θ|X)∗P(θ)

带入术语：

P(X)=∑θ*似然度 * 先验概率

简直不能更赞！但是 ∑θ 是啥意思呢？ 这是指我们所有参数值的总和。在我们抛硬币的例子中， 我们给参数 p 定义了一百个数值，所以我们需要对所有值计算它的似然度*先验概率并求和。这个值也就是贝叶斯理论的分母。因此，我们得出最后的公式：

后验概率 = 似然度 * 先验概率 / ∑θ * 似然度 * 先验概率

啰嗦了一大堆。上代码：

In [29]:

Python

def bern_post(n_params=100, n_sample=100, true_p=.8, prior_p=.5, n_prior=100):    
    params = np.linspace(0, 1, n_params)    
    sample = np.random.binomial(n=1, p=true_p, size=n_sample)    
    likelihood = np.array([np.product(st.bernoulli.pmf(sample, p)) for p in params])    
#likelihood = likelihood / np.sum(likelihood)    
    prior_sample = np.random.binomial(n=1, p=prior_p, size=n_prior)    
    prior = np.array([np.product(st.bernoulli.pmf(prior_sample, p)) for p in params])    
    prior = prior / np.sum(prior)    
    posterior = [prior[i] * likelihood[i] for i in range(prior.shape[0])]    
    posterior = posterior / np.sum(posterior)        
    fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8))    
    axes[0].plot(params, likelihood)    
    axes[0].set_title("Sampling Distribution")    
    axes[1].plot(params, prior)    
    axes[1].set_title("Prior Distribution")    
    axes[2].plot(params, posterior)    
    axes[2].set_title("Posterior Distribution")    
    sns.despine()    
    plt.tight_layout()        
    return posterior

In [25]:

example_post = bern_post()

你可能注意到了，我为先验概率和似然度设置了100个观察数据。增加了分布的方差。数据越多，分布数据表现的越集中。同样，当你有更多的数据来估计似然度时，先验分布的作用就会变小了。

In [30]:

moredata_post = bern_post(n_sample=1000)

从上图也可以看出数据量大小对结果的影响。因为我们使用更多的数据来估计似然度，从而我们的后验分布概率看起来与似然度十分接近。非常酷。

Conclusion

到这里，你对贝叶斯理论就有了一个大概的了解了。如果你曾经怀疑一枚硬币是否均匀，现在你就知道怎么来分析问题了。又或者你想估计投票结果的概率？或者其他是或否的问题。我知道你在想什么——你觉得这些都不够酷。你想要做出一些预言或者跑一个非常牛叉的算法…… 未来有可能。因为我确实想要写一个系列来介绍贝叶斯线性回归。希望这篇文章是这个系列的第一篇 :)

一些附注

1.    你可能注意到贝叶斯理论的分母是个常数。所以如果你想要获得后验概率最大化，不需要计算出这个常数的实际值。因为这个缘故，你常常可以看到后验概率正比例于（数学符号 ∝  ）似然度 * 先验概率。

2.    频率论统计关心的是似然度。或者可以说频率论就是无信息先验的贝叶斯（类似于均匀分布）。但也不要太讨厌频率论； 因为在实际应用过程中许多贝叶斯理论依赖于频率论统计。

3.   现在你知道频率论统计关注点在似然度，那么也就明白了为什么人们经常会误解频率论的置信区间。似然度是 P(X|θ) ——或者说在假设前提下得到我们数据的概率。有那么一点迷惑人，因为我们给定的是数据，而不是我们的假设。大多数频率论模型都是取似然函数分布的最大值（也称为最大似然估计）。主要是找出什么样的参数使得频率达到极值。这里有个关键点是：他们将数据视为随机的，所以频率论的置信区间意味着：如果你不断地有新的数据（可能是一些问卷调查）加入样本中，并对每个新样本来计算置信区间，那么其中 95% 的置信区间都会包含你正在估计的总体参数。