正規分布の分散の検定、カイ二乗検定とは：Juliaによる求め方

どうも、木村（@kimu3_slime）です。

今回は、正規分布の分散の検定、カイ二乗検定について、Juliaによる求め方を紹介します。

分散の検定、カイ二乗検定の原理
Juliaによる求め方
こちらもおすすめ

分散の検定、カイ二乗検定の原理

10人のテストの点数を測ったところ、分散が約83でした。点数が正規分布に従っていると仮定するとき、分散は$\sigma ^2 = 80$という仮説は正しいでしょうか？

このようなパラメータの値に関する仮説は、検定によって調べることができます。平均の検定ではt分布を用いましたが、分散の検定ではカイ二乗分布が使えます。

$X_1,\dots, X_n$を一般の正規分布に従う独立な確率変数族（サンプル）とします。$\sigma ^2$を母集団の分散、$S^2$をサンプルの不偏分散とするとき、統計量

\[Y= (n-1)\frac{S^2}{\sigma ^2}\]

は自由度$n-1$のカイ二乗分布に従うことが知られています。これはカイ二乗検定統計量（chi-squared test statistic）とも呼ばれるものです。

つまり、$F$を自由度$n-1$のカイ二乗分布の累積分布関数とするとき、$c_1,c_2$を実数として

\[P(c_1 \leq Y \leq c_2 ) =F(c_2)-F(c_1)\]

が成り立ちます。

これを変形して、サンプルの不偏分散$S ^2$の区間の形にしましょう。

左辺の中身を変形すると、

\[ c_1 \leq (n-1)\frac{S^2}{\sigma^2}\leq c_2\]

\[ c_1 \sigma ^2 \frac{1}{n-1} \leq S^2 \leq c_2 \sigma ^2 \frac{1}{n-1} \]

となります。

一方、右辺について考えましょう。与えられた有意水準$\alpha$に対し、$F(c_1) = \frac{\alpha}{2}$、$F(c_2)=1- \frac{\alpha}{2}$を満たす$c_1,c_2$を求めましょう。つまり、$c_1 = F^{-1}(\frac{1-\gamma}{2})$、$c_2 = F^{-1}(\frac{1+\gamma}{2})$とします。すると、$F(c_2)-F(c_1)=1-\alpha$ですね。

よって、サンプル分散が棄却されない区間

\[P(c_1 \sigma ^2 \frac{1}{n-1} \leq S^2 \leq c_2 \sigma ^2 \frac{1}{n-1})=1- \alpha\]

を導くことができました。

カイ二乗分布を使った分散の検定の流れは、次の通りです。

データが正規分布に従っていると仮定する。
帰無仮説を$H_0: \sigma ^2 =80$とする。対立仮説を$H_1 : \sigma^2 \neq 80$とする（両側検定）。
有意水準を決める。今回は$\alpha = 0.05$とする。
$c_1,c_2$を求め、$\sigma^2 =80$のとき（帰無仮説が正しいとき）、サンプル分散が棄却されない区間$c_1 \sigma ^2 \frac{1}{n-1} \leq S^2 \leq c_2 \sigma ^2 \frac{1}{n-1}$を求める。
データを使って、検定の考え方にもとづいた判断をする。
1. サンプル分散が棄却されない区間に含まれるなら、帰無仮説は棄却されない。
2. 棄却されない区間に含まれない（棄却域に含まれる）なら、棄却域は棄却される。

この検定は、正規分布の分散のカイ二乗検定（chi-squared test）と呼ばれます。

ただし、カイ二乗分布を使った検定には適合度検定などさまざまなものがあり、文脈によって指しているものが違うことがあるので注意しましょう。

Juliaによる求め方

より具体的に、コンピュータ、Juliaを使って分散の検定を行ってみましょう。

正規分布に従う乱数により、10個のサンプルを作りました。ここでは母集団分布の分散を与えていますが、それを未知として進めましょう。

using Distributions, Random

d1 = Normal(50,10)
Random.seed!(2022)
x = rand(d1,10)

using Distributions, Random

d1 = Normal(50,10)

Random.seed!(2022)

x = rand(d1,10)

10-element Vector{Float64}:
 50.6912301164913
 38.84399961848597
 42.84480604223728
 60.9769650429804
 54.978651946367876
 33.951339022096505
 40.96806064167559
 32.42960541746047
 49.11909337153176
 48.368593011534614

10-element Vector{Float64}:

50.6912301164913

38.84399961848597

42.84480604223728

60.9769650429804

54.978651946367876

33.951339022096505

40.96806064167559

32.42960541746047

49.11909337153176

48.368593011534614

サンプルの分散は

var(x)

var(x)

83.71338280293834

1	83.71338280293834

です。このとき、帰無仮説を$H_0 : \sigma ^2 =80$、対立仮説を$H_1 : \sigma ^2 \neq 80$として、検定を行いましょう。

自由度$n-1$のカイ二乗分布は「Chisq(n-1)」、累積分布関数の逆関数は「quantile」で求められます。

サンプルと帰無仮説を与えて、帰無仮説が棄却されないか棄却されるかを判別する関数を作りました。

function OneSampleChisqTest(x, σ_0sq)
    n = length(x)
    α = 0.05
    c_1 = quantile(Chisq(n-1), α/2)
    c_2 = quantile(Chisq(n-1), 1-α/2 )
    LL = c_1 * σ_0sq / (n-1)
    UL = c_2 * σ_0sq / (n-1)
    
    println("H_0: σ^2 = $σ_0sq")
    println("サンプルサイズ ",n)
    println("サンプル分散 ",var(x))
    println("棄却されない区間 ","($LL, $UL)")
    if LL <= var(x) <= UL
        println("fail to reject H_0")
    else
        println("reject H_0")
    end
end

function OneSampleChisqTest(x, σ_0sq)

n = length(x)

α = 0.05

c_1 = quantile(Chisq(n-1), α/2)

c_2 = quantile(Chisq(n-1), 1-α/2 )

LL = c_1 * σ_0sq / (n-1)

UL = c_2 * σ_0sq / (n-1)

println("H_0: σ^2 = $σ_0sq")

println("サンプルサイズ ",n)

println("サンプル分散 ",var(x))

println("棄却されない区間 ","($LL, $UL)")

if LL <= var(x) <= UL

println("fail to reject H_0")

else

println("reject H_0")

end

実行してみると、

H_0: σ^2 = 80
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (24.003462222047624, 169.09126932125895)
fail to reject H_0

H_0: σ^2 = 80

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (24.003462222047624, 169.09126932125895)

fail to reject H_0

$\sigma ^2 =80$という仮説は棄却されないことがわかりました。

仮説をいくつか変えて、結果を見てみましょう。

for θ in 20:100:320
    OneSampleChisqTest(x, θ)
    println("\n")
end

for θ in 20:100:320

OneSampleChisqTest(x, θ)

println("\n")

end

H_0: σ^2 = 20
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (6.000865555511906, 42.27281733031474)
reject H_0


H_0: σ^2 = 120
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (36.00519333307144, 253.63690398188842)
fail to reject H_0


H_0: σ^2 = 220
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (66.00952111063097, 465.0009906334621)
fail to reject H_0


H_0: σ^2 = 320
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (96.0138488881905, 676.3650772850358)
reject H_0

H_0: σ^2 = 20

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (6.000865555511906, 42.27281733031474)

reject H_0

H_0: σ^2 = 120

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (36.00519333307144, 253.63690398188842)

fail to reject H_0

H_0: σ^2 = 220

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (66.00952111063097, 465.0009906334621)

fail to reject H_0

H_0: σ^2 = 320

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (96.0138488881905, 676.3650772850358)

reject H_0

$\sigma ^2 =20,320$のとき、帰無仮説は棄却されます。一方、$\sigma ^2 =120,220$のときは棄却されません。サンプルサイズが小さいので、棄却されない仮説の範囲は広いですね。

検定において、棄却されない仮説は、「唯一正しい仮説」や「最も適切な仮説」を意味しないことに注意しましょう。

仮説を$\sigma ^2 =80$としたまま、同じサンプルでサンプル数を増やすとどうなるか試してみましょう。

for n in 10:100:210
    Random.seed!(2022)
    x = rand(d1,n)
    OneSampleChisqTest(x, 80)
    println("\n")
end

for n in 10:100:210

Random.seed!(2022)

x = rand(d1,n)

OneSampleChisqTest(x, 80)

println("\n")

end

H_0: σ^2 = 80
サンプルサイズ 10
サンプル分散 83.71338280293834
棄却されない区間 (24.003462222047624, 169.09126932125895)
fail to reject H_0


H_0: σ^2 = 80
サンプルサイズ 110
サンプル分散 93.75717377033588
棄却されない区間 (60.183253758042504, 102.5936862022903)
fail to reject H_0


H_0: σ^2 = 80
サンプルサイズ 210
サンプル分散 97.63422672233817
棄却されない区間 (65.39905022111957, 96.05013490145926)
reject H_0

H_0: σ^2 = 80

サンプルサイズ 10

サンプル分散 83.71338280293834

棄却されない区間 (24.003462222047624, 169.09126932125895)

fail to reject H_0

H_0: σ^2 = 80

サンプルサイズ 110

サンプル分散 93.75717377033588

棄却されない区間 (60.183253758042504, 102.5936862022903)

fail to reject H_0

H_0: σ^2 = 80

サンプルサイズ 210

サンプル分散 97.63422672233817

棄却されない区間 (65.39905022111957, 96.05013490145926)

reject H_0

サンプルサイズが110では棄却されませんが、210となると棄却されました。もしサンプルサイズが大きく取れるならば、仮説の精度を上げることができますね。

同じサンプルサイズと同じ仮説で、有意水準を変えるとどうなるか試しましょう。

function OneSampleChisqTest(x, σ_0sq, α)
    n = length(x)
    c_1 = quantile(Chisq(n-1), α/2)
    c_2 = quantile(Chisq(n-1), 1-α/2 )
    LL = c_1 * σ_0sq / (n-1)
    UL = c_2 * σ_0sq / (n-1)
    
    println("H_0: σ^2 = $σ_0sq")
    println("サンプルサイズ ",n)
    println("サンプル分散 ",var(x))
    println("有意水準α ", α)
    println("棄却されない区間 ","($LL, $UL)")
    if LL <= var(x) <= UL
        println("fail to reject H_0")
    else
        println("reject H_0")
    end
end

function OneSampleChisqTest(x, σ_0sq, α)

n = length(x)

c_1 = quantile(Chisq(n-1), α/2)

c_2 = quantile(Chisq(n-1), 1-α/2 )

LL = c_1 * σ_0sq / (n-1)

UL = c_2 * σ_0sq / (n-1)

println("H_0: σ^2 = $σ_0sq")

println("サンプルサイズ ",n)

println("サンプル分散 ",var(x))

println("有意水準α ", α)

println("棄却されない区間 ","($LL, $UL)")

if LL <= var(x) <= UL

println("fail to reject H_0")

else

println("reject H_0")

end

A = [0.05, 0.01]
for α in A
    Random.seed!(2022)
    x = rand(d1,210)
    OneSampleChisqTest(x, 80, α)
    println("\n")
end

A = [0.05, 0.01]

for α in A

Random.seed!(2022)

x = rand(d1,210)

OneSampleChisqTest(x, 80, α)

println("\n")

end

H_0: σ^2 = 80
サンプルサイズ 210
サンプル分散 97.63422672233817
有意水準α 0.05
棄却されない区間 (65.39905022111957, 96.05013490145926)
reject H_0


H_0: σ^2 = 80
サンプルサイズ 210
サンプル分散 97.63422672233817
有意水準α 0.01
棄却されない区間 (61.28026693556816, 101.59264968428663)
fail to reject H_0

H_0: σ^2 = 80

サンプルサイズ 210

サンプル分散 97.63422672233817

有意水準α 0.05

棄却されない区間 (65.39905022111957, 96.05013490145926)

reject H_0

H_0: σ^2 = 80

サンプルサイズ 210

サンプル分散 97.63422672233817

有意水準α 0.01

棄却されない区間 (61.28026693556816, 101.59264968428663)

fail to reject H_0

$\alpha =0.05$のときは$\sigma ^2 =60$という仮説が棄却されますが、$\alpha =0.01$のときは棄却されませんでした。信頼水準を厳しく（小さく）するほど、棄却されない区間は広くなっていきます。

同じサイズのサンプルを何度も検定して、どのくらいの割合で正しい仮説$H_0 : \sigma ^2 =100$が棄却されるか調べてみましょう。

棄却されないときは0、棄却されるときを1として、回数をカウントしていきます。

function OneSampleChisqTest(x)
    n = length(x)
    α = 0.05
    c_1 = quantile(Chisq(n-1), α/2)
    c_2 = quantile(Chisq(n-1), 1-α/2 )
    LL = c_1 * var(d1) / (n-1)
    UL = c_2 * var(d1) / (n-1)
    
    if LL <= var(x) <= UL
        return 0
    else
        return 1
    end
end

function OneSampleChisqTest(x)

n = length(x)

α = 0.05

c_1 = quantile(Chisq(n-1), α/2)

c_2 = quantile(Chisq(n-1), 1-α/2 )

LL = c_1 * var(d1) / (n-1)

UL = c_2 * var(d1) / (n-1)

if LL <= var(x) <= UL

return 0

else

return 1

end

k = 10^4
sum = 0
for i in 1:k
    x = rand(d1,10)
    sum += OneSampleChisqTest(x)
end
sum/k

k = 10^4

sum = 0

for i in 1:k

x = rand(d1,10)

sum += OneSampleChisqTest(x)

end

sum/k

0.0542

0.0542

1万回サンプルを取って検定を繰り返すと、そのうち約5%で帰無仮説が棄却されました。これが有意水準$\alpha = 0.05$（正しい仮説が棄却される：第一種の過誤の確率）の意味ですね。

以上、正規分布の分散の検定、カイ二乗検定について、Juliaによる求め方を紹介してきました。

検定の考え方は誤解しやすいので、理論と合わせてコンピュータで実験しながらその意味を理解していくと良いでしょう。

木村すらいむ（@kimu3_slime）でした。ではでは。

Probability and Statistics: Pearson New International Edition

posted with AmaQuick at 2022.03.24

DeGroot, Morris(著), Schervish, Mark(著)
Pearson Education Limited (2013-07-30T00:00:01Z)

￥10,792 (中古品)

Amazon.co.jpで詳細を見る

入門数理統計学

posted with AmaQuick at 2022.03.27

P.G.ホーエル(著), 浅井晃(翻訳), 村上正康(翻訳)
培風館 (1978-01-01T00:00:01Z)

￥5,280

Amazon.co.jpで詳細を見る

Advanced Engineering Mathematics

posted with AmaQuick at 2022.03.14

Kreyszig, Erwin(著)
John Wiley & Sons Inc (2011-05-03T00:00:01Z)

￥5,862 (中古品)

Amazon.co.jpで詳細を見る

こちらもおすすめ

カイ二乗分布とその性質、Juliaによる計算方法

統計的仮説検定とは：平均の検定、t検定を例に、Juliaを使って

正規分布の分散の区間推定：カイ二乗分布を使った方法