Scala

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

二項分布は試行回数が多くなると階乗周りの計算量が多くなる。一回の事象の発生確率pが十分に小さければ、ポアソン分布で近似ができる。

例えば、5枚の硬貨を投げてすべて表がでるかどうかを64回やったときの分布を二項分布とポアソン分布で比較すると、以下のように高い精度で近似できていることがわかる。

import scala.math.BigInt

def factor(n: BigInt): BigInt = if (n <= 1) 1 else n * factor(n - 1)

def combi(n: BigInt, m: BigInt): BigInt =
  factor(n) / factor(m) / factor(n - m)

def binomialDistribution(n: Int, p: Double): (Int) => Double = {
  import scala.math.pow
  def distribution(x: Int) = combi(n, x).toDouble * pow(p, x) * pow(1 - p, n - x)
  distribution _
}

def poissonDistribution(m: Double): (Int) => Double = {
  import scala.math.pow
  def distribution(x: Int) = pow(m, x).toDouble / pow(2.71828, m).toDouble / factor(x).toDouble
  distribution _
}

val n = 64
val p = 1.0 / 32.0
val binomial = binomialDistribution(n, p)
val poisson = poissonDistribution(n * p)

(0 to 10).foreach {n =>
  println("[n=%d]".format(n))
  println("binomial: %.3f".format(binomial(n)))
  println("poisson: %.3f".format(poisson(n)))
}

[n=0]
binomial: 0.131
poisson: 0.135
[n=1]
binomial: 0.271
poisson: 0.271
[n=2]
binomial: 0.275
poisson: 0.271
[n=3]
binomial: 0.183
poisson: 0.180
[n=4]
binomial: 0.090
poisson: 0.090
[n=5]
binomial: 0.035
poisson: 0.036
[n=6]
binomial: 0.011
poisson: 0.012
[n=7]
binomial: 0.003
poisson: 0.003
[n=8]
binomial: 0.001
poisson: 0.001
[n=9]
binomial: 0.000
poisson: 0.000
[n=10]
binomial: 0.000
poisson: 0.000

3囚人問題は有名だけど、こちらを見かけたのは恥ずかしながら初めて。

「プレイヤーの前に3つのドアがあって、1つのドアの後ろには景品の新車が、2つのドアの後ろにはヤギ(はずれを意味する)がいる。プレイヤーは新車のドアを当てると新車がもらえる。プレイヤーが1つのドアを選択した後、モンティが残りのドアの内ヤギがいるドアを開けてヤギを見せる。

ここでプレイヤーは最初に選んだドアを、残っている開けられていないドアに変更しても良いと言われる。プレイヤーはドアを変更すべきだろうか？」

モンティ・ホール問題

答えは、変更した方が勝率が高くなる。ベイズの定理できちんと考えれば答えは出るのだけど、非常に嘘っぽく感じてしまうのがこの問題の特徴*1。嘘っぽいなら実際やらせて確かめてみればいい。

import scala.util.Random._

def didWinGame(isChange: Boolean): Boolean = {
  val doors: Array[Boolean] = Array(false, false, false)
  doors(nextInt(3)) = true

  val playerChoice = nextInt(3)

  def tryChoosing(): Int = {
    val choice = nextInt(3)
    if (choice != playerChoice && ! doors(choice))
      choice
    else
      tryChoosing()
  }
  val montyChoice = tryChoosing()

  val finalAnswer =
    if (isChange) List(0, 1, 2).filter(
      (x: Int) => x != playerChoice && x != montyChoice
    )(0)
    else playerChoice

  return doors(finalAnswer)
}

def winningRate(doGame: () => Boolean) = {
  val total = 10000
  val winCount = (1 to total).map(_ => doGame()).count((x: Boolean) => x)

  winCount.toDouble / total.toDouble * 100
}

println("not change: %.2f %%".format(winningRate(() => didWinGame(false))))
println("change: %.2f %%".format(winningRate(() => didWinGame(true))))

実行するとほぼ理論値となる。

not change: 33.18 %
change: 66.38 %