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. (This code has been tested as part of the tests run on the Beta distribution.)
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
}
# This is for convenience: incBeta_series reduces b to < 1, so better
# if b is smaller
if {$a > $b} {
return [incBeta_series $a $b $x $tol]
} else {
set z [incBeta_series $b $a [expr {1.0 - $x}] $tol]
return [expr {1.0 - $z}]
}
}
#####################################################
#
# Series expansion for Ix(a, b)
#
# Abramowitz & Stegun formula 26.5.4
#
# This series has terms
#
# [B(a+1,n)/B(a+b,n)] * x^n
#
# This is guaranteed to converge only if
# b < 1. Use the recurrence formula 26.5.15
# from A&S to bring b below 1. Recurrence is:
#
# Ix(a,b) = (G(a+b)/(G(a+1)G(b))) * x^a * (1-x)^(b-1) +
# Ix(a+1, b-1)
#
# Also, use B(a,b) = G(a)G(b)/G(a+b) to
# rewrite coeff as
#
# Cn = [G(a+1)/G(a+b)] * G(a+b+n)/G(a+n+1)
#
#####################################################
proc incBeta_series {a b x tol} {
# a+b is invariant under recurrence formula
set aplusb [expr {$a + $b}]
set lnGapb [ln_Gamma $aplusb]
# Calculate for convenience -- these don't change
set lnx [expr {log($x)}]
set ln1mx [expr {log(1.0 - $x)}]
set retval 0.0
while {$b > 1} {
# Convenient to increment a here, because expression a+1 appears in Gamma
set a [expr {$a + 1}]
# Pack everything in one big exp, since individual terms can over/underflow
set retval [expr {$retval + \
exp($lnGapb - [ln_Gamma $a] - [ln_Gamma $b] + \
($a - 1) * $lnx + ($b - 1) * $ln1mx)}]
set b [expr {$b - 1}]
}
set pref_num [expr {pow($x, $a) * pow(1.0 - $x, $b)}]
set pref_denom [expr {$a * [Beta $a $b]}]
set pref [expr {$pref_num/$pref_denom}]
set term 1.0
set z 1.0
set sum 1.0
set n 0
set adjtol [expr {$tol * $pref}]
set t1 [expr {[ln_Gamma [expr {$a + 1}]] - $lnGapb}]
while 1 {
set z [expr {$z * $x}]
incr n
set abn [expr {$a + $n + $b}]
set an1 [expr {$a + $n + 1}]
set factor [expr {exp($t1 + [ln_Gamma $abn] - [ln_Gamma $an1])}]
set nexterm [expr {$z * $factor}]
if {abs($nexterm - $term) < $adjtol} {
break
}
set term $nexterm
set sum [expr {$sum + $term}]
}
# Add previous cumulative sum to the series expansion
expr {$retval + $pref * $sum}
}