正確検定への近似法―漸近検定―

　これまで紹介してきたMann-WhitneyのＵ検定などのノンパラメトリック検定（そのほかの例1、例2、例3、例4）では、しばしば正確検定Exact testとか漸近検定Asymptotic testと呼ばれる方法が登場してきた。基本的な概念としては、条件を満たさないとき（大体はサンプルサイズが小さいとき）は組み合わせの数え上げによって正確な帰無分布を算出し、それを使ってp値を計算する方法が正確検定である。一方、条件を満たしたときに、正確な帰無分布が別のより簡単に計算可能な確率分布に近似できることを利用して、p値を計算する方法が漸近検定となる。

　このような概念とノンパラメトリック検定が関連づいている理由は、ノンパラメトリック検定が計測値の順位の並べ替えとそれによって計算される検定統計量によって帰無分布を作っているからである。組み合わせさえ分かれば、全場合の数と特定の条件を満たす場合の数の割合から、確率を計算できる。正確なp値を算出するという点から、一見すれば正確検定のほうがより優れている、という印象を抱くかもしれない。しかし、正確検定は組み合わせ計算に依存しているので、サンプルサイズが大きくなれば計算量が膨大になり、検定を終えるまでに時間がかかる、ひいては有限時間に終わらないという事態が起こる。漸近検定は、そのような正確検定の弱点を補う手法である。正確な帰無分布を計算せずとも、サンプルサイズが大きいときには標本から計算できる少数のパラメータ（平均や分散など）を使って、他の分布に帰着することで、大幅に計算時間を軽減できるのである。例えば、正規分布に帰着できれば、標本の平均と分散だけで確率分布を表現できる。正確検定は純朴かつ、確実な方法であるのだが、その裏には膨大な計算が潜んでおり、必ずしもより少数のパラメータで帰無分布を表そうとする漸近検定より優れている、というわけではない。本稿では二項検定を例に、正確検定と漸近検定を解説する。

二項分布と二項検定

　以前、二項検定については紹介した。二項検定は、標本が二者択一の状態を示すときに、注目している現象の出現割合が、指定した確率と異なるかを判定する検定である。二項検定における検定統計量は、試行回数と注目している現象が出現した回数（成功回数）である。試行回数と期待される成功確率がわかれば、帰無分布を作ることができて、その分布は二項分布と呼ばれている。この成功回数が帰無分布のどの位置に存在するかで、p値を計算して、検定を行うわけである。この帰無分布は組み合わせの数え上げなので、理論上、「簡単に」計算することができる。

　例えば、コイントスを10回行うことを考える。まず、コイン10枚の表と裏の出現の場合の数は、2^10 = 1024通りである。この10枚のコインのうち、表の枚数が0枚の場合の数は、1通りなので、表の枚数を成功回数とすると表0枚は1/1024となる。表1枚なら、10枚のコインのうちどれか1枚のコインが表なので、表にするコインは10通り選べる。なので、表1枚は10/1024となる。10枚のコインのうち、どれか2枚を表にする選び方は(10×9)/(2×1) = 45通りなので、表2枚は45/1024……という風に、表10枚までそれぞれの確率は計算可能である。これは表と裏が1/2の確率で出現する仮定で成り立っており、表になる枚数をkとしたときにその確率は[10!/{k!(10 - k)!}] ×{(1/2)^k}{(1/2)^(10 - k)}と一致する。このk = 0～10までを分布としてまとめたのが二項分布であり、帰無分布である。実際の分布の形は、次のようになる。

------------------------------------------------------

library(plotn)

pos <- barplotn(dbinom(0:10, size = 10, prob = 1/2), 

                xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率", space = 0)

overdraw(axis(1, at = pos, labels = 0:10, cex.axis = 1.2)) #図1の描画

------------------------------------------------------

表裏が平等に出現するコインを10回投げればこのような帰無分布になることが期待されるわけで、10回投げて1回も表が出なかったコインがあった時に、いかさまコインかを判定するのに使える。二項検定では、10回投げて1回も表が出なかったとき、p値として1/1024 + 1/1024 = 1/512を算出する（計算法は両側検定と片側検定参照）。この場合、有意水準5%で帰無仮説、つまり表の出現確率1/2が棄却されるので、このコインはいかさまコインである可能性があるわけだ。このp値の算出法は、場合の数の数え上げであり、正確検定にあたる。

二項分布の正規近似

　ところが、試行回数（サンプルサイズ）が多くなると、非常に大きな値を扱う必要が出てきて、確率計算は計算式を書き下すことはできても、値を出すのは大変になってくる。試行回数が大きいとき、試行回数をn、成功確率をpとすると、二項分布は平均np、分散np(1 - p)の正規分布で近似できることが知られている。例えば、試行回数20、成功確率1/3なら、平均20/3、分散40/9の正規分布とほぼ同じになる、ということである。Rでは次のように描画する。

------------------------------------------------------

total <- 20

p <- 1/3

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab), max(lab), length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率", space = 0) #図2の描画

xx <- seq(min(pos), max(pos), length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

確かに図2では、棒グラフで示される二項分布と赤線で示される正規分布が近い値をとっているように見える。繰り替えずが、この正規近似は試行回数が大きいときに有効であり、条件を満たさないときは近似精度が良くない。例えば、試行回数が小さいとき、他にも同じ試行回数であっても成功確率が1/2から大きく外れるときには、あまり近似がうまくいかない。

------------------------------------------------------

total <- 10

p <- 1/3

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab), max(lab), length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率", space = 0) #図3の描画

xx <- seq(min(pos), max(pos), length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

total <- 20

p <- 1/10

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab), max(lab), length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率", space = 0) #図4の描画

xx <- seq(min(pos), max(pos), length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

上記の場合では、特に試行回数0周りでの近似精度が悪くなっている。では、試行回数や成功確率は、どんな時に近似精度が良くなるのか、近似精度の程度はどのように変化するのかを見てみることにしよう。今次のような図5を考える。

------------------------------------------------------

total <- 10

p <- 1/10

xx <- 0:total

yy1 <- dbinom(xx, size = total, prob = p)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

d <- data.frame(x = xx, binom = yy1, norm = yy2)

plotn(d[,1],d[,2:3],mode = "m", xlab = "成功回数k", ylab = "確率",

      line = T, legend = T, lty = 2, leg.sp = 5, leg.lab = c("二項","正規")) #図5の描画

for(i in 1:nrow(d)){

  overdraw(

    segments(d[i,1], d[i,2], d[i,1], d[i,3])

           )

}

------------------------------------------------------

図5において、黒線で示した部分が、二項分布と正規分布の分布の違いを表す偏差である（実際は、二項分布と正規分布は確率計算の性質が違っており、単純比較できないが、「見た目」の比較としては十分である）。この偏差の平方和を分布間の違いを示す値とすれば、偏差平方和が小さいほどよく近似が成り立っている、と判断することができる。次のように、まずは成功確率を連続的に変化させて、どの点でもっともよりよく近似が成り立つかを探索する。

------------------------------------------------------

SqSum <- NULL

for(j in c(10, 20, 50, 100)){

  total <- j

  tmp <- NULL

  for(i in seq(0,1,length = 200)){

    p <- i

    bin <- dbinom(0:total, size = total, prob = p)

    nom <- dnorm(0:total, mean = total * p, sd = sqrt(total * p * (1 - p)))

    tmp <- c(tmp, sum((bin - nom)^2))

  }

  SqSum <- cbind(SqSum, tmp)

}

SqSum <- as.data.frame(SqSum)

colnames(SqSum) <- paste0("Total_", c(10, 20, 50, 100))

plotn(seq(0,1,length = 200), SqSum, mode = "m", 

      ylim = c(0, 0.05), xlab = "確率p", ylab = "偏差平方和",

      line = T, leg.sp = 7, pch = "", legend = T) #図6の描画

------------------------------------------------------

図6からまず、成功確率が1/2に近づくほど、偏差平方和が小さくなる、つまり近似精度が良くなることがわかる。また、試行回数が大きくなれば、成功確率が1/2から多少離れてていても、近似精度が保たれることもわかる。次は、試行回数を連続的に変化させた場合を考えてみる。

------------------------------------------------------

SqSum <- NULL

for(j in c(1/10, 1/5, 1/3, 1/2)){

  p <- j

  tmp <- NULL

  for(i in 1:100){

    total <- i

    bin <- dbinom(0:total, size = total, prob = p)

    nom <- dnorm(0:total, mean = total * p, sd = sqrt(total * p * (1 - p)))

    tmp <- c(tmp, sum((bin - nom)^2))

  }

  SqSum <- cbind(SqSum, tmp)

}

SqSum <- as.data.frame(SqSum)

colnames(SqSum) <- paste0("Prob_", c("1/10", "1/5", "1/3", "1/2"))

plotn(1:100, SqSum, mode = "m", 

      ylim = c(0, 0.05), xlab = "確率p", ylab = "偏差平方和",

      line = T, leg.sp = 7, pch = "", legend = T) #図6の描画

------------------------------------------------------

図7から、試行回数が大きくなるほど、近似精度が良くなることがわかり、成功確率が1/2から離れているほど、近似精度のよくなり方が遅いことがわかる。

　以上のように、試行回数と成功確率によって、近似精度が変わることがわかったと思う。二項分布は性質上、0から試行回数nまでの分布である。試行回数が小さいときや、成功確率が1/2から離れているとき、というのは、二項分布が0付近、あるいはn付近に分布している、ということを表しており、図3、4でも見た結果と同じく、正規分布とのずれが大きくなっているのである。

　二項検定を行うとき、帰無仮説として期待される成功確率を指定する。では、成功確率を決定した時、どれくらいの試行回数があれば正規近似をもとにp値を計算をしてもよいか、できればそれを簡単に算出できないかを考察しよう。成功確率が1/2から離れており、試行回数も小さいときは（図3、4）、二項分布は0付近や、n付近に分布している。そのまま正規近似すると精度が悪いだけでなく、正規分布は0より小さい値、nより大きい値をとるため、その分だけ正規近似した時にp値の計算に誤差が生じてしまう。そこで、試行回数を大きくすることで、二項分布を0付近や、n付近から離し、正規近似した時に、0以下やn以上の値が計算される割合を小さくすることを考える。

つまり、上図において、正規近似した分布が、0以下やn以上の領域にはみ出た部分の面積が小さくなるように、試行回数を指定する、ということである。正規分布は-∞から∞まで定義されているので、はみ出した部分の面積を0にすることはできないが、p値の計算の要領と同じく、全体の面積の何%以下にする、ということは可能である。二項分布の近似する正規分布の平均はnp、分散はnp(1 - p)である。ここで、正規分布の性質を思い出すと、平均±標準偏差の範囲の面積が約68%、平均±2×標準偏差なら95%、平均±3×標準偏差なら99%であった。そこで平均が、√10×標準偏差以上、0やnから離れているような条件を考えてみる（√10に必ずしもする必要はないが、あとの計算が簡単になるため）。√10 > 3なので、平均±√10×標準偏差の範囲は99.843%の面積に相当し、0以下やn以上にはみ出す部分の面積は(1 - 0.99843)/2 = 0.00078、0.1%以下である。念のため、以下のように計算して確認する。

------------------------------------------------------

1 - 2 * pnorm(sqrt(10), mean = 0, sd = 1, lower.tail = F) #√10×標準偏差の範囲の面積

## [1] 0.9984346

pnorm(sqrt(10), mean = 0, sd = 1, lower.tail = F) #√10×標準偏差の範囲外の片側面積

## [1] 0.0007827011

1 - 2 * pnorm(3, mean = 0, sd = 1, lower.tail = F) #3×標準偏差の範囲の面積

## [1] 0.9973002

1 - 2 * pnorm(2, mean = 0, sd = 1, lower.tail = F) #2×標準偏差の範囲の面積

## [1] 0.9544997

1 - 2 * pnorm(1, mean = 0, sd = 1, lower.tail = F) #標準偏差の範囲の面積

## [1] 0.6826895

------------------------------------------------------

下記のように、nとpの関係を計算していく。

p ≦ 1/2のときは、二項分布は0に近いので平均 - √10×標準偏差 > 0となるように、p > 1/2のときは、二項分布はnに近いので、平均 + √10×標準偏差 < nとなるようなnを求める。最終的な計算結果は、p ≦ 1/2でも、p > 1/2でも似たような結果になり、p ≦ 1/2なら n > 10 (1 - p)/p ≧ 10 p/(1 - p)、p > 1/2なら n > 10 p/(1 - p) ≧ 10 (1 - p)/pとなるので、pが1/2より大きいか小さいかに関わらず、試行回数nは、p/(1 - p)と(1 - p)/pのうち、より大きい方の値に10を掛けた値より大きく設定すればよいことがわかる。ここで、p/(1 - p)あるいは(1 - p)/pはオッズ比Odds ratioと呼ばれている。オッズ比は現象の出現比ともいえる。例えば、コインのように表1/2、裏1/2で出現するなら、オッズ比は(1/2)/(1 - 1/2) = 1であり、表と裏が平等に出現する様子を表している。サイコロで1の出目に注目し、1とそれ以外の出目、と考えるなら、平等なサイコロであれば(1の出る確率)/(他の出目が出る確率) = (1/6)/(1 - 1/6) = 1/5となり、1の出目より他の出目のほうが5倍出やすい、様子を示している。このオッズ比の概念がわかっていれば、帰無仮説の成功確率に関して、オッズ比とその逆数を求めて、より大きな方に10を掛けた値よりも大きい試行回数をおこなえば、正規近似でp値を求めてもよいことになる（以降、オッズ比とその逆数の最大値を断りなくオッズ比ということにする）。

　この関係式は、これまでの議論と同じく、p = 1/2のときに最も小さな試行回数で正規近似が成り立つことを示している。 p = 1/2であれば、オッズ比=1となり、その逆数も1となる。p = 1/3となれば、オッズ比=1/2、その逆数は2であり、採用されるのは2のほうとなる。このようにpが1/2から離れれば、オッズ比は1より大きな値を示すようになり、その分だけ正規近似を成り立たせるために、試行回数が必要になるのである。また、この関係式は、もし正規近似で定義域外にはみ出す面積を0.1%以下に抑えたいなら、試行回数10以下では不可能であることを示している。そのような小標本であれば、正規近似に頼らず、正確検定を行うべき、ということだ。

　図1～4の状況で、それぞれが正規近似を行う上で、妥当な試行回数となっているかを判断してみよう。図1であれば、p = 1/2でオッズ比 = 1、なので試行回数は10×1 = 10以上は欲しい。試行回数10回は正規近似するうえでギリギリの値だ。図2はp = 1/3でオッズ比 = 2なので、試行回数は20回欲しい。これも試行回数20回はギリギリ、ということだ。図3もオッズ比 = 2なので、試行回数10回では正規近似は難しい。正確検定を行おう。図4はオッズ比 = 9となり、試行回数20回では、正規近似は無理だ。試しに、試行回数を91回にして帰無分布を描いてみよう。

------------------------------------------------------

total <- 91

p <- 1/10

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab), max(lab), length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率", space = 0) #図8の描画

xx <- seq(min(pos), max(pos), length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

確かに、これなら問題なく近似ができそうである。このように、帰無仮説で期待される成功確率から、正規近似に必要な試行回数を計算することができた。

Rで行う二項検定の正確検定と漸近検定

　では、実際にRの二項検定に関して、正確検定と漸近検定を実行してみよう。まずはコイントスを例にとり、データを見てみる。

------------------------------------------------------

#Coin

xx <- c("F", "F", "F", "B", "F", "F", "B", "F")

yy <- 1:8

d <- data.frame(try = yy, coin = xx)

d

##   try coin

## 1   1    F

## 2   2    F

## 3   3    F

## 4   4    B

## 5   5    F

## 6   6    F

## 7   7    B

## 8   8    F

table(d$coin)

## 

## B F 

## 2 6

------------------------------------------------------

まず、試行回数は8回である。そして、表が出た回数に着目すると、6回である。このような標本が得られた時、二項検定の正確検定と漸近検定を行い、p値にどれくらい違いが出るかを見てみることにする。帰無仮説はコインなので表が出る確率は1/2であると期待される。次のようなコマンドを実行する。

------------------------------------------------------

total <- sum(nrow(d))

front <- table(d$coin)[2]

p <- 1/2

binom.test(x = front, n = total, p = p) #正確検定

## 

##  Exact binomial test

## 

## data:  front and total

## number of successes = 6, number of trials = 8, p-value = 0.2891

## alternative hypothesis: true probability of success is not equal to 0.5

## 95 percent confidence interval:

##  0.3491442 0.9681460

## sample estimates:

## probability of success 

##                   0.75

------------------------------------------------------

正確検定時、p = 0.2891であり、帰無仮説は棄却できない。では、もし、二項分布を正規分布で近似できると考えてp値を算出したらどうなるだろうか。

------------------------------------------------------

if(front > total/2){ #表の数が期待値（試行回数×成功確率）より大きいかどうかで、計算する領域を変える。

  lower <- F

  correct <- -1/2 #連続値補正

} else {

  lower <- T

  correct <- 1/2 #連続値補正

}

pnorm(front + correct, mean = total * 1/2, sd = sqrt(total * 1/2 * 1/2), lower.tail = lower) * 2 #連続値補正ありの漸近検定

##         F 

## 0.2888444

------------------------------------------------------

連続値補正とは、離散分布を連続分布に近似するときに行う補正で、私の経験上、試行回数が小さく、成功確率 が1/2に近いときに行うと、二項検定の正規近似の精度が上がる。試行回数が大きくなるほど、連続値補正は重要な補正ではなくなる。pnorm()で、近似される正規分布からp値を計算できて、p = 0.2888となった。試行回数10以下であるが、いい近似精度である。しかし、成功確率が1/2から離れると途端に近似精度が悪くなるので注意が必要である。では、どんな領域がp値として計算されているか図示してみよう。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab) - 1, max(lab) + 1, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(front, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図9の描画

dif <- pos[2] - pos[1]

if(front > total - front){

  minor <- total - front

} else {

  minor <- front

}

axx1 <- seq(min(lab) - 1, minor + 1/2, length = 200)

ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx1 <- c(axx1, axx1[length(axx1)])

ayy1 <- c(ayy1, 0)

axx2 <- seq(total - minor - 1/2, max(lab) + 1, length = 200)

ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx2 <- c(axx2[1], axx2)

ayy2 <- c(0, ayy2)

xx <- seq(min(pos) - dif, max(pos) + dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = pos[total + 1 - minor] - dif/2 - min(lab), lty = 2),

         abline(v = pos[minor + 1] + dif/2 - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

図9を見ると、成功確率1/2の時は、二項分布が左右対称となり、正規分布との近似精度が非常に良いので、正確検定と漸近検定において、（連続値補正に注意すれば）p値に大きな差がないのだと考えられる。では、試行回数を増やしたときの挙動を見てみよう。試行回数10回の標本だと次のようになる。

------------------------------------------------------

xx <- c("F", "B", "B", "F", "B", "F", "B", "F", "B", "B")

yy <- 1:10

d <- data.frame(try = yy, coin = xx)

table(d$coin)

## 

## B F 

## 6 4

total <- sum(nrow(d))

front <- table(d$coin)[2]

p <- 1/2

binom.test(x = front, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  front and total

## number of successes = 4, number of trials = 10, p-value = 0.7539

## alternative hypothesis: true probability of success is not equal to 0.5

## 95 percent confidence interval:

##  0.1215523 0.7376219

## sample estimates:

## probability of success 

##                    0.4

if(front > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(front + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##         F 

## 0.7518296

------------------------------------------------------

正確検定ではp = 0.7539、漸近検定ではp = 0.7518であり、漸近検定の精度は良いようだ。同じく、図示を行う。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab) - 1, max(lab) + 1, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(front, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図10の描画

dif <- pos[2] - pos[1]

if(front > total - front){

  minor <- total - front

} else {

  minor <- front

}

axx1 <- seq(min(lab) - 1, minor + 1/2, length = 200)

ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx1 <- c(axx1, axx1[length(axx1)])

ayy1 <- c(ayy1, 0)

axx2 <- seq(total - minor - 1/2, max(lab) + 1, length = 200)

ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx2 <- c(axx2[1], axx2)

ayy2 <- c(0, ayy2)

xx <- seq(min(pos) - dif, max(pos) + dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = pos[total + 1 - minor] - dif/2 - min(lab), lty = 2),

         abline(v = pos[minor + 1] + dif/2 - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

コインの例の最後は、試行回数40回の標本を見てみる。

------------------------------------------------------

xx <- c("B", "B", "B", "B", "F", "B", "B", "F", "B", "B", 

        "B", "B", "F", "B", "F", "B", "B", "F", "B", "F",

        "B", "F", "B", "F", "B", "B", "B", "F", "F", "B", 

        "F", "B", "F", "B", "B", "F", "F", "B", "F", "F")

yy <- 1:40

d <- data.frame(try = yy, coin = xx)

table(d$coin)

## 

##  B  F 

## 24 16

total <- sum(nrow(d))

front <- table(d$coin)[2]

p <- 1/2

binom.test(x = front, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  front and total

## number of successes = 16, number of trials = 40, p-value = 0.2682

## alternative hypothesis: true probability of success is not equal to 0.5

## 95 percent confidence interval:

##  0.2486500 0.5667329

## sample estimates:

## probability of success 

##                    0.4

if(front > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(front + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##         F 

## 0.2683816

------------------------------------------------------

正確検定ではp = 0.2682、漸近検定ではp = 0.2683であり、漸近検定の精度は良いようだ。試行回数が増えるほど、徐々に近似精度が上昇していることがわかるだろう。同じく、図示を行う。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

xx <- seq(min(lab) - 1, max(lab) + 1, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(front, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図11の描画

dif <- pos[2] - pos[1]

if(front > total - front){

  minor <- total - front

} else {

  minor <- front

}

axx1 <- seq(min(lab) - 1, minor + 1/2, length = 200)

ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx1 <- c(axx1, axx1[length(axx1)])

ayy1 <- c(ayy1, 0)

axx2 <- seq(total - minor - 1/2, max(lab) + 1, length = 200)

ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

axx2 <- c(axx2[1], axx2)

ayy2 <- c(0, ayy2)

xx <- seq(min(pos) - dif, max(pos) + dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = pos[total + 1 - minor] - dif/2 - min(lab), lty = 2),

         abline(v = pos[minor + 1] + dif/2 - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

　次は、サイコロの例を考える。今、1の出目の出現率に注目し、1が出ればTRUE、それ以外ならFALSEという値をとして記録した。試行回数20として、二項検定を行おう。このとき、いかさまのないサイコロであれば1の出現率は1/6である。

------------------------------------------------------

xx <- c("F", "T", "F", "F", "T", "F", "T", "F", "T", "F", 

        "T", "F", "F", "T", "T", "F", "F", "F", "F", "F")

yy <- 1:20

d <- data.frame(try = yy, dice = xx)

table(d$dice)

## 

##  F  T 

## 13  7

total <- sum(nrow(d))

One <- table(d$dice)[2]

p <- 1/6

binom.test(x = One, n = total, p = p) #正確検定

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 7, number of trials = 20, p-value = 0.03714

## alternative hypothesis: true probability of success is not equal to 0.1666667

## 95 percent confidence interval:

##  0.1539092 0.5921885

## sample estimates:

## probability of success 

##                   0.35

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2 #連続値補正ありの漸近検定

##          T 

## 0.05743312

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2 #連続値補正なしの漸近検定

##         T 

## 0.0278069

------------------------------------------------------

正確検定ではp = 0.03714である。漸近検定では、連続値補正ありではp = 0.05743、なしではp = 0.02781である。むしろ連続値補正がない方が、よりp値は正確検定の時に近い。それゆえ、成功確率が1/2ではないときは、連続値補正が必ずしも有効に働くとは限らない。今回の場合、正確検定と連続値補正なしの漸近検定では、有意水準5％で帰無仮説を棄却できるが、連続値補正をした場合は棄却できない。このように、補正の有無によって判断が変わりうるので、特に試行回数が小さく、成功確率がゆがんでいるときの漸近検定は注意しよう。試行回数が少ないときは、正確検定が無難である。この時のp値として計算される領域を図示しよう。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

min_c <- 2.5

max_c <- 1

xx <- seq(min(lab) - min_c, max(lab) + max_c, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(One, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図12の描画

dif <- pos[2] - pos[1]

if(One > total*p){

  axx2 <- seq(One, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  s <- abs(One - total * p)

  axx1 <- seq(min(lab) - min_c, total * p - s, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  lwr <- total * p - s + 1/2

  upr <- One + 1/2

} else {

  axx1 <- seq(min(lab) - min_c, One, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  s <- abs(One - total * p)

  axx2 <- seq(total * p + s, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  lwr <- One + 1/2

  upr <- total * p + s + 1/2

}

xx <- seq(min(pos) - min_c*dif, max(pos) + max_c*dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = lwr - min(lab), lty = 2),

         abline(v = upr - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

次は試行回数40の時を考える。

------------------------------------------------------

xx <- c("T", "F", "F", "F", "T", "F", "T", "F", "T", "F", 

        "F", "F", "F", "F", "F", "F", "T", "F", "F", "F",

        "T", "F", "F", "T", "T", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "T", "F", "F", "T", "F", "F")

yy <- 1:40

d <- data.frame(try = yy, dice = xx)

table(d$dice)

## 

##  F  T 

## 30 10

total <- sum(nrow(d))

One <- table(d$dice)[2]

p <- 1/6

binom.test(x = One, n = total, p = p) #正確検定

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 10, number of trials = 40, p-value = 0.1985

## alternative hypothesis: true probability of success is not equal to 0.1666667

## 95 percent confidence interval:

##  0.1269148 0.4119620

## sample estimates:

## probability of success 

##                   0.25

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2 #連続値補正ありの漸近検定

##         T 

## 0.2293319

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2 #連続値補正なしの漸近検定

##         T 

## 0.1572992

------------------------------------------------------

正確検定ではp = 0.1985である。漸近検定では、連続値補正ありではp = 0.2293、なしではp = 0.1572である。今回は、連続値補正ありのほうがわずかに精度が良い。図示も行おう。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

min_c <- 3

max_c <- 1

xx <- seq(min(lab) - min_c, max(lab) + max_c, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(One, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図13の描画

dif <- pos[2] - pos[1]

if(One > total*p){

  axx2 <- seq(One, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  s <- abs(One - total * p)

  axx1 <- seq(min(lab) - min_c, total * p - s, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  lwr <- total * p - s + 1/2

  upr <- One + 1/2

} else {

  axx1 <- seq(min(lab) - min_c, One, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  s <- abs(One - total * p)

  axx2 <- seq(total * p + s, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  lwr <- One + 1/2

  upr <- total * p + s + 1/2

}

xx <- seq(min(pos) - min_c*dif, max(pos) + max_c*dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = lwr - min(lab), lty = 2),

         abline(v = upr - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

サイコロの例の最後は試行回数60の時を考える。

------------------------------------------------------

xx <- c("F", "F", "F", "F", "T", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "T", "F", "F", "F",

        "T", "F", "F", "F", "T", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "F", "F", "F", "F",

        "F", "T", "F", "T", "F", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "F", "F", "F", "F")

yy <- 1:60

d <- data.frame(try = yy, dice = xx)

table(d$dice)

## 

##  F  T 

## 54  6

total <- sum(nrow(d))

One <- table(d$dice)[2]

p <- 1/6

binom.test(x = One, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 6, number of trials = 60, p-value = 0.2233

## alternative hypothesis: true probability of success is not equal to 0.1666667

## 95 percent confidence interval:

##  0.03759127 0.20505774

## sample estimates:

## probability of success 

##                    0.1

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##         T 

## 0.2253457

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##         T 

## 0.1658567

------------------------------------------------------

正確検定ではp = 0.2233である。漸近検定では、連続値補正ありではp = 0.2253、なしではp = 0.1659である。今回は、連続値補正ありのほうが精度が良い。図示も行おう。

------------------------------------------------------

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

min_c <- 3

max_c <- 1

xx <- seq(min(lab) - min_c, max(lab) + max_c, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(One, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図14の描画

dif <- pos[2] - pos[1]

if(One > total*p){

  axx2 <- seq(One, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  s <- abs(One - total * p)

  axx1 <- seq(min(lab) - min_c, total * p - s, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  lwr <- total * p - s + 1/2

  upr <- One + 1/2

} else {

  axx1 <- seq(min(lab) - min_c, One, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  s <- abs(One - total * p)

  axx2 <- seq(total * p + s, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  lwr <- One + 1/2

  upr <- total * p + s + 1/2

}

xx <- seq(min(pos) - min_c*dif, max(pos) + max_c*dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = lwr - min(lab), lty = 2),

         abline(v = upr - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

もう一回試行回数60の時を考える。今度の例は、補正しないほうが近似精度が良くなる例だ。

------------------------------------------------------

xx <- c("F", "F", "F", "F", "F", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "T", "F", "F", "F",

        "T", "F", "F", "F", "T", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "F", "F", "F", "F",

        "F", "T", "F", "F", "F", "F", "F", "F", "F", "F", 

        "F", "F", "F", "F", "F", "F", "F", "F", "F", "F")

yy <- 1:60

d <- data.frame(try = yy, dice = xx)

table(d$dice)

## 

##  F  T 

## 56  4

total <- sum(nrow(d))

One <- table(d$dice)[2]

p <- 1/6

binom.test(x = One, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 4, number of trials = 60, p-value = 0.03664

## alternative hypothesis: true probability of success is not equal to 0.1666667

## 95 percent confidence interval:

##  0.01846178 0.16198676

## sample estimates:

## probability of success 

##             0.06666667

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##          T 

## 0.05674682

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

##          T 

## 0.03766692

threshold <- 0.00001

tmp <- dbinom(0:total, size = total, prob = p)

yy1 <- tmp[tmp > threshold]

lab <- 0:total

lab <- lab[tmp > threshold]

min_c <- 3

max_c <- 1

xx <- seq(min(lab) - min_c, max(lab) + max_c, length = 200)

yy2 <- dnorm(xx, mean = total * p, sd = sqrt(total * p * (1 - p)))

cols <- col_genelator(alpha = 0.5)

fils <- sapply(yy1 <= dbinom(One, size = total, prob = p),

              FUN = function(x){if(x){NA} else {cols[1]}})

pos <- barplotn(yy1, xaxt = "n", mar = c(3.8, 3.8, 1, 1),

                xlab = "成功回数k", ylab = "確率",

                col.fill = fils, space = 0) #図15の描画

dif <- pos[2] - pos[1]

if(One > total*p){

  axx2 <- seq(One, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  s <- abs(One - total * p)

  axx1 <- seq(min(lab) - min_c, total * p - s, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  lwr <- total * p - s + 1/2

  upr <- One + 1/2

} else {

  axx1 <- seq(min(lab) - min_c, One, length = 200)

  ayy1 <- dnorm(axx1, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx1 <- c(axx1, axx1[length(axx1)])

  ayy1 <- c(ayy1, 0)

  s <- abs(One - total * p)

  axx2 <- seq(total * p + s, max(lab) + max_c, length = 200)

  ayy2 <- dnorm(axx2, mean = total * p, sd = sqrt(total * p * (1 - p)))

  axx2 <- c(axx2[1], axx2)

  ayy2 <- c(0, ayy2)

  lwr <- One + 1/2

  upr <- total * p + s + 1/2

}

xx <- seq(min(pos) - min_c*dif, max(pos) + max_c*dif, length = 200)

overdraw(points(xx, yy2, type = "l", col = "red"),

         abline(v = lwr - min(lab), lty = 2),

         abline(v = upr - min(lab), lty = 2),

         polygon(axx1 + dif/2 - min(lab), ayy1, col = cols[4], border = NA),

         polygon(axx2 + dif/2 - min(lab), ayy2, col = cols[4], border = NA),

         axis(1, at = pos, labels = lab, cex.axis = 1.2))

------------------------------------------------------

正確検定はp = 0.03664、連続値補正ありでp = 0.05675、連続値補正なしでp = 0.03767である。このように、成功確率がゆがんでいるときは、連続値補正が不安定であり注意を要する。図12～15を見比べると、成功回数の近似正規分布における位置関係は4パターン存在しており、下記に示したようなことを判定すると、補正したほうが良いか否かを判断できる。

上記の関係を判定する、次のような関数を定義することで、補正をしたほうが良いか否かを簡易的に判定することができる。

------------------------------------------------------

correct_judge <- function(total, success, prob){

  s <- abs(success - total * prob)

  round2 <- function(x){

    fx <- floor(x)

    s <- x - fx

    if(s >= 0.5){

      fx + 1

    } else {

      fx

    }

  }

  if(success > total*prob){

    x <- round2(total * prob - s)

  } else {

    x <- round2(total * prob + s)

  }

  if(x < 0 | x > total){

    dist <- "Check distribution"

    q <- "?"

  } else {

    dist <- "No problem"

    q <- ""

  }

  judge <- dbinom(x, size = total, prob = prob) <= dbinom(success, size = total, prob = prob)

  if(judge){

    paste0(dist, ", Correction better", q)

  } else {

    paste0(dist, ", Non-correction better", q)

  }

}

------------------------------------------------------

この関数は、標本の試行回数total、標本の成功回数success、帰無仮説における成功確率probを引数として与えることで、正規分布による漸近検定は、連続値補正ありのほうが良い（"Correction better"）か、ない方が良い（"Non-correction better"）かを判定してくれる。また、もし、二項分布が0や試行回数n付近に分布しており、近似が怪しそうなら、分布を確認してください（"Check distribution"）を出力する。下記のように使う。

------------------------------------------------------

correct_judge(8, 6, 1/2)

## [1] "No problem, Correction better"

correct_judge(10, 4, 1/2)

## [1] "No problem, Correction better"

correct_judge(40, 16, 1/2)

## [1] "No problem, Correction better"

correct_judge(20, 7, 1/6)

## [1] "No problem, Non-correction better"

correct_judge(40, 10, 1/6)

## [1] "No problem, Correction better"

correct_judge(60, 6, 1/6)

## [1] "No problem, Correction better"

correct_judge(60, 4, 1/6)

## [1] "No problem, Non-correction better"

------------------------------------------------------

出力がこれまでの議論してきたことと一致していることを確認してほしい。例えば、今まで議論していなかった状況、試行回数45、成功回数39、成功確率7/10のときを判定すると、以下のようになる。併せて、本当に判定が正しいかを確認してみる。

------------------------------------------------------

correct_judge(45, 39, 7/10)

## [1] "No problem, Non-correction better"

total <- 45

One <- 39

p <- 7/10

binom.test(x = One, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 39, number of trials = 45, p-value = 0.01393

## alternative hypothesis: true probability of success is not equal to 0.7

## 95 percent confidence interval:

##  0.7320750 0.9494576

## sample estimates:

## probability of success 

##              0.8666667

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

## [1] 0.02278024

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

## [1] 0.01469742

------------------------------------------------------

確かに補正無しのほうが漸近検定の精度が良い。また、試行回数10、成功回数4、成功確率1/6のように二項分布が0付近に分布していそうな状況の出力は以下となる。

------------------------------------------------------

correct_judge(10, 4, 1/6)

## [1] "Check distribution, Correction better?"

total <- 10

One <- 4

p <- 1/6

binom.test(x = One, n = total, p = p)

## 

##  Exact binomial test

## 

## data:  One and total

## number of successes = 4, number of trials = 10, p-value = 0.06973

## alternative hypothesis: true probability of success is not equal to 0.1666667

## 95 percent confidence interval:

##  0.1215523 0.7376219

## sample estimates:

## probability of success 

##                    0.4

if(One > total*p){

  lower <- F

  correct <- -1/2

} else {

  lower <- T

  correct <- 1/2

}

pnorm(One + correct, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

## [1] 0.1197949

pnorm(One, mean = total * p, sd = sqrt(total * p * (1-p)), lower.tail = lower) * 2

## [1] 0.04771488

------------------------------------------------------

"Check distribution"が出力されたので、判定はあてにならないが、無理やり計算してみると、予測は外れているし、また近似精度は全く良くない。こんな時は正確検定を行うべきだ。

おわりに

　今回は二項検定を一例として、正確検定と漸近検定の考え方を紹介した。実際のところ、手計算していな時代とは異なり、計算機の計算速度が高くなっているので、二項検定の計算程度であれば、試行回数が多少大きくても問題とならないことが多い。それゆえ、漸近検定も出番が減りがちであるが、もっと重たい計算（Fisherの正確確率検定など）を行うときにはとても重要である。これを機に理解を深めてもらえたら幸いである。