In double precision mode, the following ties up the calculator for a painfully long period of time, and ultimately displays "Domain Error":

t^(-1)(.5001, 5), i.e. the t-quantile for the probability given with 5 df

However, the following equivalent problem readily returns the correct result to 34 digits:

sqrt(F^(-1)(.0002, 1,5)) ~= 2.6343e-4 

The calculator will work a long time yet return the result in single precision, but if the argument is much closer to 0.5, say 0.500000001, the Domain Error gets thrown in single precision too when one attempts to compute the t-quantile directly.

This occurs with the most recent version of the firmware (2810).

Oddly, it does NOT occur with Qt emulator on my Mac, which reports a firmware version of 2783.

Could it be that some of those changes I instigated regarding the inverse F CDF ended up breaking the inverse t CDF for probabilities near 0.5 (where the quantiles would be close to zero)?

Hope this was worth reporting.

Les

Edited: 17 Apr 2012, 7:50 p.m.

Yes this will have been caused by the recent changes :-(

I'm pretty sure I know what is happening but will have to investigate properly later.

- Pauli

The next build addresses this partially.

For single precision, it should be good.
For double precision, it can still fail with a domain error.

The issue was to do with the Newton steps wandering around very close to the solution but never reaching a point considered good enough, the cdf is returning the same value at multiple points in this neighbourhood. I suspect the real underlying issue is some loss of precision either in the T cdf or the subtraction from (almost) one half of another (almost) one half value.

A further down side is I had to increase the iteration limit to allow the bisection steps, which are now taken if this condition is detected, to converge sufficiently. In other words, maximum run time is now longer than before.

I've no real idea how to fix this properly for double precision mode :-(

- Pauli

I think the solution may be to compute the quantile as a special case of the inverse F distribution, as I did in my opening post:

For p = 0.5, t^(-1)(p, df) = 0

For p > 0.5, t^(-1)(p, df) = +sqrt(F^(-1)(2p-1, 1, df))

For p < 0.5, t^(-1)(p, df) = -sqrt(F^(-1)(1-2p, 1, df))

This turns an inverse-t-close-to-0.5 problem into an inverse-F-close-to-zero problem, where the inverse F actually fares rather well.

I wonder what Dieter thinks? If computing the inverse t CDF as a special case of the inverse F CDF saves some complexity and bytes while preserving most of the accuracy, ditching the dedicated inverse t code altogether could be an option.

Les

I do not think it's a good idea to evaluate 1-2p for p<0.5 (or 2p-1 for p>0.5), since in both cases may round to 1 for small p. Also, the inverse F estimate requires both df to be larger than 1 as it computes 1/(n1-1) and 1/(n2-1).

I think the Student quantile function should return correct results within reasonable time. To me it sounds more like a problem with the Newton solver that will not find a final 34-digit result. If the problem really is in the Student CDF for t close to 0, this should not be too difficult to fix.

Dieter

Fair enough. But for arguments close to 0.5, shunting things through the inverse F CDF with numerator df of 1 and a suitably transformed probability is a good way to go.

For example, in double precision the smallest number just above .5 that the calc can represent is (0.5 + 1e-34}. This corresponds to computing the inverse F CDF for p = 2e-34, which is returned readily and to 34 digits within 1 or 2 ULP.

Here is what I would propose (in very rough pseudocode):

If abs(p-0.5) < something smallish (e.g, .01)

Then

Compute the inverse F CDF for 2p-1 or 1-2p as appropriate and return the positive or negative sqrt as appropriate

Else

Do it directly with the existing inverse t CDF code and return that

End

Indeed, I may write myself a FOCAL wrapper to the inverse t that does just this very thing. As we know, computing the Student t quantile to 34 digits precision for probabilities close to 1/2 is a highly relevant and practical problem in the daily life of statisticians everywhere. ;)

Les

There are plans afoot to eventually rework the statistical code into FOCAL programs as a space saving measure. I'd welcome code donations :-)

- Pauli

Where would they live in such a case? XROM? Or would they be demoted to "optional" routines in the standard wp34s-lib.dat library? Would the "PROB" menu be no more?

Just wondering since I usually flash only with calc.bin and use my own wp34s-lib.dat, and not the one provided.

Les

XROM definitely. We've no plans on removing any of the internal functions.

- Pauli

Thanks, Pauli. I really wish I had a better grasp of the layout of the calculator's "brain"--where stuff goes, etc.

The layout isn't that difficult.

From the user's perspecitve:

      RAM -- 2kb for everything.  Most of this can be registers and programs.

      Flash -- circa 8kb of library.

Internally, there are two major divisions:

      XROM -- keystroke programs with a few special abilities.

      the rest -- raw machine code stuff.  Some of this is implementing functions, some display and some the keyboard handler.

We're migrating mathematical functions from "the rest" to XROM. This saves circa 50% of the space used for these. We're also running out of things to move -- very low level functions can't easily move over and most of the higher level ones have been moved already. The statistical distributions are one big area that hasn't really migrated yet. The other is portions of the solver.

When I get time and motivation, I'll address these but it is a huge task for either.

- Pauli