home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / amiga / misc / icalc.lha / ICalc / Scripts / gamma.ic < prev    next >
Encoding:
Text File  |  1992-08-03  |  3.2 KB  |  142 lines
  1. # Various 'special' functions of mathematics and statistics
  2. # Most are based on details and algorithms given in
  3. # "Numerical Recipes - Art Of Scientific Computing" by Press, et al.
  4. #
  5. #    gamma, beta, factorial functions
  6. #    incomplete gamma functions
  7. #    error function and it's complement
  8. #    standard and generalised univariate normal distribution
  9. #
  10. # nb: routines assume array base is 1
  11. # mws, February, 1992.
  12.  
  13. echo "Defining gamma, factorial and related functions"
  14. silent
  15.  
  16. ROOT2PI = sqrt(2*PI)
  17. array gammcoef[6] = {
  18.     76.18009173, -86.50532033,
  19.     24.01409822, -1.231739516,
  20.     0.120858003e-2, -0.536382e-5
  21. }
  22.  
  23. # ln(gamma(z))
  24. func gammaln(z) = {
  25.     local tmp, ser, j
  26.     z -= 1
  27.     tmp = z+5.5; tmp = (z+0.5)*ln(tmp)-tmp
  28.     ser = 1
  29.     for (j = 1; j < 7; j += 1) {
  30.         z += 1
  31.         ser += gammcoef[j]/z
  32.     }
  33.     tmp+ln(ROOT2PI*ser)
  34. }
  35.  
  36. func gamma(z) = exp(gammaln(z))
  37.  
  38. array facttab[33] = {1}        #table for speed
  39. facttop = 0            #filled as required
  40. func fact(n) = {
  41.     local j
  42.     if (n < 0) error(1)
  43.     if (n <= facttop) {
  44.         facttab[n+1]    # implicit return (last expr evaluated)
  45.     } else if (n <= 32) {
  46.         for (j = facttop+1; j <= n; j += 1)
  47.             facttab[j+1] = j*facttab[j]
  48.         facttop = n
  49.         facttab[n+1]    # implicit return
  50.     } else gamma(n+1)
  51. }
  52.  
  53. func beta(z,w) = exp(gammaln(z) + gammaln(w) - gammaln(z+w))
  54.  
  55.  
  56. GITMAX = 100    # max iterations to calculate gamma functions
  57. GEPS = 3e-7    #
  58.  
  59. func gamser(a,x) = {        # evaluate P(a,x) by series
  60.     local gln,ap,sum,del,n
  61.  
  62.     gln = gammaln(a)
  63.     if (x == 0) return 0 else {
  64.         ap = a
  65.         sum = 1/a        
  66.         del = sum
  67.         for (n = 1; n <= GITMAX; n += 1) {
  68.             ap += 1
  69.             del *= x/ap
  70.             sum += del
  71.             if (abs(del) < abs(sum)*GEPS)
  72.                 return sum*exp(-x+a*ln(x)-gln)
  73.         }
  74.         error(2)    # can't reach desired accuracy
  75.     }
  76. }
  77.  
  78. func gamcf(a,x) = {        # evaluate Q(a,x) by continued fraction
  79.     local gln,gold,a0,a1,b0,b1,factor,n,na,nf,g
  80.  
  81.      gln = gammaln(a)
  82.     gold = 0
  83.     a0 = 1; a1 = x; b0 = 0; b1 = 1
  84.     factor = 1
  85.     for (n = 1; n <= GITMAX; n += 1) {
  86.         na = n-a
  87.         a0 = (a1+a0*na)*factor
  88.         b0 = (b1+b0*na)*factor
  89.         nf = n*factor
  90.         a1 = x*a0+nf*a1
  91.         b1 = x*b0+nf*b1
  92.         if (a1 != 0) {
  93.             factor = 1/a1
  94.             g = b1*factor
  95.             if ((abs(g-gold)/g) < GEPS)
  96.                 return exp(-x+a*ln(x)-gln)*g
  97.             gold = g
  98.         }
  99.     }
  100.     error(3)    # can't reach desired accuracy
  101. }
  102.  
  103. # incomplete gamma function P(a,x) (by standard symbology)
  104. func gammap(a,x) = {
  105.     if (x < 0 || a <= 0) error(4)    # invalid parameters
  106.     if (x < a+1) gamser(a,x) else 1-gamcf(a,x)
  107. }
  108.  
  109. # and the other, Q(a,x) = 1 - P(a,x)
  110. func gammaq(a,x) = 1-gammap(a,x)
  111.  
  112. # and the error function erf(x) and its complement erfc(x)
  113. func erf(x) = x < 0 ? -gammap(0.5,sqr(x)) : gammap(0.5,sqr(x))
  114. func erfc(x) = x < 0 ? 1+gammap(0.5,sqr(x)) : gammaq(0.5,sqr(x))
  115.  
  116. # BUT a (much) quicker calculation of erfc (and hence erf), with fractional
  117. # error less than 1.2e-7 is:
  118. func Erfc(x) = {
  119.     local z,t,rv
  120.     z = abs(x)
  121.     t = 1/(1+0.5*z)
  122.     rv = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+ \
  123.         t*(0.09678418+t*(-0.18628806+t*(0.27886807+t*(-1.13520398+ \
  124.         t*(1.48851587+t*(-0.82215223+t*0.17087277)))))))))
  125.  
  126.     x >= 0 ? rv : 2-rv
  127. }
  128.  
  129. func Erf(x) = 1-Erfc(x)
  130.  
  131. # cumulative normal distribution functions
  132. ROOT2 = sqrt(2)
  133.  
  134. # standard normal - mean 0, variance 1
  135. func snorm(x) = 1 - 0.5*Erfc(x/ROOT2)
  136.  
  137. # generalised normal - mean mu, variance sigma
  138. func gnorm(x,mu,sigma) = 1 - 0.5*Erfc((x-mu)/(sigma*ROOT2))
  139.  
  140. verbose        # turn messages back on
  141. echo "Done."
  142.