ついに11月の更新は1回
12月は最低でも3つ書くから勘弁
はじめに
エンジニアにはAdvent Calenderなる文化があり、Adventカレンダーのテーマに則ったブログ記事を書いて投稿していこうというやつです。ちなみに去年もやっていて、
去年はPythonのアドベントカレンダーのトップバッターを掴んでました。まぁ内容は中々怖いですが。
今年は一応
- 沖縄 Advent Calender
- Python Advent Calender
- 琉大 Advent Calender
の3つをやろうと思っているのでその第1弾になる沖縄AdventCalenderをやっていこうというわけです。
ほんへ
今年は突然の翁長知事の急逝によって、急遽沖縄で知事選があり玉城デニー沖縄新知事が誕生いたしました。メディアはこぞって開票から10分足らずで当選確実を出していたのですが、「ほんとか?」と疑問に思っていたわけです。
ちょうど統計学の勉強もしていたので、自分でも当選確実を出してみようと思い立った感じです。というわけでお得意のPythonつかって俺も当選確実を出します!
モジュールのインポート
今回は以下のモジュール使います!
import pandas as pd from scipy import stats import numpy as np import matplotlib.pyplot as plt %matplotlib inline plt.style.use("ggplot") import seaborn as sns from tqdm import tqdm_notebook as tqdm from bs4 import BeautifulSoup from urllib.parse import urljoin from urllib.request import urlopen
県知事選のデータをクローリング
正確な県知事選のデータを使いたいのでさっとクローリングします。
url = "https://www.nhk.or.jp/senkyo/database/local/2018/0930_okinawa/"
governer_page = urlopen(url).read()
soup = BeautifulSoup(governer_page, "lxml")
計算に必要な情報をスクレイピング
投票数
voters, voter_turnout = soup.find(class_="data").find_all("dd")[1:3]
voters = int(voters.text[:-1].replace(",", "")) voter_turnout = float(voter_turnout.text.split("%")[0]) votes_all = round(voters*voter_turnout/100)
print("有権者数: {}".format(voters)) print("投票率: {}".format(voter_turnout)) print("投票数: {}".format(votes_all))
有権者数: 1146815 投票率: 63.24 投票数: 725246
各候補者の投票数
candidates = soup.find_all(class_="profile-data") vote_denny, vote_sakima, vote_kaneshima, vote_hatsumi = \ [int(_.text.split("\n")[0].replace(",", "")) for _ in candidates] vote_all_ok = vote_denny + vote_sakima + vote_kaneshima + vote_hatsumi
print("玉城デニー: {}人".format(vote_denny)) print("佐喜真淳: {}人".format(vote_sakima)) print("兼島俊: {}人".format(vote_kaneshima)) print("渡口初美: {}人".format(vote_hatsumi)) print("総投票数 {}人".format(vote_all_ok))
玉城デニー: 396632人 佐喜真淳: 316458人 兼島俊: 3638人 渡口初美: 3482人 総投票数 720210人
# おそらく無効票 invalid_votes = round(voters*voter_turnout/100) - vote_all_ok invalid_votes
5036 票あるみたいです。
データの作成
votes_data = [] for i, votes in enumerate([invalid_votes, vote_denny, vote_sakima, vote_kaneshima, vote_hatsumi]): votes_data += [i+1 for _ in range(votes)] votes_data = np.random.permutation(np.array(votes_data))
統計的に誰が当選するのか見てみる
ここからが勉強の成果の見せ所で3つのアプローチで見ていきたいと思います。
前提知識
- 標本平均
- 母集団(ここでは投票結果)から任意のサイズを指定して取り出したn回のサンプリングの平均
- つまりサンプリング回数が多い・サイズが大きいほど母平均(母集団の平均)に近づく
# ヒストグラム可視化用関数 def plot_votes(x, y): plt.barh(x, y) plt.title("Okinawa governer selection in 2018") plt.xlabel("frequency") plt.ylabel("votes").set_rotation(0) plt.yticks([_ for _ in range(1,6)], ["invalid", "denny", "sakima", "kaneshima", "hatsumi"])
少ない値で開票してみる
# 標本(サンプルサイズ100)で取り出す np.random.seed(98) sample_100 = pd.Series(np.random.choice(votes_data, 100)) plot_votes(sample_100.value_counts().index, sample_100.value_counts().values)
# サンプル1000で取り出してみる np.random.seed(98) sample_100 = pd.Series(np.random.choice(votes_data, 1000)) plot_votes(sample_100.value_counts().index, sample_100.value_counts().values)
この時点でもおそらくデニーさんが当選しそうな雰囲気がしていますが、次はサンプリング回数を増やしてより、統計的に詰めてみます。
標本平均を見てみる
上記では1回の試行なので、より多くサンプリングした平均(=標本平均)をとってみます。
sample_num = round(votes_all * 0.01) sample_num
サンプリング数は全体の1%である7252枚にします。
np.random.seed(98) sample_means = [] for i in tqdm(range(10000)): sample = pd.Series(np.random.choice(votes_data, sample_num)).value_counts().sort_index() sample_means.append(sample.values.tolist())
なうサンプリング(10000回)
# 標本のデータフレームをつくる sample_means_df = pd.DataFrame(sample_means, columns=["invalid", "denny", "sakima", "kaneshima", "hatsumi"]) sample_means_df.head()
invalid | denny | sakima | kaneshima | hatsumi | |
---|---|---|---|---|---|
0 | 43 | 3957 | 3162 | 45 | 45 |
1 | 52 | 4034 | 3100 | 36 | 30 |
2 | 45 | 3947 | 3176 | 37 | 47 |
3 | 43 | 4029 | 3108 | 41 | 31 |
4 | 54 | 3928 | 3210 | 33 | 27 |
10000回行ったサンプリングで基礎統計量を見てみます。
sample_means_df.describe()
invalid | denny | sakima | kaneshima | hatsumi | |
---|---|---|---|---|---|
count | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 |
mean | 50.461700 | 3965.033800 | 3165.353500 | 36.329600 | 34.821400 |
std | 7.071487 | 42.009493 | 41.928331 | 6.079643 | 5.914508 |
min | 26.000000 | 3819.000000 | 3007.000000 | 15.000000 | 15.000000 |
25% | 46.000000 | 3937.000000 | 3136.000000 | 32.000000 | 31.000000 |
50% | 50.000000 | 3965.000000 | 3165.000000 | 36.000000 | 35.000000 |
75% | 55.000000 | 3994.000000 | 3194.000000 | 40.000000 | 39.000000 |
max | 81.000000 | 4122.000000 | 3313.000000 | 62.000000 | 65.000000 |
# 7252票のサンプリングを10000回やった平均 plot_votes([_+1 for _ in range(5)], sample_means_df.mean().values)
サンプルサイズ1%, 試行回数10000回の開票でもデニーさんが当選しそうですね。
標本の信頼区間を見てみる
最後の砦である信頼区間を見ることで判断しましょう。信頼度95%(ある人が95%で当選する確証がある)で投票率を計算してみます。rをサンプリング数/全体数とすると以下の式で示せます。
def conf_interval_left(r, sample_num): return r - 1.96 * np.sqrt(r * (1-r) / sample_num) def conf_interval_right(r, sample_num): return r + 1.96 * np.sqrt(r * (1-r) / sample_num)
つまりデニーさんの投票数を見ると
sample_means_df.denny.mean()
平均が 3965.0338 みたいですね。
denny_r = sample_means_df.denny.mean() / sample_num
# 左辺
conf_interval_left(denny_r, sample_num)
0.5352928708428588
# 右辺
conf_interval_right(denny_r, sample_num)
0.5582079013579135
つまり上記から1%の開票時点でデニーさんは全体のうちの約53%から約56%の支持率を95%の確率で得ていることになるので、この時点で統計的には過半数を超える投票率のため当選確実を出すことができます。
ちなみに佐喜真さんの投票率を見ると、
sample_means_df.sakima.mean()
平均は 3165.3535 票ですね。
sakima_r = sample_means_df.sakima.mean() / sample_num
# 左辺
conf_interval_left(sakima_r, sample_num)
0.42506538703700086
# 右辺
conf_interval_right(sakima_r, sample_num)
0.44789476188743377
95%の確率で42%から44%なので、ギリギリ届かずという感じになります。
という感じで簡単な統計でしたが、1%わかるだけでもほぼ確実にデニーさんが当選するということがわかるみたいです。統計素晴らしい。とはいえ、既に全体の票数がわかっているから卑怯なことをしているのでちょっと良くない。
沖縄を題材にしたよい記事がかけたんじゃないかということで、次回のAdvent Calenderの更新もお楽しみに。