home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Special Functions < prev   
Encoding:
Text File  |  1992-09-03  |  6.6 KB  |  211 lines  |  [TEXT/NCII]
  1. ################
  2. # Special Functions
  3. #   (Also see the file Fitting Functions to Arrays, which
  4. #     has definitions for Legendre Polynomials.)
  5. # 
  6. # 
  7. # function Gamma(x) # Returns  integral 0->∞ { exp(-t) t^(x-1) dt .
  8. # function IncompleteGamma(a,x)  # Returns P(a,x)= 1/gamma(x)  integral 0->x { exp(-t) t^(a-1) dt }
  9. # 
  10. # function ErrorFunction(x)  # Return erf(x).
  11. # function ChiSqProbability(ChiSq,Nu)  # Returns prob. of ChiSqaured fit.
  12. # 
  13. # function Bessel(L,x) # Returns the spherical Bessel functions j_l(x), l=1,2,3,…
  14. # 
  15. # 
  16. #  This text file explains and gives examples of
  17. #  some of the routines in the file All Library Routines, which
  18. #  should be Opened before trying any of these examples.
  19. #################
  20.  
  21. # Here are the xCOD's that are used.
  22.  
  23. xGAMMA
  24.   # xCOD xGAMMA(N,x,Answer:num), finds Answer[1…N]=Gamma(x[1…N])=integral 0->∞ of t^(x-1) e^(-t) dt, with Gamma(J)=(J-1)! for integer J ; an external program.
  25.  
  26. xINCOMPLETE_GAMMA
  27.   # xCOD xINCOMPLETE_GAMMA(N,a,X,Answer,Err:num), calculates Answer[1…N]=1/gamma(a) integral 0 to x[1…N] {exp(-t) t^(a-1)}. Err: 0=> success, -1=> hit iteration limit, -2=> x≤0.; an external program.
  28.  
  29. xBESSEL
  30.   # xCOD xBESSEL(L,N,X,Answer,err:num); computes a spherical Bessel function Answer[1…N]=j_L(X[1…N]). Err: 0=> success, -1=> L < 0.; an external program.
  31.  
  32.  
  33. ####
  34. # Interface routines:
  35. #
  36.  
  37. # The gamma function is Gamma(x) = integral 0->∞ { exp(-t) t^(x-1) dt},
  38. # which is Gamma(n) = (n-1)! for integer n.
  39. #
  40. function Gamma(x) # Returns  integral 0->∞ { exp(-t) t^(x-1) dt .
  41. .   var n,y
  42. . #   Input:     x = number or array
  43. . #   Output:   Gamma(x) = number or array = G(x)
  44. . #                     = integral 0->∞ { exp(-t) t^(x-1) dt
  45. .   begin
  46. .     n = size(x)
  47. .     y = x # save space for answer
  48. .     xGAMMA(n,x,y)
  49. .     Gamma = y
  50. .   end
  51. .   
  52.  
  53. # This incomplete gamma function (sometimes called P(a,x)) is
  54. #    1/gamma(x)  integral 0->x { exp(-t) t^(a-1) dt }
  55. # and has IncompleteGamma(a,0) = 0 and 
  56. # IncompleteGamma(1,∞)=1. 
  57. #   The error function and the goodness of a Chi Squared 
  58. # fit may be calculated in terms of this.
  59. #
  60. function IncompleteGamma(a,x)  # Returns P(a,x)= 1/gamma(x)  integral 0->x { exp(-t) t^(a-1) dt }
  61. .   var n,y, err
  62. . #     Input:     a = number > 0
  63. . #                   x = number or array
  64. . #     Output:
  65. . #                   IncompleteGamma = P(a,x) = number or array
  66. . #                     = 1/gamma(x)  integral 0->x { exp(-t) t^(a-1) dt }
  67. .   begin
  68. .     n=size(x)
  69. .     y=x  # save space for answer
  70. .     xINCOMPLETE_GAMMA(n,a,x,y,err)
  71. .     if err=0 then
  72. .       IncompleteGamma = y
  73. .     else if -2 then
  74. .       IncompleteGamma = " ERROR in xINCOMPLETE_GAMMA, x<0 "
  75. .     else
  76. .       IncompleteGamma = " ERROR in xINCOMPLETE_GAMMA, did not converge"
  77. .   end
  78. .  
  79.  
  80. # the Error Function is erf(x) = 2/sqrt(π) integral 0->x { exp(-t^2) dt }
  81. # which is the integrated probality of a Gaussian normal distribution 
  82. # with sigma=1/sqrt(2) from -x to x.
  83. #
  84. #
  85. function ErrorFunction(x)  # Return erf(x).
  86. .   #    Input:  x = number or array
  87. .   #    Output  ErrorFunction = erf(x) = number or array
  88. .   #                   = 2/sqrt(π) integral 0->x { exp(-t^2) dt }
  89. .   begin
  90. .     if x<0 then
  91. .       ErrorFunction = -IncompleteGamma(0.5, x^2)
  92. .     else
  93. .       ErrorFunction = IncompleteGamma(0.5, x^2)
  94. .     end
  95. .   end
  96. .   
  97.  
  98. #
  99. # This returns the probability that a certain model with
  100. # Nu degress of freedom and a given ChiSq = c2 could happen
  101. # by chance. 
  102. #
  103. function ChiSqProbability(ChiSq,Nu)  # Returns prob. of ChiSqaured fit.
  104. .  # Input: ChiSq = sum [(yData[1…N]-yModel(xData))/sigma]^2
  105. .  #          Nu = N data points - J model params = degrees of freedom
  106. .  # Output: ChiSqProbability = prob of a fit this bad being due
  107. .  #              to Gaussian normal errors with stnd. dev. Sigma[1…N].
  108. .  #     (ChiSqProbabilitity << 1 means the fit is bad.)
  109. .  begin
  110. .      ChiSqProbability = 1 - IncompleteGamma(Nu/2, ChiSq/2);
  111. .   end
  112. .   
  113.  
  114.  
  115.  
  116. function Bessel(L,x) # Returns the spherical Bessel functions j_l(x), l=1,2,3,…
  117. .   var n,y, err
  118. . #   Input:    
  119. . #                     L = positive integer
  120. . #                     x = number or array
  121. . #   Output:
  122. . #                     Bessel = j_l(x) = number or array
  123. .   begin
  124. .     n = size(x)
  125. .     y = x # save space for answer
  126. .     xBESSEL(L,n,x,y,err)
  127. .     if err<>0 then
  128. .       Bessel = "ERROR: L must be >0 in Bessel."
  129. .     else
  130. .       Bessel = y
  131. .   end
  132. .   
  133.  
  134.  
  135. #######
  136. # Examples
  137. ########
  138.  
  139. # A Bessel function:
  140. # For a single value,
  141.   Bessel(2, 3.2)
  142.   0.30536678    #  this is j_2(3.2)
  143. #
  144. #     The easiest way to see the whole function is to just
  145. #     plot   [x, Bessel(L,x)]   for L=1,2,3, etc.
  146. # Or, for a short little table,
  147. x = [0, 0.5, 1, 1.5, 2];   j_1 = Bessel(1,x);
  148. #  table(x,j_1)       # evaluate the "table" command to get:
  149. # x                  j_1
  150. # - - - - - - - - - - - - - - - 
  151.    0                  0
  152.    0.5               0.16254
  153.    1                  0.30117
  154.    1.5               0.39617
  155.    2                  0.4354
  156.  
  157.  
  158. # The Gamma function:
  159. Gamma(5) 
  160.   24 
  161. # This is the same as 4!=4*3*2*1
  162. 4!
  163.   24 
  164. #  whilel
  165. Gamma(0.5)
  166.   1.77245385 
  167. # is the same as sqrt(π)
  168. sqrt(π)
  169.   1.77245385 
  170. # Again, you can just plot
  171. #      [x, Gamma(x)] 
  172. #  to see it.  Note that it blows up at x=0, -1, -2, ….
  173.  
  174.  
  175. # The error function is usually used to describe the probability that
  176. # a normal variable is within some number of standard deviations
  177. # from the mean.  This is often used in scientific data to describe
  178. # how accurate a result is : "one sigma" means a result which, 
  179. # on the average, woule happen as often as a normal variable
  180. # would fall within 1 standard deviation of the mean.
  181. #
  182. # The sigma for the error function is 1/sqrt(2):
  183. sigma = 1/sqrt(2)
  184.   sigma = 0.70710678 
  185. #
  186. #  Then
  187. ErrorFunction(1 sigma)
  188.   0.68268948                   # 68% (~2/3) chance of being ± 1 sigma from mean.
  189. ErrorFunction(2 sigma)
  190.   0.95449972                   # 95% chance of being with ± 2 sigma
  191. ErrorFunction(3 sigma)
  192.   0.9973002                     # 99.7% chance of being with ± 3 sigma
  193.  
  194.  
  195.  
  196. # Usually the agreement between a model G(x) with Nu degrees of freedom 
  197. # and a data set xData, yData[1…n] with standard deviations is given by 
  198. # ChiSquared = sum over n { [yData - G(xData)]/sigma }^2.
  199. # Then the probability that that the ChiSquared errors could be as large 
  200. # as they are by pure chance, assuming that they are distributed 
  201. # normally with standard deviations sigma[1…n], is given by
  202. #  ChiSqProbability(ChiSq, Nu), given above.
  203. # 
  204. # If , for instance, a certain data model with 3 free parameters is 
  205. # supposed to model a certain data set which has a ChiSq errors of
  206. # 8.3, then the probability that this could have happened by chance is
  207. ChiSqProbability(8.3, 3)
  208.   0.0402019 
  209. # or 4%.  Thus this model is not a good fit to the data.
  210.  
  211.