[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].

''[EKB] 11 Jan 2007: Significantly improved on the code by replacing the series expansion with a continued fraction implementation. Running times are shorter and more consistent with different parameters.''

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]:

 23 tests run
 23 tests passed
 100.0% pass rate

 All tests:
  Call: ::beta::pdf-beta 1.3 2.4 0.2
    Error: PASSED
    Result: 1.689031804741256
    Expected: 1.68903180472449 ± 1.0e-9

  Call: ::beta::pdf-beta 1 1 0.5
    Error: PASSED
    Result: 0.9999999999900844
    Expected: 1.0 ± 1.0e-9

  Call: ::beta::pdf-beta 3.7 0.9 0.0
    Error: PASSED
    Result: 0.0
    Expected: 0.0 ± 1.0e-9

  Call: ::beta::pdf-beta 1.8 4.2 1.0
    Error: PASSED
    Result: 0.0
    Expected: 0.0 ± 1.0e-9

  Call: ::beta::pdf-beta 320 400 0.4
    Error: PASSED
    Result: 1.1819237678887626
    Expected: 1.18192376783860 ± 1.0e-9

  Call: ::beta::pdf-beta 500 1 0.2
    Error: PASSED
    Result: 0.0
    Expected: 0.0 ± 1.0e-9

  Call: ::beta::pdf-beta 1000 1000 0.50
    Error: PASSED
    Result: 35.678022292091086
    Expected: 35.6780222917086 ± 1.0e-9

  Call: ::beta::cdf-beta 2.1 3.0 0.2
    Error: PASSED
    Result: 0.16220409276377096
    Expected: 0.16220409275804 ± 1.0e-9

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

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

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

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

  Call: ::beta::cdf-beta 0.1 30 0.1
    Error: PASSED
    Result: 0.998641008647038
    Expected: 0.998641008671625 ± 1.0e-9

  Call: ::beta::cdf-beta 0.01 0.03 0.9
    Error: PASSED
    Result: 0.7658650057014063
    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.00012523707557178366
    Expected: 0.000125237075575121 ± 1.0e-9

  Call: ::beta::cdf-beta 1000 1000 0.4
    Error: PASSED
    Result: 8.23161135455908e-20
    Expected: 8.23161135486914e-20 ± 1.0e-9

  Call: ::beta::cdf-beta 1000 1000 0.499
    Error: PASSED
    Result: 0.46436944398866614
    Expected: 0.464369443974288 ± 1.0e-9

  Call: ::beta::cdf-beta 1000 1000 0.5
    Error: PASSED
    Result: 0.4999999999893452
    Expected: 0.5 ± 1.0e-9

  Call: ::beta::cdf-beta 1000 1000 0.7
    Error: PASSED
    Result: 1.0
    Expected: 1.0 ± 1.0e-9

  Call: ::beta::cdf-beta 2 3 0.6
    Error: PASSED
    Result: 0.8207999999940695
    Expected: 0.8208 ± 1.0e-9

  Average time/1000 iterations for ::beta::cdf-beta 2 3 0.9999:
  51.041 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 250 760 0.47:
  87.181 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 249.9999 759.99999 0.47:
  106.473 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 2.1 3.2 0.7:
  74.872 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 4.3 9.2 0.3:
  91.824 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 2.5 1.9 0.1:
  63.872 microseconds per iteration
  Average time/1000 iterations for ::beta::cdf-beta 5 7 0.3:
  75.041 microseconds per iteration

The running times are average times over 1000 runs (using [time]).

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
        }
        
        # Rearrange if necessary to get continued fraction to behave
        if {$x < 0.5} {
            return [beta_cont_frac $a $b $x $tol]
        } else {
            set z [beta_cont_frac $b $a [expr {1.0 - $x}] $tol]
            return [expr {1.0 - $z}]
        }
    }

    #####################################################
    #
    # Continued fraction for Ix(a,b)
    #
    # Abramowitz & Stegun 26.5.9
    #
    #####################################################
    proc beta_cont_frac {a b x {tol 1.0e-9}} {
        set max_iter 512
        
        set aplusb [expr {$a + $b}]
        set amin1 [expr {$a - 1}]
        set lnGapb [ln_Gamma $aplusb]
        set term1 [expr {$lnGapb- [ln_Gamma $a] - [ln_Gamma $b]}]
        set term2 [expr {$a * log($x) + ($b - 1.0) * log(1.0 - $x)}]
        set pref [expr {exp($term1 + $term2)/$a}]
        
        set z [expr {$x / (1.0 - $x)}]
        
        set v 1.0
        set h_1 1.0
        set h_2 0.0
        set k_1 1.0
        set k_2 1.0
        
        for {set m 1} {$m < $max_iter} {incr m} {
            set f1 [expr {$amin1 + 2 * $m}]
            set e2m [expr {-$z * double(($amin1 + $m) * ($b - $m))/ \
                double(($f1 - 1) * $f1)}]
            set e2mp1 [expr {$z * double($m * ($aplusb - 1 + $m)) / \
                double($f1 * ($f1 + 1))}]
            set h_2m [expr {$h_1 + $e2m * $h_2}]
            set k_2m [expr {$k_1 + $e2m * $k_2}]
            
            set h_2 $h_2m
            set k_2 $k_2m
            
            set h_1 [expr {$h_2m + $e2mp1 * $h_1}]
            set k_1 [expr {$k_2m + $e2mp1 * $k_1}]
            
            set vprime [expr {$h_1/$k_1}]
            
            if {abs($v - $vprime) < $tol} {
                break
            }
            
            set v $vprime
            
        }
        
        if {$m == $max_iter} {
            error "beta_cont_frac: Exceeded maximum number of iterations"
        }
        
        set retval [expr {$pref * $v}]
        
        # 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] |%
!!!!!!