(35S) Statistical Distributions Functions

[Dieter]

(11-15-2015 03:04 PM)PedroLeiva Wrote:  Much I would like to try the Fisher F function with the emulator (* .EP), please also send me the instructions for the incoming and outgoing data.

I have attached a zip file that contains both the emulator .ep file and a PDF with instructions for Student's t distribution. The four distributions (XEQ N, S, C and F) are handled much the same, so instructions for one of them should do. Only the two asymmetric distributions (Chi² and F) of course do not return the two-sided symmetric integrals A(x) resp. B(x) = 1–A(x).

EDIT: the original attachment has been removed, please see note at the end of this post.

The programs run either in cdf mode (flag 1 clear) or in quantile mode (flag 1 set). In cdf mode, various cumulative distribution functions and the density function are calculated and returned on the stack and in registers A...D. In quantile mode, the value of the random variable is determined for which the (lower tail) cdf equals the entered probability ("inverse cdf").

The quantile functions use Halley's method for the Chi² and Student case and a very effective extrapolation scheme due to Abramovitz & Stegun for the Normal distribution. The Fisher quantile is handled by a simple Newton iteration, it is intended for both degress of freedom larger than 2. Here and there the iteration may diverge, some code has been added to handle this, so you may also try 1 or 2 degrees of freedom, but this may require some time until the iteration converges.

Please note: all this is experimental code! The programs may, and probably will, contain any arror you can imagine, and even those you can't. ;-) So try this and see what you get. In the meantime I am waiting for your error reports.

Some remarks:

• LBL E evaluates a rational approximation for a rough estimate of the Standard Normal quantile. It expects a probability 0<q≤0.5 in register Q and returns a positive result. It handles q values down to Q(50)≈1E–545. The error of this approximation is taylored to the mentioned extrapolation scheme, the absolute error is about 2,5E–4 in the center and less than 5E–5 for results > 20. This way one single extrapolation step should yield a result with approx. 10 valid significant digits.
• LBL A is an input routine for x resp. p. For the latter there is a range check and min(p, 1–p) is stored in Q, along with the appropriate sign for the Student and Normal distribution, which is stored in E (+1 if p>0.5 resp. –1 if p<0.5).
• LBL J is an input routine for positive integers, i.e. the degrees of freedom j and k as required by the Student, Fisher and Chi² distribution.
• LBL C uses two different ways of calculating the Chi² cdf. For large x the incomplete Gamma function is evaluated, while for smaller x (below df) a short and fast algorithm similar to the one in the HP67/97 stat pack used.
• LBL I gives access to the Regularized Incomplete Beta function Ix(a,b) which is used for the Student and Fisher distribution. It implements a continued fraction method with a modified Lentz algorithm. On exit, X and Y hold Ix resp. 1–Ix.
• LBL G gives access to the Regularized Incomplete Gamma functions P(a,x) and Q(a,x) which are used for the Normal and Fisher distribution. It evaluates either a series expansion or a continued fraction (modified Lentz algorithm). On exit, X and Y contain P(a,x) resp. Q(a,x).
• LBL B allows the calculation of Euler's Beta function B(a,b).
  

Finally, here is a test case for the Fisher distribution:

For 10 and 15 degrees of freedom, determine the probability for an F-value of 2,7.

Code:

Just to be sure, clear flag 1:
[FLAGS] CF 1

XEQ F [ENTER]   J?
   10 [R/S]     K?
   15 [R/S]     X?
  2,7 [R/S]
                0,040332463   Q(x)
                0,959667537   P(x)

So the probability that F is less than 2,7 is about 96%.

For a two-sided F-test, find the critial F-values for a significance level of 5%.
This means we are looking for the two F-values that yield a probability of 0,05/2 resp. 1–0,05/2.

Code:

Set quantile mode:
[FLAGS] SF 1
      [R/S]     P?
0,025 [R/S]
                0,282416352
                0,283964407
      [R/S]
                0,283964407   // so F~0,2839...
                0,283955930   // usually you can stop at this point
      [R/S]
                0,283955930
                0,283955930

Start over:
   0  [R/S]     P?
0,975 [R/S]
                3,063246171
                3,060190176
      [R/S]
                3,060190176   // so F~3,0602
                3,060196851   // which usually is sufficiently accurate
      [R/S]
                3,060196851
                3,060196851

Now try and see what you get. I am waiting for your error reports and other suggestions. ;-)

Dieter

EDIT and WARNING: this post originally contained an attachment with an emulator file (.ep) and instructions for Student's t distribution. However, there were a few errors in the program that now have been corrected. Please see post #59 below for the latest version along with a complete set of instructions. Please do not use older program versions before Nov 17, 2015.