[wp34s] Incomplete Gamma on the wp34s  Printable Version + HP Forums (https://archived.hpcalc.org/museumforum) + Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum1.html) + Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum2.html) + Thread: [wp34s] Incomplete Gamma on the wp34s (/thread206268.html) 
[wp34s] Incomplete Gamma on the wp34s  Les Wright  12052011 Thought there would be some interest in this article: For reasons given at the end of the article, this routine may not be terribly useful to most, but I must say it is a programming job I have been pretty pleased with over the years and I really enjoyed porting it to the wp34s.
Les
Re: [wp34s] Incomplete Gamma on the wp34s  Paul Dale  12052011 The 34S already includes a regularised incomplete gamma functions for real arguments. The difference between this and the incomplete gamma function is a scaling by a normal gamma function evaluation.
Re: [wp34s] Incomplete Gamma on the wp34s  Les Wright  12052011 Dang, I missed that one completely...
The stuff you guys have packed in... :D
Re: [wp34s] Incomplete Gamma on the wp34s  Les Wright  12052011 I should point out that the documentation of the IGamma function is actually incorrect. It returns IGamma(y,x). The arguments are are reversed.
Les
Re: [wp34s] Incomplete Gamma on the wp34s  Les Wright  12052011 Quote: Likewise with the incomplete beta function. It actually returns IBeta(x,z,y). Les
Re: [wp34s] Incomplete Gamma on the wp34s  Paul Dale  12052011 The documentation will be updated in due course I suspect :) I must admit, I never really tested this function in isolation. It is used by some of the probability distributions and was tested as part of them.
Maybe Not Completely Redundant  Les Wright  12052011 Actually, the wp34s provides only the lower regularized incomplete gamma and cdf's (lower tail) for phi and chisquare. But what if one wants high precision in the upper tail of the regularized incomplete gamma or its specific relatives? For phi it's easy, since the symmetry of the standard normal distribution means that 1  phi(x) = phi(x), so if I want the uppertail probability associated with the normal statistic x = 3 I would compute phi(3) instead. But there is no such reflection relation for chisquare. A highvalue of the statistic will give a probability close to one, and digits (and precision) will be lost due to subtraction from 1 to get the complementary uppertail probability. It is probably better to compute the uppertail chisquare probability by way of the underlying continued fraction expansion for the rightsided incomplete gammaIGL in my program. Of course, this only matters if one really cares about the lost digits. In practical applications, most likely wouldn't.
I notice the builtin erfc function gives all or almost all 16 digits accurate for higher values of the argument in cases where subtracting erf from 1 would give zero. I am surmising that erfc has it's own routine and is not derived from the builtin regularized incomplete gamma.
Re: Maybe Not Completely Redundant  Walter B  12052011 Anybody doing statistics with more than a few significant digits shall check whether his mathematical model describes the reality observed to this precision  and if true then please check whether it's necessary. Mathematical statistics provides tools for better estimations, nothing more. I've written this earlier already, presumably.
Re: [wp34s] Incomplete Gamma on the wp34s  Namir  12052011 Les, Can you please provide a link to the algorithms you used to calculate the two versions of the incomplete gamma functions?
Namir
Re: Maybe Not Completely Redundant  Les Wright  12052011 True dat, as the young people like to say today. Obviously, my interest in computing the uppertail functions to high precision has no practical indication at all. It is more of an aesthetic thing :) That said, the 48series (48G(X), 49g+, 50g) do return the uppertail probabilities to full displayed precisionUTPN, UTPC, UTPT, and UTPF.
Les
Re: [wp34s] Incomplete Gamma on the wp34s  Les Wright  12052011 Sure, Namir. I just ported the routines from an older version of Numerical Recipes. The authors of those books have put older versions online for free perusal, and the relevant section can be found here starting at page 216: Numerical Recipes in C, Second Edition In my code, the first few steps compute the series, the next several (the bulky middle) compute the continued fraction, and the last few steps where both routines end up compute and multiply in factors common to both.
Les Edited: 5 Dec 2011, 12:56 p.m.
Re: Maybe Not Completely Redundant  Paul Dale  12052011 Rather than guessing you could look at the source code ;)
Re: [wp34s] Incomplete Gamma on the wp34s  Paul Dale  12052011 The Numerical Recipes books don't seem to be very well respected in numerical mathematics circles. The code I've used from them had issues which needed to be addressed. They are easy for a layman to read however. Far more so than any of the other numeric texts I've got.
Re: Maybe Not Completely Redundant  Les Wright  12052011 Still working my way around the files. I must admit I can comprehend much of the code, but the special number types and functions acting on them make for some slow going.
Les
Re: [wp34s] Incomplete Gamma on the wp34s  Ángel Martin  12062011 I guess you've already seen JeanMarc Baillard's page: http://hp41programs.yolasite.com/gamma.php section 4 within that document. He uses two approaches, one with continued fractions and another with the Kummer function (a particular form of the Generalized Hypergeometric Function, HGF+) as an accessory to calculate the Incomplete Gamma  and many others BTW. I found that a very clever and efficient scheme, once the HGF+ function is tackeld down in MCODE that provides a slew of many others for the same price of admission. FWIW, these functions are included in JeanMarc's own JMB_MATH and in the SandMath modules.
/AM
Re: Maybe Not Completely Redundant  Paul Dale  12062011 These functions are all implemented in decn.c and stats.c. The code there mostly uses the decNumber type which is a real. The function calls should be self explanatory > fn(a, b, c) means a = fn(b, c). I've reused variables a lot for stack space reasons which does confuse things a bit.
Re: [wp34s] Incomplete Gamma on the wp34s  Namir  12062011 The HP Museum HP41C Software library pages also has an upper incomplete gamma function program ported by Tony Duell from the HP67 High Level Math Solution Book. The code is standalone and is still relatively short. I am a bit puzzled, because Duell only worked with the code for the upper incomplete gamma (whose definition is an integral between one of the function's arguments x and infinity) which seems, to me, to be harder to code than the lower incomplete function (whose definition is an integral from 0 to x) which should be easier to code. Namir
Edited: 6 Dec 2011, 10:25 a.m.
Re: [wp34s] Incomplete Gamma on the wp34s  Les Wright  12062011 Namir, I am familiar with that code you are noting, and I believe it computes the continued fraction in a very old fashioned way that is costly for steps and prone to overflow. I have always liked the NR routine for its brevity and its use of a more modern algorithm. I am particularly proud of a SysRPL port I did of it for my HP49g+ Angel, I am very familiar with JMB's work, and particularly his approach to the problem with the hypergeometric and Kummer routines. Mine is just one more way to skin that cat :)
Les
Re: [wp34s] Incomplete Gamma on the wp34s  Namir  12062011 Thanks for letting know about the vulnerability of the HP67 program ported to the HP41C. Nobody likes bad surprises when computing .... especially if you are in a space capsule going to the Moon, or worse yet, to Mars! Does Baillard implement routines for both the lower AND upper incomplete gamma functions? PS: I can see a cartoon where two astronauts are traveling to Mars and the display calculating the incomplete gamma function gives a mission critical error message! One astronaut turns to the other one and says, "I told you fifty million times use the Baillard code ... but noooooooooo .. you had to use the code translated from the HP67 library. We are so dead now!!"
Edited: 6 Dec 2011, 11:15 a.m.
