10章

1 前書き.

正規分布しないときに行われる検定方法が「ノンパラメトリック検定」です。
ここでは、データが正規分布に従わないとき、通常の検定にどのような不具合が生じるかをシミュレーションによって確かめ、代替方法としてノンパラメトリック検定を紹介します。

２標本t検定 は、独立２群の平均差を調べるパラメトリックな方法で、「データの正規性」「２群の等分散性」を前提とします。
説明の流れとしては、

• 前提条件を満たす２標本t検定を実施し、
  
• 等分散性を満たさない場合、どんな問題が生ずるか.
  
• 正規性を満たさない場合、どんな問題が生ずるか.

これらを確かめ、パラメトリック検定の限界を知ったうえでノンパラメトリック検定へ進みます。

Wilcoxonの符合順位和検定 は、独立２群の差を順位の和から判断するノンパラメトリックな方法です。
順位は離散変量であり数列を学んだ高校生であれば、期待値や分散などの統計量を数理的に導出することができます。
高校でデータサイエンスが学ばれるようになれば、高校生が統計的検定を学ぶ格好の題材となるでしょう。

2 ２標本 t 検定.

２標本t検定は、独立２群から得られた２標本の平均差から、元の２群の平均に差があるかどうかを調べる方法です。

＜例＞.
２つの学校A,Bの学力に差があるかを調べたいとします。
各校から大きさ40の標本を抽出し比較します。
シミュレーションでは、同じ平均と標準偏差の正規乱数を用いてデータを出力します。

x1 <- rnorm(40,50,10);x2 <- rnorm(40,50,10)
list=list(A=x1, B=x2)
boxplot(list,ylim=c(20,80),horizontal=T)

母平均と母分散が同じであっても、誤差変動によりboxplot（箱ヒゲ図）にはズレが生じています。
しかし、平均点に差があったとしても、それは誤差変動のためかもしれません。

今、平均差は.

mean(x1)-mean(x2)

## [1] -1.203986

同じ母集団からとった２群ですが、平均に差があります。
これが誤差変動の範囲であるか、本当に差があるのかを判断する尺度は平均差の標準偏差です。
差を測る尺度としての標準偏差について整理しておきます。

2.1 マハラノビスの距離.

http://hyoccori-toukei.la.coocan.jp/id-10/id-2/id10.html

マハラノビスの距離は、
　　集団においてデータが集団の中心から確率的にどのくらい離れているか.
を示す値です。（中心は集団の平均）

通常の距離は２点間の距離を表すのに対して、マハラノビスの距離は集団の中心からの確率的距離であることが特徴です。
中心に近ければ起こりやすく、中心から離れると起こりにくいと考えます。

2.1.1 ＜標準正規分布＞.

平均０, 分散１の確率変数の理想的な分布です。
標準正規分布に従うデータを1000個出力しヒストグラムを描きます。

x <- rnorm(1000,0,1)
hist(x,prob=T,breaks = "Scott")
curve(dnorm(x,0,1),add=T)

中心（平均）に近いところでは起こりやすく、中心から遠いところでは起こりにくいことが分かります。

曲線はガウスの発見した理論曲線で、次の式で表されます。

　　　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\sim N(0,1)\).
　　　中心 0 から距離 \(|z|\) にあるデータが出力される確率は \(\phi(z)\) .

中心極限定理では、一般の確率変数でもサンプル数が大きくなれば、標本平均を標準化した確率変数は標準正規分布に漸近的に従うことが知られています。

確率は中心からの距離によって決まります。（中心＝平均のこと）

• 中心に近いところでは確率が高く.
• 中心から遠いところでは確率は低くなります.
  この中心からの距離 \(|z|\) をマハラノビスの距離といいます。

注意.
　z そのものは、マハラノビスの座標と考えればよいでしょう。

2.1.2 ＜正規分布＞.

一般の正規分布では、尺度を標準偏差にすることで統一的な扱いが可能になります。
つまり中心から標準偏差いくつぶん離れているかを確率的距離と考えるわけです。
　　　\(z=\dfrac{x-\mu}{\sigma}\)

　　　標準化によって得られるzの絶対値が xのマハラノビス距離に相当します。

\(x\sim N(\mu,~\sigma^2)\) のとき、

　　　\(z=\dfrac{x-\mu}{\sigma} \sim N(0,1)\)　　であり、

• |z| が小さければ、xは平均から近く高確率、
• |z| が大きければ、xは平均から遠く低確率、

|z| はマハラノビスの距離であり、x が平均から標準偏差いくつぶん離れているかを表します。

• マハラノビスの確率距離が小さければ確率的に起こりやすく、誤差変動の範囲内.
  
• マハラノビスの確率距離が大きければ確率的に起こりにくく、誤差変動を超える範囲.

と考えます。

2.1.3 ＜検定との関係＞.

統計的検定では、帰無仮説のもとでの確率分布を考えます。

有意水準 5 %のZ検定では、得られた統計量と分布の中心とのマハラノビス距離が、

• 1.96以内ならば、統計量は確率95%の範囲内にあり、帰無仮説は棄却されず、
  
• 1.96を超えていれば、統計量は誤差を超える範囲と判断し、帰無仮説を棄却します。

検定では偏差を標準偏差で割った z 値の絶対値を確率的距離として、差が有意かどうかを判断しています。

標準偏差が未知の場合は標準偏差の推定値を使用します。この場合、
偏差を標準偏差の推定値で割った t 値の絶対値を確率的距離とし、t分布表を用いて差が有意かどうかを判断します。

注意.
　z 値はマハラノビスの座標と考えましたが、t 値は確率的座標と考えると良いでしょう。.

2.2 平均差の確率分布.

本節のテーマは「２群の平均差の検定」ですが、差を測るには２群の平均差の標準偏差が必要です。
平均差の期待値・分散について確認しておきましょう。
まず、基本的な性質から確認します。

2.2.1 ＜期待値・分散の性質＞.

確率変数 X,Y が独立であるとき、

• \(E(X-Y)=E(X)-E(Y)\)　　　　\(E(kX)=kE(X)\).

• \(V(X-Y)=V(X)+V(Y)\)　　　　\(V(kX)=k^2V(X)\).

が成立します。（期待値E( )の性質は、独立でなくとも成立します）
期待値の性質は線形性と呼ばれる自然な性質です。
分散の性質は意外に思われるかもしれません。

分散の性質についてシミュレーションで確かめてみましょう。
　　　x1：平均 3 、標準偏差 1 の正規乱数からの1000個のデータ.
　　　x2：平均 2 、標準偏差 1 の正規乱数からの1000個のデータ.
　　　x3=x1-x2.

x1 = rnorm(1000,mean=3,sd=1)
x2 = rnorm(1000,mean=2,sd=1)
x3 = x1-x2
hist(x1,prob=T,col="pink",xlim=c(-2,12))
hist(x2,prob=T,col="lightblue",add=T)
hist(x3,prob=T,col="lightgreen",add=T)

\(x3=x1-x2\) のヒストグラムの幅は広がり、分散は大きくなっているようです。 分散の性質を確認します。

v <- c(var(x3),var(x3/2))
v

## [1] 2.0375635 0.5093909

• \(V(X1-X2)=V(X1)+V(X2)=2\)　　　 \(V\left(\dfrac{1}{2}X3\right)=\dfrac{1}{4}V(X3)=0.5\)

となっていることが分かります。

2.2.2 ＜標本平均の分布＞.

データが \(X\sim N(\mu,\sigma^2)\) に従い、
\(n\)個のデータ\(x_1,~x_2,~\cdots,x_{n}\) が独立に得られているとします。

\(x_1,~x_2,~\cdots,x_{n}\) は独立な \(n\) 個の確率変数であるので、

　　　\(E(x_1+x_2+\cdots +x_{n})=E(x_1)+E(x_2)+\cdots +E(x_{n})\)
　　　\(V(x_1+x_2+\cdots +x_{n})=V(x_1)+V(x_2)+\cdots +V(x_{n})\) 　　が成立します。

平均 \(\mu\),　分散 \(\sigma^2\) は共通なので、
　　　\(E(x_1+x_2+\cdots +x_{n})=n\mu\)
　　　\(V(x_1+x_2+\cdots +x_{n})=n\sigma^2\)

\(E(kX)=kE(X),~V(kX)=k^2V(X)\) であるので、

　　　\(E\left(\dfrac{x_1+x_2+\cdots +x_{n}}{n}\right)=\mu\)
　　　\(V\left(\dfrac{x_1+x_2+\cdots +x_{n}}{n}\right)= \dfrac{\sigma^2}{n}\)

\(n\)個データの標本平均　\(\overline{X}=\dfrac{x_1+x_2+\cdots +x_{n}}{n}\)　について、

　　　\(\overline{X}\sim N\left(\mu,\dfrac{\sigma^2}{n}\right)\)　　が成立します。

注意.
一般の分布に対しても標本数が大きければ、中心極限定理により標本平均の分布は正規近似できることが知られています。
標本数は30以上あれば正規近似は妥当とされていますが、シミュレーションによって確かめておく必要があるでしょう。

2.2.3 ＜標本平均の差の分布＞.

\(X_1\sim N(\mu_1,\sigma_1^2)\:,\:X_2\sim N(\mu_2,\sigma_2^2)\) 　は独立な２群を表す確率変数とします。

各群から標本を抽出し、その標本数を \(n_1,n_2\) 、不偏分散を \(s_1^2,s_2^2\) とします。

各抽出群の標本平均を \(\overline{X_1},~~\overline{X_2}\) とすると、

　　\(\overline{X_1}\sim N\left(\mu_1,\dfrac{\sigma_1^2}{n_1}\right)\:,\:\overline{X_2}\sim N\left(\mu_2,\dfrac{\sigma_2^2}{n_2}\right)\)　より

　　　　　\(\overline{X_1}-\overline{X_2}\sim N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma_1^2}{n_1}+\dfrac{\sigma_2^2}{n_2}\right)\)

母分散が共通：\(\sigma^2=\sigma_1^2=\sigma_2^2\) という前提のもとでは、

　　　　　\(\overline{X_1}-\overline{X_2}\sim N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)

母分散が未知の場合、次のように母分散の推定値 \(s^2\) を構成します。

• ２群をプールした偏差平方和は 　　\(s_1^2(n_1-1)+s_2^2(n_2-1)\)

• 自由度は 　　\((n_1-1)+(n_2-1)=n_1+n_2-2\)

• ２群をプールした不偏分散は 　\(s^2=\dfrac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}\)

これが母分散の推定値です。

注意.
一般の分布に対しても標本数が大きければ、中心極限定理により標本平均の分布は正規近似できるので、標本平均の差の分布は正規分布で近似できます。

2.2.4 ＜２群の平均差の検定＞.

２群が等分散であることを前提として、

帰無仮説：「２群の平均に差はない」　を検定します。

　　　\(\overline{X_1}-\overline{X_2}\sim N\left(0\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)　であるので、

＜母分散 \(\sigma^2\) が既知の場合＞、

　　統計量　\(z=\dfrac{\overline{X_1}-\overline{X_2}}{\sigma\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim N(0\:,\:1)\) をマハラノビスの座標と考え、

　　得られた統計量 \(\overline{X_1}-\overline{X_2}\) が有意な差かどうかを判断します。

＜母分散 \(\sigma^2\) が未知の場合＞.
　母分散の推定量 \(s^2\) を用い、

　　統計量　\(t=\dfrac{\overline{X_1}-\overline{X_2}}{s\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim t(n_1+n_2-2)\) 　を確率的座標と考え、

　　t分布表で、統計量 \(\overline{X_1}-\overline{X_2}\) が有意な差かどうかを判断します。

注意.
マハラノビスの距離も確率的距離も本来は正の数ですが、上の説明では、z値とt値は標準偏差を尺度とした座標と考えています。つまり、中心から左へ行けば、マイナスの距離と考えています。

2.3 シミュレーション①.

目的は、２群の平均差を適切に検定することです。
適切でない検定であっても、検定結果は出力されてしまいます。結果が信頼に値するためには、事前に十分なシミュレーションを行い、検定の可能性と限界について調べておく必要があります。 以下のシミュレーションにおいても、データ数や分布などの前提条件を変え、納得いくまで確かめてください。

2.3.1 母分散が既知.

母分散を用いた尺度を用いると、Z検定になります。
簡単のため、２群の標本数は同じとします。
また、検定は右側検定のみを実施します。

2.3.1.1 標本数40.

＜帰無仮説＞.

• 群G1（平均5, 標準偏差1）からの大きさ40の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ40の標本 \(x_2\)

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0,\dfrac{1^2}{40}+\dfrac{1^2}{40}\right)\)
\(\overline{X_1}-\overline{X_2}\) の母標準偏差を尺度としたマハラノビスの距離は、
標準正規分布 \(N(0,1)\) に従うはずです。
標準正規分布の理論曲線を赤で重ね書きして確かめましょう。

X=c()
XZ=c()
for (i in 1:1000){
  x1 <- rnorm(40,5,1);x2 <- rnorm(40,5,1)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  xz <- x/sqrt(1^2/40+1^2/40)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dnorm(x, 0, 1),add=T,col="red")

＜対立仮説＞.

• 群G1（平均5.5, 標準偏差1）からの大きさ40の標本x1.
  
• 群G2（平均5, 標準偏差1）からの大きさ40の標本x2.

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0.5,\dfrac{1^2}{40}+\dfrac{1^2}{40}\right)\)

このように、２群に平均差があるとしたら、\(\overline{X_1}-\overline{X_2}~~\) の 0 からのマハラノビスの距離はどうなるでしょうか。

帰無仮説と対立仮説の中心間の距離は 0.5 ですが、
マハラノビス距離は \(\dfrac{0.5}{\sqrt{\frac{1}{40}+\frac{1}{40}}}\fallingdotseq 2.24~\)

Z検定では、有意差が検知される可能性が高いようです。
ヒストグラムで視覚化してみましょう。

• 帰無仮説のもとでの理論曲線を青で示します。
  
• 赤の縦線は帰無仮説のもとでの棄却境界であり、
  ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。 　　

X=c()
Y=c()
XZ=c()
YZ=c()
for (i in 1:1000){
  y1 <- rnorm(40,5.5,1);y2 <- rnorm(40,5,1)
  y <- mean(y1)-mean(y2)
  yz <- y/sqrt(1^2/40+1^2/40)
  YZ <- c(YZ,yz)
}

hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dnorm(x, 0, 1),col="blue",add=T)
abline(v = qnorm(0.95,0,1),col="red")

ヒストグラムの赤い縦線より右の範囲では帰無仮説が棄却されることになります。
\(\overline{X_1}-\overline{X_2}\) の０からのマハラノビス距離の期待値は 2.24 ですが、誤差変動のため棄却域に入らないこともあるようです。

2.3.1.2 標本数10.

• 標本数が小さくても、母集団が正規分布に従えば、標本平均も正規分布に従います.
  
• 母分散を用いた標準化は、標準正規分布に従います。

帰無仮説のもとでのマハラノビス座標の分布は、標準正規分布に従うはずです。確かめておきましょう。

X=c()
XZ=c()
for (i in 1:1000){
  x1 <- rnorm(10,5,1);x2 <- rnorm(10,5,1)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  xz <- x/sqrt(1^2/10+1^2/10)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dnorm(x, 0, 1),add=T,col="red")

マハラノビスの座標が標準正規分布に従うことは、問題ないようです。

＜対立仮説＞.

• 群G1（平均5.5, 標準偏差1）からの大きさ10の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ10の標本 \(x_2\)

のもとで、
　　標本平均の差 　 \(\overline{X_1}-\overline{X_2}~~\sim N\left(0.5,\dfrac{1^2}{10}+\dfrac{1^2}{10}\right)\)

このように、２群に平均差があるとしたら、\(\overline{X_1}-\overline{X_2}~~\) の 0 からのマハラノビスの座標はどうなるでしょうか。

帰無仮説と対立仮説の中心間の距離は 0.5 ですが、
マハラノビス距離は \(\dfrac{0.5}{\sqrt{\frac{1}{10}+\frac{1}{10}}}\fallingdotseq 1.12~\)

標本数が小さいと、尺度である標準偏差が大きくなり、Z検定では棄却されない可能性が高いようです。
ヒストグラムで視覚化してみましょう。

• 帰無仮説のもとでの理論曲線を青で示します。
  
• 赤の縦線は帰無仮説のもとでの棄却境界であり、
  ヒストグラムで赤の縦線の右側の値が得られたとき帰無仮説は棄却されます。 　　

X=c()
Y=c()
XZ=c()
YZ=c()
for (i in 1:1000){
  y1 <- rnorm(10,5.5,1);y2 <- rnorm(10,5,1)
  y <- mean(y1)-mean(y2)
  yz <- y/sqrt(1^2/10+1^2/10)
  YZ <- c(YZ,yz)
}

hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dnorm(x, 0, 1),col="blue",add=T)
abline(v = qnorm(0.95,0,1),col="red")

ヒストグラムの赤い縦線より右の範囲は少なく、有意差が検出される可能性は低いと言えます。 平均の有意差を検出するには、標本数が大きい方が有利です。

2.3.2 母分散が未知.

乱数モデルによるシミュレーションでは母分散が分かっていますが、実際の検定では母分散が未知として進めなければなりません。
以下では、母分散の推定値を用いて確率距離を計算します。
この場合、確率的座標はt分布に従います。

2.3.2.1 標本数20.

標本数が大きければ標本平均は正規分布で近似され分散の推定値も母分散に近似されるので、t分布と標準正規分布との差はあまりないはずです。シミュレーションで確かめましょう。

＜帰無仮説＞.

• 群G1（平均5, 標準偏差1）からの大きさ20の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ20の標本 \(x_2\)

母標準偏差は未知であるとし、母標準偏差の推定値 s を用いて標準化します。
t分布の確率密度関数を青で、標準正規分布を赤で描き比較します。

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,1)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

t検定の棄却境界を青の縦線で、z検定の棄却境界を赤の縦線で表示していますが、ほとんど差がないようです。（標本数40にすると完全に重なって見えます）.

＜対立仮説＞.

• 群G1（平均5.5, 標準偏差1）からの大きさ20の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ20の標本 \(x_2\)

20個の標本で、平均差 0.5 を検知できるでしょうか。
対立仮説のもとでの出力をヒストグラムで、帰無仮説のもとでの理論曲線を青で表示します。

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 20
n2 <- 20
for (i in 1:1000){
  y1 <- rnorm(n1,5.5,1);y2 <- rnorm(n2,5,1)
  y <- mean(y1)-mean(y2)
  s <- (var(y1)*(n1-1)+var(y2)*(n2-1))/(n1+n1-2)
  yz <- y/sqrt(s/n1+s/n2)
  YZ <- c(YZ,yz)
}

hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="red")

検出力（ヒストグラムで赤縦線の右の割合）は 50% ほどで低いですが、これは標本数が20 で小さいため有意差が検出しにくいからです。
標本数を40とすれば、Z検定とほぼ同じ検出力になります。

つまり、母分散を推定値で代用し t検定を使用しても、検定の有効性に変わりはないと考えられます。
しかし、標本数が少ない場合は検出力が低下するので、注意深く観察する必要があります。

2.3.2.2 標本数5.

同じ t検定ですが、標本数を 5 に減らしてみましょう。

＜帰無仮説＞.

• 群G1（平均5, 標準偏差1）からの大きさ5の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ5の標本 \(x_2\)

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 5
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,1)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

標本数が小さいと、t分布の裾は標準正規分布よりも広がり、棄却境界もZ検定から離れてきます。
この状態でZ検定を行えば検定は甘くなります。推定や検定を精密にするには、t検定がより適切と考えられます。

＜対立仮説＞.

• 群G1（平均5.5, 標準偏差1）からの大きさ5の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）からの大きさ5の標本 \(x_2\)

小標本で、平均差 0.5 が検出できるでしょうか？

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 10
n2 <- 10
for (i in 1:1000){
  y1 <- rnorm(n1,5.5,1);y2 <- rnorm(n2,5,1)
  y <- mean(y1)-mean(y2)
  s <- (var(y1)*(n1-1)+var(y2)*(n2-1))/(n1+n2-2)
  yz <- y/sqrt(s/n1+s/n2)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dt(x,(n1+n2-2)),col="blue",add=T)
abline(v = qt(0.95,n1+n2-2),col="red")

t分布への当てはまりは良いですが、小標本のため有意差の検出が難しくなっています。
小標本で有意差を検出するには、母平均の差がかなり大きい必要があります。

3 ヴェルチの検定.

3.1 ヴェルチの標準化.

　標本平均から母平均を推定する際、高校では母標準偏差が既知であることが前提となっています。
高校生の頃、「母標準偏差は既知」という前提に違和感を覚えました。分散は母平均が分からなければ求まらないからです。 母分散が未知の場合の母平均推定については、母分散の推定値を用いたゴセットの標準化がその方法を与えてくれます。
（拙著「高校生の統計学」では、不偏分散を用いた標本平均の標準化をゴセットの標準化と呼んでいます）.

　次に、独立２群の平均差を検定する際には、２群の等分散性が必要になります。 しかし、母集団から抽出した標本の分散変動は大きく、抽出標本による分散の推定値が異なっているからといって、２群の等分散性が否定されるわけではありません。
例えば、「母平均5、母分散100の母集団から大きさ10の標本を抽出し不偏分散を算出する」という操作を1000回繰り返し、不偏分散の変動を見ると、

XV=c()
n1 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,10)
  XV <- c(XV,var(x1))
}
mv <- c(mean(XV),sd(XV))
mv

## [1] 98.80381 45.37832

不偏分散の平均は100に近いですが、変動幅である標準偏差は 50 に近く、かなり大きな変動となっています。そこで、
　　　「等分散という前提を設けない検定を考えてはどうか」　というのがヴェルチのアイディアです。

方法としてはゴセットの標準化の類比を考えます。

＜標本平均から母平均を推定する際＞.

• 母分散が既知のときは.
  \(\dfrac{\bar{X}-\mu}{\left. \sigma \right/ \sqrt{n}}\) が標準正規分布に従うことから区間推定が可能でした。

• 母分散が未知のとき.
  \(\dfrac{\bar{X}-\mu}{\left. s \right/ \sqrt{n}}\) が t分布に従うことから区間推定を行います。
  （本書ではゴセットの標準化と呼んでいます）.
  母分散を推定値で代用し t分布につなげています。

t分布へのつなげ方を、次の場合にも適用します。

＜独立２群の平均差を検定する際＞.

\(X_1\sim N(\mu_1,\sigma_1^2)\:,\:X_2\sim N(\mu_2,\sigma_2^2)\) 　は独立な２群を表す確率変数とします。

各群から標本を抽出し、その標本数を \(n_1,n_2\) 、不偏分散を \(s_1^2,s_2^2\) とします。

各群抽出サンプルの標本平均を \(\overline{X_1},~~\overline{X_2}\) とすると、

　　　\(\overline{X_1}~\sim N\left(\mu_1,\dfrac{\sigma_1^2}{n_1}\right)\:,\:\overline{X_2}~\sim N\left(\mu_2,\dfrac{\sigma_2^2}{n_2}\right)\)　より

　　　　　\(\overline{X_1}-\overline{X_2}~\sim N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma_1^2}{n_1}+\dfrac{\sigma_2^2}{n_2}\right)\)

• 母分散 \(\sigma_1^2,\sigma_2^2\) が既知のときは.

　　　\(z=\dfrac{\overline{X_1}-\overline{X_2}}{\sqrt{\dfrac{\sigma_1^2}{n_1}+\dfrac{\sigma_2^2}{n_2}}}\sim N(0\:,\:1)\)　により２群の平均差をZ検定します。

• 母分散が未知のときは、母分散を推定値 \(s_1^2,s_2^2\) で代用し.
  \(t=\dfrac{\overline{X_1}-\overline{X_2}}{\sqrt{\dfrac{s_1^2}{n_1}+\dfrac{s_2^2}{n_2}}}\sim t(f)\) 　により２群の平均差をt検定します。
  これを本書ではヴェルチの標準化と呼ぶことにします。
  （自由度 \(f\fallingdotseq n_1+n_2-2\)と考えられますが、以下の補正が必要になります）
  　補正自由度：\(f=\dfrac{(s_1^2/n_1 + s_2^2/n_2)^2 }{ \dfrac{(s_1^2 / n_1)^2}{ (n_1 - 1) }+ \dfrac{(s_2^2 / n_2)^2}{(n_2 - 1) }}\).
  補正自由度は２群の標本数と標本の不偏分散に影響されます。標本数が同数の場合と大きく異なる場合では、自由度も大きく変わることが予想されます。シミュレーションでは、標本数の状況を変えて調べる必要があるでしょう。
  自由度が変動する確率変数となることで、シミュレーションが困難になりますが、自由度の平均を用いて t分布への近似の状態を把握します。

数理的導出については、
@tabintone(komatsu) 様「Welchのt検定ってどこから来たの？」に詳しい説明があります。
かなり難しいです。
ここではシミュレーションを用いて妥当性を確かめることにしましょう。

参考サイト
https://ikuro-kotaro.sakura.ne.jp/koramu/toukei2.htm

3.1.1 不等分散・母分散既知

母分散が既知であれば、２群の分散が大きく異なっていても大丈夫でしょうか？
標準化したz値が標準正規分布に適合しているかに注目します。

＜\(n_1=5,\sigma_1=1、 n_2=15,\sigma_2=10\)＞.

XW=c()
n1 <- 5;s1 <- 1
n2 <- 15;s2 <- 10
for (i in 1:1000){
  x1 <- rnorm(n1,5,s1);x2 <- rnorm(n2,5,s2)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(s1^2/n1+s2^2/n2)
  XW <- c(XW,xw)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

母分散が既知であれば、

• 分散が大きく異なっていても、
  
• 標本数が大きく異なっていても、

標本差の標準化は標準正規分布に従うと考えて良さそうです。

3.1.2 不等分散・母分散未知

大標本であれば、母分散の推定値は母分散にほぼ等しく、標本平均は漸近的に正規分布に従うので、ヴェルチの標準化は標準正規分布に漸近的に従うと考えられます。検討が必要になるのは、

• 小標本.
• 不等分散.
• 異なる標本数.

のときです。

母分散が未知の場合、母分散の推定には次の２つの方法が考えられます。

＜２群の標本変動をプールした分散で代用＞

次のように母分散の推定値 \(s^2\) を構成します。

• ２群をプールした偏差平方和は 　　\(s_1^2(n_1-1)+s_2^2(n_2-1)\)
  
• 自由度は 　　\((n_1-1)+(n_2-1)=n_1+n_2-2\)

　２群をプールした不偏分散　\(s^2=\dfrac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}\) 　を母分散の推定量として標準化します。
　この標準化は２群の母分散が等しいとき 自由度 \((n_1+n_2-2)\)の t分布に従いますが、不等分散の場合はその保証はありません。

＜２群の不偏分散で代用＞

　　　\(s^2=\dfrac{s_1^2}{n_1}+\dfrac{s_2^2}{n_2}\)　を母分散の推定量とするヴェルチの標準化です。
この場合も t分布に従う保証はありません。ヴェルチは自由度を
　　　\(f=\dfrac{(s_1^2/n_1 + s_2^2/n_2)^2 }{ \dfrac{(s_1^2 / n_1)^2}{ (n_1 - 1) }+ \dfrac{(s_2^2 / n_2)^2}{(n_2 - 1) }}\)　と補正することを主張しています。
　　　 　

この２つの方法をシミュレーションで比較検討し、t分布からどの程度乖離するか確かめます。.

3.1.3 ヴェルチの自由度.

ヴェルチの自由度が確率変数であることにも注意する必要があります。
つまり近似されるt分布が常に変化するので、分布の当てはまりの判断ができなくなります。
ヴェルチの自由度がどのような値を取るか調べてみましょう。

XW=c()
FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(FF,prob=T)
abline(v = n1+n2-2,col="red")
abline(v = mean(FF),col="red")

XW=c()
FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,3)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(FF,prob=T)
abline(v = n1+n2-2,col="red")
abline(v = mean(FF),col="red")

補正自由度はプール分散の自由度より小さく、２群の分散の差が大きくなければ、プール分散の自由度を補正して小さくしたものと考えることができそうです。シミュレーションでは、変動する自由度の平均を自由度としたt分布を考えます。

3.2 シミュレーション②.

3.2.1 不等分散かつ等平均.

• プール分散による標準化をxp.
• ヴェルチの標準化をxw.

として1000回の繰り返しデータをヒストグラムで並置し比較します。

3.2.1.1 \(n_1=4,\sigma_1=1、 n_2=7,\sigma_2=5\).

XW=c()
XP=c()
FF=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,mean(FF)),add=T,col="green")
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,mean(FF)),add=T,col="green")

#補正自由度、プール自由度
free <- c(mean(FF),n1+n2-2)
free

## [1] 6.945462 9.000000

プール分散を推定値として用いた平均差の標準化は、理論分布から乖離しています。
プール分散は、２群の分散の加重平均です。
分散の大きい群の標本数が大きければ、プール分散は過大になるので、標準化は小さな値になりがちです。 分散の小さい群の標本数が大きければ、プール分散は過小になるので、標準化は大きな値になると予想されます。

（標本数が等しい場合は、乖離度は小さい：平均分散になるため？）

3.2.1.2 \(n_1=10,\sigma_1=1、 n_2=4,\sigma_2=5\).

XW=c()
XP=c()
FF=c()
n1 <- 10
n2 <- 4
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-8,8))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,mean(FF)),add=T,col="green")
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-8,8))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,mean(FF)),add=T,col="green")

#補正自由度、プール自由度
free <- c(mean(FF),n1+n2-2)
free

## [1]  3.258067 12.000000

標本数が変わると、補正自由度・プール自由度の差も大きくなっています。
注目したいのは、ヴェルチの標準化は標本数の影響が小さいことです。
母分散が既知でない場合にヴェルチの方法が有効とされるのも納得がいきます。

課題.
２群の標本数・分散を変えてシミュレーションしましょう。
標本数が同じ場合は、プール分散による標準化もヴェルチの標準化もあまり差がないと予想されます。

3.3 検定の有効性.

２群の分散が異なる場合、 ヴェルチの標準化が有効であることが分かりましたが、標準正規分布や通常のt分布とどのくらい乖離しているでしょうか。

ヴェルチの標準化のヒストグラムと理論曲線を比べてみます。
赤は標準正規分布、青はプール分散の自由度によるt分布、緑はヴェルチの自由度によるt分布です。

＜帰無仮説＞.

• 群G1（平均5, 標準偏差1）からの大きさ15の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差3）からの大きさ5の標本 \(x_2\)

XW=c()
FF=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,5,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

free <- c(mean(FF),n1+n2-2)
free

## [1]  4.621646 18.000000

ヒストグラムへの当てはまりは、ヴェルチの自由度による曲線のほうが良く、棄却境界もかなり異なるので、ヴェルチの検定の有効性が高いと考えられます。

通常は、

• ２群の等分散検定を行い.
• 棄却されなければ、プール分散を用いたt検定.
• 棄却されれば、ヴェルチの検定.

とされていました。しかし、検定を繰り返すことで有意差が出やすくなり、第１種の過誤が増加する恐れがあります。シミュレーションの結果を見れば、直ちにヴェルチの検定を実施するほうが有効性が高いようです。

＜対立仮説＞.

• 群G1（平均8, 標準偏差1）からの大きさ15の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差3）からの大きさ5の標本 \(x_2\)

XW=c()
FF=c()
n1 <- 15
n2 <- 5
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
free <- c(mean(FF),n1+n2-2)
free

## [1]  4.567134 18.000000

hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

検出力は高いとは言えないようです。 　

分散が大きい群の標本を大きくした対立仮説でシミュレーションします。

＜対立仮説＞.

• 群G1（平均8, 標準偏差1）からの大きさ5の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差3）からの大きさ15の標本 \(x_2\)

XW=c()
FF=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
#自由度
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
free <- c(mean(FF),n1+n2-2)
free

## [1] 16.38312 18.00000

hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-3,3))
curve(dt(x,mean(FF)),add=T,col="green")
abline(v = qt(0.95,f),col="green")

分散が大きい群の標本数が大きいほうが、検出力は高くなるようです。

参考.
層化抽出法では、分散の大きい群からの抽出数を大きくすると、集団の状況が標本に反映されやすいということが知られており、配分の比率をネイマン配分といいます。

3.4 P値の分布の比較

ヴェルチの検定の有効性を確かめるため、以下のシミュレーションで、対立仮説のもとでのp値を算出し、プール分散によるt検定のp値の分布と比較します。
p値の分布を調べる理由は、ヴェルチのt検定では補正自由度が抽出ごとに変化するため、 どの自由度のt分布に従うのか固定できないためです。

ここでのシミュレーションは、分散が大きい群の標本数を大きくして行います。

＜\(n_1=5,\sigma_1=1、 n_2=15,\sigma_2=3\)＞.

PP=c()
PW=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,8,1);x2 <- rnorm(n2,5,3)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean

## [1] 0.09548768 0.01204315

PP=c()
PW=c()
n1 <- 5
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,4)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean

## [1] 0.2552570 0.1130431

＜\(n_1=15,\sigma_1=1、 n_2=5,\sigma_2=2\)＞.

PP=c()
PW=c()
n1 <- 10
n2 <- 15
for (i in 1:1000){
  x1 <- rnorm(n1,7,1);x2 <- rnorm(n2,5,2)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  pw <- 1- pt(xw,f)
  PW <- c(PW,pw)
}
par(mfrow=c(1,2))
par(family = "HiraKakuProN-W3") 
hist(PP,prob=T,breaks = "Scott",main="等分散仮定でのP値",xlab="",ylim=c(0,12))
hist(PW,prob=T,breaks = "Scott",main="ヴェルチt検定のP値",xlab="",ylim=c(0,12))

mean <- c(mean(PP),mean(PW))
mean

## [1] 0.02528680 0.01240895

ここでは、分散の大きい群の標本が大きい場合に、ヴェルチのt検定の有効性を説明しました。
標本数が同じ場合や分散の大小が逆の場合はどうなるでしょうか。
研究課題として、シミュレーションで確かめてみてください。

4 正規性の問題.

z検定やt検定などのパラメトリック検定では、データの正規性が前提となっています。

大標本の場合、中心極限定理により母集団の分布によらず標本平均は正規分布に近似されるので、検定の有効性に問題は生じません。（頑強性と言われます）.

しかし、実際の調査研究では、小標本が用いられることがしばしばです。
ゴセットは発酵タンクのビール酵母の数を推定する過程でt分布を発見したと言われています。どのような推定方法であったかは想像するしかありませんが、標本調査による商品ロスを少なくするためには、小標本で推定することが重要だったと考えられます。当時は推定に正規近似を用いていましたが、ゴセットの出発点は「小標本に正規近似を適用すると推定精度が良くない」ことに気づいたことではないかと想像します。

小標本の場合、標本平均が正規分布に従う保証はありません。

独立２群の差の検定において、正規性の前提が崩れた場合の影響を調べます。

• １群が正規分布に、もう１群がカイ二乗分布に従う場合.
  母平均が同じモデルで平均差の分布を調べます。

• ２群が外れ値を持つ場合
  ２つの正規分布に外れ値を付加し、平均は同じになるように設定します。

いずれも等分散ではないので、ヴェルチのt検定を用います。.

4.1 非正規かつ大標本.

標本数**のカイ二乗分布.

＜帰無仮説＞.

• 群G1（平均5, 標準偏差1）のカイ二乗分布からの大きさ5の標本 \(x_1\)
  
• 群G2（平均5, 標準偏差1）のカイ二乗分布からの大きさ10の標本 \(x_2\)

等分散.

XZ=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,4,4)
  x <- mean(x1)-mean(x2)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dt(x,n1+n1-2),add=T,col="blue")#自由度

検出力は低い.

P値は.

PP=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,4,4)
  x <- mean(x1)-mean(x2)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xp <- x/sqrt(s/n1+s/n2)
  pp <- 1- pt(xp,n1+n2-2)
  PP <- c(PP,pp)  
}
hist(PP,prob=T,breaks = "Scott",main="",xlab="",col=0)

mean(PP)

## [1] 0.1373345

不等分散.

XW=c()
FF=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,8);x2 <- rnorm(n2,8,4)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dt(x,mean(FF)),add=T,col="blue")#自由度

t値の分布は、理論上のt分布や標準正規分布にほぼ一致しており、検定の有効性に問題はなさそうです。

帰無仮説のもとでの分布の当てはまりは良いが、平均差があるときの検出力は良くない。

XW=c()
FF=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,5);x2 <- rnorm(n2,3,5)
  x <- mean(x1)-mean(x2)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
  f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
  FF <- c(FF,f)
}
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dt(x,mean(FF)),add=T,col="blue")#自由度

4.2 非正規かつ小標本.

標本数5のカイ二乗分布.

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- rchisq(n1,5);x2 <- rchisq(n1,5)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

理論上のt分布との大きな乖離はありません。

＜対立仮説＞.

• 群G1（平均5, 標準偏差1）からの大きさ5の標本 \(x_1\)
  
• 群G2（平均4, 標準偏差1）からの大きさ5の標本 \(x_2\)

小標本で、平均差 1 が検出できるでしょうか？

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 5
n2 <- 5
for (i in 1:1000){
  y1 <- rchisq(n1,5);y2 <- rchisq(n2,4)
  y <- mean(y1)-mean(y2)
  sey <- mean(var(y1),var(y2))
  yz <- y/sqrt(sey/n1+sey/n2)
  YZ <- c(YZ,yz)
}
hist(YZ,prob=T,breaks = "Scott",main="",xlab="")
curve(dt(x,(n1+n2-2)),col="blue",add=T)
abline(v = qt(0.95,n1+n2-2),col="red")

検出力は低下しますが、思ったほどではないようです。

4.3 外れ値の影響.

X=c()
Y=c()
XZ=c()
YZ=c()
n1 <- 5
n2 <- 10
for (i in 1:1000){
  x1 <- c(rchisq(n1-1,5),20);x2 <- rchisq(n1,5)
  x <- mean(x1)-mean(x2)
  X <- c(X,x)
  s <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)
  xz <- x/sqrt(s/n1+s/n2)
  XZ <- c(XZ,xz)
}
hist(XZ,prob=T,breaks = "Scott",main="",xlab="",col=0)
curve(dnorm(x, 0, 1),add=T,col="red")
abline(v = qnorm(0.95,0,1),col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
abline(v = qt(0.95,n1+n2-2),col="blue")

XW=c()
XP=c()
n1 <- 4
n2 <- 7
for (i in 1:1000){
  x1 <- c(rnorm(n1-1,5,1),20);x2 <- rnorm(n2,5,5)
  x <- mean(x1)-mean(x2)
  sp <- (var(x1)*(n1-1)+var(x2)*(n2-1))/(n1+n1-2)#プール分散
  xp <- x/sqrt(sp/n1+sp/n2)
  XP <- c(XP,xp)
  xw <- x/sqrt(var(x1)/n1+var(x2)/n2)
  XW <- c(XW,xw)
}
f <- ((var(x1)/n1+var(x2)/n2))^2/((var(x1)^2/(n1^2*(n1-1))+var(x2)^2/(n2^2*(n2-1))))
f

## [1] 4.969163

par(mfrow=c(1,2))
hist(XP,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")
hist(XW,prob=T,breaks = "Scott",main="",xlab="",col=0,xlim=c(-4,4))
curve(dnorm(x, 0, 1),add=T,col="red")
curve(dt(x,(n1+n2-2)),add=T,col="blue")
curve(dt(x,f),add=T,col="green")

5 まとめ.

独立２標本の差の検定では、

• 共通の平均と標準偏差のもとで、出力される平均点差がどのように変動するか分布を調べます. （差がないという帰無仮説のもとでの平均点差の分布です）.
• 平均点に差がある前提のもとで、出力される平均点差がどのように変動するか分布を調べます.
  （一定の差あるという対立仮説のもとでの平均点差の分布です）.

帰無仮説と対立仮説の分布に大きな差があれば、平均点の差を統計的に検出できることになります。

２群の等分散性、正規性、独立性を前提とする。（等分散が仮定できなければ welch検定）
実際には、次の手順で行うことが多い。
　　　\(X_1\sim N(\mu_1,\sigma_1^2)\:,\:X_2\sim N(\mu_2,\sigma_2^2)\::i.i.d.\) 　　は独立な２群を表す確率変数とする。
独立な２群の標本数を \(n_1,n_2\) 、母分散を \(\sigma_1^2,\sigma_2^2\) 、不偏分散を \(s_1^2,s_2^2\) とする。

（1）等分散性の検定を行う。

　　＜帰無仮説＞　\(\sigma_1^2=\sigma_2^2\)　のもとでは、　\(F=\dfrac{s_1^2}{s_2^2}\sim F(n_1-1,n_2-1)\)
　　　　　　　（市原p75によると、不偏分散の大きい方を分子とした片側検定である）

（2）等分散が棄却されなければ、２群の平均に差があるかどうか t 検定を行う。
　　　\(X_1\sim N(\mu_1,\sigma^2)\:,\:X_2\sim N(\mu_2,\sigma^2)\)　 、ただし 　\(\sigma^2=\sigma_1^2=\sigma_2^2\)

　　　\(\overline{X_1}\sim N\left(\mu_1,\dfrac{\sigma^2}{n_1}\right)\:,\:\overline{X_2}\sim N\left(\mu_2,\dfrac{\sigma^2}{n_2}\right)\)　より

　　　　　　　　　　　　　　　\(\overline{X_1}-\overline{X_2}\sim N\left(\mu_1-\mu_2\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)
　 　＜帰無仮説＞　２群の平均に差はないのもとでは、
　　　　　　　　　　　　　　　\(\overline{X_1}-\overline{X_2}\sim N\left(0\:,\:\dfrac{\sigma^2}{n_1}+\dfrac{\sigma^2}{n_2}\right)\)
　　　２群をプールした偏差平方和は 　　\(s_1^2(n_1-1)+s_2^2(n_2-1)\)
　　　自由度は 　　\((n_1-1)+(n_2-1)=n_1+n_2-2\)
　　　２群をプールした不偏分散は 　\(s^2=\dfrac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}\)　となるので
　　　検定量　\(z=\dfrac{\overline{X_1}-\overline{X_2}}{\sigma\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim N(0\:,\:1)\)　　\(t=\dfrac{\overline{X_1}-\overline{X_2}}{s\sqrt{\dfrac{1}{n_1}+\dfrac{1}{n_2}}}\sim t(n_1+n_2-2)\)

　　　母分散は未知なので、自由度 (\((n_1+n_2-2)\) の t検定を行う。（両側検定）