質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Java

Javaは、1995年にサン・マイクロシステムズが開発したプログラミング言語です。表記法はC言語に似ていますが、既存のプログラミング言語の短所を踏まえていちから設計されており、最初からオブジェクト指向性を備えてデザインされています。セキュリティ面が強力であることや、ネットワーク環境での利用に向いていることが特徴です。Javaで作られたソフトウェアは基本的にいかなるプラットフォームでも作動します。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

プログラミング言語

プログラミング言語はパソコン上で実行することができるソースコードを記述する為に扱う言語の総称です。

Q&A

解決済

2回答

1503閲覧

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

tapico.8190

総合スコア5

Java

Javaは、1995年にサン・マイクロシステムズが開発したプログラミング言語です。表記法はC言語に似ていますが、既存のプログラミング言語の短所を踏まえていちから設計されており、最初からオブジェクト指向性を備えてデザインされています。セキュリティ面が強力であることや、ネットワーク環境での利用に向いていることが特徴です。Javaで作られたソフトウェアは基本的にいかなるプラットフォームでも作動します。

アルゴリズム

アルゴリズムとは、定められた目的を達成するために、プログラムの理論的な動作を定義するものです。

プログラミング言語

プログラミング言語はパソコン上で実行することができるソースコードを記述する為に扱う言語の総称です。

0グッド

0クリップ

投稿2020/08/03 14:49

前提・実現したいこと

( 実現したいこと )

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

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

  • 試行回数nは最大で10,000程度になる
  • 確率pは1/4〜1/65536の値をとる
  • 一度に最大で50パターンの計算をしたいため1パターンあたりの処理速度は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先のサンプルコード通りに書いてみたところです。

java

1public class Pascal { 2 public static void main(String[] args) { 3 // 配列の準備 4 int c[][] = new int [2005][2005];//n=2000段作りきれるように、少し余裕を持って作成します。 5 6 //パスカルの三角形作成 7 int mod = 1000000007;// オーバーフロー対策:問題で指定されます。 8 c[0][0]=1;// 初期値として、最上部の1を与えます。 9 for(int i=0;i<2003;i++) {//2000段作りきるようにループを回します。細かい事を考えないで済むように少し余裕を持たせて回しています。 10 for(int j=0;j<2003;j++) { 11 long tmp = c[i][j]%mod;// オーバーフロー対策あり 12 //long tmp = c[i][j];// オーバーフロー対策なし 13 c[i+1][j]+=tmp; 14 c[i+1][j+1]+=tmp; 15 } 16 } 17 System.out.println(c[5][2]); 18 System.out.println(c[1000][10]); 19 } 20}

処理結果

//オーバーフロー対策ありの結果 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)があればそれを教えて頂きたくお願いします。

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

よろしくお願いします。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

swordone

2020/08/03 15:57

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

2020/08/03 16:38

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

2020/08/05 09:10

swordoneさん ありがとうございます。 競技向けのお話しなのですね。 >実用で出てくるのは余り存じ上げません。 >よろしければ、その調べたサイトを提示していただけないでしょうか? おそらく私が競技プログラミングの記事だときちんと理解できていなかった。 上記コードの参考元などで1000000007で割っているのを見てそれで対応できるものと誤解しておりました。 https://qiita.com/naru7sh/items/7f3f47fbf2e415c0ec86
tapico.8190

2020/08/05 09:15

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

回答2

0

ベストアンサー

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

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

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

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

python

1from decimal import * 2# 素数列挙(エラトステネスの篩) 3def eratosthenes(n): 4 primeNumbers = [] 5 isprime = [True]*(n+1) 6 for i in range(2, n+1): 7 if(isprime[i]): 8 primeNumbers.append(i) 9 temp = i 10 while(temp <= n): 11 isprime[temp] = False 12 temp += i 13 return primeNumbers 14 15# n!を素因数分解(ルジャンドルの定理) 16def primeFactorizeKaijo(n, primeNumbers): 17 assert(n <= max_int) 18 res = {} 19 for i in primeNumbers: 20 if(i > n): 21 break 22 temp = i 23 restemp = 0 24 while(temp <= n): 25 restemp += int(n/temp) 26 temp *= i 27 res[i] = restemp 28 return res 29 30# nを素因数分解する 31def primeFactorize(n, primeNumbers): 32 assert(n <= max_int) 33 res = {} 34 temp = n 35 for i in primeNumbers: 36 if(i > temp): 37 break 38 while(temp % i == 0): 39 if(i in res): 40 res[i] += 1 41 else: 42 res[i] = 1 43 temp /= i 44 return res 45def calculatePowMul(ans, fact, power): 46 for k, v in fact.items(): 47 ans[k] += v*power 48 49# nCk * p^k * q^{n-k} 50n = 10000 51k = 5000 52 53# p = [pの分子, pの分母] 54p = [1, 65535] 55max_int = max([n, k, p[0], p[1]]) 56 57# 素因数分解 58primeNumbers = eratosthenes(max_int) 59pbunbo_fact = primeFactorize(p[1], primeNumbers) 60pbunshi_fact = primeFactorize(p[0], primeNumbers) 61qbunshi_fact = primeFactorize(p[1]-p[0], primeNumbers) 62nkaijo_fact = primeFactorizeKaijo(n, primeNumbers) 63kkaijo_fact = primeFactorizeKaijo(k, primeNumbers) 64nmkkaijo_fact = primeFactorizeKaijo(n-k, primeNumbers) 65 66# 指数同士の足し算をおこない、掛け算する 67ans = [0]*(max_int+1) 68calculatePowMul(ans, pbunbo_fact, -n) 69calculatePowMul(ans, pbunshi_fact, k) 70calculatePowMul(ans, qbunshi_fact, n-k) 71calculatePowMul(ans, nkaijo_fact, 1) 72calculatePowMul(ans, kkaijo_fact, -1) 73calculatePowMul(ans, nmkkaijo_fact, -1) 74 75# 指数の正負で分ける(0は入れない) 76plus = {} 77minus = {} 78for i in primeNumbers: 79 if(ans[i] > 0): 80 plus[i] = ans[i] 81 if(ans[i] < 0): 82 minus[i] = abs(ans[i]) 83 84# 実際に掛け算を行う( 85powplus = 1 86for i in plus: 87 powplus *= pow(i, plus[i]) 88powminus = 1 89for i in minus: 90 powminus *= pow(i, minus[i]) 91 92# 有効桁数を50桁とする 93getcontext().prec = 50 94 95# 実際に割って小数にする 96print(Decimal(powplus)/Decimal(powminus))

text

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

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

投稿2020/08/03 16:03

編集2020/08/03 22:48
Penpen7

総合スコア698

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

tapico.8190

2020/08/05 09:20

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

2020/08/10 04:13

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

0

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

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

投稿2020/08/03 18:43

swordone

総合スコア20649

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

swordone

2020/08/03 18:46

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問