Yup, it's me again obsessing about the statistical distributions.

Dreaded digit loss due to subtraction rears its head again.

In the case of the upper F, the the half-df parameters are swapped correctly before passage to IBETA. However, the main argument is computed in a way that begs for digit loss.

The code computes this as:

1 - d1*x/(d1*x+d2)

This is mathematically correct but computationally challenged. In situations where d1*x is huge compared to d2, the second term of this goes to 1 and the whole thing goes to zero. Indeed, it may very well BE zero within the calculator's 34-digits. Zero gets passed to IBETA, which returns zero, even if the true F score is a tiny non-zero result representable by the calculator.

This is easily solved by elementary algebra. Arrange said term to

d2/(d1*x+d2)

and pass that. Indeed, this is precisely what I do in my own HP41 routine that computes BETA, IBETA, and the upper and lower tail t and F distributions. For example, for the absurdly unrealistic F-score of 7e19 at 5 numerator and 10 denominator df, my HP41 (EMU41, actually) returns 2.233444399e-97, which is just 1ULP off of the true value. WP34S's F-sub-u returns zero, even though this tiny upper tail probability is easily represented by the calculator in both single and double precision. Indeed, when I compute the parameter "my" way and pass the correct arguments manually to the calculator's IBETA function, I get a 34-digit result that is off only 4ULP.

That is an easy enough fix, I think. The upper-tail chi-square is more challenged. In that case, it only has access to the left-sided IGAMMA. So, if the statistic is much larger than the df, when the adjusted parameters are passed to IGAMMA, it is possible to get back unity. The upper-tail complement of this is zero as far as the calculator is concerned, even though the true value may be a tiny non-zero result that the calculator can easily handle.

The challenge here is that the calculator does not offer a RIGHT-sided incomplete gamma exposed to the XROM (or user) for use. There is only IGAMMA, which calls gammap, which (correctly) uses the series (gser) when the x parameter is less than the a parameter and takes the complement of the continued fraction result (from gcf) in the opposite case. But there is no right-sided incomplete gamma analogue, say gammaq. Whether or not it is user-accessible (and I vote it should be), gammaq should be easy enough to code and expose to XROM since it does the OPPOSITE of gammap--i.e., for a < x use the gcf and return that, if a > x, use gser and return the complement 1 - gser. To compute Chi-square-sub-u, one simply calls gammaq as one does IGAMMA, and thus gammap, for the cumulative CDF.

As an alternative, one needn't write a new gammaq routine, but rather permit XROM access to the inner gcf routine. For the upper CDF, if the statistic is less than the degrees of freedom, one simply calls IGAMMA (gammap) as before and takes the complement. BUT, if the statistic is greater than the df, in the absence of a gammaq wrapper one needs to call gcf directly. Indeed, this is what happens when one uses IGAMMA when the statistic is much greater than the df: the C code gammap routine is called, which calls gcf, subtracts that small result from unity, returns a result close unity (or even unity itself) to the stack, and again that result is subtracted from unity. Subtraction loss galore due to two subtractions that basically cancel each other out and shouldn't have to happen in the first place.

For example, with my user code which offers a right-sided incomplete gamma, I get an upper-tail chi-square probability for 1000 at 5 df that is approximately 6.01e-214, accurate within 2 ULP in the 34-digit result. The built-in Chi-square-sub-u returns zero, even through the true result and all intermediate calculations are well within range.

I know there is statistical distribution burnout out there, so this is back burner stuff, and I just offer this as food for thought for future reference.

Les

*Edited: 12 May 2012, 10:35 p.m. after one or more responses were posted*