home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Frozen Fish 2: PC
/
frozenfish_august_1995.bin
/
bbs
/
d07xx
/
d0713.lha
/
ICalc
/
Scripts
/
gamma.ic
< prev
next >
Wrap
Text File
|
1992-08-19
|
3KB
|
142 lines
# Various 'special' functions of mathematics and statistics
# Most are based on details and algorithms given in
# "Numerical Recipes - Art Of Scientific Computing" by Press, et al.
#
# gamma, beta, factorial functions
# incomplete gamma functions
# error function and it's complement
# standard and generalised univariate normal distribution
#
# nb: routines assume array base is 1
# mws, February, 1992.
echo "Defining gamma, factorial and related functions"
silent
ROOT2PI = sqrt(2*PI)
array gammcoef[6] = {
76.18009173, -86.50532033,
24.01409822, -1.231739516,
0.120858003e-2, -0.536382e-5
}
# ln(gamma(z))
func gammaln(z) = {
local tmp, ser, j
z -= 1
tmp = z+5.5; tmp = (z+0.5)*ln(tmp)-tmp
ser = 1
for (j = 1; j < 7; j += 1) {
z += 1
ser += gammcoef[j]/z
}
tmp+ln(ROOT2PI*ser)
}
func gamma(z) = exp(gammaln(z))
array facttab[33] = {1} #table for speed
facttop = 0 #filled as required
func fact(n) = {
local j
if (n < 0) error(1)
if (n <= facttop) {
facttab[n+1] # implicit return (last expr evaluated)
} else if (n <= 32) {
for (j = facttop+1; j <= n; j += 1)
facttab[j+1] = j*facttab[j]
facttop = n
facttab[n+1] # implicit return
} else gamma(n+1)
}
func beta(z,w) = exp(gammaln(z) + gammaln(w) - gammaln(z+w))
GITMAX = 100 # max iterations to calculate gamma functions
GEPS = 3e-7 #
func gamser(a,x) = { # evaluate P(a,x) by series
local gln,ap,sum,del,n
gln = gammaln(a)
if (x == 0) return 0 else {
ap = a
sum = 1/a
del = sum
for (n = 1; n <= GITMAX; n += 1) {
ap += 1
del *= x/ap
sum += del
if (abs(del) < abs(sum)*GEPS)
return sum*exp(-x+a*ln(x)-gln)
}
error(2) # can't reach desired accuracy
}
}
func gamcf(a,x) = { # evaluate Q(a,x) by continued fraction
local gln,gold,a0,a1,b0,b1,factor,n,na,nf,g
gln = gammaln(a)
gold = 0
a0 = 1; a1 = x; b0 = 0; b1 = 1
factor = 1
for (n = 1; n <= GITMAX; n += 1) {
na = n-a
a0 = (a1+a0*na)*factor
b0 = (b1+b0*na)*factor
nf = n*factor
a1 = x*a0+nf*a1
b1 = x*b0+nf*b1
if (a1 != 0) {
factor = 1/a1
g = b1*factor
if ((abs(g-gold)/g) < GEPS)
return exp(-x+a*ln(x)-gln)*g
gold = g
}
}
error(3) # can't reach desired accuracy
}
# incomplete gamma function P(a,x) (by standard symbology)
func gammap(a,x) = {
if (x < 0 || a <= 0) error(4) # invalid parameters
if (x < a+1) gamser(a,x) else 1-gamcf(a,x)
}
# and the other, Q(a,x) = 1 - P(a,x)
func gammaq(a,x) = 1-gammap(a,x)
# and the error function erf(x) and its complement erfc(x)
func erf(x) = x < 0 ? -gammap(0.5,sqr(x)) : gammap(0.5,sqr(x))
func erfc(x) = x < 0 ? 1+gammap(0.5,sqr(x)) : gammaq(0.5,sqr(x))
# BUT a (much) quicker calculation of erfc (and hence erf), with fractional
# error less than 1.2e-7 is:
func Erfc(x) = {
local z,t,rv
z = abs(x)
t = 1/(1+0.5*z)
rv = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+ \
t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+ \
t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))))
x >= 0 ? rv : 2-rv
}
func Erf(x) = 1-Erfc(x)
# cumulative normal distribution functions
ROOT2 = sqrt(2)
# standard normal - mean 0, variance 1
func snorm(x) = 1 - 0.5*Erfc(x/ROOT2)
# generalised normal - mean mu, variance sigma
func gnorm(x,mu,sigma) = 1 - 0.5*Erfc((x-mu)/(sigma*ROOT2))
verbose # turn messages back on
echo "Done."