2024-12-21 21:28:05

By Long Luo

2024 年很快就要过去了，就在今天，我们脚下的地球已经以每秒约 30公里的速度飞快越过了近日点，飞驰在围绕着太阳的椭圆轨道上。当 2025年新年钟声敲响时，对于位于银河系第三旋臂边缘的这颗蓝色行星来说，不过是围绕着一颗黄矮星完成了一次再普通不过的公转，正如之前的40多亿次一样。而对于这颗行星上的碳基生物来说，不同的生物感受大不一样，这一刻却意味着对过去366 个日夜 的告别与总结，也承载着对未来的期待与梦想。

和 2023 年不一样的是，我们在 2024 年要经历 366 个日夜，因为 2024年是闰年。小学时，课本和老师都告诉我们一年有365 天，但是闰年却是 366天，这多出来的一天就是 2 月 29 日，在英语中叫做 \(\textit{Leap Day}\)。四年一闰，百年不闰，四百年再闰，这是闰年的规则。后来学习编程时，判断某一年是不是闰年也是常见的编程练习题。在学习之余，你有没有想过，为什么闰年的规则要这么奇怪呢？背后的原因是什么呢？

小学读书时，就想 2 月份很委屈， 1 月 和 3 月都是 31 天，2 月却只有 28或者 29 天，为什么 2 月这么特别呢？为什么有的月份是 31 天，有的月份是 30天呢？但我的老师并没有讲清楚为什么，因为当时我的老师们也不清楚为什么，我也一直到大学里读了一本天文学书才知道了这个问题的答案。

为什么 2024 年 2 月有 29天？要回答这个问题，我们需要穿越历史的迷雾，回顾人类文明史，才能找到答案。

逝者如斯夫，不舍昼夜！

《鲁滨逊漂流记》中的鲁滨逊流落到荒岛之后第一时间就是竖起了一根大柱子，用刀子在立柱上刻上凹口当作日历。一方面是因为鲁滨逊是基督徒需要做礼拜，另外一方面也是为了记录时间。我们的现代生活是离不开钟表的，如果没有钟表来量化时间的话，我们的工作生活将失去秩序。当然鲁滨逊在荒岛上也只能过着“日出而作，日落而息”的农业社会生活，无法精确的安排工作和生活。

“朝菌不知晦朔，蟪蛄不知春秋”，我们智人的寿命足够长，不像朝生暮死的蜉蝣，也不似春生夏死的寒蝉，可以目睹很多生命的诞生、成长以及消亡，感受时间的流逝。“逝者如斯夫，不舍昼夜！”，时间的洪流永远奔涌向前，然而虽然以我们生命的长度可以跨越四季与年轮，但是依然无法触及时间的尽头。

寄蜉蝣于天地，渺沧海之一粟。哀吾生之须臾，羡长江之无穷。当然我们现在知道时间并不是均匀流逝的，也不能脱离物质而存在，但在足够宏大的尺度上，时间在均匀的流逝着。

虽然你可能没有意识到，但是我们一直都在用着天然的时钟，它们就位于我们头顶的星空。这些时钟都足够精准，地球自转的每日误差在毫秒级别，月球公转和地球公转上百年也仅有几毫秒的差别。

“日出东南隅，照我秦氏楼。”，新的一天又开始了；残月如弓，新月如眉，满月如镜，周而复始；“未觉池塘春草梦，阶前梧叶已秋声。”，四季轮回时，我们知道新的一年又来临了。

日长篱落无人过

“或许明日太阳西下倦鸟已归时，你将已经踏上旧时的归途。”每天清晨，太阳都会从东方的地平线上升起来，慢慢爬高，在正午时达到最高点，而后又慢慢落下，在傍晚时分又在西方的地平线落下，千百年来亘古不变。因为地球太大，在古人眼里，我们的大地显然是静止的，不是我们在自西向东转，而是太阳自东向西围绕着我们运动。

日出日落作为一个天然的钟表，人类直到现在依然是依照这个而作息。但一天的时间显然太长，颗粒度不够小，人们迫切需要更小的时间计量单位。在白天观察太阳，初升时太阳很低，影子很长，在正午之前太阳越来越高，影子也随之越来越短。在正午之后太阳越来越矮，影子也随之越来越长，直至日落。而古代文明中心都位于北回归线以北，一年中任何时刻影子都不会消失。

“日长篱落无人过，唯有蜻蜓蛱蝶飞。”，随着季节变化，太阳直射点北移，白天变得越来越长，太阳高度越来越高，影子也随之变短。春分这天昼夜平分，春分过后影子就会产生这样的变化；夏至这天影子最短，白天最长。然后太阳直射点逐渐南移，白天越来越短，太阳高度越来越低，影子也随之变长。冬至这天影子最长，白天最短。如此周而复始，循环往复。

尽管规律复杂，但经过我们的观察记录，我们仍然可以发现太阳的运动轨迹是有迹可循的。太阳在天空中是沿着圆形轨道运动，在一年里日出点和日落点都在变化，但每年都会重复同样的变化，如下图1 所示：

上图位于回归线以北，但纬度是比较低的。越往北走，纬度越高，太阳高度越低，也更冷，如下图2 所示为北纬 \(50^{\circ}\)太阳路径：

如果我们在一个圆盘中心垂直插上一根棍子，然后调整圆盘的方向使得棍子的影子都落在圆盘上，连续记录影子的运动可以发现，太阳的角位移与时间之间存在固定的比例关系，如下图3 所示，于是我们就发明了计时工具：日晷。

世界上很多古文明都独立地发明了日晷这种计时工具，由于不同地区的古代文明计时单位不一样，我们这里就统一为太阳旋转一周需要\(24\) 小时，圆周为 \(360\) 度，于是我们可以得到下列关系：

\[\frac{角位移}{时间} = \frac{360度}{24小时} = \frac{15度}{1小时} =\frac{1度}{4分钟}\]

我们的一天是来自地球自转，日出把我们的时间均匀分割成一个个 \(24\) 小时。

月有阴晴圆缺

当太阳落在地平线之下时，收起了它的光芒，天空变暗，星星和月亮便显现出来。在夜晚，我们很难忽视月亮的存在，各个人类文明都流传着关于月亮的神话故事。我们所在的地球很幸运，和太阳的距离不远不近，又有一颗大卫星—月亮。月亮离我们足够近，在满月时很亮，在日落时升起，日出时落下，这样几乎一整天都有足够的照明时间。

月这个词当然起源于月相变化周期，新月、上弦月、满月、下弦月，周而复始。我们今天使用的公历中月份天数从28 天到 31天不等，月份长度已经和月相周期没有直接关系了。小时候听说“十五的月亮十六圆”，看到日历上的十五，怎么纳闷月亮怎么还是弯弯一枚？后来妈妈告诉我十五是阴历的十五。

在古代，掌握月亮周期是非常重要的，因为如果要在夜间赶路，需要尽量选择有满月照亮道路的夜晚。辛弃疾的词《西江月·夜行黄沙道中》里有一句“明月别枝惊鹊，清风半夜鸣蝉。”说的就是在满月附近那几天赶路的场景。18世纪时在英国的伯明翰诞生了一个日后改变世界的社团：月光社 ( \(\textit{Lunar Society}\))，而这个组织之所以用月亮命名，是因为当时还没有路灯，所以成员们在每月最临近月圆的星期日之夜时聚会，这样散会之后回家也安全方便。而对于生活在海边的人来说，满月和新月时有大潮，上弦月和下弦月时是小潮。

月相周期很容易观察，如果一个月相周期的长度与一天的长度之比是整数的话，那么对于计时来说是非常方便的。但直到太阳潮汐锁定地球之前，地球自转和月亮按照各自的规律运转，两个周期之比不是整数，因为月球一直在远离地球，地球自转也会发生改变，所以两个值时刻都在发生着微小的变化。我们不可能仅仅凭借逻辑分析把地球自转周期和月亮周期联系起来，只能通过仔细的观察和计算发现自然规律的安排，最终得到一个平均值。

伊斯兰中最著名的标志就是星月，因为阿拉伯人因为身处沙漠，白天太热，因而相比其他文明更需要依靠月光晚上在沙漠中穿行，所以月相循环对他们来说非常重要。在8 世纪或 9 世纪时，阿拉伯人已经得到一个相当精确的结果:

\[1个月 \approx 29.5306天 = 29天12小时44分2.9秒\]

这个值精确到了小数点后 4位，我看到书里的资料这么说一直很好奇古代阿拉伯人是如何得到这个精确的数字的，但目前还没有找到答案。这个精度非常重要，伊斯兰历法从公元约622 年开始执行，迄今为止，一直严格地追随月亮的运动，而误差从未超过 1整天。

因为月亮周期约为 29.5 天，所以不管是农历还是伊斯兰历，都分为有 30天的大月和 29 天的小月。

那么为什么目前现行历法中一年都分为 \(12\) 个月呢？这就要从年说起了！

年年岁岁花相似

离离原上草，一岁一枯荣。野火烧不尽，春风吹又生。春天，草木吐绿，万物复苏；夏日，生长旺盛，绿意盎然。秋风渐起，草木渐黄，硕果累累，丰收的季节；冬雪覆盖，大地沉寂，等待新一轮的生命复生。

当远古人类从渔猎生活转变到农业社会，需要具备的能力就是了解季节的变换，尤其是种植谷物的时机。人误地一时，地误人一年，太早或太晚都会影响谷物生长和收成。但是在历法还没发明前，人类是如何判断耕种时机呢？最简单的方法是观察植物的生长周期，比如把草木变绿开花当作一年的开始。不过每年植物开花的时间跟气候有关，降雨、气温和阴晴日照都可能会影响花开的时间，用这种方法并不准确。

我们现在知道季节变迁来自于我们地球公转以及地轴倾角所示，如下图 7所示。地球公转如此稳定，在上千年的尺度上，误差也就 2 分钟左右。

地球公转一周的时间我们称为一年，但一年长度与一天的长度之比并不是整数，实际关系是:

\[1回归年 \approx 365.2421990741日 = 365天5小时48分46秒\]

这个比值其实是个无理数，所以我们只能得到一个大致精确值，永远得不到具体值。令人惊讶的是古埃及早在大约公元前2650 年也就是 \(4650\) 年前就得到了一年为 \(365.25\)天，简直就是个奇迹。古埃及历法把一年定为严格的 \(365\) 天，分为 \(12\) 个月，每个月 \(30\) 天，剩余 \(5\) 天作为神的节日。而古埃及人得到 \(365.25\)天这个值并不是通过逻辑计算，而是通过观察天上的星星得到的。

天狼星，是夜空中最明亮的一颗星，这让它在各地的古文明中都充满有趣的故事，冬季大三角中最低最亮的那颗星就是天狼星。

由于地球公转，恒星每天在天上的位置都会有点不一样。因为太阳光实在太亮了，所以当恒星在天上的位置与太阳太近的时候，就会一整天都看不到那颗星。而恒星在一年当中刚好赶在日出前夕从东方升起的日子称为偕日升。在这一天之后，天狼星每天都比太阳早4 分钟从东边升起，每天都能在夜空中多闪耀 4分钟，所以偕日升是一颗恒星能被看到最短的日子。

在古埃及时代，当天狼星在拂晓出现在东方地平线上的话，就代表尼罗河将要泛滥，带来大地肥沃的泥土，提供农作物丰富的养分，让人们得以耕作。这一天大概是现在的每年的7 月 19日左右，也因此天狼星被用来作为古代埃及历法的标准，作为新一年的开始。

古埃及人发现大约 \(4\)年之后，天狼星就会晚一天偕日升，认识到了每 \(4\)年应该加一个闰年。但是古埃及人将天狼星偕日升这个日子会每 \(4\) 年推后一天，这样在 \(365 \times 4 + 1 = 1461\)个埃及年之后，太阳会再次在年初和天狼星偕日升，也就得到了地球公转周期为：\(365 \times \frac {1461}{1460} =365.25\) 。

这些偕日升的星星，其实跟地球的季节变化没有因果关系，它们只是在对的时间，出现在对的位置。不过为什么不同地方的人会发展出类似的方法呢？主要是这个方法简单好用，不需要精确的算出太阳运行周期，只要早起看星星就行。

古代科技不发达，天文学成为了人们理解自然和预测未来的重要工具。会看星星的人类祖先才能建立农业社会，可以比其他人占优势。农耕稳定生产粮食，足够粮食能养活更多人，这群人再建立社会制度，发展成强大的文明。其实，我们都受益于会看星星的祖先，会看星星的人比较容易存活下去。

一年为什么分为 12 个月？

大自然为人类提供了三个各自独立的计时标准：日、月和年，但把这三个标准统一起来却并不是件容易的事。这三者两两比例都是无理数，我们只能得到一个大致的精确值。

我们已经知道一年约为 \(365.24\)天，一月约为 \(29.5306\) 天，那么 \(\dfrac {365.24}{29.5306} = 12.3682\)，但是余数 \(0.3682\)怎么处理呢？直接忽略的话，几十年的尺度里可能无所谓，但是几百年下来就会偏差特别大，必须修正。那么如何修正这个值呢？

从数学角度来看，一年约为 \(365.24\)天，而 \(365 = 5 \times 73\) ， \(366 = 2 \times 3 \times 61\)，但遗憾的是， \(73\) 和 \(61\) 都是质数，而且是大质数。人类天生就不喜欢大质数，因为不好均分，要是一年是 \(360\) 天就好了：

\[360 = 2 \times 2 \times 2 \times 3 \times 3 \times 5 = 12 \times 30\]

这样可以被 \(2, 3, 4, 5, 6, 8, 9, 10, 12\dots\) 整除，这也是古巴比伦人为什么将一个圆周分为 \(360\) 度的原因。

不过我们生活地球并不是这样，但是几乎所有历法 1 年都是分为 12个月，大概还是因为月份是质数的话，一年就无法均分了。对于余数部分则打上不同的补丁来修正，比如各种历法在某些月份增加天数或者增加某个月份来修正。

为什么有闰年？

闰，本意是余数，闰年，英语中叫 Leap Year ，而 Leap是跳跃的意思，这个含义是怎么来的呢？

之前我们已经知道我们定义地球绕太阳运动一圈的时间就是一年，这个回归年约为\(365.24219\)天。地球公转极其稳定，其实每一年都是一模一样的时间，本质上闰年也不会多一天。但是显然把一年的天数设置为整数更加方便，不然谁记得住呢！四舍五入取整之后的一年便是\(365\) 天，而实际上每年都会多出 \(0.24219\) 天， \(4\) 年下来就会多出 \(0.24219 \times 4 = 0.96876\)天，约等于多出一天，看下面几张图会理解的更清楚一些。

第 \(1\) 年，我们会延后 \(0.24\) 天：

第 \(2\) 年，我们会延后 \(0.48\) 天：

第 \(3\) 年，我们会延后 \(0.73\) 天：

第 \(4\) 年，我们会延后 \(0.97\) 天：

当到了第 \(4\)年时，为了和真实时间看齐，我们要把这一天补回来，于是制定了闰年加一天的规则。这时可能有人会问，既然是为了提高生活效率，那直接不补每年天数都一样多岂不是更加方便。但是他们明显小看了每\(4\)年差一天的威力，积少成多之下会产生巨大影响。在不补的情况下每过 \(100\) 年，就会相差 \(100 \times 0.24219 = 24.219天\)，这已经接近一个月的差距了！农作物差上一个月，收成就别指望了！

所以闰年必须有，哪怕每四年只差不到一天，也必须补回来。但是为什么又会\(100\) 年不闰， \(400\) 年再闰的规则？为什么 \(100\) 年时的那一天为什么又不用补了呢？

根据之前的数据，我们可以得知每 \(4\)年就会比实际时间少了 \(0.96876\)天，但是闰年是直接补充一整天，这就导致了在出现了闰年规则后，每 \(4\) 年反而多出了 \(1 - 0.96876 = 0.03124\) 天。这样每过 \(128\) 年就会多出了 \(0.03124 \times 128 \div 4 = 0.99968\)天，为了方便计算，在每一百年的时候再减掉一天，也就是整百年份的时候把闰年变回普通年。

同样再做一次更加精确的计算，每一百年就会比实际的时间少了 \(1 - 0.03124 \times 100 \div 4 = 0.219\)天，那么每四百年就少了 \(0.219 \times 4 =0.876\) 天，那就再把这一天补回来，也就是年份为 \(400\)的倍数的时候又需要再一次变成闰年。所以目前闰年的完整规则为什么会这么复杂，年份不能被\(4\) 整除时不是闰年，年份能被 \(4\) 整除但不能被 \(100\) 整除时是闰年；年份能被 \(100\) 整除但不能被 \(400\) 整除时不是闰年；年份能被 \(400\) 整除时是闰年。

其实由上面的计算可以知道其实每 \(400\) 年还是多了 \(1 - 0.876 = 0.124\) 天，所以每 \(4000\) 年还会多出 \(1.24\)天，也许到时候会再打上更细致的补丁。今年才公元 2024 年，离公元 4000年还有 1000 多年，也许人类不久之后换了历法也不一定。

回到开头的问题，因为我们通过额外添加天数，所以使得日期产生了跳跃，这就是Leap Year 的由来。

为什么 2 月有 29 天？

目前国际通用的历法是公历，也叫格里历，大家也可能叫阳历，但是公历只是阳历的一种。小时候不太懂什么阳历阴历的，阳历是根据地球呈现出太阳直射点的周期性变化所制定的历法，而阴历则是根据月亮的周期。

之前月相的章节提过，对于阿拉伯人来说，月亮非常重要，伊斯兰历是一种纯粹阴历，也是现今仍在使用的历法中唯一的纯阴历。中国现行的农历则兼顾了太阳周期和月亮周期，属于阴阳合历，是400多年前明朝由欧洲传教士汤若望和中国第一批天主教徒合力编撰的。犹太历也兼顾了太阳周期和月亮周期，和农历一样通过闰月实现和实际的太阳年一致。

一年有 365 天，平均分为 30 天，那么只需要 5 个月是大月，也就是 31天，而现在的历法是 1 月，3 月， 5月，7月， 8月， 10月，12月都是 31天，而平年 2 月只有 28 天？为什么 7 月 和 8 月都是 31天呢？要回答这个问题，我们需要回到回顾过去，回到古罗马时期，才能找到答案。

旧罗马历法时一年只有 10个月，比现行的公历少了2个月，冬天是没有月份名称的。旧罗马历的月份算法从春天的March（3月）开始，持续到年底的December（12月）。（古罗马认为偶数不祥，为了避免偶数不吉利，罗穆卢斯为大部分月份设定了奇数天数（29或31天），但为了全年天数平衡，2月被赋予28天，成为唯一的偶数月份。2月的特殊地位还和罗马的宗教传统有关：2月被视为祭祀亡灵和赎罪的月份，因此是一个“缩短”的月份。

古罗马时期，罗马历法在政治目的的影响下经常遭到任意篡改。由于教皇通常是政客，而且罗马行政长官的任期与日历年相对应，这种权力很容易被滥用：教皇可以延长他或其政治盟友在任的年份，或者拒绝延长其对手在任的年份。到了凯撒大帝执政时，历法已被改得面目全非，为了解决这个问题，他采纳埃及托勒密王朝亚历山大的希腊数学家兼天文学家索西琴尼计算的历法，按照太阳周期制定了新的历法，也就是儒略历，在公元前45年1月1日起执行，以取代旧罗马历法。新历法改为12个月，但仍然保留了之前历法月份天数不同的规则。七月Julius，原名Quintilis，是“第五”的意思。由于凯撒在该月出生，经元老院一致通过，将此月改为凯撒的名字“Julian”。

之后凯撒的继承人屋大维，是最伟大的罗马皇帝之一。8月，原名Sextilis，旧罗马历中第六个月，屋大维在建立罗马帝国的过程中打的几次辉煌胜仗都发生在8月，其中也包括攻陷亚历山大港、征服埃及之役。Augustus本是元老院赐封给他的称号，意为“神圣伟大”，从此屋大维也被称为奥古斯都Augustus。后来屋大维死于此月，元老院于是将此月改为他的称号“奥古斯都”。

因为旧的儒略历每四年都设置了一个闰年，通过上一节的计算我们知道时间足够长的话，每100 年会偏差达到将近一天，1500年之后偏差就会达到 10天。16世纪时，春分已经向前移动了10天。春分实际上在3月21日，在日历上却变成了3月11日。对于基督徒来说，这是个严重的问题，因为当时一年中最重要的节日是复活节，而复活节的日期是依靠春分确定的：春分后的第一个满月后的第一个星期天。

为了使复活节回到正确的日期，教皇格里高利十三世于1582年规定，把日历中的日期后移10天，同时把闰年的规则修订为现行的规则，以确保这样的问题不再出现，这就是我们现行公历格里历的由来。

总结

2024 年即将结束，这一年里我们经历了 366天。2024年2月有29天，这并不是偶然的结果，而是历经历史演变和历法改革的产物，是人类对自然规律更深刻的认识，是人类智慧和文明进步的体现。

参考文献

1. 天文学Astronomy
2. LeapDay
3. 闰年 LeapYear
4. 季节 Season
5. Daniel Defoe:《鲁宾逊漂流记》
6. EricChaisson/Steve McMillan: Astronomy: 《天文学与生活》
7. Michael Hoskin:《The Cambridge Illustrated History of Astronomy剑桥插图天文学史》
8. CamilleFlammarion: 《大众天文学》
9. 太阳 Sun
10. 太阳路径 SunPath
11. 计时工具的历史Timekeeping devices
12. 日晷Sundial
13. 日晷的历史History of Sundials
14. 月亮 Moon
15. 月相 MoonPhase
16. 月光社Lunar Society
17. 潮汐 Tide
18. 星月
19. 伊斯兰历Islamic calendar
20. 年 Year
21. 天狼星Sirius
22. 天狼星周期Sothic cycle
23. 偕日升Heliacal rising
24. AReal Scorcher! — Sirius At Heliacal Rising
25. 儒略年Julian calendar
26. 公历Gregorian calendar
27. 阳历 Solarcalendar
28. 回归年Tropical year
29. 农历Chinese calendar
30. 利玛窦 MatteoRicci
31. 汤若望Adam Schall
32. 犹太历Hebrew calendar
33. 尤利乌斯·恺撒Julius Caesar
34. 7月 July
35. 8月 August
36. 屋大维 Augustus38 Interplanetary:This is why we have Leap Years
37. 2月February
38. WhyFebruary has only 28 days? Nothing more scientific than Romansuperstition…

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}x \mathrm{d}y = r\mathrm{d}r\mathrm{d}\theta\) ？

在直角坐标系中，面积微元 \(\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}\]

3D 钟形曲面下的体积

转换为极坐标法是常用的求高斯积分的方法，尤其是对高斯积分升维更是神来之笔，那么升维之后的高斯函数是什么样的呢？

\(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}\]

沿着 \(y\) 轴方向进行切片

我们通过对高斯积分进行升维得到了一个 \(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\)维钟形曲线等概率事件是在同一等高线形成的圆上，这就是圆的由来。

参考文献

1. 高斯积分Gaussian Integral
2. 高斯函数Gaussian function
3. 圆周率 Pi
4. 笛卡尔坐标系
5. 极坐标
6. 直角坐标与极坐标的互化中，为什么dxdy=rdrdθ？
7. Why π is inthe normal distribution (beyond integral tricks)
8. TheGaussian Integral, a Geometrically Annotated Proof

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}\)）把已知太阳系中行星的数据列成如下一张表，如下表所示：

他发现行星与太阳之间的距离大致满足以下经验公式：

\[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\)个数据就计算出了谷神星的轨道。

他的最终计算是指向了一个完全不同的天空区域，之前被其他科学家所忽视的区域。从上图2也可以看出谷神星并不在黄道平面上，而是和黄道平面有个较大的夹角，所以谷神星才不容易寻找。

1801 年 12 月 31 日夜，德国天文爱好者奥伯斯( \(\text{Heinrich Olbers}\) ,1758-1840)，在高斯预言的时间里，用望远镜对准了这片天空。果然不出所料，谷神星出现了！

高斯计算谷神星轨道的方法比较复杂，可以写一篇上百页的论文，这里提供一篇论文HowGauss Determined The Orbit of Ceres 供大家参考。

不同于其他科学家把皮亚齐的观测数据当作真实数据，高斯一开始就考虑到皮亚齐的观测数据误差问题，毕竟当时的观测设备精度有限，所以求得的轨道数据应该是一个范围。如下图3 所示就是高斯计算出谷神星轨道的原始手稿：

高斯使用的数据分析方法，就是以正态误差分布为基础的最小二乘法，如下图4 所示：

高斯认为这些数据点使用最小二乘法计算出的回归曲线，如下图 5所示，等价于这些点沿着这条回归曲线呈正态分布。

那 么高斯是如何推导出误差是正态分布的呢？我们将在下一节进行推导。

误差是怎么分布的？

天文学是人类历史上最古老的科学之一，其起源可以追溯到史前时代。古代文明通过观察天空中的天体运动来制定历法、导航和进行农业活动。天文学家不光对天体进行了长期的系统观测，记录了大量的天文现象，而且试图预测不同行星和恒星的轨迹和位置。

由于观察者、仪器和许多其他因素，天文学家的测量数据存在误差。即使对于同一个物体，他们也有不同的观察结果。这里误差指的是观测值与实际值之间的差异。千百年以来，对于有误差的测量数据，多次测量取算术平均是比较好的处理方法。虽然这种方法缺乏理论上的论证，也不断的受到一些人的质疑，但取算术平均作为一种异常直观的方式，在多年积累的数据的处理经验中也得到相当程度的验证，被认为是一种良好的数据处理方法。

伽利略在 1632年的著作《关于托勒密和哥白尼两个主要世界体系的对话》中，非正式地讨论了他所谓的“观测误差”，也就是今天所谓的“随机误差分布”，对误差的分布做过一些定性的描述，主要是以下2 个假设：

1. 观测值围绕真实值对称分布，也就是说，误差对称分布在零左右；
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}\]

总结

高斯通过对观测误差的分析，寻找误差的概率分布密度函数。大胆假设最大似然估计的结果是算数平均，而在所有的概率密度函数中，唯一满足这个性质的就是正态分布。

参考文献

1. StackExchange:How was the normal distribution derived?
2. 哲学家孔德
3. 棣莫弗
4. 拉普拉斯
5. 博得法则：TheTitius-Bode Formula
6. 小行星带
7. 谷神星Ceres
8. 天文学家朱塞普·皮亚齐
9. 高斯Carl Friedrich Gauss
10. 高斯计算小行星的方法
11. Gauss,Least Squares, and the Missing Planet
12. 最小二乘法
13. 最小二乘法 (LeastSquares)
14. 最小二乘法的本质是什么？
15. 误差
16. 概率密度函数Probability density function
17. 最大似然估计
18. 如何通俗地理解“最大似然估计法”?
19. HowDid Gauss Derive The Normal Distribution
20. 靳志辉：正态分布的前世今生(上)
21. 靳志辉：正态分布的前世今生(下)
22. 米勒：《普林斯顿概率论读本》
23. 杰恩斯：《概率论沉思录》
24. 格林伯格：《普林斯顿数学分析读本》

2024-04-27 13:16:57

By Long Luo

相信大家或多或少都听过六西格玛( \(\text{6Sigma}\) ) ¹ 这个词，六西格玛是指生产的产品中，\(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}\) ) ² 中， \(68.27\%\) 的数据位于平均值的一个标准差内，\(95.45\%\) 位于两个标准差内， \(99.73\%\) 位于三个标准差内，这也是著名的68-95-99.7 Rule ³ ，如下图 1 所示：

什么是正态分布？

数据可以用不同的方式“分布”，比如数据可以向左散布的多一些，也可以向右散布的多一些，或者分布的乱七八糟，如下图2 - 4 所示，

但数据经常会集中在一个中心值的附近，而不向左或右偏斜，像一个钟形，如下图5 所示。

正态分布，又称高斯分布（ \(\text{GaussianDistribution}\) ），是一种重要的概率分布，数学王子高斯 ⁴在正态分布的研究和应用上做出了巨大贡献。有很多日常现象都符合这种分布，如人的身高、考试成绩等。正因为它几乎无处不在，所以叫\(\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}\) ） ⁵ 在 1850 年给出的。3Blue1Brown 的视频 Why π is in thenormal distribution (beyond integral tricks)中详细介绍了这种方式。不过视频中有一些不够严谨的地方，下面会先介绍视频中的推导方法，然后再介绍严谨的数学分析法。

考虑向一个镖盘投掷飞镖，过镖盘中心作 \(x\) 轴和 \(y\)轴。每次投掷飞镖都会受到随机因素而偏离目标，故每次飞镖在镖盘的落点 \((x, y)\) 是 \(2\) 维随机变量。

假设满足以下 \(2\) 个条件：

1. 落点的 \(x\) 轴和 \(y\) 轴坐标是相互独立的；
2. 落点的概率密度函数仅与其到原点的距离有关，即分布在空间上具有旋转对称性。

3Blue1Brown Herschel 推导方法

如下图 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}\) 是一个函数方程 ⁶，熟悉函数方程的同学可能一眼就知道满足公式 \(\eqref{6}\)的函数解是指数函数 ⁷ 。

\[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}\) ） ⁸ ，如何求解高斯积分可以参考这篇文章：数学之美：几何视角下的高斯积分(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}\) ） ⁹有关，我会在下一篇文章详细解释，敬请期待！

参考文献

2024-04-21 08:08:32

By Long Luo

世界那么大，我想去看看！随着科技的发展，我们早已做到无需离开家就能领略世界各地的风景和文化。通过高清视频、高清直播、社交媒体、VR技术、各种图片或者视频分享平台，我们不仅可以体验世界名胜古迹、自然美景和各地的风土人情，还能与当地居民互动，了解他们的日常生活和传统文化。互联网丰富了人们的生活，缩小了地域的界限，真正实现了让世界触手可及，足不出户便可周游世界的梦想。但正所谓“百闻不如一见”，“读万卷书，不如行万里路！”，尽管互联网让我们可以虚拟游览世界，但亲自出行的体验无可替代。自由行不仅提供了前所未有的自由和灵活性，还能让我们亲身感受到大自然的美妙、城市的活力。这种身临其境的体验，远非屏幕前的感受可比。

当你打开地图软件时，你会看到如图 1 所示的道路标志，

当你在道路出行时，你也会看到看到如下图 2所示的路牌，但你可能并未真正留意过这些标志。因为现在我们只需要有一部联网的智能手机，在地图类软件里，设定出发地和目的地，自然有导航会指引我们到达目的地。

这些编号肯定不是随机的，那么这些道路编号到底有什么用呢？出于好奇心你可能会去寻找答案，你很容易轻松找到中国国家高速的编号密码 ¹ 和 高速公路是怎样命名和编号 ²这样的文章。但这些文章只是告诉我们是什么，并没有告诉我们为什么。

国内现行的高速公路命名是由交通部从 2005 年启动的 ³，之前道路都以起始地和终点地命名。你可能会想，之前那种命名方式不是更合理吗？用了数字不是更加不清晰易懂吗？如果图1 不是路牌上写了济广高速，谁知道 G35 ⁴ 是哪里到哪里呢？

如果你更进一步的话，如果你去了解世界其他国家的高速公路命名的话，你会惊讶的地发现为什么居然全世界各主要大国都选择了类似的编号系统，这背后的原因是什么呢？

要回答这个问题，我们需要把时钟拨回几十年前，回到高速公路诞生的时期，那个没有GPS，没有手机，只有纸质地图的时代，我们才能知道这种编号系统的重要意义和实用性，以及背后的数学密码。

世界各国如何对高速公路进行编号？

我们已经了解了国内高速公路编号 ⁵，让我们看看其他国家的高速公路系统编号是什么样的。这些国家需要国土面积足够大，高速公路系统足够发达，国土疆域长宽比例没有太夸张，人口分布比较均匀。

德国 （Germany）

首先我们把目光投向德国 ⁶，因为世界上第一条被广泛认可的现代意义上的高速公路就是德国的“Autobahn”⁷ ，于 1932年建成，连接科隆与波恩，完全是为高速行驶的机动车辆设计和建造的。

在纳粹上台之后，高速公路建设被视为国家的重要工程，是当时德国最具标志性的基础设施项目之一，以展示德国的工程实力和国家威严⁸。高速公路建提高了交通运输效率和便利性，同时也为了支持军事行动和纳粹的宣传活动。这些公路被设计成宽敞平坦、直线且无障碍的道路，可以实现高速行驶。

这些高速公路被称为“Autobahn”，德语的意思就是“汽车道”，命名方式则采用了简单的字母和数字组合，例如\(A1, A2, \cdots , A9\) 等 ⁹。这种命名方式简洁明了，便于驾驶员识别和导航。

德国现在高速公路 (Bundesautobahn，BAB) 的编号系统是在 1974年开始使用的，所有的高速公路都以 \(A\)开头后接一个空白与数字编号，如 \(A8\)。穿越德国全国的东西向主要高速公路以偶数编号，南北向的道路则以单数编号。用来连接区域性重点城市的较短高速公路则以两位数字来编号，具体如下图3 所示：

• \(A1 - A9\) :主要南北向高速公路，数字越小越西。
• \(A10 - A99\) :主要东西向高速公路，数字越小越北。
• \(A100 - A999\) :环路、支线和次要高速公路。两位数和三位数的编号通常用于城市环线和连接支线。

当然德国也接入了欧盟标准的“E”编号系统与其他欧洲国家的公路系统协调，如下图 4 所示：

法国 （France）

法国的高速公路系统 ¹⁰被称为“Autoroute”，通常以字母“A”开头，后跟一个数字。编号体系如下：

• \(A1 - A16\)沿着首都巴黎的放射性高速公路；
• \(A100 - A999\) :区域性或次要高速公路。

由于巴黎长期是法国的中心 ¹¹，所以法国的高速系统沿巴黎放射性分布，如下图 5所示。当然法国是欧洲国家，所以也使用了欧盟标准的“E”编号系统与其他欧洲国家的公路系统协调。

俄罗斯 （Russia）

俄罗斯的高速公路系统 ¹² 被称为 Avtomagistral，俄语里就是高速公路的意思，编号使用字母“M”、“A”或“P”开头，后跟数字：

• \(M\) 系列:连接首都莫斯科与其他重要城市和边界。
• \(A\) 系列:连接各主要城市，通常为区域性主干道。
• \(P\) 系列:其他重要的区域性公路。

俄罗斯和法国类似，精华都在莫斯科圣彼得堡的东欧平原区，西伯利亚和远东地区人口稀少，高速分布如下图6 所示：

欧洲 （Europe）

欧洲的高速公路系统 ¹³使用“E”字母（Europe）作为前缀，后接一个数字，代表着欧洲主要道路网。

欧洲路线编号的决定是由 1975 年的欧洲经济委员会作出的，并在 1992年进行了修正。作为干线的A级道路编号是两位数，而作为地方线的B 级的道路编号则是三位数，主要规则如下所示：

• 南北方向的道路按两位数编号，从西向东依次递增，以 \(5\) 结尾。
• 东西方向的道路按两位数编号，从北向南依次递增，以 \(0\) 结尾。

当然还有其他一些规则，欧洲高速公路分布如下图 7 所示：

由于加拿大 ¹⁴ 和 巴西 ¹⁵人口分布不均匀，没有什么借鉴价值。当然有同学看到这会问，欧洲高速公路系统和国内编号规则几乎一致，那么国内是借鉴了欧洲的公路编号系统吗？

其实世界各大国的高速公路编号体系在设计上几乎一致，这种统一编码方案可以追溯到美国的州际高速公路¹⁶系统，这个我们会在下面的章节里进行详细讲解。

我们可以先问下自己几个问题：

1. 欧洲这种编号规则有什么好处？为什么要这么编号？
2. 如果让你来对公路进行编号，你会怎么编号呢？
3. 2个城市之间可能有多条线路，又该如何编号呢？

如何给道路进行编码？

假设我们在下图 8 所示的纽约曼哈顿的 \(A\)点，我们要去古根海姆博物馆。此时我们位于 东区第 87 大街和 第 3大道的十字路口，我们要去的目标在 东区第 88 大街和 第 5大道交汇处十字路口。很显然，我们有很多条路线可以到达目的地。

即使你是个路痴，当你沿着东区第 87大街走的时候，如果你走着走着发现走到了第 4大道，你就知道至少在东西方向上你的方向对了，而当你发现走到了第 2大道，你就知道你方向走反了。同理，对于南北方向上也是如此。这样，即使你绕了很多路，但你只要知道目的地所在街区的编号，你就调整方向，多走点路多花点时间也一定能走到目的地。

当然我在写这篇文章时，想找个浅显易懂的例子，结果发现曼哈顿街区完美的符合我想找的例子。因为道路规划太好了，呈棋盘状。实际中大部分城市街区命名不一定使用数字，也不一定如此横平竖直。但是曼哈顿街区因为都使用数字编号，我们想找具体地点太容易了，这还没用上二分查找呢！

由于现实中由于不同地区地理条件不一样，城市的分布往往在地图上是随机的，比如下图9 就是我随手绘制的一个道路分布图。

如何对这些道路进行命名？你可能会说，那就以出发地和目的地命名好了，比如\(AB\) , \(BC\)等。当然这样命名也可以，但是我们可以更进一步，让命名和曼哈顿街区命名一样，让名字就自带导航功能。

那比如有 100 条道路，如何给每条道路编上正确的号码呢？从 1开始递增当然也可以，但就没有任何信息了。

让我们回顾现实中的高速公路，我们知道一条高速公路有哪些属性呢？

1. 是东西走向？南北？还是东北、东南、西南、西北？
2. 是国家级主干道还是省级，市级？
3. 连接哪些城市？在哪个位置？

还记得地球仪上的经线和纬线吗？墨卡托投影制作的地图帮助无数航海家顺利抵达目的地，航海家只要知道测量出当前所在地经纬度，沿着连接目的地的方向就不会迷失方向。

想一想，如何数字也有很多属性，比如奇偶性、大小、位数等，参考之前的编号系统，我们可以用数字属性和高速公路属性一一对应¹⁷ ，比如：

1. 奇偶 \(\leftrightarrow\)公路走向
2. 大小 \(\leftrightarrow\)公路位置
3. 数字位数 \(\leftrightarrow\)主干道和次要干道
4. \(5\) 的倍数 \(\leftrightarrow\) 重要道路

不得不说，这种设计比曼哈顿的纯数字编号要好很多，绝对是天才般的设计 ¹⁸！下面我们就来介绍这种编号系统的诞生地：美国州际高速公路系统。

州际公路（InterstateHighway System）的天才设计

也许是艾森豪威尔总统 ¹⁹二战时担任盟军统帅，战后又担任美国驻德国占领军总司令的缘故，我想当时德国领先世界的高速公路系统肯定给他留下了深刻印象。在他1953 - 1961年担任总统期间，提出了州际高速公路系统的构想，并在1956年签署了《国家公路与交通法案》，旨在建立一套横贯美国的高速公路网络，以提升国家交通运输效率和安全性。这个宏伟的工程由联邦政府主导，并联合各州政府共同筹资和建设。经过数十年的不懈努力，美国州际高速公路系统于1992年正式完成，涵盖了数万英里的公路，连接了各个州份和城市，成为美国国家交通运输体系的重要组成部分，也促进了经济的发展和人口的流动。

州际公路系统每条线路和编号，可以参考这个在线网站：Interstate Highway Map

下面我们来介绍下州际高速公路的编号规则 ²⁰：

2 位数字

1. 偶数且为 \(5\)的倍数的 \(2\)位数字表示东西向主干道，从南到北依次递增，如下图 10 所示：

其实也很好记， \(\text{Even}\) 对应\(\text{East}\) ，而且 \(0\)不偏不倚，表示水平方向，很容易记吧！

2. 奇数且为 \(5\)的倍数的 \(2\)位数字表示南北向主干道，从西到东依次递增，如下图 11 所示：

其实也很好记，偶数表示东西向，那奇数就是南北向了，这是主要干道，肯定不能依次递增，所以每次递增\(5\) 。注意西海岸的 \(5\) 要当成 \(05\) 。

如何理解数字增加方向呢？

由于美国本土大致呈矩形，如果我们在地图左下角作为坐标原点，建立一个笛卡尔坐标系，如下图12 所示，想必你忘不了：

由于中国的国土疆域并不像美国接近一个矩形，所以坐标点实际是放在靠近北京的右上角，如下图13 所示：

如果你掌握了这个规则，那么下面这个路标，你能猜出大概在美国哪里吗？可以先猜测然后自寻答案验证下你的猜测。

3 位数字

相比 \(2\)位数字高速公路是跨州线路(注：不完全如此)， \(3\)位数的州际公路是服务于各个都会区的较短路线，它们连接到更长的两位数路线，并且充当环城公路、支线或连接线路。

\(3\) 位数编号规则如下：

• 第 \(1\) 位数字反映了道路的用途，后\(2\)位数字反映路线连接到的任何两位数的州际公路。例如，I-395 连接到I-95，I-270 连接到 I-70 ；
• 第 \(1\)位数字如果是偶数，说明和州际公路会有 \(2\) 个交点，奇数则仅相交一次；
• 对于 \(3\)位数的州际公路，只要同一数字不在同一州内重复，就可以根据需要重复相同的数字。

国内公路编号系统则大致相同，可以参考 主干道 ²¹和 辅道 ²² ，具体可参考 中国高速公路网 。

上面就是美国州际高速公路的大致规则，当然也会有一些例外。

这套编号系统不仅方便导航，有效帮助驾驶员辨别方向，还提供了关于道路走向和地理位置的信息，提高了道路网络的整体效率和安全性，它是如此的成功，以致于后来陆续被其他国家所采用。

总结

生活处处皆学问，正如王建硕老师写的系列文章：无用的冷知识。如果我们对日常事务多一点好奇心，我们可以发现后面更多的原理！

参考文献

2024-04-16 11:13:48

By Long Luo

阿里巴巴达摩院 从2018年开始每年都会举办一届全球数学竞赛，之前一方面自己数学水平比较弱，另外一方面也没有报名，但一直很仰慕那些数学大神的风采。今年是第一次报名参加2024阿里巴巴全球数学竞赛，上周末参加了预选赛，但遗憾的是，全部 \(7\) 道题中只有第 \(1, 2, 6\) 题会做，这里分享下我的解答：

Problem 1

几位同学假期组成一个小组去某市旅游. 该市有 \(6\) 座塔，它们的位置分别为 \(A, B, C, D, E, F\)。同学们自由行动一段时间后，每位同学都发现，自己在所在的位置只能看到位于\(A, B, C, D\) 处的四座塔，而看不到位于\(E\) 和 \(F\) 的塔。已知：

1. 同学们的位置和塔的位置均视为同一平面上的点，且这些点彼此不重合；
2. 塔中任意 \(3\) 点不共线；
3. 看不到塔的唯一可能就是视线被其它的塔所阻挡，例如，如果某位同学所在的位置\(P\) 和 \(A, B\) 共线，且 \(A\) 在线段\(PB\) 上，那么该同学就看不到位于 \(B\) 处的塔。

（5 分） 请问 这个旅游小组最多可能有多少名同学？

\(A. 3\)
\(B. 4\) \(C.6\) \(D. 12\)

Solution

这道题选 \(C\) ，最多只能有 \(6\) 名同学。

1. 假设有 \(1\) 座塔 \(A\) ，那么很显然有无数多同学可以看到塔\(A\) ，也可以有无数多同学看不到塔\(A\)​ ；

2. 假设有 \(2\) 座塔 \(A, B\) ，那么只有以 \(A\) 为起点的射线 \(AB\) 且位于 \(B\) 之后的同学无法看到塔 \(A\) ；

3. 假设有 \(3\) 座塔 \(A, B, C\)，同理可知存在无数位同学至少可以看见 \(2\) 座塔；

4. 假设有 \(4\) 座塔 \(A, B, C, D\)，同理可知存在无数位同学至少可以看见 \(2\) 座塔；

5. 假设有 \(6\) 座塔 \(A, B, C, D, E, F\) ，如果每位同学都无法看见\(E, F\) 塔，如下图1 所示：

所以至多有 \(6\) 位同学位于 \(M, N, O, P, R, Q\) 处，无法看到塔 \(E, F\) 。

Problem 2

小明玩战机游戏。初始积分为 \(2\)。在游戏进行中，积分会随着时间线性地连续减少 (速率为每单位时间段扣除\(1\) )。游戏开始后，每隔一个随机时间段(时长为互相独立的参数为 \(1\)的指数分布)，就会有一架敌机出现在屏幕上。当敌机出现时，小明立即进行操作，可以瞬间击落对方，或者瞬间被对方击落。如被敌机击落，则游戏结束。如小明击落敌机，则会获得\(1.5\)个积分，并且可以选择在击落该次敌机后立即退出游戏，或者继续游戏。如选择继续游戏，则须等待到下一架敌机出现，中途不能主动退出。游戏的难度不断递增：出现的第\(n\) 架敌机，小明击落对方的概率为\((0.85)^n\) ，被击落的概率为 \(1 - (0.85)^n\)，且与之前的事件独立。在任何时刻，如果积分降到 \(0\) ，则游戏自动结束。

第 1 问

小问 1 (5分)如果游戏中，小明被击落后，其之前的积分保持。那么为了游戏结束时的累积积分的数学期望最大化，小明应该在其击落第几架敌机后主动结束游戏？

\(A. 1\) \(B. 2\) \(C.3\) \(D. 4\)

Solution

根据小概率事件原理，我们可以把打飞机事件视为泊松过程¹，那么事件分布就是泊松分布² 。其 Python代码可以直接调用 API numpy.random.exponential。虽然是 指数分布，但在 Java 中 需要使用 -Math.log(1 - random.nextDouble())而不是 Math.exp(double a) 。模拟代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

public class WarPlaneGame {

    private static Random random = new Random();

    private static double getExpectScore(int planes) {
        double score = 2;

        for (int i = 0; i < planes; i++) {
            double waitTime = -Math.log(1 - random.nextDouble());

            if (waitTime >= score) {
                score = 0;
                break;
            }

            score -= waitTime;

            double possonPr = random.nextDouble();

            double shootDownPr = Math.pow(0.85, i + 1);

            if (possonPr < shootDownPr) {
                score += 1.5;
            } else {
                break;
            }
        }

        return score;
    }

    private static double simulate(int planes) {
        int simulateTimes = 50000;

        double scoreSum = 0.0;

        for (int i = 0; i < simulateTimes; i++) {
            scoreSum += getExpectScore(planes);
        }

        return scoreSum / simulateTimes;
    }

    public static void main(String[] args) {

        for (int i = 1; i < 5; i++) {
            double result = simulate(i);
            System.out.println("Shoot down " + i + " planes, Expect Score: " + result);
        }
    }
}

输出如下：

1
2
3
4

Shoot down 1 planes, Expect Score: 2.2344824425425442
Shoot down 2 planes, Expect Score: 2.290207012168609
Shoot down 3 planes, Expect Score: 2.2653361024420042
Shoot down 4 planes, Expect Score: 2.187342196221392

可以看出击落第 \(2\)架敌机后主动结束游戏，期望积分最大，所以答案选 \(B\) 。

第 2 问

小问 2 (5分)如果游戏中，小明被击落后，其之前积累的的积分会清零。那么为了游戏结束时的期望积分最大化，小明也会选择一个最优的时间主动结束游戏。请问在游戏结束时（小明主动结束游戏、或积分减到\(0\)），下列哪一个选项最接近游戏结束时小明的期望积分？

\(A. 2\)
\(B. 4\) \(C.6\) \(D. 8\)

Solution

继续写代码模拟其过程，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

public class WarPlaneGame {

    private static Random random = new Random();

    private static double getExpectScore2(int planes) {
        double score = 2;

        for (int i = 0; i < planes; i++) {
            double waitTime = -Math.log(1 - random.nextDouble());

            if (waitTime >= score / 2) {
                score /= 2;
                break;
            }

            score -= waitTime;

            double possonPr = random.nextDouble();

            double shootDownPr = Math.pow(0.85, i + 1);

            if (possonPr < shootDownPr) {
                score += 1.5;
            } else {
                score = 0;
                break;
            }
        }

        return score;
    }

    private static double simulate(int planes) {
        int simulateTimes = 50000;

        double scoreSum = 0.0;

        for (int i = 0; i < simulateTimes; i++) {
            scoreSum += getExpectScore2(planes);
        }

        return scoreSum / simulateTimes;
    }

    public static void main(String[] args) {

        for (int i = 1; i <= 8; i++) {
            double result = simulate(i);
            System.out.println("Shoot down " + i + " planes, Expect Score: " + result);
        }
    }
}

输出如下：

1
2
3
4
5
6
7
8

Shoot down 1 planes, Expect Score: 2.0230245088842116
Shoot down 2 planes, Expect Score: 1.7658505471404586
Shoot down 3 planes, Expect Score: 1.42047193787333
Shoot down 4 planes, Expect Score: 1.0837167796697962
Shoot down 5 planes, Expect Score: 0.8877915996251725
Shoot down 6 planes, Expect Score: 0.7556905685107955
Shoot down 7 planes, Expect Score: 0.7055539105268976
Shoot down 8 planes, Expect Score: 0.6896372679317954

可以看出最大期望积分是 \(2.023\)左右，所以答案选 \(A\) 。

Problem 3

对于实数 \(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\) .

1. 证明题(10分) 假设 \(\operatorname{tr}(A)=0\) . 证明存在 \(C > 0\) ，使得对任意正整数 \(n\)​ ，集合

\[A^{n} \mathbb{Z}^{2}:=\left\{A^{n} v: v \in \mathbb{Z}^{2}\right\}\]

是 \(C|\operatorname{det}(A)|^{n /2}\) -稠密的.

2. 证明题 (10分) 假设 \(A\) 的特征多项式在有理数域上不可约.证明与(1)相同的结论.

注: 这里 \(\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}\)​ 中的非零向量.)

Solution

Problem 4

设 \(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\).

1. 证明题 (10分) 证明: \(f\) 的全部特征值为 \(-d,-d+1, \cdots, d\).

2. 问答题 (5分) 记 \(W\) 是从属于特征值 \(-d+2 k(0 \leq k \leq d)\) 的 \(f\) 的特征子空间的和. 求 \(W \cap U_0\) 的维数.

3. 问答题 (5分) 对任一整数 \(j\left( 1 \leq j \leq \dfrac{d}{2}\right)\) ，求 \(W \cap U_j\)​​​​的维数.

Solution

Problem 5

证明题 (20分) 对于 \(\mathbb{R}^3\) 中的任何中心对称的凸多面体\(V\) ，证明可以找到一个椭球面 \(E\) ，把凸多面体包在内部，且 \(E\) 的表面积不超过 \(V\) 的表面积的 \(3\) 倍。

Solution

Problem 6

第 1 问

假设有一枚硬币，投掷得到正面的概率为 \(\dfrac{1}{3}\) 。独立地投掷该硬币 \(n\) 次，记 \(X_n\) 为其中得到正面的次数。试求 \(X_n\) 为偶数的概率在 \(n\) 趋于正无穷时的极限。

Solution

由题意可知，设随机变量 \(X_n\)表示在 \(n\)次独立投掷中正面出现的次数，每次出现正面的概率为 \(p = \dfrac{1}{3}\) ，则 \(X_n\) 服从参数为 \(\operatorname{B}(n, p)\) 的二项分布³ ，那么 \(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}\]

由 二项式定理⁴ 可知：

\[\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}\]

第 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 )\) 。

Solution

由多项式定理⁵ ：

\[\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}\]

Problem 7

有这么一个音乐盒，它上面有一个圆形的轨道，轨道上的一点处还有一棵开花的树。当音乐盒处于开启模式时，音乐盒会发出音乐，轨道会按照顺时针匀速转动。

你可以在轨道上放置象征恋人的两颗棋子，我们不妨称它们为小红和小绿。当小红和小绿没有到达树下时，它们就会在轨道上各自移动。当某一颗棋子到达树下时，它就会在树下原地等待一段时间。此段时间内，如果另外一颗棋子也达到了树下，那么两颗棋子就会相遇，之后在它们将随即一起顺着轨道移动，不再分开；否则，等待时间结束，两颗棋子将各自顺着轨道继续移动。

考虑这个音乐盒的数学模型。我们把这个圆形轨道参数化成一个周长为 \(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)\) 的形式。

1. 证明题 (10分) 我们设 \(f:\mathbb{R} \to \mathbb{R}\) 是一个光滑函数， 满足

\[f^{\prime} > 0, \quad f^{\prime \prime} < 0, \quad f(0)=0, \quadf(1) = 1 .\]

并设 \(\varepsilon\)是一个充分小的正的常数。我们定义等待时间函数

\[K(\phi ) = f^{-1}(f(\phi ) + \epsilon ) - \phi .\]

证明除了唯一的例外（特定的初始距离）之外，无论小红和小绿的初始距离如何，他们最终会相遇的。

2. 问答题 (10分) 我们考虑一个如下形式的 \(f\) 函数

\[f(\phi ) = \frac {1}{b} \ln \left (1 + \left (e^b - 1 \right ) \phi\right )\]

这里 \(b>0\) 是一个常数。当 \(b \ll 1, \varepsilon \ll 1\)时，请估算出相遇之前小红小绿走过的圈数的数量级。

Solution

由题设条件，我们知道当小红和小绿都没有在树下时，都随着圆形轨道顺时针匀速转动，也因此小红和小绿只可能在树下相遇。

不妨设最初条件为小红在小绿的前方，两者的距离为 \(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}\) 可知：

1. 如果 \(f^{-1}(f(d_k) + \varepsilon) \ge1\)，则小红将在等待时间结束之前到达树下，小红和小绿相遇，结束分析。

2. 如果 \(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\) 种情况：

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)\]

参考文献