javaで二項分布を用いた確率計算を高速に行う方法を教えてください

Question

前提・実現したいこと

( 実現したいこと )

javaでandroidアプリの開発をしています。
アプリの中で「コインを5回投げた時に表が2回出る確率」のような「確率pをn回試行した時にk回当たる確率」を求めるプログラムを作りたいと考えています。
下記の前提条件を満たした上で高速に計算できるようにするにはどのようなコードを組めば良いのでしょうか。

( 作成したいアプリでの前提条件 )

• 試行回数nは最大で10,000程度になる
• 確率pは1/4〜1/65536の値をとる
• 一度に最大で50パターンの計算をしたいため１パターンあたりの処理速度は0.1秒以下にしたい

発生している問題・エラーメッセージ

調べたところ二項分布を使えば目的の計算ができることが分かりましたので以下のようなプログラムを組もうと考えいます

• nCrの組み合わせを求める
• 二項分布の計算式 P(X=k)= nCr * (p**k) * (1-p)**(n-k) に当てははめて確率を計算する

組み合わせの計算には、nが大きいことから処理速度の観点で「パスカルの三角形」を利用した計算をしようとしましたが、そこで困っています。

nが大きい場合、いくつかのサイトでオーバーフロー対策として計算結果を素数1000000007で割った余りを取ることを推奨していましたが、1000000007で割った余りでは上記二項分布の計算式に当てはめられず上手く計算できません。
かと言って、オーバーフロー対策をしないと正しい結果が得られない状態で困っています。

該当のソースコード

こちらの記事を参考に以下のようなコードを書きました。
https://qiita.com/naru7sh/items/7f3f47fbf2e415c0ec86

米nは最大で10,000となりますが、まずは試しにURL先のサンプルコード通りに書いてみたところです。

public class Pascal {
    public static void main(String[] args) {
        // 配列の準備
        int c[][] = new int [2005][2005];//n=2000段作りきれるように、少し余裕を持って作成します。

        //パスカルの三角形作成
        int mod = 1000000007;// オーバーフロー対策：問題で指定されます。
        c[0][0]=1;// 初期値として、最上部の1を与えます。
        for(int i=0;i<2003;i++) {//2000段作りきるようにループを回します。細かい事を考えないで済むように少し余裕を持たせて回しています。
            for(int j=0;j<2003;j++) {
                long tmp =  c[i][j]%mod;// オーバーフロー対策あり
                //long tmp =  c[i][j];// オーバーフロー対策なし
                c[i+1][j]+=tmp;
                c[i+1][j+1]+=tmp;
            }
        }
        System.out.println(c[5][2]);
        System.out.println(c[1000][10]);
    }
}

処理結果

//オーバーフロー対策ありの結果
10 　　　　　← 正しい5C2
302505679 　← 1000000007で割った余りらしいがこれをどう使えばいいのかが分からない

//オーバーフロー対策なしの結果
10 　　　　　← 正しい5C2 (nが小さいと正しい答えになる)
1774896272 ← 正しくない1000C10 (nが大きいと値がおかしい。本来は263,409,560,461,970,212,832,400)

試したこと

他の方法がないか以下のことを試しましたがともに上手く行きませんでした。

• nの値が大きいために単純なforループでの計算速度が遅く使えませんでした。
• Apacheのmath3ライブラリの CombinatoricsUtils.binomialCoefficient を試しましたが、こちらも計算速度が遅い上にnが大きすぎるために上手く計算できませんでした。

環境

• java
• android8以上

質問

• 上記パスカルの三角形を利用したコードで1000000007で割った余りを利用して、二項分布を使い「確率pをn回試行した時にk回当たる確率」を求めるにはどうすれば良いでしょうか
• そもそものやり方がよくない場合、別の手法(計算方法、コードの書き方、別のライブラリetc)があればそれを教えて頂きたくお願いします。

※正しく、かつ、高速であれば今試しているパスカルの三角形を利用したやり方にはこだわりません。

よろしくお願いします。

Comments

[swordone]

そもそも、1000000007で割った余りを取った時点で、確率の計算としては成り立たなくなります。
通常、1000000007で割った余りを取るのは、競技プログラミングで途中で巨大な数が出てくる場合に行う処置であり、実用で出てくるのは余り存じ上げません。
よろしければ、その調べたサイトを提示していただけないでしょうか？

[tiitoi]

そもそも、アプリの仕様として、そのぐらい n が大きい場合において、正確な値を計算する必要はあるのでしょうか？

[tapico.8190]

swordoneさん

ありがとうございます。
競技向けのお話しなのですね。

>実用で出てくるのは余り存じ上げません。
>よろしければ、その調べたサイトを提示していただけないでしょうか？

おそらく私が競技プログラミングの記事だときちんと理解できていなかった。
上記コードの参考元などで1000000007で割っているのを見てそれで対応できるものと誤解しておりました。
https://qiita.com/naru7sh/items/7f3f47fbf2e415c0ec86

[tapico.8190]

tiitoiさん

多少の誤差(数%)は問題ないと考えていますが、なるべく正確な値を計算したいと考えています。

Answer 1

上記パスカルの三角形を利用したコードで1000000007で割った余りを利用して、二項分布を使い「確率pをn回試行した時にk回当たる確率」を求めるにはどうすれば良いでしょうか

競技プログラミングで組み合わせを求めるときはオーバーフローを防ぐため余りを使うことはやりますが、確率で余りをつかうことはできないと思います。

二項分布にはほかの分布で近似されることが知られています。正規分布(np>5 np(1-p)>5のとき) あるいはポアソン分布(nが大きくpが十分に小さい場合)で近似できるため、こころへんの性質をつかうといいのではないでしょうか。

追記) 愚直に計算を行うのではなく、前処理で素因数分解を行うことで高速(マイクロ秒くらい)に解けることを確認しました。
まずは、エラトステネスの篩で素数列挙を行い, 通常の整数では素数で順番に割り算してどこまで割れるか試していき、階乗ではルジャンドルの定理を使って、素因数分解を行います。
指数の足し算を行ったあと、実際に掛け算を行って、有理数(powplus/powminus)を求めてやります。
(なお、掛け算時に累乗を計算する必要がありますが、これは繰り返し二乗法を使って高速化します。)
ここまでは整数同士の掛け算割り算なので誤差はありませんが、最後に割り算を行い小数に直します。
(試作のためpythonですが, 流れさえわかれば, javaでも組めると思います。javaでも書こうと思いましたが、pythonで力尽きました)

from decimal import *
# 素数列挙(エラトステネスの篩)
def eratosthenes(n):
    primeNumbers = []
    isprime = [True]*(n+1)
    for i in range(2, n+1):
        if(isprime[i]):
            primeNumbers.append(i)
            temp = i
            while(temp <= n):
                isprime[temp] = False
                temp += i
    return primeNumbers

# n!を素因数分解(ルジャンドルの定理)
def primeFactorizeKaijo(n, primeNumbers):
    assert(n <= max_int)
    res = {}
    for i in primeNumbers:
        if(i > n):
            break
        temp = i
        restemp = 0
        while(temp <= n):
            restemp += int(n/temp)
            temp *= i
        res[i] = restemp
    return res

# nを素因数分解する
def primeFactorize(n, primeNumbers):
    assert(n <= max_int)
    res = {}
    temp = n
    for i in primeNumbers:
        if(i > temp):
            break
        while(temp % i == 0):
            if(i in res):
                res[i] += 1
            else:
                res[i] = 1
            temp /= i
    return res
def calculatePowMul(ans, fact, power):
    for k, v in fact.items():
        ans[k] += v*power

# nCk * p^k * q^{n-k}
n = 10000
k = 5000

# p = [pの分子, pの分母]
p = [1, 65535]
max_int = max([n, k, p[0], p[1]])

# 素因数分解
primeNumbers  = eratosthenes(max_int)
pbunbo_fact   = primeFactorize(p[1],      primeNumbers)
pbunshi_fact  = primeFactorize(p[0],      primeNumbers)
qbunshi_fact  = primeFactorize(p[1]-p[0], primeNumbers)
nkaijo_fact   = primeFactorizeKaijo(n,    primeNumbers)
kkaijo_fact   = primeFactorizeKaijo(k,    primeNumbers)
nmkkaijo_fact = primeFactorizeKaijo(n-k,  primeNumbers)

# 指数同士の足し算をおこない、掛け算する
ans = [0]*(max_int+1)
calculatePowMul(ans, pbunbo_fact,   -n)
calculatePowMul(ans, pbunshi_fact,   k)
calculatePowMul(ans, qbunshi_fact, n-k)
calculatePowMul(ans, nkaijo_fact,    1)
calculatePowMul(ans, kkaijo_fact,   -1)
calculatePowMul(ans, nmkkaijo_fact, -1)

# 指数の正負で分ける(0は入れない)
plus  = {}
minus = {}
for i in primeNumbers:
    if(ans[i] > 0):
        plus[i] = ans[i]
    if(ans[i] < 0):
        minus[i] = abs(ans[i])

# 実際に掛け算を行う(
powplus = 1
for i in plus:
    powplus *= pow(i, plus[i])
powminus = 1
for i in minus:
    powminus *= pow(i, minus[i])

# 有効桁数を50桁とする
getcontext().prec = 50

# 実際に割って小数にする
print(Decimal(powplus)/Decimal(powminus))

CPU times: user 3 µs, sys: 0 ns, total: 3 µs
Wall time: 5.25 µs
6.3420873578683450119531499686198773734249989811395E-21075

ほぼ無視できるような確率ですね。

Comments

[tapico.8190]

サンプルプログラムまで記載して頂きありがとうございます。
とても参考になります。
ご説明頂いた内容とコード内容を理解してjavaで実装できるか工夫してみたいと思います。

[tapico.8190]

内容をきちんと理解するのに少し時間がかかりました。
計算方法を工夫すればここまで早く計算できるんですね。
大変参考になりました。ありがとうございます。

Answer 2

この計算をする、このようなサイトがあります。
二項分布 - 高精度計算サイト

このサイトにおいても、n=1000程度の計算で4秒前後かかります。n=10000でもやってみましたが、結果が返ってきませんでした。
nが小さいことが保証されているならともかく、n=10000まであり得る状況で「1つの処理を0.1秒以内」は到底不可能です。
高速に出したいというなら、nCrの全パターンのリストをデータベースのようなものに予め作っておく必要があるでしょう。
n=10000まで全部作っても5000パターン程度ですし、nCr=nC(n-r)を利用すればさらに半分程度になります。
ただ、通常の数値型では無理なので、文字列で保存しておいて、計算する際にはBigDecimalを使うなどする必要があるでしょう。

Comments

[swordone]

なお、最も大きくなる10000C5000は、値そのものではなく桁数を調べてみたところ、3009桁あるようです。