利用贝叶斯最优实验设计预测美国总统选举结果¶
在本教程中，我们将探讨使用最优实验设计技术创建最优抽样策略，以预测美国总统选举的结果。在之前的教程中，我们探讨了使用贝叶斯最优实验设计来学习单个人的工作记忆容量。在这里，我们将相同的概念应用于研究整个国家。
首先，我们需要一个关于选举获胜者 w 的贝叶斯模型，以及我们可能计划进行的任何抽样的结果 y。实验设计是在每个州进行抽样的人数 \(n_i\)。为了建立我们的探索性模型，我们将进行一些简化假设。我们将使用 1976-2012 年的历史选举数据构建一个合理的先验，并使用 2016 年的选举作为我们的测试集：我们设想我们正在 2016 年选举前进行抽样。
选择先验¶
在我们的模型中，我们包含了一个 51 维的隐变量 alpha。对于 50 个州加上哥伦比亚特区的每个州，我们定义
我们假设所有其他选民投票给共和党。选举前夕，\(\alpha\) 的值是未知的，我们希望通过对州 \(i\) 的 \(n_i\) 个人进行抽样来估计它，其中 \(i=1, ..., 51\)。选举的获胜者 \(w\) 由选举人团制度决定。民主党在州 \(i\) 获得的选举人票数是
（这是对真实系统的一个粗略近似）。所有其他选举人票都归共和党。这里 \(k_i\) 是分配给州 \(i\) 的选举人票数，这些票下列在以下数据框中。
[3]:
# Data path
BASE_URL = "https://d2hg8soec8ck9v.cloudfront.net/datasets/us_elections/"
[4]:
import pandas as pd
import torch
from urllib.request import urlopen
electoral_college_votes = pd.read_pickle(urlopen(BASE_URL + "electoral_college_votes.pickle"))
print(electoral_college_votes.head())
ec_votes_tensor = torch.tensor(electoral_college_votes.values, dtype=torch.float).squeeze()
Electoral college votes
State
AL 9
AK 3
AZ 11
AR 6
CA 55
选举的获胜者 \(w\) 是
在代码中，表达如下
[5]:
def election_winner(alpha):
dem_win_state = (alpha > 0.).float()
dem_electoral_college_votes = ec_votes_tensor * dem_win_state
w = (dem_electoral_college_votes.sum(-1) / ec_votes_tensor.sum(-1) > .5).float()
return w
我们对能够帮助我们预测 \(w\) 的抽样策略感兴趣，而不是预测更复杂的逐州结果 \(\alpha\)。
为了建立一个完整的贝叶斯模型，我们需要 \(\alpha\) 的先验。我们将基于一些历史总统选举的结果来建立先验。具体来说，我们将使用以下包含 1976-2012 年总统选举逐州结果的数据集。请注意，除了民主党和共和党以外的投票已被忽略。
[7]:
frame = pd.read_pickle(urlopen(BASE_URL + "us_presidential_election_data_historical.pickle"))
print(frame[[1976, 1980, 1984]].head())
1976 1980 1984
Democrat Republican Democrat Republican Democrat Republican
State
AL 659170 504070 636730 654192 551899 872849
AK 44058 71555 41842 86112 62007 138377
AZ 295602 418642 246843 529688 333854 681416
AR 499614 268753 398041 403164 338646 534774
CA 3742284 3882244 3083661 4524858 3922519 5467009
仅基于这些数据，我们将把 \(\alpha\) 的先验均值完全基于 2012 年的选举结果。我们的模型将基于逻辑回归，因此我们将使用 logit 函数转换投票给民主党的概率。具体来说，我们将按如下方式选择先验均值
[8]:
results_2012 = torch.tensor(frame[2012].values, dtype=torch.float)
prior_mean = torch.log(results_2012[..., 0] / results_2012[..., 1])
我们对 \(\alpha\) 的先验分布将是具有均值 prior_mean 的多变量正态分布。唯一需要决定的就是协方差矩阵。由于 alpha 值经过 logit 转换，协方差也将定义在 logit 空间中。
旁注：先验协方差在许多方面都很重要。如果我们允许方差过大，先验将对每个州的结果不确定，并需要在各地进行抽样。如果我们允许方差过小，我们可能会因意外的选举结果而措手不及。如果我们假设各州是独立的，那么我们将无法跨州汇集信息；但如果假设相关性过高，我们可能会过于忠实地根据一个州的抽样结果来预测另一个州。
我们通过计算 1976 年至 2012 年选举结果的经验协方差，并在对角线上加上一个小的数值 0.01 来选择先验协方差。
[9]:
idx = 2 * torch.arange(10)
as_tensor = torch.tensor(frame.values, dtype=torch.float)
logits = torch.log(as_tensor[..., idx] / as_tensor[..., idx + 1]).transpose(0, 1)
mean = logits.mean(0)
sample_covariance = (1/(logits.shape[0] - 1)) * (
(logits.unsqueeze(-1) - mean) * (logits.unsqueeze(-2) - mean)
).sum(0)
prior_covariance = sample_covariance + 0.01 * torch.eye(sample_covariance.shape[0])
建立模型¶
我们现在可以定义我们的模型了。从高层次上看，模型工作如下
\(\alpha\) 是多变量正态分布
\(w\) 是 \(\alpha\) 的一个确定性函数
\(y_i\) 是 Binomial(\(n_i\), sigmoid(\(\alpha_i\)))，因此我们假设人们对抽样的回应方式与他们在选举日的投票方式完全相同
在 Pyro 中，这个模型看起来如下
[10]:
import pyro
import pyro.distributions as dist
def model(polling_allocation):
# This allows us to run many copies of the model in parallel
with pyro.plate_stack("plate_stack", polling_allocation.shape[:-1]):
# Begin by sampling alpha
alpha = pyro.sample("alpha", dist.MultivariateNormal(
prior_mean, covariance_matrix=prior_covariance))
# Sample y conditional on alpha
poll_results = pyro.sample("y", dist.Binomial(
polling_allocation, logits=alpha).to_event(1))
# Now compute w according to the (approximate) electoral college formula
dem_win = election_winner(alpha)
pyro.sample("w", dist.Delta(dem_win))
return poll_results, dem_win, alpha
理解先验¶
在我们进一步之前，我们将研究该模型，以检查它是否符合我们对美国总统选举的直觉。
首先，让我们看看每个州投票给民主党的选民比例的置信上限和下限。
[11]:
std = prior_covariance.diag().sqrt()
ci = pd.DataFrame({"State": frame.index,
"Lower confidence limit": torch.sigmoid(prior_mean - 1.96 * std),
"Upper confidence limit": torch.sigmoid(prior_mean + 1.96 * std)}
).set_index("State")
print(ci.head())
Lower confidence limit Upper confidence limit
State
AL 0.272258 0.517586
AK 0.330472 0.529117
AZ 0.321011 0.593634
AR 0.214348 0.576079
CA 0.458618 0.756616
\(\alpha\) 上的先验隐式定义了我们对 w 的先验。我们可以通过多次从先验中进行模拟来研究这个先验。
[12]:
_, dem_wins, alpha_samples = model(torch.ones(100000, 51))
prior_w_prob = dem_wins.float().mean()
print("Prior probability of Dem win", prior_w_prob.item())
Prior probability of Dem win 0.6799200177192688
由于我们的先验基于 2012 年，并且民主党在 2012 年获胜，因此我们偏向于民主党在 2016 年获胜是有道理的（这是在我们看到任何抽样数据或纳入任何其他信息之前）。
我们还可以初步研究哪些州是最边缘化的。
[13]:
dem_prob = (alpha_samples > 0.).float().mean(0)
marginal = torch.argsort((dem_prob - .5).abs()).numpy()
prior_prob_dem = pd.DataFrame({"State": frame.index[marginal],
"Democrat win probability": dem_prob.numpy()[marginal]}
).set_index('State')
print(prior_prob_dem.head())
Democrat win probability
State
FL 0.52501
NC 0.42730
NH 0.61536
OH 0.61997
VA 0.63738
这是一个健全性检查，似乎与我们的直觉相符。佛罗里达州经常是一个重要的摇摆州，并且在先验下处于我们边缘州列表的首位。我们还可以看到宾夕法尼亚州和威斯康星州也靠近列表顶部——我们知道这些州在 2016 年选举中起着重要作用。
最后，我们仔细看看我们的先验协方差。具体来说，我们检查我们预计相关性更强或更弱的州。让我们先看看新英格兰地区的州
[14]:
import numpy as np
def correlation(cov):
return cov / np.sqrt(np.expand_dims(np.diag(cov.values), 0) * np.expand_dims(np.diag(cov.values), 1))
new_england_states = ['ME', 'VT', 'NH', 'MA', 'RI', 'CT']
cov_as_frame = pd.DataFrame(prior_covariance.numpy(), columns=frame.index).set_index(frame.index)
ne_cov = cov_as_frame.loc[new_england_states, new_england_states]
ne_corr = correlation(ne_cov)
print(ne_corr)
State ME VT NH MA RI CT
State
ME 1.000000 0.817323 0.857351 0.800276 0.822024 0.825383
VT 0.817323 1.000000 0.834723 0.716342 0.754026 0.844140
NH 0.857351 0.834723 1.000000 0.871370 0.803803 0.873496
MA 0.800276 0.716342 0.871370 1.000000 0.813665 0.835148
RI 0.822024 0.754026 0.803803 0.813665 1.000000 0.849644
CT 0.825383 0.844140 0.873496 0.835148 0.849644 1.000000
显然，这些州倾向于以相似的方式投票。我们还可以检查一些我们同样预计会相似的南方州。
[15]:
southern_states = ['LA', 'MS', 'AL', 'GA', 'SC']
southern_cov = cov_as_frame.loc[southern_states, southern_states]
southern_corr = correlation(southern_cov)
print(southern_corr)
State LA MS AL GA SC
State
LA 1.000000 0.554020 0.651511 0.523784 0.517672
MS 0.554020 1.000000 0.699459 0.784371 0.769198
AL 0.651511 0.699459 1.000000 0.829908 0.723015
GA 0.523784 0.784371 0.829908 1.000000 0.852818
SC 0.517672 0.769199 0.723015 0.852818 1.000000
这些协方差矩阵表明，正如预期的那样，逻辑分组的州倾向于具有相似的投票趋势。我们现在看看组之间的相关性（例如，缅因州和路易斯安那州之间）。
[16]:
cross_cov = cov_as_frame.loc[new_england_states + southern_states, new_england_states + southern_states]
cross_corr = correlation(cross_cov)
print(cross_corr.loc[new_england_states, southern_states])
State LA MS AL GA SC
State
ME 0.329438 0.309352 -0.000534 0.122375 0.333679
VT -0.036079 0.009653 -0.366604 -0.202065 0.034438
NH 0.234105 0.146826 -0.105781 0.008411 0.233084
MA 0.338411 0.122257 -0.059107 -0.025730 0.182290
RI 0.314088 0.188819 -0.066307 -0.022142 0.186955
CT 0.139021 0.074646 -0.205797 -0.107684 0.125023
现在，我们看到新英格兰州和南方州之间的相关性比这些分组内部的相关性要弱。这再次是预料之中的。
衡量抽样策略的预期信息增益¶
我们建立的先验似乎至少在近似程度上符合直觉。然而，我们现在想添加第二种信息来源：抽样。我们的目标是使用我们的先验来选择最有信息量的抽样策略，以了解我们的目标 \(w\)。在这种简化设置中，抽样策略是在每个州抽样的人数。（我们忽略了任何其他协变量，例如州内的区域差异、人口统计学等）。我们可能认为在佛罗里达州（最边缘的州）对 1000 人进行抽样会比在哥伦比亚特区（最不边缘的州）对 1000 人进行抽样有效得多。这是因为哥伦比亚特区的结果基于我们的先验已经相当可预测，而佛罗里达州的结果则完全是待定的。
事实上，我们的模型基于设计 \(d\) 进行抽样并获得结果 \(y\) 所获得的关于 \(w\) 的信息可以用数学描述如下
由于抽样结果目前未知，我们考虑预期信息增益 [1]
EIG 的变分估计器¶
在工作记忆教程中，我们使用“边缘”估计器来找到 EIG。这涉及到估计边缘密度 \(p(y|d)\)。在此实验中，这将相对困难：\(y\) 是 51 维的，并且具有一些相当棘手的约束，使得对其密度进行建模变得困难。此外，“边缘”估计器要求我们解析地知道 \(p(y|w)\)，而我们并不知道。
幸运的是，存在其他 EIG 的变分估计器：更多细节请参见 [2]。其中一个变分估计器是“后验”估计器，基于以下表示
在这里，\(H(p(w))\) 是 \(w\) 上的先验熵（我们可以很容易地计算出来）。重要的项涉及到变分近似 \(q(w|y)\)。这个 \(q\) 可以用于执行摊销变分推断。具体来说，它将 \(y\) 作为输入，并输出一个关于 \(w\) 的分布。当 \(q(w|y) = p(w|y)\) 时，上界达到最大值 [2]。由于 \(w\) 是一个二元随机变量，我们可以将 \(q\) 视为一个分类器，它试图根据抽样结果来决定最终谁将赢得选举。在本 Jupyter Notebook 中，\(q\) 将是一个神经网络分类器。训练一个神经网络分类器比学习 \(y\) 的边缘密度要容易得多，因此我们在本教程中采用这种方法来估计 EIG。
[17]:
from torch import nn
class OutcomePredictor(nn.Module):
def __init__(self):
super().__init__()
self.h1 = nn.Linear(51, 64)
self.h2 = nn.Linear(64, 64)
self.h3 = nn.Linear(64, 1)
def compute_dem_probability(self, y):
z = nn.functional.relu(self.h1(y))
z = nn.functional.relu(self.h2(z))
return self.h3(z)
def forward(self, y_dict, design, observation_labels, target_labels):
pyro.module("posterior_guide", self)
y = y_dict["y"]
dem_prob = self.compute_dem_probability(y).squeeze()
pyro.sample("w", dist.Bernoulli(logits=dem_prob))
现在我们将使用它来计算几种可能的抽样策略的 EIG。首先，我们需要计算上述公式中的 \(H(p(w))\) 项。
[18]:
prior_entropy = dist.Bernoulli(prior_w_prob).entropy()
让我们考虑四种简单的抽样策略。1. 仅在佛罗里达州抽样 1000 人 2. 仅在哥伦比亚特区抽样 1000 人 3. 在美国各地均匀抽样 1000 人 4. 使用侧重于摇摆州的抽样分配
[19]:
from collections import OrderedDict
poll_in_florida = torch.zeros(51)
poll_in_florida[9] = 1000
poll_in_dc = torch.zeros(51)
poll_in_dc[8] = 1000
uniform_poll = (1000 // 51) * torch.ones(51)
# The swing score measures how close the state is to 50/50
swing_score = 1. / (.5 - torch.tensor(prior_prob_dem.sort_values("State").values).squeeze()).abs()
swing_poll = 1000 * swing_score / swing_score.sum()
swing_poll = swing_poll.round()
poll_strategies = OrderedDict([("Florida", poll_in_florida),
("DC", poll_in_dc),
("Uniform", uniform_poll),
("Swing", swing_poll)])
现在我们将计算每个选项的 EIG。由于这需要训练网络（四次），可能需要几分钟。
[20]:
from pyro.contrib.oed.eig import posterior_eig
from pyro.optim import Adam
eigs = {}
best_strategy, best_eig = None, 0
for strategy, allocation in poll_strategies.items():
print(strategy, end=" ")
guide = OutcomePredictor()
pyro.clear_param_store()
# To reduce noise when comparing designs, we will use the precomputed value of H(p(w))
# By passing eig=False, we tell Pyro not to estimate the prior entropy on each run
# The return value of `posterior_eig` is then -E_p(w,y)[log q(w|y)]
ape = posterior_eig(model, allocation, "y", "w", 10, 12500, guide,
Adam({"lr": 0.001}), eig=False, final_num_samples=10000)
eigs[strategy] = prior_entropy - ape
print(eigs[strategy].item())
if eigs[strategy] > best_eig:
best_strategy, best_eig = strategy, eigs[strategy]
Florida 0.3088182508945465
DC 0.000973820686340332
Uniform 0.30316317081451416
Swing 0.3254041075706482
运行实验¶
我们现在已经评估了我们的四种候选设计，并可以选择最优的一种来实际收集新数据。在此 Jupyter Notebook 中，我们将使用 2016 年选举的结果模拟新数据。具体来说，我们将假设抽样结果来自我们的模型，其中我们将 alpha 的值设置为与 2016 年的实际结果相对应。
首先，我们使用选定的抽样策略重新训练 \(q\)。
[21]:
best_allocation = poll_strategies[best_strategy]
pyro.clear_param_store()
guide = OutcomePredictor()
posterior_eig(model, best_allocation, "y", "w", 10, 12500, guide,
Adam({"lr": 0.001}), eig=False)
[21]:
tensor(0.3653, grad_fn=<DivBackward0>)
由 2016 年结果推断出的 \(\alpha\) 值，其计算方法与我们计算先验的方法相同。
[23]:
test_data = pd.read_pickle(urlopen(BASE_URL + "us_presidential_election_data_test.pickle"))
results_2016 = torch.tensor(test_data.values, dtype=torch.float)
true_alpha = torch.log(results_2016[..., 0] / results_2016[..., 1])
[24]:
conditioned_model = pyro.condition(model, data={"alpha": true_alpha})
y, _, _ = conditioned_model(best_allocation)
让我们看看我们抽样的结果。
[25]:
outcome = pd.DataFrame({"State": frame.index,
"Number of people polled": best_allocation,
"Number who said they would vote Democrat": y}
).set_index("State")
print(outcome.sort_values(["Number of people polled", "State"], ascending=[False, True]).head())
Number of people polled Number who said they would vote Democrat
State
FL 191.0 89.0
NE 66.0 30.0
NJ 41.0 22.0
OH 40.0 19.0
VT 35.0 22.0
分析数据¶
收集数据后，我们现在可以在模型下进行推断，以获得民主党获胜的后验概率。进行推断的方法有很多种：例如，我们可以使用 Pyro 的 SVI 进行变分推断，或使用 MCMC（例如 Pyro 的 NUTS）。使用这些方法，我们将计算 alpha 的后验，并使用此结果获得 w 的后验。
然而，一种快速分析数据的方法是使用我们已经训练好的神经网络。在收敛时，我们期望网络能很好地近似真实的后验，即 \(q(w|y) \approx p(w|y)\)。
[26]:
q_w = torch.sigmoid(guide.compute_dem_probability(y).squeeze().detach())
[27]:
print("Prior probability of Democrat win", prior_w_prob.item())
print("Posterior probability of Democrat win", q_w.item())
Prior probability of Democrat win 0.6799200177192688
Posterior probability of Democrat win 0.40527692437171936
这就完成了实验设计、数据获取和数据分析的整个过程。通过使用我们的分析模型和先验来设计抽样策略，我们能够选择产生最大预期信息增益的设计。
结论¶
在本教程中，我们展示了如何将最优实验设计原则应用于选举抽样。如果我们可以选择首先使用我们的先验来设计抽样策略，那么一旦我们收集了数据，这将有助于获得更好的预测。
我们的模型可以通过多种方式增强：[3] 提供了一份有用的指南，介绍了政治学中使用的模型类型。我们的最优设计策略的质量永远取决于我们先验的质量，因此投入时间选择一个好的先验肯定有很多用处。
还可以通过更复杂的方式搜索可能的抽样策略。例如，模拟退火 [4] 是一种可用于高维离散空间优化的技术。
我们选择抽样策略时只有一个目标：预测最终结果 \(w\)。另一方面，如果我们想进行更精细的预测，我们可以使用相同的过程，但将 \(\alpha\) 视为目标而不是 \(w\)。
最后，值得注意的是，这与贝叶斯模型选择有关联：具体来说，在不同的模型中，我们可能将 \(w\) 作为二元随机变量来选择不同的子模型。虽然概念上不同，但我们的选举设置和模型选择设置共享许多数学特征，这意味着类似的实验设计方法可以在这两种情况下使用。
参考文献¶
[1] Chaloner, K. and Verdinelli, I., 1995. 贝叶斯实验设计：一篇综述. Statistical Science, pp.273-304.
[2] Foster, A., Jankowiak, M., Bingham, E., Horsfall, P., Teh, Y.W., Rainforth, T. and Goodman, N., 2019. **变分贝叶斯最优实验设计.** Advances in Neural Information Processing Systems 2019 (待出版).
[3] Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A. and Rubin, D.B., 2013. **贝叶斯数据分析.** Chapman and Hall/CRC.
[4] Kirkpatrick, S., Gelatt, C.D. and Vecchi, M.P., 1983. **模拟退火优化.** Science, 220(4598), pp.671-680.