机器学习项目中不可忽视的一个密辛 – 大数定理、中心极限定理[通俗易懂]

1. 前言

0x1：发展历史

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

在概率论发展的初期，人们主要还是以实验归纳为主，人们发现，当大量重复某一相同的实验的时候，其最后的实验结果可能会稳定在某一数值附近。就像抛硬币一样，当我们不断地抛，抛个上千次，甚至上万次，我们会发现，正面或者反面向上的次数都会接近一半。

除了抛硬币，现实中还有许许多多这样的例子，像掷骰子，最著名的实验就是泊松抛针实验。这些实验都像我们传达了一个共同的信息，那就是大量重复实验最终的结果都会比较稳定。那稳定性到底是什么？怎样去用数学语言把它表达出来？这其中会不会有某种规律性？是必然的还是偶然的？

这一系列问题其实就是大数定律要研究的问题。很早的时候，人们其实就发现了这一规律性现象，也有不少的数学家对这一现象进行了研究，这其中就包括伯努利（后来人们为了纪念他，都认为他是第一个研究这一问题的人，其实在他之前也早有数学家研究过）。

伯努利在1713年在他的《猜度数》中提出了一个极限定理，当时这个定理还没有名称，后来人们称这个定理为伯努利大数定律。因此概率论历史上第一个有关大数定律的极限定理是属于伯努利的，它是概率论和数理统计学的基本定律，属于弱大数定律的范畴。

当大量重复某一实验时，最后的频率无限接近事件概率。而伯努利成功地通过数学语言将现实生活中这种现象表达出来，赋予其确切的数学含义。他让人们对于这一类问题有了新的认识，有了更深刻的理解，为后来的人们研究大数定律问题指明了方向，起到了引领作用，其为大数定律的发展奠定了基础。

除了伯努利之外，还有许许多多的数学家为大数定律的发展做出了重要的贡献，有的甚至花了毕生的心血，像德莫佛—拉普拉斯，李雅普诺夫，林德伯格，费勒，切比雪夫，辛钦等等。这些人对于大数定律乃至概率论的进步所起的作用都是不可估量的。
1733年，德莫佛—拉普拉斯经过推理证明，得出了二项分布的极限分布是正态分布的结论，后来他又在原来的基础上做了改进，证明了不止二项分布满足这个条件，其他任何分布都是可以的，为中心极限定理的发展做出了伟大的贡献。在这之后大数定律的发展出现了停滞。直到20世纪，李雅普诺夫又在拉普拉斯定理的基础上做了自己的创新，他得出了特征函数法，将大数定律的研究延伸到函数层面，这对中心极限定理的发展有着重要的意义。到1920年，数学家们开始探讨中心极限定理在什么条件下普遍成立，这才有了后来发表的林德伯格条件和费勒条件，这些成果对中心极限定理的发展都功不可没。
经过几百年的发展，大数定律体系已经很完善了，也出现了更多更广泛的大数定律，例如切比雪夫大数定律，辛钦大数定律，泊松大数定律，马尔科夫大数定律等等。正是这些数学家们的不断研究，大数定律才得以如此迅速发展，才得以完善。

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

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年）

最早的中心极限定理是由法国数学家棣莫弗（1677年-1754年）于1733年发现的，他使用正态分布取估计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：棣莫弗-拉普拉斯定理应用举例

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

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

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

0x1：数学公式定义

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 

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

重构TODO LIST：

频率和概率的关系：从极限角度出发，基于极限收敛的理论得到等价近似关系

大数定律的应用：1）算术平均法则；2）频率估计概率；3）数理统计中参数的点估计思想；4）

在大数n的情况下，所有概率分布函数都会收敛到正态分布这个中心极限的形式上。同时各个概率分布函数之间也存在互相推导的关系。

这也是为什么GMM从理论上可以拟合任意复杂的概率分布函数的原因。

讨论一下保险行业的赚钱原理

中心极限定理和大数定律的区别
（1）大数定律是说，n只要越来越大，把这n个独立同分布的数加起来去除以n得到的这个样本均值（也是一个随机变量）会依概率收敛到真值u，但是样本均值的分布是怎样的我们不知道。

（2）中心极限定理是说，n只要越来越大，这n个数的样本均值会趋近于正态分布，并且这个正态分布以u为均值，sigma^2/n为方差。

（3）综上所述，这两个定律都是在说样本均值性质。随着n增大，大数定律说样本均值几乎必然等于均值。中心极限定律说，它越来越趋近于正态分布，并且这个正态分布的方差越来越小。

直观上来讲，想到大数定律的时候，脑海里浮现的应该是一个样本；
想到中心极限定理的时候脑海里应该浮现出很多个样本。

1. 大数定理是什么？

大数定理是一个非常底层的基础性原理，大量的机器学习理论和算法实际上都建立在其之上。只是因为其太底层，也太合理的，合理到我们感觉这就是一个毫无争议的真理而不需要去单独讨论。我们本文和读者来一起从贝叶斯推断的角度来尝试对大数定理展开一些讨论。

0x1：大数定理基本数学定义

设?? 为服从特定概率分布的?个独立采样样本，按照大数定理理论，只要?的期望?[?]是有限的，则就有如下关系：

用文字来描述：即来自同一分布的一组随机变量，其均值收敛域该分布的期望。

以上等式只有在极限条件下成立，但是随着参加平均计算的样本越来越多，均值将越来越接近?的期望值?[?]。除了一些极少数的特殊情况外，大数定理对于大多数分布都是正确的。

0x2：泊松随机变量的收敛 – 大数定理普遍存在的一个例证

下面是三个不同的泊松随机变量在大数定理上的实例，用三个不同的曲线表示。 三个泊松分布的参数? = 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）通过公式方式来描述收敛速率

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

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

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

0x3：连接概率和频率之间的桥梁 – 大数定理

在开始这个小节讨论之前，我们先抛出一个公式：

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

在期望值和估计的概率之间存在一种隐性关系。

首先定义如下函数：

，这种0-1分布有非常好的特性，可以很清晰地表达出期望和概率的关系。

根据大数定理，如果有样本?? ，可以估计出事件 A 的概率，用?(?)表示

由上面只是函数的定义可知，只有事件发生时1?(?)才取 1 其它情况都取 0， 所以只要用试验中事件发生的总次数除以试验的次数就是概率（回想一下平时我们是怎 么样用频率近似概率的）

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

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

2. 小数据的无序性 – 当 N 很小时大数定理的表现如何？

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

0x1：实例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 轴就变成一个点了，红色圆圈的点就会在一条竖线上了，所以就反映出最高人口的县也会有最低人口，最低人口的县也会有最高人口)。

1. 在小数据集情况下大数定理失效了吗？

但是这种推断结果真的是对的吗？

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

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

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

0x2：实例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 很小的时候，后验概率的摆动是很剧烈的。

0x3：实例3 – 如何对Reddit网站上的评论进行排序

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

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

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

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

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

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

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

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

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

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

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

1）热度

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

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

2）差数

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

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

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

3）时间调节

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

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

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

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

4）好评率

对某一评论的好评数除以对该评论的好评和差评的总和，得到好评率， 作为排名权重。这解决了 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. 给定两个票数一样的评论，我们会选择好评数更多的。

3. 大数定理的一些更多思考

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

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