多項分布

1 = 黒、2 = 赤、... としてそれぞれの生起確率 \(p_k\) を調整できる多項分布シミュレーション。+1 ボタンを押すと試行を 1 回行い観測値 \(\hat{\vector{x}}\) が得られる。\(K\) の調整は重みを 0 にすることで行う。

定義と性質

独立した \(K\) 個の事象 \(k \in \{1, 2, \ldots, K\}\) のいずれか一つが観測される確率をそれぞれ \(\vector{p} = (p_1, p_2, \ldots, p_K)\) とする。この試行を \(n\) 回繰り返したときに各事象の観測された回数を確率変数 \(x_k\) とする。

このとき、事象ごとの観測回数が \(\vector{x} = (x_1, x_2, \ldots, x_K)\) となる確率は \(\vector{p}, n\) をパラメータとする多項分布 (multinomial distribution) に従う。\[ {\rm Multi}(\vector{x} ; \vector{p}, n) = \frac{n!}{\prod_{k=1}^K x_k!} \prod_{k=1}^K p_k^{x_k} \] ここで \(\sum_{k=1}^K p_k = 1\), \(\sum_{k=1}^K x_k = n\) である。\(x_k\) が離散値でなく実数の場合はガンマ関数を用いて表すことができる。

以下は確率変数 \(\vector{x}\) と生起確率 \(\vector{p}\) を入力するとその事象が起きうる確率を算出する。

多項分布の最尤推定、MAP 推定、ベイズ推定は多項分布の推定を参照。

確率変数 \(\vector{x}\) の組み合わせ

確率変数 \(\vector{x}\) は「合計で \(n\) となる \(K\) 個の非負整数の数列」、つまり \(\sum_{i=1}^K x_i = n\) を満たす数列と考える事ができる。これは重複組み合わせとみなすことができるため、\(\vector{x}\) が取り得る組み合わせ数は式 (\(\ref{combinations}\)) で求めることができる。\[ \begin{equation} {}_K H _n = {}_{K+n-1} C_n = \frac{(K+n-1)!}{n! (K-1)!} \label{combinations} \end{equation} \] 例えば \(K=5\), \(n=6\) の場合 \({}_5H_6=\frac{(5+6-1)!}{6! (5-1)!}=210\) の組み合わせが存在する (Wolfram|Alpha による結果)。通常の組み合わせと同様に \(K\) と \(n\) の増加によって容易に計算が非現実的となることに注意。例えば \(K=1000\), \(n=25\) では \(8.7 \times 10^{49}\) の組み合わせが存在し、一つの確率計算が 1nsec で終わるとしても全てに \(2.8 \times 10^{33}\) 年を必要とする。

実装例

以下は \({}_K H _n\) 個のすべての確率変数の組み合わせとその生起確率を出力するプログラムの実装例である。

// 階乗
def factrial(n:Int):Double = if(n <= 1) 1.0 else factrial(n - 1) * n

// 総乗
def product(min:Int, max:Int)(generator:Int => Double):Double = (min to max).map(generator).product

// 多項分布に従う確率を求める
def multi(x:Array[Int], p:Array[Double]):Double = {
  assert(p.length == p.length && math.abs(p.sum - 1.0) <= 1e-15)
  val k = x.length
  val n = x.sum
  factrial(n) / product(0, k - 1) { i => factrial(x(i)) } * product(0, k - 1) { i => math.pow(p(i), x(i)) }
}

// 確率変数の組み合わせを生成する (総数が膨大になる可能性があるためコールバックする)
def generateRandomVariables(k:Int, n:Int)(callback:Array[Int] => Unit):Unit = {
  def _generate(x:Array[Int], n:Int, i:Int, callback:Array[Int] => Unit):Unit = {
    if(i == x.length) {
      if(n == 0) {
        callback(x)
      }
    } else for(value <- 0 to n) {
      x(i) = value
      _generate(x, n - value, i + 1, callback)
    }
  }

  _generate(new Array[Int](k), n, 0, callback)
}

val p = Array(0.1, 0.2, 0.4, 0.2, 0.1)
var combinations = 0
var totalProbability = 0.0
generateRandomVariables(k = 5, n = 6) { x =>
  val probability = multi(x, p)
  println(f"${x.mkString("[", ", ", "]")}%s: $probability%.6f")
  combinations += 1
  totalProbability += probability
}
println(f"$combinations%,d combinations, $totalProbability%f probability total")
// [0, 0, 0, 0, 6]: 0.000001
// [0, 0, 0, 1, 5]: 0.000012
// [0, 0, 0, 2, 4]: 0.000060
// ...
// [5, 1, 0, 0, 0]: 0.000012
// [6, 0, 0, 0, 0]: 0.000001
// 210 combinations, 1.000000 probability total 

共通性と近似

カテゴリカル分布との共通性

多項分布はカテゴリカル分布に従う試行を \(n\) 回繰り返し、各事象の観測された回数を \(x_k\) としている。言い換えれば \(n=1\) とした多項分布はカテゴリカル分布と等価である。カテゴリカル分布が「どの事象が観測されるか」を論点としていたのに対して、多項分布は「どの事象が何回観測されるか」までを論点としていることから、考え方としてカテゴリカル分布を多試行の方向に一般化している。

二項分布との共通性

\(K=2\) とした場合、多項分布の事象 \(k \in \{1,2\}\) はベルヌーイ試行における \(b \in \{0, 1\}\) と見なすことができる。ベルヌーイ試行を \(n\) 回行うとは二項分布と等価である。二項分布が「1が何回観測されるか (=0が何回観測されるか)」を論点としていたのに対して、多項分布は「どの事象が何回観測されるか」までを論点としていることから、考え方として二項分布を多事象の方向に一般化している。

適用例

「サンプルとして抽出するものが有限個のどれかに属している」というモデルに適応できる。よく例示されるのはサイコロで、\(K=6\) のサイコロを \(n\) 回振ったときにそれぞれの目の出た回数が \(\vector{x}\) となる確率分布、といった説明がなされる。サイコロ以外では宝くじの抽選で回転する的に矢を当て \(K=\{0, 1, \ldots, 9\}\) の番号を選択するのも同じモデルだろう。

多項分布で注目するのは確率変数の出現回数のみで、出現する順序は無視されることに注意。

Example 1

一様に落ちてくる \(n\) 個のボールが \(K\) 個に区切られた区画のいずれかに確率 \((p_1, p_2, \ldots, p_K)\) で入る。

上記の例でボールがそれぞれの区画に入る確率を \(\vector{p} = (0.10, 0.38, 0.14, 0.30, 0.08)\)、落とすボールの数 \(n=10\) とした場合に \(\vector{x} = (1, 4, 2, 2, 1)\) となる確率 (つまり区画1に1個、区画2に4個、区画3に2個... となる確率) は: \[ \frac{10!}{1! \times 4! \times 2! \times 2! \times 1!} (0.10^1 \times 0.38^4 \times 0.14^2 \times 0.30^2 \times 0.08^1) = 0.0111 \] となる。

Example 2

サイコロは \(K=6\) 全ての目の出現確率が等しく \(\vector{p} = (1/6, 1/6, ..., 1/6)\) である。これを \(n=3\) 回振って (あるいは \(3\) 個同時に振って) 2, 3, 6 が出る \(\vector{x} = (0, 1, 1, 0, 0, 1)\) に対する確率は: \[ \frac{3!}{0! \times 1! \times 1! \times 0! \times 0! \times 1!} \left(\frac{1}{6^0} \times\frac{1}{6^1} \times \frac{1}{6^1} \times \frac{1}{6^0} \times \frac{1}{6^0} \times\frac{1}{6^1}\right) = 0.0278 \] となる。

Example 3

ある駅で列車に乗る人数の偏りを考える。\(K=4\) 両編成の車両それぞれに乗車する確率を \(\vector{p} = (0.1, 0.2, 0.4, 0.3)\) とする。この駅で \(n=20\) 人が乗車するとき \(\vector{x} = (9, 6, 4, 1)\) の配分で人が乗る確率は: \[ \frac{20!}{9! \times 6! \times 4! \times 1!} (0.1^9 \times 0.2^6 \times 0.4^4 \times 0.3^1) = 1.907 \times 10^{-7} \] となる。確率 \(\vector{p}\) が 3 両目に偏っているのに対して (たぶん 3 両目付近に改札があるのだろう)、この確率変数 \(\vector{x}\) は 1 両目に偏っているため、偶然起きうる可能性はほとんどないと言える。

Example 4

100個に1個の不良品が混じる生産レーンを考える。つまり \(K\) は正常か不良品かの \(2\)、確率は \(\vector{p} = (99/100, 1/100)\) となる。このレーンから任意の \(n=3\) 個を取り出したときに不良品が混じっている \(\vector{x} = (2, 1) \mbox{ or } (1, 2) \mbox{ or } (0, 3)\) から確率は: \[ \frac{3!}{2! \times 1!}\left\{(\frac{99}{100})^2 \times (\frac{1}{100})^1\right\} + \frac{3!}{1! \times 2!}\left\{(\frac{99}{100})^1 \times (\frac{1}{100})^2\right\} + \frac{3!}{0! \times 3!}\left\{(\frac{99}{100})^0 \times (\frac{1}{100})^3\right\} = 0.03 \] となる。これは「任意の 1 個取り出して不良品を引くか」を \(3\) 回試行するのと等価であることが分かる。

Example 5

ビタミンC、カルシウム、鉄の 3 種類の錠剤を 3:2:2 の比率で十分な量を混ぜ合わせた。この中から任意の 4 錠を取り出したときにビタミンCが1錠、カルシウムが2錠、鉄が1錠となる確率を求める。

この場合、確率変数が取り得る値は VC, Ca, Fe であるから \(K=3\) である。それぞれの錠剤の出現確率は配合の割合 (ディリクレ分布の最頻度) から \(\vector{p}=(3/7, 2/7, 2/7)\) となる (錠剤を取り出すことによる若干の配合比率の変化は無視する)。試行回数は \(n=4\)、\(\vector{x}=(1, 2, 1)\) であることから: \[ \frac{4!}{1! \times 2! \times 1!}\left\{ \left(\frac{3}{7}\right)^1 \times \left(\frac{2}{7}\right)^2 \times \left(\frac{2}{7}\right)^1 \right\} = 0.1200 \] となる。

参照

1. 多項分布の推定 (MOXBOX)