通知

備用：如果以后讀者無法再閱讀到本公號的更新，可以去網易搜索同名的網易號。海外的讀者也可以下載Telegram , 然后再手機瀏覽器中打開鏈接“https://t.me/joinchat/MvXTABj7X6uQxRjnNxaHmg”，點擊“ Join Group”加入“北美新藥科普歷史網”的讀者群。

2020年底，距離新年鐘聲的敲響還有三天的時候，洛杉磯郊區(qū)帕薩迪納小鎮(zhèn)的一家老人院內，一位93歲高齡的老婦人因為新冠肺炎靜靜地走完了一生。

她的名字叫艾麗亞娜（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）。

這是他們當年發(fā)表的論文，按照姓氏排序艾麗亞娜的名字排在共同作者的第二位，和所有做科研的人一樣，她以為這篇論文問世后就埋藏在故紙堆中，很快被人遺忘，就像99.99%的科技論文一樣的命運。

但是，你今天如果隨便谷歌一下諸如“20世紀10大科學算法的話”，這個以艾麗亞娜為第二作者的Metropolis Algorithm for MCMC經常是排名第一的。

這個算法神奇在哪里？這個MCMC又代表了什么？這篇文章有五位作者，包括艾麗亞娜在內，他們分別都做出了什么貢獻呢？

如果想略微談得深入一些的話，我們可能不得不讓時光倒流300年，去認識一位叫布馮的法國數(shù)學家和博物學家（Georges-Louis Leclerc, Comte de Buffon）。

對于一些比較艱深的數(shù)學問題，除了使用理論推導，演算和證明之外，布馮認為可以通過大量的重復試驗而去逼近客觀真理。舉一個簡單的例子， 圓周率pi究竟是多少？當時數(shù)學家已經證明 pi是一個無理數(shù)，也就是說沒法用兩個整數(shù)的相除，或者其他一個精準的公式來概括表達它。于是布馮提出了著名的“布馮投針”試驗來估算pi，如下圖所示，把一把大頭針隨機投向畫有幾條平行線的紙板，針的一半長度和平行線之間距離的比例，乘以所有這些大頭針中跨越平行線者的比例，就是pi的估計值。大頭針的數(shù)量越多，或者投擲計數(shù)的次數(shù)越多，計算值就越逼近pi的真實值。這個術語叫做“模擬”，simulation。

再讓時間快進到第二次世界大戰(zhàn)，在美國研發(fā)原子彈的曼哈頓計劃中，科學家們發(fā)現(xiàn)中子可以轟擊原子核引發(fā)核裂變，而中子的隨機運動非常適合于使用模擬計算的方法。但是這個東西的計算量實在太大了。中國兩彈一星研制過程中的一個著名軼事就是，為了驗證蘇聯(lián)專家留下的一個重要參數(shù)，鄧稼先領著上百的中國科學家硬是打了半年的算盤。

曼哈頓計劃中沒有算盤，用布馮投大頭針的方法去模擬核裂變中的中子軌跡也不現(xiàn)實，必須使用計算機。

于是Nicholas Metropolis這個名字就走進了歷史舞臺，他是一位計算物理學家，在芝加哥大學期間結識了費米和泰勒兩位現(xiàn)代物理學的宗師（費米人稱原子彈之父，泰勒是美國氫彈之父）。在他們的引薦之下，曼哈頓工程的總負責人奧本海默把Metropolis請到了新墨西哥的洛斯阿莫斯國家實驗室，在這個地方，Metropolis又和兩位大牛有了交集，一個是費曼，這是一位在人類所有物理學家排名中位居前10的人物，甚至超過了狄拉克和薛定諤。還有一位是馮諾伊曼，這是一個大概所有的理工男都熟悉的名字，因為迄今為止所有的計算機都可被名命為馮諾伊曼型。

在費米與馮諾伊曼的啟發(fā)下，Metropolis設計出了洛斯阿莫斯的第一臺用于熱力學模擬的計算機，它具有馮諾伊曼型計算機的最主要的特征：具有一個內存，這相當于計算機的大腦，它依靠讀入內存的計算機程序來指揮計算機電子管中的電流，進行科學計算。這樣Metropolis的計算機給未來工作奠定了物質基礎，但是這個大名鼎鼎的算法的問世，還需要另外兩個重要條件，一個理論框架，和一個碼工來把這個理論思想寫成代碼。

40年代曼哈頓計劃剛剛正式立項的時候，泰勒就來到了洛斯阿莫斯。廣島長崎核爆之后，蘇聯(lián)緊跟美國腳步也成功爆炸原子彈。于是美國政府的注意力立即轉移到了氫彈項目上，泰勒開始組織團隊，第一個想到的當然是自己在芝大培養(yǎng)的研究生，但不是他最著名的學生楊振寧，而是一個相對默默無聞的學生，一個叫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的計算機為依托，用數(shù)值模擬的方法做出一套基本粒子在熱力學分布場中的運動模型。這當然也包括用來轟擊原子核引發(fā)核裂變的中子的運動分布圖。

Marshal Rosenbluth很快拿出了初步結果。他不依賴于高深的量子力學理論去構建粒子的運動分布，而是讓粒子做隨機的運動。正如在布馮投針實驗中，大頭針的一端在接觸紙面后，針體以相等的幾率向任何一個方向隨機倒下。

不過，如果基本粒子真的可以不受限制地到處轉悠的話，這個熱力學分布場就了均勻的一鍋粥了，和實際對不上號。Marshal最聰明的地方在下一步。

根據粒子場的能量越低就越穩(wěn)定的基本物理定律，他提出，如果隨機模擬把假想中的粒子帶向一個能量更高的位置，那就減少了體系的穩(wěn)定性，這樣的隨機運動會被拒絕；而如果粒子的目的地由于低能級而增強系統(tǒng)的穩(wěn)定性，那些這個隨機運動就被許可，粒子會以下一個點為起始開始新一輪的隨機運動，被拒絕，或者被接受，這個過程周而復始，運行多少次取決于程序員設定的循環(huán)次數(shù)。

我對Metropolis算法的粗淺解釋可以歸納為上面這簡單的幾句話，但是大量虛擬粒子的隨機運動和每步定位，需要的是Metropolis計算機的模擬計算，和驅動這個這個運算大腦的海量計算機代碼。

這完全是艾麗亞娜的任務。

當現(xiàn)代科技史家們在回顧20世紀中葉時尚在襁褓之中的計算科學的時候，他們或許會吃驚地發(fā)現(xiàn)，女性在早期的計算機編程領域曾經扮演了多么關鍵和主導性的角色。比如當年計算機語言的霸主，Cobol語言的創(chuàng)始人Grace Hopper，和NASA以計算登月軌道而出名的黑人女數(shù)學家Katherine Johnson。和這些彪炳史冊的女杰們相比，艾麗亞娜只是一個更加不出名的幕后英雄罷了。

也許有科技史家們推測，在當今這個以男性為主導的計算機編程領域，女性其實早在最初階段就先走了一步，但是很多當年計算領域的女杰被長期埋沒了，所以現(xiàn)在的編程領域是鳩占鵲巢，男人取代了女人，這個歷史現(xiàn)象很可能是性別歧視作祟的結果。

這個話只對了一半。

事實是，當年的算法編程被認為是類似于文書文秘一類的簡單重復性工作，反而被認為是最適合女性的。

Metropolis算法論文署名最后的，也是最有名的作者，是氫彈之父泰勒，倒數(shù)第二位叫Augusta Teller, 她是泰勒的妻子，昵稱叫Mici，也是一位女科學家。他們夫婦是最早入駐洛斯阿莫斯國家實驗室的核物理學家，在那個白手起家的時代，科學家的家庭就類似一個創(chuàng)業(yè)的手工作坊，男主外女主內，男人在科研攻關中撕殺，沉浸在理論研究和公式推導中，而這個“內”則是各種行政瑣事，其中居然還包括了被認為是“瑣屑”的計算編程工作。

這是因為，在曼哈頓計劃的前期，大量的計算工作是靠女人，很多是科學家的太太們，她們操作手搖計算器完成的。帶有內存的馮諾伊曼型計算機問世之后，才慢慢有了編程這個工種，自然而然地也被太太們捷足先登了。泰勒夫人就扮演了這樣的角色，她的編程能力曾在原子彈計算中展露頭角。

Marshal Rosenbluth夫婦來了。作為老板的泰勒給了一個理論的大框架大方向，全部的細節(jié)論證和推導都是Marshal的任務，同時泰勒夫人也順理成章地把編程的工作交接給了學生的妻子艾麗亞娜，這個哈佛大學物理學博士。

是她寫出了運行這個算法的全部程序。

遺憾的是，我們今天無從得知艾麗亞娜這個工作的任何細節(jié)。2003年，物理界有一個慶祝Metropolis算法誕生50周年的研討會，在這篇文章的五位作者中，Metropolis和泰勒夫人已經作古，95歲高齡的泰勒中風后失去了思維能力，只有76歲身患癌癥的Marshal Rosenbluth拖著病體堅持來到大會，給這個算法的問世留下了一段珍貴的口述歷史，從大會回來之后，他很快也駕鶴西行了。

令人不解的是，身體最好的艾麗亞娜卻沒有得到邀請，只是由大會召集人給她打了一個電話，當時距Rosenbluth這對科學伉儷的勞燕分飛，也已經25年了。

我們都知道一個好的程序員需要有比較強的邏輯思維能力，現(xiàn)代的編程語言基本還以英語為基礎的，比如if then, do until, for (i in 1:100)。而當年Metropolis研發(fā)的計算機內存只懂0和1這樣最簡單基本的信號。也就是說，艾麗亞娜必須先要把丈夫的算法理論吃透了，然后把他們轉換成計算機內存指揮電子管開關的邏輯流程，然后再轉換成無數(shù)行0和1這樣只有計算機才懂的天書。

僅憑這點粗淺的推理，我感覺她的這個工作比今天灣區(qū)FLAG動輒年薪半個米的高級碼工要難多了。

在規(guī)定了假想粒子在模擬場中的運動規(guī)則之后，Marshal Rosenbluth進而推導出這篇論文中最重要的理論結果：只要粒子的每步的運動軌跡符合能級最低規(guī)律，在經過大量的模擬運動之后，整體的粒子運動規(guī)律將無限向統(tǒng)計力學中經典的麥克斯韋-玻爾茲曼分布收斂。

如果打個簡單比方的話，這就好像是在布馮的投針實驗中，隨著投針數(shù)量的無限增加，pi的模擬值將無限趨近于3.1415926和3.1415927之間。

也許我和讀者的數(shù)學能力都不足以把這個證明和推導講清楚和理解明白。那么就讓我們把艾麗亞娜的程序可視化，產生類似下面的一條基本粒子的運動軌跡，這在數(shù)學上有個學名叫馬可夫鏈（Markov Chain)，這也就是之前說的MCMC算法當中的第一個MC，它說的是粒子在每個位點之間的移動符合固定的概率。但是運動軌跡的大趨勢卻是斗折蛇行，好像毫無規(guī)律，不是嗎？

那么就讓我們在下圖中的左側加快粒子位移的速率，然后在右邊統(tǒng)計粒子在每個位點出現(xiàn)的總頻率并作柱狀圖。 這樣一個驚人的規(guī)律就開始浮現(xiàn)了，在虛擬粒子剛剛開始位移的時候，這個柱狀圖的形狀是非常不穩(wěn)定的，各個位點的頻率此起彼落一會一變。 但是隨著步驟的增多，比如在進入400步之后，各個點的頻率分布柱狀圖就趨向穩(wěn)定了，形成了一個以位點D為峰頂（最高頻率）的中間高，兩邊低的“鐘型曲線”。

一提“鐘型曲線”，人們也許意識到這可以是一個概率分布。Metropolis算法的初始目的是用計算機模擬粒子運動的玻爾茲曼分布，但是這個技術可以推廣到任意的統(tǒng)計分布中。也就是說，在實踐中，人們可以通過特定的規(guī)則構建一條隨機運動的鏈條，這個鏈條終點的位置就代表了從所需的任意概率分布中的隨機取樣。

這就決定了Metropolis算法的普世價值。

在完成Metropolis算法之后，艾麗亞娜和丈夫雙雙離開了洛斯阿莫斯，當時她還不到30歲，卻從此放棄了職業(yè)生涯，直到93歲高齡去世，大概當了60年的家庭婦女。

艾麗亞娜的科學成就仿佛是彗星般地靈光閃現(xiàn)，然后就沉寂了，直現(xiàn)在人們才知道她是一位相當特立獨行的女性，她在上高中期間就是一位幾乎達到職業(yè)水平擊劍選手，有時參加男子比賽的。可惜壯志未酬，第二次世界大戰(zhàn)讓她失去了首次參加奧運會的機會，1948年的倫敦奧運會她沒錢參加。聊以自慰的是，她在第二年拿在哈佛到了物理學博士。

當?shù)貓蠹垐蟮佬∨悂喣鹊膿魟Τ删?/p>

也沒有人知道她早早退隱江湖的原因，是對科學厭倦了嗎？應該不是，因為她女兒小的時候曾記得母親閑著沒事推導數(shù)學公式，權當是大腦散步。也許是她認為業(yè)余搞數(shù)學可以平衡家庭和事業(yè)的雙重壓力？因為她兒子回憶到母親在持家之余，曾鉆研究紐結理論，這是高端的拓撲數(shù)學，不過一直沒有發(fā)表什么成果。

1978年，50歲的Rosenbluth夫婦結束了維系了20多年的婚姻，兩年后Marshal娶了加州一位著名的藝術家，艾麗亞娜后來一直沒有再婚，直到在93歲的高齡上因新冠去世。

不過Metropolis算法的命運，和自己的作者相比卻仿佛是截然相反。

在問世的30多年里，Metropolis算法雖然一直默默無聞，但是進入90年代之后就開始忽然火爆起來，在今天成為統(tǒng)計分析和數(shù)據科學不可或缺的利器。 澳大利亞的一位物理學家給艾麗亞娜做了維基介紹網頁，他說你今天在斯坦福或伯克隨便扔一石子，大幾率能砸到一個正在用Metropolis跑MCMC的。 這種流行度，要歸功于它和數(shù)理統(tǒng)計學的聯(lián)系，特別是貝葉斯統(tǒng)計學中的地位，以及現(xiàn)代計算機技術的大躍進，下面讓我們一一道來。

在對這個事件的歷史回顧中，我們或許可以看到，計算物理學家Metropolis主要貢獻其實就是提供了計算機，給這個課題奠定了物質基礎，成為發(fā)表文章的首席作者，也拿到了這個經典成就的冠名權，其實他并未做出真正具體的貢獻。

但是Metropolis也并非浪得虛名，他和另一個曼哈頓計劃的物理學家Stanislaw Ulam合作，在1949年寫了一篇經典文獻《蒙特卡羅方法》Monte Carlo Method，MC，此文的核心是，如果某個特定概率無法用數(shù)學或物理方法推導，比如一個形狀奇怪重心未知的篩子，就可以使用反復投篩子的方法估算其每面朝上的概率。Ulam把這個經驗方法名命為蒙特卡洛，這是法國南部一個小國摩納哥的城市，以博彩業(yè)出名。據說Ulam小時候總是看到叔叔去蒙特卡洛賭錢，就用這個擲色子的城市名命了這個擲色子的方法。這就是MCMC算法中第二個MC的含義。

當然實際的問題要比擲色子復雜太多了，計算量也大太多了，當時只有Metropolis的計算機能夠勝任。同時這個方法也蘊含了從特定的概率分布中隨機取樣的思想，所以這篇經典文獻發(fā)表在《美國統(tǒng)計學會雜志》JASA，這預示著這一系列不平凡的物理學思想將在統(tǒng)計學領域發(fā)揚光大。

也許有人要打破沙鍋問到底，從概率分布中隨機取樣有何用處？

就舉一個最熱門的例子吧，大家都知道輝瑞RNA疫苗在臨床實驗中的有效率高達90%以上，那么大家在激動的同時也關心，這個被輝瑞力捧上天的偉大有效率，它的誤差是多少？

疫苗有效率是由兩個數(shù)字來決定的：試驗疫苗組中的新冠感染數(shù)，和對照組的感染數(shù)。前者越小，疫苗的效力就越高。那么計算疫苗保護率的誤差就有一個很簡單直觀的辦法：疫苗組和對照組感染數(shù)是兩個隨機變量，如果我們更夠找到它們的統(tǒng)計分布的話，就可以從這兩個分布中大量隨機取樣，讓兩者的差值除以對照組的隨機感染數(shù)，得到的就是疫苗保護率的估計值，然后對得到的這一大堆估計值進行掐頭去尾，就得到了所謂了95%的置信區(qū)間，這是一種在臨床試驗中最常見的統(tǒng)計量。

輝瑞疫苗有效率的95%的置信區(qū)間

再回到古老的布馮投針試驗，那也是一種隨機取樣，只不過是占了針末端向四周各個方向以均等概率倒下的便宜，所以這個取樣不需要計算機，只需要牛頓萬有引力自個干活就成了。

但并不是任何的隨機取樣都是如此容易得到的。舉個最簡單的例子，怎樣獲取一系列成年人類身高的隨機數(shù)呢？

如果采用最簡單的均勻取樣的方法，類似布馮投針那樣，你也許能得到一列這樣的“身高”樣本（以厘米為單位）：

140，146，152，158，165，169，172，179，185，189，194，197， 202， 209， 215

明眼人都能看出，這列數(shù)固然“隨機”，但并不能代表真實世界中人類的身高，因為人的身高屬于正態(tài)分布，中等身材的人多，特高特矮的人少，象上面那樣的2米以上的人和170左右的人一樣多，不可能。真實的人類身高隨機取樣應該是象下面這樣的，中間多，兩邊少。

157， 162， 165， 168， 171， 172， 173，175， 177， 179， 184， 192

怎樣才能得到符合正態(tài)分布的隨機數(shù)？在沒有計算機的年代，人們會去下苦功夫在實際世界中進行大量的測量。比如在本文的上集，《借茅臺院士的熱度，科普這樣一位啤酒總工》一文中，我記錄了三位有師承關系的統(tǒng)計學大師的事跡，他們是皮爾森，戈塞特和費舍爾。其中戈塞特是著名的Student’s t test的發(fā)明人，皮爾森是他的老師，其拿手好戲就是派學生在民間大量采樣，比如監(jiān)獄犯人或者軍人的身高體重胸圍等生理參數(shù)，然后依此繪制出完備的統(tǒng)計曲線和圖表。

現(xiàn)在有了計算機，產生隨機數(shù)就太簡單了。因為所有統(tǒng)計學軟件都儲存了完備的概率分布數(shù)學解析式，一切照公式而行即可。

但關鍵是，你要知道公式才行。

在所有的學科分支中最流行的正態(tài)分布，它有一個簡潔而優(yōu)美的數(shù)學公式：

這個公式是高斯推導出來的初始樣子，和現(xiàn)代課本中的略有不同。為了紀念先賢，特保留原樣，所以正態(tài)分布也叫高斯分布

即便是這個最常用的公式，也是歷經了從伽利略，棣莫弗， 到拉普拉斯，高斯，幾代數(shù)學家歷時300多年的摸索才拿到的。

給一個統(tǒng)計變量的分布找到解析解非常困難，我們可以再舉那個新冠疫苗有效率的例子。

我在前面提到疫苗臨床試驗中，疫苗組和對照組的感染數(shù)都符合某個特定的概率分布，當時賣了個關子，沒提究竟是什么樣的分布，其實簡單的很，就是常用程度僅次于正態(tài)分布的二元分布， binomial distribution。因為它的結果是不連續(xù)的，不是0就是1，感染或者不感染。疫苗有效率是1減去兩個感染率的比值，這是一個簡單的算術操作。

可讓人不可思議的，疫苗有效率，這個對兩個最簡單統(tǒng)計變量進行的最簡單的算術操作，任憑你是牛頓還是高斯的數(shù)學天才，硬是給它找不到一個解析解，也就是一個象正態(tài)曲線那樣的公式，術語叫closed-form。沒有解析解就無法直接按公式進行直接取樣，而只能采用數(shù)值模擬的辦法（相比之下，兩個正態(tài)分布變量之比例，就有一個解析解，它在物理學上叫洛倫斯分布，在數(shù)學上叫柯西分布）。

另一個概率分布公式之難求的例子，就是我們上文提到的“啤酒”總工，威廉戈塞特，愛爾蘭吉尼斯啤酒集團的總工，也是Student’s t 檢驗的發(fā)明者。他在啤酒釀造過程中，發(fā)現(xiàn)小樣本實驗的均值和標準方差之間的比值是一個有規(guī)律的統(tǒng)計變量分布，這被后世名命為Student’s t 分布，取自當年戈塞特發(fā)表論文的筆名：Student。

戈塞特是數(shù)理統(tǒng)計鼻祖皮爾森的學生，也是牛津數(shù)學系的高材生，但就是他也沒有能力給出t分布的解析公式。是戈塞特的晚輩費舍爾（R.A.Fisher）這位數(shù)學天才，在幾年后解決了這個問題，也許是直接的代數(shù)推導不易，他采用的是高維幾何的方法。對此戈塞特不明覺厲，但直覺上知道晚輩是對的。費舍爾得到的t分布的解析解是長這個樣子的：

明顯比她的母分布正態(tài)曲線要復雜太多了。

作為實踐大師的戈塞特，在費舍爾的答案出來之前，他的解決方案是不厭其煩地從大量實驗中測出各種型態(tài)t分布的概率和相應的關鍵值，以便啤酒車間的實驗人員按圖索驥依數(shù)查表，做出合乎科學的決定：這罐啤酒母液要不要倒掉？

從某種意義上，戈塞特和他的工作人員們類似于布馮投針實驗中人肉大頭針。在大規(guī)模的科學計算中，這是一種不可持續(xù)的玩法。

由此可見，做一個“蒙特卡洛”實驗，從一個缺乏解析解的統(tǒng)計分布中隨機取樣，是不容易的。如果我們還記得的話，Marshal和艾麗亞娜夫婦的工作，恰恰通過馬可夫鏈的數(shù)值方式，用計算機從未知分布中隨機取樣，這就預示了他們這個工作在未來的意義。

我們也會看到，在解析解未知的t分布中取樣，或者是在根本就沒有解析解的疫苗有效率的分布中取樣，他們的難度和另外一類的問題相比，根本就不在一個層次上。

就在法國人布馮提出投針實驗這個天才構想的15年前，有人在英國皇家學會年會上宣讀了一篇文章，有些不同尋常的是，這篇文章的作者已經在兩年前仙去了，他的名字叫托馬斯貝葉斯（ Thomas Bayes）。在今天的統(tǒng)計學，工程和醫(yī)學制藥中，這個名字已經是如雷貫耳了，因為他被認為是貝葉斯統(tǒng)計學派的開山鼻祖。

• 你如果有耐心通讀輝瑞疫苗臨床試驗的protocol的話，你會發(fā)現(xiàn)其中的統(tǒng)計學部分幾乎全部是用貝葉斯的算法和語言寫成；

• 美國知名的政治觀察家Nate Silver，他的成名作就是使用貝葉斯模型成功預測美國2008、2012的大選結果；2020的美國大選，著名的《經濟學人》雜志特邀了當今貝葉斯計算界的權威Andrew Gelman，全盤使用貝葉斯方法進行預測，取得了成功（雖然高估了拜登取勝的盤面，不過這個是民調數(shù)據而不是概率模型的問題）；

• 2009年，法航447在從巴西飛回巴黎的途中墜落大洋，2年后，搜尋者借助貝葉斯方法在茫茫南大西洋4000米的海底找到了飛機殘骸，找到了墜機原因；

• 早在第二次世界大戰(zhàn)期間，盟軍就在貝葉斯名家Edward Simpson的幫助下，使用該方法成功破譯了納粹德國的密電碼神器Enigma，破壞了德軍的重大軍事行動。

• 就連不才小編如我，在這個以為生命科學和制藥為主的公共號上，也有數(shù)篇文章在暢談貝葉斯概率。

幾個星期之前，網上有一篇奇文在流行，號稱是用貝葉斯模型證明了新冠病毒來源是人工合成并泄露的概率高達98%。此文大家讀來都不明覺厲，因為生物學家覺得作者是一個統(tǒng)計高手，而統(tǒng)計學家覺得他們病毒專家。

后來讀者請貝葉斯算法的權威Andrew Gelman評論，他說此文的生物學他不懂，但是其貝葉斯部分是胡扯。作者把貝葉斯分析摻和進來，分明是用來拉大旗作虎皮嚇唬人的。

難道使用了貝葉斯就是高大上了嗎？歷史上并非如此，實際上在很長一個歷史時期內，貝葉斯概率思想被認為是離經叛道的。因為，就在250年前在英國皇家學會上宣布的那篇論文，托馬斯貝葉斯第一次提出了“反向概率”（reverse probability）的思想。

再用新冠疫苗當例子，有一個我們所感興趣的參數(shù)是：接種疫苗后依然感染的概率，P（感染 | 接種疫苗）。接種疫苗會誘發(fā)人體的免疫性，從而降低感染率和死亡率，這兩者有直接的因果次序，所以這個概率是符合常理的。但是貝葉斯在擺弄各種條件下的概率換算時，卻遇到了這么一個問題：在所有的已感染病例中，有多少是在之前接種了疫苗的？P（接種了疫苗 | 感染）。

在貝葉斯死后才見天日的這篇經典文獻中，他給出了一個被后世尊稱為貝葉斯定理的概率換算，如果換成新冠疫苗的例子就是這樣的：

P（感染 | 之前接種了疫苗）= P（之前接種了疫苗 | 感染）X 群體感染率 / 群體接種率

顯然，接種疫苗會降低感染率，但是今天感染不會影響昨天的接種行為，所以在等號右邊的第一個概率，在表達次序上有本末倒置之嫌，有人叫它“反向概率”。這個概念在認知哲學上的真正含義，直到今天還有爭議，就遑論200多年前的人了。

同時，這也是一個極容易引發(fā)誤導的概念。舉個例子：即使在疫苗普及后，依然也會有極小部分人拒絕接種。這個時候，你肯定希望使用公式左邊的那個概率去說服他們，因為接種疫苗后的新冠感染率和死亡率都會變得很低很低。但是它的“反向概率”卻可能把你嚇一跳，因為在確診病例中，會有很大比例也都是接種過疫苗的。

這個能作為反對疫苗的理由嗎，絕對不能。具體原因自己去想，因為這不在本文主題之內。

在歷史上，貝葉斯本人只是給后世的貝葉斯學派開了個頭，這個領域內第一位真正的大師是法國的拉普拉斯，Laplace，他的名字在前文出現(xiàn)過一次，他在高斯之前為推導正態(tài)分布的解析解做出過開拓性貢獻。但是貝葉斯體系在當時的爭議太大了，在拉普拉斯去世之后，他的朋友建議在訃告悼詞中免談死者在貝葉斯領域的工作，因為“何必給逝者臉上抹黑呢”？

在20世紀的大部分時間，貝葉斯技術雖然偶露崢嶸，比如幫忙破解了納粹德國的電報密碼，但是其應用基本處于冬眠狀態(tài)，只有少數(shù)幾個理論家在默默耕耘。這主要的原因是著名的費舍爾，R.A.Fisher，他不是貝葉斯的粉絲，無法接受這個違反常識的所謂“反向概率”，再加上此人是數(shù)學天才，我們也許記得他超越了前輩戈塞特而獨立推導出了t分布的解析解。后來費舍爾幾乎用一己之力奠定了以常規(guī)概率為藍圖的統(tǒng)計方法之數(shù)學基礎，他的學派史稱概率學派（Frequentist），把貝葉斯學派整整壓制了半個世紀。

你很快就會看到，是Marshal和艾麗亞娜夫婦的Metropolis算法讓貝葉斯門派滿血復活。

貝葉斯定理即便是在今天也非常有實用價值的。假如你想用接種疫苗后的極低感染率來給輝瑞Moderna做宣傳，你會驚訝地發(fā)現(xiàn)這個精確的統(tǒng)計數(shù)字是找不到的，因為這要求全國每一個人的接種情況和感染狀態(tài)都得到精準掌控，至少美國政府沒有能力做到這一點。

最簡單的辦法是用貝葉斯定律換算。

由于感染人口遠遠小于全體人口，所以相對容易統(tǒng)計確診者之前的疫苗接種情況（在核酸檢測表中加這么一個問題就夠了）。至于貝葉斯定律中的總確診率和總接種率呢，這個大概在每天的新聞中都能找到。這樣接種后的感染危險就算出來了（希望是大大地降低了）。

也許更重要的是，貝葉斯定理中蘊含了統(tǒng)計學的核心思想。

讓我們再復習一下新冠中的貝葉斯定律：

P（感染 | 接種過疫苗）= P（接種過疫苗 | 感染）X 群體感染率 / 群體接種率

在疫情控制中，我們最關心的是感染率的概率分布，因為我們需要這個參數(shù)來評估未來醫(yī)院的承受力，和經濟重新開放的政策，顯然這個概率是和疫苗接種情況相關的，所以公式的左邊被叫做后驗分布（Posterior Distribution），因為它是在獲知疫苗接種情況之“后’才獲得的有條件的概率分布。

在研究一個函數(shù)的變化趨勢的時候，最直接的方法是對其求導（微積分的范疇）。而在貝葉斯公式的分母中，群體接種率中不包含感染率這個參數(shù)，因此它的導數(shù)就是一個常數(shù)。所以貝葉斯在這里，把一個概率分布的核心寫成了兩個分布式的乘積。其中，整體感染率是在獲知疫苗接種情況之前的分布，被稱之為“先驗分布”（Prior Distribution），在某些場合下也可被稱為人類在沒有客觀數(shù)據時對未知事物的主觀判斷。另一個概率叫做似然函數(shù)。

如此以來，人類就獲得了幾乎無限之多的獲取概率分布的方法。

但是老問題依然沒有解決，我們或許還記得，兩個二元分布變量的簡單除法沒有解析解；從正態(tài)分布變量的簡單算術中拿到的t分布，聰明如戈塞特的都沒有拿到解析解，最后是靠費舍爾這個數(shù)學天才用高維解析幾何的方法解決的。

貝葉斯公式也面臨類似的困境，把一個未知概率分解為兩個已知分布的乘積，并不意味著容易拿到閉合形式的數(shù)學解析解，而沒有公式就無法從中取樣，對未來的疫情做推斷和預測。更何況，在真實的貝葉斯模型中參數(shù)可以是層層嵌套的，因此后驗分布的表達寫成幾十個概率分布的乘積也有可能，那樣的話就是一百個費舍爾來也是無能為力了。

這個難題在1953年看到了第一束希望的曙光。

我們也許還記得，在Metropolis算法中， Marshal Rosenbluth給粒子模擬運動設立的規(guī)則是，允許粒子移到能級較低的位點，但如果目的地的能級變高就說明這個運動給體系帶來不穩(wěn)定因素，那么就得打回原點重選方向。

Marshal也許沒有想到，這個規(guī)則簡直就是給貝葉斯的后驗概率度身定做的。艾麗亞娜寫的計算機程序既可以模擬中子的運動，也可以構建數(shù)值取樣的馬可夫鏈，其尾端的移動方向，可以取決于各個候選值之間后驗分布的比較，這和玻爾茲曼場中能級來決定運動是一個道理。

讀者也許要問了，我們不是還不知道這個后驗分布的真面目嗎？的確如此，但是，未知后驗分布和兩個或數(shù)個已知概率的乘積成正比， 所以馬可夫鏈的運動方向可以通過貝葉斯定律的換算而確定。

正如Marshal Rosenbluth證明了，遵循能量最低原則的粒子計算機隨機模擬運動，最終無限向麥克斯韋玻爾茲曼分布收斂；而在先驗概率和似然函數(shù)乘積指引下馬可夫鏈，在貝葉斯定律這只無形之手的牽引之下，也會無限地向真正的后驗分布收斂，這個蒙著面紗的美女就真相大白了。

十幾年后統(tǒng)計學科班出身的人又對這個算法做了必要的補充和發(fā)展，所以這個算法現(xiàn)在叫Metropolis-Hasting算法，簡稱M-H算法。

進入了90年代后，廉價高性能的個人電腦開始走入千家萬戶。以MCMC方法為基礎的貝葉斯學派才第一次走出了象牙塔，它不再是數(shù)學家們的專利，人們不必在無窮無盡的求導和積分中窮經皓首，而是可以借助強大的計算機，靠馬可夫鏈的隨機取樣，去解決實際中的問題。

貝葉斯學派終于咸魚翻身了。

在接近本文結尾的時候，我們再來欣賞一下M-H方法暴力的計算美學。今天的貝葉斯玩家們，也許大部分都讀不懂Marshal Rosenbluth證明和推導了，但是人人都站在巨人的肩膀上，也就是艾麗亞娜當年用無數(shù)的0和1堆砌起來的大廈。

這個例子很簡單，讓我們觀察到的數(shù)據的似然函數(shù)是一個正態(tài)分布，它的期望值和方差的先驗分布分別又是正態(tài)分布和均態(tài)分布，我們的任務是從這些數(shù)據和假設中推導出兩個正態(tài)參數(shù)的后驗分布。如前所述，即使是這樣的簡單例子也沒有解析解，唯有使用MCMC的數(shù)值模擬辦法。

左圖顯示的是三條獨立的馬可夫鏈，可以看出他們的起始位點頗為不同，但是在M-H算法的拉動下，三條鏈很快就收斂了。右圖是在馬可夫鏈上對均值和方差的隨機取樣，繪制成三維的概率密度圖，引人矚目的是，在模擬的初始階段，峰型非常不穩(wěn)而且四處位移，但是隨著馬可夫鏈走到了五千步以上，這個二維的分布開始收斂成一個穩(wěn)定漂亮的峰型。

?

https://blog.revolutionanalytics.com/2013/09/an-animated-peek-into-the-workings-of-bayesian-statistics.html

事就這樣成了。

在MCMC計算中，維度更高的復雜概率問題也是以此類推，我們也許永遠不能給他們寫出一個優(yōu)美的數(shù)學表達式，但是卻可以通過隨機采樣的方式給它畫一個逼真的肖像。

這就是概率和模擬的優(yōu)美之處，我們今天要講的這個科學史上的傳奇，就講完了。

當然，當時Metropolis算法文章的五位共同作者，都沒有意識到這個技術所蘊含的普世意義，所以他們日后都沒有在研究中再使用過這個方法。

Metropolis后來繼續(xù)從事計算科學在物理中的應用， 他獲得了該算法的冠名權，這一榮譽足以讓他名垂青史。

作為學界權威的泰勒，在氫彈成功后，大概不愿意繼續(xù)留在洛斯阿莫斯試驗室了，這畢竟是奧本海默的勢力范圍。泰勒夫婦搬到了加州伯克利，創(chuàng)立了勞倫斯利沃莫國家實驗室，和洛斯阿莫斯分庭抗禮。到了麥卡錫時代，泰勒去國會作證，還把老戰(zhàn)友奧恩海默給賣了，這是題外話。

Marsha Rosenbluth后來也搬到了加州，專業(yè)方向變成了凝聚態(tài)物理，也和這個他們點燃了第一簇篝火的領域拜拜了。直到2003年，在學術界聚集慶祝Metropolis算法問世50周年的大會上，他強撐病體前往，做了一個對歷史的回顧，厘清了理論工作貢獻的來龍去脈。人們開始意識到，這個算法的真名或許應該叫Rosenbluth-Teller算法。

艾麗亞娜的貢獻是最不為人知的。

幾年前，有人為了紀念這個劃時代算法，給這篇文章的五位作者做了一個“全家?！保鳛楣诿叩腗etropolis自然是畫面的中心，以藝術氣質見長的泰勒正在彈鋼琴，他的太太坐在一旁聽得入迷，年輕的Marsha Rosenbluth也是一副春風得意的氣度，唯有艾麗亞娜的影像是一個黑影，因為當時連網上都找不到她的影像信息。

除了Metropolis算法，統(tǒng)計界還特別推崇那篇1949年Metropolis和Stanislaw Ulam發(fā)表的《蒙特卡洛方法》。由于Metropolis已經備極哀榮了，貝葉斯界把他們最新創(chuàng)立的一門專用于馬可夫鏈隨機取樣的計算機語言名命為Stan， 為了紀念當年的共同作者之一Stanislaw Ulam，但是Stanislaw寫過程序嗎？無人知曉。

回首當年，Stan語言的創(chuàng)立者之一Andrew Gelman不無遺憾地說：也許我們應該把它名命為艾麗亞娜。

（圖片來自網絡）

參考資料：

A History of the Metropolis-Hastings Algorithm.

https://www.jstor.org/stable/30037292?seq=1

Marshall Rosenbluth and the Metropolis algorithm

https://aip.scitation.org/doi/10.1063/1.1887186

Arianna Rosenbluth Dies at 93; Pioneering Figure in Data Science

https://www.nytimes.com/2021/02/09/science/arianna-wright-dead.html

https://cosx.org/2013/12/1st-bayesian-statistics-meeting/

（如果喜歡本文，請點擊最下方六角形的“在看” Wow）

END

通知

備用：如果以后讀者無法再閱讀到本公號的更新，可以去網易搜索同名的網易號。海外的讀者也可以下載Telegram , 然后再手機瀏覽器中打開鏈接“https://t.me/joinchat/MvXTABj7X6uQxRjnNxaHmg”，點擊“ Join Group”加入“北美新藥科普歷史網”的讀者群。

Email 聯(lián)系：scimedhis@163.com