Calculating binomial distribution

Probability !

計算機を使って確率を考えるのは楽しいことです。特に、実装が考察した通りの結果であったときは嬉しいものです。

しかし、確率を計算するのに予期しない工夫を求められることもあります。例えば二項分布の計算はそれほど簡単ではないことが知られています。今回のエントリでは二択問題を題材にこれを考えてみることにします。


さて、100問の二択問題にまったく適当に答えを記入していった場合、何問正解することができるでしょうか?

この問題はn = 100、p = 1 / 2の二項分布の問題として考えることができます。正解する回数をkとすると、二項分布b(k; n, p)は以下のように分布します。

二項分布b(k; n, p)はある実験をn回独立に行い、ちょうどk回成功しn - k回失敗するのは如何ほどの確率であるかを表しています。

成功する回数を確率変数Xで表すとき、期待値E[X]はnpとして求めることができます。

つまり、100問の二択問題に全く適当に答えを記入していったとしても50問正解することが期待できます。もちろん試行によって結果にばらつきがありますが、この値から遠からぬ結果を期待することができる、ということです。

では100問の二択問題に同じ方法で答えを記入していった場合、60問以上正解する確率は如何ほどでしょうか?

この問いには二項分布b(k; 100, 1 / 2)を[60, 100]の範囲のkで計算した結果を足し合わせることによって求めることができます。

計算すると大体0.028444であることが分かります。およそ2.844%です。つまり、全く適当に答えを記入する方法で60問以上の問題に正解することはあまり起こらないのです。


さて、二項分布b(k; n, p)の式をもう一度眺めてみましょう。

この式は二つの計算から成り立っていることが分かります。二項係数をx、残りをyと置いてみましょう。

xは二項係数を計算します。「n回の試行の中でどのk回を成功とするか」はx通りあることを表しています。

yは確率を計算します。「成功確率がpである実験をn回独立に行い、ちょうどk回成功しn - k回失敗するときの確率」を表しています。

二項分布を計算する上での問題は、x、yのそれぞれが極端な値を取ることです。例えばn = 100のとき、xはk = 50(もしくはx = 51)のときに最大の値となり、その値は30桁の数値となります。

また二項分布b(k; n, p)は確率の値です。確率の公理の一つは確率の値は0以上1以下であることを求めています。つまりxとyの積は必ず[0, 1]の範囲に収まる値となります。

xが0を含まない自然数であることに注意すると、

これより、n = 100のとき、yは小数点の後に少なくとも0が29個並ぶ小さな値であることが分かります。How to calculate binomial probabilities ― The Endeavourは一般的な計算機ではnが170を超えるとxとyはオーバーフローやアンダーフローを起こし、計算不可能になると注意を促しています。


さてこのように極端な大きさの値を伴った計算を行うとき、対数を利用することで実用上十分な近似値を計算できることがあります。

二項分布b(k; n, p)は以下のようなものでした。

右辺を以下のように置いてみましょう。

ただし、x、y、z、u、vは以下のように定義します。

ここで、別の正の数a ≠ 1を用意して、

と置くと、xyzuvは以下のように計算できます。

さらに指数法則により、xyzuvはα、β、γ、δ、εの和を用いて以下のように求めることができます。

さて、

であることと、

に注意すると、α、β、γ、δ、εは以下のように求めることができます。

δとεを計算するのに極端な値は使われないこと、またその結果も極端な値にならないことに注目しましょう。

残るはα、β、γの計算ですが、対数の中の階乗をそのまま計算してオーバーフローを起こしてしまっては元も子もありません。そこで、階乗の対数を直接計算するlgamma関数を使います。lgamma関数はC99とPOSIXで規定された数学ライブラリ関数で、ガンマ関数の絶対値の自然対数を返します。

ガンマ関数は自然数nについて、以下が成立する関数です。

ではガンマ関数を使って、α、β、γを書き換えてみましょう。

aにネイピア数eを選ぶと、α、β、γ、δ、εは以下のように計算できます。

これを用いると、二項分布は以下のように求めることができます。

ただし、


この計算方法を使って、二項分布を計算する関数b(n, p, k)をGnuplotで実装すると以下のようになります(引数の順番に注意してください)。

b(n, p, k) = \
  exp(                                      \
             lgamma(n     + 1)              \
    -        lgamma(n - k + 1)              \
    -        lgamma(    k + 1)              \
    +      k                   * log(    p) \
    + (n - k)                  * log(1 - p) \
  )

まとめ

確率を計算するときはこのエントリで記載したような極端な値を扱わなければならない局面が多々あります。そのようなときに対数の使用を検討してみるのはよいアイデアです。


このエントリをまとめるにあたり、以下のページを大いに参考にしました。氏が記載しているエントリはどれも興味深く、何より面白いので購読してみることをお勧めします。

また、以下の書籍を参考にしました。

確率と計算 ―乱択アルゴリズムと確率的解析―

確率と計算 ―乱択アルゴリズムと確率的解析―

Probability and Computing: Randomized Algorithms and Probabilistic Analysis

Probability and Computing: Randomized Algorithms and Probabilistic Analysis

アルゴリズムイントロダクション 第1巻の付録Cでは数え上げと確率を明瞭、簡潔に解説しています。二項分布の式をb(k; n, p)と記すのにこの書籍を参考にしました。

アルゴリズムイントロダクション 第1巻 数学的基礎とデータ構造

アルゴリズムイントロダクション 第1巻 数学的基礎とデータ構造

  • 作者: T.コルメン,R.リベスト,C.ライザーソン,Thomas H. Cormen,Ronald L. Rivest,Charles E. Leiserso,浅野哲夫,梅尾博司,和田幸一,岩野和生,山下雅史
  • 出版社/メーカー: 近代科学社
  • 発売日: 1995/12
  • メディア: 単行本
  • 購入: 3人 クリック: 144回
  • この商品を含むブログ (34件) を見る

また、このエントリを記載するにあたり、以下のページを参考にしました。

確率に関しては以下のような記事を記載しています。