95%の信頼区間とは

10000本のクジにx本の当たりが入っているとする。ここから100本のクジを引いた結果を元に、当たりクジが何本入っているか区間推定したい。

非復元抽出なので厳密にいえば100本中の当たりクジの本数yは超幾何分布に従うが、抽出する数と比較して全体数が大きいので二項分布で近似して良いだろう。さらに引いたクジの本数も十分に大きいので、正規分布として近似できる。

当たりくじの割合がp = x / 10000でそこから100本クジを引くのだから、確率変数yが二項分布だと考えれば平均値は100p、分散は100p(1-p)となる。これを元に正規化した、(y - 100p) / 10√(p(1-p))は標準正規分布に従うと言える。標準正規分布であれば-1.96から1.96の間の値をとる確率が95%であるから、-1.96 < (y - 100p) / 10√(p(1-p)) < 1.96 をpについて解けば*1、pがその範囲内に入る確率は95%であるということになる。これが95%の信頼区間。

見方を変えれば、この方法で推定すると、5%の確率で外れるってことになる。本当かどうか、シミュレーションしてみる。100個を無作為抽出してそこから95%の信頼区間を求めるという試行を1000回やって、実際の確率pが信頼区間に入った回数、入らなかった回数をカウントしてみる。

import scala.math._
import scala.util.Random._

val n: Int = 10000
val x: Int = 4321
val p: Double = x.toDouble / n.toDouble
val sampleN: Int = 100
val count: Int = 1000

val population: Array[Boolean] =
  (1 to n).map(n => if (n <= x) true else false).toArray

def square(x: Double) = x * x

def randomSampling[T](x: Seq[T], n: Int) = {
  val shuffled: Seq[T] = shuffle(x)
  shuffled.take(n)
}

def solveQuadratic(
  a: Double = 1, b: Double = 0, c: Double = 0
): (Double, Double) = {
  val hanbetu = sqrt(square(b) - 4 * a * c)
  ((- b - hanbetu) / 2.0 / a, (- b + hanbetu) / 2.0 / a)
}

def estimate (samples: Seq[Boolean]): (Double, Double) = {
  val n = samples.length
  val x = samples count {x => x}

  solveQuadratic(n + square(1.96), - 2 * x - square(1.96), square(x) / n)
}

val results = (1 to count).map { _ =>
  val sample = randomSampling(population, sampleN)
  val (min, max) = estimate(sample)
  if (min <= p && p <= max) true else false
}

println(
  "hit: %d, miss: %d" format (
    results count {x => x}, results count {x => ! x}
  )
)

結果は以下の通り、やはり5%は推定に失敗するという結果になる。

hit: 952, miss: 48

*1:二乗して二次不等式にして解く