母平均の検定 -- p 値を使おう

# 母平均の検定

検定理論の細々した所は省略する。


# 2. 検定手順

例題1

「31 人の身長の平均値が 157.8 cm，不偏分散が 24.6 であった。母平均が 156.2 cm であるかどうかを，有意水準 5% で検定しなさい。」

帰無仮説 $H_0$：「$\mu=\mu_{0}$（母平均は $\mu_{0}$ である）」。

対立仮説 $H_1$：「$\mu \ne \mu_{0}$（母平均は $\mu_{0}$ ではない）」。

有意水準 $\alpha$ で両側検定を行う（片側検定も定義できる）。

母分散が既知か未知かでいずれかの方法をとる。

### 2.1. 母分散が既知の場合

母分散が既知であるというのは，ほとんどありえない。実際に使用することも稀有であろう。

なので，記述は必要ないので省略する。

### 2.2. 母分散が未知の場合

母平均，ケース数，標本平均，標本不偏分散を，それぞれ $\mu_{0}$，$n$，$\bar{X}$，$U$ とする。

例題では，$\mu_{0} = 156.2$，$n = 31$，$\bar{X} = 157.8$，$U = 24.6$ である（母分散 $\sigma^{2}$ は不明である）。

以下の実行例は Julia であるが，R でも Python でも，その他の言語でも移植は容易であろう。

```julia
mu_0 = 156.2
n = 31
bar_X = 157.8
U = 24.6;
```

次式で $t_{o}$ 検定統計量を計算する。

$\hspace{1cm}t_o = \displaystyle \frac{\left | \bar{X}-\mu_0 \right |} {\displaystyle \sqrt{\frac{U}{n}}}$



```julia
t_o = abs(bar_X - mu_0) / sqrt(U / n)
```




    1.7961114275463794



$t_{o}$ は，自由度が $n - 1$ の $t$ 分布に従う（導出過程は不要）。


```julia
df = n - 1
```




    30



例題では，自由度が $30$ の $t$ 分布に従う。

有意確率（$p$ 値）を $P = \Pr\{｜t｜ \geqq t_{o}\}$ とする。


```julia
using Distributions
p_value = 2 * ccdf(TDist(df), t_o)
```




    0.08255481102990395




```julia
using Distributions, Plots
pyplot(grid=false, tick_direction=:out, linecolor=:black, label="")
x = range(-4, 4, length=500)
y = pdf.(TDist(df), x)
plot(x, y, xlabel=raw"$t$ value", ylabel=raw"$f(t, df)$", size=(340, 230))
x2 = range(t_o, 4, length=500)
y2 = pdf.(TDist(df), x2)
plot!(vcat(x2, 4, t_o, t_o...), vcat(y2, 0, 0, y2[1]...),
    seriestype=:shape, fillcolor=:violetred1)
plot!(-vcat(x2, 4, t_o, t_o...), vcat(y2, 0, 0, y2[1]...),
    seriestype=:shape, fillcolor=:violetred1)
hline!([0], linecolor=:black)
plot!([2.3, -2.2], [0.25, 0.02], arrow=1)
plot!([2.3, 2.2], [0.25, 0.02], arrow=1)
annotate!(2.3, 0.25, text(raw" $p$ value", :left, 10))
```




    
![output_8_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/086f20d3-497b-265c-974f-4aeb40536990.png)

    



$p \gt 0.05$ ゆえ，有意水準 5% のもとで，帰無仮説は棄却できない（すなわち，母平均が 156.2 cm でないとはいえない）。

検定関数を定義しよう。

母平均 $\mu_0$，ケース数 $n$，標本平均 $\bar{x}$，標本不偏分散 $U$ を与えれば，検定統計量 $t$，自由度 $df$，$p$ 値を返す。


```julia
function testofmu(mu_0, bar_X, U, n)
    t_o = abs(bar_X - mu_0) / sqrt(U / n)
    df = n - 1
    p_value = 2 * ccdf(TDist(df), t_o)
    return t_o, df, p_value
end

testofmu(mu_0, bar_X, U, n)
```




    (1.7961114275463794, 30, 0.08255481102990395)



# 3. 応用例

例題2

あるパンの重量は 100g と規定されている。20 個のパンの重量を実際に測ったところ，以下のようであった。

 99.8, 95.7,  94.9, 101.9, 105.5, 107.9, 103.4, 110.7,  96.3, 105.1
106.0, 99.2, 102.6,  99.7, 106.4, 102.4, 104.0,  97.9, 105.5, 100.3

本当に 100g という規定は守られているのだろうか?

平均値，不偏分散，標準偏差は以下のようになっていた。


```julia
bread = [ 99.8, 95.7,  94.9, 101.9, 105.5, 107.9, 103.4, 110.7,  96.3, 105.1,
         106.0, 99.2, 102.6,  99.7, 106.4, 102.4, 104.0,  97.9, 105.5, 100.3]
bar_X = mean(bread)
```




    102.26




```julia
U = var(bread)
```




    18.17726315789474




```julia
std(bread)
```




    4.263480169755073




```julia
testofmu(100, bar_X, U, 20)
```




    (2.3706049649293197, 19, 0.028491369743615328)



$p$ 値は 0.02849 なので，5％ の有意水準の両側検定では，パンの重量の母平均は 100g とはいえないという結論になる。

平均で 2g 強しか違わないのにねえ。

ちなみに，上のデータは次のようにして作られたものだった。

母平均=103，不偏分散=25（標準偏差=5）の 20 個の正規乱数（小数点以下 1 桁）である。


```julia
using Random; Random.seed!(123)
bread = rand(Normal(103, 5), 20)
round.(bread, digits=1)|>println
```

    [99.8, 95.7, 94.9, 101.9, 105.5, 107.9, 103.4, 110.7, 96.3, 105.1, 106.0, 99.2, 102.6, 99.7, 106.4, 102.4, 104.0, 97.9, 105.5, 100.3]

### 3.1. R で t.test() を使うと

```R
> bread = c(99.8, 95.7,  94.9, 101.9, 105.5, 107.9, 103.4, 110.7,  96.3, 105.1,
+          106.0, 99.2, 102.6,  99.7, 106.4, 102.4, 104.0,  97.9, 105.5, 100.3)
> t.test(bread,mu=100)

	One Sample t-test

data:  bread
t = 2.3706, df = 19, p-value = 0.02849
alternative hypothesis: true mean is not equal to 100
95 percent confidence interval:
 100.2646 104.2554
sample estimates:
mean of x 
   102.26 
```

`p-value = 0.02849` と同じ結果になった（ホッ）

### 3.2. Python では，scipy.stats.ttest_1samp() を使う

```Python
>>> import scipy.stats as stats
>>> import numpy as np
>>> bread = np.array([99.8, 95.7,  94.9, 101.9, 105.5, 107.9, 103.4, 110.7,  96.3, 105.1,
... +          106.0, 99.2, 102.6,  99.7, 106.4, 102.4, 104.0,  97.9, 105.5, 100.3])
>>> stats.ttest_1samp(bread, 100)
Ttest_1sampResult(statistic=2.3706049649293197, pvalue=0.028491369743615335)
```

`pvalue=0.028491369743615335`　と同じ結果になった（ホッ）

### 3.3.  Julia では HypothesisTests.OneSampleTTest() が既にある

```Julia
julia> using HypothesisTests

julia> bread = [99.8, 95.7,  94.9, 101.9, 105.5, 107.9, 103.4, 110.7,  96.3, 105.1,
                106.0, 99.2, 102.6,  99.7, 106.4, 102.4, 104.0,  97.9, 105.5, 100.3];

julia> OneSampleTTest(bread, 100)
One sample t-test
-----------------
Population details:
    parameter of interest:   Mean
    value under h_0:         100
    point estimate:          102.26
    95% confidence interval: (100.3, 104.3)

Test summary:
    outcome with 95% confidence: reject h_0
    two-sided p-value:           0.0285

Details:
    number of observations:   20
    t-statistic:              2.3706049649293193
    degrees of freedom:       19
    empirical standard error: 0.9533431480294685
```

当然ながら，他の言語のパッケージを使用した場合の結果と同じになる

`two-sided p-value:           0.0285`

### 3.4. オンラインで計算できるサイトもある

https://keisan.casio.jp/exec/user/1490619220

やってみた。

![スクリーンショット 2022-05-07 20.49.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2611739/16b03010-a7bf-2d69-32a6-7df5526bee94.png)

# 4. 結論

どの言語のパッケージでも，臨界値ではなく，$p$ 値を求めていることに留意すべし。

コンピュータを使う我々にとって，「近似」とか「簡単な場合を考える」とか「単純化する」などは無縁のものだ。

- $t$ 分布は $z$ 分布で近似する必要はない。
- Fisher の正確検定は $2\times2$ 分割表に限らない。
- 独立二標本の平均値の差の検定は，等分散性の検定などせず,Welch の方法を直接使えばよい。
- 独立 $k$ 標本の平均値の差の検定も，Welch の方法を使う。
- ノンパラメトリック検定も正確な $p$ 値を求めることができる。
- 母比率の検定は二項検定をすればよい。
- サンプルサイズが大きいからといって，ノンパラメトリック検定を行わないという選択肢はない。
- パラメトリック検定であっても，並べ替え検定で正確な $p$ 値を求めることができる。
- 自前の検定手順（関数化したとしても）をわざわざ使わない（検定の意味を知るために労力を使うべし。数式を導出したり，それを理解したりする必要はさらさらない）。