統計学をpythonにやってもらおう#2

初めに

前回は、正規分布とt分布をpythonで表現する方法を解説しましたが、今回は区間推定について扱っていきます。前回の記事
特に平均を中心に区間推定を行っていきます。

目次
　　１，区間推定とは
　　２，母分散が既知の場合の母平均の区間推定
　　３，母分散が未知の場合の母平均の区間推定
　　４，母比率の区間推定

１，区間推定とは

区間推定とは、『真の値θが区間[L,U]に入る確率が(1-α)以上になるように保証する』という推定方法です。
簡単に言うと、真の値を１点で求める点推定とは反対に、真の値をある区間内(信頼区間)にいるという推定をするということです。
図のpは母比率ですが真の値と考えてもらって大丈夫です。

参考（区間推定について）

定義や理論だと私はわかりずらかったので、具体例を考えていきます。

２，母分散が既知の場合の母平均の区間推定

例えば10人の男子学生の身長データを取ってきた結果、以下の数値のようになったとします。

175,170,169,165,164,179,172,168,175,170

母集団は正規分布に従うとします。この時サンプルの平均は170.7cmでした。
もし、母集団の分散があらかじめ25とわかっている場合、信頼区間95%に設定すると正規分布表の上側の確率が0.025になるパーセント点Z=1.96を求めます。これをpythonで計算していきます。

#母分散が既知の場合の母平均の区間推定
import numpy as np
import scipy.stats as st
import math

sample = np.array([175,170,169,165,164,179,172,168,175,170])
population_stdev = math.sqrt(25.0)  #母分散既知
confidence_factor = 0.95 #信頼区間
q = 1 - (1-confidence_factor)/2   #q=0.975となるように計算
z_critical = st.norm.ppf(q=q)     #0.025での正規分布表の値を求める
sample_mean = sample.mean()
margin_or_error = z_critical*(population_stdev/math.sqrt(len(sample)))
confidence_interval = (sample_mean - margin_or_error,sample_mean + margin_or_error)
print('信頼区間[',confidence_interval[0].round(4),',',confidence_interval[1].round(4),']')

実行結果
信頼区間[ 167.601 , 173.799 ]
つまり、この男子学生の身長の母平均は信頼区間95%の時、167.6~173.8cmの間にいることがわかりました。
コードについては、z_critical = st.norm.ppf(q=q) で上側パーセント点を求め、
z_critical*(population_stdev/math.sqrt(len(sample)))で標準化の計算を行っています。

３，母分散が未知の場合の母平均の区間推定

次に、母分散が未知の場合を先ほどの例で考えてみます。
先ほどの母分散がわかっていないとします。この時、母分散をサンプルの不偏分散で置き換え、t分布を用いて推定します。今回はサンプル数が10個なので自由度(n-1)は9になります。
t分布については前回の記事を参考にしてください。正規分布とt分布

#母分散が未知の場合の母平均の区間推定
import numpy as np
import scipy.stats as st
import math

sample = np.array([175,170,169,165,164,179,172,168,175,170])

confidence_factor = 0.95 
q = 1 - (1-confidence_factor)/2
t_critical = st.t.ppf(q=q,df=len(sample)-1)     #t分布を用いる　自由度df=n-1
sample_mean = sample.mean()
sample_stdev = sample.std(ddof=1)  #母分散未知よりサンプルの不偏分散を用いる
margin_or_error = t_critical*(sample_stdev/math.sqrt(len(sample)))
confidence_interval = (sample_mean - margin_or_error,sample_mean + margin_or_error)
print('信頼区間[',confidence_interval[0].round(4),',',confidence_interval[1].round(4),']')

実行結果
信頼区間[ 167.3608 , 174.0392 ]

近似値として不偏分散を採用していますが、先ほどの既知の場合と比べてそれほどブレがないことがわかります。
なぜかというと、不偏分散 u² を代入することで、新しく出てきた標準化変数 T は、数学的な議論によって標準正規分布にではなく、自由度 n - 1 の t 分布に従うことが証明されているからです。
ただ、標本数が大きい場合は、標準正規分布と t 分布のグラフの形はほぼ一致するため、約 20 以上の標本では、自由度 n - 1 の t 分布、もしくは標準正規分布のどちらを用いても良いです。

４，母比率の区間推定

最後に、母比率の推定についてです。
母比率とは、母集団の中である事象が起こる確率のことで、製品の不良率や番組の視聴率などの様々な場面で必要となります。

推定では、二項分布の枠組みを用います。
二項分布は、一回の試行で二者択一の結果が確立pで生じるベルヌーイ試行をn回行うとき、成功回数Xが従う確率分布B(n,p)であり、今回はpを母比率に対応させます。その後標準化し信頼区間を求めますが本記事では省略します。

具体例を考えていきましょう。
サンプルから製品の不良率を区間推定します。ランダムにとった製品100個から不良品が3個見つかったとき、全体の不良率を(1-α)=0.95、0.90、0.80でそれぞれ区間推定する。この時、p=0.03なのでpythonでは以下のように記述します。

#母比率の区間推定
import numpy as np
import scipy.stats as st
import math

n = 100
defective = 3
phat = defective/n
alpha = [0.025,0.05,0.10]

for i in alpha:
  Z = -norm.ppf(i)  # 上限パーセント点を求める
  confidence_coef = Z * math.sqrt((phat*(1-phat))/n)  # 標準化(信頼係数を求める)
  print('信頼係数=',confidence_coef.round(4))
  print('信頼区間[',(phat-confidence_coef).round(4),',',(phat+confidence_coef).round(4),']')

実行結果
信頼係数= 0.0334
信頼区間[ -0.0034 , 0.0634 ]
信頼係数= 0.0281
信頼区間[ 0.0019 , 0.0581 ]
信頼係数= 0.0219
信頼区間[ 0.0081 , 0.0519 ]

結果から、αの値が大きくなると信頼係数は小さくなり、信頼区間はせまくなることが分かった。
信頼係数95%で推定すると約6%まであり得る結果になったのは意外でしたね。

まとめ

今回は、区間推定の基礎をpythonを通じて学習しました。数学的な証明などの説明は内容が複雑になってしまうのを避けるため、なるべくコードを走らせて何となくやっていることを理解できればなと思い書きました。
あくまで統計学とpythonを学ぶために行っていますが、データが大きくなってこれば実際に活用場面も出てきそう、、