区間推定2 選挙速報の当選確実の出し方

#選挙速報の当選確実出し方
テレビの選挙報道番組では、開票率0%や1%で当選確実が出る事がある。これは事前のアンケート結果や当日の出口調査、開票情報、選挙中のニュースによる票の流れの読みなど、テレビ局各社が様々なノウハウと技術を駆使して当選確実を出している。
テレビ局としては、他社を出し抜いて早く当選確実を報道したいが、証拠不十分で当選確実を報道すれば、別候補者が当選してしまう事もある。実際に過去の報道では、当選確実の報道を覆すケースが発生している。
では報道機関はどの様な計算をして当選確実を出しているのだろうか。今回は出口調査だけを使用して、母集団の得票率を推定してみよう。

#比率の区間推定：選挙速報の場合
都知事の選挙を例に考える。

##比率の区間推定の背景・条件
テレビ局ディレクターである私は、都知事選でのA氏の得票率に興味がある。投票所の出口調査の結果から母集団の得票率が50%を超える事がわかれば、番組開始と同時に当選確実を報道しようと考えている。
出口調査は、ランダムに時刻と場所を変えて合計100人に対してを行い、60人がA氏に投票している事が分かっている。

##比率の区間推定の流れ
私は東京都民の有権者全員の投票を知り得ない。手元にある情報は出口調査100人分の投票結果だけである。
A氏の真の得票率を知る為に、これから100人分の投票結果の得票率(=平均値)を元に、中心極限定理を利用して、標本の得票率(=平均値)の標準正規分布を導く。標本の得票率(=平均値)の標準正規分布が分かれば、区間推定が可能になる。

##比率の区間推定の計算
■標本比率 $\hat{p}$ と母比率 $p$
今回の100人の出口調査（標本抽出）から得たA氏の得票率（標本比率）$\hat{p}$は、$\hat{p}=60/100=0.6$　である。しかし、出口調査の対象とした100人の標本が変われば、$\hat{p}$も変わるので、$\hat{p}$は確率変数と考える。

出口調査のi番目の人がA氏に投票したと回答している場合は1、それ以外は0を取る確率変数を$X_i　(i=1,2,...,n)$とすると、次の式の確率変数$X$は、出口調査でA氏に投票した人数（得票数）を表す。
　　$X=X_1+X_2+...+X_n$　　($X_i$は0または1)
ちなみに$X$は確率$\hat{p}$で試行回数n=100の二項分布となる。
ところで $\hat{p}=X/n$ なので、$\hat{p}$の期待値は次の式になる。
　　$E[\hat{p}] = E[X/n] = \frac{1}{n}E[X_1+X_2+...+X_n] = \frac{1}{n}(E[X_1]+E[X_2]+...+E[X_n])$

標本$X_i$は無作為に選ばれたなら$E[X_i]=p$、すなわち$X_i$の期待値は母集団の得票率$p$になるので、次の式となる。
　　$E[\hat{p}]=\frac{1}{n}(np)=p$

この式から、標本数の大小にかかわらず、得票率$\hat{p}$の期待値は母集団の得票率$p$と同値になる事が分かる。今回の例に限らず、標本比率の期待値は母集団比率に一致する。
　

■標本分散$V[\hat{p}]$　（Z値の計算方法1を採用する場合に使用する）
出口調査の得票数$X=n\hat{p}$なので得票数の分散は$V[X] = V[n\hat{p}]= n^2V[\hat{p}]$であり、また[二項分布][2]の公式より$V[X] = np(1-p)$なので、
　　$V[\hat{p}] = p(1-p)/n$

■標準誤差 $SE$　（Z値の計算方法2を採用する場合に使用する）
まず母分散を計算しよう。確率変数がXだと上に記載した確率変数と混乱しそうなので、ここではYで考えている。

分散の定義は次の通りである。
　$V[Y]=\displaystyle\sum_{i} (y_i-\mu)^2 p_i$
比率とは確率変数Xが0または1の2種類を取る場合の平均値であるから、比率の分散を定義すると次の様になる。
　$V[Y]=\sum_{y=0}^1 (y-\mu)^2 p(y)　　(Y=\{0,1\})$
ここで、yの期待値であるμは、1となる確率を表しているので、p(y=1)に等しい。p(y=1)=p=μとおくと、p(y=0)=1-pとなるので、分散は次の様になる。
　$V[Y]=(0-p)^2(1-p) + (1-p)^2(p) = p(1-p)$

上の様に、確率pで1を、確率1−pで0をとる離散確率分布は[ベルヌーイ分布][2]と呼ぶ。
ちなみに、二項分布の一回だけの試行の分布がベルヌーイ分布となる。

標準誤差SEは標本平均の標準偏差であるから、中心極限定理より次の式が成り立つ。
　$SE=\sqrt{σ^2/n}$
今回の得票の例では、標準誤差算出の為には$σ^2$に[ベルヌーイ分布][2]の分散　$V[Y]=p(1-p)$を入れる。
　$SE=\sqrt{V[Y]/n}=\sqrt{p(1-p)/n} 　　(Y=\{0,1\})$

■標準正規分布に従うZ値
Z値の求め方を3種類示す。どれかを採用すれば区間推定可能である。

【Zの計算方法１】
標本比率$\hat{p}$は$\hat{p}=(X_1+X_2+...+X_n)/n$と、0か1を取る[ベルヌーイ分布][2]に従う確率変数$Xi$の和の形をしているので、nがある程度大きければ$\hat{p}$の分布は中心極限定理によって正規分布に近似できる。そこで$\hat{p}$から期待値を引き、標本標準偏差で割る事で、標準正規分布に従うZ値に変換できる。
　$Z=(\hat{p}-E[\hat{p}])/\sqrt{V[\hat{p}]}\\　=(\hat{p}-p)/\sqrt{p(1-p)/n}$

【Zの計算方法２】
標本比率$\hat{p}$は$\hat{p}=(X_1+X_2+...+X_n)/n$と、0か1を取る[ベルヌーイ分布][2]に従う確率変数$Xi$の和の形をしているので、nがある程度大きければ$\hat{p}$の分布は中心極限定理によって正規分布に近似できる。そこで$\hat{p}$から期待値を引き、標準誤差で割る事で、標準正規分布に従うZ値に変換できる。
　$Z=(\hat{p}-E[\hat{p}])/SE\\　=(\hat{p}-p)/\sqrt{p(1-p)/n}$

【Zの計算方法３】
nがある程度大きい時は、中心極限定理によって二項分布B(n,p)は正規分布N(np,np(1-p)) に近似できる。よって、Xが二項分布に従う場合、Xを標準化したZ値もnが十分に大きければZは標準正規分布N(0,1)に従う。
　$Z=(X-np)/\sqrt{np(1-p)}$
標本比率$\hat{p}$ は、$\hat{p}=X/n$（出口調査の得票数を調査数で割ったもの）なので、上の式の分母と分子をnで割り、$X/n$を$\hat{p}$とおくと、次の式となる。
　$Z=(X/n-np/n)/\sqrt{np(1-p)/n^2}\
　=(\hat{p}-p)/\sqrt{p(1-p)/n}$
　

■Z値の計算簡略化
上記Z値の式の分母に区間推定したい母集団の得票率$p$が残っていると、その後の計算が手間となる。そこで、nが大きければ標本から推定した$\hat{p}$は母集団の$p$にほぼ一致すると考えられるので($\hat{p}$は$p$の一致推定量)、分母の$p$を$\hat{p}$で置き換える。
　$Z=(\hat{p}-p)/\sqrt{\hat{p}(1-\hat{p})/n}$

※置き換え無しの厳密な計算は後で説明する。

■標準正規分布を使用した区間推定
今回は信頼係数95%の信頼区間を求めてみよう。標準正規分布表を参照し、信頼度95%(有意水準0.05)のZ値の臨界値を求めると1.96なので、
　$\begin{align}
　P[-1.96<Z<1.96] &=0.95\
　P[-1.96<Z<1.96]
&=P[-1.96<(\hat{p}-p)/\sqrt{\hat{p}(1-\hat{p})/n}<1.96　]\
&=P[　\hat{p}-1.96\sqrt{\hat{p}(1-\hat{p})/n} < p < \hat{p}+1.96\sqrt{\hat{p}(1-\hat{p})/n}　]
\end{align}$
得票率の下限は次の値となる。
　$\hat{p}-1.96\sqrt{\hat{p}(1-\hat{p})/n} = 0.6 -1.96\sqrt{0.6(1-0.6)/100}=0.6-0.096=0.504$

■結論
95%の信用区間の下限が得票率50%を超えているので、A氏の当選は濃厚と言える。もっとも、実際の報道では、間違いがあってはメディアの信用問題になるので、信用区間は99%や99.9%を使い、さらに他の情報も総動員して当選確実を出しているはずである。
　

##比率の区間推定の厳密な計算
先ほどは、次の式の分母の$p$を$\hat{p}$で置き換えて計算した。
　$Z=(\hat{p}-p)/\sqrt{p(1-p)/n}$
厳密には置き換え無しで計算するべきであろうが、どの程度結果に差がでるか見てみよう。

　$-1.96 <(\hat{p}-p)/\sqrt{p(1-p)/n}<1.96$
この不等式は次の様に書き換え可能である。
　$\begin{align}
1.96\sqrt{p(1-p)/n} &>｜(\hat{p}-p)｜ \
1.96^2p(1-p)/n &>\hat{p}^2 - 2\hat{p}p + p^2\
1.96^2p -1.96^2p^2 &>n\hat{p}^2 - 2n\hat{p}p + np^2\
(1.96^2 +n)p^2 -(1.96^2 +2n\hat{p})p +n\hat{p}^2 &> 0 \
(1.96^2/2n+1/2)p^2 -(1.96^2/2n +\hat{p})p +\hat{p}^2/2 &> 0 \
\end{align}$

pの二次方程式として解くと
　$\begin{align}
p &= \frac{1.96^2/2n + \hat{p} \pm \sqrt{(1.96^2/2n + \hat{p})^2-(1.96^2/n +1) \hat{p}^2} }{(1.96^2/n +1)}\
&=\frac{0.6192 \pm \sqrt{(0.6192)^2-(1.0384)0.36} }{1.0384}\
&=\frac{0.6192 \pm 0.0981}{1.0384}\
&=0.5018,　0.6908\
\end{align}$

∴0.5018 <p<0.6908

この計算は、正規分布への近似誤差は残るものの、$p$を$\hat{p}$で近似した誤差は無いので、より厳密な推定と言えよう。しかし結果の信頼区間を見ると、先の$p$を$\hat{p}$で近似した推定結果とほぼ同じであった。ある程度の標本数を揃えるのであれば、この厳密な計算を行う必要性は無いだろう。

ここでは得票率60%のランダムな投票データを作成し、区間推定の結果が妥当か検証してみよう。

区間推定の精度

import numpy as np
import matplotlib.pyplot as plt
from   scipy import stats

loop_count = 10000
n = 100 #サンプルサイズ
#ppf (Percent point function)
score = stats.norm.ppf( 0.975) #片側0.975信頼区間(両側0.95と同様)の値
print("Z値{:.3f}".format(score))
np.random.seed(0)#同じ疑似乱数を生成する為、シード値を指定。



#母集団の配列ndarrayを作成　
arr_a=np.array([])
arr_a=np.hstack((arr_a,np.full(int(4),0)))#A氏以外に投票
arr_a=np.hstack((arr_a,np.full(int(6),1)))#A氏に投票

#母集団の平均
mu_population =np.mean(arr_a)

p_h        =np.array( [] )
p_l        =np.array( [] )
p_h_strict =np.array( [] )
p_l_strict =np.array( [] )
for i in range(0,loop_count):
    #母集団からのランダムにサンプル取得
    sample = np.random.choice(arr_a, n, replace=True)#replace=False抽出するデータに重複無し
    mu = np.mean(sample) #標本比率pはX(0か1)の標本平均mu
    #std= np.std(sample, ddof = 1)
    std= (mu*(1-mu))**0.5 #標本平均(=標本比率)の標準偏差
    se = std / len(sample)**0.5

    p_h = np.append( p_h, mu + score * se )
    p_l = np.append( p_l, mu - score * se )

    p_h_strict = np.append( p_h_strict, \
                ( score**2 /(2*n) +mu \
                  +( (score**2 /(2*n) +mu)**2  -  (score**2/n +1)*mu**2    )**0.5\
                )\
                / (score**2 /n + 1)\
            )
    p_l_strict = np.append( p_l_strict, \
                ( score**2 /(2*n) +mu \
                  -( (score**2 /(2*n) +mu)**2  -  (score**2/n +1)*mu**2    )**0.5\
                )\
                / (score**2 /n + 1)\
            )

print("信頼区間が母平均を含む割合(一般的な方法):")
print(1 - (len(p_h[p_h < mu_population]) \
          +len(p_l[p_l > mu_population]) ) /loop_count)

print("信頼区間が母平均を含む割合(厳密な方法)）:")
print(1 - (len(p_h_strict[p_h_strict < mu_population]) \
          +len(p_l_strict[p_l_strict > mu_population]) )/loop_count)

結果

Z値1.960
信頼区間が母平均を含む割合:
0.9497
信頼区間が母平均を含む割合(厳密）:
0.9497

一般的な方法による区間推定も、厳密な方法による区間推定も、算出した信頼区間が母比率（つまり母平均）を含む割合は、ほぼ95%となった。

[区間推定1][1]
[1]:https://qiita.com/keymc021/items/4e969f85268c6498bc34
[2]:https://qiita.com/keymc021/items/301378111dd47128385b