EKB This is an implementation of the incomplete Beta function Ix(a, b), defined as:
/ x
1 | a-1 b-1
Ix(a,b) = ------- | dt t (1 - t)
B(a, b) |
/ 0where B(a,b) is the Beta function. The incomplete Beta function is the cumulative probability function for the Beta distribution.
EKB 11 Jan 2007: Significantly improved on the code by replacing the series expansion with a continued fraction implementation. Running times are shorter and more consistent with different parameters.
This code has been tested as part of the tests run on the Beta distribution. Here are the results of the tests, where the expected results were calculated using R:
23 tests run
23 tests passed
100.0% pass rate
All tests:
Call: ::beta::pdf-beta 1.3 2.4 0.2
Error: PASSED
Result: 1.689031804741256
Expected: 1.68903180472449 ± 1.0e-9
Call: ::beta::pdf-beta 1 1 0.5
Error: PASSED
Result: 0.9999999999900844
Expected: 1.0 ± 1.0e-9
Call: ::beta::pdf-beta 3.7 0.9 0.0
Error: PASSED
Result: 0.0
Expected: 0.0 ± 1.0e-9
Call: ::beta::pdf-beta 1.8 4.2 1.0
Error: PASSED
Result: 0.0
Expected: 0.0 ± 1.0e-9
Call: ::beta::pdf-beta 320 400 0.4
Error: PASSED
Result: 1.1819237678887626
Expected: 1.18192376783860 ± 1.0e-9
Call: ::beta::pdf-beta 500 1 0.2
Error: PASSED
Result: 0.0
Expected: 0.0 ± 1.0e-9
Call: ::beta::pdf-beta 1000 1000 0.50
Error: PASSED
Result: 35.678022292091086
Expected: 35.6780222917086 ± 1.0e-9
Call: ::beta::cdf-beta 2.1 3.0 0.2
Error: PASSED
Result: 0.16220409276377096
Expected: 0.16220409275804 ± 1.0e-9
Call: ::beta::cdf-beta 4.2 17.3 0.5
Error: PASSED
Result: 0.998630771122973
Expected: 0.998630771123192 ± 1.0e-9
Call: ::beta::cdf-beta 500 375 0.7
Error: PASSED
Result: 0.9999999999999996
Expected: 1.0 ± 1.0e-9
Call: ::beta::cdf-beta 250 760 0.2
Error: PASSED
Result: 0.000125234318663519
Expected: 0.000125234318666948 ± 1.0e-9
Call: ::beta::cdf-beta 43.2 19.7 0.6
Error: PASSED
Result: 0.07288812920992815
Expected: 0.0728881294218269 ± 1.0e-9
Call: ::beta::cdf-beta 500 640 0.3
Error: PASSED
Result: 2.9987254761180083e-23
Expected: 2.99872547567313e-23 ± 1.0e-9
Call: ::beta::cdf-beta 400 640 0.3
Error: PASSED
Result: 3.070566962863235e-9
Expected: 3.07056696205524e-09 ± 1.0e-9
Call: ::beta::cdf-beta 0.1 30 0.1
Error: PASSED
Result: 0.998641008647038
Expected: 0.998641008671625 ± 1.0e-9
Call: ::beta::cdf-beta 0.01 0.03 0.9
Error: PASSED
Result: 0.7658650057014063
Expected: 0.765865005703006 ± 1.0e-9
Call: ::beta::cdf-beta 2 3 0.9999
Error: PASSED
Result: 0.9999999999960003
Expected: 0.999999999996 ± 1.0e-9
Call: ::beta::cdf-beta 249.9999 759.99999 0.2
Error: PASSED
Result: 0.00012523707557178366
Expected: 0.000125237075575121 ± 1.0e-9
Call: ::beta::cdf-beta 1000 1000 0.4
Error: PASSED
Result: 8.23161135455908e-20
Expected: 8.23161135486914e-20 ± 1.0e-9
Call: ::beta::cdf-beta 1000 1000 0.499
Error: PASSED
Result: 0.46436944398866614
Expected: 0.464369443974288 ± 1.0e-9
Call: ::beta::cdf-beta 1000 1000 0.5
Error: PASSED
Result: 0.4999999999893452
Expected: 0.5 ± 1.0e-9
Call: ::beta::cdf-beta 1000 1000 0.7
Error: PASSED
Result: 1.0
Expected: 1.0 ± 1.0e-9
Call: ::beta::cdf-beta 2 3 0.6
Error: PASSED
Result: 0.8207999999940695
Expected: 0.8208 ± 1.0e-9
Average time/1000 iterations for ::beta::cdf-beta 2 3 0.9999:
51.041 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 250 760 0.47:
87.181 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 249.9999 759.99999 0.47:
106.473 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 2.1 3.2 0.7:
74.872 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 4.3 9.2 0.3:
91.824 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 2.5 1.9 0.1:
63.872 microseconds per iteration
Average time/1000 iterations for ::beta::cdf-beta 5 7 0.3:
75.041 microseconds per iterationThe running times are average times over 1000 runs (using time).
Here's the code:
package require math
namespace import ::math::ln_Gamma
namespace import ::math::Beta
#
# Implement the incomplete beta function Ix(a, b)
#
proc incompleteBeta {a b x {tol 1.0e-9}} {
if {$x < 0.0 || $x > 1.0} {
error "Value out of range in incomplete Beta function: x = $x, not in \[0, 1\]"
}
if {$a <= 0.0} {
error "Value out of range in incomplete Beta function: a = $a, must be > 0"
}
if {$b <= 0.0} {
error "Value out of range in incomplete Beta function: b = $b, must be > 0"
}
if {$x < $tol} {
return 0.0
}
if {$x > 1.0 - $tol} {
return 1.0
}
# Rearrange if necessary to get continued fraction to behave
if {$x < 0.5} {
return [beta_cont_frac $a $b $x $tol]
} else {
set z [beta_cont_frac $b $a [expr {1.0 - $x}] $tol]
return [expr {1.0 - $z}]
}
}
#####################################################
#
# Continued fraction for Ix(a,b)
#
# Abramowitz & Stegun 26.5.9
#
#####################################################
proc beta_cont_frac {a b x {tol 1.0e-9}} {
set max_iter 512
set aplusb [expr {$a + $b}]
set amin1 [expr {$a - 1}]
set lnGapb [ln_Gamma $aplusb]
set term1 [expr {$lnGapb- [ln_Gamma $a] - [ln_Gamma $b]}]
set term2 [expr {$a * log($x) + ($b - 1.0) * log(1.0 - $x)}]
set pref [expr {exp($term1 + $term2)/$a}]
set z [expr {$x / (1.0 - $x)}]
set v 1.0
set h_1 1.0
set h_2 0.0
set k_1 1.0
set k_2 1.0
for {set m 1} {$m < $max_iter} {incr m} {
set f1 [expr {$amin1 + 2 * $m}]
set e2m [expr {-$z * double(($amin1 + $m) * ($b - $m))/ \
double(($f1 - 1) * $f1)}]
set e2mp1 [expr {$z * double($m * ($aplusb - 1 + $m)) / \
double($f1 * ($f1 + 1))}]
set h_2m [expr {$h_1 + $e2m * $h_2}]
set k_2m [expr {$k_1 + $e2m * $k_2}]
set h_2 $h_2m
set k_2 $k_2m
set h_1 [expr {$h_2m + $e2mp1 * $h_1}]
set k_1 [expr {$k_2m + $e2mp1 * $k_1}]
set vprime [expr {$h_1/$k_1}]
if {abs($v - $vprime) < $tol} {
break
}
set v $vprime
}
if {$m == $max_iter} {
error "beta_cont_frac: Exceeded maximum number of iterations"
}
set retval [expr {$pref * $v}]
# Because of imprecision in underlying Tcl calculations, may fall out of bounds
if {$retval < 0.0} {
set retval 0.0
} elseif {$retval > 1.0} {
set retval 1.0
}
return $retval
}