如何在java中计算逆累积beta分布函数

Question

我正在寻找一个 java 库/实现，它支持以合理的精度计算 beta 分布的逆累积分布函数（也称为分位数估计）。

当然我已经尝试过apache commons math，但在版本 3 中似乎还有一些issues with the precision。下面对导致这个问题的问题进行了广泛的描述。

---

假设我想通过大量试验计算 beta 分布的可信区间。在 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-包stats

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

为了支持 R 的结果，这是 Wolfram Alpha 告诉我的

• InverseBetaRegularized[0.025,10007+1,161750-10007+1] => 0.06070354631...
• InverseBetaRegularized[0.975,10007+1,161750-10007+1] => 0.06305170794...

关于要求的最后说明：

• 我需要运行很多这样的计算。因此，任何解决方案所花费的时间都不应超过 1 秒（与 41 毫秒（尽管是错误的）apache commons 数学相比，这仍然很多）。
• 我知道可以在 java 中使用 R。由于我不会在这里详细说明的原因，如果其他任何东西（纯 java）失败，这是最后一个选项。

更新 21.08.12

It seems 表示该问题已在 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

虽然乍一看这个问题及其回答可能过于本地化，但我认为它很好地说明了一些数字问题无法（有效地）通过首先想到的黑客方法来解决。所以我希望它保持开放。

Answer 1

该问题已在 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 版本的广泛应用也没有造成任何问题。

Answer 2

我已经找到并尝试了库 JSci（版本 1.2 27.07.2010）

代码sn-p：

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 ？我更喜欢拥有更大用户群的那个...

Answer 3

这个问题很可能无法通用解决，因为如果累积分布函数的图形非常平坦（通常朝向分布的尾部），则需要在垂直轴上具有非常高的精度在水平轴上达到合理的精度。

因此，使用直接计算分位数的函数总是比从累积分布函数导出分位数更好。

如果您不担心精度，您当然可以数值求解方程 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]

备注：我不清楚为什么精度应该是一个问题，特别是对于 beta 分布，因为 beta 分布存在于区间 [0;1] 上，而且图表相当陡峭间隔的结束。

第二句话：你计算的上分位数有误；它应该是

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

第三次修改：算法已更正。