Julia 初體驗の立法會選舉勝算 DIY

我慣用 C++ 與 Matlab/Octave，偶爾也用 Python 及 R。近年眼見不少有趣語言出現，我尤為喜歡 Ruby, D（兩者其實都不算新）, Scala 及 Chapel，可是除了 Ruby 我依然計劃會抽時間學之外，其餘都只得三十秒熱度。直至最近偶然碰到 Julia，覺得應該先與她（實在無法說「它」呀）把臂同遊。

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. 郭家麒 (1.0)
2. 譚耀宗 (1.0)
3. 李卓人 (0.999993)
4. 田北辰 (0.99991)
5. 李永達 (0.99954)
6. 梁耀忠 (0.999528)
7. 陳偉業 (0.996816)
8. 麥美娟 (0.996756)
9. 陳樹英 (0.918242)
10. 余若薇 (0.821678)
11. 陳恒鑌 (0.716208)
12. 梁志祥 (0.716037)
13. 陳一華 (0.338804)
14. 何君堯 (0.338146)
15. 龍瑞卿 (0.066209)
16. 曾健成 (0.065039)
17. 譚駿賢 (0.004338)
18. 麥業成 (0.002589)
19. 陳強 (0.002561)
20. 張慧晶 (0.000565)

... 等等。一百萬次實驗，在我的電腦上需時共廿九秒。至於何謂有「極高機會」當選，何謂「較高機會」、「機會均等」、「機會較低」，就見仁見智了。

若只想知道那九位候選人有最高機會當選，那其實毋須搞甚麼模擬實驗，因為以上當選概率的高低名次，基本上與民調得出的支持度的高低排列相同。因此，模擬實驗結果中勝算最高的九位候選人，就是民調結果中支持度最高那九位。這類勝算計算的目的，其實並非要找出最有機會當選的是誰，而是要反映這些勝算較高的候選人，與其他候選人的差距。

這個方法也有它的毛病，當中最嚴重的，是它沒有考慮政黨配票的情形。以上例來說，民建聯梁志祥與陳恒鑌的支持度一直低企，約四、五個巴仙左右，可是到了選舉日，他們的得票率比起民調結果大幅上升，相反，譚耀宗的得票率就比民調結果低許多（變成約 8%）。要考慮配票，就要考慮各政黨互相鬥法，結果可能要計算隨機博奕下的 Nash equilibrium。除了較複雜之外，均衡點是否存在，是否唯一，亦造成很大的技術困難。