上一篇初步显示「初选民调结果」不是小概率事件。 本篇进一步探讨： [一] 模拟可否进一步简化？ [二] 假如同时考虑 (蔡韩柯) 与 (赖韩柯) 的民调结果，又该如何着手？ 先讨论 [一]: 注意到一家民调抽样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