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) | / 0
where 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 are the results of the tests, where the expected results were calculated using R:
Call: ::beta::cdf-beta 2.1 3.0 0.2 Error: PASSED Result: 0.16220409270890568 Expected: 0.16220409275804 ± 1.0e-9 Call: ::beta::cdf-beta 4.2 17.3 0.5 Error: PASSED Result: 0.998630771122991 Expected: 0.998630771123192 ± 1.0e-9 Call: ::beta::cdf-beta 500 375 0.7 Error: PASSED Result: 1.0 Expected: 1.0 ± 1.0e-9 Call: ::beta::cdf-beta 250 760 0.2 Error: PASSED Result: 0.0001252342758440994 Expected: 0.000125234318666948 ± 1.0e-9 Call: ::beta::cdf-beta 43.2 19.7 0.6 Error: PASSED Result: 0.07288812943897244 Expected: 0.0728881294218269 ± 1.0e-9 Call: ::beta::cdf-beta 500 640 0.3 Error: PASSED Result: 0.0 Expected: 2.99872547567313e-23 ± 1.0e-9 Call: ::beta::cdf-beta 400 640 0.3 Error: PASSED Result: 3.040258600428558e-9 Expected: 3.07056696205524e-09 ± 1.0e-9 Call: ::beta::cdf-beta 0.1 30 0.1 Error: PASSED Result: 0.9986410086716231 Expected: 0.998641008671625 ± 1.0e-9 Call: ::beta::cdf-beta 0.01 0.03 0.9 Error: PASSED Result: 0.7658650057031781 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 Time for ::beta::cdf-beta 2 3 0.9999: 98.56 microseconds per iteration Time for ::beta::cdf-beta 250 760 0.2: 3343.65 microseconds per iteration
The time for "::beta::cdf-beta 250 760 0.2" is not good. it would take 3 seconds to run 1000 values, so performance at larger values of a & b does not appear to be very good. More efficient algorithms are available, such as the one used in DCDFLIB (http://people.scs.fsu.edu/~burkardt/f_src/dcdflib/dcdflib.html ).
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 } # Function will converge faster if x 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 z [expr {$x / (1.0 - $x)}] set retval 0.0 # Pack everything in one big exp, since individual terms can over/underflow set factor [expr {exp($lnGapb - [ln_Gamma [expr {$a + 1}]] - [ln_Gamma $b] + \ $a * $lnx + ($b - 1) * $ln1mx)}] while {$b > 1} { set retval [expr {$retval + $factor}] set a [expr {$a + 1}] set b [expr {$b - 1}] set factor [expr {$factor * $z * double($b)/double($a)}] } 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 factor 1.0 while 1 { set z [expr {$z * $x}] set abn [expr {$a + $n + $b}] set an1 [expr {$a + $n + 1}] set factor [expr {$factor * double($abn)/double($an1)}] incr n 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 set retval [expr {$retval + $pref * $sum}] # 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 }