[R] 坦克問題與平行運算


「真想叫發明平行運算的人去死。」
(嗚喔是設計對白)

前幾天在噗浪上有噗友討論到統計學上一個名為「德國坦克問題」的玩意兒,頗有趣的,所以就來寫個CODE模擬看看吧!(你真無趣。)

文章非常之長,想睡覺的就來吧嘿!


簡單地說呢,各位可以想像以下情境:

盟軍在一次交戰中擊破德軍四台坦克,發現坦克車身上有序號,懷疑這個序號是有順序的生產序號,因此想要藉由這幾輛摧毀的戰車序號來估計德軍的戰車數量。這幾台戰車的序號分別是2、4、8、14。

這當然是一個簡化的問題。首先針對「估計」這個字眼,我覺得值得多談幾句先!當統計學家說他要「估計」一個東西,這代表什麼呢?代表:

  1. 他不知道那個東西是多少。
  2. 他很想知道。
  3. 他決定猜猜看

好像是廢話但其實很重要真的啦。

估計,說穿了就是猜,只是很優雅地猜。既然是猜,當然可能出錯,任何人都能猜——估計——上面這個問題的答案,同時鮮少有人可以證明他絕對正確。但統計學家的猜測特別之處在於,使用數學(換言之是一種科學)的方式來進行猜測,試圖從假設、機率、風險的數學模型角度來解決一個「沒有人知道答案的問題」。

對盟軍來說,這是一個很重要的問題,所以即使是用猜的,也比什麼都不知道要來的好。
但不是亂猜。
不是亂猜。
從假設出發,依邏輯推演出數學的結果,這樣的結果在假設為真的情況下是可信的,那麼根據假設偏離現實的程度,我們也可以大概感受一下自己猜測的偏誤,心裡總會有個底。
以上是白話說詞,接下來的一切就交給數學吧!

數學能力不夠高明怎辦?
那就交給模擬吧!
下面這段R的程式碼,是一個對這個問題的粗略模擬,想先看看它的輪廓。



模擬的思維與建立數學模型其實相當接近,但卻更貼近自然語言,所以用高階語言程式碼來表達是最優雅不過了。具體來說,我們一定要先有假設,這個情境我們第一個遭遇的難題並不是德軍他媽的到底有幾台坦克,而是我們到底是怎麼擊破那四台坦克的。或許這次的會戰德軍使用的是交戰點附近兵工廠剛好生產出來的一小批坦克,而其實他大本營還有更多?或許這次來犯的坦克都是比較早生產、狀況比較不穩定的坦克,所以他們才會被擊破?或許……或許,太多的或許!所以重點是你擁有什麼情報?你能排除什麼情報?該做的都做了,剩下的就只能交給假設。

如果除了盟軍擊破四台坦克之外我什麼都不知道,那麼對我而言最直觀的一個假設就是:盟軍擊破德軍每一台坦克的機會都差不多。來點騙小孩的數學,假設在任一次會戰中盟軍擊破一台德軍坦克的序號是X,一個隨機變數,那麼我可以把它表示成:


意思就是X服從間斷均勻分配。
(更精良一點的說法:擊破坦克的序號是一串隨機變數,抽後不放回,因此前後並非完全無關,但每一次抽的時候剩餘的個體都有同樣的機會被抽中。)
有了這個假設——你可以有其他不同的假設,從而導致不同的估計——這個得軍坦克問題的數學描述就是

我們對一間斷均勻分配的隨機變數抽了一個大小為四的樣本,想要估計未知母數N——母體的大小。

已知的樣本告訴了我們什麼?還記得抽到的編號是2、4、8、14。換言之,根據假設,我們馬上知道德軍至少有十四台坦克。然後咧?如果只是停在這裡就不夠優雅了。首先,十四台坦克這個猜測絕對有問題,第一,它不可能是一個「不偏估計量」——因為樣本最大序號只會比母體最大序號小而不會比較大,所以這個估計量的算數期望值必定小於母數真值;第二,呃,好吧我沒有想到第二,因為當我知道一個猜測「平均而言總會低估(=偏了)」,我就知道它不是一個好的猜測。如何解決這個問題?Roger W. Johnson(1994)這篇文章針對這個問題談到了一個不偏估計量,內文解釋十分淺顯易懂,但不是我這邊想談的,因為它不是來自一個模擬情境的直覺所造就的猜測。

不妨先看看上面的程式碼畫出的這張圖吧!


這張圖的意思是,以橫軸為假設的N(德軍坦克總數),縱軸則是每發生一千次隨機抽樣,且每次都抽(擊破!)四台,那麼這樣的樣本裡頭觀察到最大序號是14的比例有多少。有點繞口,但其實很自然。我們的問題就是,我們只有一次會戰,但有了前述的假設,我們就可以建立機率模型,有機率模型,我們就可以產生模擬。在這個完全服從我們假設的虛構世界裡,我們可以任意地讓「隨機擊破四台坦克」這個事件不斷發生,從而去觀察這個事件的統計性質。

數學不錯的讀者大概很早就發現,這個問題何必模擬,因為如果N=14,而我們要觀察到一個樣本是最大序號m=14的機率,不就是:


顯然你可以把它想像是一個排列組合問題。順帶一提在R裡頭我們可以用下面這行



來計算它的極精確數值解。

所以這張圖就是告訴我當德軍坦克總是實際上等於各種(可能的)數字時,我們觀察到手上這個樣本(最大序號是14)的機率,的近似值。我們是用一千次重複事件來逼近這個機率真值,而它的真值我們已經知道有公式可以算得出來。事實上,我們現在要做的事情就是,考慮所有可能的N,並考慮這些所有可能的N導致我們觀察到手頭上這個樣本的機率,然後我們就可以求出一個期望值、一個機率加權平均數,作為我們的猜測。

事實上,這個猜測被稱為貝氏猜測(Baysian Estimate)。

既然每個N都有解析解,我們實際上要運算的就是一個N趨近於無限大的加總數量的加權平均。
我沒有要算的意思。
不過有興趣的人可以去參照維基百科頁面,裡面有包括了我前面提到的一個不偏估計量的推導過程,以及貝氏估計式的求導。如果你真的去看那些算式,多少會發現明明貼近自然的東西突然就變得很曲折,那是因為平均而言我們對數學語言的理解遠遠不及我們對情境語言、自己日常每天接觸的語言的理解能力。

有些人可以很自然地讀通數學語言,我不是這種人。
所以我才想要寫程式來回答這個問題。

維基百科裡有一張貝氏估計的間斷機率分配函數圖形,長得和我們上面畫出來的非常像。事實上那就是我們想要逼近的目標,貝氏機率的真值,透過模擬的方式我們產生了一個近似的貝氏模型。

接下來就是模擬的問題了!
記得我們要算的是一個無窮數列的加總加權平均,如何模擬無窮?上面的程式碼我只模擬到N=50而已,我的加權平均當然也就只考慮到這裡。聽起來模擬一個無限大的東西是根本辦不到的,實際上呢,要看問題是什麼。

再次觀察上面這張圖,很明顯,隨著N變大我們觀察到四個樣本其最大值是14的機率急劇下降,這非常直覺:「如果德軍坦克可以堆得跟山一樣高,憑甚麼我們隨便看見的那四台最大編號只到14?」我們的模擬於是出現一絲曙光,也就是收斂。顯然的,這張圖在告訴我們下面這件事情:


意思是你考慮很大很大的N是沒有意義的啦!因為很大的N——根據我們的假設——你幾乎不可能觀察到你手上的這個樣本,換言之這個無窮數列加總的後面全部都是零,你有算跟沒算是一樣滴。只是我沒有辦法告訴你「後面」到底是有多後面,但這個問題可以透過模擬來回答。

所以,上面那段程式碼給我的猜測到底是多少?差不多18.45。



從維基百科我們已經得知此題的貝氏猜測是19.5,換言之,如果我們的模擬夠好,應該要能逼近這個數字。不曉得各位怎麼看,但我覺得,顯然不太妙。為了進一步追蹤我的模擬逼近真值的情形,我把程式碼稍微改寫:



這次我把N的最大可能值調高到100了,並且N每增加一次我就紀錄一次我重新計算的貝次猜測,把這些結果丟到outcome這個向量裡頭。細心的讀者還可以注意我多用了system.time()這個函數想看看運算的速度,這也是為了接下來的平行運算埋個伏筆。(該睡了)

我們來看看模擬出來的貝氏猜測隨著N逐漸增加有無更貼近真值(標上紅線):


有嗎?有嗎?
我們來print一下outcome這個向量的最後幾項:



看起來,還是不太妙。它告訴我N最後六次增加(從95到100),對我的猜測值沒有任何影響。花黑噴?絕對不是程式碼有蟲,這個我可以保證,雖然說這句話的時候聲音在抖也一樣。

其實問題在於我們的蔬菜不是是基數不夠。由於樣本太小(k=4),對於每一個N,我們都需要相當大的重複次數才能夠逼近到真值,這些真值彼此的加權平均才是最終目標,所以毫無疑問我們會離目標有那麼一段無法彌補的距離。

解決方案?提高事件模擬次數。目前是一千次,假設五千次呢?在那之前,我們可以先來驗證一下前面的想法,將情境暫時改成k=8,看看結果如何,我就不重覆秀出來了,模擬出來的貝氏猜測是15.1772,而真值呢?15.1667。這是模擬的技術面問題,你必須知道需要多大的基數,才能有良好的近似。現在這個問題做為一個練習,我們總是有真值可以互相比較,但當面臨真值未明的情況,就得靠個人功力了。不過也不完全算是個人功力,因為我們總是可以採用一個相當大的基數來進行模擬,問題在於程式的效能。

讓我們回到k=4的情形,我要一次將模擬次數調高到五千,而且考慮最大的N到五百!雖然這仍然不是一個非常暴力的解決方案,但我想為了未來的發展空間(是要多暴力啦),也是時候該導入更優雅的編程技術了。這裡我想用的是單機多核平行運算,在R 2.14版之後已經內建了一個名為parallel的套件,是基於snow與multicore這兩個套件的平行運算框架,不過這裡我想直接使用snow這個套件,因為它的名稱是如此優雅(WTF)。

除此之外,記得每次模擬我們總是希望能夠給定隨機數字產生器的種子,這樣我們才有辦法重製結果。但在平行運算的框架底下,這會有點問題。由於各個worker(相對於master)的session啟動時間就關連到預設的隨機種子,可能產生部分平行session隨機種子相同、部分又不同的問題。這邊我裝了rlecuyer這個平行運算隨機數字產生器的套件來解決這個問題,它是snow內建函式仰賴的兩個平行隨機數字產生器之一。

同樣地,為了能夠追蹤貝氏猜測隨著N膨脹而修正的軌跡,我的程式碼會需要一些額外的步驟,所以會變得稍微複雜一些:



為了進行平行運算,必須把模擬的整個程序寫成函數。要記得一個關鍵:每個worker session並不會分享master session的變數域!換言之,你在master session(就是你現在正在看的這整塊程式碼)該函式外定義的物件,worker都無法直接取用!你有兩個方法可以解決這個問題,其一是將會用到的物件全部都寫進函式的引數(我在這裡採用此法);其二是使用clusterExport()套件函式將物件傳遞給worker。

makeCluster()負責定義你的工作群組,我使用Type="SOCK"的意思是採用最原始的TCP/IP協定來進行master-worker之間的溝通,由於是單機平行,所以我指派的四個worker都是我的本機電腦(請量力而為)。

將演算法平行化的關鍵在於:辨認出哪個部位是可以而且值得被平行化的。我們這邊面臨的平行化問題,在平行運算的領域被稱之為Embarrassingly Parallel,是的,意思就是真他媽簡單平行運算。這在裡我到底做了什麼呢?我把我們的迴圈拆開來了!事實上,拆成四分,分給四個worker同時計算,然後再把結果收集起來(這部分已經由snow套件內部封裝實作完成,我們不必去煩惱),就這樣。

由於我的迴圈是針對每個N(怕你忘記,德軍坦克總數)去做,這個計算對每個N基本上都是獨立的,這是非常標準的真他媽簡單平行運算的範疇,但別看這好像不需要什麼技術力,他提升的速度是非常驚人的。我上面的程式碼只需要跑十秒鐘,如果不採用平行運算,嗯……因為我根本不想等,所以也沒去實測!

好,來收割結果吧。
下面這張圖是模擬貝氏猜測的逼近情形:


最終模擬的猜測值(max(N)=500, k=4, j=5000):



SWEET!
我們還可以再挖深一點:



結果顯示,在N=282+13=295之後,我們的貝氏猜測就再也沒有修正了。(對了之所以要加13是因為過程中我把N <= 13直接拿掉了,畢竟它們的結果必然為零。) 還記得嗎?這叫作收斂。

大功告成。

等等,其實,不該在這裡就停下來。
如果你真的想要在沒有真值的情況下,做出可信的模擬結果,我們通常需要更多。
在這個例子裡,我除了想知道隨著N越大,是否猜測值有收斂行為之外,我還想知道隨著試行次數j越大,我們的猜測是否同樣收斂。

這次我們需要真正暴力運算的模擬了。

這次我想要看的是,隨著不同的j值,我們最終的貝氏猜測值的「漸近行為」。我的預期是,結果應該要向真值收斂。但我不知道要多嚴格的模擬參數設定才能讓我清楚看到這樣的收斂。更重要的是,由於每一次的模擬——不論你設定了多大的基數——都受到隨機的影響,因此任何一次模擬的結果都會包含隨機雜訊,所以真正可靠的模擬,就是採用大量模擬的平均值作為我們的猜測,這個取平均的行為能夠把許多雜訊排除。

上面的模擬結果非常接近真值,但我懷疑這只是我們運氣好,所以這次我要做好幾次更嚴格的模擬然後看平均的結果。以我們的問題來說,這就是真正的暴力運算了,以下是我要做的:



這次我要把N考慮到一萬!(德軍坦克可能會有一萬台嗎?讓我們繼續看下去~)
然後對於模擬試行次數j,我要j=1000做一次,j=1250做一次……一路做到j=10000!

This is tough.
確定你的電腦有那個能耐再嘗試執行上面的程式。
以下是我畫出來的圖:


或許讀者覺得他不夠漂亮,看起來似乎沒有「收斂」的特性。
事實上,這是因為隨機干擾,我們很快就發現每一次的結果總是環繞著真值(紅線)上下飄動。如果我們對所有結果取平均:



非常接近真值。
從數學上我們根本可以證明我們的模擬貝氏猜測值會「機率收斂」到真值,這也是為什麼我如此有信心的原因。具體來說,對於每個不同的N的抽樣結果,都是「弱大數法則」的力量在運作,而我們要逼近的真值基本上是這許許多多不同的N對應結果的一個函數(還記得就是把無窮的它們加總起來算期望值),有一個名為Continuous Mapping Theorem的玩意兒,確保我們這麼做是可行的。



最後附帶一提,除了均勻分配之外,我們其實還有一個真正賴以為生(?)的假設,就是相信戰車編號反映了得軍坦克實際的序號。你說德軍兵工廠序號都是隨機出貨?抱歉那我們就GAME OVER惹。這又牽涉到另一個有趣的機率議題,叫做「主觀機率」(事實上,貝氏機率就是一種主觀機率)。你覺得德軍兵工廠序號竟然是隨機出貨的可能性有多高?隨便你說,反正我是不信。(喂)

最後的最後,大概會有人覺得我在搞笑,既然都有貝氏猜測的精確數學解析解了,那我寫這麼多程式碼來「逼近一個猜測」的意義何在?
意義就在於,並不是所有問題都能讓你找到清楚、可以操作的解析解。
問題的複雜度可能使你的數學模型無法產生解析解,甚至根本不知道如何建立一個模型,這個時候你能怎麼辦?事實上在當代總體經濟學領域,有非常多的模型都沒有學者慣稱之為Closed Form Solution的結果,這使得計算統計學(Computational Statistics)變得分外重要。

意義就在於,你總是可以回歸自然,建立模擬。
年輕人令人稱羨之處就在於,他們擁抱的,是貼近自然的語言,同時卻又享有摩登的效能!

1 comment(s):

林威成 | 05 February, 2016 09:51

覺得厲害~~

Post a Comment

回應文章前請注意下列三勿原則:

1)勿拍照;(→會有靈異的照片從你的相機裡跑出來...
2)勿餵食;(→會有飢渴的猛獸從我的網誌裡跑出來...
3)勿告白。(→會有奇怪的東西從站長體內裡跑出來...

謝謝大家的配合。
( > ー <)b