C++ 的發明人 Bjarne Stroustrup 將 C++ 形容為 "a general purpose programming language with a bias towards systems programming"。依此說法，Julia 大概就是 "a general purpose programming language with a bias towards scientific computing" 了。若不計 Julia 語法上對矩陣的直接支援，我想一般 programmers 應該不會將她聯想成 domain-specific language 吧。
我昨天才下載 Julia，約會了一日，她給我的第一印象，是她絕對有潛力成為 Matlab 殺手或 R 殺手。無論是語法的簡潔程度、data structures 的數量、語法上對 functional programming, generic programming 及 parallel computing 的支援，抑或程式的執行速度，Julia 都明顯超越對手，網上不少 Matlab 與 R 用家亦對她頗為讚賞。她也借用了 Ruby, Python, Matlab 與 C/C++ 語法當中一些優良部份，我學習時倍感親切。
Julia 今年一月才出 1.0 版本，算係有女初長成，距離亭亭玉立還有一段日子，現在仍有不少未成熟的地方，不過已經夠足我做練習用。前文提過，計算立法會選舉各競選名單的勝算，可用多項式分佈的常態逼近當成投票的分佈。利用蒙地卡羅模擬實驗，就可以計算出各名單的勝算。文友電鋸於選舉前已經做過類似的計算，此處只是當成我第一次的 Julia 編程練習。
先說明計算細節。設 $\mathbf{p} = (p_1,\ldots,p_n)^\top = $ 民調所得各名單的支持度 ($\sum_i p_i = 1$)，而 $n$ 是樣本數，譬如 NOW 新聞台於九月七日報道（調查時窗為九月二至六日）的結果為
$$\mathbf{p}=\frac{1}{101} (4, 7, 5, 7, 1, 2, 16, 9, 3, 3, 8, 4, 1, 8, 1, 12)^\top,\ n=503.$$
（由於 NOW 新聞台四捨五入，以上向量內各數字的總和為 101 而非 100。）我們的做法，是模擬多次投票實驗。每次實驗，均由 n=503 位選民，每人隨機投一張名單一票。投票的概率由 $\mathbf{p}$ 決定。換句話說，每一名選民都會有 $p_1=\frac4{101}$ 的機會投票予第一張名單、$p_2=\frac7{101}$ 的機會予第二張名單，餘此類推。當 503 人都投完票，就可以按比例代表制查出名單上各人是否當選，查核完畢，就完成了一次實驗。重覆同樣實驗許多次 ── 譬如 1,000,000 次，就完成了整個模擬過程。整個過程當中，若排在七號名單第二順位的余若薇當選了300,000 次，她當選的概率就估計為 300,000/1,000,000 = 0.3，其他人的當選概率也用同一方式估計。
n = 503 位選民每人按 $\mathbf{p}$ 的概率來投票，即是說各名單得票 $(X_1,\ldots,X_m)$ （今屆新界西有 m=16 張競選名單）的分佈為 $\textrm{Multinomial}(n, \mathbf{p})$。多項式分佈的常態逼近公式，可參考本網誌前文，當中用到的正交矩陣 Q 的構作方法，則見我另一篇網誌。總括來說，設
$$
\begin{eqnarray}
v &=& \frac12 \left(
\begin{bmatrix}0\\ \vdots\\0\\1\end{bmatrix} -
\begin{bmatrix}\sqrt{p_1}\\ \vdots\\\sqrt{p_m}\end{bmatrix}
\right),\\
M &=&
\begin{bmatrix}\sqrt{\frac{p_1}n}\\ &\ddots\\&&\sqrt{\frac{p_m}n}\end{bmatrix}
\left(I_m - 2\frac{vv^\top}{\|v\|^2}\right).
\end{eqnarray}
$$若每次投票實驗，我們皆能夠生成 $m-1$ 個服從標準常態分佈的隨機數字 $Z_1, \ldots, Z_{m-1}$（重申，m 是競選名單數目，n 為投票人數），則每次實驗各名單的得票率可模擬為：
$$
\begin{bmatrix}\frac{X_1}n\\ \vdots\\ \frac{X_m}n\end{bmatrix}\approx \mathbf{p} + M\begin{bmatrix} Z_1\\ \vdots\\ Z_{m-1}\\ 0\end{bmatrix}.
$$
生成得票率之後，將它正規化為 $m$ 倍：
$$\mathbf{f} = (f_1, \ldots, f_m)^\top = m\left(\frac{X_1}n, \ldots, \frac{X_m}n\right)^\top,$$
之後就可以點票。正規化後的黑爾數額為 1（原本黑爾數額為 1/m，乘以 m 倍就變成 1）。每個 $f_k$ 都是一個實數，其整數部份 $\lfloor f_k\rfloor$ 代表第 k 張名單因超過黑爾數額而取得的議席數目，小數部份 $r_k = f_k - \lfloor f_k\rfloor$ 代表餘額。比較各餘額的大小，就知道餘下 $m - \sum_k \lfloor f_k\rfloor$ 個議席落入誰家。最後程式如下：
我和 Julia 還不是很熟。上面有些迴圈，相信可以用比較 functional programming 的方式寫得精簡一些。
按上述 NOW 新聞台的調查結果，由 Julia 所計算，十六張名單中各候選人（是候選人，不是候選名單）的當選概率依次為：
- 郭家麒 (1.0)
- 譚耀宗 (1.0)
- 李卓人 (0.999993)
- 田北辰 (0.99991)
- 李永達 (0.99954)
- 梁耀忠 (0.999528)
- 陳偉業 (0.996816)
- 麥美娟 (0.996756)
- 陳樹英 (0.918242)
- 余若薇 (0.821678)
- 陳恒鑌 (0.716208)
- 梁志祥 (0.716037)
- 陳一華 (0.338804)
- 何君堯 (0.338146)
- 龍瑞卿 (0.066209)
- 曾健成 (0.065039)
- 譚駿賢 (0.004338)
- 麥業成 (0.002589)
- 陳強 (0.002561)
- 張慧晶 (0.000565)
若只想知道那九位候選人有最高機會當選，那其實毋須搞甚麼模擬實驗，因為以上當選概率的高低名次，基本上與民調得出的支持度的高低排列相同。因此，模擬實驗結果中勝算最高的九位候選人，就是民調結果中支持度最高那九位。這類勝算計算的目的，其實並非要找出最有機會當選的是誰，而是要反映這些勝算較高的候選人，與其他候選人的差距。
這個方法也有它的毛病，當中最嚴重的，是它沒有考慮政黨配票的情形。以上例來說，民建聯梁志祥與陳恒鑌的支持度一直低企，約四、五個巴仙左右，可是到了選舉日，他們的得票率比起民調結果大幅上升，相反，譚耀宗的得票率就比民調結果低許多（變成約 8%）。要考慮配票，就要考慮各政黨互相鬥法，結果可能要計算隨機博奕下的 Nash equilibrium。除了較複雜之外，均衡點是否存在，是否唯一，亦造成很大的技術困難。
沒有留言:
發佈留言