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

```py3
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本で同じ結果の場合の区間はどうなるか。

```py3
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]

```py3
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

```py3
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]

```py3
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)$