上一篇初步顯示「初選民調結果」不是小概率事件。 本篇進一步探討： [一] 模擬可否進一步簡化？ [二] 假如同時考慮 (蔡韓柯) 與 (賴韓柯) 的民調結果，又該如何著手？ 先討論 [一]: 注意到一家民調抽樣3000人求候選人支持度的過程， 其實是求 multinomial 這個分配的 mean, 根據 CLT, mean 的分配會漸進多元常態分配。 由於 multinomial 的收斂速度較快， 三千樣本的均值應該很近似多元常態分配了。 如此，等價於問底下的近似問題： 給定一個 Multi-Normal(mu, A), 抽樣 t 次, (本例 t = 5, A 與抽樣的數目有關, 每個元素都反比於人數) 可以得到每個變數的全距。 求每個全距都不超過給定值 diff (本例 diff = 0.025) 的可能性。 我們首先檢查這樣個近似算法是否準確， 先套用上文的模型2 的參數，與模型2 的結果作比較。 好消息是，兩者模擬出來的結果非常接近： 模型2 與近似算法跑出來的結果都在 0.571 左右，後者數字更穩定。 壞消息就有點麻煩了...orz 用近似算法跑的時間竟然沒有比模型2 本身快， 事實上，花了超過四倍的時間！也就是說，在這個例子裡， 跑一次 Multi-normal 費時竟然還比跑 multinomial 12000 次還久！ 嗯，這是理論感知的效率和實測結果的差距了。 再討論 [二]: 對比式民調涉及每一個受訪者對候選人的偏好排序， 由於涉及策略性的回報 (不誠實反映偏好) 顯得更為複雜： 比方說，某受訪者的真實偏好是 蔡 > 賴 > 柯 > 韓， 此人按照真實偏好回答，會在 (蔡韓柯) 中答蔡，在 (賴韓柯) 中答賴， 然而，他可能為了讓蔡在初選民調中勝出， 會在 (賴韓柯) 中答韓或柯，或者不表態 ("唯一支持蔡英文") 考量個體的策略行為，這會讓分析變得複雜， 何況我們不知道這種策略性投票的比例 (但可以合理假定比例不大)， 估計這些參數將會更顯麻煩。 然而，我們看的是總體的比例，這可以避免個體困境。 假設在 (蔡韓柯) 情境中，蔡、韓、柯、不表態分別是狀態1, 2, 3, 4， 表態率分別是 v1, v2, v3, v4, 並記 v = (v1, v2, v3, v4)'; 在 (賴韓柯) 情境中，賴、韓、柯、不表態分別是狀態1, 2, 3, 4， 表態率分別是 u1, u2, u3, u4, 並記 u = (u1, u2, u3, u4)'. 從總體上考量，只需考慮狀態轉移的 transition matrix: 記從 (蔡韓柯) 情境的狀態j 轉成 (賴韓柯) 情境的狀態i 的機率為 p_{ij}. 令轉移矩陣 P = (p_{ij}), 那麼我們有 (*) P * v = u. P 有 16 個變數, 然而, 各行的數字和皆為 1, 因此自帶 4 條限制式. 另外，可以做以下觀察: (或者稱之為合理的假定) 若 s 取自 {韓, 柯, 不表態}, 在 (蔡韓柯) 情境中投 s 的理性受訪者，在 (賴韓柯) 情境中只會投賴或 s. 如此 p_{23} = p_{p24} = p_{32} = p_{34} = p_{42} = p_{43} = 0, 又多了 6 條限制式. 分兩種觀點： (1) 若已經事先給定 v 和 u, 由 (*) 又多了 3 條限制式，那麼只剩下 3 個自由度; (因為 v 和 u 各自分量和均為1, 所以 (*) 只有成立3 條，剩下那條自動成立。) (2) 若我們不預先給定 v 和 u, (比如向上篇的最後，考量全部合理參數空間的最低可能性)， 那麼要考慮 6 個自由度。 底下處理觀點 (1). 為了節省討論，我們就用民進黨中央給的數字做為 v 和 u, 並考量底下的可能性： (蔡韓柯) 情境不超過 diff_v = 0.25, (賴韓柯) 情境不超過 diff_u = 0.22. 我們適當取值 p_{11}, p_{31}, p_{41}, 並讓 p_{22}, p_{33}, p_{44} 都接近1. 由 (*) 的第 4 條, 注意 p_{44} 不超過 1, 會發現 p_{41} 至少是 0.1278, 由 (*) 的第 3 條, 注意 p_{33} 不超過 1, 會發現 p_{31} 至少是 0.1312, 因此，由於 p_{11} = 1- p_{21} - p_{31} - p_{41}, 這表示：不超過 74.1% 的蔡支持者會在 (賴韓柯) 選擇賴！ 隨便取 p_{11} = 0.69, p_{31} = 0.14, p_{41} = 0.15, 那麼矩陣 P 可以完全確定。 (其他數值的測試以後再做了...) 每次抽樣， 在 (蔡韓柯) 情境中抽出的狀態 s_u, 在 (賴韓柯) 會透過矩陣 P 轉移到 s_v, 由線性關係知道: (蔡韓柯) 大樣本抽出來的均數 (各候選人支持度) 會透過矩陣 P 轉移到 (賴韓柯) 大樣本抽出來的均數。 因此，由 [一]， 可以令 Multi-Normal(mu, A) 的抽出值作為 (蔡韓柯) 的各候選人支持度， 經由 P 算出 (賴韓柯) 的各候選人支持度，然後檢查： 前者各分量全距同時小於 diff_u, 後者各分量全距同時小於 diff_v 的可能性。 模擬結果： 對於這組參數，完整抽樣模擬和近似算法可能性得到的數值皆為 0.517, 也是大概率事件。 [一] 的模擬 python code: import numpy as np diff = 0.026 sample_size = 3000 agency = 5 v1, v2, v3 = 0.3508, 0.2451, 0.2270 v4 = 1 - v1 - v2 - v3 v1, v2, v3 = 0.3508, 0.2451, 0.2270 v4 = 1 - v1 - v2 - v3 mu = np.array([v1 ,v2, v3, v4]) mu_ = mu.reshape(-1,1) cov = np.diag(mu) - mu_.dot(mu_.transpose()) cov /= sample_size round = 1000000 result = 0 for i in range(round): poll = np.random.multivariate_normal(mu, cov, agency) q = np.all((np.ptp(poll[:,:-1], axis = 0)) < diff) result += q prob = result / round # prob is around 0.571, cf. the result of model 2 [二] 的模擬 python code: 完整抽樣模擬： ''' model 3: model 2 with transition matrix ''' import numpy as np diff_v, diff_u = 0.025, 0.022 u1, u2, u3 = 0.2748, 0.2347, 0.2738 u4 = 1 - u1 - u2 -u3 #u = np.array([u1, u2, u3, u4]) p11, p31, p41 = 0.69, 0.14, 0.15 p23 = p24 = p32 = p34= p42 = p43 = 0 p21 = 1 - p11 - p31 - p41 p22 = (u2 - v1 * p21) / v2 p33 = (u3 - v1 * p31) / v3 p44 = (u4 - v1 * p41) / v4 p12 = 1 - p22 p13 = 1 - p33 p14 = 1 - p44 P = np.array([[p11,p12,p13,p14], [p21,p22,p23,p24], [p31,p32,p33,p34], [p41,p42,p43,p44]]) result_3 = 0 import time start = time.time() for i in range(round): poll_v = np.random.multinomial(sample_size, para, agency) / sample_size q_v = np.all((np.ptp(poll_v[:,:-1], axis = 0)) < diff_v) poll_u = ( P.dot(poll_v.transpose()) ).transpose() q_u = np.all((np.ptp(poll_u[:,:-1], axis = 0)) < diff_u) result_3 += q_v * q_u end = time.time() print(end - start) prob_3 = result_3 / round # prob_3 is around 0.517 近似算法模擬： import numpy as np diff_v = 0.025 diff_u = 0.022 sample_size = 3000 agency = 5 v1, v2, v3 = 0.3568, 0.2451, 0.2270 v4 = 1 - v1 - v2 - v3 v = np.array([v1, v2, v3 ,v4]) u1, u2, u3 = 0.2748, 0.2347, 0.2738 u4 = 1 - u1 - u2 -u3 u = np.array([u1, u2, u3, u4]) p11, p31, p41 = 0.69, 0.14, 0.15 p23 = p24 = p32 = p34= p42 = p43 = 0 p21 = 1 - p11 - p31 - p41 p22 = (u2 - v1 * p21) / v2 p33 = (u3 - v1 * p31) / v3 p44 = (u4 - v1 * p41) / v4 p12 = 1 - p22 p13 = 1 - p33 p14 = 1 - p44 P = np.array([[p11,p12,p13,p14], [p21,p22,p23,p24], [p31,p32,p33,p34], [p41,p42,p43,p44]]) v1, v2, v3 = 0.3508, 0.2451, 0.2270 v4 = 1 - v1 - v2 - v3 mu = np.array([v1 ,v2, v3, v4]) mu_ = mu.reshape(-1,1) cov = np.diag(mu) - mu_.dot(mu_.transpose()) cov /= sample_size round = 1000000 result = 0 import time start = time.time() for i in range(round): poll_v = np.random.multivariate_normal(v, cov, agency) q_v = np.all((np.ptp(poll_v[:,:-1], axis = 0)) < diff_v) poll_u = ( P.dot(poll_v.transpose()) ).transpose() q_u = np.all((np.ptp(poll_u[:,:-1], axis = 0)) < diff_u) result += q_v * q_u end = time.time() print(end - start) prob = result / round # prob is around 0.517 --

※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 101.136.218.198 (臺灣) ※ 文章網址: https://webptt.com/m.aspx?n=bbs/Statistics/M.1560741875.A.85B.html ※ 編輯: raiderho (101.137.38.68 臺灣), 06/17/2019 22:35:58