[WP34S] Improving the Upper-Tail F and Chi-Square


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


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.


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


Thanks for that. Too much beer tonight to fix this. I leave it to Pauli. :-)


Fixed the F problem. I looked at that equation last night and completely missed the obvious rewrite :-(

I agree we need the gammaQ function to complement the gammaP I implemented. Time as usual...

- Pauli


I've implemented the other half of the regularised incomplete gamma function and used them in the cumulative functions.

I get something near your df=5, x=1000 Chi2 result.

Firmware and emulators building....

- Pauli


Thank you!



With the implementation of both sides of the incomplete gamma function, I've freed up some space.

Any suggestions for more statistical distributions to include?? Hypergeometric? Negative Binomial? Pareto? Bernoulli? Beta? Gamma? Erlang?

Any other functions missing from the 34S function suite? Bessel functions won't fit sorry.

Not that we're planning on changing anything, but some idea of the shortcomings (hah!) of the device would be interesting.

- Pauli


Do I sense a tongue firmly lodged in cheek? ;)

Seriously, having nearly full DP results out of the top four distributions (normal, t, F, chi-square) and their inverses is impressive. The others--Poisson, Cauchy, Weibull, etc.--are gravy. The ones you mention are special application esoterica. Anyone with abundant time and overweening nerdiness can regard the boggling contributions of J-M Baillard to the HP41 library on this site and port things for possible inclusion in the user library.

Baillard's programs (many of which have been ported to MCODE and special modules by Angel Martin, I believe) include Bessel functions, the simplex linear programming algorithm (a personal favourite of mine), and, very impressively, the Bullrisch-Stoer differential equation solver.

Why would we want to use the WP34S to numerically solve an ODE in a few minutes that takes Mathematica or Maple mere seconds? Why, just because. ;)


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


Do I sense a tongue firmly lodged in cheek? ;)

Only half in cheek :-)

I am tempted to add the uniform and discrete uniform distributions, not so much because they are statistically very useful but more because they provide useful functions as a byproduct.

Want to linearly interpolate between two points? Use the uniform QF. Want to know how far between two points a third is? Use the uniform CDF. Want to Roll a dice? RAN# followed by a discrete uniform QF call. All handy functions to have built in.

Baillard's programs (many of which have been ported to MCODE and special modules by Angel Martin, I believe) include Bessel functions, the simplex linear programming algorithm (a personal favourite of mine), and, very impressively, the Bullrisch-Stoer differential equation solver.

I've used Jean-Marc Baillard's HP-41 program site often. A few of the 34S built in functions were based on his implementations.

- Pauli

Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP-Prime CAS] "Warning, ^ (Command) Is ambiguous on non square matrices"?? CompSystems 1 1,213 12-07-2013, 07:15 PM
Last Post: CompSystems
  HP PRIME: PISERIES without h tail on binaries Joseph Ec 0 493 11-21-2013, 11:14 AM
Last Post: Joseph Ec
  HP Prime Error: Summation Upper Bound > 1000 HP Pioneer 2 858 10-25-2013, 11:32 AM
Last Post: steindid
  A brand new calculator benchmark: "middle square method seed test" Pier Aiello 25 4,320 09-13-2013, 01:58 PM
Last Post: Pier Aiello
  Square Root Simplifier for HP39gII Mic 4 1,305 03-11-2013, 08:25 AM
Last Post: Eddie W. Shore
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 742 10-04-2012, 04:59 PM
Last Post: Paul Dale
  HP41 Upper Screw Well Busted Colin Verrilli 6 1,331 07-16-2012, 08:45 PM
Last Post: matt kernal
  [WP34S] Improving the Incomplete Gammas Les Wright 19 3,032 05-16-2012, 05:15 PM
Last Post: Marcus von Cube, Germany
  [WP34S] Chi-Square Functions Broken in 2935 Les Wright 7 1,503 05-03-2012, 11:10 AM
Last Post: Les Wright
  [wp34s] Incomplete Gamma on the wp34s Les Wright 18 3,096 12-06-2011, 11:07 AM
Last Post: Namir

Forum Jump: