频率主义 vs 贝叶斯主义 II：当结果不同时

这个notebook出自Pythonic Perambulations的博文 . The content is BSD licensed.

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

在上一篇我论述了频率主义和贝叶斯主义在科学计算方面的差异。其中，讨论了两者基础理论的差异，表明在简单问题的处理方面两者基本等价。

尽管如此，在处理复杂问题时两者会有巨大分歧。通过我的实践表明，这种分歧主要出现在两个方面：

冗余参数(nuisance parameters)的处理
频率论的置信区间(confidence intervals)和贝叶斯可信区间(credible intervals)的小差异(经常被忽略)。

第二点有点复杂，我打算后面再讲，这篇内容先讲第一点：冗余参数的差异。尽管我会尽力保持公正，但你还是会觉得我偏向贝叶斯方法。不过，这些都只是解决第二点之前的热身赛，因为那会引起更多口水战(get downright polemical)。

什么是冗余参数？

冗余参数是与最终目标无关的任何变量，但是对确定一些相关的变量却不可或缺。例如，在上一篇博文的第二个案例中，我们估计了光量子的均值

我打算用两个例子来论述冗余参数。下面我们看看频率主义和贝叶斯主义解决这些问题差异：

例 #1：贝叶斯台球(The Bayesian Billiard Game)

我打算从Thomas Bayes1763的论文里介绍的游戏出发来讨论冗余参数。这里我选一种改进版，详细介绍在Eddy 2004里面。

背景是两个人Alice和Bob赌输赢，但是他们看不到游戏的过程：

在一个房间里Alice和Bob前面挂着帘子，后面是个台球桌，他们看不到。但是他们的朋友Carol可以，Carol掷出一个球，然后记下球的位置。之后，他继续掷，如果这个球在第一个球的左边，Alice得一分。否则Bob得一分。这里假设Carol是公正无私的：也就是说他掷出的球呈均匀分布。先得6分者胜，也就是11局6胜。

第一个球的位置就是一个冗余参数：它是不确定的，也非结果，但是很明显后面的球都受其影响。如果这个球很靠右，那么Alice就占便宜。否则这个球很靠左，那么Bob就占便宜。

那么我的问题是：

在某次游戏中，当Alice赢5分时Bob才赢3分，那么Bob最后获胜的概率多大？

乍一看8个球里面Alice赢了5个球，标记位置可能更偏向她。那么，后面的球也更可能偏向她。而且她只要得一个球赢了，比Bob多3个机会；看来她是赢定了。但是定量分析Bob赢的概率是多少呢？

简单频率论方法

持频率方法论的人可能会如此解释：

要想获得结果，我们得先间接的估计一下标记的位置，假设概率

p^= 5 / 8

(结果服从二项分布的最大似然估计binomial likelihood)。假设Bob获胜的最大似然估计为

P (B) = (1 - p^) 3

也就是说，Bob要连续赢3次。我们可以算出结果：

In [1]:

p_hat = 5. / 8.
freq_prob = (1 - p_hat) ** 3
print(u"Bob获胜的频率论概率为: {0:.2f}".format(freq_prob))

Bob获胜的频率论概率为: 0.05

我们改为赔率(odds)的方式显示：

In [2]:

print(u"Bob获胜的赔率为: {0:.0f} to 1".format((1. - freq_prob) / freq_prob))

Bob获胜的赔率为: 18 to 1

可见，Bob每赢一次，Alice要赢18次。让我们再试试贝叶斯方法。

贝叶斯方法

我们再从贝叶斯的角度看看这个问题。有点小复杂，得先定义一些变量：

具体定义如下：

我们要计算

通常，贝叶斯方法处理冗余参数的方法是边缘概率法(marginalization)，即获得冗余参数全部可能的联合概率分布(joint probability)。也就是说，我们得先计算联合概率分布

P (B, p | D)

然后用下面的等式边缘化

P (B | D) \equiv \int \infty - \infty P (B, p | D) d p

这个等式源自条件概率(conditional probability)和总概率(total probability)定义，源自概率率公理的基本推论，不容置疑。频率主义者也认同这点，即使他们可能不同意我们把

建立贝叶斯表达式

要计算结果，我们就要把上面的

我们用条件概率(conditional probability)来展开

P (B | D) = \int P (B | p, D) P (p | D) d p

然后再用贝叶斯定理(Bayes' rule)改写

P (B | D) = \int P (B | p, D) P ( D | p ) P ( p ) P ( D ) d

最后，用前面同样的概率恒等式，可以将

P (B | D) = \int P ( B | p , D ) P ( D | p ) P ( p ) d p \int P (

现在的公式就可以计算了，我解释一下：

把这些放在一起，如下所示：

P (B | D) = \int 1 0 ( 1 - p ) 6 p 5 d p \int 1 0 ( 1 - p ) 3 p 5 d p

其中，积分范围介于0和1之间。

积分看着有点复杂，如果看出是贝塔函数(Beta Function)就好办了：

β (n, m) = \int 1 0 (1 - p) n - 1 p m - 1

贝塔函数可以进一步表示成伽马函数(Gamma function)，如阶乘(factorials)，但是我们可以直接用Scipy的beta函数简化计算：

In [3]:

from scipy.special import beta
bayes_prob = beta(6 + 1, 5 + 1) / beta(3 + 1, 5 + 1)

print("P(B|D) = {0:.2f}".format(bayes_prob))

P(B|D) = 0.09

我们改为赔率的方式显示：

In [4]:

print("Bob获胜的贝叶斯赔率为: {0:.0f} to 1".format((1. - bayes_prob) / bayes_prob))

Bob获胜的贝叶斯赔率为: 10 to 1

可以看出贝叶斯的结果是10比1，而频率论方法是18比1。究竟哪个是对的呢？

穷举/蒙特卡罗方法(Brute Force/Monte Carlo Approach)

由于本例内容太简单，实际上很容易通过蒙特卡罗方法仿真来获得结果。通过穷举获得所有的可能：随机产生大量的结果，然后计算Bob最后获胜的比例。本例很多相关的随机变量都是均匀分布的，可以通过 numpy 包就可以搞定：

In [5]:

import numpy as np
np.random.seed(0)

# play 100000 games with randomly-drawn p, between 0 and 1
# 玩1000万次游戏，p介于0到1之间
p = np.random.random(1000000)

# each game needs at most 11 rolls for one player to reach 6 wins
# 每次游戏最多需要投掷11次保证一个玩家赢6次，
# 最坏的情况是第10次5比5，那么第11次一定见分晓，
# 就像乒乓球比赛中的11球6胜
rolls = np.random.random((11, len(p)))

# count the cumulative wins for Alice and Bob at each roll
# 分别计算Alice和Bob赢的次数。
Alice_count = np.cumsum(rolls < p, 0)
Bob_count = np.cumsum(rolls >= p, 0)

# sanity check: total number of wins should equal number of rolls
# 合理性检验：总获胜次数应该等于投掷数
total_wins = Alice_count + Bob_count
assert np.all(total_wins.T == np.arange(1, 12))
print("(Sanity check passed)")

# determine number of games which meet our criterion of (A wins, B wins)=(5, 3)
# this means Bob's win count at eight rolls must equal 3
# 确定符合初始条件是(5, 3)的游戏次数
# 也就是说Bob在8次后获胜3次
good_games = Bob_count[7] == 3
print("Number of suitable games: {0}".format(good_games.sum()))

# truncate our results to consider only these games
# 只考虑满足初始条件的游戏
# Alice_count = Alice_count[:, good_games]
Bob_count = Bob_count[:, good_games]

# determine which of these games Bob won.
# to win, he must reach six wins after 11 rolls.
# 确定最终Bob获胜的的游戏次数
# 也就是11次后获胜6次
bob_won = np.sum(Bob_count[10] == 6)
print("Number of these games Bob won: {0}".format(bob_won.sum()))

# compute the probability
# 计算Bob获胜概率
mc_prob = bob_won.sum() * 1. / good_games.sum()
print("Monte Carlo Probability of Bob winning: {0:.2f}".format(mc_prob))
print("MC Odds against Bob winning: {0:.0f} to 1".format((1. - mc_prob) / mc_prob))

(Sanity check passed)
Number of suitable games: 111467
Number of these games Bob won: 10079
Monte Carlo Probability of Bob winning: 0.09
MC Odds against Bob winning: 10 to 1

蒙特卡罗仿真和贝叶斯的结果都是10比1。表面上看，简单频率论方法挂了。

讨论

本例比较了冗余参数p的不同处理方式。蒙特卡罗仿真给了一个真实可能的穷举方法(假设之前的假设都靠谱)，表明贝叶斯的结果是对的。而频率论方法用一个最大似然估计值，看来是错了。

随着我的评论和跟帖不断扩散，我在Reddit和Hacker News上被网友撕成了碎片，我的声明一点并不是说频率主义是错的。结果的错误源自“简单”而非“频率率”。当然，频率论中有一堆方法处理冗余参数——如通过数据变换和训练来消除对

另外，本例对经典的频率论来说可能不太公平。频率论者可能希望用假设检验或置信区间(null tests or confidence intervals)给出结果：也就是设计一个方法来构建边界，类似实验中正确的结果为百分之

这两条频率论观点有一个共同之处：都需要一定程度的努力和/或专业知识；可能一个靠谱的频率论方法对一个统计学博士是显而易见的，但是对一个普通人来说未必那么简单。总之，我认为贝叶斯主义是解决这类问题的良方：通过几条贝叶斯公理做简单的代数运算，就可以直接获得正确的结果。

下面我们研究一个复杂点的例子。

例 #2: 带奇点(Outliers)的线性拟合(Linear Fit)

冗余参数适合处理数据奇点。考虑下面的数据集合

In [6]:

x = np.array([ 0,  3,  9, 14, 15, 19, 20, 21, 30, 35,
              40, 41, 42, 43, 54, 56, 67, 69, 72, 88])
y = np.array([33, 68, 34, 34, 37, 71, 37, 44, 48, 49,
              53, 49, 50, 48, 56, 60, 61, 63, 44, 71])
e = np.array([ 3.6, 3.9, 2.6, 3.4, 3.8, 3.8, 2.2, 2.1, 2.3, 3.8,
               2.2, 2.8, 3.9, 3.1, 3.4, 2.6, 3.4, 3.7, 2.0, 3.5])

我们画个图：

In [7]:

%matplotlib inline
import matplotlib.pyplot as plt

plt.errorbar(x, y, e, fmt='.k', ecolor='gray');

我们的任务是找数据的最佳拟合线。很明显，数据里面有一些奇点，但是让我们用一个不稳定(non-robust)最大似然估计方法开始。如前所述，下面的最大似然估计结果要么概率论要么贝叶斯(有均匀分布前提)：在这类问题里面，这些方法基本等价。

我们先做一个简单的线性模型，有一个斜率和一个截距的参数矢量

y^(x | θ) = θ 0 + θ 1 x

用这个模型，我们可以计算正态似然估计：

p (x i, y i, e i | θ) \propto exp [- 1 2 e 2 i ( y i - y ^ (

总样本的似然估计是个体似然估计的乘积。两边取对数：

log L (D | θ) = c o n s t - \sum i 1 2 e 2 i ( y i - y ^ (

如果看过上一节的内容，会很明白。最终的式子是数据模型的对数似然估计(log-likelihood)，可以找到

l o s s = \sum i 1 2 e 2 i ( y i - y ^ ( x i | θ ) ) 2

这个表达式就是平方损失(squared loss)；我们简单的演示了平方损失源自正态对数似然估计。

标准似然估计法

按照前面的逻辑，我们用频率论中求最大似然估计(或者最小损失)的方式来算

简单起见，我们用scipy的optimize包来获得最小损失(如果用矩阵方法，平方损失的计算会更困难，但是我们使数值最小化，就比较容易处理)。

In [8]:

from scipy import optimize

def squared_loss(theta, x=x, y=y, e=e):
    dy = y - theta[0] - theta[1] * x
    return np.sum(0.5 * (dy / e) ** 2)

theta1 = optimize.fmin(squared_loss, [0, 0], disp=False)

xfit = np.linspace(0, 100)
plt.errorbar(x, y, e, fmt='.k', ecolor='gray')
plt.plot(xfit, theta1[0] + theta1[1] * xfit, '-k')
plt.title('Maximum Likelihood fit: Squared Loss');

由平方损失函数的本性，奇点对拟合线的影响不成比例。如果有一个奇点离拟合线有10个标准差的距离，那么它对损失函数的贡献比25个只有2个标准差远的奇点们重要。

很明显，平方损失过度依赖奇点，对我们的拟合任何造成很大影响。频率论模式解决这个问题的手段就是通过调整平方损失函数，使它更稳定。

奇点的频率论纠正方法：Huber损失(Huber Loss)

损失函数的可能看起来有无限多种，但是相对靠谱的选择就是Huber loss。Huber损失引入一个临界值(critical value)使得损失曲线从平方变为线性。让我们画个图来比较一下Huber损失和标准损失的临界值

In [9]:

t = np.linspace(-20, 20)

def huber_loss(t, c=3):
    return ((abs(t) < c) * 0.5 * t ** 2
            + (abs(t) >= c) * -c * (0.5 * c - abs(t)))

plt.plot(t, 0.5 * t ** 2, label="squared loss", lw=2)
for c in (10, 5, 3):
    plt.plot(t, huber_loss(t, c), label="Huber loss, c={0}".format(c), lw=2)
plt.ylabel('loss')
plt.xlabel('standard deviations')
plt.legend(loc='best', frameon=False);

Huber损失对模型中拟合很好的点作用是等同于标准损失的，但是会降低奇点的影响。如一个20标准差的点标准损失200，但是c=3时的Huber损失只有55。我们来看看Huber损失的最佳拟合结果与标准损失做个对比(灰色线表示)。

In [10]:

def total_huber_loss(theta, x=x, y=y, e=e, c=3):
    return huber_loss((y - theta[0] - theta[1] * x) / e, c).sum()

theta2 = optimize.fmin(total_huber_loss, [0, 0], disp=False)

plt.errorbar(x, y, e, fmt='.k', ecolor='gray')
plt.plot(xfit, theta1[0] + theta1[1] * xfit, color='lightgray')
plt.plot(xfit, theta2[0] + theta2[1] * xfit, color='black')
plt.title('Maximum Likelihood fit: Huber loss');

直观的看，会觉得好像Huber损失拟合是我们要的结果。

但一个贝叶斯主义者可能会说这种新的损失函数的动机有点可疑：如前所述，平方损失源自正态似然估计，但Huber损失看着有点(随意)ad hoc，这是哪儿来的?

贝叶斯方法处理奇点：冗余参数

贝叶斯方法处理奇点一般通过调整模型(modifying the model) ，这样奇点就被考虑到。很明显，线性函数对这个模型不是一个好结果。因此，让我们设计一个更复杂的模型来考虑奇点。一种方式就是把信号和噪声(background)混合在一起：

p ({x i}, {y i}, {e i} | θ, {g i}, σ, σ b) =

我们要做的就是把模型加上冗余参数：

现在我们的模型足够复杂了：有22个参数，但是大多数可作为冗余参数，最终将被边缘化(marginalized-out)。恰似之前我们在台球游戏里面边缘化

为了算出结果，我们用函数开始定义先验概率，似然估计和后验概率：

In [11]:

# theta will be an array of length 2 + N, where N is the number of points
# theta[0] is the intercept, theta[1] is the slope,
# and theta[2 + i] is the weight g_i
# theta是长度2 + N的数组，其中N是点的数量
# theta[0]是截距，theta[1]是斜率
# theta[2 + i]是权重g_i

def log_prior(theta):
    #g_i needs to be between 0 and 1
    if (all(theta[2:] > 0) and all(theta[2:] < 1)):
        return 0
    else:
        return -np.inf  # recall log(0) = -inf

def log_likelihood(theta, x, y, e, sigma_B):
    dy = y - theta[0] - theta[1] * x
    g = np.clip(theta[2:], 0, 1)  # g<0 or g>1 leads to NaNs in logarithm
    logL1 = np.log(g) - 0.5 * np.log(2 * np.pi * e ** 2) - 0.5 * (dy / e) ** 2
    logL2 = np.log(1 - g) - 0.5 * np.log(2 * np.pi * sigma_B ** 2) - 0.5 * (dy / sigma_B) ** 2
    return np.sum(np.logaddexp(logL1, logL2))

def log_posterior(theta, x, y, e, sigma_B):
    return log_prior(theta) + log_likelihood(theta, x, y, e, sigma_B)

现在我们运行MCMC样本来探索参数空间：

In [12]:

# Note that this step will take a few minutes to run!

ndim = 2 + len(x)  # number of parameters in the model
nwalkers = 50  # number of MCMC walkers
nburn = 10000  # "burn-in" period to let chains stabilize
nsteps = 15000  # number of MCMC steps to take

# set theta near the maximum likelihood, with 
np.random.seed(0)
starting_guesses = np.zeros((nwalkers, ndim))
starting_guesses[:, :2] = np.random.normal(theta1, 1, (nwalkers, 2))
starting_guesses[:, 2:] = np.random.normal(0.5, 0.1, (nwalkers, ndim - 2))

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

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

-c:18: RuntimeWarning: divide by zero encountered in log
-c:19: RuntimeWarning: divide by zero encountered in log

一旦有了这些样本，我们就可以获得马尔可夫链(Markov chains)的一个很好的属性。因为这些样本的分布模型是后验的，所以我们可以通过忽略它们来整合出(如，边缘化)冗余参数。

我们通过样本的前两列来看斜率和截距的(边缘化)分布：

In [13]:

plt.plot(sample[:, 0], sample[:, 1], ',k', alpha=0.1)
plt.xlabel('intercept')
plt.ylabel('slope');

可见，斜率分布在

通过分析MCMC样本会发现，冗余参数的选择是完全对称的(completely symmetric)：就像我们可以把

In [14]:

plt.plot(sample[:, 2], sample[:, 3], ',k', alpha=0.1)
plt.xlabel('$g_1$')
plt.ylabel('$g_2$')

print("g1 mean: {0:.2f}".format(sample[:, 2].mean()))
print("g2 mean: {0:.2f}".format(sample[:, 3].mean()))

g1 mean: 0.64
g2 mean: 0.41

这两个参数不是很相关，但是我们可以看到

让我们来确认一下，画出原始数据的边缘化测试模型，奇点用红圈标识。

In [15]:

theta3 = np.mean(sample[:, :2], 0)
g = np.mean(sample[:, 2:], 0)
outliers = (g < 0.5)

plt.errorbar(x, y, e, fmt='.k', ecolor='gray')
plt.plot(xfit, theta1[0] + theta1[1] * xfit, color='lightgray')
plt.plot(xfit, theta2[0] + theta2[1] * xfit, color='lightgray')
plt.plot(xfit, theta3[0] + theta3[1] * xfit, color='black')
plt.plot(x[outliers], y[outliers], 'ro', ms=20, mfc='none', mec='red')
plt.title('Maximum Likelihood fit: Bayesian Marginalization');

看！上面黑线显示的结果更符合我们的直觉。而且，自动被识别的奇点也是我们可以手工识别出来的。为了比较，我把之前的方法用灰线标识：简单的最大似然估计和Huber损失的频率论方法。

讨论

今天，我们讨论线性回归(linear regression)中的奇点发现问题。一个经典的正态最大似然估计方法不能解决，但是在频率论中我们能够调整损失函数来纠正，在贝叶斯理论中通过引入一堆冗余参数的混合模型来调整。

两种方法都有优缺点：频率论方法相对直接而且计算高效，但基于难以解释的损失函数。贝叶斯方法容易建立并能产生好的结果，但是要求一堆主观前提。在写代码和计算时间上耗费更多。

最终结论

在上一篇中，我讨论了两种理论的差异，并且说明在简单问题下，忽略彼此差异，二者等价。在这一篇中，我们探讨了差异的一种：冗余参数处理。

在贝叶斯台球的例子中，我们论述频率论方法不靠谱，而贝叶斯方法靠谱。这并非说频率主义是错的，但是告诉我们用它的时候得当心。

在线性回归的例子中，我们论述了一种用频率主义和贝叶斯主义解决数据奇点问题的可行方法。用稳定的频率论损失函数相对的短平快，但是其动机混乱而难以解释。用贝叶斯混合模型更费力且需要大量计算，但是能够算出完美的结果，使多个问题被一次解决：也就是，边缘化一种方式找出最佳拟合模型，边缘化另种方式找出数据奇点。

哪个更好呢？频率主义和贝叶斯主义。

答案依赖于你对两者的掌握程度，也有赖于问题的复杂度，以及可用的计算资源。我，和很多物理背景的人一样，倾向于贝叶斯方法，因为它们能够从任何基本原理派生出我想要的结果。贝叶斯主义，基于一些概率公理的代数计算，我们可以构建非常复杂的方法来解决一大堆问题。就如质-能转换可以用到任何场合，从抛物运动(projectile motion)到恒星结构(stellar structure)，贝叶斯原理和贝叶斯概率解释能用来解决几乎任何统计问题：从计算赌博赔率(gambling odds)到探索噪声数据中的外星信息(exoplanet transits in noisy photometric data)。在贝叶斯模式中，你不需要花大量时间来记忆和理解频率论技巧和术语(p值，假设检验，置信区间，断点等等)，就为了解决一个小问题。我认为这是一种误解：人们总觉得贝叶斯很难。相反，很多科学家觉得它更容易。

简单和优美之外，我坚定贝叶斯主义立场的深层理由是必须比较科学研究中频率论置信区间(confidence intervals)和贝叶斯可信区域(credibility regions)。这篇已经够长了，我打算下一篇(http://jakevdp.github.io/blog/2014/06/12/frequentism-and-bayesianism-3-confidence-credibility/)来讨论。(有费马之嫌~~)

这篇博客是IPython notebook写的。可以[下载](http://jakevdp.github.io/downloads/notebooks/FreqBayes2.ipynb)，或者看[nbviewer](http://nbviewer.ipython.org/url/jakevdp.github.io/downloads/notebooks/FreqBayes2.ipynb)静态网页。

菜鸡一枚

frequentism-and-bayesianism-chs-ii