|
她的名字叫艾丽亚娜（Arianna），她还在使用着40年前和她离婚的前夫的姓，Rosenbluth。在几乎所有人看来，在这家老人院所有失去自理能力的老人当中，她是普通得不能再普通的一位。 Dr. Arianna Wright Rosenbluth in 2013. She helped create what has become one of the most important algorithms of all time. Credit...via Rosenbluth family （NYT） 多年来艾丽亚娜也不认为自己有何杰出之处。所以，在76岁的时候她曾接到洛斯阿莫斯国家实验室的一位物理学家的电话，对方恭维她说您在一项叫做Metropolis Algorithm的科学算法中做出了杰出贡献啊，我们想请您做一个口述回顾。艾丽亚娜说Metropolis算法是什么啊，从来没听说过。直到对方反复提示甚至列出参考文献，她才恍然大悟，意识这个东西原来就是自己在50年代曾经参与过的一个课题：“原来你说的就是那个”（oh, that thing）。 这是他们当年发表的论文，按照姓氏排序艾丽亚娜的名字排在共同作者的第二位，和所有做科研的人一样，她以为这篇论文问世后就埋藏在故纸堆中，很快被人遗忘，就像99.99%的科技论文一样的命运。 但是，你今天如果随便谷歌一下诸如“20世纪10大科学算法的话”，这个以艾丽亚娜为第二作者的Metropolis Algorithm for MCMC经常是排名第一的。 这个算法神奇在哪里？这个MCMC又代表了什么？这篇文章有五位作者，包括艾丽亚娜在内，他们分别都做出了什么贡献呢？  如果想略微谈得深入一些的话，我们可能不得不让时光倒流300年，去认识一位叫布冯的法国数学家和博物学家（Georges-Louis Leclerc, Comte de Buffon）。 对于一些比较艰深的数学问题，除了使用理论推导，演算和证明之外，布冯认为可以通过大量的重复试验而去逼近客观真理。举一个简单的例子， 圆周率pi究竟是多少？当时数学家已经证明 pi是一个无理数，也就是说没法用两个整数的相除，或者其他一个精准的公式来概括表达它。于是布冯提出了著名的“布冯投针”试验来估算pi，如下图所示，把一把大头针随机投向画有几条平行线的纸板，针的一半长度和平行线之间距离的比例，乘以所有这些大头针中跨越平行线者的比例，就是pi的估计值。大头针的数量越多，或者投掷计数的次数越多，计算值就越逼近pi的真实值。这个术语叫做“模拟”，simulation。 再让时间快进到第二次世界大战，在美国研发原子弹的曼哈顿计划中，科学家们发现中子可以轰击原子核引发核裂变，而中子的随机运动非常适合于使用模拟计算的方法。但是这个东西的计算量实在太大了。中国两弹一星研制过程中的一个著名轶事就是，为了验证苏联专家留下的一个重要参数，邓稼先领着上百的中国科学家硬是打了半年的算盘。 曼哈顿计划中没有算盘，用布冯投大头针的方法去模拟核裂变中的中子轨迹也不现实，必须使用计算机。 于是Nicholas Metropolis这个名字就走进了历史舞台，他是一位计算物理学家，在芝加哥大学期间结识了费米和泰勒两位现代物理学的宗师（费米人称原子弹之父，泰勒是美国氢弹之父）。在他们的引荐之下，曼哈顿工程的总负责人奥本海默把Metropolis请到了新墨西哥的洛斯阿莫斯国家实验室，在这个地方，Metropolis又和两位大牛有了交集，一个是费曼，这是一位在人类所有物理学家排名中位居前10的人物，甚至超过了狄拉克和薛定谔。还有一位是冯诺伊曼，这是一个大概所有的理工男都熟悉的名字，因为迄今为止所有的计算机都可被名命为冯诺伊曼型。 在费米与冯诺伊曼的启发下，Metropolis设计出了洛斯阿莫斯的第一台用于热力学模拟的计算机，它具有冯诺伊曼型计算机的最主要的特征：具有一个内存，这相当于计算机的大脑，它依靠读入内存的计算机程序来指挥计算机电子管中的电流，进行科学计算。这样Metropolis的计算机给未来工作奠定了物质基础，但是这个大名鼎鼎的算法的问世，还需要另外两个重要条件，一个理论框架，和一个码工来把这个理论思想写成代码。 40年代曼哈顿计划刚刚正式立项的时候，泰勒就来到了洛斯阿莫斯。广岛长崎核爆之后，苏联紧跟美国脚步也成功爆炸原子弹。于是美国政府的注意力立即转移到了氢弹项目上，泰勒开始组织团队，第一个想到的当然是自己在芝大培养的研究生，但不是他最著名的学生杨振宁，而是一个相对默默无闻的学生：一个叫Marshal Rosenbluth犹太人。于是Marshal从四季如春的加州搬到了荒漠一般的洛斯阿莫斯实验室，和他一块来的还有他的新婚妻子艾丽亚娜。 在科学史上，艾丽亚娜比她的丈夫更加默默无闻，但是她的学术背景比Marshal只强不弱。从世俗一点的观点看，丈夫的研究生导师是费米和泰勒，虽然他们是美国核武器的奠基人，但都和诺贝尔奖的殊荣擦肩而过；而艾丽亚娜年仅21岁就拿到了哈佛大学物理学博士，导师是诺贝尔奖得主John Van Vleck，她是哈佛历史上第五位拿到理学博士的女性，她在斯坦福做博士后研究的时候，和未来的丈夫相识。 Dr. Rosenbluth in an undated family photo. As a young woman she was a champion fencer.Credit...via Rosenbluth family （NYT) 泰勒给学生布置的任务是：以Metropolis的计算机为依托，用数值模拟的方法做出一套基本粒子在热力学分布场中的运动模型。这当然也包括用来轰击原子核引发核裂变的中子的运动分布图。 Marshal Rosenbluth很快拿出了初步结果。他不依赖于高深的量子力学理论去构建粒子的运动分布，而是让粒子做随机的运动。正如在布冯投针实验中，大头针的一端在接触纸面后，针体以相等的几率向任何一个方向随机倒下。 不过，如果基本粒子真的可以不受限制地到处转悠的话，这个热力学分布场就了均匀的一锅粥了，和实际对不上号。Marshal最聪明的地方在下一步。 根据粒子场的能量越低就越稳定的基本物理定律，他提出，如果随机模拟把假想中的粒子带向一个能量更高的位置，那就减少了体系的稳定性，这样的随机运动会被拒绝；而如果粒子的目的地由于低能级而增强系统的稳定性，那些这个随机运动就被许可，粒子会以下一个点为起始开始新一轮的随机运动，被拒绝，或者被接受，这个过程周而复始，运行多少次取决于程序员设定的循环次数。 我对Metropolis算法的粗浅解释可以归纳为上面这简单的几句话，但是大量虚拟粒子的随机运动和每步定位，需要的是Metropolis计算机的模拟计算，和驱动这个这个运算大脑的海量计算机代码。 这完全是艾丽亚娜的任务。 当现代科技史家们在回顾20世纪中叶时尚在襁褓之中的计算科学的时候，他们或许会吃惊地发现，女性在早期的计算机编程领域曾经扮演了多么关键和主导性的角色。比如当年计算机语言的霸主，Cobol语言的创始人Grace Hopper，和NASA以计算登月轨道而出名的黑人女数学家Katherine Johnson。和这些彪炳史册的女杰们相比，艾丽亚娜只是一个更加不出名的幕后英雄罢了。 也许有科技史家们推测，在当今这个以男性为主导的计算机编程领域，女性其实早在最初阶段就先走了一步，但是很多当年计算领域的女杰被长期埋没了，所以现在的编程领域是鸠占鹊巢，男人取代了女人，这个历史现象很可能是性别歧视作祟的结果。 这个话只对了一半。 事实是，当年的算法编程被认为是类似于文书文秘一类的简单重复性工作，反而被认为是最适合女性的。 Metropolis算法论文署名最后的，也是最有名的作者，是氢弹之父泰勒，倒数第二位叫Augusta Teller, 她是泰勒的妻子，昵称叫Mici，也是一位女科学家。他们夫妇是最早入驻洛斯阿莫斯国家实验室的核物理学家，在那个白手起家的时代，科学家的家庭就类似一个创业的手工作坊，男主外女主内，男人在科研攻关中厮杀，沉浸在理论研究和公式推导中，而这个“内”则是各种行政琐事，其中居然还包括了被认为是“琐屑”的计算编程工作。 这是因为，在曼哈顿计划的前期，大量的计算工作是靠女人，很多是科学家的太太们，她们操作手摇计算器完成的。带有内存的冯诺伊曼型计算机问世之后，才慢慢有了编程这个工种，自然而然地也被太太们捷足先登了。泰勒夫人就扮演了这样的角色，她的编程能力曾在原子弹计算中展露头角。 Marshal Rosenbluth夫妇来了。作为老板的泰勒给了一个理论的大框架大方向，全部的细节论证和推导都是Marshal的任务，同时泰勒夫人也顺理成章地把编程的工作交接给了学生的妻子艾丽亚娜，这个哈佛大学物理学博士。 是她写出了运行这个算法的全部程序。 遗憾的是，我们今天无从得知艾丽亚娜这个工作的任何细节。2003年，物理界有一个庆祝Metropolis算法诞生50周年的研讨会，在这篇文章的五位作者中，Metropolis和泰勒夫人已经作古，95岁高龄的泰勒中风后失去了思维能力，只有76岁身患癌症的Marshal Rosenbluth拖着病体坚持来到大会，给这个算法的问世留下了一段珍贵的口述历史，从大会回来之后，他很快也驾鹤西行了。 令人不解的是，身体最好的艾丽亚娜却没有得到邀请，只是由大会召集人给她打了一个电话，当时距Rosenbluth这对科学伉俪的劳燕分飞，也已经25年了。 我们都知道一个好的程序员需要有比较强的逻辑思维能力，现代的编程语言基本还以英语为基础的，比如if then, do until, for (i in 1:100)。而当年Metropolis研发的计算机内存只懂0和1这样最简单基本的信号。也就是说，艾丽亚娜必须先要把丈夫的算法理论吃透了，然后把他们转换成计算机内存指挥电子管开关的逻辑流程，然后再转换成无数行0和1这样只有计算机才懂的天书。 仅凭这点粗浅的推理，我感觉她的这个工作比今天湾区FLAG动辄年薪半个米的高级码工要难多了。 在规定了假想粒子在模拟场中的运动规则之后，Marshal Rosenbluth进而推导出这篇论文中最重要的理论结果：只要粒子的每步的运动轨迹符合能级最低规律，在经过大量的模拟运动之后，整体的粒子运动规律将无限向统计力学中经典的麦克斯韦-玻尔兹曼分布收敛。 如果打个简单比方的话，这就好像是在布冯的投针实验中，随着投针数量的无限增加，pi的模拟值将无限趋近于3.1415926和3.1415927之间。 也许我和读者的数学能力都不足以把这个证明和推导讲清楚和理解明白。那么就让我们把艾丽亚娜的程序可视化，产生类似下面的一条基本粒子的运动轨迹，这在数学上有个学名叫马可夫链（Markov Chain)，这也就是之前说的MCMC算法当中的第一个MC，它说的是粒子在每个位点之间的移动符合固定的概率。但是运动轨迹的大趋势却是斗折蛇行，好像毫无规律，不是吗？ 一提“钟型曲线”，人们也许意识到这可以是一个概率分布。Metropolis算法的初始目的是用计算机模拟粒子运动的玻尔兹曼分布，但是这个技术可以推广到任意的统计分布中。也就是说，在实践中，人们可以通过特定的规则构建一条随机运动的链条，这个链条终点的位置就代表了从所需的任意概率分布中的随机取样。 这就决定了Metropolis算法的普世价值。 在完成Metropolis算法之后，艾丽亚娜和丈夫双双离开了洛斯阿莫斯，当时她还不到30岁，却从此放弃了职业生涯，直到93岁高龄去世，大概当了60年的家庭妇女。 艾丽亚娜的科学成就仿佛是彗星般地灵光闪现，然后就沉寂了，直现在人们才知道她是一位相当特立独行的女性，她在上高中期间就是一位几乎达到职业水平击剑选手，有时参加男子比赛的。可惜壮志未酬，第二次世界大战让她失去了首次参加奥运会的机会，1948年的伦敦奥运会她没钱参加。聊以自慰的是，她在第二年拿在哈佛到了物理学博士。 当地报纸报道小女孩艾丽亚娜的击剑成就 也没有人知道她早早退隐江湖的原因，是对科学厌倦了吗？应该不是，因为她女儿小的时候曾记得母亲闲着没事推导数学公式，权当是大脑散步。也许是她认为业余搞数学可以平衡家庭和事业的双重压力？因为她儿子回忆到母亲在持家之余，曾钻研究纽结理论，这是高端的拓扑数学，不过一直没有发表什么成果。 1978年，50岁的Rosenbluth夫妇结束了维系了20多年的婚姻，两年后Marshal娶了加州一位著名的艺术家，艾丽亚娜后来一直没有再婚，直到在93岁的高龄上因新冠去世。 不过Metropolis算法的命运，和自己的作者相比却仿佛是截然相反。 在对这个事件的历史回顾中，我们或许可以看到，计算物理学家Metropolis主要贡献其实就是提供了计算机，给这个课题奠定了物质基础，成为发表文章的首席作者，也拿到了这个经典成就的冠名权，其实他并未做出真正具体的贡献。 但是Metropolis也并非浪得虚名，他和另一个曼哈顿计划的物理学家Stanislaw Ulam合作，在1949年写了一篇经典文献《蒙特卡罗方法》Monte Carlo Method，MC，此文的核心是，如果某个特定概率无法用数学或物理方法推导，比如一个形状奇怪重心未知的筛子，就可以使用反复投筛子的方法估算其每面朝上的概率。Ulam把这个经验方法命名为蒙特卡洛，这是法国南部一个小国摩纳哥的城市，以博彩业出名。据说Ulam小时候总是看到叔叔去蒙特卡洛赌钱，就用这个掷色子的城市命名了这个掷色子的方法。这就是MCMC算法中第二个MC的含义。 当然实际的问题要比掷色子复杂太多了，计算量也大太多了，当时只有Metropolis的计算机能够胜任。同时这个方法也蕴含了从特定的概率分布中随机取样的思想，所以这篇经典文献发表在《美国统计学会杂志》JASA，这预示着这一系列不平凡的物理学思想将在统计学领域发扬光大。 也许有人要打破沙锅问到底，从概率分布中随机取样有何用处？ 就举一个最热门的例子吧，大家都知道辉瑞RNA疫苗在临床实验中的有效率高达90%以上，那么大家在激动的同时也关心，这个被辉瑞力捧上天的伟大有效率，它的误差是多少？ 疫苗有效率是由两个数字来决定的：试验疫苗组中的新冠感染数，和对照组的感染数。前者越小，疫苗的效力就越高。那么计算疫苗保护率的误差就有一个很简单直观的办法：疫苗组和对照组感染数是两个随机变量，如果我们更够找到它们的统计分布的话，就可以从这两个分布中大量随机取样，让两者的差值除以对照组的随机感染数，得到的就是疫苗保护率的估计值，然后对得到的这一大堆估计值进行掐头去尾，就得到了所谓了95%的置信区间，这是一种在临床试验中最常见的统计量。 辉瑞疫苗有效率的95%的置信区间 再回到古老的布冯投针试验，那也是一种随机取样，只不过是占了针末端向四周各个方向以均等概率倒下的便宜，所以这个取样不需要计算机，只需要牛顿万有引力自个干活就成了。 但并不是任何的随机取样都是如此容易得到的。举个最简单的例子，怎样获取一系列成年人类身高的随机数呢？ 如果采用最简单的均匀取样的方法，类似布冯投针那样，你也许能得到一列这样的“身高”样本（以厘米为单位）： 140，146，152，158，165，169，172，179，185，189，194，197， 202， 209， 215 明眼人都能看出，这列数固然“随机”，但并不能代表真实世界中人类的身高，因为人的身高属于正态分布，中等身材的人多，特高特矮的人少，象上面那样的2米以上的人和170左右的人一样多，不可能。真实的人类身高随机取样应该是象下面这样的，中间多，两边少。 157， 162， 165， 168， 171， 172， 173，175， 177， 179， 184， 192 怎样才能得到符合正态分布的随机数？在没有计算机的年代，人们会去下苦功夫在实际世界中进行大量的测量。比如在本文的上集，《借茅台院士的热度，科普这样一位啤酒总工》一文中，我记录了三位有师承关系的统计学大师的事迹，他们是皮尔森，戈塞特和费舍尔。其中戈塞特是著名的Student’s t test的发明人，皮尔森是他的老师，其拿手好戏就是派学生在民间大量采样，比如监狱犯人或者军人的身高体重胸围等生理参数，然后依此绘制出完备的统计曲线和图表。
现在有了计算机，产生随机数就太简单了。因为所有统计学软件都储存了完备的概率分布数学解析式，一切照公式而行即可。 但关键是，你要知道公式才行。 在所有的学科分支中最流行的正态分布，它有一个简洁而优美的数学公式：
这个公式是高斯推导出来的初始样子，和现代课本中的略有不同。为了纪念先贤，特保留原样，所以正态分布也叫高斯分布 即便是这个最常用的公式，也是历经了从伽利略，棣莫弗， 到拉普拉斯，高斯，几代数学家历时300多年的摸索才拿到的。 给一个统计变量的分布找到解析解非常困难，我们可以再举那个新冠疫苗有效率的例子。 我在前面提到疫苗临床试验中，疫苗组和对照组的感染数都符合某个特定的概率分布，当时卖了个关子，没提究竟是什么样的分布，其实简单的很，就是常用程度仅次于正态分布的二元分布， binomial distribution。因为它的结果是不连续的，不是0就是1，感染或者不感染。疫苗有效率是1减去两个感染率的比值，这是一个简单的算术操作。 可让人不可思议的，疫苗有效率，这个对两个最简单统计变量进行的最简单的算术操作，任凭你是牛顿还是高斯的数学天才，硬是给它找不到一个解析解，也就是一个象正态曲线那样的公式，术语叫closed-form。没有解析解就无法直接按公式进行直接取样，而只能采用数值模拟的办法（相比之下，两个正态分布变量之比例，就有一个解析解，它在物理学上叫洛伦斯分布，在数学上叫柯西分布）。 另一个概率分布公式之难求的例子，就是我们上文提到的“啤酒”总工，威廉戈塞特，爱尔兰吉尼斯啤酒集团的总工，也是Student’s t 检验的发明者。他在啤酒酿造过程中，发现小样本实验的均值和标准方差之间的比值是一个有规律的统计变量分布，这被后世名命为Student’s t 分布，取自当年戈塞特发表论文的笔名：Student。 戈塞特是数理统计鼻祖皮尔森的学生，也是牛津数学系的高材生，但就是他也没有能力给出t分布的解析公式。是戈塞特的晚辈费舍尔（R.A.Fisher）这位数学天才，在几年后解决了这个问题，也许是直接的代数推导不易，他采用的是高维几何的方法。对此戈塞特不明觉厉，但直觉上知道晚辈是对的。费舍尔得到的t分布的解析解是长这个样子的：
明显比她的母分布正态曲线要复杂太多了。 作为实践大师的戈塞特，在费舍尔的答案出来之前，他的解决方案是不厌其烦地从大量实验中测出各种型态t分布的概率和相应的关键值，以便啤酒车间的实验人员按图索骥依数查表，做出合乎科学的决定：这罐啤酒母液要不要倒掉？ 从某种意义上，戈塞特和他的工作人员们类似于布冯投针实验中人肉大头针。在大规模的科学计算中，这是一种不可持续的玩法。 由此可见，做一个“蒙特卡洛”实验，从一个缺乏解析解的统计分布中随机取样，是不容易的。如果我们还记得的话，Marshal和艾丽亚娜夫妇的工作，恰恰通过马可夫链的数值方式，用计算机从未知分布中随机取样，这就预示了他们这个工作在未来的意义。 我们也会看到，在解析解未知的t分布中取样，或者是在根本就没有解析解的疫苗有效率的分布中取样，他们的难度和另外一类的问题相比，根本就不在一个层次上。  就在法国人布冯提出投针实验这个天才构想的15年前，有人在英国皇家学会年会上宣读了一篇文章，有些不同寻常的是，这篇文章的作者已经在两年前仙去了，他的名字叫托马斯贝叶斯（ Thomas Bayes）。在今天的统计学，工程和医学制药中，这个名字已经是如雷贯耳了，因为他被认为是贝叶斯统计学派的开山鼻祖。
几个星期之前，网上有一篇奇文在流行，号称是用贝叶斯模型证明了新冠病毒来源是人工合成并泄露的概率高达98%。此文大家读来都不明觉厉，因为生物学家觉得作者是一个统计高手，而统计学家觉得他们病毒专家。 后来读者请贝叶斯算法的权威Andrew Gelman评论，他说此文的生物学他不懂，但是其贝叶斯部分是胡扯。作者把贝叶斯分析掺和进来，分明是用来拉大旗作虎皮吓唬人的。 难道使用了贝叶斯就是高大上了吗？历史上并非如此，实际上在很长一个历史时期内，贝叶斯概率思想被认为是离经叛道的。因为，就在250年前在英国皇家学会上宣布的那篇论文，托马斯贝叶斯第一次提出了“反向概率”（reverse probability）的思想。 再用新冠疫苗当例子，有一个我们所感兴趣的参数是：接种疫苗后依然感染的概率，P（感染 | 接种疫苗）。接种疫苗会诱发人体的免疫性，从而降低感染率和死亡率，这两者有直接的因果次序，所以这个概率是符合常理的。但是贝叶斯在摆弄各种条件下的概率换算时，却遇到了这么一个问题：在所有的已感染病例中，有多少是在之前接种了疫苗的？P（接种了疫苗 | 感染）。 在贝叶斯死后才见天日的这篇经典文献中，他给出了一个被后世尊称为贝叶斯定理的概率换算，如果换成新冠疫苗的例子就是这样的： P（感染 | 之前接种了疫苗）= P（之前接种了疫苗 | 感染）X 群体感染率 / 群体接种率 显然，接种疫苗会降低感染率，但是今天感染不会影响昨天的接种行为，所以在等号右边的第一个概率，在表达次序上有本末倒置之嫌，有人叫它“反向概率”。这个概念在认知哲学上的真正含义，直到今天还有争议，就遑论200多年前的人了。 同时，这也是一个极容易引发误导的概念。举个例子：即使在疫苗普及后，依然也会有极小部分人拒绝接种。这个时候，你肯定希望使用公式左边的那个概率去说服他们，因为接种疫苗后的新冠感染率和死亡率都会变得很低很低。但是它的“反向概率”却可能把你吓一跳，因为在确诊病例中，会有很大比例也都是接种过疫苗的。 这个能作为反对疫苗的理由吗，绝对不能。具体原因自己去想，因为这不在本文主题之内。 在历史上，贝叶斯本人只是给后世的贝叶斯学派开了个头，这个领域内第一位真正的大师是法国的拉普拉斯，Laplace，他的名字在前文出现过一次，他在高斯之前为推导正态分布的解析解做出过开拓性贡献。但是贝叶斯体系在当时的争议太大了，在拉普拉斯去世之后，他的朋友建议在讣告悼词中免谈死者在贝叶斯领域的工作，因为“何必给逝者脸上抹黑呢”？ 在20世纪的大部分时间，贝叶斯技术虽然偶露峥嵘，比如帮忙破解了纳粹德国的电报密码，但是其应用基本处于冬眠状态，只有少数几个理论家在默默耕耘。这主要的原因是著名的费舍尔，R.A.Fisher，他不是贝叶斯的粉丝，无法接受这个违反常识的所谓“反向概率”，再加上此人是数学天才，我们也许记得他超越了前辈戈塞特而独立推导出了t分布的解析解。后来费舍尔几乎用一己之力奠定了以常规概率为蓝图的统计方法之数学基础，他的学派史称概率学派（Frequentist），把贝叶斯学派整整压制了半个世纪。 你很快就会看到，是Marshal和艾丽亚娜夫妇的Metropolis算法让贝叶斯门派满血复活。 贝叶斯定理即便是在今天也非常有实用价值的。假如你想用接种疫苗后的极低感染率来给辉瑞Moderna做宣传，你会惊讶地发现这个精确的统计数字是找不到的，因为这要求全国每一个人的接种情况和感染状态都得到精准掌控，至少美国政府没有能力做到这一点。 最简单的办法是用贝叶斯定律换算。 由于感染人口远远小于全体人口，所以相对容易统计确诊者之前的疫苗接种情况（在核酸检测表中加这么一个问题就够了）。至于贝叶斯定律中的总确诊率和总接种率呢，这个大概在每天的新闻中都能找到。这样接种后的感染危险就算出来了（希望是大大地降低了）。 也许更重要的是，贝叶斯定理中蕴含了统计学的核心思想。 让我们再复习一下新冠中的贝叶斯定律： P（感染 | 接种过疫苗）= P（接种过疫苗 | 感染）X 群体感染率 / 群体接种率 在疫情控制中，我们最关心的是感染率的概率分布，因为我们需要这个参数来评估未来医院的承受力，和经济重新开放的政策，显然这个概率是和疫苗接种情况相关的，所以公式的左边被叫做后验分布（Posterior Distribution），因为它是在获知疫苗接种情况之“后’才获得的有条件的概率分布。 在研究一个函数的变化趋势的时候，最直接的方法是对其求导（微积分的范畴）。而在贝叶斯公式的分母中，群体接种率中不包含感染率这个参数，因此它的导数就是一个常数。所以贝叶斯在这里，把一个概率分布的核心写成了两个分布式的乘积。其中，整体感染率是在获知疫苗接种情况之前的分布，被称之为“先验分布”（Prior Distribution），在某些场合下也可被称为人类在没有客观数据时对未知事物的主观判断。另一个概率叫做似然函数。 如此以来，人类就获得了几乎无限之多的获取概率分布的方法。 但是老问题依然没有解决，我们或许还记得，两个二元分布变量的简单除法没有解析解；从正态分布变量的简单算术中拿到的t分布，聪明如戈塞特的都没有拿到解析解，最后是靠费舍尔这个数学天才用高维解析几何的方法解决的。 贝叶斯公式也面临类似的困境，把一个未知概率分解为两个已知分布的乘积，并不意味着容易拿到闭合形式的数学解析解，而没有公式就无法从中取样，对未来的疫情做推断和预测。更何况，在真实的贝叶斯模型中参数可以是层层嵌套的，因此后验分布的表达写成几十个概率分布的乘积也有可能，那样的话就是一百个费舍尔来也是无能为力了。 这个难题在1953年看到了第一束希望的曙光，因为Marshal和艾丽亚娜夫妇来了。 我们也许还记得，在Metropolis算法中， Marshal Rosenbluth给粒子模拟运动设立的规则是，允许粒子移到能级较低的位点，但如果目的地的能级变高就说明这个运动给体系带来不稳定因素，那么就得打回原点重选方向。 Marshal也许没有想到，这个规则简直就是给贝叶斯的后验概率度身定做的。艾丽亚娜写的计算机程序既可以模拟中子的运动，也可以构建数值取样的马可夫链，其尾端的移动方向，可以取决于各个候选值之间后验分布的比较，这和玻尔兹曼场中能级来决定运动是一个道理。 读者也许要问了，我们不是还不知道这个后验分布的真面目吗？的确如此，但是，未知后验分布和两个或数个已知概率的乘积成正比， 所以马可夫链的运动方向可以通过贝叶斯定律的换算而确定。 正如Marshal Rosenbluth证明了，遵循能量最低原则的粒子计算机随机模拟运动，最终无限向麦克斯韦玻尔兹曼分布收敛；而在先验概率和似然函数乘积指引下马可夫链，在贝叶斯定律这只无形之手的牵引之下，也会无限地向真正的后验分布收敛，这个蒙着面纱的美女就真相大白了。 十几年后统计学科班出身的人又对这个算法做了必要的补充和发展，所以这个算法现在叫Metropolis-Hasting算法，简称M-H算法。 进入了90年代后，廉价高性能的个人电脑开始走入千家万户。以MCMC方法为基础的贝叶斯学派才第一次走出了象牙塔，它不再是数学家们的专利，人们不必在无穷无尽的求导和积分中穷经皓首，而是可以借助强大的计算机，靠马可夫链的随机取样，去解决实际中的问题。 贝叶斯学派终于咸鱼翻身了。 在接近本文结尾的时候，我们再来欣赏一下M-H方法暴力的计算美学。今天的贝叶斯玩家们，也许大部分都读不懂Marshal Rosenbluth证明和推导了，但是人人都站在巨人的肩膀上，也就是艾丽亚娜当年用无数的0和1堆砌起来的大厦。 这个例子很简单，让我们观察到的数据的似然函数是一个正态分布，它的期望值和方差的先验分布分别又是正态分布和均态分布，我们的任务是从这些数据和假设中推导出两个正态参数的后验分布。如前所述，即使是这样的简单例子也没有解析解，唯有使用MCMC的数值模拟办法。 左图显示的是三条独立的马可夫链，可以看出他们的起始位点颇为不同，但是在M-H算法的拉动下，三条链很快就收敛了。右图是在马可夫链上对均值和方差的随机取样，绘制成三维的概率密度图，引人瞩目的是，在模拟的初始阶段，峰型非常不稳而且四处位移，但是随着马可夫链走到了五千步以上，这个二维的分布开始收敛成一个稳定漂亮的峰型。
https://blog./2013/09/an-animated-peek-into-the-workings-of-bayesian-statistics.html 事就这样成了。 在MCMC计算中，维度更高的复杂概率问题也是以此类推，我们也许永远不能给他们写出一个优美的数学表达式，但是却可以通过随机采样的方式给它画一个逼真的肖像。 这就是概率和模拟的优美之处，我们今天要讲的这个科学史上的传奇，就讲完了。 当然，当时Metropolis算法文章的五位共同作者，都没有意识到这个技术所蕴含的普世意义，所以他们日后都没有在研究中再使用过这个方法。 Metropolis后来继续从事计算科学在物理中的应用， 他获得了该算法的冠名权，这一荣誉足以让他名垂青史。 作为学界权威的泰勒，在氢弹成功后，大概不愿意继续留在洛斯阿莫斯试验室了，这毕竟是奥本海默的势力范围。泰勒夫妇搬到了加州伯克利，创立了劳伦斯利沃莫国家实验室，和洛斯阿莫斯分庭抗礼。到了麦卡锡时代，泰勒去国会作证，还把老战友奥恩海默给卖了，这是题外话。 Marshal Rosenbluth后来也搬到了加州，专业方向变成了凝聚态物理，也和这个他们点燃了第一簇篝火的领域拜拜了。直到2003年，在学术界聚集庆祝Metropolis算法问世50周年的大会上，他强撑病体前往，做了一个对历史的回顾，厘清了理论工作贡献的来龙去脉。人们开始意识到，这个算法的真名或许应该叫Rosenbluth-Teller算法。 艾丽亚娜的贡献是最不为人知的。 几年前，有人为了纪念这个划时代算法，给这篇文章的五位作者做了一个“全家福”，作为冠名者的Metropolis自然是画面的中心，以艺术气质见长的泰勒正在弹钢琴，他的太太坐在一旁听得入迷，年轻的Marsha Rosenbluth也是一副春风得意的气度，唯有艾丽亚娜的影像是一个黑影，因为当时连网上都找不到她的影像信息。
除了Metropolis算法，统计界还特别推崇那篇1949年Metropolis和Stanislaw Ulam发表的《蒙特卡洛方法》。由于Metropolis已经备极哀荣了，贝叶斯界把他们最新创立的一门专用于马可夫链随机取样的计算机语言名命为Stan， 为了纪念当年的共同作者之一Stanislaw Ulam，但是Stanislaw写过程序吗？无人知晓。 回首当年，Stan语言的创立者之一Andrew Gelman不无遗憾地说：也许我们应该把它命名为艾丽亚娜。 |
|
|
