1 混合分布.

連続分布の混合分布について考えます。

A：平均５、標準偏差１の正規乱数３０個.
B：平均５、標準偏差１の正規乱数３０個.

a=rnorm(30,5,1)
A <- rep(a,1)
b=rnorm(30,5,1)
B <- rep(b,1)
U <- c(A,B)
nA <- length(A)
nB <- length(B)
nU <- length(U)
vA <- var(A)*(nA-1)/nA
vB <- var(B)*(nB-1)/nB
vU <- var(U)*(nU-1)/nU#合併集合の分散

　　　vA,vB,vU.

v <- c(vA,vB,vU)
v

## [1] 0.4587298 1.1922383 0.8403222

par(mfrow=c(2,2))
hist(A,prob=T,breaks=c(seq(-0.5,10.5,by=0.5)),col="#ff00ff40",ylim=c(0,1))
hist(B,prob=T,breaks=c(seq(-0.5,10.5,by=0.5)),col="#0000ff40")
hist(U,prob=T,breaks=c(seq(-0.5,10.5,by=0.5)),col="grey88",ylim=c(0,1.5))

1.1 分散の分解.

#各群分散の加重平均
sE <- (nA/nU)*vA+(nB/nU)*vB;sE

## [1] 0.825484

#群間分散の加重平均
sR <-(nA/nU)*(mean(A)-mean(U))^2+(nB/nU)*(mean(B)-mean(U))^2;sR

## [1] 0.01483815

sE+sR

## [1] 0.8403222

vU

## [1] 0.8403222

分散が、各群分散の加重平均と群間分散の加重平均に分解されている。

1.2 平方和分解

#群内平方和
SE <- var(A)*(nA-1)+var(B)*(nB-1)
#群間平方和
SR <-(nA)*(mean(A)-mean(U))^2+(nB)*(mean(B)-mean(U))^2
SE+SR

## [1] 50.41933

#全平方和
vT <- var(U)*(nU-1)
vT

## [1] 50.41933

1.3 群内分散と群間分散.

SE/(nA-1+nB-1)

## [1] 0.853949

SR/1

## [1] 0.8902892

2 分散分析へ.

２群A,B が同じ母集団U からの標本である場合、

V(A),V(B)から母分散を推定することが可能.
群間分散からも母分散を推定することが可能.

　　\(m_A\sim N\left(m,\dfrac{\sigma^2}{n_A}\right)\)　　\((m_A-m)\sim N\left(0,\dfrac{\sigma^2}{n_A}\right)\).

　　\(m_B\sim N\left(m,\dfrac{\sigma^2}{n_B}\right)\)　　\((m_B-m)\sim N\left(0,\dfrac{\sigma^2}{n_B}\right)\).

難しくなりそうなのですが、\(n_A=30,~n_B=30\) の場合について考えてみます。

不偏分散とカイ二乗分布.
偏差平方和を母分散で割ったものは、自由度n-1のカイ二乗分布に従います。

\(\dfrac{\sum_{i=1}^n(x_i-\bar{x})}{\sigma^2}\sim\chi(n-1)\)より、

　　　　　　　\(\dfrac{\sum_{i=1}^n(x_i-\bar{x})}{n-1}\sim \sigma^2\cdot\dfrac{\chi(n-1)}{n-1}\).

自由度n-1のカイ二乗分布の期待値はn-1なので、右辺の期待値は\(\sigma^2\)であり、偏差平方和をn-1で割った値は母分散の不偏推定量となっています。

気を付けること.
カイ二乗分布の形状には歪みがあることを注意しておく必要があります。平均すると\(\sigma^2\)であっても、値の変動が大きいことがありえます。特に、自由度が小さい場合は偏りが大きくなります。

具体的にみていきます。

SE2 = c()
for (i in 1:1000) {
  x2<-rnorm(2,5,1)
  se2 <- var(x2)#群内変動
  SE2 <- c(SE2,se2)
}

SE3 = c()
for (i in 1:1000) {
  x3<-rnorm(10,5,1)
  se3 <- var(x3)#群内変動
  SE3 <- c(SE3,se3)
}

SE5 = c()
for (i in 1:1000) {
  x5<-rnorm(10,5,1)
  se5 <- var(x5)#群内変動
  SE5 <- c(SE5,se5)
}
par(mfrow=c(2,2))
hist(SE2,prob=T,xlim=c(0,10))
hist(SE3,prob=T,xlim=c(0,10))
hist(SE5,prob=T,xlim=c(0,10))

不偏分散の確率分布は、

標本数２では偏りが大きく.
標本数３以上は必要である.

と思います。
この意味するところは、
分散分析では帰無仮説のもとで群平均の分散 \(SR/(a-1)\)を考えますが、群数は３以上が望ましいということです。

2.1 群間分散・群内分散.

東京大学大森宏先生　統計特論２によれば、

２群の分散分析モデル.

　　　\(x1\sim N(\mu+\alpha_1,\sigma^2)~,~~~~x2\sim N(\mu+\alpha_2,\sigma^2)\).

各群の標本数を\(n_1,~n_2\)とすると、
　　　全変動　\(S_T=S_A+S_e\)　（群間変動と群内変動の和）.

ここで、「群間変動\(S_A\)と群内変動\(S_e\)がどのような分布するか」理解しておくことが、分散分析の理解に役立ちます。

\(E(S_A)=n_1\alpha_1^2+n_2\alpha_2^2+\sigma^2\).
\(E(S_e)=(n-2)\sigma^2\).

帰無仮説のもとでは \(\alpha_1=\alpha_2=0\)であるので、

群間分散の期待値：\(E(S_A/1)=\sigma^2\).
群内分散の期待値：\(E(S_e/(n-2))=\sigma^2\).

2.2 分散の挙動.

群間分散：\(S_A/1\).
群内分散：\(S_e/(n-2)\).

これらがどんな挙動をするか具体例で見てみましょう。

2.2.1 群平均に差がないとき.

SE = c()
SR = c()
X=c()
for (i in 1:1000) {
  x1<-rnorm(30,5,1);x2<-rnorm(30,5,1)
  se <- var(x1)*29+var(x2)*29 #群内変動
  SE <- c(SE,se)
  sr <- 30*(mean(x1)-mean(c(x1,x2)))^2+30*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
  SR <- c(SR,sr)
}
par(mfrow=c(1,2))
hist(SR/1,prob=T)
hist(SE/58,prob=T)

v <- c(mean(SR/1),mean(SE/58),sd(SE/58))
v

## [1] 1.0174412 0.9986896 0.1873226

2.2.2 群平均に差があるとき.

SE = c()
SR = c()
X=c()
for (i in 1:1000) {
  x1<-rnorm(30,5,1);x2<-rnorm(30,6,1)
  se <- var(x1)*29+var(x2)*29 #群内変動
  SE <- c(SE,se)
  sr <- 30*(mean(x1)-mean(c(x1,x2)))^2+30*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
  SR <- c(SR,sr)
}
par(mfrow=c(1,2))
hist(SR/1,prob=T)
hist(SE/58,prob=T)

v <- c(mean(SR/1),mean(SE/58),sd(SE/58))
v

## [1] 15.6700990  0.9938099  0.1814883

　　　\(30(0.5^2+0.5^2))+1=15\).

2.2.3 観察のポイント.

群内分散は.

どちらも差がなく、\(1\pm 0.2\) 程度です.

群間分散は.

群平均に差がないとき、大きな値を取ることはほとんどありません.
群平均に差があるとき. 小さな値を取ることはほとんどありません.
（これは標本数が大きいほど、平均差が大きいほど顕著になります）.

根拠：\(E(S_A)=n_1\alpha_1^2+n_2\alpha_2^2+\sigma^2\).

ということは、

２群の差が大きいときには、標本数が小さくても検出可能であり.
２群の差が小さいときには、標本数を大きくしなければならない.

これらのことは、調査をする際に気をつけておかなければなりません。

F検定.
分散分析では、群間分散と群内分散の比であるF分布を用います。
　　　\(F=\dfrac{S_A/1}{S_e/(n-2)}\)　の値が１より有意に大きくなるか判断します。

しかし、いきなりF値を計算しても、どんな意味があるのか分かりづらいので、このようなシミュレーションで十分理解しておくことは有用でしょう。

2.3 F分布.

帰無仮説のもとでは、群間分散と群内分散の比はF分布に従います。
これをシミュレーションで確かめます。

SE = c()
SR = c()
X=c()
for (i in 1:1000) {
  x1<-rnorm(30,5,1);x2<-rnorm(30,5,1)
  se <- var(x1)*29+var(x2)*29 #群内変動
  SE <- c(SE,se)
  sr <- 30*(mean(x1)-mean(c(x1,x2)))^2+30*(mean(x2)-mean(c(x1,x2)))^2 #群間変動
  SR <- c(SR,sr)
  f <- (sr)/(se/58)#F値
  X <- c(X,f)
}
mean(X)

## [1] 1.121241

hist(X,prob=T,breaks = "Scott")
curve(df(x, 1, 58),add=T)

mean(SE/58)

## [1] 1.001361

hist(SE/58)

mean(SR/1)

## [1] 1.082314

hist(SR/1)

cor(SE,SR)

## [1] -0.0006254883

ヒストグラムと理論上のグラフは一致しています。

詳しくは、当サイトF分布を参照してください。

3 探究.

群平均の差を変えると、群間分散がどのように変わるか.
標本数を変えると、群間分散がどのように変わるか.

調べてみましょう。

4 3群の分散分析.

３群で、群平均に差がある場合を下に示します。

SE = c()
SR = c()
X=c()
for (i in 1:1000) {
  x1<-rnorm(30,4,1);x2<-rnorm(30,5,1);x3<-rnorm(30,6,1)
  se <- var(x1)*29+var(x2)*29+var(x3)*29 #群内変動
  SE <- c(SE,se)
  sr <- 30*(mean(x1)-mean(c(x1,x2,x3)))^2+30*(mean(x2)-mean(c(x1,x2,x3)))^2+30*(mean(x3)-mean(c(x1,x2,x3)))^2 #群間変動
  SR <- c(SR,sr)
}

hist(SR/2,prob=T)

mean(SR)

## [1] 62.14034

　　　\(30(1^2+1^2))+2\times 1=62\).

群平均の差を変えると、群間分散がどのように変化するか探究しましょう.

混合分布から分散分析へ