Posts: 1,368
Threads: 212
Joined: Dec 2006
Thought there would be some interest in this article:
Incomplete Gamma for wp34s
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
Posts: 3,229
Threads: 42
Joined: Jul 2006
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.
 Pauli
Posts: 1,368
Threads: 212
Joined: Dec 2006
Dang, I missed that one completely...
The stuff you guys have packed in... :D
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
Posts: 1,368
Threads: 212
Joined: Dec 2006
Quote:
I should point out that the documentation of the IGamma function is actually incorrect. It returns IGamma(y,x). The arguments are are reversed.
Likewise with the incomplete beta function. It actually returns IBeta(x,z,y).
Les
Posts: 3,229
Threads: 42
Joined: Jul 2006
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.
We've also got a regularised incomplete beta function, again for real arguments only. And again, in support of the statistical functions.
 Pauli
Posts: 1,368
Threads: 212
Joined: Dec 2006
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.
Posts: 4,587
Threads: 105
Joined: Jul 2005
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.
Posts: 2,247
Threads: 200
Joined: Jun 2005
Les,
Can you please provide a link to the algorithms you used to calculate the two versions of the incomplete gamma functions?
Namir
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
Posts: 1,368
Threads: 212
Joined: Dec 2006
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.
Posts: 3,229
Threads: 42
Joined: Jul 2006
Rather than guessing you could look at the source code ;)
Yes, the upper tail probabilities will suffer from cancellation  we ran out of space to implement both properly.
 Pauli
Posts: 3,229
Threads: 42
Joined: Jul 2006
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.
 Pauli
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
Posts: 1,253
Threads: 117
Joined: Nov 2005
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
Posts: 3,229
Threads: 42
Joined: Jul 2006
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.
 Pauli
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
Posts: 2,247
Threads: 200
Joined: Jun 2005
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.
