区間推定#
区間推定（interval estimation）は推定したいパラメータ
の
なお、この
正規母集団の母平均の区間推定#
確率変数
これを標準化して
とすると、これは平均0、標準偏差1の標準正規分布に従う。
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
from scipy.stats import norm
z = np.linspace(-4, 4, 100)
y = norm.pdf(z)
fig, ax = plt.subplots(dpi=100, figsize=[6, 3])
ax.set(title="標準正規分布", xlabel="Z", xlim=(-4, 4))
ax.plot(z, y, color="dimgray")
ax.axhline(y=0, color="dimgray", linewidth=1)
alpha = 0.05 / 2
for a in [alpha, (1 - alpha)]:
x = norm.ppf(a)
ax.axvline(x=x, color="steelblue")
if x < 0:
ax.text(x - 0.1, norm.pdf(x) + 0.01, f"Z = {x:.2f}", color="steelblue", horizontalalignment="right")
ax.fill_between(z, 0, y, where = z <= x, color="steelblue")
else:
ax.text(x + 0.1, norm.pdf(x) + 0.01, f"Z = {x:.2f}", color="steelblue", horizontalalignment="left")
ax.fill_between(z, 0, y, where = z >= x, color="steelblue")
標準正規分布は
一般化してこのような値を
となり、これを
であり、信頼区間は
となる。
pythonでの実装#
母集団が
Show code cell source
np.random.seed(0)
n = 100
x = np.random.uniform(low=0, high=1, size=n)
mu = 0.5 # 母平均
fig, ax = plt.subplots()
ax.hist(x)
ax.set(title="Histogram of Data")
ax.axvline(x.mean(), color="darkorange")
ax.text(x.mean() + 0.02, 1, f"標本平均: {x.mean():.3f}", color="darkorange", size=14)
ax.axvline(mu, color="limegreen")
ax.text(mu + 0.02, 3, f"母平均: {mu:.3f}", color="limegreen", size=14)
fig.show()
式をpythonに落とし込んで計算すると次のようになる
# ※ddof=1: 不偏標準偏差にするためのオプション
std_error = x.std(ddof=1) / np.sqrt(n)
# 信頼区間
alpha = 0.05
z = norm.ppf(1 - alpha / 2)
[x.mean() - z * std_error, x.mean() + z * std_error]
[0.4160030960878336, 0.5295845829372018]
statsmodelsのemplike.DescStatで計算することもできる
# statsmodelsを使う場合
from scipy.stats import norm
import statsmodels.api as sm
el = sm.emplike.DescStat(x)
el.ci_mean()
(0.4166429503391242, 0.5294769896986984)
scipy.stats の norm.interval() で計算することもできる。
from scipy.stats import norm, sem
norm.interval(confidence=0.95, loc=x.mean(), scale=sem(x))
# ※ sem は standard error of mean、つまり x.std(ddof=1) / np.sqrt(n)と等しい
(0.4160030960878336, 0.5295845829372018)
サンプルをとって95%信頼区間を計算する作業を100回繰り返したものが以下の図である。
Show code cell source
np.random.seed(0)
fig, ax = plt.subplots(dpi=90, figsize=[6, 3.5])
ax.set(title="信頼区間", xlabel="平均", ylabel="試行回数")
ax.axvline(mu, color="limegreen", label="母平均μ")
is_first_success = True
is_first_fail = True
n_trial = 100
for i in range(n_trial):
x = np.random.uniform(size=n)
lower, upper = norm.interval(confidence=0.95, loc=x.mean(), scale=sem(x))
args = dict(y=i, xmin=lower, xmax=upper)
if mu < lower or upper < mu:
# 凡例表示のため初回だけlabelをつける
if is_first_fail:
ax.axhline(**args, color="crimson", label="信頼区間（失敗）")
is_first_fail = False
else:
ax.axhline(**args, color="crimson")
else:
if is_first_success:
ax.axhline(**args, color="steelblue", label="信頼区間（成功）")
is_first_success = False
else:
ax.axhline(**args, color="steelblue")
ax.legend()
fig.show()
母分散が未知の場合、
は自由度
となる。