5

妥当な精度でベータ分布の逆累積分布関数の計算 (別名、分位数の推定) をサポートする Java ライブラリ/実装を探しています。

もちろん、 apache commons mathも試しましたが、バージョン 3 ではまだ精度に問題があるようです。この質問につながる問題を以下に詳しく説明します。

試行回数が多いベータ分布の信頼区間を計算したいとします。Apacheコモンズ数学では...

final int trials = 161750;
final int successes = 10007;
final double alpha = 0.05d;

// the supplied precision is the default precision according to the source code
BetaDistribution betaDist = new BetaDistribution(successes + 1, trials - successes + 1, 1e-9);

System.out.println("2.5 percentile :" + betaDist.inverseCumulativeProbability(alpha / 2d));
System.out.println("mean: " + betaDist.getNumericalMean());
System.out.println("median: " + betaDist.inverseCumulativeProbability(0.5));
System.out.println("97.5 percentile :" + betaDist.inverseCumulativeProbability(1 - alpha / 2d));

配信する

2.5 percentile :0.062030402074808505
mean: 0.06187249616697166
median: 0.062030258659508855
97.5 percentile :0.06305170793994147

問題は、2.5 パーセンタイルと中央値が同じである一方で、どちらも平均よりも大きいことです。

比較すると、Rパッケージbinomは以下を実現します。

binom.confint(10007+1,161750+2,methods=c("agresti-coull","exact","wilson"))
         method     x      n      mean      lower      upper
1 agresti-coull 10008 161752 0.0618725 0.06070873 0.06305707
2         exact 10008 161752 0.0618725 0.06070317 0.06305756
3        wilson 10008 161752 0.0618725 0.06070877 0.06305703

およびRパッケージ統計

qbeta(c(0.025,0.975),10007+1,161750-10007+1)
[1] 0.06070355 0.06305171

R の結果に次ぐために、Wolfram Alphaが私に言ったことは次のとおりです。

要件に関する最終的な注意事項:

  • これらの計算をたくさん実行する必要があります。したがって、どのソリューションも 1 秒よりも長くかかるべきではありません ((間違っているとはいえ) Apache Commons の数学の 41 ミリ秒と比較すると、これはまだ多くの時間です)。
  • Java内でRを使用できることを認識しています。ここでは詳しく説明しませんが、これは他の方法 (純粋な Java) が失敗した場合の最後のオプションです。

更新 21.08.12

この問題は apache-commons-math の 3.1-SNAPSHOT で修正または少なくとも改善されたようです。上記のユースケースの場合

2.5 percentile :0.06070354581340706
mean: 0.06187249616697166
median: 0.06187069085946604
97.5 percentile :0.06305170793994147

更新 23.02.13

一見すると、この質問とその回答は局所的すぎるかもしれませんが、最初に頭に浮かぶハッカーのアプローチでは、いくつかの数値問題を (効率的に) 解決できないことを非常によく示していると思います。ですので、このまま開いていてほしいです。

4

3 に答える 3

2

この問題は、 apache commons math 3.1.1で修正されました。

上記のテストケースが配信されました

2.5 percentile :0.06070354581334864
mean: 0.06187249616697166
median: 0.06187069085930821
97.5 percentile :0.0630517079399996

これは、r-package stats の結果と一致します。3.1-SNAPSHOT + x バージョンの広範なアプリケーションも問題を引き起こしませんでした。

于 2013-02-23T18:04:34.777 に答える
0

累積分布関数のグラフが非常に平坦な場合 (通常は分布の裾に向かっている場合)、妥当な値に到達するには縦軸で非常に高い精度が必要になるため、この問題を一般的に解決することはおそらく不可能です。横軸は精度。

したがって、累積分布関数から分位数を導出するよりも、分位数を直接計算する関数を使用する方が常に優れています。

もちろん、精度を気にしなければ、方程式 q = F (x) を数値的に解くことができます。F は増加しているので、それは難しくありません。

   double x_u = 0.0;
   double x_l = 0.0;

   // find some interval quantile is in
   if ( F (0.0) > q) {
      while ( F (x_l) > q) {
         x_u = x_l;
         x_l = 2.0 * x_l - 1.0;
      }
   } else {
      while ( F (x_u) < q) {
         x_l = x_u;
         x_u = 2.0 * x_u + 1.0;
      }
   }

   // narrow down interval to necessary precision
   while ( x_u - x_l > precision ) {
      double m = (x_u - x_l) / 2.0;
      if ( F (m) > q ) x_u = m; else x_l = m;
   }     
   // quantile will be within [x_l; x_u]

注意:ベータ分布は区間 [0;1] に存在し、グラフは区間の終わりに向かってかなり急勾配であるため、特にベータ分布で精度が問題になる理由は明確ではありません。

2 番目の注意:上位分位数の計算が間違っています。それは読むべきです

System.out.println( "97.5 percentile :" + betaDist.inverseCumulativeProbability( 1 - alpha / 2d ) );

3 番目の編集:アルゴリズムが修正されました。

于 2012-08-20T09:10:16.587 に答える
0

ライブラリJSci（バージョン1.2 27.07.2010）を見つけて試しました

コードスニペット：

final int trials = 162000;
final int successes = 10000;
final double alpha =0.05d;

BetaDistribution betaDist = new BetaDistribution(successes + 1, trials - successes + 1);
long timeSum = 0;
for(double perc : new double[]{alpha/2,0.5,1-alpha/2}){
    long time = System.currentTimeMillis();
    System.out.println((perc*100) + " percentile :" + betaDist.inverse(perc));
    timeSum += System.currentTimeMillis()-time;
}
System.out.println("Took ~" + timeSum/3 + " per call");

戻った

2.5 percentile :0.060561615036184686
50.0 percentile :0.06172659147924378
97.5 percentile :0.06290542466617127
Took ~2ms per call

内部的には、JohnB によって提案されているように、ルート検索アプローチが使用されます。ProbabilityDistribution#inverseを拡張して、より高い精度を要求できます。残念ながら、大量の反復 (100k) と要求された精度 10^-10 を使用しても、アルゴリズムは依然として戻ります

2.5 percentile :0.06056698485628473
50.0 percentile :0.06173200221779383
97.5 percentile :0.06291087598052053
Took ~564ms per call

今: 誰のコードのほうが間違っていませんか? R または JSci ? ユーザーベースが大きい方がいいと思います...

于 2012-08-20T11:36:00.787 に答える