2024-05-11 15:34:38
By Long Luo
从之前的文章 正态分布(NormalDistribution)公式为什么长这样？ 和 从最小二乘法到正态分布：高斯是如何找到失踪的谷神星的？，我们使用了 \(2\)种不同的方法最终得到了如下公式 \((1)\)所示的误差的概率密度函数 ( \(\text{ProbabilityDensity Function}\) ) ：
\[f(x) = \mathrm{e}^{-cx^2}, \, c > 0 \tag{1} \label{1}\]
其函数图像如下图 1 所示的钟形曲线 ( \(\text{Bell Curve}\) ) ：
在概率论中，我们需要保证上图 1 中 \(f(x)\) 和 \(x\) 轴围成的面积是 \(1\) , 即：
\[\int_{- \infty}^{+ \infty} f(x) \mathrm{d}x = 1 \tag{2}\label{2}\]
最终我们得到了正态分布 ( \(\text{NormalDistribution}\) ) 的公式如下所示：
\[f(x) = {\frac {1}{\sigma {\sqrt {2 \pi }}}}\;e^{-{\frac {\left(x - \mu\right)^{2}}{2 \sigma ^{2}}}} \tag{3} \label{3}\]
上式中有一个 \(\pi\) ，用费曼( \(\text{Richard Feynman}\))的话来说，当我们看到一个公式中存在 \(\pi\) 时，我们都要问自己“Where is thecycle?”。我们知道公式 \(\eqref{3}\)中的归一化系数 \(\dfrac {1}{\sigma {\sqrt {2\pi }}}\) 是为了保证 \(f(x)\)下的面积为 \(1\) ，出现 \(\pi\) 是因为高斯积分 ( \(\text{Gaussian Integral}\) ) 的结果为 \(\sqrt{\pi}\) 。
那么什么是高斯积分呢？高斯积分和圆有什么关系呢？
高斯积分是数学王子高斯 ( \(\text{CarlFriedrich Gauss}\) ) 名字命名，和正态分布密切相关，对高斯函数 (\(\text{Gauss function}\) )，也就是正态分布的概率密度函数在整个实数域进行积分，即 \(\int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2}\,\mathrm{d}{x}\) ，如下图 2 所示：
初看 \(\int e^{-x^2} \mathrm{d}{x}\)，你可能会觉得用分部积分可以求解原函数吧！一顿操作猛如虎之后，发现根本就计算不出来。
因为高斯积分在初等函数范围内是不可积的，其原函数无法用初等函数来表示，如下公式所示：
\[F(x) = \int_{0}^{x} e^{-t^2} \mathrm{d}t \tag{4}\label{4}\]
好在早有天才数学家用天才的解法得到了高斯积分的值为 \(\sqrt{\pi}\) ，如下图 4 所示：
\[\int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2} \,\mathrm{d}{x} =\sqrt{\pi} \tag{5} \label{5}\]
高斯积分看起来和圆没有任何关系，但结果中却出现了\(\pi\) ，那圆在哪里呢？
求解高斯积分有很多种解法，均来自天才数学家的天才想法。这里我们将介绍\(3\)种方法，第一种是常见的直角坐标系转为极坐标系法，另外 \(2\) 种来自 3Blue1Brown 的视频 Why π is in thenormal distribution (beyond integral tricks)，都非常优雅且直观，下面就让我们沉浸在数学的海洋里吧，感受数学的美！
我们设 \(I\) 表示高斯积分的值：
\[I = \int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2}\,\mathrm{d}{x} \tag{6} \label{6}\]
因为定积分是一个数，与被积函数的自变量无关，简单换下元则有：
\[I = \int_{-\infty}^{+\infty} \mathrm{e} ^{-y^2}\,\mathrm{d}{y} \tag{7} \label{7}\]
二者相乘，则对高斯积分进行了升维，得到 \(I^2\) ：
\[\begin{aligned}I^2 & = \int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2} \,\mathrm{d}{x}\int_{-\infty}^{+\infty} \mathrm{e} ^{-y^2} \,\mathrm{d}{y} \\& = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{e}^{-(x^2 + y^2)} \,\mathrm{d}{x} \,\mathrm{d}{y}\end{aligned}\]
下面我们从直角坐标系 ( \(\text{Cartesiancoordinates}\) ) 转换为极坐标系 ( \(\text{Polar coordinates}\) ) ，则有：
\[\begin{cases}x = r \cos \theta \\y = r \sin \theta\end{cases}\]
在直角坐标系中，我们的积分区域是：
\[\begin{cases}-\infty < x < \infty \\-\infty < y < \infty\end{cases}\]
转换为极坐标系之后，积分区域变为：
\[\begin{cases}0< r < \infty \\0< \theta < 2 \pi\end{cases}\]
在直角坐标系中，面积微元 \(\mathrm{d}\sigma = \mathrm{d}x \mathrm{d}y\)，而在极坐标系中，其面积微元如下图 7 所示：
计算其面积：
\[\begin{aligned}\mathrm{d} \sigma & = \frac{1}{2}\left(r + \mathrm{d}r\right)^2\cdot \mathrm{d}\theta - \frac{1}{2}r^2\cdot \mathrm{d}\theta\\& = \frac{1}{2}\left[r^2 + 2r \mathrm{d}r + \left(\mathrm{d}r\right)^2-r^2\right] \mathrm{d}\theta \\& = \frac{1}{2}\left[2r \mathrm{d} r + \left(\mathrm{d} r \right)^2\right] \cdot \mathrm{d}\theta\end{aligned}\]
因为 \((\mathrm{d}r)^2\) 是 \(2r\mathrm{d}r\) 的高阶无穷小，故有：
\[\mathrm{d} \sigma = \frac{1}{2}(2r \mathrm{d}r) \mathrm{d}\theta = r\mathrm{d}r \mathrm{d}\theta\]
所以：
\[\mathrm{d}x \mathrm{d}y = r \mathrm{d}r\mathrm{d}\theta\]
在极坐标系中， \(r^2 = x^2 + y^2\)，则有：
\[\begin{aligned}I^2 & = \int_0^{2\pi} \int_0^{+\infty} \mathrm{e} ^{-r^2} r\,\mathrm{d}{r} \,\mathrm{d}{\theta} \\& = \int_{\theta = 0}^{2\pi} \mathrm{d}{\theta} \int_{r=0}^{\infty}\mathrm{e} ^{-r^2} r \,\mathrm{d}{r} \\& = 2 \pi \int_0^{+\infty} r \mathrm{e} ^{-r^2} \,\mathrm{d}{r}\end{aligned}\]
令 \(u = r^2\) ，则有：
\[\mathrm{d}u = 2r \mathrm{d}{r}\]
\[r \mathrm{d}r = \frac{1}{2}\mathrm{d}u\]
根据微积分基本定理易得：
\[\begin{aligned}I^2 & = \pi \int_0^{+\infty} \mathrm{e} ^{-u} \,\mathrm{d}{u} \\& = \pi \left (-\mathrm{e} ^{-u}) \right\rvert _0^{+\infty} \\& = \pi\end{aligned}\]
开方，最终我们求出了高斯积分的值：
\[I = \int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2} \,\mathrm{d}{x} =\sqrt{\pi}\]
转换为极坐标法是常用的求高斯积分的方法，尤其是对高斯积分升维更是神来之笔，那么升维之后的高斯函数是什么样的呢？
\(1\) 维的高斯函数为：
\[f_1(x) = \mathrm{e}^{-x^2}\]
那么 \(2\) 维的高斯函数则为：
\[f_2(x, y) = \mathrm{e}^{-x^2} \mathrm{e}^{-y^2} = \mathrm{e}^{- (x^2 +y^2)}\]
\(1\) 维的高斯函数表现为钟形曲线，\(2\) 维的高斯函数则是 \(1\) 维的高斯函数沿着 \(z\) 轴旋转，如下图 8 所示：
钟形曲线绕 \(z\) 轴旋转一周会形成\(3\) 维的钟形曲面，如下图 9 所示：
对于钟形曲面上的任意点 \(P(x, y)\)，其距离 z 轴的距离 \(r = x^2 + y^2\)， 如下图 10 所示：
很明显， \(3\)维钟形曲面满足径向对称性，到了这里圆出现了！
\(2\)维的高斯函数表达式也可以写成：
\[f_2(x, y) = \mathrm{e}^{-r^2}\]
求解 \(2\)维的高斯函数，很明显就是求解这个 \(3\)维钟形体积。考虑如下图 12 所示的同心圆筒，其周长为 \(2 \pi r\) ，高为 \(\mathrm{e}^{-r^2}\) ，厚度为 \(\mathrm{d}r\) ，则其体积为：
\[\mathrm{d}V = 2 \pi r \mathrm{e}^{-r^2}\mathrm{d}r\]
对其积分可得：
\[\begin{aligned}\mathbb{V} & = \int_0^{+\infty} \mathrm{d}V \\& = \int_0^{+\infty} 2 \pi r \mathrm{e}^{-r^2}\mathrm{d}r \\& = \pi \int_0^{+\infty} 2r \mathrm{e}^{-r^2}\mathrm{d}r \\& = \pi \left (-e^{-\infty^2} -(-e^0) \right ) \\& = \pi\end{aligned}\]
我们通过对高斯积分进行升维得到了一个 \(2\) 维高斯曲面，如果我们沿着平行 \(y\) 轴方向对其进行切片，如下图 13所示：
我们可以得到了一系列切片，而每个切片都是钟形曲线，如下图 14 所示就是\(y = 0\)时形成的钟形曲线，此时就是：
\[f_2(x, y) = \mathrm{e}^{-x^2} \mathrm{e}^{-y^2} = \mathrm{e}^{-x^2}\]
而所有的切片都可以理解为钟形曲线乘以某个比例常数进行放缩，如下图 15所示：
那么每个切片的面积就是：
\[S = \mathrm{e}^{-x^2} \int_{-\infty}^{+\infty} \mathrm{e} ^{-y^2}\,\mathrm{d}{y}\]
每个切片的体积微元则为：
\[\mathrm{d}V = \mathrm{e}^{-y^2} \int_{-\infty}^{+\infty} \mathrm{e}^{-x^2} \,\mathrm{d}{x}\]
对体积进行积分则可得：
\[\mathbb{V} = \int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2} \,\mathrm{d}{x}\int_{-\infty}^{+\infty} \mathrm{e} ^{-y^2} \,\mathrm{d}{y}\]
由于 \(x\) 和 \(y\) 轴的对称性，易得：
\[\mathbb{V} = \left ( \int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2}\,\mathrm{d}{x} \right )^2\]
而我们已经知道体积 \(\mathbb{V} =\pi\) ，所以：
\[\int_{-\infty}^{+\infty} \mathrm{e} ^{-x^2} \,\mathrm{d}{x} = \sqrt{\pi}\]
高斯积分是最美积分之一，其结果是出人意料的 \(\sqrt {\pi}\) ，而 \(\pi\)正是来自旋转对称性。通过之前的文章我们知道满足旋转对称性和各方向互相独立\(2\) 个性质可以推导出正态分布，在\(2\)维钟形曲线等概率事件是在同一等高线形成的圆上，这就是圆的由来。
2024-05-02 10:35:08
By Long Luo
在上一篇文章 正态分布(NormalDistribution)公式为什么长这样？中，我们使用了投掷飞镖的模型，推导出了正态分布( \(\text{Normal Distribution}\))的表达式。这种方法既优雅又直观，所以常被用于科普视频或者文章中。那么这个例子是怎么来的呢？我们知道这个方法是天文学家赫歇尔（\(\text{John Herschel}\) ）在 1850年给出的，难道他在投掷飞镖时想到的吗？
答案是否定的，原因是因为赫歇尔作为一个天文学家，需要精确的测量天体的位置，而在观测星星时，必须要考虑误差的影响。星星在天球中的位置误差是二维的，考虑到误差大家不太好理解，所以用了投掷飞镖这个更通俗易懂的例子。
正如法国著名哲学家孔德( \(\text{AugusteComte}\) ，1798-1857)所说“To understand a science, it isnecessary to know its history.”，只有了解这个学科的发展历史，了解这个学科的重要概念是如何建立起来的，才能真正理解这个学科。不同于我们在课本中学习顺序，科学是用来解决实际问题的，科学是由一个个问题所驱动发展的。正如仅次牛顿和爱因斯坦的伟大物理学家麦克斯韦(\(\text{James Clerk Maxwell}\) )曾说过“It is of great advantage to the student of any subject to readthe original memoirs on that subject, for science is always mostcompletely assimilated when it is in the nascentstate…”，我们学习历史上科学家是如何解决这些问题，用了什么方法，才能获取某个概念的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{JohannesKepler}\) ，1571-1630)的行星运动三定律、伽利略( \(\text{Galileo Galilei}\)，1564-1642)的望远镜观测以及牛顿( \(\text{Isaac Newton}\)，1643-1727)的万有引力定律共同构成了现代天文学的基础。这一时期的科学家们不仅改变了人类对宇宙的理解，还为后续的科学研究提供了重要的方法和工具。
1766 年德国物理学家提丢斯（ \(\text{JohannDaniel 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 ElertBode}\)），告知博得自己这一发现。博得看到之后，于 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)，在高斯预言的时间里，用望远镜对准了这片天空。果然不出所料，谷神星出现了！
高斯计算谷神星轨道的方法比较复杂，可以写一篇上百页的论文，这里提供一篇论文HowGauss 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 \cdotf(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 \cdote^{\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}\]
高斯通过对观测误差的分析，寻找误差的概率分布密度函数。大胆假设最大似然估计的结果是算数平均，而在所有的概率密度函数中，唯一满足这个性质的就是正态分布。
2024-04-27 13:16:57
By Long Luo
相信大家或多或少都听过六西格玛( \(\text{6Sigma}\) ) 1 这个词，六西格玛是指生产的产品中，\(99.99966\%\)的产品是没有质量问题的，即只有 \(3.4ppm\) 的不良率。
假如一家工厂生产某型号零件，零件的长度要求是 \(100mm\) ，允许的标准差是 \(0.1mm\) 。根据 \(6 \sigma\) 原则，零件规格允许的偏差范围是：\(100 \pm 6 \times 0.1 = 100 \pm 0.6\)。
这意味着，零件长度超过 \(100.6mm\)或低于 \(99.4mm\)的概率是非常低的，约为 \(0.00034\%\)。如果工厂每天生产 100 万个零件，只允许有 \(3.4\) 个零件会超出 \(6 \sigma\)的范围，几乎可以忽略不计。因此，生产过程是极其稳定和可靠的，达到了六西格玛水平。
那么 \(6 \sigma\) 中 \(3.4ppm\) 的不良率来自哪里呢？
学过中学数学都知道，在正态分布( \(\text{Normal Distribution}\) ) 2 中， \(68.27\%\) 的数据位于平均值的一个标准差内，\(95.45\%\) 位于两个标准差内， \(99.73\%\) 位于三个标准差内，这也是著名的68-95-99.7 Rule 3 ，如下图 1 所示：
数据可以用不同的方式“分布”，比如数据可以向左散布的多一些，也可以向右散布的多一些，或者分布的乱七八糟，如下图2 - 4 所示，
但数据经常会集中在一个中心值的附近，而不向左或右偏斜，像一个钟形，如下图5 所示。
正态分布，又称高斯分布（ \(\text{GaussianDistribution}\) ），是一种重要的概率分布，数学王子高斯 4在正态分布的研究和应用上做出了巨大贡献。有很多日常现象都符合这种分布，如人的身高、考试成绩等。正因为它几乎无处不在，所以叫\(\text{Normal Distribution}\)。德国曾经发行的一款 10 马克的纸币上就印着高斯和正态分布曲线，如下图 6所示。
这个曲线的数学公式大家在中学里都早已见过，如下所示：
\[f(x; \mu ,\sigma ) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac{\left(x - \mu \right)^{2}}{2 \sigma ^{2}}}} \tag{1} \label{1}\]
\(\mu\)是正态分布的数学期望值，可解释为位置参数，决定了分布的位置，表示曲线中心在哪里；方差\(\sigma ^{2}\)为尺度参数，决定了分布的幅度，表示曲线的扁平情况。均值\(\mu\) 和方差 \(\sigma ^{2}\)不同，曲线形状也有所不同，如下图 7 所示：
正态分布的公式看起来非常复杂，里面有 \(\pi\) 、 \(e\) 、 \(\mu\) 、 \(\sigma\)，组合起来非常复杂。在学习时，课本介绍正态分布时就直接给出这个公式，却从来不说明这个概率密度函数是怎么推导来的，来龙去脉是什么。最近看了3Blue1Brown 关于 概率论的系列视频，我知道了正态分布曲线公式为什么是这样，我们将在下一章节中推导出这个公式。
有很多种方法都可以推导出正态分布公式，这里将介绍一种既优雅又直观的推导方式，由天文学家赫歇尔（\(\text{John Herschel}\) ） 5 在 1850 年给出的。3Blue1Brown 的视频 Why π is in thenormal distribution (beyond integral tricks)中详细介绍了这种方式。不过视频中有一些不够严谨的地方，下面会先介绍视频中的推导方法，然后再介绍严谨的数学分析法。
考虑向一个镖盘投掷飞镖，过镖盘中心作 \(x\) 轴和 \(y\)轴。每次投掷飞镖都会受到随机因素而偏离目标，故每次飞镖在镖盘的落点 \((x, y)\) 是 \(2\) 维随机变量。
假设满足以下 \(2\) 个条件：
如下图 9 所示，设箭头所示落点区域 \(P\) 的概率密度函数为 \(f_2(x, y)\) ， \(f_2\) 表示有 \(2\) 个输入参数，落点的坐标 \(x\) 和 \(y\) 。
由条件 \(1\) 我们知道每个落点区域\(P\) 的概率密度，可以表示为 \(x\) 轴方向上概率密度函数与 \(y\)轴方向上概率密度函数的乘积，每个方向上的概率密度函数只有对应方向上一个参数。
设 \(x\) 轴方向上概率密度函数为\(g(x)\) ， \(y\) 轴方向上概率密度函数为 \(h(y)\) ，则有以下公式：
\[f_2(x, y) = g(x) h(y) \tag{2} \label{2}\]
又因为条件 \(2\)我们知道可以通过旋转对称性，可知 \(x\)轴和 \(y\)轴方向上概率密度函数相同，所以 \(g(y) =h(y)\) 。同时落点 \(P\)距离原点的距离为 \(r = \sqrt{x^2 +y^2}\) ，如下图 10 所示：
则 \(f_2(x, y)\)概率密度函数表示为半径（距原点的距离）的单变量函数 \(f(r)\) ，则有：
\[f_2(x, y) = f(r) = f(\sqrt {x^2 + y^2}) \tag{3} \label{3}\]
综合 \(\eqref{2}\) 和 \(\eqref{3}\) 可得：
\[f(\sqrt {x^2 + y^2}) = g(x) g(y) \tag{4} \label{4}\]
假设 \(x\) 轴上距离原点距离为 \(r\) 的点 \(P\) ，坐标为 \((r, 0)\) ，如下图 11 所示：
点 \(P(r, 0)\)处的概率密度函数可以写成：
\[f_2(r, 0) = f(r) = g(r) g(0) \tag{5} \label{5}\]
从上式 \(\eqref{5}\) 可知 的 \(f(r)\) 等于 \(g(r)\) 乘以某个常数 \(g(0)\) ，所以 \(f(r) = C g(r)\) ，其中 \(C\) 为某个常数。
由于最终都需要进行归一化，这里不妨设 \(C =1\) ，则公式 \(\eqref{4}\)可以写成：
\[f(\sqrt {x^2 + y^2}) = f(x) f(y) \tag{6} \label{6}\]
至此我们得到了最重要的关系式，问题转变为如何求解函数\(f\) 。公式 \(\eqref{6}\) 是一个函数方程 6，熟悉函数方程的同学可能一眼就知道满足公式 \(\eqref{6}\)的函数解是指数函数 7 。
\[f({r}) = e^{-{r}^2} \tag{7} \label{7}\]
如果我们不知道这个函数方程的解怎么办呢？下面我们就来求解下。
设函数 \(h(x) = f(\sqrt {x})\)，则有 \(h(x^2) = f(x)\) ，那么公式\(\eqref{6}\) 可以写成：
\[h(x^2 + y^2) = h(x^2) h(y^2) \tag{8} \label{8}\]
从公式 \(\eqref{8}\) 我们知道任取\(2\)个正整数，先相加再代入函数 \(h\) ，结果等于先分别代入函数 \(h\) 再相乘。这意味着，函数 \(h\)把加法变成了乘法。
实际上到了这里，还记得中学时学过的指数函数吗？
对于以某个常数 \(c\)为基数的指数函数 \(h(x) = c^x\)，可以将加法变成乘法。
但这个只是我们的猜测，我们还需要严谨的证明我们的猜测。
如果这个性质适用于任意 \(2\)个正整数，那么很容易扩展至任意 \(n\)个正整数：
\[h(x_1 + x_2 + \cdots + x_n) = h(x_1)h(x_2) \cdotsh(x_n) \tag{9} \label{9}\]
不妨代入 \(x = 5\) ，则有：
\[\begin{aligned}h(5) & = h(1 + 1 + 1 + 1 + 1) \\& = h(1)h(1)h(1)h(1)h(1) \\& = h(1)^5\end{aligned}\]
这里 \(x\) 可以为任意整数 \(n\) ，则有：
\[h(n) = h(1 + \cdots + 1) = h(1) \cdots h(1) =h(1)^n \tag{10} \label{10}\]
由公式 \(\eqref{6}\)，我们知道：
\[f(0) = f(0)f(0)\]
因为 \(f(0) \ne 0\) ，因此 \(f(0) = 1\) ，那么有：
\[h(0) = h(0)h(0) = h(0)^2\]
对于 \(p \in \mathbb{Z}^+\)，有：
\[1 = h(0) = h(p - p) = h(p)h(-p) = h(1)^p h(-p)\]
易得对于任意负整数 \(-p\) ，有：
\[h(-p) = h(p)^{-1} = h(1)^{-p}\]
这样我们就将函数 \(h\) 推广至整数域\(\mathbb{Z}\) 。
同理，对于任意正整数 \(q \in\mathbb{Z}^+\) ，可以得到以下表达式
\[h(q) = \prod_{k=1}^q h(1) = h(1)^q\]
对公式适当变形，可得：
\[h(1) = h(\frac {q}{q}) = h(\frac {1}{q})^q\]
可以得到：
\[h(\frac {1}{q}) = h(1)^{\frac{1}{q}}\]
综合上述可得，对于任意有理数 \(\dfrac{p}{q} \in \mathbb{Q}\)，都满足：
\[h(\frac{p}{q}) = h(1)^\frac{p}{q}\]
同时由于概率密度函数连续，那么对于无理数也成立。则根据上述推导我们知道公式\(\eqref{10}\) 在实数域 \(\mathbb{R}\) 都成立。
不妨令 \(a = h(1)\)，则我们得到了最终表达式：
\[h(x) = a^x \tag{11} \label{11}\]
为了方便和统一，任意指数函数都可以写成以 \(e\) 为底的指数函数：
\[h(x) = e^{\ln a x} = e^{c x} \tag{12} \label{12}\]
至此我们就求出了满足函数方程 \(\eqref{6}\) 函数 \(f\) 为：
\[f(x) = e^{cx^2} \tag{13} \label{13}\]
那么常数 \(c\)是多少呢？根据概率论我们知道：
\[\int f(x) \mathrm{d}x = \int e^{cx^2} \mathrm{d}x =1 \tag{14} \label{14}\]
上式就是大名鼎鼎的高斯积分（ \(\text{Gaussian Integral}\) ） 8 ，如何求解高斯积分可以参考这篇文章：数学之美：几何视角下的高斯积分(GaussianIntegral) 。
最终我们得出了 \(2\)维情况下的未归一化的概率密度函数：
\[f_2(x, y) = f_1(x)f_1(y) = e^{-x^2} e^{-y^2} = e^{-(x^2 +y^2)} \tag{15} \label{15}\]
上一节我们使用了不那么严谨的方法得到了正态分布的概率密度函数，下面我们使用另外一种方法求出正态分布的概率密度函数。
由落点分布分布在空间上具有旋转对称性，我们可知 \(x\) 轴和 \(y\)轴具有相同且连续的概率密度函数。
设落点 \(P\) 的概率密度函数为 \(\rho (x, y)\) ， \(x\) 轴方向上概率密度函数为 \(f(x)\) ，则 \(y\) 轴方向上的概率密度函数为 \(f(y)\) ，那么考虑如下图 12所示的一个充分小的黄色区域 \(\Box ABCD\) ：
飞镖落在黄色区域 \(\Box ABCD\)的概率为：
\[\rho (x, y) \mathrm{d}x \mathrm{d}y = f(x) \mathrm{d}x \cdot f(y)\mathrm{d}y \tag{16} \label{16}\]
将等式左边转换为极坐标形式，
\[\begin{cases}x = r \cos \theta \\y = r \sin \theta\end{cases}\]
在极坐标下的概率密度函数设为 \(g(r, \theta)\) , 则有：
\[\rho (x, y) = \rho (r \cos \theta , r \sin \theta) = g(r, \theta) \]
由条件 \(2\) , \(g(r, \theta )\) 具有旋转对称性，也就是和\(\theta\) 无关，所以
\[\frac{\mathrm{d} g(r, \theta)}{\mathrm{d} \theta} =0 \]
对公式 \(\eqref{16}\) 两边对 \(\theta\) 求导，可得：
\[\frac{\mathrm{d} f(x)}{\mathrm{d} \theta } f(y) + f(x) \frac{\mathrm{d}f(y)}{\mathrm{d} \theta } = 0 \]
利用链式法则，有：
\[-r \sin \theta f'(x)f(y) + f(x) f'(y) r \cos \theta =0 \]
上式移项可得：
\[\frac{f'(x)}{f(x) x} = \frac{f'(y)}{f(y) y} =C \]
我们则有：
\[\frac{f'(x)}{f(x)} = C x \]
对上式进行积分，可得：
\[\int \frac{f'(x)}{f(x)} \mathrm{d}x = \int C x\mathrm{d}x \]
求解上式可得：
\[\ln f(x) = \frac{C}{2}x^2 + C' \]
则求得函数 \(f(x)\) 为：
\[f(x) = A e^{\frac{1}{2}Cx^2} \tag{17}\label{17}\]
同理 \(f(y)\) 为：
\[f(y) = A e^{\frac{1}{2}Cy^2} \tag{18}\label{18}\]
由概率论我们知道 \(C < 0\) ，同时\(\int_{0}^{\infty } f(x) \mathrm{d}x =\dfrac{1}{2}\) ，则：
\[\int_{0}^{\infty } e^{\frac{C}{2} x^2} \mathrm{d}x =\frac{1}{2A} \]
对高斯核函数进行升维，对于 \(2\)维正态分布，从笛卡尔坐标系转换为极坐标系则有：
\[\begin{aligned}\int_{0}^{\infty} \int_{0}^{\infty} e^{\frac{C}{2}(x^2 + y^2)}\mathrm{d}x \mathrm{d}y & = \int_{0}^{\infty} \int_{0}^{\infty}e^{\frac{C}{2}r^2} r \mathrm{d}r \mathrm{d}\theta \\& = \int_{0}^{\frac{\pi }{2}} \int_{0}^{\infty} e^{\frac{C}{2}r^2} r\mathrm{d}r \mathrm{d} \theta \\& = \frac{1}{4A^2}\end{aligned}\]
令 \(u = r^2\) ，则有：
\[\int_{0}^{\infty} e^{\frac{C}{2}r^2} r \mathrm{d}r = \int_{0}^{\infty}e^{\frac{C}{2}u} \frac{\mathrm{d}u}{2} = \frac{1}{2} \left [ \frac{2}{C}e^{\frac{C}{2} u} \right ]_{0}^{\infty } =-\frac{1}{C} \]
可得：
\[-\frac{1}{C} \int_{0}^{\frac{\pi }{2}} \mathrm{d}\theta = -\frac{1}{C}\frac{\pi}{2} = \frac{1}{4A^2} \]
所以可求得：
\[A = \sqrt{\frac {-C}{2 \pi}} \]
至此，我们只剩下一个未知参数 \(C\) 就得到所求公式 \(\eqref{17}\) 。
考虑方差 \(\sigma^2\) 定义，对于期望\(\mu = 0\) ，则有：
\[\sigma ^2 = \int_{-\infty }^{\infty }x^2 f(x) \mathrm{d}x = 2\sqrt{\frac {-C}{2 \pi}} \int_{0}^{\infty }x^2 e^{\frac{C}{2} x^2}\mathrm{d}x \tag{19}\label{19}\]
令 \(u = x\) ， \(v = \dfrac{1}{C} e^{\dfrac{C}{2} x^2}\)，则有：
\[\mathrm{d}v = x e^{\frac{C}{2} x^2}\]
根据分部积分公式：
\[\int u \frac{\mathrm{d}v}{\mathrm{d}x}\,\mathrm{d}x = uv - \int\frac{\mathrm{d}u}{\mathrm{d}x}v\,\mathrm{d}x\]
则对公式 \(\eqref{19}\)进行分部积分求解可得：
\[\begin{aligned}\int_{0}^{\infty } u \frac{\mathrm{d}v}{\mathrm{d}x}\,\mathrm{d}x &= \lim_{x \to \infty} \frac{x}{C} e^{\frac{C}{2} x^2} - \frac{0}{C}e^{\frac{C}{2} 0} - \frac{1}{C} \int_{0}^{\infty } e^{\frac{C}{2} x^2}\mathrm{d}x \\& = 0 - 0 - \frac{1}{C} \frac{\sqrt {2 \pi }}{2 \sqrt {-C}} \\& = - \frac{1}{C} \frac{\sqrt {2 \pi }}{2 \sqrt {-C}}\end{aligned} \]
所以我们求得 \(\sigma^2\) 为：
\[\sigma^2 = 2 \sqrt{\frac {-C}{2 \pi}} \frac{-1}{C} \frac{\sqrt {2 \pi}}{2 \sqrt {-C}} = -\frac{1}{C} \]
即：
\[C = -\frac{1}{\sigma^2} \]
将求得 \(A\) 和 \(C\) 代入公式 \(\eqref{17}\)，最终我们求得概率密度函数为：
\[f(x) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac {x^{2}}{2\sigma^{2}}}} \tag{20}\label{20}\]
公式 \(\eqref{20}\) 是期望 \(\mu = 0\) 的特殊情况，当期望 \(\mu \ne 0\) 时，更一般的公式为：
\[f(x) = {\frac {1}{\sigma {\sqrt {2\pi }}}}\;e^{-{\frac {\left(x-\mu\right)^{2}}{2\sigma^{2}}}} \tag{21}\label{21}\]
特别地，当 \(\mu = 0\) 且 \(\sigma = 1\)，这个分布被称为标准正态分布：
\[f(x)={\frac {1}{\sqrt {2\pi }}}\, e^{-{\frac{x^{2}}{2}}} \tag{22}\label{22}\]
通过 \(\text{Herschel}\)给出的优雅直观方法，仅仅依靠那 \(2\)个假设条件，我们居然最终求出了正态分布的公式。有没有感觉到数学的美感？
最初看到 3Blue1Brown的这个视频，感觉非常美，正态分布那么复杂的公式居然有这么优雅直观的方式自然而然的出来了！
分析正态分布公式，公式中的 \(\pi\)意味着空间上的对称性，即点分布距离中心是对称的。而\(e\)的出现意味着取了时间上的极限，而这和中心极限定理（\(\text{Central limit theorem}\) ） 9有关，我会在下一篇文章详细解释，敬请期待！
2024-04-21 08:08:32
By Long Luo
世界那么大，我想去看看！随着科技的发展，我们早已做到无需离开家就能领略世界各地的风景和文化。通过高清视频、高清直播、社交媒体、VR技术、各种图片或者视频分享平台，我们不仅可以体验世界名胜古迹、自然美景和各地的风土人情，还能与当地居民互动，了解他们的日常生活和传统文化。互联网丰富了人们的生活，缩小了地域的界限，真正实现了让世界触手可及，足不出户便可周游世界的梦想。但正所谓“百闻不如一见”，“读万卷书，不如行万里路！”，尽管互联网让我们可以虚拟游览世界，但亲自出行的体验无可替代。自由行不仅提供了前所未有的自由和灵活性，还能让我们亲身感受到大自然的美妙、城市的活力。这种身临其境的体验，远非屏幕前的感受可比。
当你打开地图软件时，你会看到如图 1 所示的道路标志，
当你在道路出行时，你也会看到看到如下图 2所示的路牌，但你可能并未真正留意过这些标志。因为现在我们只需要有一部联网的智能手机，在地图类软件里，设定出发地和目的地，自然有导航会指引我们到达目的地。
这些编号肯定不是随机的，那么这些道路编号到底有什么用呢？出于好奇心你可能会去寻找答案，你很容易轻松找到中国国家高速的编号密码 1 和 高速公路是怎样命名和编号 2这样的文章。但这些文章只是告诉我们是什么，并没有告诉我们为什么。
国内现行的高速公路命名是由交通部从 2005 年启动的 3，之前道路都以起始地和终点地命名。你可能会想，之前那种命名方式不是更合理吗？用了数字不是更加不清晰易懂吗？如果图1 不是路牌上写了济广高速，谁知道 G35 4 是哪里到哪里呢？
如果你更进一步的话，如果你去了解世界其他国家的高速公路命名的话，你会惊讶的地发现为什么居然全世界各主要大国都选择了类似的编号系统，这背后的原因是什么呢？
要回答这个问题，我们需要把时钟拨回几十年前，回到高速公路诞生的时期，那个没有GPS，没有手机，只有纸质地图的时代，我们才能知道这种编号系统的重要意义和实用性，以及背后的数学密码。
我们已经了解了国内高速公路编号 5，让我们看看其他国家的高速公路系统编号是什么样的。这些国家需要国土面积足够大，高速公路系统足够发达，国土疆域长宽比例没有太夸张，人口分布比较均匀。
首先我们把目光投向德国 6，因为世界上第一条被广泛认可的现代意义上的高速公路就是德国的“Autobahn”7 ，于 1932年建成，连接科隆与波恩，完全是为高速行驶的机动车辆设计和建造的。
在纳粹上台之后，高速公路建设被视为国家的重要工程，是当时德国最具标志性的基础设施项目之一，以展示德国的工程实力和国家威严8。高速公路建提高了交通运输效率和便利性，同时也为了支持军事行动和纳粹的宣传活动。这些公路被设计成宽敞平坦、直线且无障碍的道路，可以实现高速行驶。
这些高速公路被称为“Autobahn”，德语的意思就是“汽车道”，命名方式则采用了简单的字母和数字组合，例如\(A1, A2, \cdots , A9\) 等 9。这种命名方式简洁明了，便于驾驶员识别和导航。
德国现在高速公路 (Bundesautobahn，BAB) 的编号系统是在 1974年开始使用的，所有的高速公路都以 \(A\)开头后接一个空白与数字编号，如 \(A8\)。穿越德国全国的东西向主要高速公路以偶数编号，南北向的道路则以单数编号。用来连接区域性重点城市的较短高速公路则以两位数字来编号，具体如下图3 所示：
当然德国也接入了欧盟标准的“E”编号系统与其他欧洲国家的公路系统协调，如下图 4 所示：
法国的高速公路系统 10被称为“Autoroute”，通常以字母“A”开头，后跟一个数字。编号体系如下：
由于巴黎长期是法国的中心 11，所以法国的高速系统沿巴黎放射性分布，如下图 5所示。当然法国是欧洲国家，所以也使用了欧盟标准的“E”编号系统与其他欧洲国家的公路系统协调。
俄罗斯的高速公路系统 12 被称为 Avtomagistral，俄语里就是高速公路的意思，编号使用字母“M”、“A”或“P”开头，后跟数字：
俄罗斯和法国类似，精华都在莫斯科圣彼得堡的东欧平原区，西伯利亚和远东地区人口稀少，高速分布如下图6 所示：
欧洲的高速公路系统 13使用“E”字母（Europe）作为前缀，后接一个数字，代表着欧洲主要道路网。
欧洲路线编号的决定是由 1975 年的欧洲经济委员会作出的，并在 1992年进行了修正。作为干线的A级道路编号是两位数，而作为地方线的B 级的道路编号则是三位数，主要规则如下所示：
当然还有其他一些规则，欧洲高速公路分布如下图 7 所示：
由于加拿大 14 和 巴西 15人口分布不均匀，没有什么借鉴价值。当然有同学看到这会问，欧洲高速公路系统和国内编号规则几乎一致，那么国内是借鉴了欧洲的公路编号系统吗？
其实世界各大国的高速公路编号体系在设计上几乎一致，这种统一编码方案可以追溯到美国的州际高速公路16系统，这个我们会在下面的章节里进行详细讲解。
我们可以先问下自己几个问题：
假设我们在下图 8 所示的纽约曼哈顿的 \(A\)点，我们要去古根海姆博物馆。此时我们位于 东区第 87 大街和 第 3大道的十字路口，我们要去的目标在 东区第 88 大街和 第 5大道交汇处十字路口。很显然，我们有很多条路线可以到达目的地。
即使你是个路痴，当你沿着东区第 87大街走的时候，如果你走着走着发现走到了第 4大道，你就知道至少在东西方向上你的方向对了，而当你发现走到了第 2大道，你就知道你方向走反了。同理，对于南北方向上也是如此。这样，即使你绕了很多路，但你只要知道目的地所在街区的编号，你就调整方向，多走点路多花点时间也一定能走到目的地。
当然我在写这篇文章时，想找个浅显易懂的例子，结果发现曼哈顿街区完美的符合我想找的例子。因为道路规划太好了，呈棋盘状。实际中大部分城市街区命名不一定使用数字，也不一定如此横平竖直。但是曼哈顿街区因为都使用数字编号，我们想找具体地点太容易了，这还没用上二分查找呢！
由于现实中由于不同地区地理条件不一样，城市的分布往往在地图上是随机的，比如下图9 就是我随手绘制的一个道路分布图。
如何对这些道路进行命名？你可能会说，那就以出发地和目的地命名好了，比如\(AB\) , \(BC\)等。当然这样命名也可以，但是我们可以更进一步，让命名和曼哈顿街区命名一样，让名字就自带导航功能。
那比如有 100 条道路，如何给每条道路编上正确的号码呢？从 1开始递增当然也可以，但就没有任何信息了。
让我们回顾现实中的高速公路，我们知道一条高速公路有哪些属性呢？
还记得地球仪上的经线和纬线吗？墨卡托投影制作的地图帮助无数航海家顺利抵达目的地，航海家只要知道测量出当前所在地经纬度，沿着连接目的地的方向就不会迷失方向。
想一想，如何数字也有很多属性，比如奇偶性、大小、位数等，参考之前的编号系统，我们可以用数字属性和高速公路属性一一对应17 ，比如：
不得不说，这种设计比曼哈顿的纯数字编号要好很多，绝对是天才般的设计 18！下面我们就来介绍这种编号系统的诞生地：美国州际高速公路系统。
也许是艾森豪威尔总统 19二战时担任盟军统帅，战后又担任美国驻德国占领军总司令的缘故，我想当时德国领先世界的高速公路系统肯定给他留下了深刻印象。在他1953 - 1961年担任总统期间，提出了州际高速公路系统的构想，并在1956年签署了《国家公路与交通法案》，旨在建立一套横贯美国的高速公路网络，以提升国家交通运输效率和安全性。这个宏伟的工程由联邦政府主导，并联合各州政府共同筹资和建设。经过数十年的不懈努力，美国州际高速公路系统于1992年正式完成，涵盖了数万英里的公路，连接了各个州份和城市，成为美国国家交通运输体系的重要组成部分，也促进了经济的发展和人口的流动。
州际公路系统每条线路和编号，可以参考这个在线网站：Interstate Highway Map
下面我们来介绍下州际高速公路的编号规则 20：
其实也很好记， \(\text{Even}\) 对应\(\text{East}\) ，而且 \(0\)不偏不倚，表示水平方向，很容易记吧！
其实也很好记，偶数表示东西向，那奇数就是南北向了，这是主要干道，肯定不能依次递增，所以每次递增\(5\) 。注意西海岸的 \(5\) 要当成 \(05\) 。
如何理解数字增加方向呢？
由于美国本土大致呈矩形，如果我们在地图左下角作为坐标原点，建立一个笛卡尔坐标系，如下图12 所示，想必你忘不了：
由于中国的国土疆域并不像美国接近一个矩形，所以坐标点实际是放在靠近北京的右上角，如下图13 所示：
如果你掌握了这个规则，那么下面这个路标，你能猜出大概在美国哪里吗？可以先猜测然后自寻答案验证下你的猜测。
相比 \(2\)位数字高速公路是跨州线路(注：不完全如此)， \(3\)位数的州际公路是服务于各个都会区的较短路线，它们连接到更长的两位数路线，并且充当环城公路、支线或连接线路。
\(3\) 位数编号规则如下：
国内公路编号系统则大致相同，可以参考 主干道 21和 辅道 22 ，具体可参考 中国高速公路网 。
上面就是美国州际高速公路的大致规则，当然也会有一些例外。
这套编号系统不仅方便导航，有效帮助驾驶员辨别方向，还提供了关于道路走向和地理位置的信息，提高了道路网络的整体效率和安全性，它是如此的成功，以致于后来陆续被其他国家所采用。
生活处处皆学问，正如王建硕老师写的系列文章：无用的冷知识。如果我们对日常事务多一点好奇心，我们可以发现后面更多的原理！
2024-04-16 11:13:48
By Long Luo
阿里巴巴达摩院 从2018年开始每年都会举办一届全球数学竞赛，之前一方面自己数学水平比较弱，另外一方面也没有报名，但一直很仰慕那些数学大神的风采。今年是第一次报名参加2024阿里巴巴全球数学竞赛，上周末参加了预选赛，但遗憾的是，全部 \(7\) 道题中只有第 \(1, 2, 6\) 题会做，这里分享下我的解答：
几位同学假期组成一个小组去某市旅游. 该市有 \(6\) 座塔，它们的位置分别为 \(A, B, C, D, E, F\)。同学们自由行动一段时间后，每位同学都发现，自己在所在的位置只能看到位于\(A, B, C, D\) 处的四座塔，而看不到位于\(E\) 和 \(F\) 的塔。已知：
（5 分） 请问 这个旅游小组最多可能有多少名同学？
\(A. 3\)
\(B. 4\) \(C.6\) \(D. 12\)
这道题选 \(C\) ，最多只能有 \(6\) 名同学。
这道题的解题思路是从假设只有 \(1\) 座塔开始，一直到 \(6\) 座塔，找到思路。
假设有 \(1\) 座塔 \(A\) ，那么很显然有无数多同学可以看到塔\(A\) ，也可以有无数多同学看不到塔\(A\) ；
假设有 \(2\) 座塔 \(A, B\) ，那么只有以 \(A\) 为起点的射线 \(AB\) 且位于 \(B\) 之后的同学无法看到塔 \(A\) ；
假设有 \(3\) 座塔 \(A, B, C\)，同理可知存在无数位同学至少可以看见 \(2\) 座塔；
假设有 \(4\) 座塔 \(A, B, C, D\)，同理可知存在无数位同学至少可以看见 \(2\) 座塔；
假设有 \(6\) 座塔 \(A, B, C, D, E, F\) ，如果每位同学都无法看见\(E, F\) 塔，如下图1 所示：
所以至多有 \(6\) 位同学位于 \(M, N, O, P, R, Q\) 处，无法看到塔 \(E, F\) 。
小明玩战机游戏。初始积分为 \(2\)。在游戏进行中，积分会随着时间线性地连续减少 (速率为每单位时间段扣除\(1\) )。游戏开始后，每隔一个随机时间段(时长为互相独立的参数为 \(1\)的指数分布)，就会有一架敌机出现在屏幕上。当敌机出现时，小明立即进行操作，可以瞬间击落对方，或者瞬间被对方击落。如被敌机击落，则游戏结束。如小明击落敌机，则会获得\(1.5\)个积分，并且可以选择在击落该次敌机后立即退出游戏，或者继续游戏。如选择继续游戏，则须等待到下一架敌机出现，中途不能主动退出。游戏的难度不断递增：出现的第\(n\) 架敌机，小明击落对方的概率为\((0.85)^n\) ，被击落的概率为 \(1 - (0.85)^n\)，且与之前的事件独立。在任何时刻，如果积分降到 \(0\) ，则游戏自动结束。
小问 1 (5分)如果游戏中，小明被击落后，其之前的积分保持。那么为了游戏结束时的累积积分的数学期望最大化，小明应该在其击落第几架敌机后主动结束游戏？
\(A. 1\) \(B. 2\) \(C.3\) \(D. 4\)
这道题考察的就是泊松过程，数学好的同学推出其表达式，然后计算可得。
根据小概率事件原理，我们可以把打飞机事件视为泊松过程1，那么事件分布就是泊松分布2 。其 Python代码可以直接调用 API numpy.random.exponential。虽然是 指数分布，但在 Java 中 需要使用 -Math.log(1 - random.nextDouble())而不是 Math.exp(double a) 。模拟代码如下所示：
1 |
public class WarPlaneGame { |
输出如下：
1 |
Shoot down 1 planes, Expect Score: 2.2344824425425442 |
可以看出击落第 \(2\)架敌机后主动结束游戏，期望积分最大，所以答案选 \(B\) 。
小问 2 (5分)如果游戏中，小明被击落后，其之前积累的的积分会清零。那么为了游戏结束时的期望积分最大化，小明也会选择一个最优的时间主动结束游戏。请问在游戏结束时（小明主动结束游戏、或积分减到\(0\)），下列哪一个选项最接近游戏结束时小明的期望积分？
\(A. 2\)
\(B. 4\) \(C.6\) \(D. 8\)
通过第一问，我们知道期望积分是随着次数逐渐递减的。
继续写代码模拟其过程，如下所示：
1 |
public class WarPlaneGame { |
输出如下：
1 |
Shoot down 1 planes, Expect Score: 2.0230245088842116 |
可以看出最大期望积分是 \(2.023\)左右，所以答案选 \(A\) 。
对于实数 \(T > 0\) ，称欧式平面\(\mathbb{R}^2\) 的子集 \(\Gamma\) 为 \(T\) -稠密的，如果对任意 \(v \in \mathbb{R}^{2}\) ，存在 \(w \in \Gamma\) 满足 \(\|v-w\| \leqslant T\) . 设 \(2\) 阶整方阵 \(A\in \mathrm{M}_{2}(\mathbb{Z})\) 满足 \(\operatorname{det}(A) \neq 0\) .
\[A^{n} \mathbb{Z}^{2}:=\left\{A^{n} v: v \in \mathbb{Z}^{2}\right\}\]
是 \(C|\operatorname{det}(A)|^{n /2}\) -稠密的.
注: 这里 \(\mathbb{R}^{2}\) 和 \(\mathbb{Z}^{2}\) 中的向量约定为列向量,\(\mathbb{R}^{2}\) 中的内积为标准内积,即 \(\langle v, w\rangle=v^{t} w\) .(提示: 在对(2)的证明中, 可使用如下 \(\text{Minkowski}\) 凸体定理的特殊情形：\(\mathbb{R}^{2}\)中以原点为中心且面积为 \(4\)的任意闭平行四边形中总包含 \(\mathbb{Z}^{2}\) 中的非零向量.)
先挖坑，等我看懂了大神的解答再来填坑！
设 \(d \geq 0\) 是整数， \(V\) 是 \(2d+1\) 维复线性空间，有一组基
\[\left\{v_1, v_2, \cdots, v_{2 d+1}\right\} \text {. }\]
对任一整数 \(j\left(0 \leq j \leq\dfrac{d}{2}\right)\) ，记 \(U_j\) 是
\[v_{2 j+1}, v_{2 j+3}, \cdots, v_{2 d-2 j+1}\]
生成的子空间. 定义线性变换 \(f: V\rightarrow V\) 为
\[f \left( v_i \right) = \frac{(i-1)(2d + 2 - i)}{2} v_{i-1} + \frac{1}{2}v_{i+1}, 1 \leq i \leq 2 d+1 .\]
这里我们约定 \(v_0=v_{2d+2}=0\).
证明题 (10分) 证明: \(f\) 的全部特征值为 \(-d,-d+1, \cdots, d\).
问答题 (5分) 记 \(W\) 是从属于特征值 \(-d+2 k(0 \leq k \leq d)\) 的 \(f\) 的特征子空间的和. 求 \(W \cap U_0\) 的维数.
问答题 (5分) 对任一整数 \(j\left( 1 \leq j \leq \dfrac{d}{2}\right)\) ，求 \(W \cap U_j\)的维数.
先挖坑，等我看懂了大神的解答再来填坑！
证明题 (20分) 对于 \(\mathbb{R}^3\) 中的任何中心对称的凸多面体\(V\) ，证明可以找到一个椭球面 \(E\) ，把凸多面体包在内部，且 \(E\) 的表面积不超过 \(V\) 的表面积的 \(3\) 倍。
先挖坑，等我看懂了大神的解答再来填坑！
假设有一枚硬币，投掷得到正面的概率为 \(\dfrac{1}{3}\) 。独立地投掷该硬币 \(n\) 次，记 \(X_n\) 为其中得到正面的次数。试求 \(X_n\) 为偶数的概率在 \(n\) 趋于正无穷时的极限。
当 \(n \to \infty\)，直觉告诉我们，偶数次正面出现的概率和奇数次正面出现的概率是一样的，而奇数和偶数是均匀分布的，答案应该是\(\dfrac{1}{2}\)。但这道题不是选择题也不是填空题，我们需要严谨证明这个结论！
由题意可知，设随机变量 \(X_n\)表示在 \(n\)次独立投掷中正面出现的次数，每次出现正面的概率为 \(p = \dfrac{1}{3}\) ，则 \(X_n\) 服从参数为 \(\operatorname{B}(n, p)\) 的二项分布3 ，那么 \(n\) 次独立投掷中正面出现 \(k\) 次的概率是：
\[\begin{equation}\operatorname{Pr}(X_n = k) = \binom{n}{k} p^k (1−p)^{n−k} \tag{6.1.1}\label{6.1.1}\end{equation}\]
要求 \(X_n\) 为偶数的概率，即：
\[\begin{aligned}\operatorname{Pr}(X_n \text { is even}) & = \operatorname{Pr}(X_n =0) + \operatorname{Pr}(X_n = 2) + \cdots + \operatorname{Pr}(X_n = 2k, k= \left \lfloor \frac{n}{2} \right \rfloor ) \\& = \sum_{k=0}^{\left \lfloor \frac{n}{2} \right \rfloor}\binom{n}{2k} p^{2k} (1 − p)^{n − 2k}\end{aligned}\]
带入 \(p = \dfrac{1}{3}\)，可得：
\[\begin{equation}\operatorname{Pr}(X_n \text { is even}) = \sum_{k=0}^{\left \lfloor\frac{n}{2} \right \rfloor} \binom{n}{2k} (\frac{1}{3})^{2k}(\frac{2}{3})^{n−2k} \tag{6.1.2} \label{6.1.2}\end{equation}\]
由 二项式定理4 可知：
\[\begin{equation}(x + y)^n = \sum_{k=0}^{n} \binom{n}{k} x^k y^{n−k} \tag{6.1.3}\label{6.1.3}\end{equation}\]
那么易得共轭表达式：
\[\begin{align}(x + y)^n + (x − y)^n & = 2 \sum_{k=0}^{\left \lfloor \frac{n}{2}\right \rfloor} x^{2k}y^{n - 2k} \tag{6.1.4} \label{6.1.4} \\(x + y)^n - (x − y)^n & = 2 \sum_{k=0}^{\left \lfloor \frac{n}{2}\right \rfloor} x^{2k + 1}y^{n - 2k - 1} \tag{6.1.5} \label{6.1.5}\end{align}\]
可得：
\[\begin{aligned}\operatorname{Pr}(X_n \text { is even}) & = \sum_{k=0}^{\left\lfloor \frac{n}{2} \right \rfloor} \binom{n}{2k} (\frac{1}{3})^{2k}(\frac{2}{3})^{n−2k} \\& = \frac{1}{2} \left [ \left (\frac{1}{3} + \frac{2}{3} \right )^n+ \left (\frac{1}{3} - \frac{2}{3} \right )^n \right ] \\& = \frac{1}{2} \left [1 + \frac{1}{3^n} \right ]\end{aligned}\]
故答案为：
\[\lim_{n \to \infty} \operatorname{Pr}(X_n \text { is even}) = \lim_{n\to \infty} \frac{1}{2} \left (1 + \frac{1}{3^n} \right ) = \frac{1}{2}\]
这道题也可以用 马尔可夫链来做，构建递推关系式，感兴趣的同学可以试试！
某人在过年期间参加了集五福活动，在这项活动中此人每扫描一次福字，可以随机地得到五张福卡中的一张。假设其每次扫福得到五福之一的概率固定，分别为\(p_i \in (0, 1) , i = 1, 2, \cdots ,5\) 且 \(\sum_{i = 1}^{5} p_i =1\) ，并假设其每次扫描得到的结果相互独立。在进行了 \(n\) 次扫福之后，记 \(X^{i}_n, i =1, 2, \cdots, 5\)为其得到每种福卡的张数。那么求极限 \(\lim _{n\to \infty} \operatorname{P} \left ( X^{(i)}_{2n}, i = 1, 2, \cdots, 5\text { 全部为偶数} \right )\) 。
直觉告诉我们，当 \(n \to \infty\)时，五种福卡每种都是偶数的事件是相互独立的。通过第一问，我们已经知道答案是\(\dfrac{1}{2}\)，那么五种福卡每种福卡的张数都是偶数的概率就是 \(\dfrac{1}{2^5} = \dfrac{1}{32}\) ，而 \(2n\) 次扫福卡的概率就是 \(\dfrac{1}{16}\)。这个猜测对不对呢？下面我们就来证明下。
由多项式定理5 ：
\[\left (x_{1} + x_{2} + \cdots + x_{m} \right)^{n} = \sum_{\begin{array}{c} \alpha _{1} + \alpha _{2} + \cdots + \alpha _{m} = n\\ \alpha _{1},\alpha _{2},\cdots ,\alpha _{m} \geq 0 \end {array}}{\frac {n!}{\alpha _{1}! \dots \alpha _{m}!}}x_{1}^{\alpha _{1}} \dotsx_{m}^{\alpha _{m}} \tag{6.2.1} \label{6.2.1}\]
设 \(k_i, i = 1,2, \cdots ,5\)表示是 \(n\) 次独立扫描福卡中得到第\(i\) 种福卡的张数，则其概率为：
\[\begin{align}\operatorname{Pr} \left (X_{n}^{(i)} = k_i, i=1,2, \cdots ,5 \right)&= \binom {n}{k_1, k_2, \cdots, k_5} p_1^{k_1} p_2^{k_2} \cdotsp_5^{k_5} \nonumber \\& = \frac {n!}{k_1!k_2! \cdots k_5!} p_1^{k_1} p_2^{k_2} \cdotsp_5^{k_5} \tag{6.2.2} \label{6.2.2}\end{align}\]
观察上式可知，所求概率为多项式 \(\left (p_1+ p_2 + p_3 + p_4 + p_5 \right )^{n}\) 的 \(p_1^{k_1} p_2^{k_2} \cdots p_5^{k_5}\)项，\(\binom {n}{k_1, k_2, \cdots ,k_5}\) 为其系数。
和问题 \(1\)的共轭表达式类似，我们给不同福卡添加符号位，考虑如下求和表达式：
\[S_{2n} = \frac{1}{2^{5}} \sum_{\beta_{i} = \pm 1} \left (\beta_{1} p_{1}+ \beta_{2} p_{2} +\cdots + \beta_{5} p_{5}\right)^{2n} \tag{6.2.3} \label{6.2.3}\]
对上式进行多项式展开，可得：
\[\begin{aligned}S_{2n} & = \frac{1}{2^{5}} \sum_{\substack {\beta_{i} = \pm 1 \\x_{1} + x_{2} + \cdots + x_{5} = 2n}} \frac{(2n)!}{x_{1}!x_{2}! \cdotsx_{5}!} \beta{1}^{x_{1}} \beta_{2}^{x_{2}} \cdots \beta_{5}^{x_{5}}p_{1}^{x_{1}} p_{2}^{x_{2}} \cdots p_{5}^{x_{5}} \\& = \frac{1}{2^{5}} \sum_{x_{1} + x_{2} + \cdots + x_{5} = 2n}\frac{(2n)!}{x_{1}!x_{2}! \cdots x_{5}!} p_{1}^{x_{1}} p_{2}^{x_{2}}\cdots p_{5}^{x_{5}} \sum_{\beta_{i} = \pm 1} \beta_{1}^{x_{1}}\beta_{2}^{x_{2}} \cdots \beta_{5}^{x_{5}}\end{aligned}\]
考虑 \(\sum _{\substack {\beta_{i} = \pm 1\\x_{1} + x_{2} + \cdots + x_{5} = 2n}} \beta_{1}^{x_{1}}\beta_{2}^{x_{2}} \cdots \beta_{5}^{x_{5}}\) ，如果存在 $k $ 使得\(x_{k}\)为奇数的话，那么：
\[\sum_{\beta_{i} = \pm 1} \beta_{1}^{x_{1}} \beta_{2}^{x_{2}} \cdots\beta_{5}^{x_{5}} = \sum_{\substack {\beta_{i}= \pm 1 \\ i \neq k}}\left[\left (1^{x_{k}} + (-1)^{x_{k}} \right) \prod_ {i \neq k}\beta_{i}^{x_{i}} \right ]=0\]
由于奇数项最终都会消去，只有偶数项 \(x_{i}\) 才会留下来，故有：
\[\sum_{\beta_{i} = \pm 1} \beta_{1}^{x_{1}} \beta_{2}^{x_{2}} \cdots\beta_{k}^{x_{k}} = 2^{k}\]
那么求和表达式为：
\[\begin{aligned}S_{2n} & = \sum_{\substack{x_{1} + x_{2} + \cdots + x_{k} = 2 n \\x_{i} \text { is even }}} \frac{(2 n)!}{x_{1}!x_{2}! \cdots x_{k}!}p_{1}^{x_{1}} p_{2}^{x_{2}} \cdots p_{k}^{x_{k}} \\& = \operatorname{Pr} \left \{ X_{2 n}^{(i)} \text { is all even }\right\}\end{aligned}\]
因此，所求问题转化为在 \(X_{2n}^{(i)}\) 均为偶数情况下，当 $n $ 时，其极限为：
\[\lim _{n \to \infty} \operatorname{Pr}\left \{X_{2 n}^{(i)} \text { isall even } \right\} = \lim _{n \to \infty} S_{2 n}\]
因为 \(\left | \beta_{1} p_{1} + \beta_{2}p_{2} + \cdots + \beta_{5} p_{5} \right | \leq 1\) , 所以当 \(n \to \infty\) 时，只有 \(\beta_i\) 全为 \(1\) 或者 \(-1\) 情况下，
\[\left | \sum_{i=1}^{5} \beta_{i} p_{i} \right | = 1\]
因此，我们可以得到答案：
\[\lim _{n \to \infty} \operatorname{Pr} \left \{X_{2n}^{(i)} \text { isall even } \right \} = \frac{1}{2^{5}} \left [ (+1)^{2n} + (-1)^{2n}\right ] = \frac{1}{16}\]
有这么一个音乐盒，它上面有一个圆形的轨道，轨道上的一点处还有一棵开花的树。当音乐盒处于开启模式时，音乐盒会发出音乐，轨道会按照顺时针匀速转动。
你可以在轨道上放置象征恋人的两颗棋子，我们不妨称它们为小红和小绿。当小红和小绿没有到达树下时，它们就会在轨道上各自移动。当某一颗棋子到达树下时，它就会在树下原地等待一段时间。此段时间内，如果另外一颗棋子也达到了树下，那么两颗棋子就会相遇，之后在它们将随即一起顺着轨道移动，不再分开；否则，等待时间结束，两颗棋子将各自顺着轨道继续移动。
考虑这个音乐盒的数学模型。我们把这个圆形轨道参数化成一个周长为 \(1\)的圆环，我们认为棋子和树都可以用圆环上点表示。具体来说，我们用 \(X(t) \in [0, 1]\) 和 \(Y(t) \in [0, 1]\) 分别表示 \(t\)时刻小红和小绿的在轨道上的位置坐标，而树的坐标是 \(\phi = 1\) ，或者，等价地， \(\phi = 0\) 。
当他们都没有抵达树下时 (见左图) ，他们的位置变化规律满足
\[\frac{\mathrm{d}}{\mathrm{d} t} X(t)=1, \quad\frac{\mathrm{d}}{\mathrm{d} t} Y(t) = 1\]
假设在 \(t_0\)时刻，小绿到达了树下（见中图），即 \(Y \left(t_0 \right) = 1\) ，它就会至多等待
\[\tau = K \left (X \left( t_0 \right ) \right)\]
的时间，换句话说，最长等待时间依赖于小红的当时的位置。
在等待期间，小绿不动，小红继续移动。如果等待期间的某时刻 \(t^\ast \in \left(t_0, t_0+\tau\right]\)，小红也达到了树下，即 \(X\left (t^\ast\right) = 1\)，那么两棋子相遇。如果等待时间结束时（见右图），小红仍没有到达树下，那么它们俩继续移动，此时他们的位置分别是
\[X \left(t_0 + \tau \right) = X \left(t_0 \right) + \tau, \quad Y\left(t_0 + \tau \right) = 0 .\]
注意，虽然小绿的坐标被重置了，但是它在圆环上的位置并没有变。
如果在某时刻小红到达树下，它也会按照相同的规则等待，最长等待时间取决于此时小绿的位置。显然，小红小绿的命运取决于最长等待时间函数\(K(\phi)\) 的形式。
\[f^{\prime} > 0, \quad f^{\prime \prime} < 0, \quad f(0)=0, \quadf(1) = 1 .\]
并设 \(\varepsilon\)是一个充分小的正的常数。我们定义等待时间函数
\[K(\phi ) = f^{-1}(f(\phi ) + \epsilon ) - \phi .\]
证明除了唯一的例外（特定的初始距离）之外，无论小红和小绿的初始距离如何，他们最终会相遇的。
\[f(\phi ) = \frac {1}{b} \ln \left (1 + \left (e^b - 1 \right ) \phi\right )\]
这里 \(b>0\) 是一个常数。当 \(b \ll 1, \varepsilon \ll 1\)时，请估算出相遇之前小红小绿走过的圈数的数量级。
这道题考试的时候没做出来，最近几天看了知乎上关于这次考试的讨论 如何评价2024阿里巴巴数学竞赛预选赛试题?，看了大神们的解答，发现这道题不难，不要以为它是压轴题就觉得很难。这道题的关键在于找到小红和小绿的距离递推关系式，然后对这个关系式进行分析。下面的解法参考了知乎Fiddie 的解答[^6] ，喵喵 的解答[^7]，我综合他们的解题思路自己推了一遍。
由题设条件，我们知道当小红和小绿都没有在树下时，都随着圆形轨道顺时针匀速转动，也因此小红和小绿只可能在树下相遇。
不妨设最初条件为小红在小绿的前方，两者的距离为 \(d_0, \ d_0 \in (0, 1)\) ， \(d_0 \le 0\) 或 \(d_0 \ge 1\) 两者相遇无需讨论。
当小红来到树下，此时小绿距离树下还有一段距离，小红将在树下等待一段时间。假设在小红等待时间结束之前，小绿都没能赶到树下，那么在等待时间结束时，记为\(t_{k}\) 时刻。在 \(t_{k}\) 时小绿距离树下还有 \(d_k\) 的距离，即 \(X(t_k) = 0, \ Y(t_k) = 1 - d_k\) 。
小红继续出发，在 \(d_k\) 时之后，即\(t_k + d_k\)时小绿将到达树下。此时小红已经出发了 \(d_k\) 的距离，即 \(X(t_k) = d_k, \ Y(t_k) = 1\) 。
小绿将在树下等待小红 \(\tau = K \left (X\left (d_{k} \right ) \right )\) 的时间，即：
\[\begin{equation}\tau = K \left (X \left (d_{k} \right ) \right ) = f^{-1}(f(d_k) +\varepsilon) - d_k \tag{7.1.1} \label{7.1.1}\end{equation}\]
分析 \(\eqref{7.1.1}\) 可知：
如果 \(f^{-1}(f(d_k) + \varepsilon) \ge1\)，则小红将在等待时间结束之前到达树下，小红和小绿相遇，结束分析。
如果 \(f^{-1}(f(d_k) + \varepsilon)< 1\)，那么在小绿等待时间结束之前，小红没能赶到树下。
在小绿等待时间结束那一刻，我们记为 \(t_{k+1}\) 时刻，此时 \(X(t_{k+1}) = f^{-1}(f(d_k) + \varepsilon), \Y(t_{k+1}) = 0\) ，两者距离为：
\[\begin{equation}d_{k+1} = 1 - f^{-1}(f(d_k) + \varepsilon) \tag{7.1.2} \label{7.1.2}\end{equation}\]
至此我们找到了小红与小绿之间的距离递推关系式。
设两者之间距离数列 ${ d_n } $表示一个人刚要从树下出发，另外一个距离树下的距离，那么问题转化为：对于任意初值\(d_0 \in (0, 1)\)，除了某个特定的初始距离值之外，都存在某个 \(k\in \mathrm{Z^+}\) ，使得数列 \(d_k \le0\) 或者 \(d_k \ge 1\) 。
这里我们需要证明 \(2\) 种情况：
考虑函数 \(g(x)\)表示两者之间距离，函数 \(h(x)\)表示前后时刻（ \(t_{k}\) 和 \(t_{k+1}\)）两者之间的距离变化，即 \(\Delta d = d_{k+1} - d_k = 1 - f^{-1}(f(d_k) +\varepsilon) - d_k\) 。
则有：
\[\begin{equation}g(x) = 1 - f^{-1}(f(x) + \varepsilon) , \ x \in (0, 1) \tag{7.1.3}\label{7.1.3}\end{equation}\]
\[\begin{equation}h(x) = g(x) - x = 1 - f^{-1}(f(x) + \varepsilon) - x \tag{7.1.4}\label{7.1.4}\end{equation}\]
对 \(g(x) , \ h(x)\) 求导可得：
\[\begin{equation}g^{\prime}(x) = \frac{\mathrm{d} g(x)}{\mathrm{d} x} = -\frac{f^{\prime}(x)} {f^{\prime}\left(f^{-1}(f(x) +\varepsilon)\right)} \tag{7.1.5} \label{7.1.5}\end{equation}\]
\[\begin{equation}h^{\prime}(x) = \frac{\mathrm{d} h(x)}{\mathrm{d} x} = -1 -\frac{f^{\prime}(x)} {f^{\prime}\left(f^{-1}(f(x) +\varepsilon)\right)} \tag{7.1.6} \label{7.1.6}\end{equation}\]
因为 \(f^{\prime}>0 , f^{\prime \prime}< 0\) ，所以 \((f^{-1})^{\prime}>0\) , \(\left(f^{-1}\right)^{\prime \prime}>0\)，
又 \(h(0) = 1 - f^{-1}(\varepsilon) >0\) ， \(h(1) = - f^{-1}(f(1) +\varepsilon) < 0\) ，所以 \(h(x)\) 在
针对第 \(1\)种情况，需要证明：存在性和唯一性。
设距离数列 ${ d_n } $ 存在某个初值 \(d_0 =d^\ast\) 满足公式 \(\eqref{7.1.2}\) 使得 \(d_k = d^\ast , k = 0, 1, \cdots , n\)，所以：
\[\begin{equation}d_{k+1} = d_k \Leftrightarrow d^\ast = 1 - f^{-1} \left (f(d^\ast) +\varepsilon \right ) \end{equation}\]
进一步化简可得：
\[\begin{equation}f(1 - d^\ast) = f(d^\ast) + \varepsilon \end{equation}\]
设 \(x_1 = d^\ast, \ x_2 = 1 -d^\ast\) ，那么 \(x_1, \ x_2\)关于 \(\dfrac{1}{2}\) 对称，那么存在\(x_1 = \dfrac{1}{2} - \varepsilon^\ast ，x_2= \frac{1}{2} + \varepsilon^\ast, \ \varepsilon^\ast > 0\)。
根据题设条件 \(f^{\prime} > 0\)， \(f(0) = 0, \ f(1) = 1\) ， \(f\) 在 \([0,1]\) 是单调递增的光滑函数，那么 \(f(x_2) > f(x_1) > 0\)，故存在性得证。
下面来证明唯一性，由公式 \(\eqref{7.1.5}\) ，可得： \(f(1 - d^\ast) - f(d^\ast) = \varepsilon\)。
令函数 \(g(x) = f(1 - x) - f(x), \ x \in(0,1)\) ，对 \(g(x)\) 求导：
\[\begin{equation}g^{\prime}(x) = - f^{\prime}(1 - x) - f^{\prime}(x) < 0 , \ x \in(0,1) \end{equation}\]
\(g(x)\) 单调递减， \(g(0) = 1, \ g(\dfrac{1}{2}) = 0\) ，\(g(x)\)连续，故有且仅存在一个 \(x\in (0, \dfrac{1}{2})\) ，使得 \(g(x) =\varepsilon\) ，所以唯一性得证。
故存在某个特定的初始距离 \(d_0 =d^\ast\) ，使得小红和小绿永远不相遇。
下面来证明 \(d_0 \ne d^\ast\)的情况，这里也可以分为 \(2\)种情况讨论，\(0 < d_0 < d^\ast\)和 \(d^\ast < d_0 < 1\) 。
\[\phi_{n+1}=g\left(g\left(\phi_n\right)\right)=1-\left(g\left(\phi_n\right)+K\left(g\left(\phi_n\right)\right)\right)=\phi_n+K\left(\phi_n\right)-K\left(g\left(\phi_n\right)\right)<\phi_n\] , also \(g\left(\phi_{n+1}\right)>g\left(\phi_n\right)\).Then
\[\phi_{n+1}-\phi_{n+2}=K\left(g\left(\phi_{n+1}\right)\right)-K\left(\phi_{n+1}\right)>K\left(g\left(\phi_n\right)\right)-K\left(\phi_n\right)=\phi_n-\phi_{n+1}\]
, showing that
\[d_k = d_0 - (d_0 - d_1) - (d_1 - d_2) - \cdots - (d_{k-1} - d_k) \leqd_0 - k(d_0 - d_1)\]
故所求上界 \(k = \left \lceil\dfrac{d_0}{d_0 - d_1} \right \rceil\) ，使得 \(d_k \leq 0\) ，小红和小绿终将相遇。
\[d_k = d_0 + (d_0 - d_1) + (d_1 - d_2) + \cdots + (d_{k-1} - d_k) \ge d_0+ k(d_0 - d_1)\]
2024-02-16 20:32:58
By Long Luo
1 |
/** |
1 |
/** |
1 |
/* @(#)k_sin.c 1.3 95/01/18 */ |
\[\sin (x_0 + \Delta x) \approx \sin (x_0) + \sin'(x_0) \frac {\Deltax}{1!} + \sin''(x_0) \frac { \Delta x^2}{2!} +\sin'''(x_0) \frac {\Delta x^3}{3!} + \cdots\]
