我正在寻找一个java库/实现,它支持以合理的精度计算β分布的逆累积分布函数(也称为分位数估计)。
当然,我已经尝试过apachecommons数学,但在版本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-包装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-包统计
qbeta(c(0.025,0.975),10007+1,161750-10007+1)
[1] 0.06070355 0.06305171
为了验证R的结果,以下是Wolfram Alpha告诉我的
- InverseBeta正则化[0.02510007+1161750-10007+1]=>0.06070354631
- InverseBeta正则化[0.97510007+1161750-10007+1]=>0.06305170794
关于要求的最后说明:
- 我需要做很多这样的计算。因此,任何解决方案都不应该花费超过1的时间(与apachecommonsmath的41ms(尽管是错误的)相比,这仍然是一个很大的时间)
- 我知道在java中可以使用R。由于我在这里不详细说明的原因,如果其他任何东西(纯java)失败,这是最后一个选项
2012年8月21日更新
这个问题似乎已经在apachecommons数学的3.1-SNAPSHOT中得到了修复或至少得到了改进。对于以上的用例
2.5 percentile :0.06070354581340706
mean: 0.06187249616697166
median: 0.06187069085946604
97.5 percentile :0.06305170793994147
更新23.02.13
虽然乍一看,这个问题及其答案可能过于本地化,但我认为它很好地说明了一些数值问题无法用最初想到的黑客方法(有效)解决。所以我希望它仍然开放。
该问题已在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版本的广泛应用也没有造成任何问题。
最有可能的是,这个问题不能用一般的方法来解决,因为如果累积分布函数的图非常平坦(通常会朝向分布的尾部),则需要在垂直轴上达到非常高的精度,才能在水平轴上达到合理的精度。
因此,直接使用计算分位数的函数总是比从累积分布函数导出分位数更好。
如果你不担心精度,你当然可以用数值求解方程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 ) );
第三次编辑:已更正算法。
我已经找到并尝试了库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?我更喜欢用户群更大的。。。