标签归档:统计学

A/B测试的数学原理与深入理解(达观数据 陈运文)

Deep Learning Specialization on Coursera

当面对众多选择时,如何选才能最大化收益(或者说最小化我们的开销)?比如,怎么选择最优的上班的路线才能使途中花费的时间最少?假设每天上下班路线是确定的,我们便可以在账本中记下往返路线的长度。

A/B测试便是基于数据来进行优选的常用方法,在记录多次上班路线长度后,我们便会从数据中发现到一些模式(例如路线A比路线B花的时间更少),然后最终一致选择某条路线。

当A/B测试遇到非简单情况时(如分组不够随机时,或用户量不够大到可以忽略组间差异,或不希望大规模A/B测试长期影响一部分用户的收益),该怎样通过掌握理论知识来更好的指导实践呢?本文尝试通过由浅入深的介绍,希望能够帮助大家对A/B测试有更加深入的理解。

NO.1
为什么需要A/B测试

任何问题,只要它的每个选项能够被多次进行测试,并且每个选项在被测试时都能返回固定的结果,那么它就能使用A/B测试技术来进行优化。在上述例子中,每天的上下班路线是确定的,所以我们能够在账本中记下往返路线的长度。

那么什么样的路线对于用户来说才是一个好的方案呢?是考虑路线A还是B?什么时候用户才有充分的数据去确定哪条线路是最好的?测试线路好与不好的最优策略又是什么?图1用形式化概括定义了问题。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

图1 形式化定义的问题

在这个场景中,参与的用户正面临一个选择,根据他的决策会生成一个结果,而这个结果会对应一份给参与者的反馈。假设用户持续地暴露于这个决策,他应该怎么制定获得最大收益(或等效地说,最小成本)的策略?

图1中假定了用户多次处于需要进行选择的场景中,每一次进行决策都会达成一项结果,而这个结果会关联相应的反馈。在上下班这个例子中,假定他每天都需要上下班,而且他每次上下班都必须进行线路的选择,产出的结果是这次上下班中所有因素的结合体,反馈就是从这些因素中构建出来的(陈运文 达观数据)。

这是个浅显的例子,在互联网产品研发时,有大量类似的场景需要做出各种正确的选择,例如:

1

着陆页优化(Landing-page optimization)

在用户点击去往的页面(着陆页),如何获得最大的转化率(常用计算方法为有购买行为或深度网页交互行为的用户数占网站访问总用户数的比率)。决策要考虑到着陆页的形式和内容(要从可能已有的3或4个备选方案中做出选择),希望能够从候选集合中选出最好的着陆页,以能够吸引来访的用户,并让深度交互或者购买行为的概率最大化。

2

广告创意优化(Ad creative optimization)

在线广告提出了许多适合机器学习技术应用的挑战,其中之一就是如何选择广告的形式和内容。当我们决定将要进行广告展示,以及确定了广告的价格后,在这个广告位上选择放置什么广告呢?我们需要对大量的决策进行测试,选出正确的广告创意组合。

NO.2
什么是A/B测试

经常遇到的问题是,我们应该怎么评估各不相同的决策,以及应该采用哪些策略来测试我们的产出? A/B测试(A/B testing)就是其中之一的方法。A/B测试近年来很受欢迎,但大部分产品经理也许会简单地认为它只不过是一种包含两个组的实验,其实背后有更为复杂的数学统计理论知识。

具体细节
当进行A/B测试时,通常会采用两个(或多个)组:A组和B组。第一个组是对照组,第二个组会改变其中一些因素。就以着陆页优化为例,A组会展示现有的着陆页,B组会展示一个内容或者内容作了某些修改的新着陆页。A/B测试的目的就是尝试了解新的布局是否在统计上显著地改变了转化率。

特别值得注意的是,将用户分配到对应的组需要经过深思熟虑。对于A/B测试,我们可以高效地进行随机分组。当用户数量较大时,各组间用户行为可以假设是相同的(即组间没有偏差)。但是,这里有三个非常重要的关键点,是大家有必要进一步理解其数学理论原理的原因:

1
问题1
怎样验证两个组的用户的行为是无偏差、完全相同的
2
问题2

当两个组的用户行为不完全相同时(例如分组不够随机或者组内用户数量较小时),该如何设计AB测试以实现期望的验证结果

3
问题3

当用户基础行为受其他因素影响发生整体变化了呢?例如季节、时间波动、热度等因素影响下,怎样更好的剔除干扰来评估结果

NO.3
AB测试的统计理论

假设我们已经构建了两组数目较大的用户组,这些用户组的区别仅在于他们到达的着陆页。我们现在希望能测试两组间的转化率在统计上是否存在明显差异。由于样本量大,我们可以采用双样本单尾z-检验(two-sample, one-tailed z-test)。另外,对于较小的样本集合,我们可以依赖于t-检验。

z检验(z-test)是在数据是正态分布和随机抽样的假设下运行的,目的是验证测试集(B组)是否与该对照集(A组)有显著不同,但是如何执行这个测试呢?

假设有来自A组和B组中的每一组的5,000个样本。我们需要一个数学公式来说明我们的零假设(null hypothesis)——两组群体的转化率没有显著的正差异,和备择假设(或称对立假设,alternative hypothesis)——不同人群间的转化率确实存在着正差异。

我们可将采样转化率视为一个正态分布的随机变量,也就是说,采样的转化率是在正态分布下对转化率的一个观测。要了解这一点,请考虑从同一组中提取多个样本进行实验将导致略有不同的转化率。每当对某组进行抽样时,可获得群体转化率的估计,对于A组和B组都是如此。为此我们提出一个新的正态随机变量,它是A和B组的随机变量的组合,是差值的分布。让我们用X来表示这个新的随机变量,定义为:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

其中,Xe表示实验组的转化率的随机变量,Xn表示对照组的转化率的随机变量。现在我们可以写出零假设和备择假设。零假设可以表示为:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

这表示实验组和对照组是相同的。两个随机变量Xe和Xn分布在相同的群体平均值周围,所以我们的新随机变量X应该分布在0左右。我们的备择假设可以表示如下:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

实验组的随机变量的期望值大于对照组的期望值;该群体的平均值较高。

我们可以在零假设的前提下,对X的分布执行单尾z检验,以确定是否有证据支持备择假设。为了达到这个目的,我们对X进行采样,计算标准分,并测试已知的显著性水平。

X的采样等效于运行两个实验,确定它们各自的转化率,并将对照组和实验组的转化率相减。按照标准分的定义,可以写作:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

其中,P_experiment是实验组的转化率,P_control 是对照组的转化率,SE是转化率差值的标准差。

为确定标准误差,注意到转化过程是符合二项分布的,因此访问该网站可以被看作单次伯努利试验(single Bernoulli trial),而积极结果(完成转化)的可能性是未知的。

假设样本数量足够大,我们可以使用广泛采用的Wald方法(参考Lawrence D. Brown, T. Tony Cai, and Anirban DasGupta, “Confidence Intervals for a Binomial Proportion and Asymptotic Expansions,” The Annals of Statistics 30, no. 1 (2002): 160–201.)将该分布近似为正态分布。为了捕获特定转化率的不确定性,我们可以将标准误差(SE)写入实验组和对照组,其中p是转化的可能性,n是样本数量,具体如下:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

从二项分布(np(1-p))的方差得到分子,而分母表示当采用更多的样本时,转化率的误差会随之下降。请注意正面结果的概率等同于转化率,并且因为两个变量的标准误差可以通过相加来合并,得到如下结果:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

通过替换,可获得如下的z检验公式,这是一个符合二项分布的Wald(或正态)区间的公式:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

z的值越大,反对零假设的证据就越多。为了获得单尾测试的90%置信区间,我们的z值将需要大于1.28。这实际上这是指在零假设(A组和B组的人口平均值是相同的)的条件下,等于或大于这个转化率差值的偶然发生的概率小于10%。

换句话说,在对照组和实验组的转化率来自具有相同平均值的分布的假设前提下,如果运行相同的实验100次,只会有10次具有这样的极端值。我们可以通过95%的置信区间,更严格的边界和更多的证据来反对零假设,这时需要将z值增加到1.65。

研究影响z大小的因素会带来很多有用的帮助。很显然,如果在一个给定的时间点从一个实验集和一个对照集中提取两个转化率,转化率的差值越大将导致z分数越大。因此就有了更多的证据表明两个集合分别来自不同的人群,而且这些人群带有不同的均值。然而样品的数量也很重要,如你所见,大量样本将导致总体较小的标准误差。这表明运行实验的时间越长,转化率的估算越准确。

NO.4
评估效果的代码实现

设想你在负责大型零售网站,设计团队刚刚修改了着陆页。每周有约20,000用户,并可以量化用户的转化率:即购买产品的百分比。设计团队向你保证新网站将带来更多的客户。但你不太确定,希望运行A / B测试来看看效果是否真的会提高。

用户在第一次访问网站时被随机分配到A组或B组,并在实验期间始终保留在该组中,实验结束时评估两组用户的平均转化率。统计结果是,新着陆页的平均转化率是0.002,而原先的着陆页的平均转化率是0.001。在着陆页永久更改为新设计之前,你需要知道这一增长是否足够明确。下面这段代码帮你回答这个问题。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

这段代码获取实验中z的值,在上述参数条件下z值为1.827,超过了92%置信区间,但不在95%的区间内。可以说,从控制分布中抽取数据的概率小于0.08。因此在该区间内数据提升是显著的。我们应该否定零假设,接受备择假设,即组之间有差异,第二组具有较高的转化率。如果我们控制了用户组的所有其他方面,就意味着网站的新设计产生了积极的效果。

你应该能够从代码中看到转化率分布的标准误差对返回的z值有直接影响。 对给定的常数值p_experiment和p_control,两个组的SE越高,z的数值越小,结果就越不显著。还注意到由于SE的定义,z的数值与样本的数量具有直接关系,对于给定的转换概率也同样如此。图2展示了这种关系。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

图2

图2 展示了A / B组的固定转化率,以及A / B组中的用户数量和z值之间的关系。 假设转化率不会随着我们收集更多数据而改变,我们需要每个组中大约3,000个用户达到70%的置信区间。 要达到80%的置信区间时需要每组约5000个用户,达到90%时需要 7500个用户,达到95%时需要12000个用户。

图2中可见对于两个组的给定转化率,测试组中的用户越多,备择假设的证据就越充分。直观上来看这很容易理解:当收集的数据越多,我们对结果越自信!我们也可以绘制一张类似的图,保持用户数量不变,改变组之间的差异。但必须注意,对正在关注的应用,不应该期望效果的大幅度变化。

 

NO.5
A/B测试方法的副作用和处理办法

对于非常小的效果变化,往往都需要创建相当大的对照组和测试组来实现AB测试,这个的代价往往是很大的。设想下在零售商场中,每天观察到的用户数量,往往需要很久的时间才能得出明显的结论。在实际业务应用中,会遇到的问题是:当你运行测试时整体运行的效果是受到很大影响的,因为必须有一半的用户处于效果不佳的实验组,或者有一半的用户处于效果不佳的对照组,而且你必须等待测试完成才能停止这种局面。

这是被称为探索利用难题(explore-exploit conundrum)的一个经典问题。我们需要运行次优方法,以探索空间,并找到效果更好的解决方案,而一旦找到了更好的解决方案,我们还需要尽快利用它们来实现效果提升。能否可以更快地利用新的解决方案,而不必等待测试完全完成呢?答案是肯定的。下面简单介绍下多臂赌博机(multi-armed bandit,MAB)的概念。

1

多臂赌博机的定义

多臂赌博机(multi-armed bandit,MAB)的名字来源于著名的赌博游戏角子赌博机(one-armed bandit)。对那些从没去过赌场的人,我们来做下解释:角子机(又称老虎机)是一个需要你拉杠杆(或摇臂)的赌博机器,根据机器展示的数值,你可能会得到一笔奖励,也可能(更大几率)得不到任何东西。和你想的一样,这些机器的设置都对庄家有利,所以能获的奖励的几率是非常非常小的。

多臂赌博机(理论上的)扩展了这种形式,想象你面对的是一堆角子赌博机,每个赌博机都被分配按照一个独立的概率进行奖励。作为一个玩家,你不知道在这些机器后的获奖概率,你唯一可以找到获奖概率的方法是进行游戏。你的任务是通过玩这些机器,最大限度地提高所获的奖励。那么你应该使用什么策略呢?

2

多臂赌博机策略

为了更严格地定义问题,我们通过数学形式化来表达,假设现在有k个赌博机,可观察到的每台的获奖概率等于p_k。假设一次只能拉动一个摇臂,并且赌博机只会按照它关联的概率机型奖励。这是一个设置了限定局数的有限次的游戏。在游戏期间任意时间点时,水平线H被定义为允许的剩余游戏的数量。

对所有机器用户会尝试最大化的获奖回报。在游戏中的任一时间点,我们都可以通过使用称为遗憾值(regret)来度量用户的表现。遗憾值的意思是,假设用户能在每一步选择最优的赌博机,得到的奖励和目前获得的实际奖励的差值。遗憾值的数学定义为:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

其中T表示我们到目前为止进行过的步数,r_t表示在第t步获得的奖励,u_opt表示每一局从最优赌博机返回来的期望奖励。遗憾值的数值越低,策略越优。但因为这个度量值会受到偶然性的影响(奖励可能会被从最优赌博机选择中获得的期望奖励更高),我们可以选择使用遗憾值的期望值代替,定义为:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

其中μ_t是在第t步从赌博机中获得的平均奖励(不可观测的)。因为第二项是来自所选策略的期望奖励,所以它将小于或等于来自最优策略(每一步都选择最优的赌博机)的期望奖励。

3

Epsilon优先方法

Epsilon优先(Epsilon first)是MAB策略中最简单的一种方式,它被认为和事先执行A/B测试方法具有同等意义。给定ε,执行探索空间操作的次数为(1 – ε) × N,其中N是游戏中总共的局数,剩余的次数都是执行后续探索的局数。

update_best_bandit算法会持续统计记录每一个赌博机的奖励收入和游戏局数。变best_bandit会在每一局结束进行更新,记录当前具有最高获奖概率的赌博机的编号,流程如下:

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

4

Epsilon贪婪

Epsilon贪婪(epsilon-greedy)策略中,ε表示我们进行探索空间的概率,和进行利用已知最优摇臂的事件互斥

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

该方法的特点:不需要等到探索阶段完成,才能开始利用有关赌博机的奖励表现的知识。但要小心,该算法不会考虑效果数据的统计意义。因此可能发生这样的情况:个别赌博机的奖励峰值导致后续的所有局游戏都错误地选择了这个赌博机(陈运文 达观数据)。

5

Epsilon递减

Epsilon递减(epsilon-decreasing)策略在实验开始阶段,会有一个很高的ε值,所以探索空间的可能性很高。ε值会随着水平线H上升而不断递减,致使利用似然知识的可能性更高。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

需要注意这里有几种方法去来选择一个最优的速率来更新ε值,具体取决于赌博机的数量,以及他们各自进行奖励的权重。

6

贝叶斯赌博机

与A / B测试类似,贝叶斯赌博机(Bayesian bandits)假设每个赌博机的获奖概率被建模为获奖概率的分布。当我们开始实验时,每个赌博机都有一个通用的先验概率(任意赌博机的奖励比率初始都是同等的)。

在某一个赌博机上进行的局数越多,我们对它的奖励信息就了解越多,所以基于可能的奖励概率更新其获奖概率分布。当需要选择玩哪一个赌博机的时候,从获奖概率分布中采样,并选择对应样本中具有最高奖励比率的赌博机。图3提供了在给定时间内对三个赌博机所含信息的图形化表示。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

图3

使用贝叶斯赌博机策略对三个赌博机的获奖概率信息进行建模。第1、2和3个赌博机的平均获奖率分别为0.1、0.3和0.4。 第1个赌博机具有较低的平均值而且方差也比较大,第2个赌博机具有较高的平均值和较小的方差,第3个赌博机具有更高的平均值和更小的方差。

可以看到2018免费送彩金游戏赌博机的获奖概率分布的信息被编码为三个分布。每个分布具有递增的平均值和递减的方差。因此,我们不太确定奖励期望值为0.1的真实奖励率,最可靠的是奖励期望值为0.4的赌博机。因为赌博机的选择是通过对分布进行抽样来进行的,所以分布期望值是0.1的赌博机的摇臂也可能被拉动。这个事件会发生在第2个赌博机和第3个赌博机的采样样本奖励值异常小,而且第1个赌博机的采样样本异常大时,相应代码如下(陈运文 达观数据):

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

NO.6
总结

A/B测试和贝叶斯赌博机的各自的优点和局限是:两者有各自适用的场景,也验证的变量数量也各不相同,具体如下表。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

此外,两个方法的收敛速度也很不一样。在A/B测试中是指获得统计意义,在贝叶斯赌博机中是指累积遗憾值不再增加。以本章最开始的网站优化为例,首先请注意,任何行为的改变可能是微小的(<0.01),而我们已经知道贝叶斯赌博机相比大的改变提升,需要更多的收敛时间。如果加了多种选择,在同一个实验中测试多种登陆页面,将更加会影响收敛速度。假如用户变化导致的底层分布变的比模型收敛更快呢?比如,季节趋势,销售或者其他因素可能会影响。

技术干货 | 如何选择上班路线最省时间?从A/B测试数学原理说起

显然,收集的数据越多,对效果的潜在变化的把握度就越高。当2个组划分本身就存在统计差异时,通过多臂赌博机而不是A/B测试的方法可以从概率上修正我们选择的分布。本文还重点介绍了z检验(z-test)的数学知识,因为其构成了A/B测试的统计理论基础。

正态分布的前世今生(五)

Deep Learning Specialization on Coursera

(六) 开疆扩土,正态分布的进一步发展

19世纪初,随着拉普拉斯中心极限定理的建立与高斯正态误差理论的问世,正态分布开始崭露头角, 逐步在近代概率论和数理统计学中大放异彩。在概率论中,由于拉普拉斯的推动,中心极限定理发展 成为现代概率论的一块基石。而在数理统计学中,在高斯的大力提倡之下,正态分布开始逐步畅行于天下。

1. 论剑中心极限定理

先来说说正态分布在概率论中的地位,这个主要是由于中心极限定理的影响。 1776 年,拉普拉斯开始考虑一个天文学中的彗星轨道的倾角的计算问题,最终的问题涉及 独立随机变量求和的概率计算,也就是计算如下的概率值

 S_n = X_1 + X_2 + cdots + X_n

P(a < S_n < b) = ?

在这个问题的处理上,拉普拉斯充分展示了其深厚的数学分析功底和高超的概率计算技巧,他首次引入了 特征函数(也就是对概率密度函数做傅立叶变换)来处理概率分布的神妙方法,而这一方法经过几代概率学家的发展, 在现代概率论里面占有极其重要的位置。基于这一分析方法,拉普拉斯通过近似计算, 在他的1812年发表的名著《概率分析理论》中给出了中心极限定理的一般描述:

[定理 Laplace, 1812] 假设 $ e_i (i=1, cdots n)$ 为独立同分布的测量误差, 具有均值$mu$ 和方差 $sigma^2$。如果 $lambda_1, cdots, lambda_2$ 为常数,$a>0$, 则有

 displaystyle P(|sum_{i=1}^n lambda_i(e_i - mu)| le a sqrt{sum_{i=1}^n lambda_i^2})approx frac{2}{sqrt{2pi}sigma} int_0^a e^{-frac{x^2}{2sigma^2}} dx

理科专业的本科生学习《概率论与数理统计》这门课程的时候, 除了学习棣莫弗-拉普拉斯中心极限定理,通常还学习如下中心极限定理的一般形式:

[Lindeberg-Levy 中心极限定理] 设$X_1,cdots, X_n$ 独立同分布,且具有有限的均值 $mu$ 和方差 $sigma^2$ , 则在 $n rightarrow infty$ 时,有

 displaystyle frac{sqrt{n}(bar{X} - mu)}{sigma} rightarrow N(0,1)

多么奇妙的性质,随意的一个概率分布中生成的随机变量, 在序列和(或者等价的求算术平均)的操作之下,表现出如此一致的行为,统一的规约到正态分布。 概率学家们进一步的研究结果更加令人惊讶,序列求和最终要导出正态分布的条件并不需要这么苛刻, 即便$X_1,cdots, X_n$ 并不独立,也不具有相同的概率分布形式,很多时候他们求和的最终的归宿仍然是正态分布。 一切的纷繁芜杂都在神秘的正态曲线下被消解,这不禁令人浮想联翩。 中心极限定理恐怕是概率论中最具有宗教神秘色彩的定理,如果有一位牧师拿着 一本圣经向我证明上帝的存在,我是丝毫不会买账;可是如果他向我展示中心极限定理并且声称那是神迹, 我会很乐意倾听他的布道。如果我能坐着时光机穿越到一个原始部落中,我也一定带上中心极限定理,并 劝说部落的酋长把正态分布作为他们的图腾。

中心极限定理虽然表述形式简洁,但是严格证明它却非常困难。 中心极限定理就像一张大蜘蛛网,棣莫弗和拉普拉斯编织了它的雏形,可是这张网上漏洞太多,一个多世纪来, 数学家们就像蜘蛛一样前赴后继,努力想把所有的漏洞都补上。 在十九世纪,珀松(Poission)、狄利克莱(Dirichlet)、柯西(Cauchy)、贝塞尔(Bessel)这些大蜘蛛 都曾经试图对把这张网上的漏洞补上。从现代概率论来看角度, 整个十九世纪的经典概率理论并没有能输出一个一般意义下严格的证明。 而真正把漏洞补上的是来自俄罗斯的几位蜘蛛侠:切比雪夫(Chebyshev)、马尔可夫(Markov)和李雅普诺夫(Lyapunov)。 俄罗斯是一个具有优秀的数学传统的民族,产生过几位顶尖的的数学家,在现代概率论的发展中, 俄罗斯的圣彼得堡学派可以算是顶了半边天。 把漏洞补上的严格方案的雏形是从切比雪夫1887年的工作开始的,不过切比雪夫的证明存在一些漏洞。 马尔可夫和李雅普诺夫都是切比雪夫的学生,马尔科夫沿着老师的基于矩法的思路在蜘蛛网上辛勤编织,但洞还是补得不够严实; 李雅普诺夫不像马尔可夫那样深受老师的影响,他沿着拉普拉斯当年提出的基于特征函数的思路,于1901年给出了一个补洞的方法, 切比雪夫对这个方法大加赞赏,李雅普诺夫的证明被认为是第一个在一般条件下的严格证明; 而马尔科夫也不甘示弱,在1913年基于矩法也把洞给补严实了。

20世纪初期到中期,中心极限定理的研究几乎吸引了所有的概率学家,这个定理俨然成为了概率论的明珠,成为了各大概率论 武林高手华山论剑的场所。不知道大家对中心极限定理中的“中心”一词如何理解,许多人都认为"中心"这个词描述的是这个定理的 行为:以正态分布为中心。这个解释看起来确实合情合理,不过并不符合该定理被冠名的历史。 事实上,20世纪初概率学家大都称呼该定理为极限定理(Limit Theorem),由于该定理在概率论中 处于如此重要的中心位置,如此之多的概率学武林高手为它魂牵梦绕, 于是数学家波利亚(Polya)于1920年在该定理前面冠以"中心"一词,由此后续人们都称之为中心极限定理。


论剑中心极限定理

数学家们总是及其严谨苛刻的,给定了一个条件下严格证明了中心极限定理。数学家就开始 探寻中心极限定理成立的各种条件,询问这个条件是否充分必要条件,并且进一步追问序列和在该条件下以 什么样的速度收敛到正态分布。 1922年 Lindeberg 基于一个比较宽泛容易满足的条件,给中心极限定理提出了一个很容易理解的初等证明。 这个条件我们现在称之为Lindeberg 条件。然后概率学家 Feller 和 Levy 就开始追问Lindeberg 条件是充分必要的吗? 基于 Lindeberg 的工作, Feller 和 Levy 都于 1935 年独立的得到了中心极限定理成立的充分必要条件, 这个条件可以用直观的非数学语言描述如下:

[中心极限定理充要条件]  假设独立随机变量序列 $X_i$ 的中值为0, 要使序列和 $S=sum_{i=1}^n X_i$ 的分布函数逼近正态分布,以下条件是充分必要的:

  1. 如果 $X_i$相对于序列和$S$的散布(也就是标准差)是不可忽略的,则 $X_i$ 的分布必须接近正态分布
  2. 对于所有可忽略的 $X_i$, 取绝对值最大的那一项,相对于可忽略项这个子序列和的散布,这个绝对值也是可忽略的

事实上这个充分必要条件发现的优先权,Feller 和 Levy 之间还出现了一定的争论。 在 Levy 证明这个充分必要条件的过程中, Levy发现了正态分布的一个有趣的性质。 我们在数理统计中都学过,如果两个独立随机变量 $X,Y$ 具有正态分布,则$S=X+Y$ 也具有正态分布。奇妙的是这个定理的逆定理也成立:

[正态分布的血统] 如果 $X,Y$ 是独立的随机变量,且 $S=X+Y$ 是正态分布,那么 $X,Y$ 也是正态分布。

正态分布真是很奇妙,就像蚯蚓一样具有再生的性质,你把它一刀两断,它生成两个正态分布; 或者说正态分布具有及其高贵的优良血统,正态分布的组成成分中只能包含正态分布,而不可能含有其它杂质。 1928 年 Levy 就猜到了这个定理,并使用这个定理于1935年对中心极限定理的充分必要条件作了证明。 但是 Levy 却无法证明正态分布的这个看上去及其简单的再生性质。直到 1936 年 Cramer 才给出了证明。

中心极限定理成为了现代概率论中首屈一指的定理,事实上中心极限定理在现代概率论里面已经不是指一个定理, 而是指一系列相关的定理。 统计学家们也基于该定理不断的完善拉普拉斯提出的元误差理论(the hypothesis of elementary errors), 并据此解释为何世界上正态分布如此常见。而中心极限定理同时成为了现代统计学中大样本理论的基础。

正态分布的前世今生(二)

Deep Learning Specialization on Coursera

三、最小二乘法,数据分析的瑞士军刀

第二个故事的主角是欧拉(Euler), 拉普拉斯(Lapalace),勒让德Legendre) 和高斯(Gauss),故事发生的时间是十八世纪中到十九世纪初。十七、十八世纪是科学发展的黄金年代,微积分的发展和牛顿万有引力定律的建立,直接的推动了天文学和测地学的迅猛发展。当时的大科学家们都在考虑许多天文学上的问题。几个典型的问题如下:

  • 土星和木星是太阳系中的大行星,由于相互吸引对各自的运动轨道产生了影响,许多大数学家,包括欧拉和拉普拉斯都在基于长期积累的天文观测数据计算土星和木星的运行轨道。
  • 勒让德承担了一个政府给的重要任务,测量通过巴黎的子午线的长度,
  • 海上航行经纬度的定位。主要是通过对恒星和月面上的一些定点的观测来确定经纬度。

这些天文学和测地学的问题,无不涉及到数据的多次测量,数据的计算与分析;十七、十八世纪的天文观测,也积累了大量的数据需要进行分析和计算。很多年以前,学者们就已经经验性的认为,对于有误差的测量数据,多次测量取平均是比较好的处理方法,虽然缺乏理论上的论证,也不断的受到一些人的质疑。取平均作为一种异常直观的方式,已经被使用了千百年,在多年积累的数据的处理经验中也得到一定的验证,被认为是一种良好的数据处理方法。

以上涉及的问题,我们直接关心的目标量往往无法直接观测,但是一些相关的量是可以观测到的,而通过建立数学模型,最终可以解出我们关心的量。这些天文学的问题大体都可以转换为描述如下的问题:有我们想估计的量 $beta_0,cdots,beta_p$, 另有若干个可以测量的量 $x_1,cdots,x_p, y$, 这些量之间有线性关系
 y = beta_0 + beta_1x_1 + cdots + beta_px_p

如何通过多组观测数据求解出参数$beta_0,cdots,beta_p$呢? 欧拉和拉普拉斯采用的都是求解线性方程组的方法。

begin{eqnarray}
left{
begin{array}{lll}
y_1 = beta_0 + beta_1x_{11} + cdots + beta_px_{p1}
y_2 = beta_0 + beta_1x_{12} + cdots + beta_px_{p2}
vdots
y_n = beta_0 + beta_1x_{1n} + cdots + beta_px_{pn}
end{array}
right.
end{eqnarray}

但是面临的一个问题是,有 $n$ 组观测数据,$p + 1$ 个变量, 如果 $n > p + 1$, 则得到的线性矛盾方程组,无法直接求解。 所以欧拉和拉普拉斯采用的方法都是通过一定的对数据的观察,把$n$个线性方程分为 $p+1$组,然后把每个组内的方程线性求和后归并为一个方程,从而就把$n$个方程的方程组划归为$p+1$个方程的方程组,进一步解方程求解参数。这些方法初看有一些道理,但是都过于 adhoc, 无法形成统一处理这一类问题的一个通用解决框架。

以上求解线性矛盾方程的问题在现在的本科生看来都不困难,就是统计学中的线性回归问题,直接用最小二乘法就解决了,可是即便如欧拉、拉普拉斯这些数学大牛,当时也未能对这些问题提出有效的解决方案。可见在科学研究中,要想在观念上有所突破并不容易。有效的最小二乘法是勒让德在 1805 年发表的,基本思想就是认为测量中有误差,所以所有方程的累积误差为

累积误差 = $sum($ 观测值 - 理论值 $)^2$

我们求解出导致累积误差最小的参数即可。

begin{eqnarray}
label{least-square-error}
begin{array}{lll}
hat{beta}& = & displaystyle argmin_{beta} sum_{i=1}^n e_i^2
& = & displaystyle
argmin_{beta} sum_{i=1}^n [y_i - (beta_0 + beta_1x_{1i} + cdots + beta_px_{pi})]^2
end{array}
end{eqnarray}

勒让德在论文中对最小二乘法的优良性做了几点说明:

  •  最小二乘使得误差平方和最小,并在各个方程的误差之间建立了一种平衡,从而防止某一个极端误差取得支配地位
  •  计算中只要求偏导后求解线性方程组,计算过程明确便捷
  • 最小二乘可以导出算术平均值作为估计值

对于最后一点,从统计学的角度来看是很重要的一个性质。推理如下:假设真值为 $theta$, $x_1, cdots, x_n$为n次测量值, 每次测量的误差为$ e_i = x_i - theta $,按最小二乘法,误差累积为

 L(theta) = sum_{i=1}^n e_i^2 = sum_{i=1}^n (x_i - theta)^2

求解$theta$ 使得 $L(theta)$达到最小,正好是算术平均 $bar{x} = frac{sum_{i=1}^n x_i}{n} $。

由于算术平均是一个历经考验的方法,而以上的推理说明,算术平均是最小二乘的一个特例,所以从另一个角度说明了最小二乘方法的优良性,使我们对最小二乘法更加有信心。

最小二乘法发表之后很快得到了大家的认可接受,并迅速的在数据分析实践中被广泛使用。不过历史上又有人把最小二乘法的发明归功于高斯,这又是怎么一回事呢。高斯在1809年也发表了最小二乘法,并且声称自己已经使用这个方法多年。高斯发明了小行星定位的数学方法,并在数据分析中使用最小二乘方法进行计算,准确的预测了谷神星的位置。

扯了半天最小二乘法,没看出和正态分布有任何关系啊,离题了吧?单就最小二乘法本身,虽然很实用,不过看上去更多的算是一个代数方法,虽然可以推导出最优解,对于解的误差有多大,无法给出有效的分析,而这个就是正态分布粉墨登场发挥作用的地方。勒让德提出的最小二乘法,确实是一把在数据分析领域披荆斩棘的好刀,但是刀刃还是不够锋利;而这把刀的打造后来至少一半功劳被归到高斯,是因为高斯不单独自的给出了造刀的方法,而且把最小二乘这把利刀的刀刃造得无比锋利,把最小二乘打造为了一把瑞士军刀。高斯拓展了最小二乘法,把正态分布和最小二乘法联系在一起,并使得正态分布在统计误差分析中确立了自己的定位,否则正态分布就不会被称为高斯分布了。 那高斯这位神人是如何把正态分布引入到误差分析之中,打造最小二乘这把瑞士军刀的呢?看下一个故事。

正态分布的前世今生(一)

Deep Learning Specialization on Coursera

神说,要有正态分布,就有了正态分布。
神看正态分布是好的,就让随机误差就服从了正态分布。

创世纪-数理统计

一、正态分布

学过基础统计学的同学大都对正态分布非常熟悉。这个钟型的分布曲线不但形状优雅,其密度函数写成数学表达式
 displaystyle f(x)=frac{1}{sqrt{2pi}sigma}e^{-frac{{(x-mu})^2}{2sigma^2}}
也非常具有数学的美感。其标准化后的概率密度函数
 displaystyle f(x)=frac{1}{sqrt{2pi}}e^{-frac{x^2}{2}}
更加的简洁漂亮,两个最重要的数学常量 $ pi, e$ 都出现在了公式之中。在我个人的审美之中,它也属于 top-N 的最美丽的数学公式之一,如果有人问我数理统计领域哪个公式最能让人感觉到上帝的存在,那我一定投正态分布的票。因为这个分布戴着神秘的面纱,在自然界中无处不在,让你在纷繁芜杂的数据背后看到隐隐的秩序。

正态分布又通常被称为高斯分布,在科学领域,冠名权那是一个很高
的荣誉。去过德国的兄弟们还会发现,德国的钢镚和10马克的纸币上都留有高斯的头像和正态密度曲线。正态分布被冠名高斯分布,我们也容易认为是高斯发现了正态分布,其实不然,不过高斯对于正态分布的历史地位的确立是起到了决定性的作用。

正态曲线虽然看上去很美,却不是一拍脑袋就能想到的。我在本科学习数理统计的时候,课本一上来介绍正态分布就给出密度分布函数,却从来不说明这个分布函数是通过什么原理推导出来的。所以我一直搞不明白数学家当年是怎么找到这个概率分布曲线的,又是怎么发现误差服从这个奇妙的分布的。直到我读研究生的时候我的导师给我介绍了陈希儒院士的《数理统计简史》这本书,看了之后才了解了正态分布曲线从发现到被人们重视进而广泛应用,也是经过了几百年的历史。

正态分布的这段历史是很精彩的,我们通过讲几个故事来揭开她的神秘面纱。

二、邂逅,正态曲线的首次发现
第一个故事和概率论的发展密切相关,主角是棣莫弗(De Moivre) 和拉普拉斯(Laplace)。

拉普拉斯是个大科学家,被称为法国的牛顿;棣莫弗名气可能不算很大,不过大家应该应该都熟悉这个名字,因为我们在高中数学学复数的时候都学过棣莫弗定理$(costheta + i sintheta)^n = cos(ntheta) + i sin(ntheta)$。

古典概率论发源于赌博,惠更斯、帕斯卡、费马、贝努力都是古典概率的奠基人,他们那会研究的概率问题大都来自赌桌上,最早的概率论问题是赌徒梅累在1654年向帕斯卡提出的如何分赌金的问题。统计学中的总体均值之所以被称为期望(Expectation), 就是源自惠更斯、帕斯卡这些人研究平均情况下一个赌徒在赌桌上可以期望自己赢得多少钱。

有一天一个哥们,也许是个赌徒,向棣莫弗提了一个和赌博相关
的一个问题:A,B 两人在赌场里赌博,A,B各自的获胜概率是$p, q=1-p$,赌 n 局,若 A 赢的局数 $X > np$, 则 A 付给赌场 $X-np$ 元,否则B 付给赌场 $np-X$ 元。 问赌场挣钱的期望值是多少。

问题并不复杂, 本质上是一个二项分布,最后求出的理论结果是
 2npq b(n, p, np)
其中 $b(n,p,i) = binom{n}{i}p^iq^{n-i}$ 是常见的二项概率。 但是对具体的 $n$, 要把这个理论结果实际计算出数值结果可不容易, 因为其中的二项公式中有组合数.这就驱动 De Moivre寻找近似计算的方法计算。

与此相关联的另一个问题,是遵从二项分布的随机变量 $X sim B(n,p)$, 求X 落在二项分布中心点一定范围的概率 $P_d = P(|X - np| le d)$

对于 p=1/2 的情形, 棣莫弗 做了一些计算并得到了一些近似结果,但是还不够漂亮,幸运的是 棣莫弗 和 Stirling 处在同一个时代, 而且二人之间有联系,Stirling 公式是在数学分析中必学的一个重要公式(事实上Stirling 公式的形式其实是棣莫弗最先发现的,但是 Stirling 改进了公式)

 displaystyle n! sim sqrt{2pi n} (frac{n}{e})^{n}

1733 年,棣莫弗很快利用 Stirling 公式进行计算并取得了重要的进展。考虑 n 是偶数的情形,令二项概率
b(i) = b(n, frac{1}{2}, i) = binom{n}{i}(frac{1}{2})^n
通过 Stirling 公式做一些简单的计算容易得到,

 displaystyle b(frac{n}{2}) sim sqrt{frac{2}{pi n}}
 displaystyle frac{b(frac{n}{2}+d)}{b(frac{n}{2})} sim e^{-frac{2d^2}{n}}

于是有
 displaystyle b(frac{n}{2}+d) sim frac{2}{sqrt{2 pi n}}e^{-frac{2d^2}{n}}

使用上式的结果,并在二项概率累加求和的过程中近似的使用定积分代替求和,很容易就能得到

 displaystyle P(|frac{X}{n} - frac{1}{2}| le frac{c}{sqrt{n}} ) simint_{-2c}^{2c} frac{1}{sqrt{2pi}} e^{-x^2/2} dx

看,正态分布的密度函数的形式在积分公式中出现了!这也就是我们在数理统计课本上学到的二项分布的极限分布是正态分布。

以上只是讨论了 $p=1/2$ 的情形, 棣莫弗也对 $p ne 1/2$做了一些计算,后来拉普拉斯对 $p ne 1/2$ 的情况做了更多的分析,并把二项分布的正态近似推广到了任意 $p$ 的情况。 这是第一次正态密度函数被数学家勾画出来,而且是以二项分布的极限分布的形式被推导出来的。 熟悉基础概率统计的同学们都知道这个结果其实叫棣莫弗-拉普拉斯中心极限定理。

[De Moivre-Laplace 中心极限定理]
设随机变量 $X_n (n=1,2,cdots)$ 服从参数为 $p$ 的二项分布,则对任意的 $x$, 恒有
displaystylelim_{nrightarrowinfty}P{ frac{X_n - np}{sqrt{np(1-p)}} le x }=int_{-infty}^x frac{1}{sqrt{2pi}} e^{frac{-t^2}{2}}dt

我们在大学学习数理统计的时候,学习的过程都是先学习了正态分布,然后才学习中心极限定理。而学习到正态分布的时候,直接就描述了其概率密度的数学形式,虽然数学上很漂亮,但是当时很容易困惑数学家们是如何凭空就找到这个分布的。读了陈希孺的《数理统计学简史》之后,我才明白正态分布的密度形式首次发现是在棣莫弗-拉普拉斯的中心极限定理中。数学家研究数学问题的进程很少是按照我们数学课本的安排顺序推进的,现代的数学课本都是按照数学内在的逻辑进行组织编排的,虽然逻辑结构上严谨优美,却把数学问题研究的历史痕迹抹得一干二净。DNA 双螺旋结构的发现者之一 Waston 在他的名著《DNA 双螺旋》序言中说:“科学的发现很少会像门外汉所想象的一样,按照直接了当合乎逻辑的方式进行的。”

棣莫弗 出他的发现后40年(大约是 1770), 拉普拉斯建立了中心极限定理较一般的形式,中心极限定理后续又被其它数学家们推广到了其它任意分布的情形,而不限于二项分布。后续的统计学家发现,一系列的重要统计量,在样本量 N 趋于无穷的时候, 其极限分布都有正态的形式, 这构成了数理统计学中大样本理论的基础。

棣莫弗在二项分布的计算中瞥见了正态曲线的模样,不过他并没有能展现这个曲线的美妙之处。棣莫弗的这个工作当时并没有引起人们足够的重视,原因在于棣莫弗 不是个统计学家,从未从统计学的角度去考虑其工作的意义。 正态分布(当时也没有被命名为正态分布) 在当时也只是以极限分布的形式出现,并没有在统计学,尤其是误差分析中发挥作用。这也就是正态分布最终没有被冠名 棣莫弗分布的重要原因。 那高斯做了啥工作导致统计学家把正态分布的这顶桂冠戴在了他的头上呢?这先得从最小二乘法的发展说起。下回分解:-)