前回書いた結果をもとにして、ちゃんと近似できているのか試してみることにしました。
今回は Welch (1938) の (2) 式の統計量
\[
v=\frac{\bar{X}_1-\bar{X}_2}{\sqrt{\frac{V_1}{n_1}+\frac{V_2}{n_2}}}
\] を導出する過程で出てくる分布の近似を考えてみようと思います。
この統計量は、
\[
v=\frac{\frac{\bar{X}_1-\bar{X}_2}{\sqrt{\frac{{\s}_1^2}{n_1}+\frac{{\s}_2^2}{n_2}}}}{\sqrt{\left(\nu\frac{\frac{V_1}{n_1}+\frac{V_2}{n_2}}{\frac{{\s}_1^2}{n_1}+\frac{{\s}_2^2}{n_2}}\right)/\nu}}
\] この形に変形することが出来ます。
$\bar{X}_1$ と $\bar{X}_2$ が、それぞれ $\mathrm{N}(\mu_1, \s_1^2/n_1), \mathrm{N}(\mu_2,\s_2^2/n_2)$ に従うとき、$\bar{X}_1-\bar{X}_2$ は $\mathrm{N}(\mu_1-\mu_2,\s_1^2/n_1+\s_2^2/n_2)$ に従います。
それを標準偏差で割っているので、分子は $\mathrm{N}(\mu_1-\mu_2,1^2)$ に従います。
では、分母はどうなのでしょう？という問いの答えのひとつが Satterthwaite の近似です。
Satterthwaite の近似は、カイ二乗分布に従う確率変数の線形結合の分布を近似する方法なので、このヘンテコな分母を近似することができます。
ここで、$X\sim\mathrm{N}(\mu,1^2)$、$Y\sim\chi^2(\gamma)$ のとき、$X/\sqrt{Y/\gamma}\sim t(\mu, \gamma)$ となることを思い出すと、
\[
W=\left(\nu\frac{\frac{V_1}{n_1}+\frac{V_2}{n_2}}{\frac{{\s}_1^2}{n_1}+\frac{{\s}_2^2}{n_2}}\right)=\frac{c}{n_1}V_1+\frac{c}{n_2}V_2
\] が、自由度 $\nu$ のカイ二乗分布に従えば、統計量 $v$ が非心度 $\mu$、自由度 $\nu$ の非心 $t$ 分布に従うと考えられます。
一般の場合の Satterthwaite の近似を使うと、$W$ の分布はカイ二乗分布で近似でき、その自由度は近似的に、
\[
\hat{\nu}=\frac{\left(c \sum_{i=1}^2 \frac{V_i}{n_i} \right)^2}{c^2 \sum_{i=1}^2 \frac{V_i^2}{n_i^2 \nu_i}}=\frac{\left(\sum_{i=1}^2 \frac{V_i}{n_i} \right)^2}{\sum_{i=1}^2 \frac{V_i^2}{n_i^2 \nu_i}}
\] で求めることが出来ます。
$W$ に含まれていた真値 $\nu$ は、打ち消されています。
この統計量 $v$ が、帰無仮説 $H_0: \mu=0$ のもとで近似的に自由度 $\nu$ の $t$ 分布に従うとして検定を行うのが、Welch の $t$ 検定です。
今回は、$W$ の分布を本当に近似できているのかをシミュレーションで確認してみることにしました。
ただし、$W$ 自身は未知の真値 $\nu$ を含むため、これを打ち消した $W/\nu$ の分布について考えてみることにしました。
# パスの設定と関数の読み込み if (is.na(commandArgs()[6])) { path <- "C:/"; } else { path <- commandArgs()[6]; } tapplySSD <- function(x, index) { # create index indicate matrix index <- as.numeric(factor(index)) mindex <- NULL for(i in 1:length(unique(index))) { mindex <- cbind(mindex, index-(i-1)) } mindex[mindex!=1]<-0 # Sum of squared deviation ssd <- t(mindex) %*% x^2 - (t(mindex) %*% x)^2 / as.numeric(table(index)) return(ssd) } # 条件設定とシミュレーションデータの生成 set.seed(31035963) n <- c(20, 10) sd <- c( 1, 4) mean <- c(1, 1) rep <- 1000000 v_denom <- sum(sd^2 / n) cx <- do.call("rbind", mapply(function(n, mean, sd, rep) matrix(rnorm(rep * n, mean = mean, sd = sd), nrow=n, ncol=rep), n, mean, sd, rep, SIMPLIFY = FALSE) ) g <- rep(1:length(n), n) var <- tapplySSD(cx, g) / (n - 1) W_nu <- as.vector((rep(1, length(n)) %*% (var / n)) / v_denom) # Satterthwaiteの近似 sa_est <- (rep(1, length(n)) %*% (var / n))^2 / (rep(1, length(n)) %*% (var^2 / n^2 / (n-1))) # Satterthwaiteの近似式にsigma_iの真値を代入してみたもの sa_true <- (sum(sd^2 / n))^2 / (sum(sd^4 / n^2 / (n-1))) # nlm関数で nu の最尤推定値を求める mlogL <- function(df, x) { -sum(-lgamma(df/2) + (df/2) * log(df/2) + (df/2-1) * log(x) - df * x / 2) } df.start <- sa_true out <- nlm(mlogL, df.start, x = W_nu, steptol = 1e-8, gradtol = 1e-8) # 結果出力 cat("Simulation MLE =", sprintf("%4.2f", out$estimate), "\n") cat("Simulation 5% point =", sprintf("%4.2f", quantile(sa_est, 0.05)), "\n") cat("Simulation 95% point =", sprintf("%4.2f", quantile(sa_est, 0.95)), "\n") # グラフ出力 postscript( file=paste(path, "sim.eps", sep=""), width = 9, height = 9, pointsize = 25, horizontal = FALSE, onefile = FALSE, paper = "special", family="Palatino" ) density.cv <- function(df, x) { exp(-lgamma(df/2) + (df/2) * log(df/2) + (df/2-1) * log(x) - df * x / 2) } par(bty="l") par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c(2.0, 0.7, 0)) h <- hist( W_nu, probability=T, main=NA, breaks=100, xlab=expression(italic(W/nu)), ylim=c(0,1.1) ) l <- seq(min(h$breaks), max(h$breaks), 0.01) par(new=T) plot( l, density.cv(out$estimate, l), type="l", ylim=c(0,1.1), xlab=NA, ylab=NA, xaxt="n", yaxt="n", col="#FAA55C", lwd=3 ) par(new=T) plot( l, density.cv(quantile(sa_est, 0.05), l), type="l", ylim=c(0,1.1), xlab=NA, ylab=NA, xaxt="n", yaxt="n", col="#5ECD22", lwd=3 ) par(new=T) plot( l, density.cv(quantile(sa_est, 0.95), l), type="l", ylim=c(0,1.1), xlab=NA, ylab=NA, xaxt="n", yaxt="n", col="#E64B6B", lwd=3 ) legend( max(h$breaks), 1.1, c("Simulation (MLE)", "Simulation 5% point (SA)", "Simulation 95% point (SA)"), lty=c(1,1,1), lwd=c(3,3,3), col=c("#FAA55C", "#5ECD22", "#E64B6B"), cex=.65, adj = c(0, .5), xjust=1, yjust=1, x.intersp = 0.5, y.intersp = 0.9 ) dev.off()
10万組ほどデータを生成するので、結構時間がかかります。
試してみたい奇特な人がいたら、適当に調整してください。
高速な関数に取り替えました。
MLEは、10万組すべてを使って $W/\nu$ の分布を、$W$ がカイ二乗分布に従うと仮定して最尤推定した結果です。
SA は、10万組それぞれに Satterthwaite の近似を行ったものから、5% 点と 95% 点を求めた結果です。
大雑把にみても、5% 点や 95% 点でもいい感じの近似になっている事が分かります。
