推 coldeye: 林澤民教授的文章有推薦你的兩篇文章 06/30 23:19
→ raiderho: 多天後補充一下：multinomial不是真的跑3000次，因為公 07/15 12:49
→ raiderho: 式很明確，所以常態分配近似算法比multinomial慢很合理 07/15 12:50
上一篇初步顯示「初選民調結果」不是小概率事件。
本篇進一步探討：
[一] 模擬可否進一步簡化？
[二] 假如同時考慮 (蔡韓柯) 與 (賴韓柯) 的民調結果，又該如何著手？
先討論 [一]:
注意到一家民調抽樣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://www.ptt.cc/bbs/Statistics/M.1560741875.A.85B.html
※ 編輯: raiderho (101.137.38.68 臺灣), 06/17/2019 22:35:58