Pythonで数学の勉強：統計学(仮説検定の基礎と区間推定)

統計学基礎 (仮説検定)

帰無仮説と対立仮説

帰無仮説 ${\rm H}_0$: 棄却されることを前提とした仮説を立てる。
対立仮説 ${\rm H}_1$: 採択されることを前提とした仮説を立てる。

(例題) サイコロを10回振ったら1が5回出た。このサイコロは1が出やすいのだろうか。

${\rm H}_0$: このサイコロの1が出る確率は1/6である。
${\rm H}_1$: このサイコロの1が出る確率は1/6より大きい。

有意水準と信頼度

${\rm H}_0$が棄却される領域を棄却域、採択される領域を採択域という。
${\rm H}_0$が正しいときに棄却域に入る確率を有意水準という。
有意水準1%、5%がよく使われる。
有意水準の逆を信頼度99%、95%ともいう。

片側検定と両側検定

棄却域を分布の片側に取るものを片側検定という。
棄却域を分布の両側に取るものを両側検定という。
例題の場合はサイコロの1が出る確率が1/6より小さいことは問題にしないので片側検定を使います。

仮説の真偽と採択

$1.$ 統計量が棄却域に落ちた場合
${\rm H}_0$が棄却され、${\rm H}_1$が採択される。
「このサイコロの1が出る確率は1/6より大きいといえる。」
$2.$ 統計量が採択域にある場合
${\rm H}_1$が棄却され、${\rm H}_0$が採択される。
「このサイコロの1が出る確率は1/6より大きいとはいえない。」

第1種の誤り、第2種の誤り

第1種の誤りとは${\rm H}_0$が真にもかかわらず棄却されること。
高々有意水準でしか起こらない。
第2種の誤りとは${\rm H}_0$が偽にもかかわらず採択されること。
第2種の誤り確率を$\beta$、$1-\beta$を検出力という。

パラメトリック検定とノンパラメトリック検定

母集団の分布に正規分布などを仮定して検定を行う手法をパラメトリック検定という。
そのような仮定を置かずに検定を行う手法をノンパラメトリック検定という。

区間推定

import

import numpy as np
import scipy

正規分布に従う母集団の区間推定

母平均： $\mu$
母分散： $\sigma^2$
標本分散： $\bar{X}$
標本分散： $s^2$
標本数： $n$

$\bar{X}$は正規分布$N(\mu, \sigma^2/n)$に従う。
$\sqrt{n}(\bar{X}-\mu)/s$は自由度$(n-1)$の$t$分布に従う。
$\mu$の信頼係数$1-\alpha$の信頼区間は

$$\bar{X} \pm t_{\alpha/2} (n-1) \frac{s}{\sqrt{n}} $$

$(n-1)s^2/\sigma^2$は自由度$(n-1)$の$\chi^2$分布に従う。
$\sigma^2$の信頼係数$1-\alpha$の信頼区間は

$$[\frac{(n-1)s^2}{\chi^2_{\alpha/2}(n-1)}, \frac{(n-1)s^2}{\chi^2_{1-\alpha/2}(n-1)}]$$

(問) ねじ100本の長さを測定したところ平均45.625mm, 標準偏差0.188mmであった。
信頼係数95%で平均値と分散の区間を求めよ。
測定本数が10000本で同じ結果の場合の区間はどうなるか。

alpha = 0.95
n1 = 100
n2 = 10000
mean = 45.625
std = 0.188

def confidence_interval_normal(n, mean, std, alpha=0.95):
    sem = std / np.sqrt(n) # 標準偏差を自由度の平方根で割った値をscaleに渡す
    range_mean = scipy.stats.t.interval(alpha, n-1, loc=mean, scale=sem)
    (chi2_lower, chi2_upper) = scipy.stats.chi2.interval(alpha, n-1) # カイ二乗の下値、上値を求める
    temp = (n-1) * std * std
    range_var = (temp / chi2_upper, temp / chi2_lower)
    return [range_mean, range_var]

print(confidence_interval_normal(n1, mean, std, alpha))
print(confidence_interval_normal(n2, mean, std, alpha))

$[(45.587696721311637, 45.662303278688363), (0.0272465489512402, 0.047696353309908554)]$
$[(45.621314821624395, 45.628685178375605), (0.034384386015036454, 0.036344548305711967)]$

2つの正規母集団の母平均の差、母分散の比の区間推定

母平均： $\mu_1, \mu_2$
母分散： $\sigma_1^2, \sigma_2^2$
標本平均： $\bar{X}, \bar{Y}$
標本分散： $s_1^2, s_2^2$
標本数： $n, m$

ウェルチの近似式

$$t=\frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{\sqrt{\frac{s_1^2}{m}+\frac{s_2^2}{n}}}$$
は自由度が
$$\nu = \frac{(\frac{s_1^2}{m}+\frac{s_2^2}{n})^2}{\frac{s_1^4}{m^2(m-1)}+\frac{s_2^4}{n^2(n-1)}}$$
に最も近い整数$\nu*$の$t$分布に近似的に従う。

$\mu_1-\mu_2$の信頼係数$1-\alpha$の信頼区間は
$$\bar{X}-\bar{Y}\pm t_{\alpha/2}(\nu*) \sqrt{\frac{s_1^2}{m}+\frac{s_2^2}{n}}$$

母分散の比

$(\sigma_2^2s_1^2)/(\sigma_1^2s_2^2)$は自由度$(m-1, n-1)$の$F$分布に従う。
$\sigma_2^2/\sigma_1^2$の信頼係数$1-\alpha$の信頼区間は
$$[F_{1-\alpha/2}(m-1, n-1)s_2^2/s_1^2, F_{\alpha/2}(m-1, n-1)s_2^2/s_1^2]$$

(問) 以下の標本から2つの母集団の母平均の差と母分散の比を信頼係数95%で区間推定せよ。
x=[100,98,92,90,87,85,81,75,65,64,63,60,58,29,20,15,13,10,4,0]
y=[100,95,90,70,70,69,60,50,45,40,38,33,32,30,25,25,20,19,2,0]

alpha = 0.95
x = np.array([100, 98, 92, 90, 87, 85, 81, 75, 65, 64, 63, 60, 58, 29, 20, 15, 13, 10, 4, 0])
y = np.array([100, 95, 90, 70, 70, 69, 60, 50, 45, 40, 38, 33, 32, 30, 25, 25, 20, 19, 2, 0])
n = len(x)
m = len(y)
X_ = np.mean(x)
Y_ = np.mean(y)
s1 = np.var(x, ddof=1)
s2 = np.var(y, ddof=1)
# ここから母平均の差を推定する
temp = s1 / m + s2 / n
nu = temp * temp / (s1 * s1 / (m * m * (m - 1)) + s2 * s2 / (n * n * (n - 1)))
t = scipy.stats.t.ppf(alpha + (1 - alpha) / 2, round(nu))
mean_sub = (X_ - Y_ - t * np.sqrt(temp), X_ - Y_ + t * np.sqrt(temp))
# ここから母分散の比を推定する
(f_lower, f_upper) = scipy.stats.f.interval(alpha, m-1, n-1)
var_p = (s2 / (s1 * f_upper), s2 / (s1 * f_lower))

print(X_, Y_)
print(s1, s2)
print(mean_sub)
print(var_p)

$55.45 45.65$
$1192.57631579 852.871052632$
$(-10.69084424953008, 30.290844249530089)$
$(0.28306509921801942, 1.8067915978336075)$

二項分布に従う分布のpの推定

正規分布での近似

試行回数が十分多いかつpが0や1に近くないとき
二項分布$Bi(1, p), 母平均:p, 母分散:p(1-p)$
$\bar{X}: 標本平均$

$p$の信頼係数$1-\alpha$の信頼区間は近似的に
$$\bar{X}\pm Z_{\alpha/2}\sqrt{\bar{X}(1-\bar{X})/n}$$

Clopper-Pearson法

$$[B_{(1-\alpha/2)}(n-k+1, k), B_{(\alpha/2)}(n-k, k+1)]$$

(問) コインをn回投げてx回表が出たとする。以下の場合について表が出る確率pを信頼係数95%で区間推定せよ。
1. n=10, x=7
2. n=10, x=10
3. n=100, x=70
4. n=100, x=100

alpha = 0.95
x1 = 7
n1 = 10
x2 = 70
n2 = 100

def binom_norm(k, n, alpha):
    p = k / n
    alpha2 = alpha + (1 - alpha) / 2
    z = scipy.stats.norm.ppf(alpha2)
    std = np.sqrt(p * (1 - p) / n)
    return (p - z * std, p + z * std)

def clopper_pearson(k, n, alpha):
    alpha2 = (1 - alpha) / 2
    lower = scipy.stats.beta.ppf(alpha2, k, n - k + 1)
    upper = scipy.stats.beta.ppf(1 - alpha2, k + 1, n - k)
    return (lower, upper)

print(f"正規:{x1}/{n1}{binom_norm(x1, n1, alpha)}")
print(f"　CP:{x1}/{n1}{clopper_pearson(x1, n1, alpha)}")
print(f"正規:{n1}/{n1}{binom_norm(n1-1, n1, alpha)}")
print(f"　CP:{n1}/{n1}{clopper_pearson(n1, n1, alpha)}")

print(f"正規:{x2}/{n2}{binom_norm(x2, n2, alpha)}")
print(f"　CP:{x2}/{n2}{clopper_pearson(x2, n2, alpha)}")
print(f"正規:{n2}/{n2}{binom_norm(n2-1, n2, alpha)}")
print(f"　CP:{n2}/{n2}{clopper_pearson(n2, n2, alpha)}")

正規:7/10(0.41597423491067459, 0.98402576508932538)
　CP:7/10(0.3475471499400028, 0.93326048882226553)
正規:10/10(0.71406149030863153, 1.0859385096913685)
　CP:10/10(0.69150289218123928, nan)
正規:70/100(0.6101831668145794, 0.78981683318542051)
　CP:70/100(0.60018532382019585, 0.78759357951046338)
正規:100/100(0.97049860458201209, 1.0095013954179879)
　CP:100/100(0.96378330735482365, nan)

ポアソン分布に従う集団の区間推定

$P(\lambda)$は平均、分散とも$\lambda$である。
nが大きい場合中心極限定理より正規分布で近似でき、
$\lambda$の信頼係数$1-\alpha$の信頼区間は$\bar{\lambda}=\bar{X}$として
$$\bar{\lambda}\pm Z_{(\alpha/2)}\sqrt{\frac{\bar{\lambda}}{n}}$$

(問) 以下の標本の母集団がポアソン分布に従うとして$\lambda$を区間推定せよ。
[14, 22, 16, 18, 17, 15, 16, 20, 12, 10, 13, 23]

alpha = 0.95
x = np.array([14, 22, 16, 18, 17, 15, 16, 20, 12, 10, 13, 23])
n = len(x)
X_ = np.mean(x)
std = np.sqrt(X_ / n)
alpha2 = alpha + (1 - alpha) / 2
z = scipy.stats.norm.ppf(alpha2)
(X_ - z * std, X_ + z * std)

$(14.046708684703269, 18.619957981963395)$