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

前回の記事を、 元の記事の著者である江端さんにお伝えしたところ ご好評をいただきました。

そのうえで、質問として件のRコードを「50.01%のバイアス」がある場合へと拡張するにはどうしたらよいのか? ……というご質問をいただきました。(江端さんは普段はR言語以外の言語で、こういうことをなされているみたいです。)

この江端さんの質問へ回答を行いたいです。

そこでこの記事では、サンプルコードの拡張と、そもそもどの様な仕組みで動いているのかをお話いたします。

二項分布と、今回のシミュレーションの対応関係

そもそも二項分布とは?

二項分布とは、 「コイントスを何度か行ったときに、コインが表になる回数は、回数ごとに何%の確率で生じるのか」 についての分布です。

またコイントス1回分について「出る確率が表にどれだけ偏っているのか」も踏まえて計算することが出来ます。

この端的な言い方では、ちょっと分かりづらいと思います。なので少し具体的に言います。

今、10枚のコインを考えます。

〇〇〇〇〇〇〇〇〇〇

この10枚のコインは、表裏があって……

〇:オモテ
●:ウラ

コイントスをすると、このどちらかになります。

さて、今1回のコイントスで、オモテが出る確率が1/2(ということはウラが出る確率もまた、1/2ですね)としましょう。

この状況で、実際にコイントスを10回すると、以下のような順番でコインが出ました

〇●〇●●〇〇●●〇

このコイントスの結果はオモテが5回出ていますね。

もう一度、実際にコイントスを10回投げてみましょう。すると以下のような順番でコインがでました。

●〇〇●〇〇●●〇〇

今回はオモテが6回出ていますね。

さて。二項分布とは、以下のような のことを言います。(これは概算なので適当ですが。)

オモテになったコインの枚数 そうなる確率
0枚 約0%
1枚 約1%
2枚 約4%
3枚 約12%
4枚 約21%
5枚 約25%
6枚 約21%
7枚 約12%
8枚 約3%
9枚 約1%
10枚 約0%

この表に従えば、コインを10回投げて、3回「ピッタリ」出す確率は、約12%程度だと分かります。 では5回「以上」出せる確率は……? それは5~10枚の確率の合計になるので「25% + 21% + 12% + 3% + 1% + 0%」で、「62%」となります。

一般に、二項分布を表すために

$\mathbb{B}($ <何回連続でコインを投げるのか>, <表が何%の確率で出るコインなのか> $)$

という表記をします。例えば、上記の表は $\mathbb{B}(10, 1/2)$ という表(二項分布)です。

また重要な点として、二項分布は、「何番目が、何%なのか」を、簡単に求める式が知られています。 なので実際には表を書かなくても、その式で、その場所の値を計算することができます。

今回の有権者との対応関係

今回、各州の各有権者は、まさにコイントスのように 1/2(=50%) で投票を行うのでした。 もともとの考え方は、3億回のコイントスを実施する(州ごとにその人口分のコイントスを実施する)ようなものでした。

ですが各州においての投票結果は、二項分布の定義と合致した挙動を示すはずです。 したがって二項分布に従う乱数というものを作ることができれば、 3億回のコイントスを実施しなくても、まったく同じ結果を得られるはずです。

さて二項分布に従う確率変数(乱数) というのをシミュレーションすることを考えます。 その場合、結果を出力するためには、以下のような二項分布の表を作って、一番右の「ここまでの合計確率」の値を使います。

ちなみにこの表の一番右の「ここまでの合計確率」の表を「累積分布」と言います。 (累積分布も何番目が何%なのか計算する式が知られてい…たはず。/ すいません。これ確認できないでいます。 →以下の追記へ)

(数学に詳しい方へ向けた追記: この累積分布は「正則不完全ベータ関数」を使えば手に入ることを 内場先生 から教えていただけました。ただそうすると途中に積分を要するのですから、計算コストや計算精度としては微妙かもしれないです。実際には幾何分布に従う乱数列との間のうまい関係(リンク先PDF/対象は23~24ページ)を使って求めますが、ここではオーダが減らせることを中心に解説をしたいので、累積分布関数から逆関数法を解説します。これはもちろん微妙な面もありますが……。)

(さらに追記: 実際にRが使っているのはCatherine Loaderさんが2000年に公開した「Fast and Accurate Computation of Binomial Probabilities」という論文のアルゴリズムで、これは本当にめちゃくちゃ早いです……。この辺の技術は世の中、極まってますよね……。ちなみに↑のPDFのアルゴリズムだと、そもそもnの増加に対してはまだ遅くてダメでした。あやふやな内容を書いちゃダメですね。)

オモテになったコインの枚数 そうなる確率 ここまでの合計確率
0枚 約0% 約0%
1枚 約1% 約1%
2枚 約4% 約5%
3枚 約12% 約17%
4枚 約21% 約38%
5枚 約25% 約63%
6枚 約21% 約84%
7枚 約12% 約96%
8枚 約3% 約99%
9枚 約1% 約100%
10枚 約0% 約100%

まず0~1までの浮動小数点数値の乱数を振ります。 その後その乱数値が、ひとつ前の「ここまでの合計確率」から、「ここまでの合計確率」との間に入るような位置を、検索します。 そうすると、コイントスが表になる回数を求める処理は、検索アルゴリズムのオーダ分のコストで完了します。

この時、二項分布の表はソート済みと見做せます。 したがってバイナリサーチが使えて、結果を $\mathcal{O}(\lg n)$ で取得することができます。

実際今回の場合は、 $\lg(300000000) \leq 29$ です。

つまり、たった29回の検索操作で、3億回のコイントス相当の処理を実施できます。 (さらにもっと効率が良い二項分布の検索方法があるかもしれませんが、少なくともこれより早いハズです。)

Rにおける二項分布の実装と、前回のコード解説

さてRでは二項分布に従う乱数取得を行うために、以下の関数呼び出しを使います。

rbinom( "<取得したい乱数の個数>", "<何回のコイントスを行うか>", "<表が何%の確率で出るコインなのか>" )

この結果は配列(Rではベクトルといいます)で得られます。

R言語は配列(ベクトル)が基礎データになっていて、 配列に対して適当な計算を書くと、配列の中身を一括で処理する仕組みになっています。

たとえば

"<コイントスが表になった回数の配列>" > "<コイントス回数の半分>"

という比較処理を書くと、これは「行ったコイントスのうち、半分以上が表だったかどうかを真偽値で持った配列」が得られます。

そこで改めて元のコードを眺めてみると以下のようになっています

rbinom(N, population, 1/2) > (population/2)

つまり「1/2の確率で表になるコインについて、州の人口回のコイントスを行う」時に、 オモテは何枚が出るのか……という乱数値を「N個」求めて、その配列を取得しています。

この配列の1つ1つの要素は「0~population(州の人口)」の値を取ります。 その値が(population/2)よりも大きいなら真、そうでないなら偽が入る配列が最終的に得られます。

前回のサンプルでは、このような仕組みで得た配列の値に応じて、選挙人数を足していく操作を行っていました。

今回ほしいもの(バイアス付きの結果を求める方法)

さて今回欲しい物を得るための結論としては……

# このような形で、1/2の部分を50.01% (= 0.5001)にすればよいです。
rbinom(N, population, 0.5001) > (population/2)
# あるいはこれは、以下のように書いても構いません(0.0001 = 0.01%のバイアス)
rbinom(N, population, 1/2 + 0.0001) > (population/2)

とするとよいでしょう。 単にコインのオモテが出る確率が、まさに求めている党の支持率に他ならないためです。

ちなみに元の制御の仕方に合わせて、

( rbinom(N, population - (population/10000), 1/2) + (population/10000) ) > (population/2)

という方式も可能です。分布の絵面は、こっちのほうが江端さんの元のプロットに近いものになります。 この形だと確定で、州の人口の1/10000人が投票するという差異があるためです。

今回の改変版の実際のサンプル

ということで、これらの結論を取り込んで、さらにプロットまで自動化したりしたサンプルコードが以下のものになります。 確率の偏り方は P という変数にしましたので、この値をいじって動かせば、もっといろいろなシミュレーションを行えます。 今回も結果はすぐに出ます。(ただし1回目だけ ggplot2reshape2 といったパッケージのインストールには時間がかかります)

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

# 各種data.frameの操作用(TODO: dplyrに乗り換え)
if( !require(plyr) ) {
  install.packages("plyr")
}
library(plyr)

# melt関数用
if( !require(reshape2) ) {
  install.packages("reshape2")
}
library(reshape2)

# プロット処理用
if( !require(ggplot2) ) {
  install.packages("ggplot2")
}
library(ggplot2)

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

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

# 得票率の偏り
P <- 1/2 + 0.0001

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

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

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

# 投票人の数からdemocraticへ投票した人を引けばrepublicへの投票数になる
elector_count <- sum(electors)

democratics <- results$result
republics   <- elector_count - democratics

# ----------- 統計情報出力

# democraticへの投票の平均・分散・標準偏差
print(mean(democratics))
print(var(democratics))
print(sd(democratics))

# republicへの投票の平均・分散・標準偏差
print(mean(republics))
print(var(republics))
print(sd(republics))

# 互いの勝率
democratics_wins = (republics < democratics)
print(sum(ifelse(democratics_wins, 1, 0)))
print(sum(ifelse(democratics_wins, 1, 0)) / N)

republics_wins = (democratics < republics)
print(sum(ifelse(republics_wins, 1, 0)))
print(sum(ifelse(republics_wins, 1, 0)) / N)

# ----------- プロット

df <- data.frame(
  democratic = democratics,
  republic = republics
)
graph_df <- melt(df)

g <- ggplot(
  graph_df,
  aes(
    x = value,
    fill = variable
  )
)
g <- g + geom_histogram(
  alpha = 0.5,
  position = "identity",
  binwidth=1
)
g <- g + scale_x_continuous(
  limits = c(0, elector_count)
)
plot(g)

実際にP値を50.01%に動かしたところ以下の結果を得ました。 プロットの階級幅は(前回と同じく)1です。

P値+0.0001

変数名
democraticsの平均選挙人獲得数 391.895
democraticsの分散 1625.285
democraticsの標準偏差 40.31482
republicsの平均選挙人獲得数 146.105
republicsの分散 1625.285
republicsの標準偏差 40.31482
democraticsの勝利回数 9958
democraticsの勝率 0.9958
republicsの勝利回数 42
republicsの勝率 0.0042

この結果は、江端さんの結果と大きな差はないため、シミュレーションとして問題ないといえるでしょう。

まとめ

追記

この記事の内容通りのアルゴリズムで同様のシミュレータをC++でも書いてみた。 (つまり言ってる内容で、Rのより素早い方法を使わなくても実際に高速化できることを実証した。)

ソースコード(cpp)

ちなみに最初Cの範囲でやろうと思ってたので、いまいちC++ではないコードになってしまった。 C++でやるならCDFとかは std::shared_ptr<std::vector<double>> であるべきだとか、 まあそもそもクラスを作るべきだとかあるが、てかコードきたねえなおい。とかいろいろあるが、サンプルなのでご勘弁を。