频率主义和贝叶斯主义——一个实用的介绍

此notebook最初出现在博客Pythonic Perambulations的文章。BSD licensed。

这个系列共4个部分：中文版Part I Part II Part III Part IV，英文版Part I Part II Part III Part IV

科学工作者学习统计学的第一件事儿是要知道有两种不同的方法论：频率主义和贝叶斯主义。尽管这点很重要，但很多科学工作者从来没有机会认识两者的差异，及其产生结果的不同。本文的目的是综合频率主义和贝叶斯主义的理论与实践内容，让科学工作者在数据分析之前做出更好的选择。

开始我将通过理论来论述两者区别，之后用Python代码演示两者在应用中的差异。

频率主义 VS 贝叶斯主义：理论之争

频率主义与贝叶斯主义的本质之争是对概率定义的差异。

频率主义认为概率是表示有限次的重复测量 。假设我多次测量某个星星的亮度(光量子通量)(F) 时(假设星星通量不会变化)，由于测量误差每次得到不同的测量值。在多次测量的前提下，某个变量的频率指的是测量变量的概率。频率主义认为概率本质上是与频率相关的事件 。因此，一个纯粹的频率主义者会认为，讨论星星真实亮度的概率是没啥意义的，因为真实亮度其实是固定值，也就是说讨论一个固定值的频率分布是瞎扯淡。

贝叶斯理论认为，概率的概念应该被扩展为对一个命题信任的程度 。贝叶斯理论会说测量概率为

奇怪的是理论争议不大，但是在实际的数据分析中会产生极大差异。下面我会给出一个实际案例的两种理论实现，用Python演示。

频率理论与贝叶斯理论的实际应用: 光量子计量

让我们看一个简单问题，然后比较一下两种方法的差异。这儿得有点儿数学公式，不过我不会入得特别深哈。如果你想深入，你可以考虑一下练练教材4-5章的内容。

研究对象: 简易光量子计量方法

想象一下我们拿着望远镜，观察遥远的星星传来的无尽星光。首先，让我们假设星星的真实亮度是固定值，不会随时间变化。用常量

问题是有了这些测量值

(关于测量误差：如果假设测量误差是按正态分布。频率主义者的观点是

我们用Python随机生成一些数据来演示两种理论的差异。由于测量是按次计数的离散事件，我们可以使用泊松分布做演示：

In [1]:

# Generating some simple photon count data
import numpy as np
from scipy import stats
np.random.seed(1)  # for repeatability

F_true = 1000  # true flux, say number of photons measured in 1 second
N = 50 # number of measurements
F = stats.poisson(F_true).rvs(N)  # N measurements of the flux
e = np.sqrt(F)  # errors on Poisson counts estimated via square root

现在让我们为”测量“值做个简单的图形：

In [2]:

get_ipython().magic(u'matplotlib inline')
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.errorbar(F, np.arange(N), xerr=e, fmt='ok', ecolor='gray', alpha=0.5)
ax.vlines([F_true], 0, N, linewidth=5, alpha=0.2)
ax.set_xlabel("Flux");ax.set_ylabel("measurement number");

用标准平方根法则对泊松分布进行分析，可以看到每次测量值的误差都不同。上面的例子中我们是先知道真实值

下面让我们来看看两种理论的不同解决方案。

简易光量子计量的频率主义方法

我们用频率主义中经典的似然估计方法开始。对应某一次测量值

P (D i | F t r u e) = 1 2 π e 2 i----\sqrt exp [ - ( F i - F

公式应该读作“在

我们通过计算每个观测值概率的乘积来构建似然函数。

L (D | F t r u e) = \prod i = 1 N P (D i | F t r u e)

其中，

log L = - 1 2 \sum i = 1 N [ log ( 2 π e 2 i ) + ( F i - F t r u e

我们这么做是为了确定

F e s t = \sum w i F i \sum w i ; w i = 1 / e 2 i

特殊情况：当所有的测量误差

F e s t = 1 N \sum i = 1 N F i

也就是说，当误差相等时，

我们可以进一步研究下这么估计的误差。按照频率主义方法，可以通过李雅普诺夫定理将最大似然函数的曲线进行正态近似；在这个简单的例子中也可以容易解决。正态近似的标准差是：

σ e s t = (\sum i = 1 N w i) - 1 / 2

这个结果很容易计算；让我们估计一些我们的数据：

In [3]:

w = 1. / e ** 2
print("""
      F_true = {0}
      F_est  = {1:.0f} +/- {2:.0f} (based on {3} measurements)
      """.format(F_true, (w * F).sum() / w.sum(), w.sum() ** -0.5, N))

      F_true = 1000
      F_est  = 998 +/- 4 (based on 50 measurements)

可以看出，以50次测量得到的观测值为输入，得出我们的估计误差在千分之4。

简易光量子计量的贝叶斯方法

你可能认为，如果按照贝叶斯方法，通篇都得使用概率计算。我们真正想计算的其实是我们对参数的主观认识，在本例中就是：

P (F t r u e | D)

可以看出这与频率主义完成相反，要对

贝叶斯主义使用贝叶斯定理这一基本概率论公式：

P (F t r u e | D) = P ( D | F t r u e ) P ( F t r

虽然贝叶斯定理这一法则本身没有争议，贝叶斯主义的名字就从那儿来，但是贝叶斯主义对概率有争议的解释却是暗含在法则的

让我们仔细看看公式中的术语：

如果我们设置

P (F t r u e | D) \propto L (D | F t r u e)

也就是说贝叶斯主义的概率是频率主义似然估计结果的最大值！因此，遑论两种理论的差异，我们会发现(至少是简单问题)两种估计是等价的。

什么叫先验概率？

你会发现，我为了图省事儿，先验概率

当先验事实不存在时，频率主义者会说先验概率不靠谱。尽管用一个未知先验概率更简单，像前面提到那个先验概率常量，但是也会有一些意外的微妙之处。很多实验表明，真正的未知先验概率是不存在的！频率主义者认为主观选择先验概率会偏离结果，不能进行统计数据分析。

贝叶斯主义者认为频率主义解决不了问题，只想简单的避开。在确定（隐含的）先验概率这点上，频率主义经常被简单地看作贝叶斯主义理论的特殊情形：贝叶斯主义会要求最好让隐含变成明确，即使这种选择可能包含主观性。

光量子计量：贝叶斯方法

离开这些理论之争，让我们看看贝叶斯主义的实际计算结果。因为只考虑一个参数，可以将后验概率

这里不细究MCMC的理论知识，我用Dan Foreman-Mackey的Python emcee包做一个简单的应用。目标是产生一堆满足后验概率分布的点，用这些点来确定我们问题的答案。

要演示MCMC，我们用Python定义先验概率

In [4]:

def log_prior(theta):
    return 1  # flat prior

def log_likelihood(theta, F, e):
    return -0.5 * np.sum(np.log(2 * np.pi * e ** 2)
                         + (F - theta[0]) ** 2 / e ** 2)

def log_posterior(theta, F, e):
    return log_prior(theta) + log_likelihood(theta, F, e)

现在我们开始解决问题，产生随机数来猜测这些马尔可夫链上的点。

In [5]:

ndim = 1  # number of parameters in the model
nwalkers = 50  # number of MCMC walkers
nburn = 1000  # "burn-in" period to let chains stabilize
nsteps = 2000  # number of MCMC steps to take

# we'll start at random locations between 0 and 2000
starting_guesses = 2000 * np.random.rand(nwalkers, ndim)

import emcee
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F, e])
sampler.run_mcmc(starting_guesses, nsteps)

sample = sampler.chain  # shape = (nwalkers, nsteps, ndim)
sample = sampler.chain[:, nburn:, :].ravel()  # discard burn-in points

如果都能正常运行，数值sample会包含满足后验概率分布的50000个点。让我们画出来看看：

In [6]:

# plot a histogram of the sample
plt.hist(sample, bins=50, histtype="stepfilled", alpha=0.3, normed=True)

# plot a best-fit Gaussian
F_fit = np.linspace(975, 1025)
pdf = stats.norm(np.mean(sample), np.std(sample)).pdf(F_fit)

plt.plot(F_fit, pdf, '-k')
plt.xlabel("F"); plt.ylabel("P(F)")

Out[6]:

<matplotlib.text.Text at 0x7f061dec3f90>

我们从(正态)后验概率分布得到一样本。后验概率的期望和标准差是由前面的频率主义极大似然估计得来的：

In [7]:

print("""
      F_true = {0}
      F_est  = {1:.0f} +/- {2:.0f} (based on {3} measurements)
      """.format(F_true, np.mean(sample), np.std(sample), N))

      F_true = 1000
      F_est  = 998 +/- 4 (based on 50 measurements)

可以看出，在简单情况下，贝叶斯主义和频率主义两种方法的结果相同！

分析

现在，你可能认为贝叶斯主义的复杂性没啥必要，这个例子已经充分证明了。其实，用一个仿射不变MCMC综合样本发生器来模拟这样一个简单的一维正态分布，实在是杀鸡用牛刀，但是这么做的目的是要表明该方法同样适用于多维度的复杂后验概率情形，也能在复杂条件下提供似然估计不可能提供的优美的结果。

另外，你还可能注意到一个小伎俩：最后，我们用来频率主义方法来描述后验概率样本！当我们计算样本期望和标准差的时候，我们使用了一个标准的频率主义方法来刻画后验概率分布。用纯贝叶斯理论分析这个问题的结果可能要包含后验概率本身(比如，它的类似样本)。也就是说，在纯贝叶斯主义的结果中，问题不是一个简单的数值和误差，而是后验概率分布和模型参数！

增加一个维度：研究更复杂的模型

让我们看一个更复杂的情况，再比较一下频率主义和贝叶斯主义的结果。前面我们假设星星是静止的，现在假设我们观察的天体是一个随机变量；就是说它会随着时间变化(就是类星体那种)。

我们将用双参数正态分布来表示这个天体：

F t r u e \sim 1 2 π σ 2----\sqrt exp [ - ( F - μ ) 2 2 σ 2 ]

现在，我们再考虑

In [8]:

np.random.seed(42)  # for reproducibility
N = 100  # we'll use more samples for the more complicated model
mu_true, sigma_true = 1000, 15  # stochastic flux model

F_true = stats.norm(mu_true, sigma_true).rvs(N)  # (unknown) true flux
F = stats.poisson(F_true).rvs()  # observed flux: true flux plus Poisson errors.
e = np.sqrt(F)  # root-N error, as above

可变条件下的光量子计量：频率主义方法

似然估计结果是本征分布与误差分布的卷积，可以得出

L (D | θ) = \prod i = 1 N 1 2 π ( σ 2 + e 2 i )---------\sqrt exp [ - (

与前面类似，我们可以通过求微分解出似然估计的最大值，即

μ e s t = \sum w i F i \sum w i ; w i = 1 σ 2 + e 2 i

这里有个问题：参数

然而，我们能用数值分析计算来解极大似然估计值。这里我们用Scipy里面的最优化工具optimize：

In [9]:

def log_likelihood(theta, F, e):
    return -0.5 * np.sum(np.log(2 * np.pi * (theta[1] ** 2 + e ** 2))
                         + (F - theta[0]) ** 2 / (theta[1] ** 2 + e ** 2))

# maximize likelihood <--> minimize negative likelihood
def neg_log_likelihood(theta, F, e):
    return -log_likelihood(theta, F, e)

from scipy import optimize
theta_guess = [900, 5]
theta_est = optimize.fmin(neg_log_likelihood, theta_guess, args=(F, e))
print("""
      Maximum likelihood estimate for {0} data points:
          mu={theta[0]:.0f}, sigma={theta[1]:.0f}
      """.format(N, theta=theta_est))

Optimization terminated successfully.
         Current function value: 502.839505
         Iterations: 58
         Function evaluations: 114

      Maximum likelihood estimate for 100 data points:
          mu=999, sigma=19

Scipy给出了

在频率主义模式中有不同的方法来确定误差。如前所述，我们可以用正态近似去拟合最大似然估计和协方差矩阵(这里我们应该用数值分析而不是解析方法)。另外，我们可以计算统计量如

这些方法都行，但是每种方法都有各自的假设前提和特点。我们用Python的astroML包里面的basic bootstrap resampler:

In [10]:

from astroML.resample import bootstrap

def fit_samples(sample):
    # sample is an array of size [n_bootstraps, n_samples]
    # compute the maximum likelihood for each bootstrap.
    return np.array([optimize.fmin(neg_log_likelihood, theta_guess,
                                   args=(F, np.sqrt(F)), disp=0)
                     for F in sample])

samples = bootstrap(F, 1000, fit_samples)  # 1000 bootstrap resamplings

和前面的贝叶斯后验概率的MCMC方法类似，我们将计算样本期望值和标准差来确定参数误差。

In [11]:

mu_samp = samples[:, 0]
sig_samp = abs(samples[:, 1])

print " mu    = {0:.0f} +/- {1:.0f}".format(mu_samp.mean(), mu_samp.std())
print " sigma = {0:.0f} +/- {1:.0f}".format(sig_samp.mean(), sig_samp.std())

 mu    = 999 +/- 4
 sigma = 18 +/- 5

我得说明bootstrap resampling的细究有大堆资料可供参考，这里粗略地说明一下这个方法的特点。最明显的地方就是误差需要是相关的、非正态的，但是这两个特点都不能通过简单的计算每个模型参数的均值和方差反映出来。然而，我相信这给频率主义提供的基本的解题思路。

可变条件下的光量子计量：贝叶斯主义方法

贝叶斯主义对这个问题的解决方法可之前那个问题基本雷同，我们可以把上面的代码简单改改就行。

In [12]:

def log_prior(theta):
    # sigma needs to be positive.
    if theta[1] <= 0:
        return -np.inf
    else:
        return 0

def log_posterior(theta, F, e):
    return log_prior(theta) + log_likelihood(theta, F, e)

# same setup as above:
ndim, nwalkers = 2, 50
nsteps, nburn = 2000, 1000

starting_guesses = np.random.rand(nwalkers, ndim)
starting_guesses[:, 0] *= 2000  # start mu between 0 and 2000
starting_guesses[:, 1] *= 20    # start sigma between 0 and 20

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, args=[F, e])
sampler.run_mcmc(starting_guesses, nsteps)

sample = sampler.chain  # shape = (nwalkers, nsteps, ndim)
sample = sampler.chain[:, nburn:, :].reshape(-1, 2)

现在有了样本，我们可以利用astroML包方便的功能来画出轨迹图，并用等高线描述1、2标准差的情况：

In [13]:

from astroML.plotting import plot_mcmc
fig = plt.figure()
ax = plot_mcmc(sample.T, fig=fig, labels=[r'$\mu$', r'$\sigma$'], colors='k')
ax[0].plot(sample[:, 0], sample[:, 1], ',k', alpha=0.1)
ax[0].plot([mu_true], [sigma_true], 'o', color='red', ms=10);

红点显示在真实值附近(在前述问题条件下)，等高线显示了1、2个标准差的置信区间情况(68%和95%置信水平)。

这里要注意

另外，后验概率确实是非正态的：可以通过垂直方向的对称性来看。也就意味着用频率主义的正态近似不会将真实的不确定性反映到结果中。对频率主义者来说这不是啥事儿(频率主义模式有方法会考虑到非正态性)，但是频率主义方法的绝大多数使用者都直接或间接地使用了正态分布的假设。贝叶斯方法一般不需要这个假设。

(先验概率:有些观点在确定

结论

我希望我已经通过这篇文章展示了频率主义与贝叶斯主义根本的区别，导致两者本质上对同样问题有着不同的解释，然而彼此却经常可以产生相似甚至相同的结果。

两者差异总结如下：

频率主义认为概率与事实或假想事件的频率有关。
贝叶斯主义认为概率是对事件认识程度的度量。
频率主义通常使用点估计和极大似然估计方法进行分析。
贝叶斯主义通常直接的计算后验概率，或者通过一些MCMC抽样的方法计算。

处理简单问题时，两种方法会产生近似的结果。但当数据和模型越来越复杂时，两种方法会有天壤之别。在以后的文章中，我会显示一两个更复杂的案例。继续关注噢！

Update: 看下一篇文章：频率主义与贝叶斯主义 II: 当结果不同时

菜鸡一枚

frequentism-and-bayesianism-chs

frequentism-and-bayesianism-chs

频率主义和贝叶斯主义——一个实用的介绍

频率主义 VS 贝叶斯主义：理论之争

频率理论与贝叶斯理论的实际应用: 光量子计量

研究对象: 简易光量子计量方法

简易光量子计量的频率主义方法

简易光量子计量的贝叶斯方法

什么叫先验概率？

光量子计量：贝叶斯方法

分析

增加一个维度：研究更复杂的模型

可变条件下的光量子计量：频率主义方法

可变条件下的光量子计量：贝叶斯主义方法

结论

公告

菜鸡一枚

frequentism-and-bayesianism-chs

frequentism-and-bayesianism-chs

频率主义和贝叶斯主义——一个实用的介绍

频率主义 VS 贝叶斯主义： 理论之争

频率理论与贝叶斯理论的实际应用: 光量子计量

研究对象: 简易光量子计量方法

简易光量子计量的频率主义方法

简易光量子计量的贝叶斯方法

什么叫先验概率？

光量子计量：贝叶斯方法

分析

增加一个维度：研究更复杂的模型

可变条件下的光量子计量：频率主义方法

可变条件下的光量子计量：贝叶斯主义方法

结论

公告

频率主义 VS 贝叶斯主义：理论之争