如何计算在Java中的逆累积beta分布函数(How to calculate the invers

我要寻找支持的Beta分布（分位数的估计又名） 以合理的精度逆累积分布函数的计算一个java库/执行。

当然，我已经试过阿帕奇百科全书数学，但在第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 -package 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 -package 统计

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

第二来自R的结果，这里是Wolfram Alpha的告诉我

InverseBetaRegularized [0025.10007 + 1.161750至10007 + 1] => 0.06070354631 ...
InverseBetaRegularized [0975.10007 + 1.161750至10007 + 1] => 0.06305170794 ...

在要求最后要注意的：

我需要运行很多这些计算的。因此，任何解决方案不应该需要更长的时间超过1秒（这仍然是一个很大相比的41ms（尽管是错误的）阿帕奇百科全书数学）。
我知道一个可以的Java中使用R上。至于原因，我在这里就不细说了，这是最后的选择，如果别的（纯Java）失败。

更新12年8月21日

看来，该问题已得到修复或Apache的公地数学的3.1-SNAPSHOT至少提高。对于上面的用例

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

更新13年2月23日

虽然乍一看这个问题，它的反应可能过于本地化，我觉得这很好说明了一些数值问题不能（有效）与解决的是什么先来到到心灵黑客的方法。所以，我希望它保持打开状态。

Answer 1:

这个问题已被固定在Apache的百科全书数学3.1.1

测试用例上述交付

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

从R-包统计的结果相匹配。 3.1-SNAPSHOT + X版本的广泛应用，也没有造成任何问题。

Answer 2:

最可能的是，这个问题一般不能被解决，因为如果累积分布函数的曲线图是非常平坦（其中它通常将朝着分布的尾部），需要在垂直轴上一个非常高的精度以达到合理精度在水平轴上。

因此，它总是会更好地使用功能直接计算位数不是从累积分布函数导出的分位数。

如果您不担心精度，就可以了，当然，解方程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]和图形是在接近该间隔的端部相当陡峭。

第二句话：你的上分位数的计算是错误的; 应该读

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

三编辑：算法修正。

Answer 3:

我发现并试图库JSci （1.2版2010年7月27日）

代码片段：

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＃逆要求更高的精度。不幸的是，即使吨迭代（100K）和10 ^ -10的算法仍然返回所请求的精度

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

现在：其代码是少错了吗？ R或JSci？我赞成一个具有更大的用户群...

文章来源: How to calculate the inverse cumulative beta distribution function in java