大統領選のシミュレーション(in R)

大統領選に関する記事

沈黙する人工知能 ~なぜAIは米大統領選の予測に使われなかったのか という記事があります。 この記事の中では、アメリカ大統領選挙の性質を調べる為に、全てのアメリカ国民が1/2の確率で、どちらかの候補に投票を行うなら何がおきるのか? というシミュレーションを計算で行っています。(そのソースコードは ここ )

そのシミュレーションの実装では、一人ずつ1/2の確率で投票を行う操作を、全員について実際に乱数を引くことで行っています。 シミュレーションは合計1万回実施し、そこから得られる平均と分散を求めることで結果を出しています。 この「全人口(=3億人)について乱数を取る操作」がとっても大変なようです。 実際、記事では1回のシミュレーションに5秒。最終的な結果を得るのに14時間を費やしたようです……。

1つの州当たりの投票を二項分布で考えてもよい

ところでこの記事の話をそのまま理解すると、人口分の回数の乱数を取る必要はなさそうです。 というのはこれは二項分布で求められる値だからです。

つまり州$state$ごとに、人口$\mathrm{population}(state)$と、選挙人数$\mathrm{elector}(state)$があるとまず見做します。 それから二項分布$\mathbb{B}(\mathrm{population}(state), 1/2)$に従う確率変数 $X_{state}$ を考えます。 $X_{state}$ が人口の半数を超えた州 $X_{state} > \mathrm{population}(state)/2$ は、 その州に割り当てられた選挙人数$\mathrm{elector}(state)$を取得する。……と見做せばいいだろうと思ったわけです。

二項分布で考えて書いたシミュレーション

そこで二項分布を使って実際にシミュレーションを取ってみました。 これは以下のRコードになります。

# ---------- パッケージ読み込み

if( !require(plyr) ) {
  install.packages("plyr")  
}
library(plyr)

# ---------- アメリカ合衆国の州の情報

names <- c("AL","AK","AZ","AR","CA","CO","CT","DE","FL","GA","HI","ID","IL","IN","IA","KS","KY","LA","ME","MD","MA","MI","MN","MS","MO","MT","NE","NV","NH","NJ","NM","NY","NC","ND","OH","OK","OR","PA","RI","SC","SD","TN","TX","UT","VT","VA","WA","WV","WI","WY","DC")
populations <- c(4708708,698473,6595778,2889450,36961664,5024748,3518288,885122,18537969,9829211,1295178,1545801,12910409,6423113,3007856,2818747,4314113,4492076,1318301,5699478,6593587,9969727,5266214,2951996,5987580,974989,1796619,2643085,1324575,8707739,2009671,19541453,9380884,646844,11542645,3687050,3825657,12604767,1053209,4561242,812383,6296254,24782302,2784572,621760,7882590,6664195,1819777,5654774,544270,601723)
electors <- c(9,3,11,6,55,9,7,3,29,16,4,4,20,11,6,6,8,8,4,10,11,16,10,6,10,3,5,6,4,14,5,29,15,3,18,7,7,20,4,9,3,11,38,6,3,13,12,5,10,3,3)

# ---------- シミュレーション回数
N <- 10000

# ---------- シミュレーション実装

states <- data.frame(name = names, population = populations, elector = electors)

transformed <- ddply(
  states,
  .(name),
  transform,
  counter=1:N,
  elected=ifelse(
    rbinom(N, population, 1/2) > (population/2),
    elector, 0
  )
)
results <- ddply(
  transformed,
  .(counter),
  summarize,
  result=sum(elected)
)

print(mean(results$result))
print(var(results$result))
print(sd(results$result))

動かすと割と即座に結果が出て、平均 $\mu = 268.939$ 標準偏差 $\sigma = 48.90947$ でした。

江端家の電気代を安く済ませるには、二項分布が有効だったわけですね。

より多くのシミュレーション回数でも比較してみる

ところでこのシミュレーションって、本当に1万回で十分なんでしょうか。

ということで、せっかく軽い実装も作ったわけなので、シミュレーション回数をもっと増加させてみました。 そして標準偏差がどう変わるのか見てみました。

以下が結果です。

シミュレーション回数 平均 標準偏差
10000 268.939 48.90947
20000 269.115 50.83685
30000 269.112 52.72271
40000 268.778 51.43323
50000 267.167 48.54549
60000 267.827 51.43179
70000 267.822 52.07724
80000 269.005 51.09182
90000 266.911 51.9244
100000 268.782 50.56805
以下、飛んで - -
1000000 271.18 52.39613

(本当はシミュレーション回数ごとに、また大きな回数シミュレーションを行って、その平均的な値を見たほうがいいのかもしれないが。) この結果を見る限りは、シミュレーション回数を増やしたとしても、特にこれ以上分散が大きく変化することは無さそうですね。

ちなみに100000回、シミュレーションを回した時のヒストグラムは以下のようになりました。階級幅は1です。

100000回シミュレーションのヒストグラム

ところで元記事に関する私なりの考察。

も、最初書こうと思ってたんですけれど、止めておきます。 なんというか政治の話は難しいです。

その他のお話

今回書いたRコード。正直、あまりきれいに書けた気がしなくて、頭痛いです。 Rのもつ基礎データ構造の理解がどうもイマイチらしくて、 それぞれの型( data.framelist )がどう振る舞うのか理解できてない気がするんですよね……。 plyr もそれを持ち出さなくていい方法があればそっちを使いたいし、 R言語徹底解説を積読してるから、あれを読みたいですね……。