統計量は、標本抽出によって生成されます。
母集団 U から標本 S を抽出すると、
数学ではこれを、標本から実数への写像といいます。
不偏分散もカイ二乗値も標本からの写像です。
\(U \sim N(\mu,\sigma^2)\)、\(S=\{X_1,\cdots,X_n\}\) を \(U\) からの標本とすると、
これらの統計量は、標本抽出に実数を対応させる写像となっています。
フィッシャーの考えたF統計量は、
「2つの標本抽出に対して分散の推定量の比を対応させる2変量の写像」
であり、2つの標本の分散比較に用いることができます。
\(U \sim N(\mu,\sigma^2)\)、\(S_1=\{X_1,\cdots,X_n\}\)、 \(S_2=\{Y_1,\cdots,Y_m\}\) は \(U\) からの標本とすると、
は分散の推定量であり、
\(F(S_1,S_2) \rightarrow
F=\dfrac{\sigma_1^2}{\sigma_2^2}\) はF統計量です。
フィッシャーは、この統計量によって 2群 \(S_1,S_2\)
が共通の分散を持つか検定しようと考えました。
独立2群の等分散検定(F検定)です。
分散の推定量は1種類ではありません。このことは F統計量に様々な応用の可能性を与えてくれます。
不偏分散も分散の推定量です。
は分散の推定量であり、
\(F(S_1,S_2) \rightarrow
F=\dfrac{s_1^2}{s_2^2}\) もF統計量です。
ただ、平方和の自由度が異なるため、分散の除数が異なっています。
自由度は偏差平方和から分散を求めるための除数と考えて差し支えありません。
不偏分散の自由度は、(標本数−1)です。
そこで、Fに自由度を付記すると、
\(\dfrac{s_1^2}{s_2^2}\sim
F(n-1,m-1)\) \(\dfrac{\sigma_1^2}{\sigma_2^2}\sim
F(n,m)\)
と表すことができます。
Rでは F統計量の確率密度関数は、
df( x, 分子の自由度, 分母の自由度)
という関数で表されます。 この関数を用いれば、
\(\dfrac{s_1^2}{s_2^2}\sim f(x,n-1,m-1)\) \(\dfrac{\sigma_1^2}{\sigma_2^2}\sim F(x,n,m)\)
と表すことができます。
以下、特に断わらない限り、母分散の推定には不偏分散を用います。
Fがどのような確率分布に従うか、シミュレーションしましょう。
簡単のため、2群の標本数が等しい場合を扱います。
これがフィッシャーの大きな発見であり、様々な応用へと広がっていきます。
以下に示すのは、標本分散を用いたF値のシミュレーションです。
平均や分散を変えてシミュレーションすると、上記のことが確認できるでしょう。
S1 = c()
S2 = c()
F12 = c()
n <- 25;m <- 25
for (i in 1:1000) {
x <- rnorm(n,5,2)
y <- rnorm(m,5,2)
S1 <- c(S1,var(x))
S2 <- c(S2,var(y))
F12 <- c(F12,S1/S2)
}
hist(F12,prob=T,xlim=c(0,4))
curve(df(x,n-1,m-1),add=T)
mean(F12)
## [1] 1.097895
Fの期待値は1ですが、少し右に広がった分布になります。
(理論曲線の導出は数学者に任せることにしましょう)
Rでは、分布曲線は df(x,n-1,m-1) という関数で表されます。
ヒストグラムとこの関数の曲線はほぼ一致しているようです。
分散に差があるときのF値の挙動を見てみましょう。
分散の大きい方を分子にとります。
n <- 25;m <- 25
S1 = c()
S2 = c()
G12 = c()
for (i in 1:1000) {
x <- rnorm(n,5,2.2)
y <- rnorm(m,5,2)
S1 <- c(S1,var(x))
S2 <- c(S2,var(y))
G12 <- c(G12,S1/S2)
}
hist(G12,prob=T,xlim=c(0,4))
curve(df(x,n-1,m-1),add=T)
abline(v = qf(0.95,n-1,m-1),col="red")
mean(G12)
## [1] 1.308278
ヒストグラムで、赤い縦線より右の範囲では帰無仮説が棄却されることになり、この範囲の割合を検出力といいます。
分散に差があると理論曲線より大きい値をとるので、分散の差を検出することができます。しかし、この例のように分散の差が小さい場合、検出力は良くないようです。
このように、2群の分散が等しいかどうかは、分散の推定値の比の値を見ることで判断できます。
[課題].
標本数を変えたり、標準偏差の差を変えたりしてシミュレーションし、分散の差が検出可能かどうかを調べましょう。
標本平均の分散は、
\(V(\bar{X})=\dfrac{\sigma^2}{n}\) で与えられます。
標本抽出を無限回繰り返し標本平均の確率分布を調べると、もとの母分散の1/n
になっているというものです。
実際の調査で無限回の抽出を行うことは不可能ですが、乱数でのシミュレーションは可能です。
フィッシャーは、
「標本平均の実験から群間分散の発想を得たのではないか。」
というのが筆者の想像です。
\(U \sim N(\mu,\sigma^2)\)
を母集団とします。
n個の標本抽出を5回行い、得られた標本を \(S_1,S_2,S_3,S_4,S_5\) とするとき、
その標本平均 \(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\)
の分散の推定値は \(\dfrac{\sigma^2}{n}\) で与えられます。
確かめてみましょう。
標本数10の標本を抽出し標本平均を計算する操作を5回行います。
S = c()
for (i in 1:1000) {
s1 <- mean(rnorm(10,5,1))
s2 <- mean(rnorm(10,5,1))
s3 <- mean(rnorm(10,5,1))
s4 <- mean(rnorm(10,5,1))
s5 <- mean(rnorm(10,5,1))
s= var(c(s1,s2,s3,s4,s5))
S<- c(S,s)
}
mean(S)
## [1] 0.09771242
mean(S)*10
## [1] 0.9771242
標本平均の分散は \(\dfrac{\sigma^2}{n}\)
であるので、n倍すると母分散の推定値になります。
上記の例では、
mean(S)*10 が母分散 1 の推定値になっています。
また、5個のデータ \(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\) から求めているので、自由度は4です。
n個の標本抽出を2回行い、得られた標本を \(S_1,S_2\) とします。
標本数が2個でも、標本平均 \(\bar{S}_1,\bar{S}_2\) の分散の推定値は
\(\dfrac{\sigma^2}{n}\) です。
S = c()
for (i in 1:1000) {
s1 <- mean(rnorm(10,5,1))
s2 <- mean(rnorm(10,5,1))
s= var(c(s1,s2))
S<- c(S,s)
}
mean(S)
## [1] 0.08851925
mean(S)*10
## [1] 0.8851925
2つの標本平均から推定した分散は、推定精度は下がりますが、ほぼ母分散に一致しています。
また、2個のデータ \(\bar{S}_1,\bar{S}_2\) から求めているので、自由度は1です。
このように、n群の標本平均から推定した母分散を群間分散と呼びます。
これを、各群の平均の差の検定に利用しようと考えたのがフィッシャーです。
母集団からn群の標本を抽出したとき、n群の標本平均から推定した母分散の推定量が群間分散です。(群間分散の自由度は 群数−1).
これに対し、n群の分散の加重平均を群内分散と呼びます。群内分散は抽出したn群をまとめた分散のことで母分散の推定量です。プールした分散とも呼ばれます。(群内分散の自由度は データ数ー群数).
フィッシャーは、 群間分散/群内分散 もまたF分布に従うことを発見しました。
シミュレーションで確かめてみましょう。
(ここでは、各群の標本数は同じとしておきます)
標本 \(S_1,S_2\) の標本数はともに10 とすると、
分散比は自由度(1,18)のF分布に従うので、確率密度関数は df(x, 1, 18)
になるはずです。
群間分散 sr 、群内分散 se を1000回出力したものを SR,SE、分散比を F12
として期待値を求めます。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(10,5,1);x2 <- rnorm(10,5,1)
sr <- 10*var(c(mean(x1),mean(x2))) #群間分散
SR <- c(SR,sr)
se <- mean(var(x1),var(x2)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 1, 18),add=T)
VR <- c(mean(SR),mean(SE),mean(F12)) #(群間分散、群内分散、分散比)
VR
## [1] 1.079136 0.995622 1.371000
がわかります。
自由度(1,18)のF分布は \(F(1,18)\)
で表されます。
標本 \(S_1,S_2,S_3,S_4,S_5\) の標本数はともに10 とすると、
分散比は自由度(4,45)のF分布に従うので、確率密度関数は df(x, 4, 45)
になるはずです。
群間分散 sr 、群内分散 se を1000回出力したものを SR,SE、分散比を F12
として期待値を求めます。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(10,5,1);x2 <- rnorm(10,5,1);x3 <- rnorm(10,5,1)
x4 <- rnorm(10,5,1);x5 <- rnorm(10,5,1)
sr <- 10*var(c(mean(x1),mean(x2),mean(x3),mean(x4),mean(x5))) #群間分散
SR <- c(SR,sr)
se <- mean(var(x1),var(x2),var(x3),var(x4),var(x5)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 4, 45),add=T)
VR <- c(mean(SR),mean(SE),mean(F12)) #(群間分散、群内分散、F値)
VR
## [1] 1.0162532 0.9936127 1.3242028
ヒストグラムは、df(x, 4, 45) の曲線に一致しています。
自由度(4,45)のF分布は \(F(4,45)\)
で表されます。
これまでのことを整理します。
母集団からいくつかの標本を抽出するとしましょう。
各標本が同じ母集団からのものであれば、群間分散の推定値は母分散になり、群内分散との比はF分布に従います。
ある標本が別の母集団からのものであった場合、見抜くことができるでしょうか。
次の手順で判断します。
等分散検定でF値が大きくなり、等分散という帰無仮説が棄却されれば、別の母集団からの標本が混じっていると判断します。
等分散という帰無仮説が棄却されなければ、各標本の分散は共通と考えます。
各標本の分散は共通として、各標本の平均が異なっていた場合、別の母集団からの標本が混じっていることになるでしょうか。
各標本の平均の差が偶然誤差の範囲であるか、群間分散と群内分散の比を用いて判断します。
各標本の平均の差が偶然誤差の範囲を超えていれば、群間分散の値は跳ね上がり、F値が大きくなることで別の母集団からの標本が混じっていると判断されます。
まず「2群に平均差がある場合」について、シミュレーションで確かめます。
分散は共通、平均差を1とします。
平均差がなければF値は自由度(1,18)のF分布に従い、期待値は1ですが、平均差があればF値は大きくなるはずです。
df(x, 1, 18)
の理論曲線から、ヒストグラムがどの程度離れるかも見ておきましょう。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(10,6,1);x2 <- rnorm(10,5,1)
sr <- 10*var(c(mean(x1),mean(x2))) #群間変動
SR <- c(SR,sr)
se <- mean(var(x1),var(x2)) #群内変動
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 1, 18),add=T)
abline(v = qf(0.95,1,18),col="red")
VR <- c(mean(SR),mean(SE),mean(F12)) #(群間分散、群内分散、F値)
VR
## [1] 5.9166200 0.9789699 7.6588214
平均が異なると群間変動が跳ね上がりF値が大きくなります。
ヒストグラムで、赤の縦線より右の値の場合に有意差が検出されますが、赤の縦線の左の部分も大きいので、検出力が十分とはいえません。
ある程度の標本数がないと、平均差が検出できない可能性があります。
標本数や平均の差を変えてシミュレーションし、有効な調査のために必要なサンプル数を確かめておくことが重要です。
標本 \(S_1,S_2\) の標本数はともに 25 とすると、
平均差がなければF値は自由度(1,48)のF分布に従います。
分散は共通、平均差を1とします。
標本数が増えると、F値が大きくなり検出力が高くなることを確かめましょう。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(25,6,1);x2 <- rnorm(25,5,1)
sr <- 25*var(c(mean(x1),mean(x2))) #群間変動
SR <- c(SR,sr)
se <- mean(var(x1),var(x2)) #群内変動
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 1, 48),add=T)
abline(v = qf(0.95,1,48),col="red")
VR <- c(mean(SR),mean(SE),mean(F12)) #(群間分散、群内分散、F値)
VR
## [1] 13.47747 1.00800 14.49195
赤の縦線の右にあるものがほとんどになり、
検出力が高くなっていることがわかります。
標本数が大きければ、群の平均差を検出できる可能性が大きくなります。
独立2群の差の検定においてフィッシャーの画期的なところは、
ということです。
また、2群の差の検定から多群の差の検定への拡張が容易なこともフィッシャーの方法の優れているところです。
ただし、等分散の検定を多群に拡張するには、バートレット検定やハートレイの検定などが必要になります。ここでは立ち入りません。
[参考] 池田郁男先生 統計検定を理解せずに使っている人のためにIII.
簡単のため、3群の標本数は共通で 25
としてシミュレーションします。
分散は共通、平均差を1とします。
群内分散 \(\dfrac{24var(S_1)+24var(S_2)+24var(S_3)}{24\times
3}\) なので自由度は 72.
群間分散 \(\bar{S}_1,\bar{S}_2,\bar{S}_3\) の分散なので、自由度は 2
平均差がなければF値は自由度(2,72)のF分布に従い、期待値は1ですが、平均差があればF値は大きくなるはずです。 df(x, 2, 72) の理論曲線から、ヒストグラムがどの程度離れるかも見ておきましょう。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(25,6,1);x2 <- rnorm(25,5,1);x3 <- rnorm(25,5,1)
sr <- 25*var(c(mean(x1),mean(x2),mean(x3))) #群間分散
SR <- c(SR,sr)
se <- mean(var(x1),var(x2),var(x3)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 2, 72),add=T)
abline(v = qf(0.95,2,72),col="red")
mean(F12)
## [1] 10.448
3つの群に平均の異なる群があると、F値が大きくなります。
上記の例では、標本 x1 のみ平均が異なりましたが、3群とも平均を変えてみます。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(25,6,1);x2 <- rnorm(25,5,1);x3 <- rnorm(25,4,1)
sr <- 25*var(c(mean(x1),mean(x2),mean(x3))) #群間分散
SR <- c(SR,sr)
se <- mean(var(x1),var(x2),var(x3)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 2, 72),add=T)
abline(v = qf(0.95,2,72),col="red")
mean(F12)
## [1] 28.12351
検出力が大きくなります。
次に、平均差を小さくしてみます。
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(25,5.5,1);x2 <- rnorm(25,5,1);x3 <- rnorm(25,5,1)
sr <- 25*var(c(mean(x1),mean(x2),mean(x3))) #群間分散
SR <- c(SR,sr)
se <- mean(var(x1),var(x2),var(x3)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
hist(F12,prob=T,breaks = "Scott",main="",xlab="")
curve(df(x, 2, 72),add=T)
abline(v = qf(0.95,2,72),col="red")
mean(F12)
## [1] 3.474884
検出力は小さくなります。 平均差を検出するには、ある程度平均差があることが必要なようです。
この方法は分散分析と呼ばれ、多群の比較分析に活用されています。
分散分析は、各群の差が偶然変動の範囲か偶然を超えた差かを判断する重要な方法です。
フィッシャーは1910年代に、t分布に関するゴセットの論文に興味を持ち、その数理を一般化したと言われています。
ゴセットの t値の二乗は、フィッシャーから見れば分散比でした。
実際、 \(E\left((\bar{X}-\mu)^2\right)=\dfrac{\sigma^2}{n}\) \(E\left(\dfrac{s^2}{n}\right)=\dfrac{\sigma^2}{n}\) であるので、
とすれば、\(t^2\) は母分散の推定値の比になります。
定義を最後に持ってきた理由は、定義が一般化されフィッシャーの分散比とのつながりが見えにくいものになっているからです。
フィッシャーの考えたF統計量は、
2つの標本 \(S_1=\{X_1,\cdots,X_n\}\)、
\(S_2=\{Y_1,\cdots,Y_m\}\)
から求めた母分散の推定量の比でした。
分散の推定量は偏差平方和を自由度で割った、
F統計量は、\(F=\dfrac{\sigma_1^2}{\sigma_2^2}=\dfrac{\left. \sigma_1^2 \right/ \sigma^2}{\left. \sigma_2^2 \right/ \sigma^2}\) と考えれば、
\(\left. \sigma_1^2 \right/ \sigma^2 =\left. \left(\dfrac{X_1-\mu}{\sigma}\right)^2+\left(\dfrac{X_2-\mu}{\sigma}\right)^2+\cdots +\left(\dfrac{X_n-\mu}{\sigma}\right)^2\right/ n\)
\(\left. \sigma_2^2 \right/ \sigma^2 =\left. \left(\dfrac{Y_1-\mu}{\sigma}\right)^2+\left(\dfrac{Y_2-\mu}{\sigma}\right)^2+\cdots +\left(\dfrac{Y_m-\mu}{\sigma}\right)^2\right/ m\).
これらの分子はカイ二乗統計量であり、自由度 n のカイ二乗統計量を \(\chi^2(n)\) と書くことにすれば、
F統計量は、\(F=\dfrac{\sigma_1^2}{\sigma_2^2}=\dfrac{\left. \sigma_1^2 \right/ \sigma^2}{\left. \sigma_2^2 \right/ \sigma^2}=\dfrac{\left. \chi^2(n) \right/ n}{\left. \chi^2(m) \right/ m}\)
この最後の式が数理統計学でのF統計量の定義です。
抽象化によって分かりにくくなりますが、抽象化によって応用範囲が広がることが数学の不思議です。
群間分散と群内分散をカイ二乗統計量で表す変形には立ち入りません。
このように数理的導出が難しい事柄は、Rによるシミュレーションで確かめることが有効です。
群間分散と群内分散の比が、
\(F=\dfrac{\left. \chi^2(n) \right/
n}{\left. \chi^2(m) \right/
m}\) の形で表されることを確かめましょう。
[カイ二乗分布の定義].
2章「3. 母分散の推定」の課題で示したように、
\(U=\dfrac{(X_1-\bar{X})^2+(X_2-\bar{X})^2+\cdots +(X_n-\bar{X})^2}{\sigma^2}\)
は自由度 (n-1) のカイ二乗分布に従います。
つまり、偏差平方和を母分散で割ったものは、カイ二乗統計量です。
「3.群内分散との比」の抽出5回のコードを利用します。
[群間分散].
標本 \(S_1,S_2,S_3,S_4,S_5\)
の標本数はともに10 とすると、
自由度は4なので、偏差平方和は \(4
var\left(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\right)\).
これを標本平均の分散 \(\frac{\sigma^2}{10}\) で割った値、
\(\left. 4
var\left(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\right)
\right/
\frac{\sigma^2}{10}=40var\left(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\right)\).
は自由度4のカイ二乗統計量で、
\(\chi^2(4)=40var\left(\bar{S}_1,\bar{S}_2,\bar{S}_3,\bar{S}_4,\bar{S}_5\right)\).
[群内分散].
群内分散は加重平均 \(\frac{9var(S_1)+9var(S_2)+\cdots
+9var(S_5)}{9\times 5}\) なので自由度は 45.
分子は偏差平方和なので、これを 分散 \(\sigma^2=1\) で割った、
\(\frac{9var(S_1)+9var(S_2)+\cdots
+9var(S_5)}{1}\).
は自由度 45 のカイ二乗統計量で、
\(\chi^2(45)=9var(S_1)+9var(S_2)+\cdots
+9var(S_5)\).
SE = c()
SR = c()
F12= c()
for (i in 1:1000){
x1 <- rnorm(10,5,1);x2 <- rnorm(10,5,1);x3 <- rnorm(10,5,1)
x4 <- rnorm(10,5,1);x5 <- rnorm(10,5,1)
sr <- 40*var(c(mean(x1),mean(x2),mean(x3),mean(x4),mean(x5))) #群間分散
SR <- c(SR,sr)
se <- 9*(var(x1)+var(x2)+var(x3)+var(x4)+var(x5)) #群内分散
SE<- c(SE,se)
F12 <- SR/SE
}
par(mfrow=c(1,2))
hist(SR,prob=T,breaks = "Scott",main="dchisq(x,4)",xlab="")
curve(dchisq(x,4),add=T)
hist(SE,prob=T,breaks = "Scott",main="dchisq(x,45)",xlab="")
curve(dchisq(x,45),add=T)
VR <- c(mean(SR),mean(SE)) #(群間偏差平方和、群内偏差平方和)
VR
## [1] 4.087195 44.702659
分子と分母の分散が異なっていても、各々の分散で割ることでカイ二乗値になるのは驚くべきことで、数学における一般化の凄さを感じます。
自由度は理解が難しいですが、このように偏差平方和の期待値を求めれば、それが自由度になっています。
最後に、2つのカイ二乗統計量を用いて商を作ると、F統計量になります。
hist((SR/4)/(SE/45),prob=T,breaks = "Scott",main="df(x, 4, 45)",xlab="")
curve(df(x, 4, 45),add=T)
F値 <- mean((SR/4)/(SE/45))
F値
## [1] 1.076942
母集団の分散によらず、F統計量の期待値は1になります。
[課題].
抽出2回の場合について、シミュレーションしましょう。
統計検定を理解せずに使っている人のために II によると.
「データ数が少ない場合は等分散性の検定はあてにできない」とされています。
環境と品質のためのデータサイエンス様によると.
「少ないサンプル数では、極端に違う場合しか判定できない」とされています。