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