Incomplete gamma function et al. for 35s Message #1 Posted by Les Wright on 26 Aug 2007, 4:00 p.m.
Here is a long promised contribution. It is sort of long, and perhaps should eventually be reworked as an article or software library submission.
The gamma function for real numbers is, as we all know, an extension of the factorial function. Gamma(a) is the integral of exp(-t)*t^(a-1) as t goes from 0 to infinity. If we split this interval of integration at some point x, we get two incomplete gamma functions, where the left-sided gamma(a,x) takes the integral from 0 to x and the right sided Gamma(a,x) takes it from x to infinity. Obviously, Gamma(a) = gamma(a,x) + Gamma(a,x).
The incomplete gamma functions may not be of great interest to most here, but certain of their familiar special cases--the error functions and the normal and chi-square probability functions--possibly are.
The following listing computes the incomplete gamma functions gamma(a,x) and Gamma(a,x). The code is adapted from the series and continued fraction routines given in Chapter 6 of Numerical Recipes. I just give the listing here and describe use below:
G001 LBL G
G002 STO X
G003 x<>y
G004 STO A
G005 STO B
G006 1/x
G007 STO C
G008 STO D
G009 STO E
G010 1
G011 STO+ B
G012 RCL X
G013 STO* C
G014 RCL B
G015 STO/ C
G016 RCL C
G017 STO+ D
G018 STO E
G019 RCL D
G020 x#y?
G021 GTO G009
G022 GTO G067
G023 STO X
G024 x<>y
G025 STO A
G026 -
G027 1
G028 =
G029 STO B
G030 1e250
G031 STO E
G032 1/x
G033 +
G034 1/x
G035 STO C
G036 STO D
G037 CLx
G038 STO I
G039 1
G040 STO+ I
G041 STO+ B
G042 STO+ B
G043 RCL A
G044 RCL- I
G045 RCL* I
G046 STO G
G047 RCL* C
G048 RCL+ B
G049 1e-250
G050 +
G051 STO C
G052 RCL G
G053 RCL/ E
G054 RCL+ B
G055 1e-250
G056 +
G057 STO E
G058 RCL C
G059 1/x
G060 x<> C
G061 /
G062 STO* D
G063 1
G064 x#y?
G065 GTO G039
G066 RCL D
G067 RCL X
G068 RCL A
G069 y^x
G070 *
G071 RCL X
G072 +/-
G073 e^x
G074 *
G075 RTN
G076 x^2
G077 0.5
G078 x<>y
G079 XEQ G002
G080 GTO G085
G081 x^2
G082 0.5
G083 x<>y
G084 XEQ G023
G085 PI
G086 SQRT
G087 /
G088 RTN
G089 x^2
G090 2
G091 /
G092 0.5
G093 x<>y
G094 XEQ G002
G095 2
G096 /
G097 PI
G098 SQRT
G099 /
G100 0.5
G101 +
G102 RTN
G103 x^2
G104 2
G105 /
G106 0.5
G107 x<>y
G108 XEQ G023
G109 2
G110 /
G111 PI
G112 SQRT
G113 /
G114RTN
G115 XEQ G122
G116 XEQ G002
G117 GTO G120
G118 XEQ G122
G119 XEQ G023
G120 RCL/ F
G121 RTN
G122 x<>y
G123 2
G124 /
G125 ENTER
G126 ENTER
G127 1
G128 -
G129 !
G130 STO F
G131 Rv
G132 x<>y
G133 2
G134 /
G135 RTN
a ENTER x XEQ G002 gives gamma(a,x). Best used when a < x. For a > x, compute gamma(a,x) = (a-1)! - Gamma(a,x).
a ENTER x XEQ G023 gives Gamma(a,x). Use for a > x usually. For a < x, use Gamma(a,x) = (a-1)! - gamma(a,x).
x XEQ G076 gives erf(x). I use for x < 1.8, and prefer for x > 1.8 erf(x) = 1 - erfc(x).
x XEQ G081 gives erfc(x). I use it for x > 1.8, and prefer erfc(x) = 1 - erf(x) for smaller values.
z XEQ G089 gives the lower tailed normal probability associated with a standard normal variate z--i.e., the probability that a normal variate takes on a value less z. I use it for z < 2.3. For larger values, I prefer to compute 1 - the upper tailed probability.
z XEQ G103 gives the upper tailed probability associated with a standard normal variate z--i.e., the probability that variate takes on a value exceeding z. I use it for z > 2.3, and for smaller values compute 1 - lower-tailed probability.
nu ENTER x XEQ G115 gives the lower tailed probability associated with a chi-squared variate x at nu degrees of freedom. Works best for nu < x. For nu > x, compute 1 - upper-tailed probability.
nu ENTER x XEQ G118 gives the upper tailed probability associated with a chi-squared variate x at nu degrees of freedom. Use for nu > x. For nu < x compute 1 - lower-tailed probability.
Here are some examples:
5 ENTER 4 XEQ G002 gives gamma(5,4) as 8.90791255566. The last digit should be an 8.
5 ENTER 6 XEQ G023 gives Gamma(5,6) as 6.84135600761. The last digit should be a 0.
1 XEQ G076 gives erf(1) = 8.42700792945e-1. The last two digits should be 50.
2 XEQ G081 gives erfc(2) = 4.67773498102e-3. The last digit should be a 5.
1.5 XEQ G089 gives 9.33192798730e-1 as the probability that a standard normal variate will be less than 1.5. The last digit should be a 1.
2.5 XEQ G103 gives 6.20966532587e-3 as the probability that a standard normal variate will exceed 2.5. The last two digits should be 77.
4 ENTER 3 XEQ G115 gives 4.42174599629e-1 as the probability that a chisquare variate at 4 df is less than 3. All digits are correct here.
2 ENTER 7 XEQ G118 gives 3.01973834223e-2 as the probability that a chisquare variate at 2 df exceeds 7. All digits are correct here.
A few caveats are in order. First of all, I don't do sign checking, so make sure that input for all functions is positive and meaningful. I am assuming that anyone with an interest in these functions is going to be aware of the usual reflection and complementary formulae. Second, the user needs to know about the basics of the computation to judge whether one launches a routine that enters a series computation or a continued fraction computation. This is very much a matter of judgement and personal taste--I give only rough guidelines. Finally, keep in mind that due to rounding and digit loss due to subtraction in some cases it is uncommon to expect all twelve digits to be correct. Indeed, I find that in general the twelfth digit is usually off a bit, and occasionally even the eleventh. In general, I think at least 10 digits are good the vast majority of the time, at least within 1 or 2 ULP.
Sadly, these routines, while pretty fast, and certainly much faster than equivalents I have for the 41CX and 42S, are slower than on the 33S. Regrettably, they don't fit the 33S paradigm as neatly, since they would gobble up oodles of those precious labels with the various entry points and GTOs etc.
Enjoy this work, and let me know what you think.
Les
Edited: 26 Aug 2007, 4:00 p.m.
|