10月8日（水）

検出力について

数理統計の教科書を読んでいたら、検定における「検出力」に関する内容があって、すぐにはピンと来なかったので、Ｒで検出力を求める関数をつくってみた。
（ここへ来てより道しすぎな感は否めないが…。）

正規母集団 N(mu,1)からn個の無作為標本をとってくる場合を考える。
いま、mu>mu0のとき、mu=mu0を帰無仮説として有意水準95％で片側検定を行って帰無仮説が棄却される確率を求める。

とりあえずモンテカルロ法が使えると思ったのでそのバージョン。

mu0 <- 0

ken.p <- function(x){

a <- qnorm(0.975,mu0,1/x[2])

cnt <- 0
for (i in 1:10000){
if (rnorm(1,x[1],1/x[2])>a) cnt <- cnt+1
else cnt <- cnt
}
return(cnt/10000)
}

mu <- seq(0,1,length=50)
n <- 1:10

mu.n <-  as.matrix(expand.grid(mu,n))

ken <- apply(mu.n,1,ken.p)

library(rgl)

plot3d(mu.n[,1],mu.n[,2],ken,xlab="mu",ylab="n")

mu,nが大きくなるほど検出力が大きくなることが確認できた。
そういえばモンテカルロ・シミュレーションを使ったのは初めてなので、感慨深い。

ただ、モンテカルロを使う必要はなくて、

ken.p. <- function(x){

a <- qnorm(0.975,mu0,1/x[2])

b <- pnorm(a,x[1],1/x[2],lower=F)

return(b)
}

としても可能。

とりあえず検出力の概念については理解できた。

要するに、標準正規分布の表しか与えられていなくて紙と鉛筆だけでやる場合（教科書のやり方）、式を同値変形して解かなくてはいかないのでイメージとしてつかみにくかったということらしい。
確率分布とかについて考える場合、直観的に処理できるという意味でもコンピュータは優れていると感じた。

以上