[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
            }
            #tk_messageBox -message "$x : $term : $nexterm"
            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
        
    }
----
!!!!!!
%| [Category Mathematics] | [Category Statistics] |%
!!!!!!