Computes the regularized incomplete beta function, I x (a, b). The result of which is always in the range [0, 1] @param x any value in the range [0, 1] @param a any value ≥ 0 @param b any value ≥ 0 @return the result in a range of [0,1]
(double x, double a, double b)
| 563 | * @return the result in a range of [0,1] |
| 564 | */ |
| 565 | public static double betaIncReg(double x, double a, double b) |
| 566 | { |
| 567 | if(a <= 0 || b <= 0) |
| 568 | throw new ArithmeticException("a and b must be > 0, not" + a+ ", and " + b); |
| 569 | if(x == 0 || x == 1) |
| 570 | return x; |
| 571 | else if(x < 0 || x > 1) |
| 572 | throw new ArithmeticException("x must be in the range [0,1], not " + x); |
| 573 | |
| 574 | //We use this identity to make sure that our continued fraction is always in a range for which it converges faster |
| 575 | if(x > (a+1)/(a+b+2) || (1-x) < (b+1)/(a+b+2) ) |
| 576 | return 1-betaIncReg(1-x, b, a); |
| 577 | |
| 578 | |
| 579 | /* |
| 580 | * All values are from x = 0 to x = 1, in 0.025 increments |
| 581 | * a = 0.5, b = 0.5: max rel error ~ 2.2e-15 |
| 582 | * a = 0.5, b = 5: max rel error ~ 2e-15 |
| 583 | * a = 5, b = 0.5: max rel error ~ 1.5e-14 @ x ~= 7.75, otherwise rel error ~ 2e-15 |
| 584 | * a = 8, b = 10: max rel error ~ 9e-15, rel error is clearly not uniform but always small |
| 585 | * a = 80, b = 100: max rel error ~ 1.2e-14, rel error is clearly not uniform but always small |
| 586 | */ |
| 587 | |
| 588 | double numer = a*log(x)+b*log(1-x)-(log(a)+lnBeta(a, b)); |
| 589 | |
| 590 | return exp(numer)/regIncBeta.lentz(x,a,b); |
| 591 | } |
| 592 | |
| 593 | private static final Function betaIncRegFunc = new Function() { |
| 594 |