mathjsで桁あふれする規模の二項分布の計算をしてみた

やってみたら、小規模なら簡単でも大規模だと工夫が必要であった。
別に言語にはこだわりは無いのだけど、今回は JavaScript を使用した。

実装的な結論から言うと BigNumber を使うだけの話です。
数式からプログラムを実装することに慣れていない人は参考になるかもしれない。

環境

Ubuntu 18.04
node: v10.16.2

Step 1. 数式のおさらい

$n$枚のコインを投げて、$x$枚が表になるような確率は、次のように表される。
Xは確率変数である。具体的には X={0枚のコインが表(すべて裏)の事象, 1枚のコインが表になる事象, ...} である。

P_{X}(x) = {}_{n}C_{k} p^{x}(1-p)^{n-x}

答え合わせしやすいように、次のサイトと同じパラメーターで数式、およびプログラムの実装例を示す。

13-1. 二項分布 | 統計学の時間 | 統計WEB

具体的に、コインを10回投げて、表が3回となる確率を考える。
根元事象は 3枚のコインが表になる事象 で、確率変数$x$の値は 3 である。
なお、理想的なコインを考えて、$p=0.5$とする。

\begin{aligned}
P_{X}(3) &= {}_{10}C_{3} 0.5^{3}(1-0.5)^{10-3} \\
&= \dfrac{10 \cdot 9 \cdot 8}{3 \cdot 2 \cdot 1} \cdot 0.125 \cdot 0.0078125 \\
&= 0.1171875
\end{aligned}

mathjs (Node.js)を用いた実装例
数式そのままです。特に解説はいりませんね。

const math = require('mathjs')

const n = 10
const k = 3
const p = 0.5

p_n_k = math.combinations(n, k) * math.pow(p, k) * math.pow(1 - p, n - k)
console.log('P(n=%d,x=%d)=%f', n, k, p_n_k)

実行例
結果も手計算で求めたものに一致します。

$ node throw_coin_10_times_3_head.js
P(n=10,x=3)=0.1171875
$

Step 2. コインを10回投げて、表が3回以上となる確率

コインを10回投げて、表が3回以上となる確率を考える。
根元事象の集合は 3枚のコインが表になる事象 ∪ 4枚のコインが表になる事象 ∪…∪ 10枚のコインが表になる事象 で、確率変数$x$={3, 4, ..., 10} である。離散型の確率分布であることに注意(3.5とかは存在していない)。この確率、どのくらいになると思いますか？

数式はこんな感じです。

\begin{aligned}
P_{X}(3 \cup 4 \cup \cdots \cup 10)
&= P_{X}(3) +  P_{X}(4)  +  \cdots + P_{X}(10) \\
&= \sum_{x=3}^{10} {}_{10}C_{x} 0.5^{x}(1-0.5)^{10-x} \\
&= {}_{10}C_{3} 0.5^{3}(1-0.5)^{10-3} + \cdots + {}_{10}C_{10} 0.5^{10}(1-0.5)^{10-10} \\
&\simeq 0.945
\end{aligned}

mathjs (Node.js)を用いた実装例
for文でぐるぐる回しているだけです。まだ、簡単なはず。特に解説はいりませんね。

const math = require('mathjs')

const k = 3
const p = 0.5

for (let n = 10; n <= 10; n++) {
    let sum = 0.0;
    for (let i = k; i <= n; i++) {
        p_n_k = math.combinations(n, i) * math.pow(p, i) * math.pow(1 - p, n - i)
        console.log('P(n=%d,x=%d)=%f', n, i, p_n_k)
        sum += p_n_k
    }
    console.log('P(n=%d,3<=x)=%f', n, sum)
}

実行例
結果も手計算で合計したものに大体一致します。ちなみに、約95%ですね。

$ node throw_coin_10_times_3to10_head.js
P(n=10,x=3)=0.1171875
P(n=10,x=4)=0.205078125
P(n=10,x=5)=0.24609375
P(n=10,x=6)=0.205078125
P(n=10,x=7)=0.1171875
P(n=10,x=8)=0.0439453125
P(n=10,x=9)=0.009765625
P(n=10,x=10)=0.0009765625
P(n=10,3<=x)=0.9453125
$

Step 3. : 5%のあたりくじを400本引いて、あたりの本数が20回以上となる確率

Step 2.の数値を置き換えができれば、数式は簡単です。
$p= 0.05, n=400, k=\{20, 21, \cdots, 400\}$

数式はこんな感じです。

\begin{aligned}
P_{X}(20 \cup 21 \cup \cdots \cup 400) &= P_{X}(20) +  P_{X}(21)  +  \cdots + P_{X}(21) \\
&= \sum_{x=20}^{400} {}_{400}C_{x} \cdot 0.05^{x}\cdot(1-0.05)^{400-x} \\
&= \text{???}
\end{aligned}

ただ、 ₄₀₀C_x をこれまでのやり方で数値計算してみようとすると、桁あふれします。

mathjs (Node.js)を用いた実装例(順列の計算だけ)
順列の計算だけを抜き出したコードです。

const { combinations } = require('mathjs')
let n = 400;
for (let k = 20; k <= 400; k++) {
    console.log('(n=%d, k=%d) -> %d', n, k, combinations(n, k))
}

実行結果
Number型で桁あふれしてしまっています。Σの合計どころではありません。

$ node mathjs_node_combinations_sample.js
(n=400, k=20) -> 2.788360983670897e+33
(n=400, k=21) -> 5.045605589499718e+34
(n=400, k=22) -> 8.692202356456333e+35
(n=400, k=23) -> 1.4285445611915188e+37
...snip
(n=400, k=120) -> 5.7072682234758346e+104
(n=400, k=121) -> 1.3206901674158961e+105
(n=400, k=122) -> Infinity
(n=400, k=123) -> Infinity
(n=400, k=124) -> Infinity
(n=400, k=125) -> Infinity
(n=400, k=126) -> Infinity
...snip

Final Step. :プログラムの実装例

実は、桁あふれ以外にも、加算時の情報落ちなど、もとのコードは色々問題があります。

mathjs は BigNumber を使えるので、全面的にこれを使います。これで、BigNumberで扱える範囲の問題なら数値計算ができます。

mathjs (Node.js)を用いた実装例(BigNumer使用)

const math = require('mathjs')

const n = 400
const k = 20
const p = 0.05

let sum = math.bignumber(0.0);
for (let i = k; i <= n; i++) {
    let c_n_k = math.bignumber(math.combinations(math.bignumber(n), math.bignumber(i)))
    let pow_p = math.pow(math.bignumber(p), math.bignumber(i))
    let pow_np = math.pow(math.bignumber(1 - p), math.bignumber(n - i))
    p_n_k = math.bignumber(c_n_k * pow_p * pow_np)
    console.log('P(n=%d,k=%d)=%s * %s * %s -> %s', n, i, c_n_k.toPrecision(8), pow_p.toPrecision(4), pow_np.toPrecision(4), p_n_k.toPrecision(4))
    sum = math.add(sum, p_n_k)
}

console.log('P(n=%d, %d<=x) = %s', n, k, sum.toFixed(4))

実行結果
表示桁が長い部分(4桁以上)は、削っています。Infinity扱いとなっていた計算に成功しています。

$ node lottery_400_hit_over_20.js
P(n=400,k=20)=2.7883610e+33 * 9.537e-27 * 3.427e-9 -> 0.09114
P(n=400,k=21)=5.0456056e+34 * 4.768e-28 * 3.608e-9 -> 0.08680
P(n=400,k=22)=8.6922024e+35 * 2.384e-29 * 3.798e-9 -> 0.07870
P(n=400,k=23)=1.4285446e+37 * 1.192e-30 * 3.998e-9 -> 0.06808
...snip
P(n=400,k=120)=5.7072682e+104 * 7.523e-157 * 5.789e-7 -> 2.486e-58
P(n=400,k=121)=1.3206902e+105 * 3.762e-158 * 6.094e-7 -> 3.027e-59
P(n=400,k=122)=3.0202669e+105 * 1.881e-159 * 6.414e-7 -> 3.644e-60
P(n=400,k=123)=6.8262942e+105 * 9.404e-161 * 6.752e-7 -> 4.334e-61
P(n=400,k=124)=1.5249060e+106 * 4.702e-162 * 7.107e-7 -> 5.096e-62
P(n=400,k=125)=3.3669925e+106 * 2.351e-163 * 7.482e-7 -> 5.922e-63
P(n=400,k=126)=7.3485948e+106 * 1.175e-164 * 7.875e-7 -> 6.803e-64
...snip
P(n=400,k=247)=1.5240999e+114 * 4.422e-322 * 0.0003906 -> 2.618e-211
P(n=400,k=248)=9.4027133e+113 * 2.211e-323 * 0.0004111 -> 7.640e-213
P(n=400,k=249)=5.7398089e+113 * 1.105e-324 * 0.0004328 -> 0.000
P(n=400,k=250)=3.4668446e+113 * 5.527e-326 * 0.0004556 -> 0.000
...snip
P(n=400,k=397)=10586800 * 3.098e-517 * 0.8574 -> 0.000
P(n=400,k=398)=79800.000 * 1.549e-518 * 0.9025 -> 0.000
P(n=400,k=399)=400.00000 * 7.745e-520 * 0.9500 -> 0.000
P(n=400,k=400)=1.0000000 * 3.873e-521 * 1.000 -> 0.000
P(n=400, 20<=x) = 0.5320
$

所感

数値計算は奥深いですね。
実際、BigNumer の範囲では計算はできるのですが、逆に言えばできることはその範囲までです。
BigNumberを使用すると整数型や浮動小数点の計算に比べてパフォーマンス上の問題もあって、ループ回数を増やしたりすると時間がかかり、時間的にできることも縮小します。
(最近流行のDeep Learningが流行っているのは、数値計算ライブラリのアルゴリズム、計算時間、ハードウェアなどを活用して、計算が可能になったことが一つのピースとなったと言えます。)

こういったAPIなんて使用せず、自前で実装するなら、計算順序などを工夫して情報落ちなど避けることはできます。
今回で言うと、 $0.05^{122}$ とかですね。 1.880e-159 とかほぼ0です。今回は　(デカイ数) × (小さい数) × (小さい数) という構造をしていて、₄₀₀C_x が math.combinations(n,k) と math.pow(p, k) などと、数式が単位にAPIがあるのは、とても分かりやすい。ただ、分かりやすくAPIが意味的(数式的に)に分離されていることとは、必ずしも数値計算に有利だとは限らない。

こういうのは、素朴に数式をハードウェアにあわせて最適化するアプローチもあるが一方で、シミュレーテッド・アニーリングのアルゴリズムのように意味的な関連付けは困難だが、入力とそれに対する解の性質は分かっていて使い物になる……そういうアプローチもある。今回のは、実装は綺麗に見えて、計算機科学的にダサいそういう感じのものであろう。

参考

数学・計算機科学関連

13-1. 二項分布 | 統計学の時間 | 統計WEB
順列・組合せ - 高精度計算サイト
異なる n個のものから r個を選ぶ組み合わせの総数 nCr を求めます。
情報落ち、桁落ち、丸め誤差、打切り誤差の違い

プログラミング関連

その他

「世界最速・最大規模」──東芝、量子コンピュータより高速に組み合わせ最適化問題を計算するアルゴリズムを開発 - ITmedia NEWS