机器学习项目中不可忽视的一个密辛 - 大数定理、中心极限定理

1. 前言

0x1：概率论的起源

大数定律和中心极限定律是伴随着古典统计，古典概率论，函数分析，极限理论，现代概率论与数理统计这些理论一起发展起来的。科学家在提出和发展概率论理论体系的过程中，一个最大也是最根本的挑战就是，如何从理论上强有力地证明概率论是一个理论合理的数学理论，而不仅仅是基于大量的实际实验而得到的归纳假设，诚然，很多科学理论都是从科学实验里归纳总结得到的一般性推理，例如亚里士多德的重力理论，但是最终它们都必须被合理地嵌入在某个理论体系中，和其他理论互相推导成立，才能真正成为一门科学。

在概率论发展的初期，人们主要还是以实验归纳为主，人们发现，当大量重复某一相同的实验的时候，其最后的实验结果可能会稳定在某一数值附近。就像抛硬币一样，当我们不断地抛，抛个上千次，甚至上万次，我们会发现，正面或者反面向上的次数都会接近一半。除了抛硬币，现实中还有许许多多这样的例子，像掷骰子，最著名的实验就是泊松抛针实验。这些实验都像我们传达了一个共同的信息，那就是大量重复实验最终的结果都会比较稳定。那稳定性到底是什么？怎样去用数学语言把它表达出来？这其中会不会有某种规律性？是必然的还是偶然的？

17 世纪中叶，人们开始对机会性游戏的数学规律进行探讨。它的发展与数学史上一些伟大的名字相联系，如帕斯卡、费马、惠更斯、詹姆斯、伯努利、棣莫弗、拉普拉斯等。

1654 年，费马与帕斯卡的通信中关于分赌注问题的讨论被公认为是概率论诞生的标志。问题是这样的：

“两个赌徒相约赌若干局，谁先赢 s 局就算赢了，当赌徒 A 赢 a 局 (a < s)，而赌徒 B 赢 b 局(b < s) 时，赌博被迫中止，应该怎样分配赌注才合理？”

在三年后，惠根斯亦用自己的方法解决了这一问题，并写成了《论赌博中的计算》一书，这就是概率论最早的论著，他们三人提出的解法中，都首先涉及了数学期望（mathematical expectation）这一概念，并由此奠定了古典概率的基础。

经过几百年的发展，大数定律体系已经很完善了，也出现了更多更广泛的大数定律，例如切比雪夫大数定律，辛钦大数定律，泊松大数定律，马尔科夫大数定律等等。正是这些数学家们的不断研究，大数定律才得以如此迅速发展，才得以完善。

在这篇文章中，笔者将尝试从时间线的维度，介绍几个核心的大数定理和中心极限定律的发展脉轮，并讨论其应用场景。

Relevant Link: 

https://baike.baidu.com/item/%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B

 

2. 伯努利大数定律（1713年）

自然科学的发展总是从简入难，从物理、生活中的日常现象发现规律，并进行归纳总结逐步发展起来的。

历史上首个极限定理是由瑞士数学家雅各布.伯努利（1654年-1705年），在其遗著中发表了概率论中的第一个极限定理——伯努利大数定理，即 “在多次重复独立的试验中，事件发生的频率有越来越稳定的趋势”。

这正是频率稳定性的定理形式，它的出现意味着概率论由建立走向发展的阶段。

0x1：数学公式定义

设 na 是n重伯努利试验中，事件A发生的次数，p是事件A在每次事件中发生的概率，则对任意的ε>0，有：

或

在n重伯努利试验中，设：

则 X1，X2，....，Xn 是相互独立同分布的随机变量，且Xi~B(1，p)，i=1，2，...，n，其中 p=P(A)。

由伯努利大数定律知，当n充分大时，有如下等价近似公式：

A发生的频率

这是概率论第一个非常深刻的认识，当大量相互独立重复试验中，可以用某个事件A发生的频率来近似每次试验中事件A发生的概率。这就是伯努利大数定律的直观意义。

当n充分大时，频率与其概率能以任意接近的概率趋向于1。

因此实际中，只要试验次数足够多，可以用频率作为概率的估计。

同时伯努利大数定律也解释了概率存在的客观意义，即为什么“大数次”重复试验下，事件的概率是存在的，正是因为频率的这种稳定性，我们才意识到概率的存在，才有了概率论这门学科。

Relevant Link: 

https://baike.baidu.com/item/%E4%BC%AF%E5%8A%AA%E5%88%A9%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B/7922139

 

3. 棣莫弗-拉普拉斯中心极限定理（1733年）

1730 年，法国数学家棣莫弗（1677年-1754年）出版的著作《分析杂论》中包含了著名的棣莫弗─拉普拉斯定理。他使用正态分布取估计n（很大）时抛掷硬币出现正面次数的分布，即二项分布B(n，0.5)。这就是概率论中第二个基本极限定理的雏形。

将近80年后，拉普拉斯（1749年-1827年）在 1812 年出版的《概率的分析理论》中，首先明确地对概率作了古典的定义。他指出当n很大时，二项分布B(n，p)(0<p<1)都可以用正态分布逼近。所以后人称之为棣莫弗-拉普拉斯中心极限定理。

另外，拉普拉斯他又和数个数学家建立了关于正态分布及最小二乘法的理论。

0x1：数学公式定义

定理核心公式如下：

设随机变量序列X1，X2，..... 是一个独立同分布的随机变量序列（注意这是二项分布的前提假设），且Xi ~ B(1，p)，i = 1，2，....，则对任意实数x，有：

仔细看这个公式，会发现几点：

• 分子就是二项分布随机变量序列的累加和，减去随机变量的期望np
• 分布就是二项分布的标准差
• 分子除分母的整个公式，可以看成是随机变量期望值标准化的过程
• 同时括号里的<=符号，表达说标准化的期望偏离0原点的概率，也即误差概率，这个偏离误差概率是符合正态分布的

德莫弗-拉普拉斯中心极限定理不仅在概率论发展的早期起过重要的作用，而且在工程实践中依然被大量使用，在大数据前提下，基于该定理的近似等价计算方法在复杂工程问题中非常方便。

0x2：由该定理得到二项分布的正态分布近似解

由二项分布的可加性知道，因此我们有：

这就是二项分布的原始概率计算公式，精确的概率值理论上是可以精确算出的。

但是，实际问题中当n较大时，计算并不方便，这时候就需要借助极限、近似等价这2个数学工具来帮助我们简化计算。

泊松定理曾经告诉我们，当p ≤ 0.1时，可以用泊松分布作近似计算，这个近似等价其实已经很好了，但是这个假设还是有些局限，如果p比较大，则无法继续用泊松近似。

现在棣莫弗-拉普拉斯中心极限定理告诉我们，也可以用正态分布作近似计算，它的优点是不受“p≤0.1”的限制，只需n足够大。

由德莫弗 -拉普拉斯中心极限定理推得，如果随机变量Y~B(n，p)，那么，当n较大时，有如下近似等价公式：

显然，棣莫弗-拉普拉斯中心极限定理是列维-林德伯格中心极限定理的特例，因为服从二项分布，所以这个定理又称为二项分布的正态近似。

0x3：棣莫弗-拉普拉斯中心极限定理对伯努利大数定理的更细致证明

前面由伯努利大数定理知，当n充分大时，可以用作为p的近似，但是至于近似程度如何，伯努利大数定律没有给出。

中心极限定律对近似的程度进行了注释：

中心极限定理对频率趋近于概率的论证更为细致。

0x4：棣莫弗-拉普拉斯定理应用举例

1. 题目

一本20万字的长篇小说进行排版。假定每个字被错排的概率为10-5。试求这本小说出版后发现有6个以上错字的概率，假定各个字是否被错排是相互独立的。

2. 解题分析

显然，这道题符合二项分布，我们应该用二项分布公式来解题，但是注意这里n有20w，是一个比较大的数字，因此我们采用棣莫弗-拉普拉斯定义或者泊松定理来计算近似等价结果。

设错字总数为随机变量X ，X~B (200000，10^-5 )，np=2，

所以有正态分布近似等价结果：

以及泊松分布近似等价结果（λ=2）：

Relevant Link: 

https://en.wikipedia.org/wiki/De_Moivre%E2%80%93Laplace_theorem 

 

4. 列维-林德伯格中心极限定理（1920年）

我们将棣莫弗-拉普拉斯定理推广到一般化的相互独立同部分随机变量中，即列维-林德伯格中心极限定理。列维（1886年-1971年）是法国数学家，对极限理论和随机过程理论做出了杰出的贡献，林德伯格（1876年-1932年）是芬兰数学家因中心极限定理而闻名于世。

0x1：数学公式定义

设X1，X2，... 是任意一个独立同分布的随机变量序列，且均值和方差存在（随机变量序列收敛），即：

则对任意一个x，-∞ < x < ∞，总有：

其中，Φ(x)是 N(0，1) 的正态分布函数。

上式也被称为列维 -林德伯格(Levy-Lindberg)中心极限定理。

由于：

因此，定理中的概率实际上是的标准化随机变量的分布函数值，任意的概率分布函数都可以进行标准化。

列维-林德伯格中心极限定理告诉我们，不论X1，X2，... 原来服从什么分布，当n足够大时，总可以近似地认为：

或者

在实际问题中，若n较大，可以利用正态分布近似求得概率：

这在工程实践中特别有用。

0x2：公式的合理性分析

1. 微观少数样本的随机构成宏观统计意义上的有序

我们在讨论正态分布的时候，曾经阐述过正态分布存在的合理性。当一个变量，受到大量微小的，独立的随机因素（即X1，X2，... Xn）的累计影响时，这个变量一般服从正态分布，也可以反过来理解为，正是因为这些大量微小的，独立的随机因素（即X1，X2，... Xn）的累计影响，导致最终从宏观上呈现出正态分布的结果。中心极限定理正是这种直观经验的严格数学表达。

定理的条件要求随机变量相互独立并且服从同一分布。

这里相互独立意味着随机变量之间不相互影响；

同分布是指每个随机变量在随机变量序列的前n项部分和中的地位相同，也即每个随机变量对前n项部分和的影响都是微小的。

这就解释了自然界中一些现象受到许多相互独立且微小的随机因素影响，总的影响就可以看作服从或近似服从正态分布。

2. 高尔顿钉板实验

为了更好地阐明“微观少数样本的随机性”以及“宏观大数统计意义上的正态有序”这2个概念，我们来看一个历史上著名的例子，高尔顿钉板实验。

弗朗西斯·高尔顿（Francis Galton，1822年2月16日－1911年1月17日），英国科学家和探险家。他曾到西南非洲探险，因树立功绩而知名并被选为英国皇家地理学会会员，三年后又入选英国皇家学会，晚年受封为爵士。
他的学术研究兴趣广泛，包括人类学、地理、数学、力学、气象学、心理学、统计学等方面。他是查尔斯·达尔文的表弟，深受其进化论思想的影响，把该思想引入到人类研究。
他着重研究个别差异，从遗传的角度研究个别差异形成的原因，开创了优生学。他关于人类官能的研究开辟了个体心理和心理测验研究的新途径。
他在统计学方面也有突出的贡献，高尔顿在1877年发表关于种子的研究结果，指出回归到平均值（regression toward the mean）现象的存在，这个概念与现代统计学中的“回归”并不相同，但是却是回归一词的起源。

下图是高尔顿钉板实验的示意图，实验过程如下：

每一圆圈点表示钉在板上的一颗钉子，它们彼此的距离均相等，上一层的每一颗的水平位置恰好位于下一层的两颗正中间，总共有n排钉子，相当于做n次伯努利实验。

从入口处放进一个直径略小于两颗钉子之间的距离的小圆玻璃球，当小圆球向下降落过程中，碰到钉子后皆以1/2的概率向左或向右滚下，于是又碰到下一层钉子。如此继续下去，直到滚到底板的一个格子内为止。

把许许多多同样大小的小球不断从入口处放下，只要球的数目相当大，它们在底板将堆成近似于正态 的密度函数图形（即：中间高，两头低，呈左右对称的古钟型）。

来稍微分析一下小球下落的结果，小球堆积的形态取决于小球最终下落在底部隔板的位置的分布。设随机变量X为“小球最终下落在底部隔板中的位置”，X的概率分布由一系列的随机变量序列组成，即：

• Xi = -1：小球碰到第 i 排钉子后向左下落
• Xi = 1：小球碰到第 i 排钉子后向右下落

显然对单次实验（投掷小球）来说，小球会落在哪个方格Xi完全是随机的，几乎看不出什么规律，这就是所谓的“微观少数样本的随机性”，这也是概率论的精髓思想，表面随机的事物背后都蕴含着某种确定的统计规律性。

从概率分布的角度来看，显然有，随机变量X表示最终小球下落的结果，是由所有随机变量序列累加得到的，和的分布计算很复杂。但是高尔顿经过试验发现，随着试验次数的增加，小球的堆积形态呈现出正态分布的形态。这就是所谓的“宏观统计意义上的有序性”。

经过叠加，原本取值任意一点的可能性由相同变为了向中心位置聚拢，这和我们的直觉是多么不同。 

Relevant Link: 

https://baike.baidu.com/item/%E9%AB%98%E5%B0%94%E9%A1%BF%E9%92%89%E6%9D%BF/6765470
https://baike.baidu.com/item/%E6%9E%97%E5%BE%B7%E4%BC%AF%E6%A0%BC%E5%88%97%E7%BB%B4%E5%AE%9A%E7%90%86/10388983?fr=aladdin 

 

5. 辛钦大数定律（相互独立同部分大数定律）

辛钦大数定律继承并发展了伯努利大数定律。

辛钦大数定律面向的是一般化的相互独立同分布的随机变量情况，并不局限于伯努利随机试验。伯努利大数定律是相互独立同分布大数定律的特例，这点我们接下来会讨论。

0x1：数学公式定义

则随机变量序列 X1，X2，....，Xn 相互独立同分布（任意两个随机变量线性无关），若均值存在，即

则对任意ε>0，有：

也可以表示成：

在许多实际问题中，方差存在不一定满足，苏联数学家辛钦（1894年-1959年）证明了相互独立同分布情形下，仅期望存在、方差不存在时结论仍然成立，因此相互独立同分布大数定律又称为辛钦大数定律。

0x2：算数平均值法则的理论依据

相互独立同分布（辛钦）大数定律是我们日常生活中经常使用的算术平均值法则的理论依据。 

为了精确称量物体的质量μ，可在相同条件下重复称n次，结果可记为x1，x2，....，xn。可看做n个相互独立同分布的随机变量（X1，X2，....，Xn）的一次观测值。

X1，X2，....，Xn服从同一分布，它们共同的期望记为物体的真实质量μ。

由相互独立同分布（辛钦）大数定律可知，当n充分大时有：

。

这意味着逐渐趋向于μ，也即随机变量的算数平均值具有稳定性。

在物理实验中我们就是采用这种方法测得物体质量的，例如，为得到一颗钻石的真实质量，我们测n次取其算法平均值即可。

算法平均值法则提供了一条切实可操作的途径来得到物体的真实值。大数定律从理论上给出了这个结论的严格证明，而不是仅仅靠直觉。

Relevant Link: 

https://baike.baidu.com/item/%E8%BE%9B%E9%92%A6%E5%A4%A7%E6%95%B0%E5%AE%9A%E5%BE%8B/7922188

 

6. 切比雪夫大数定律

切比雪夫大数定律是对伯努利大数定律和辛钦大数定律的继承和发展。是由俄国数学家切比雪夫（1821年-1894年）发表的。

切比雪夫大数定律是面向所有两两不相关且方差一致有界的随机变量，并不要求随机变量独立同分布。相互独立同分布的大数定律是切比雪夫大数定律的特例，因为切比雪夫大数定律中的在相互独立同分布大数定律的条件下为μ。

切比雪夫大数定律比较理论化，在讨论具体的数学公式之前，我们需要先介绍几个数学概念。

0x1：切比雪夫（Chebyshev）不等式

随机变量X的取值总是围绕着其期望变动，若X的分布已知时，可以计算事件的概率。

切比雪夫不等式给出了任意概率分布X时，的上限计算公式。

设随机变量X的数学期望E(X)及方差D(X)存在，则对于任意的ε>0，有：

事件理解为：“随机变量X关于其期望发生了较大偏差”，切比雪夫不等式给出了此事件的概率上界，它与方差成正比：

• 方差越大，此上界就越大
• 方差越小，X在其期望附近取值的密集程度就越高，那么远离期望的区域的概率上界就越小

我们对切比雪夫公式进行一个变换使其更好地体现出“概率和偏差的反比关系”，

，令，则有，

从上式中，可以很清楚看到，X偏离E[X]两个标准差距离的概率小于1/4。 

切比雪夫不等式进一步说明了方差的概率意义，即方差是随机变量取值与其中心位置的偏离程度的一种度量指标。

值得注意的是，显然利用切比雪夫不等式估计“随机变量X关于其期望μ发生了较大偏差”的概率是粗糙的，引入切比雪夫不等式的另一个目的是，它是证明大数定律的工具之一。

0x2：依概率收敛

在前面第二小节讨论到伯努利大数定律的时候，曾经提到频率的稳定性。

设随机事件A的概率P(A)=p，在n重伯努利试验中事件A发生的频率为fn(A)，当n很大时，fn(A)将与p非常接近。

自然会想到，应该用极限概念来描述这种稳定性，但是不能简单地使用高等数学中数列的极限，因为fn(A)本质上是一个随机变量，它随着不同的n次试验可能取不同的值，这就需要对随机变量序列引进新的收敛性定义。

设X1，X2，... 是一个随机变量序列，如果存在一个常数c，使得对任意一个ε>0，总有：

那么，称随机变量序列X1，X2，... 依概率收敛于c，记作

依概率收敛性的直观意义是，当n足够大时，随机变量Xn几乎总是取接近于常数c的值。

利用求对立事件的概率计算公式，依概率收敛性也可以等价地表示成：

一般地，Xi 不一定相互独立，也不一定服从0-1分布，把具有这种形式的依概率收敛性的结论统称为大数定律，即一般化的大数定律。

0x3：切比雪夫大数定律数学公式

设 X1，X2，... 是两两不相关（注意两两不相关并不是独立同分布，只是两两不线性相关）的随机变量序列。如果存在常数c，使得 D(Xi)≤c，i=1，2，...，那么有：

由切比雪夫不等式推得，对任意一个ε>0，当n → ∞时，有：

要特别注意的是，切比雪夫大数定律的主要条件是“方差有界”。

在特殊条件下，即随机变量相互独立同分布，即E(Xi) = μ，i=1，2，... ，则切比雪夫大数定律等价于辛钦大数定律，即：

0x4：三个大数定律之间的条件关系

三个大数定律条件是不同的：

• 切比雪夫大数定律不要求随机变量序列同分布，甚至不要求相互独立，只要两两不线性相关、方差一致有界即可
• 辛钦大数定律要求随机变量相互独立且同分布，但不要求方差存在，仅期望存在即可
• 努利大数定律要求随机变量相互独立且同分布，但限定于伯努利两点分布

 

7.  大数定律和中心极限定理的异同

0x1：中心极限定理和大数定律的区别

中心极限定理是随机变量和的分布收敛到正态分布的一类定理，而随机变量的和又和随机变量的均值有密切的联系，而大数定律论证的主要部分就是随机变量均值的收敛性特点，因此，中心极限定理和大数定律之间有千丝万缕的联系。

• 大数定律是说，n只要越来越大，把这n个独立同分布的数加起来去除以n得到的这个样本均值（也是一个随机变量）会依概率收敛到真值u，但是样本均值的分布是怎样的我们不知道
• 中心极限定理是说，n只要越来越大，这n个数的样本均值会趋近于正态分布，并且这个正态分布以u为均值，sigma^2/n为方差
• 综上所述，这两个定律都是在说样本均值性质。随着n增大，大数定律说样本均值几乎必然等于均值。中心极限定律说，它越来越趋近于正态分布，并且这个正态分布的方差越来越小

0x2：中心极限定理在概率论中的核心地位

不同的中心极限定理的差异就在于对随机变量序列做出了不同的假设，由于中心极限定理的有力支撑使正态分布在概率论与数理统计中占据了独特的核心地位，它是20世纪初概率论研究的中心内容。

在大数n的情况下，所有概率分布函数都会收敛到正态分布这个中心极限的形式上。同时各个概率分布函数之间也存在互相推导的关系，这也是为什么GMM从理论上可以拟合任意复杂的概率分布函数的原因。

每当这时，笔者就非常喜欢祭出一张非常出名且有趣的图：

 

8. 大数定律的理论应用 

大数定律支撑了很多应用理论的发展，例如：

• 1）算术平均法则
• 2）频率估计概率方法
• 3）数理统计中参数的点估计思想，这也是很多传统机器学习模型训练的理论支撑
• 4）中心极限定理在数理统计的区间估计与假设检验问题中的应用  

 

9. 大数定律存在的一个实验例证

0x1：泊松随机变量的收敛 - 在大数N时，大数定理普遍存在的一个例证

下面是三个不同的泊松随机变量在大数定理上的实例，用三个不同的曲线表示。 三个泊松分布的参数𝜆 = 4.5，样本个数 sample_size=100000 个。

泊松分布本身也是一个随机过程，参数 𝜆 一样，只是表明它们符合同一个概率分布，但是泊松分布本身每次的取值依然是一个随机过程。这里用三个泊松分布来进行实验，只是想说明不管过程如何随机，但是其期望依然会收敛到一个具体的值上，就像有一只魔法之手在牵引着每次的取值。

下面代码计算𝑛 = 1到 sample_size 时样本的均值：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

figsize(12.5, 5)
import pymc as pm

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)
print N_samples

for k in range(3):

    samples = poi(lambda_, size=sample_size)

    partial_average = [samples[:i].mean() for i in N_samples]   # 从[1:sample_size]不同样本数下，三个泊松分布的均值

    plt.plot(N_samples, partial_average, lw=1.5, label="average \
of  $n$ samples; seq. %d" % k)


plt.plot(N_samples, expected_value * np.ones_like(partial_average),
         ls="--", label="true expected value", c="k")

plt.ylim(4.35, 4.65)
plt.title("Convergence of the average of \n random variables to its \
expected value")
plt.ylabel("average of $n$ samples")
plt.xlabel("# of samples, $n$")
plt.legend()

plt.show()

从上图中，我们可以清晰地看出，当 n 很小时，样本均值的变动很大，随着 n 的增大都在逐渐趋向平缓，最终在 4.5 水平线附近小幅摆动。数学家和统计学家将这种“小幅摆动”称为收敛。

我们所熟知的定理：泊松分布的均值等于其参数𝜆。这里所谓的泊松分布指的是一个真实概率分布，但是真实概率分布只有上帝知道我们并不知道，但是当样本数足够大时，根据大数定理，我们可以用一定数量的样本均值来近似代替真实分布的均值，进而得到真实分布的期望。

1. 收敛到期望值的速度

明白了分布的收敛特性之后，另一个与此相关的问题是：收敛到期望值的速度有多快？

1）通过实验方式来观察收敛速度

我们选择一个特定的 n，然后重复实验数千次，并计算与实际期望值的距离的均值：

以上公式可以理解为：前 n 个采样的样本均值与真实值的距离（平均意义上的）

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)

N_Y = 250  # use this many to approximate D(N)
N_array = np.arange(1000, 50000, 2500)  # use this many samples in the approx. to the variance.
D_N_results = np.zeros(len(N_array))

def D_N(n):
    """
    This function approx. D_n, the average variance of using n samples.
    """
    Z = poi(lambda_, size=(n, N_Y))
    average_Z = Z.mean(axis=0)
    return np.sqrt(((average_Z - expected_value) ** 2).mean())


for i, n in enumerate(N_array):
    D_N_results[i] = D_N(n)


plt.xlabel("$N$")
plt.ylabel("expected squared-distance from true value")
plt.plot(N_array, D_N_results, lw=3,
         label="expected distance between\n\
expected value and \naverage of $N$ random variables.")
plt.plot(N_array, np.sqrt(expected_value) / np.sqrt(N_array), lw=2, ls="--",
         label=r"$\frac{\sqrt{\lambda}}{\sqrt{N}}$")
plt.legend()
plt.title("How 'fast' is the sample average converging? ")

plt.show()

如我们预期，随着 N 增大，样本均值与实际值间距的期望逐渐减小。但也注意到，这个收敛速率也在降低。

样本点从 10000 增加到 20000 时，均值和期望的距离从 0.020 降低到 0.015，减少了 0.005；当样本点从 20000 增加到 40000 时，样本均值和期望之间的距离从 0.015 降低到 0.010，仅仅减少了 0.005。 说明收敛的速率在逐渐变小。

2）通过公式方式来描述收敛速率

实际上收敛的速率是可以通过公式衡量的，在上图的橙色曲线中，是通过下面的公式生成的：

该公式模拟了样本的均值和真实期望之间的距离。

这个特性是非常有趣的，因为对于一个给定的𝑁，就可以知道估计出的样本平均值和真实的期望值之间的误差大小。也很容易直观地想象，如果当前方差很小了，说明摆动幅度已经很小了，基本上肯定是快要进入收敛状态了，所以收敛速率自然也减小了。

2. 连接概率和频率之间的桥梁 - 大数定理

频率统计 <---->(大数定理) <----> 期望 <----> 概率

在这里，大数定理充当了一个桥梁的作用，将频率统计和概率连接了起来。而连接的水泥就是期望。

同时，我们也要明白，大数定理在 N 很小的时候会失效，在 N 很小的时候，此时得到的近似结果是不可靠的。

0x2：小数据的无序性 - 当 N 很小时大数定理会失效

只有在 N 足够大时，大数定理才成立，然后数据量并非总是足够大的。如果任意运用这一定理，就会导致很荒谬的错误，我们这章来讨论这个问题。

1. 实例1 - 地理数据的聚合

有些数据会以某一属性进行分组。例如有些数据按照国家，县区或者城市进行分组。当然人口数量也会根据地理区域的不同而不同。如果要对某一地理区域中的一些特征进行统计平均时，必须要注意大数定理适用范围，并了解在人口数较少时它是如何失效的。
假设在数据集 中包含有 5000 个县的数据。除此之外人口数量在整个国家的分布符合 100 到 1500 的均匀分布。 我们关心的是每个县中人口的平均高度。

我们假设不论生活在哪里，每个人的身高都服从相同的分布：height~Normal(150,15)

把数据按照县级别进行分组，所以只能对每个县里的数据做平均。数据如下图：

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)
std_height = 15
mean_height = 150

n_counties = 5000
pop_generator = pm.rdiscrete_uniform
norm = pm.rnormal

# generate some artificial population numbers
population = pop_generator(100, 1500, size=n_counties)

average_across_county = np.zeros(n_counties)
for i in range(n_counties): # 5000个县
    # generate some individuals and take the mean
    average_across_county[i] = norm(mean_height, 1. / std_height ** 2,  # 每个县都服从相同的分布: height~Normal(150,15)
                                    size=population[i]).mean()　　# 不同的县人口

# located the counties with the apparently most extreme average heights.
i_min = np.argmin(average_across_county)
i_max = np.argmax(average_across_county)

# plot population size vs. recorded average
plt.scatter(population, average_across_county, alpha=0.5, c="#7A68A6")
plt.scatter([population[i_min], population[i_max]],
            [average_across_county[i_min], average_across_county[i_max]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="extreme heights")

plt.xlim(100, 1500)
plt.title("Average height vs. County Population")
plt.xlabel("County Population")
plt.ylabel("Average height in county")
plt.plot([100, 1500], [150, 150], color="k", label="true expected \
height", ls="--")
plt.legend(scatterpoints=1)

plt.show()

从上图中可以观察到什么现象？如果不考虑人口数量大小，推断出的结果就可能存在错误。

例如，如果不考虑人口数大小，我们可能会说有最高身高人口的县也有最低身高的人口，反之亦然（不考虑人口数量，图中的 x 轴就变成一个点了，红色圆圈的点就会在一条竖线上了，所以就反映出最高人口的县也会有最低人口，最低人口的县也会有最高人口)。

通过这个例子，我们似乎得出了一个结论，小数据集情况下大数定理失效了！但是这种推断结果真的是对的吗？

从贝叶斯不确定性理论的角度出发，准确的答案应该是：不确定！

理论上来说，确实存在可能在说图中红圈所在的 x 轴对应的县确实存在极值的身高，但是在这种小数据情况下，不确定性是非常高的，高到我们几乎不能将其采纳作为一种推断结果。换句话说，这种高风险的推断对我们来说是没有意义的。

从大数定理的角度来说，这个县不一定有最极端高度的人存在。错误的原因是人口数量少时计算的均值并不能反映出人口真实的期望值(真实是 μ = 150)。样本的数量 N 越小，大数定理就会失效的越厉害。

2. 实例2 - Kaggle的美国人口普查反馈比例预测

下面是 2010 年美国人口普查数据，它不是按照县对数据分区，而是按照街区进行分区。

本例中的数据来自于 Kaggle 机器学习竞赛中使用的数据集。竞赛的目的是为了预测各个街区对人口普查邮件的回复率，回复率的取值范围为 0 到 100，使用的预测特征包括:人口平均收入，女性数量，停车场的个数，小孩的平均个数等。

下面画出人口调查邮件的回复率和街区中人口数量的关系图（多维自变量不方便在一张图上展示）

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 6.5)
data = np.genfromtxt("./data/census_data.csv", skip_header=1,
                     delimiter=",")
plt.scatter(data[:, 1], data[:, 0], alpha=0.5, c="#7A68A6")
plt.title("Census mail-back rate vs Population")
plt.ylabel("Mail-back rate")
plt.xlabel("population of block-group")
plt.xlim(-100, 15e3)
plt.ylim(-5, 105)

i_min = np.argmin(data[:, 0])
i_max = np.argmax(data[:, 0])

plt.scatter([data[i_min, 1], data[i_max, 1]],
            [data[i_min, 0], data[i_max, 0]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="most extreme points")

plt.legend(scatterpoints=1)

plt.show()

上图反映了典型的统计现象。这里的典型指的是散点图的形状。散点图是一个典型的三 角形，随着样本数量的增加，这个三角形变得越紧实（因为越来越符合大数定理的要求）。

可以这么理解，随着 N 不断增大，后验概率越来越收敛到某一个固定值，而在 N 很小的时候，后验概率的摆动是很剧烈的。

3. 实例3 - 如何对Reddit网站上的评论进行排序

上面的两个例子都旗帜鲜明地表达了一个观点：大数据！big data！概率统计和机器学习只有在大数据时才能体现出威力，如果你没有大数据，那就别痴心妄想做机器学习项目了。

但是世上之事，不如意者十之，我们不可能永远在项目中能拿到大数据，就算你身在BAT、google这样的大数据公司，你所拥有的也不一定就是大数据。

参考Andrew Gelman的一句名言：样本从来都不是足够大的。如果 N 太大不足以进行足够精确的估计，你需要获得更多的数据。但当 N “足够大”，你可以开始通过划分数据研究更多的问题，例如在民意调查中，当你已经对全国的民意有了较好的估计，你可以开始分性别、地域、年龄进行更多的统计。N 从来都无法做到足够大，因为当它一旦大了，你总是可以开始研究下一个问题从而需要更多的数据。

回到我们这个小节的问题，考虑在线商品的评分问题，这在生活中很常见，各种商城类app里都会有针对某个商品的历史购买者评论。

你如何评价只有一个人，两个人或者三个人对一个商品的五星级评价结果。我们多少都会明白，这样少的评分结果无法真正的反映出商品的真实价值。这种评价信息不充分的情况给排序带来了困难。很多人都知道在线搜索图书， 视频时，如果用它们的评分或者评论对返回结果进行排序，最终用户得到的结果并不会很好。通常情况下搜索结果中排在最前面的视频或评论往往有很好的评分，而这些好的评分都是由粉丝评出的。还有一些真正高质量的视频或者评论排在了稍微靠后的位置， 它们的评分在 4.8 左右。我们如何对这种情况进行修正？

1）我们如何对评论进行排序

Reddit是一个流行的网站。它主要提供一些小说或者图片的链接，这个网站最流行的部分是它对每个链接的评论。Reddit可以对每个评论进行投票(被称为赞成票或否决票)。默认状态下Reddit将把最好的评论排在前面。

现在假设我们是Reddit后台的算法工程师，我们会如何设计这种排序算法呢？

需要明白的是这是一个非常复杂且牵一发而动全身的事情，一个简单的例子是消费者公众是具有从众效应的，即大多数时候我们会优先看大多数的意见，且很容易接受所谓的大多数人的意见。这种现象表现在很多方面，例如百度搜索的首页的访问量会远大于后面的页面，且这种差距一旦建立，会越来越大；一个商品的评论一旦上了淘宝的爆款热搜，它的购买热度会越来越大，因为公众会优先购买它们最先看到的所谓爆款。

而评论排序这个问题也是一样的，一旦一个评论被排到靠前的位置，用户会优先看到这些评论，而有更大的几率接受该评论的观点，从而形成从众效应，最终导致靠前的评论的点赞数会越来越多。

基于这样的复杂性原因，我们接下来对评论排序这件事展开多维度的讨论。

热度

认为点赞数越多的评论越好。但是一个问题在于，当一评论有几百个赞，但是被反对了上千次时，虽然热度很高，但是把这条评论当成最佳评论是有争议的。

但是这里的问题也非常复杂，一个观点有上百文赞同，有上千人反对，那么这个观点就一定是对，或者一定就是错的吗？都不一定，很多时候，这取决于下游的商业和市场策略。

差数

用点赞数减去反对数来打分。这解决了以上热度指标的问题。但是依然有缺陷，它不能解决评论的时间属性导致的问题。

因为一个链接在线上会存在很多，不断有用户来点赞或反对，更老的评论会随着时间积累更多的点赞数，但并不意味着他们是最佳评论。这个问题的根源显然是没有进行归一化。

但是反过来想，这个问题也十分复杂，为什么一定要归一化呢？试想，一个评论如果在一个相当长的时间区间内，持续不断地被其他用户所接受并点赞，这不就意味着这个评论的观点接受度很高吗？

时间调节

考虑用差数除以评论的寿命（相当于归一化）。这样得到了一种类似每秒差数或每分钟差数这样的速率值。

看起来足够合理了，但是同样存在缺陷，比如一个一秒前发布的评论，只需要一票，就能击败一百秒前的有用九十九票的评论。出现这种问题的原因是没有把持续时间的跨度考虑进去。

避免这一问题的一种解法是只针对 t 秒前的评论进行统计（过滤掉刚刚出现的新鲜的评论）。但是 t 怎么选呢？而且这样是不是意味着 t 秒内的评论都是不好的呢？

读者看出来了吗？我们是在对不稳定量和稳定量进行比较（新评论 PK 老评论）

好评率

对某一评论的好评数除以对该评论的好评和差评的总和，得到好评率， 作为排名权重。这解决了 3 中的问题，如果用好评率作为排名结果，即使一个新的评论也会和老的评论一样得到较高的排名得分。

但是依然存在缺陷，在这种度量标准下，如果对一个评论的好评个数只有一个，即好评率为 1.0，而另一个评论的好评个数为 999，差评个数为 1，此时这个评论的好评率为 0.999，显然后一个评论的质量要比前一个高，这“可能”不太客观。

我说“可能”是有理由的，前者虽然只有一个赞成票，但还是有可能由于有 999 个赞成票的后者，只是说这个可能性有多少罢了，我们不能武断地下结论。

2）真实的好评率如何得到？ - 贝叶斯推断来帮忙

我们真正目标就是估计真实的好评率。注意，真实的好评率并不是观测到的好频率， 真实的好评率隐藏起来了，我们能观测到的仅仅是每个评论的好评和差评。

当有 999 个 好评，1 个差评时，我们确实可以说真实的好评率接近于 1，我们能这样说的原因是因为有大数定理的存在；

反之，当只有一次好评时，我们无法断言真实的好评率就是 1，因为在小数据下，后验估计时不确定的，这听起来像个贝叶斯问题，我们显然无法得到一个对好评率的准确估计，我们只能对好评率可能的概率分布进行一个贝叶斯推断。

贝叶斯方法需要好评率的先验知识，获得先验知识的一个方法是基于好评率的历史分布。可以用好评率的历史分布作为定好评率的先验分布。只要对Reddit的评论进行挖掘 就可以得到好评率的历史分布。计算好评率的历史分布时还应注意以下问题：

1. 倾斜数据: 大部分的评论都没有投票数据，因此有许多评论的好评率都非常小 (极端值)，这就使分布大多数集中在极端值附件。可以为投票的个数设个阈值，大于 该阈值的情况才予以考虑。阈值的选择还要考虑其对好评率精度的影响。
2. 偏差数据: Reddit有许多不同的子页面，称为subreddits。例如有两个页面，第一个为:r/aww 页面粘贴许多可爱动物的照片，第二个为:r/politics 讨论政治问题。用户对这两个页面的评论行为是不 相同的。浏览第一个页面的人的性格会更友好，更有爱心，也就会有更多的好评(与第二个页面相比)，而第二个页面的评论可能会有更多的争论或不同的观点出现。因此并不是所有的评论都是一致的。

鉴于此，我觉得用均匀分布作为评论好评率的先验分布是最佳的做法。

假设共有𝑁个投票数，好评率为𝑝，好评的个数就是二项随机变量（好评或差评），参数为 𝑝 和𝑁。 构造一个函数对 𝑝 进行贝叶斯推断。

import pymc as pm


def posterior_upvote_ratio(upvotes, downvotes, samples=20000):
    """
    This function accepts the number of upvotes and downvotes a particular submission received, 
    and the number of posterior samples to return to the user. Assumes a uniform prior.
    """
    N = upvotes + downvotes
    upvote_ratio = pm.Uniform("upvote_ratio", 0, 1)
    observations = pm.Binomial("obs", N, upvote_ratio, value=upvotes, observed=True)
    # do the fitting; first do a MAP as it is cheap and useful.
    map_ = pm.MAP([upvote_ratio, observations]).fit()
    mcmc = pm.MCMC([upvote_ratio, observations])
    mcmc.sample(samples, samples / 4)
    return mcmc.trace("upvote_ratio")[:]

下面是后验分布的结果：

figsize(11., 8)
posteriors = []
colours = ["#348ABD", "#A60628", "#7A68A6", "#467821", "#CF4457"]
for i in range(len(submissions)):
    j = submissions[i]
    posteriors.append(posterior_upvote_ratio(votes[j, 0], votes[j, 1]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.9,
             histtype="step", color=colours[i % 5], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )

plt.legend(loc="upper left")
plt.xlim(0, 1)
plt.title("Posterior distributions of upvote ratios on different submissions")

从上图中可以看出，某些分布很窄，其他分布表现为长尾。这体现了我们对真实好评率的不确定性。

3）如何针对评论的后验分布的进行排序

回到我们这个小节的主题，我们的目标是：如何对评论进行从好到坏的排序？

当然，我们是无法对分布进行排序的，排序只能是标量值。有很多种方法能够从分布中提取出标量值，用期望或均值来表示分布就是一个方法，但是均值并不是一个好办法，因为它没有考虑到分布的不确定性。

我们这里建议使用 95%最小可信值，定义为真实参数只有 5% 的可能性低于该值（即 95% 置信区间的下界），下图是根据 95%最小可信值画的后验分布。

N = posteriors[0].shape[0]
lower_limits = []

for i in range(len(submissions)):
    j = submissions[i]
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.9,
             histtype="step", color=colours[i], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )
    v = np.sort(posteriors[i])[int(0.05 * N)]
    # plt.vlines( v, 0, 15 , color = "k", alpha = 1, linewidths=3 )
    plt.vlines(v, 0, 10, color=colours[i], linestyles="--", linewidths=3)
    lower_limits.append(v)
    plt.legend(loc="upper left")

plt.legend(loc="upper left")
plt.title("Posterior distributions of upvote ratios on different submissions");
order = np.argsort(-np.array(lower_limits))
print(order, lower_limits)

为何基于这个量进行排序是个好主意？

因为基于 95% 最小可信值进行排序，是对最佳评论的一个最为保守的估计。即便在最差情况下，也就是说我们明显高估了好评率时，也能确保最佳的评论排在顶端。在这一排序思路中，我们利用了以下很自然的特性：

1. 给定两个好评率相同的评论，我们会选择票数最多的作为更加评论，因为票数更多的评论在贝叶斯后验概率分布中更集中，落在 95%置信区间的概率也更大，也即我们更确定其好评率更好。
2. 给定两个票数一样的评论，我们会选择好评数更多的。

 

10. 大数定律在工程领域的具体应用 - 保险行业的赚钱原理

0x1：保险行业的业务背景

保险是对风险的保障，它提供一份帮助人们分散危险、分摊损失的机制，这就是保险的本质。

本质是投保人通过不同数额的保费，将不同类型的风险打包转卖给了保险公司，所以保险费是建立保险基金的来源。

0x2：保费制定的理论支撑 - 大数定律

大数定律是保险业保险费计算的科学理论基础，当投保人（保单）的数量足够大时，由切比雪夫大数定律知，被保人缴纳的纯保费与其能获得的赔款的期望值是相同的。这个结论反过来，也可以说明保险公司应该如何合理的收取保费。

假设有n个投保人购买了n份相互独立的保险，每个人出险的概率为p，每个人可能获得的赔偿金为a，应缴纳的保费为b。

当n足够大时，设购买该险种的投保人出险的人数为X，由已知得 X~B(n，p)。则：

nb = a * E(X) = a * np

=> b = ap

即保费的制定是每个人可能获得赔偿金的期望。

 

11. 我们怎样更好应用大数定理

尽管大数定理非常酷，但是正如其名，它只适用于足够大的数据量。我们已经看到，如果不考虑数据的构造，那么估计结果会受到很大影响。

1. 通过（低成本地）获得大量后验样本，可以确保大数定理适用于期望的近似估计
2. 在样本数量很小的贝叶斯推断中，后验分布包括更多的随机性。从后验分布中将能看出这些随机性的存在。后验分布越平坦随机性越大，越陡峭随机性越小。因此推断的结果是需要调校的。
3. 不考虑样本量会带来很大的影响，此时排序的依据往往是不稳定的，并会导致非常病态的排序结果。