[wp34s] Incomplete Gamma on the wp34s

[Les Wright]

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

[Paul Dale]

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

[Les Wright]

Dang, I missed that one completely...

The stuff you guys have packed in... :D

[Les Wright]

I should point out that the documentation of the I-Gamma function is actually incorrect. It returns I-Gamma(y,x). The arguments are are reversed.

Les

[Les Wright]

Quote:

---

I should point out that the documentation of the I-Gamma function is actually incorrect. It returns I-Gamma(y,x). The arguments are are reversed.

---

Likewise with the incomplete beta function. It actually returns I-Beta(x,z,y).

Les

[Paul Dale]

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

[Les Wright]

Actually, the wp34s provides only the lower regularized incomplete gamma and cdf's (lower -tail) for phi and chi-square.

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 upper-tail probability associated with the normal statistic x = 3 I would compute phi(-3) instead.

But there is no such reflection relation for chi-square. A high-value 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 upper-tail probability. It is probably better to compute the upper-tail chi-square probability by way of the underlying continued fraction expansion for the right-sided incomplete gamma--IGL 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 built-in 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 built-in regularized incomplete gamma.

[Walter B]

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.

[Les Wright]

True dat, as the young people like to say today.

Obviously, my interest in computing the upper-tail functions to high precision has no practical indication at all. It is more of an aesthetic thing :)

That said, the 48-series (48G(X), 49g+, 50g) do return the upper-tail probabilities to full displayed precision--UTPN, UTPC, UTPT, and UTPF.

Les

[Paul Dale]

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

[Les Wright]

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

[Paul Dale]

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

[Namir]

Les,

Can you please provide a link to the algorithms you used to calculate the two versions of the incomplete gamma functions?

Namir

[Les Wright]

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.

[Paul Dale]

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

[Ángel Martin]

I guess you've already seen Jean-Marc 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 Jean-Marc's own JMB_MATH and in the SandMath modules.

/AM

[Namir]

The HP Museum HP-41C Software library pages also has an upper incomplete gamma function program ported by Tony Duell from the HP-67 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.

[Les Wright]

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

[Namir]

Thanks for letting know about the vulnerability of the HP-67 program ported to the HP-41C. 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 HP-67 library. We are so dead now!!"

Edited: 6 Dec 2011, 11:15 a.m.