[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.16220409276911796
    Expected: 0.16220409275804 ± 1.0e-9

  Call: ::beta::cdf-beta 4.2 17.3 0.5
    Error: PASSED
     Result: 0.9986307712891124
    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.00012523431867230727
    Expected: 0.000125234318666948 ± 1.0e-9

  Call: ::beta::cdf-beta 43.2 19.7 0.6
    Error: PASSED
    Result: 0.07288812920336896
    Expected: 0.0728881294218269 ± 1.0e-9

  Call: ::beta::cdf-beta 500 640 0.3
    Error: PASSED
    Result: 2.998725475754239e-23
    Expected: 2.99872547567313e-23 ± 1.0e-9

  Call: ::beta::cdf-beta 400 640 0.3
    Error: PASSED
    Result: 3.0705669621484744e-9
    Expected: 3.07056696205524e-09 ± 1.0e-9

  Call: ::beta::cdf-beta 0.1 30 0.1
    Error: PASSED
    Result: 0.9986410087106836
    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

  Call: ::beta::cdf-beta 249.9999 759.99999 0.2
    Error: PASSED
    Result: 0.00012523707558061254
    Expected: 0.000125237075575121 ± 1.0e-9

 Time for ::beta::cdf-beta 2 3 0.9999:
 77.121 microseconds per iteration

 Time for ::beta::cdf-beta 250 760 0.2:
 1285.309 microseconds per iteration

 Time for ::beta::cdf-beta 249.9999 759.99999 0.2:
 1500.846 microseconds per iteration

These are average times over 1000 runs (using [time]). The results are OK, but not great. I've shortened them from significantly longer values in earlier tests. However, it may make sense to use more efficient algorithms, 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
        }
        
        # Adjust to get into regime where convergence is faster
        if {$x <= 0.5} {
            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 adjtol [expr {$tol * $pref}]
        if {$adjtol == 0.0} {
            # Is pref so small that it is evaluated at exactly zero?
            set sum 0.0
        } elseif {[string is integer $b] && $b == 1} {
            # In the case of integer b, this reduces to 1 + x + x^2 + ... = 1/(1-x)
            # From routine above, b must be 1, but check
            set sum [expr {1.0/(1.0 - $x)}]
        } else {
            set term 1.0
            set z 1.0
            set sum 1.0
            set n 0
            
            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
        
    }
----
!!!!!!
%| [Category Mathematics] | [Category Statistics] |%
!!!!!!