[EKB] I've implemented the incomplete gamma function (and submitted it to tcllib). Here's the code: package require math package require math::statistics # Adapted from Fortran code in the Royal Statistical Society's StatLib # library (http://lib.stat.cmu.edu/apstat/), algorithm AS 32 (with # some modifications from AS 239) # # Calculate normalized incomplete gamma function # # 1 / x p-1 # P(p,x) = -------- | dt exp(-t) * t # Gamma(p) / 0 # # Tested some values against R's pgamma function proc incompleteGamma {x p {tol 1.0e-9}} { set overflow 1.0e37 if {$x < 0} { error "x must be positive" return } if {$p <= 0} { error "p must be greater than or equal to zero" return } # If x is zero, incGamma is zero if {$x == 0.0} { return 0.0 } # Use normal approx is p > 1000 if {$p > 1000} { set pn1 [expr {3.0 * sqrt($p) * (pow(1.0 * $x/$p, 1.0/3.0) + 1.0/(9.0 * $p) - 1.0)}] # pnorm is not robust enough for this calculation (overflows); cdf-normal could also be used return [::math::statistics::pnorm_quicker $pn1] } # If x is extremely large compared to a (and now know p < 1000), then return 1.0 if {$x > 1.e8} { return 1.0 } set factor [expr {exp($p * log($x) -$x - [::math::ln_Gamma $p])}] # Use series expansion (first option) or continued fraction if {$x <= 1.0 || $x < $p} { set gin 1.0 set term 1.0 set rn $p while {1} { set rn [expr {$rn + 1.0}] set term [expr {1.0 * $term * $x/$rn}] set gin [expr {$gin + $term}] if {$term < $tol} { set gin [expr {1.0 * $gin * $factor/$p}] break } } } else { set a [expr {1.0 - $p}] set b [expr {$a + $x + 1.0}] set term 0.0 set pn1 1.0 set pn2 $x set pn3 [expr {$x + 1.0}] set pn4 [expr {$x * $b}] set gin [expr {1.0 * $pn3/$pn4}] while {1} { set a [expr {$a + 1.0}] set b [expr {$b + 2.0}] set term [expr {$term + 1.0}] set an [expr {$a * $term}] set pn5 [expr {$b * $pn3 - $an * $pn1}] set pn6 [expr {$b * $pn4 - $an * $pn2}] if {$pn6 != 0.0} { set rn [expr {1.0 * $pn5/$pn6}] set dif [expr {abs($gin - $rn)}] if {$dif <= $tol && $dif <= $tol * $rn} { break } set gin $rn } set pn1 $pn3 set pn2 $pn4 set pn3 $pn5 set pn4 $pn6 # Too big? Rescale if {abs($pn5) >= $overflow} { set pn1 [expr {$pn1 / $overflow}] set pn2 [expr {$pn2 / $overflow}] set pn3 [expr {$pn3 / $overflow}] set pn4 [expr {$pn4 / $overflow}] } } set gin [expr {1.0 - $factor * $gin}] } return $gin } ---- !!!!!! %| enter categories here |% !!!!!!