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

Question

我正在寻找一个java库/实现,它支持以合理的精度计算β分布的逆累积分布函数(也就是分位数的估计).

当然我已经尝试过apache commons math,但是在版本3中,似乎仍然存在一些精度问题.下面引出这个问题的问题被广泛描述.

---

假设我想通过大量试验来计算β分布的可信区间.在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- 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 [0.025,10007 + 1,161750-10007 + 1] => 0.06070354631 ...
• InverseBetaRegularized [0.975,10007 + 1,161750-10007 + 1] => 0.06305170794 ......

关于要求的最后说明:

• 我需要运行很多这些计算.因此任何解决方案都不应该花费超过1秒的时间(与41ms(虽然错误的)apache commons math相比仍然很多).
• 我知道可以在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

虽然乍一看这个问题和它的回答可能过于局部化,但我认为它很好地说明了一些数值问题无法通过首先出现在头脑中的方法来解决(有效).所以我希望它仍然开放.

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