母平均の検定

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

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 でも，その他の言語でも移植は容易であろう。

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}}}$

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

1.7961114275463794

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

df = n - 1

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

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

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

0.08255481102990395

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))

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

検定関数を定義しよう。

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

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 という規定は守られているのだろうか?

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

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

U = var(bread)

18.17726315789474

std(bread)

4.263480169755073

testofmu(100, bar_X, U, 20)

(2.3706049649293197, 19, 0.028491369743615328)

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

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

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

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

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() を使うと

> 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() を使う

>>> 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> 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. オンラインで計算できるサイトもある

やってみた。

4. 結論

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

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

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