作者raiderho (冷颜冷雨)
看板Statistics
标题[讨论] 绿营初选民调分析2
时间Mon Jun 17 11:24:31 2019
上一篇初步显示「初选民调结果」不是小概率事件。
本篇进一步探讨:
[一] 模拟可否进一步简化?
[二] 假如同时考虑 (蔡韩柯) 与 (赖韩柯) 的民调结果,又该如何着手?
先讨论 [一]:
注意到一家民调抽样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/cn.aspx?n=bbs/Statistics/M.1560741875.A.85B.html
※ 编辑: raiderho (101.137.38.68 台湾), 06/17/2019 22:35:58
1F:推 coldeye: 林泽民教授的文章有推荐你的两篇文章 06/30 23:19
3F:→ raiderho: 多天後补充一下:multinomial不是真的跑3000次,因为公 07/15 12:49
4F:→ raiderho: 式很明确,所以常态分配近似算法比multinomial慢很合理 07/15 12:50