ある集団の平均値を求めるときに、そのすべてを調べるのは困難な場合は一部のサンプルの平均値を求めて元の集団の平均値を求めることがよく行われます。今回はそのような場合の、一部のサンプルの平均値(=標本平均)から元の集団の平均値(=母平均)の範囲を推定する方法をPythonを用いて説明します。
開発環境
- Python 3.7.7
- scipy 1.4.1
Contents
具体的な状況を考えてみる ― 全社員の身長の平均値は？
ある会社の社員の身長の分布は正規分布に従っているとします。この社員から無作為に10人を抽出し、以下のような身長を得ました。この時、その会社の社員全員の身長の平均値はどの程度の範囲に収まると推定できるでしょうか？
なお、母平均の区間推定とは、「その母平均が〇〇%の確率でその範囲内に入っている」という範囲を求めることであり、その範囲のことを〇〇%信頼区間と呼びます。
母分散が既知の場合の母平均の区間推定
ここでは先ほどの例で、その会社の社員全員の身長の分布の分散(母分散)が33だと分かっているものと仮定します。
数学的背景
正規母集団\(N(μ, σ^2)\)から得た標本\(\{X_1, X_2, …, X_n\}\)の標本平均\( \bar{X}\)は正規分布\(N(μ, \frac{σ^2}{n})\)に従うので、
$$Z = \frac{\bar X – μ}{\sqrt{\frac {σ^2}{n}}}$$
は標準正規分布\(N(0, 1)\)に従います。
これを用いて、信頼水準\(1-α\)で母平均\(μ\)が存在する範囲は
$$- z(\frac{α}{2}) \le \frac{\bar X – μ}{\sqrt{\frac {σ^2}{n}}} \le z(\frac{α}{2})$$
$$\bar{X} – z(\frac{α}{2}) \sqrt{\frac {σ^2}{n}} \le μ \le \bar{X} + z(\frac{α}{2})\sqrt{\frac {σ^2}{n}}$$
と表されます。ただし、\( z(\frac{α}{2}) \)は標準正規分布で確率\( \frac{α}{2} \)を与える数値とします。
以上から、標準正規分布をX軸方向に\( \bar{X}\)だけ平行移動して、Y軸方向に\(\sqrt{\frac {σ^2}{n}}\)倍した分布の〇〇%信頼区間が、母平均が〇〇%の確率で存在している範囲ということになります。
Python プログラム
それではPythonを用いて、信頼水準95%で母平均の区間推定を行ってみましょう。標準正規分布をX軸方向に\( \bar{X}\)だけ平行移動して、Y軸方向に\(\sqrt{\frac {σ^2}{n}}\)倍した分布の95%信頼区間を求めればよいので、scipyで正規分布を表すnormオブジェクトのintervalメソッドを用いて信頼区間を求めてみます。
import math
import pandas as pd
from scipy.stats import norm
df = pd.DataFrame([168.2, 172.3, 165.4, 170.8, 178.9, 173.5, 169.5, 167.1, 176.4, 174.5])
p_var = 33 # 母分散
s_mean = df.mean() # 標本平均
n = len(df) # 標本数
bottom, up = norm.interval(0.95, loc=s_mean, scale=math.sqrt(p_var/n))
print(bottom)
print(up)
[168.09954861] [175.22045139]
よって、母平均μの区間推定の95%信頼区間は
$$168.09954861 \le μ \le 175.22045139$$
となりました。
この計算には母分散の値を必要としますが、現実的には「母平均が分からないのに母分散のみ分かっている」という状況はあまり考えにくいですよね。今回の例でいえば、「その会社の社員全体の身長の平均値は分からないけれど、分散は33と分かっている」ということはどう考えてもあり得ませんよね。
そのため実際には次の「母分散が未知の場合の母平均の区間推定」を使うことがほとんどと思われます。
母分散が未知の場合の母平均の区間推定
数学的背景
正規母集団\(N(μ, σ^2)\)から得た標本\(\{X_1, X_2, …, X_n\}\)の標本平均\( \bar{X}\)と不偏分散\( U^2 \)に対して、
$$T=\frac{\bar X – μ}{\sqrt {\frac{U^2}{n}}}$$
は自由度\(n-1\)のt分布に従うことが知られています。
これを用いて、信頼水準\(1-α\)で母平均\(μ\)が存在する範囲は
$$- t(n-1, \frac{α}{2}) \le \frac{\bar X – μ}{\sqrt {\frac{U^2}{n}}} \le t(n-1, \frac{α}{2})$$
$$\bar X – t(n-1, \frac{α}{2}) \sqrt{\frac{U^2}{n}} \le μ \le \bar X + t(n-1, \frac{α}{2}) \sqrt{\frac{U^2}{n}}$$
となります。ただし、\(t(n-1, \frac{α}{2})\)は、自由度\(n-1\)のt分布で確率\( \frac{α}{2} \)を与える数値とします。
以上から、自由度\(n-1\)のt分布をX軸方向に\( \bar{X}\)だけ平行移動して、Y軸方向に\( \sqrt{\frac{U^2}{n}} \)倍した分布の〇〇%信頼区間が、母平均が〇〇%の確率で存在している範囲ということになります。
Python プログラム
それでは先ほどと同様にPythonを用いて、信頼水準95%で母平均の区間推定を行ってみましょう。自由度\(n-1\)のt分布をX軸方向に\( \bar{X}\)だけ平行移動して、Y軸方向に\( \sqrt{\frac{U^2}{n}} \)倍した分布の95%信頼区間を求めればよいので、scipyでt分布を表すtオブジェクトのintervalメソッドを用いて信頼区間を求めてみます。
以下のように、intervalメソッドの第1引数に信頼水準の0.95を、第2引数に自由度を指定して、第3引数・第4引数にもとのt分布からの変形度合を指定します。
import math
import pandas as pd
from scipy.stats import t
df = pd.DataFrame([168.2, 172.3, 165.4, 170.8, 178.9, 173.5, 169.5, 167.1, 176.4, 174.5])
u_var = df.var() # 不偏分散
s_mean = df.mean() # 標本平均
n = len(df) # 標本数
deg_of_freedom = len(df)-1 # 自由度
bottom, up = t.interval(0.95, deg_of_freedom, loc=s_mean, scale=math.sqrt(u_var/n))
print(bottom)
print(up)
[168.60535249] [174.71464751]
よって、母平均μの区間推定の95%信頼区間は
$$168.60535249 \le μ \le 174.71464751$$
となりました。