ごはんと飲み物は紙一重

あんまり更新できてないです

俺でも沖縄県知事戦で当選確実を出したい 【沖縄 Advent Calender 2018 1日目】

ついに11月の更新は1回

12月は最低でも3つ書くから勘弁

はじめに

エンジニアにはAdvent Calenderなる文化があり、Adventカレンダーのテーマに則ったブログ記事を書いて投稿していこうというやつです。ちなみに去年もやっていて、

twdlab.hatenablog.com

去年は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)

f:id:ST_ha1cyon:20181130203548p:plain

# サンプル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)

f:id:ST_ha1cyon:20181130203601p:plain

この時点でもおそらくデニーさんが当選しそうな雰囲気がしていますが、次はサンプリング回数を増やしてより、統計的に詰めてみます。

標本平均を見てみる

上記では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)

f:id:ST_ha1cyon:20181130203745p:plain

サンプルサイズ1%, 試行回数10000回の開票でもデニーさんが当選しそうですね。

標本の信頼区間を見てみる

最後の砦である信頼区間を見ることで判断しましょう。信頼度95%(ある人が95%で当選する確証がある)で投票率を計算してみます。rをサンプリング数/全体数とすると以下の式で示せます。

{ r - 1.96 \sqrt{\frac{r(1 - r)}{全体数}} \leq R \leq r + 1.96 \sqrt{\frac{r(1 - 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の更新もお楽しみに。