By Long Luo
在上一篇文章 正态分布(Normal Distribution)公式为什么长这样？ 中，我们使用了投掷飞镖的模型，推导出了正态分布( \(\text{Normal Distribution}\) )的表达式。这种方法既优雅又直观，所以常被用于科普视频或者文章中。那么这个例子是怎么来的呢？我们知道这个方法是天文学家赫歇尔（ \(\text{John Herschel}\) ）在 1850 年给出的，难道他在投掷飞镖时想到的吗？
答案是否定的，原因是因为赫歇尔作为一个天文学家，需要精确的测量天体的位置，而在观测星星时，必须要考虑误差的影响。星星在天球中的位置误差是二维的，考虑到误差大家不太好理解，所以用了投掷飞镖这个更通俗易懂的例子。
正如法国著名哲学家孔德( \(\text{Auguste Comte}\) ，1798-1857)所说“To understand a science, it is necessary to know its history. ”，只有了解这个学科的发展历史，了解这个学科的重要概念是如何建立起来的，才能真正理解这个学科。不同于我们在课本中学习顺序，科学是用来解决实际问题的，科学是由一个个问题所驱动发展的。正如仅次牛顿和爱因斯坦的伟大物理学家麦克斯韦( \(\text{James Clerk Maxwell}\) ) 曾说过“It is of great advantage to the student of any subject to read the original memoirs on that subject, for science is always most completely assimilated when it is in the nascent state…”，我们学习历史上科学家是如何解决这些问题，用了什么方法，才能获取某个概念的 insight ，建立 intuition 。
正态分布，又被称为高斯分布（ \(\text{Gaussian Distribution}\) ），人们可能会以为正态分布是由高斯发现的，但事实并非如此！
正态分布最早是由法国数学家棣莫弗( \(\text{Abraham de Moivre}\) , 1667-1754)在 1718 年左右发现的。他为了解决朋友提出的一个赌博问题，而去认真研究了二项分布。他发现当实验次数增大时，二项分布( \(p=0.5\) )趋近于一个看起来呈钟形的曲线，如下图 1 所示。后来著名法国数学家拉普拉斯( \(\text{Pierre-Simon Laplace}\) , 1749-1827)对此作了更详细的研究，并证明了 \(p \ne 0.5\) 时二项分布的极限也是正态分布。之后人们便将此称为棣莫弗 - 拉普拉斯中心极限定理（ \(\text{Central limit theorem}\) ）。
16 和 17 世纪是天文学发展的黄金时期，这一时期的科学革命彻底改变了人类对宇宙的理解。哥白尼( \(\text{Nicolaus Copernicus}\) ，1473-1543)的日心说、开普勒( \(\text{Johannes Kepler}\) ，1571-1630)的行星运动三定律、伽利略( \(\text{Galileo Galilei}\) ，1564-1642)的望远镜观测以及牛顿( \(\text{Isaac Newton}\) ，1643-1727)的万有引力定律共同构成了现代天文学的基础。这一时期的科学家们不仅改变了人类对宇宙的理解，还为后续的科学研究提供了重要的方法和工具。
1766 年德国物理学家提丢斯（ \(\text{Johann Daniel Titius}\) ）把已知太阳系中行星的数据列成如下一张表，如下表所示：
| Planet | Actual | Formula |
|---|---|---|
| Mercury | 0.387 | 0.4 |
| Venus | 0.723 | 0.7 |
| Earth | 1.00 | 1.0 |
| Mars | 1.52 | 1.6 |
| ??? | 2.8 | |
| Jupiter | 5.20 | 5.2 |
| Saturn | 9.55 | 10.0 |
| Uranus | 19.2 | 19.6 |
| Neptune | 30.1 | 38.8 |
他发现行星与太阳之间的距离大致满足以下经验公式：
\[ a(n) = 0.3 \times 2^n + 0.4 \text{AU}, \, n = -\infty , 0, 1, 2, 3, 4, \dots \]
其中 \(a_n\) 是第 \(n\) 颗行星离太阳的平均距离（以天文单位 \(\text{AU}\) 表示）， \(n\) 是行星的顺序编号，但水星的 \(n\) 值为 \(-\infty\) 。
之后提丢斯写信联系了德国天文学家博得（\(\text{Johann Elert Bode}\)），告知博得自己这一发现。博得看到之后，于 1772 年进一步推广和公开发表了这一公式。
从表格中我们也可以看到，火星和木星的轨道之间距离大得不正常。根据公式，在 \(n = 5\) ，距离太阳 \(2.8\) \(\text{AU}\) 处应该存在一颗行星，于是博得呼吁天文学家都来寻找这颗行星。
最终， 1801 年 1 月，意大利天文学家朱塞普·皮亚齐( \(\text{Giuseppe Piazzi}\) , 1746-1826)发现了第一颗小行星谷神星( \(\text{Ceres}\) ) ，其距离与公式预测的 \(a_5 = 2.8\) 天文单位非常接近，其轨道如下图 2 所示：
皮亚齐发现谷神星之后，对其进行观测了 40 天，之后谷神星就进入了太阳的耀眼的光芒之中，受当时观测技术所限而无法继续追踪观测，这样谷神星就“丢失”了！因为观测数据太少，天文学家无法从收集到的少量数据中推断出它的位置，因为要用不到总轨道 \(1 \%\) 的数据去求解开普勒的椭圆轨道复杂非线性方程，这是一项艰巨的数学挑战。
但是好不容易发现一颗新行星，这么重大的天文发现，科学家们怎么能放过这样的机会呢！于是由 24 位当时知名天文学家组成了一个“失踪星球探测协会”，他们被亲切地称为“天体警察”，因为他们负责去“逮捕”失踪的谷神星。同时也向欧洲学术界上百名数学家求助，希望他们能够提供谷神星的轨道线索。当时法国著名数学家拉普拉斯就认为，用这么少的数据是不可能解决这个问题的。
这时，天才高斯出场了，当时他只有 24 岁，但他已经在数学领域展现出非凡的才华。高斯 17 岁时发现了质数分布定理和最小二乘法，19 岁时，证明了正十七边形可以用尺规作图，解决了这一流传 2000 年的数学难题。在看到这个问题之后，高斯只用了开普勒的行星运动三定律，和他新发现的误差分布规律和最小二乘法 ( \(\text{Least Squares Method}\) ) ，仅使用了下表中的 \(3\) 个数据就计算出了谷神星的轨道。
| right ascension | declination | Time | |
|---|---|---|---|
| Jan. 2 | \(51^{\circ} 47^{\prime} 49^{\prime \prime}\) | \(15^{\circ} 41^{\prime} 5^{\prime \prime}\) | \(8 \mathrm{h} 39\mathrm{m} 4.6\mathrm{s}\) |
| Jan. 22 | \(51^{\circ} 42^{\prime} 21^{\prime \prime}\) | \(17^{\circ} 3^{\prime} 18^{\prime \prime}\) | \(7 \mathrm{h} 20 \mathrm{m} 21.7 \mathrm{s}\) |
| Feb. 11 | \(54^{\circ} 10^{\prime} 23^{\prime \prime}\) | \(18^{\circ} 47^{\prime} 59^{\prime \prime}\) | \(6 \mathrm{h} 11 \mathrm{m} 58.2 \mathrm{s}\) |
他的最终计算是指向了一个完全不同的天空区域，之前被其他科学家所忽视的区域。从上图 2 也可以看出谷神星并不在黄道平面上，而是和黄道平面有个较大的夹角，所以谷神星才不容易寻找。
1801 年 12 月 31 日夜，德国天文爱好者奥伯斯( \(\text{Heinrich Olbers}\) , 1758-1840)，在高斯预言的时间里，用望远镜对准了这片天空。果然不出所料，谷神星出现了！
高斯计算谷神星轨道的方法比较复杂，可以写一篇上百页的论文，这里提供一篇论文 How Gauss Determined The Orbit of Ceres 供大家参考。
不同于其他科学家把皮亚齐的观测数据当作真实数据，高斯一开始就考虑到皮亚齐的观测数据误差问题，毕竟当时的观测设备精度有限，所以求得的轨道数据应该是一个范围。如下图 3 所示就是高斯计算出谷神星轨道的原始手稿：
高斯使用的数据分析方法，就是以正态误差分布为基础的最小二乘法，如下图 4 所示：
高斯认为这些数据点使用最小二乘法计算出的回归曲线，如下图 5 所示，等价于这些点沿着这条回归曲线呈正态分布。
那 么高斯是如何推导出误差是正态分布的呢？我们将在下一节进行推导。
天文学是人类历史上最古老的科学之一，其起源可以追溯到史前时代。古代文明通过观察天空中的天体运动来制定历法、导航和进行农业活动。天文学家不光对天体进行了长期的系统观测，记录了大量的天文现象，而且试图预测不同行星和恒星的轨迹和位置。
由于观察者、仪器和许多其他因素，天文学家的测量数据存在误差。即使对于同一个物体，他们也有不同的观察结果。这里误差指的是观测值与实际值之间的差异。千百年以来，对于有误差的测量数据，多次测量取算术平均是比较好的处理方法。虽然这种方法缺乏理论上的论证，也不断的受到一些人的质疑，但取算术平均作为一种异常直观的方式，在多年积累的数据的处理经验中也得到相当程度的验证，被认为是一种良好的数据处理方法。
伽利略在 1632 年的著作《关于托勒密和哥白尼两个主要世界体系的对话》中，非正式地讨论了他所谓的“观测误差”，也就是今天所谓的“随机误差分布”，对误差的分布做过一些定性的描述，主要是以下 2 个假设：
为了系统地分析误差，人们希望找到误差分布曲线来描述误差。因为误差显然是连续的，我们希望找到误差的概率密度函数( \(\text{Probability Density Function}\) ) 。那么误差的 PDF 是什么样的呢？
粗看这 2 个假设，看起来和赫歇尔（ \(\text{John Herschel}\) ）的 2 个假设差不了太多，那么仅凭这 2 个假设，能推导出误差是正态分布吗？
答案是否定的！
因为满足上述假设的曲线有很多，满足上述 2 个条件的函数有很多个，如下图 6 就展示了其中的几种可能函数形状：
历史上很多科学家尝试寻找误差分布曲线，比如辛普森( Thomas Simpson , 1710-1761) 、拉普拉斯等，但他们找到的误差分布曲线都不是正态分布。而高斯最终找到了误差概率密度函数是正态分布，如下公式所示：
\[ \varphi \Delta = \frac{h}{\sqrt {\pi}} \, e^{-hh \Delta \Delta} \]
其中 \(\Delta\) 是误差，也就是 \(x - \mu\) ， \(h\) 是一个常量，也就是观测精度的度量。
下面就让我们跟上高斯的思路去推导下，看看高斯是如何推导出误差呈正态分布的呢？
我们使用大写字母 \(X, Y, Z\) 表示随机变量，小写字母 \(x, y, z\) 表示随机变量具体的值。
设观测对象的真实位置是 \(m\) ，进行多次重复测量之后，我们得到测量数据 \(X_1, X_2, \cdots, X_{n-1}, X_n\) 。这些测量值都存在误差，用希腊字母 \(\epsilon\) 来表示误差， \(\epsilon_i = X_i - m\) 。我们认为这些误差满足一个未知的概率密度函数 PDF ，设为： \(f_{\epsilon}(\varepsilon)\) 。
基于上面提到的 2 个假设，我们已经知道了 PDF 的一些属性。
假设 1 告诉我们 \(f_{\epsilon}\) 是偶函数，如下公式所示：
\[ f_{\epsilon}(\varepsilon ) = f_{\epsilon }(-\varepsilon ) \]
假设 2 则告诉我们，小误差比大误差更容易发生，则 \(\varepsilon\) 越接近 \(0\) , 发生的可能性就越大。
现在我们有了一系列重复的测量值 \(x_1, x_2, \cdots , x_{n-1}, x_n\) ，同时也有对应的一系列误差值 \(\varepsilon_i = x_i - m\) 。我们可以使用这组数据构造似然函数 ( \(\text{Likelihood function}\) )，如下所示：
\[ \mathcal{L} = \prod_{i=1}^n f_{\epsilon}(\varepsilon_i) \]
由于对数函数是单调递增的，在极大化求解时更方便，所以上式的对数似然函数是：
\[ \ell = \sum_{i=1}^{n} \ln f_{\epsilon}(\varepsilon_i) \]
我们上述对数似然函数 \(\ell\) 取最大值，则需要 \(\ell^{\prime} = 0\) ，则有：
\[ \ell^{\prime} = \frac{\mathrm{d}\ell}{\mathrm{d} \varepsilon} = 0 \]
使用链式法则，可得：
\[ \left[ \ln f(x) \right ]^{\prime} = \frac{f^{\prime}(x)}{f(x)} \]
对对数似然函数求导可得：
\[ \ell^{\prime} = \left[ \sum_{i=1}^n \ln f_{\epsilon}(\varepsilon_i) \right]^{\prime} = \sum_{i=1}^{n} \frac{f_{\epsilon}^{\prime}(\varepsilon_i)}{f_{\epsilon}(\varepsilon_i)} \]
为了方便后续计算，令 \(g_{\epsilon}(\varepsilon) = \dfrac{f_{\epsilon}^{\prime}(\varepsilon)}{f_{\epsilon}(\varepsilon)}\) ，则 \(\ell^{\prime}\) 可以重写为：
\[ \begin{aligned} \ell^\prime & = \sum_{i=1}^{n}g_{\epsilon}(\varepsilon_i) \\ & = g_{\epsilon}(\varepsilon_1) + g_{\epsilon}(\varepsilon_2) + \cdots + g_{\epsilon}(\varepsilon_{n-1}) + g_{\epsilon}(\varepsilon_n) \\ & = g_{\epsilon}(x_1 - m) + g_{\epsilon}(x_2 - m) + \cdots + g_{\epsilon}(x_{n-1} - m) + g_{\epsilon}(x_n - m) \end{aligned} \]
现在我们求极大似然估计( \(\text{Maximum Likelihood Estimation}\) ) ，也就是需要求 \(\ell^\prime = 0\) 的解。这种情况下，要么我们知道 PDF ，要么知道 \(m\) ，这样才可以求解 \(\ell^\prime = 0\) 。但是目前我们既不知道 PDF 也不知道 \(m\) ，该如何进行下一步工作呢？
到这里似乎已经陷入死胡同了，但是高斯假设极大似然估计 \(\ell^\prime = 0\) 的解就是算术平均，也就是说上述公式的解 \(m\) 是：
\[ \hat{m} = \frac{1}{n}\sum_{i=1}^n x_i \]
这里高斯把整个问题的思考模式倒过来：既然千百年来大家都认为算术平均是一个好的估计，那我就认为极大似然估计推导出的就应该是算术平均！
这确实只是一种猜测，并不是严谨的证明。但著名物理学家马赫说：“科学的功能是代替经验。这样，科学一方面必须依然停留在经验的范围之内，另一方面必须加速超越经验”。优秀的科学家往往能够通过敏锐的直觉从海量的可能性中找到正确的方向，再通过实验观测与数学推理进行验证，从而得到其中运行的规律。
令 \(\ell^\prime = 0\) ，用测量数据的算数平均值 \(\bar{x}\) 代入 \(m\) ，可得：
\[ g_{\epsilon}(x_1 - \bar{x}) + g_{\epsilon}(x_2 - \bar{x}) + \cdots + g_{\epsilon}(x_{n-1} - \bar{x}) + g_{\epsilon}(x_n - \bar{x}) = 0 \tag{3} \]
上述公式显然对于任意 \(n\) 个测量数据都成立，因此不妨设 \(x_1 = m\) ， \(x_2 = x_3 = x_4 = \cdots = x_{n-1} = x_n = m - nN\) ，这里 \(N\) 是一个常数，因此我们有：
\[ \begin{aligned} \bar{x} & = \frac{x_1 + x_2 + \cdots + x_{n-1} + x_n}{n} \\ & = \frac{m + (m - nN) + (m - nN) + \cdots + (m - nN)}{n} \\ & = \frac{m + (n-1) \cdot (m - nN) }{n} \\ & = m - (n - 1)N \end{aligned} \]
将上述值代入公式 \((3)\) 中，则有：
\[ \begin{aligned} g_{\epsilon}[m - (m-(n-1)N)] + g_{\epsilon}[(m - nN) - (m-(n-1)N)] + \cdots \\ + g_{\epsilon}[(m - nN) - (m-(n-1)N)] + g_{\epsilon}[(m - nN) - (m-(n-1)N)] = 0 \end{aligned} \]
化简可得：
\[ g_{\epsilon}[(n - 1)N] + g_{\epsilon}(-N) + \cdots + g_{\epsilon}(-N) + g_{\epsilon}(-N) =0 \]
注意到 \(g_{\epsilon}(x_2 - \bar{x}) = g_{\epsilon}(x_3 - \bar{x}) = \cdots = g_{\epsilon}(x_n - \bar{x}) = g_{\epsilon}(-N)\) ，共有 \(n-1\) 个，进一步化简可得：
\[ g_{\epsilon}[(n - 1)N] + (n - 1)g_{\epsilon}(-N) = 0 \tag{4} \]
根据假设 1 我们知道 \(f_{\epsilon}(\varepsilon) = f_{\epsilon}(-\varepsilon)\) ，则有：
\[ f_{\epsilon}^{\prime}(\varepsilon) = -f_{\epsilon}^{\prime}(-\varepsilon) \]
上式等号左边除以 \(f_{\epsilon}(\varepsilon)\) ，等式右边除以 \(f_{\epsilon}(-\varepsilon)\) ，易得：
\[ \frac{f_{\epsilon}^{\prime}(\varepsilon)}{f_{\epsilon}(\varepsilon)} = - \frac{f_{\epsilon}^{\prime}(-\varepsilon)}{f_{\epsilon}(-\varepsilon)} \]
也就是说 \(g_{\epsilon}(\varepsilon) = -g_{\epsilon}(-\varepsilon)\) ，即 \(g_{\epsilon}(-N) = -g_{\epsilon}(N)\) ，代入公式 \((4)\) ，可得：
\[ g_{\epsilon}[(n - 1)N] - (n - 1)g_{\epsilon}(N) = 0 \]
上式移项可得：
\[ g_{\epsilon}[(n - 1)N] = (n - 1)g_{\epsilon}(N) \]
上式是一个 \(f(k \cdot x) =k \cdot f(x)\) 形式的函数方程，要满足此函数方程，意味着：
\[ \frac{\Delta f(x)}{\Delta x} = C \]
这里 \(C\) 是一个常数。
故 \(g_{\epsilon}\) 是一个线性表达式：
\[ g_{\epsilon}(\varepsilon) = k \cdot \varepsilon \]
这里 \(k\) 是某个常数。
因为 \(g_{\epsilon}(\varepsilon) = \dfrac{f_{\epsilon}^{\prime}(\varepsilon)}{f_{\epsilon}(\varepsilon)}\) ，则有：
\[ \frac{f_{\epsilon}^{\prime}(\varepsilon)}{f_{\epsilon}(\varepsilon)} = k \cdot \varepsilon \tag{5} \]
对上式公式 \((5)\) 两边对 \(\varepsilon\) 进行积分，可得：
\[ \int \frac{f_{\epsilon}^{\prime}(\varepsilon)}{f_{\epsilon}(\varepsilon)}\, \mathrm{d} \varepsilon = \int k \cdot \varepsilon \mathrm{d}\varepsilon \]
可得：
\[ \ln f_{\epsilon}(\varepsilon) = \frac{k}{2} \, \varepsilon^2 + c \]
这里 \(k, c\) 均是常数，根据积分法则，最终我们得到了误差的 PDF 表达式：
\[ f_{\epsilon}(\varepsilon) = e^{\frac{k}{2}\varepsilon^2 + c} = e^c \cdot e^{\frac{k}{2}\varepsilon^2} = Ae^{\frac{k}{2}\varepsilon^2} \]
下一步就是求解 \(A, k\) 这 2 个常数。
根据假设 2 告诉我们 \(\varepsilon\) 越小，函数值越大，那么指数部分必然小于 0 。设 \(\dfrac{k}{2} = -h^2\) ， \(h\) 是某个常数，则得到了误差的概率密度函数：
\[ f_{\epsilon}(\varepsilon) = Ae^{-h^2 \varepsilon^2} \tag{6} \]
根据 PDF 满足积分为 1 ，则有：
\[ \int_{-\infty}^{+\infty} f_{\epsilon}(\varepsilon) \mathrm{d} \varepsilon = 1 \]
代入公式有：
\[ A \int_{-\infty}^{+\infty} e^{-h^2 \varepsilon^2} \mathrm{d} \varepsilon = 1 \]
上式中又出现了高斯积分，我们知道高斯积分结果为：
\[ \int_{-\infty}^{+\infty} e^{-x^2} \mathrm{d}x = \sqrt{\pi} \]
非标准高斯积分结果为：
\[ \int_{-\infty}^{+\infty} e^{-a^2x^2} \mathrm{d}x = \frac{\sqrt{\pi}}{a} \]
因此我们最终得到 \(A = \dfrac{h}{\sqrt{\pi}}\) ，于是我们也就推导出了误差分布曲线，也就是概率密度函数的形式为：
\[ \varphi (\epsilon ) = \frac{h}{\sqrt{\pi}}\,e^{-h^2\varepsilon^2} \]
高斯通过对观测误差的分析，寻找误差的概率分布密度函数。大胆假设最大似然估计的结果是算数平均，而在所有的概率密度函数中，唯一满足这个性质的就是正态分布。
此内容由惯性聚合(RSS阅读器)自动聚合整理，仅供阅读参考。 原文来自 — 版权归原作者所有。