如何在Java中计算反累积Beta分布函数

Question

我正在寻找一个支持计算Beta分布的反累积分布函数（也称为量化估计）的Java库/实现，具有合理的精度。当然，我已经尝试过Apache Commons Math，但在版本3中仍存在一些精度问题。下面详细描述了导致这个问题的情况。假设我想计算具有大量试验的Beta分布的可信区间，在apache commons math中...

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 语言中的 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秒（尽管与（虽然错误的）apache commons math的41ms相比还是很多）。
• 我知道可以在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

虽然乍一看这个问题及其回复可能过于具体，但我认为它很好地说明了有些数值问题不能用头脑风暴式的方法来高效地解决。因此，我希望它能保持开放。

- steffen

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包stats的结果相匹配。广泛应用3.1-SNAPSHOT + x版本也没有引起任何问题。

- steffen

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]

备注：我不清楚为什么精度在beta分布中应该是一个问题，因为beta分布存在于区间[0;1]上，而且图形在区间的两端非常陡峭。

第二个备注：您计算的上分位数是错误的；应该是

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

第三次编辑：算法已经修正。

- JohnB

---

非常感谢您的回复。请问您能否提供 a) 您方法的参考资料（名称、链接等均可）和 b) 一个调用示例（使用易于访问的 F 库）？ - steffen

它简单地使用嵌套区间来数值地找到函数的零点。 - JohnB

谢谢您指出错误。我已经更新了问题，并提供了另一个使用情况，其中Apache Commons Math出现错误。 - steffen

尝试通过向BetaDistribution的构造函数添加第三个参数来提高准确性，例如：BetaDistribution betaDist = new BetaDistribution(successes + 1, trials - successes + 1, 0.000001); - JohnB

谢谢您的评论。然而，根据源代码，缺省精度为1e-9并没有多大帮助。不幸的是，将其增加到10^-20也没有帮助。 - steffen

我尝试了您的方法，但对于精度=10^-5，q=0.025并使用colt库中的Beta-cdf实现，该算法在30分钟内未能终止。我可能做错了什么，因此我谦虚地请求一段代码片段。 - steffen

Answer 3

我已经找到并尝试了库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？我会支持用户群更大的那一个...

- steffen