[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'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}
        
    }
----
!!!!!!
%| [Category Mathematics] | [Category Statistics] |%
!!!!!!