java - 如何计算java中的逆累积贝塔分布函数

Question

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

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

---

假设我想通过大量试验计算 beta 分布的可信区间。在 apache commons 数学 ...

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 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 数学的 41 毫秒相比，这仍然很多)。
• 我知道可以在 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 统计的结果相匹配。 3.1-SNAPSHOT + x版本的广泛应用也没有造成任何问题。