Pythonを使った区間推定の問題を解いてみる | 数理モデルと社会科学、時々AI美女

区間推定とは何かを、まずは教科書的な言葉で説明してみます。

例題で学ぶ統計入門よると、

区間推定とは、母数の存在する範囲をデータから構成された区間で示す方法であり、点推定よりも、母数に関するより詳しい情報が得られる。

これでわかったら、はっきりって天才だと思いますね。

なのでもう少しわかりやすく説明してみましょう。

中心極限定理の記事でも説明しましたが、母集団と標本集団は同じではない可能性があります。

しかし、母集団が100万件あったとして、そこからサンプル10件を採取した場合の標本平均と、100件を採取した場合標本平均では、中心極限定理においてサンプル数が多いほど母平均に近づいていきます。

そのため、サンプル数の多い標本をもとに、母集団の特性を推定するのが推測統計といいます。

通常、母集団の特性を推定するには、母集団の平均値 $\mu$ が用いられます。

信頼区間とは

しかし、例えば100件のサンプルの平均が53だとして、それが100万件の母集団の平均といえるでしょうか。

言えないことはないのですが、裏付けのない不確実な数字ですね。

そこで、標本平均に幅を持たせる多ために用いられるのが信頼区間というわけです。

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

信頼区間の計算で、母分散がすでに分かっていることは少ないのですが、その場合に使われるのが以下の式。

$\bar{X} -\frac{\sigma }{\sqrt{n}}Z_{1} < \mu < \bar{X} +\frac{\sigma }{\sqrt{n}}Z_{1}$

この公式に数値を当てはめてもいいのですが、もう少し簡単にプログラムを使って計算してみましょう。

例題は演習で身につける統計学入門から、

殻付き卵の1個当たりの重量の分散は16。本日取れた鶏卵から25個を無作為に抜き取り図った平均は63g。本日取れた鶏卵について信頼水準 $\gamma$ = 0.95の母平均 $\mu$ の信頼区間（g）を求める。

この問題では母分散は16、サンプル数25、標本平均は63gとわかっていますので、Pythonのscipy を使って信頼区間を求めてみましょう。

コードは以下になります。

import numpy as np
from scipy import stats

egg_a = 0.95 # 信頼水準
egg_d = 25-1 # 自由度
egg_m = 63 # 標本平均
egg_v = 16 # 母分散

t_egg = stats.t(loc = egg_m, scale = np.sqrt(egg_v/25), df = egg_d)
bottom, up = t_egg.interval(alpha = egg_a)
print('95% 信頼区間: {:.2f} < x < {:.2f}'. format(bottom, up))

信頼水準95%　信頼区間: 61.35 < x < 64.65

つまり母平均は95％の確率で、61.35g から 64.65g の間に存在するということになります。

ではもう1問解いてみましょう。

同じように、演習で身につける統計学入門からの出題。

ある農園で出荷するリンゴの1個当たりの重量(g)の分散は64と分かっています。本日出荷するリンゴから25個を無作為に取り出し、その重量を計った結果、平均は295でした。本日出荷するリンゴの重要について、信頼水準 $\gamma$ = 0.9で母平均 $\mu$ の信頼区間を求めなさい。

前の問題と同じように、それぞれの数値をプログラムに当てはめていきましょう。

import numpy as np
from scipy import stats

apple_a = 0.9 # 信頼水準
apple_d = 25-1 # 自由度
apple_m = 295 # 標本平均
apple_var = 64 # 母分散

t_dist = stats.t(loc = apple_m, scale = np.sqrt(apple_var/25), df = apple_d)
bottom, up = t_dist.interval(alpha=apple_a)
print('90% 信頼区間: {:.2f} < x < {:.2f}'.format(bottom, up))

信頼水準90% 信頼区間: 292.26 < x < 297.74

この場合、母平均は95％の確率で、292.26ｇと297.74ｇの間と推定できます。

母分散がわからない場合

通常は母分散がわからない場合が多いので、その場合の信頼区間を求める公式を見てみましょう。

$\bar{X}-\frac{ut_{1}}{\sqrt{n}} < \mu < \bar{X} + \frac{ut_{1}}{\sqrt{n}}$

この公式を使って計算すれば、母分散が不明でも標本分散を使って母平均を求めることができますが、それは面倒なので、Python のプログラムに当てはめて計算してみましょう。

例題は以下演習で身につける統計学入門からの出題。

ある高校1年生の男子生徒から31人を無作為に選び、その体重を測定した結果、平均が56、不偏標本分散が25でした。この高校の1年生男子生徒全体の平均体重 $\mu$ の信頼区間を信頼水準 $\gamma$ = 0.95で求めなさい。

この問題でも、プログラムに数値を当てはめて計算してみます。

import numpy as np 
from scipy import stats

std_a = 0.95 # 信頼区間
std_d = 31-1 # 自由度
std_m = 56 # 標本平均
std_v = 25 # 不偏標本分散

t_std = stats.t(loc = std_m, scale = np.sqrt(std_v/31), df = std_d)
bottom, up = t_std.interval(alpha = std_a)
print('95% interval: {:.2f} < x < {:.2f}'.format(bottom, up))
print('信頼水準95% 信頼区間: {:.2f}kg < x < {:.2f}kg'.format(bottom, up))

信頼水準95% 信頼区間: 54.17kg < x < 57.83kg

95％の確率で、1年生全体の体重の平均は54.17kg と 57.83kgの間と推定できます。

母比率の推定

データの種類には「重さや長さ」などだけではなくコインの裏表のように「有か無しか」のようなデータもあります。

前者を計量データといい、後者を計数データといいます。

例えば、ウイルスの感染者数などは、感染しているかしていないかという２者択一の計数データになりますね。

このような場合に使用するのがコインを投げて表が出る確率と裏が出る確率を表すベルヌーイ分布であり、それを複数回繰り返した時の確率分布である二項分布になります。

以下は母比率の信頼区間を求める公式。

$\bar{X}-\sqrt{\frac{\bar{X}(1-\bar{X})}{n}}\cdot Z_{1} < P < \bar{X}+\sqrt{\frac{\bar{X}(1-\bar{X})}{n}}\cdot Z_{1}$

では、実際の例題をもとにPythonを使って母比率を求めてみます。

例題は以下演習で身につける統計学入門からの出題。

A県の県民120人を任意に選び、ウイルスCの検査をした結果、23人が陽性でした。このとき、A県でのウイルス感染率を信頼水準 $\gamma$ ＝ 0.95 で推定しなさい。

Pythonを使って母比率を計算すると、以下のコードになります（参照：たぬハックさん）。

import numpy as np
import scipy as sp

alpha = 0.95 # 信頼水準95%
n = 120 # サンプル数
pos = 23 # 陽性数
p = pos/n # 標本比率

bottom, up = sp.stats.binom.interval(alpha= alpha, n = n, p= p, loc=0)
print('母比率Pの信頼水準95%: {:.3f} < p < {:.3f}'.format(bottom/n, up/n))

母比率Pの信頼水準95%: 0.125 < p < 0.267

つまりウイルスの感染率は、0.125％から0.267％と推定できます。