[R] 數學是隱喻,程式能說書


我覺得程式之所以有趣的地方,在於它是「語言」。
編程(Programming)就是在說故事。
機器語言當然相對無趣,因為那是給機器去理解的,但是對於高階語言——通用目的的C++、Java、Python,又或者是計量目的的Stata、SAS、R……——他們的特色就是更接近自然語言,也就是人類熟悉的語言。

來個腦筋急轉彎吧。

某間電影院裡頭有十排座位,A君每次去看電影都由電腦簡單隨機分配到其中一排的座位,如果A君這個樂看了八次電影,試問她(A君是個年輕女孩)會有至少一次坐跟之前一樣排數的機率是多少?

不要跟我說前兩排不坐人啦齁

這看起來是個蠢問題,因為如果我是A君當然是買中間靠走道而且絕對不買爆米花。不過這不是重點,因為這只是把一個問題生活化而已。根據心理學領域的研究,把數學或邏輯問題生活化後,人們的平均正解率就顯著增加了!準此,我也要來生活化一下。

重點是這是一個我最近正在面對的問題的簡單版本,算是複雜機械裡的一顆螺絲吧!螺絲鎖不緊,工程一定會崩垮,國家也一樣。所以說這個問題的答案是多少呢!

好吧生活化對於解決這個機率問題,或者你要說排列組合問題也好,似乎一點幫助都沒有。真的是這麼一回事嗎?

為什麼寫程式有趣,因為它可以解決問題,而且是用一種自然的方式。數學太優雅了,不是我這種低俗的人能夠領悟貫通的,但是無論如何總而言之這題的答案是這個:


好吧,看起來不算很友善。
直覺上來說,我想要找的是每一次都坐不一樣排數的事件發生機率,從集合論的角度我當然就知道它的相反即為所求。但是排列組合的數學實在離生活太遙遠了。為什麼我們不乾脆實際去這個虛擬電影院坐坐看?程式辦得到,因為在程式的世界,你就是神。(好吧我很少真的有此感受。但這個簡單的問題,我可以。

所謂機率,不高談闊論的話,指的就是「長期而言的故事」。每天都買,長期而言,你能中大樂透;長期而言,即使是大顆的粒子也能憑空穿牆(或者是穿到你的內褲裡)——但可能需要超過宇宙年齡不知道多久的時間!那麼對於電影院的問題來說,長期是什麼?長期就是,發生非常多次。非常非常多次。

從簡開始,假設只發生一次,我們如何用程式來描述電影院問題?
R這個語言之所以好用,在於它讓我們用矩陣思考。
電影院座位就是一個矩陣。下面這個程式以R的語言來描述電影院問題:



我們從數字1到10給予隨機排序,假設A君是1。(我喜歡1,因為它讓我聯想到TRUE這個邏輯變數。)那麼在矩陣draw這個物件裡,每一欄都代表來到電影院一次的座位分配。我把1獨立出來做成矩陣filter,那麼接下來我只要針對每一列進行加總,就可以得知問題的「一次性」答案。對於已經發生的事件,不存在機率,不是有(1)就是無(0)。上面這道程式執行完成時,八次看電影的事件也結束了,機率是不存在的,A君要馬就是「有至少一次坐跟之前同一排」、要馬就是「每次都坐不同排」。

看樣子,我們還無法回答這個問題。
沒關係,讓我們多試幾次,假設上面這個事件發生了一萬次呢?
如果你丟一次銅板,你要馬得到正面,要馬得到反面。如果你丟一萬次銅板,你可能有大約五千次得到正面,其餘是反面。

下面的程式小小改寫上面的故事,我們就更接近答案了:



在一萬次的「看八場電影」的事件中,有多少次是「有坐過同一排」?這個比例,聽起來似乎就是我們的答案。問題是,在隨機的環境底下,你每次進行這一萬場的模擬,得到的比例一定都會不一樣,因為這個比例,同樣也不是機率,它只是很多很多既成事實的敘述統計而已。

然而就是這樣的東西,讓我們得以逼近機率的原本。在這背後的數學理論,被叫做「弱大數法則」。
雖然我們無法獲得精確的機率數值,正如


所代表的一個實數。

它的極精確數值解同樣可以用R來算出:



用比較鬆散的方式來描述的話,這個理論上的解,我們會稱之為分析解(Analytical solution)。
在使用程式進行思想實驗來說故事的時候,我們尋找的則是數值解(Numerical solution)。
(嚴謹的分析解/數值解定義其實我沒有仔細研究過……)
我們可以把上面的程式碼再做一層推廣,模擬一百次、兩百次、三百次、……乃至於一萬次的結果所計算出來的比例,看看它的「漸進行為」長成什麼模樣。

下面的程式這麼做:



我把幾個不同模擬次數的事件比例都算出來,然後把它在一個笛卡兒坐標系統上標示出來,同時我把分析解的極精確數值近似值用紅線標示:




基本上,故事到此大致告一段落。

這張圖告訴我們,隨著模擬的規模越大(直觀上你可以理解成時間越長、或是發生越多次),我們的數值解變得越逼近理論的分析解。事實上,數學上我們也可以證明這件事情。

好吧,整個過程的啟示到底是啥咧?啟示就是,即使我的數學不是很好,我也可以用程式設法找出十分逼近解答的結果。事實上,用個數學的術語,我可以找到任意近似的解答——只要我不斷擴大模擬規模,而且運算效能允許。

這篇文章裡用到的程式碼是非常接近自然語言的東西,沒有高深的物件導向、繼承與封裝(嘿丟考試最喜歡考這個了)。它雖然在這裡只解決一個看起來單純的問題,但是,毫無疑問我也可以用它來解決複雜許多的困難問題。倒是不敢宣稱可以解決任意困難(arbitrarily difficult???)的問題就是了哈。



[UPDATE]
一開始上傳的圖片有錯,抱歉沒經檢查就丟上來,主要錯誤的地方是在第三個模擬中,我把result <- numeric(iter)這個assignment不小心拿到迴圈外頭去定義了。看來就算是這麼簡單的問題,也是可以出錯的。(默

下面這張是模擬次數更密集的版本:





最後同場加映,我把模擬的程式碼寫成函數,將座位排數與看電影次數都參數化,如此一來就可以隨便地挑選任何一個情境來看看近似的機率。



Functional Programming是R的一大特色,哪天心血來潮再來寫篇介紹吧。[狀態顯示為欠加班]

0 comment(s):

Post a Comment

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

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

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