(35S) Statistical Distributions Functions

[Dieter]

(10-30-2015 07:43 PM)PedroLeiva Wrote:  My apologies, I expressed badly. The system don´t let me send the * .EP files of the emulator programs that I type while testing what you were designing. Try your self; what the system warned me is that do not accept that file format for sending.

Yes, the forum software allows just a limited number of file types. But AFAIK you can simply zip the .ep file and send the resulting .zip-file.

Dieter

[PedroLeiva]

[/quote]
Excellent

[Dieter]

(10-29-2015 07:18 AM)Dieter Wrote:  I am not quite happy with the initial guess of the Chi² quantile. This is a very stripped down version of what's implemented in the 34s. But a more sophisticated version that usually converges within 3...5 iterations would require an estimate for the Normal quantile, which is also true for the Student inverse. So this can be done in a unified package that includes all three (or four) distributions.

Status report: I am currently working on such a package, and I think it looks not too bad. The package can evaluate the Normal, Student, Chi² and partially F distribution's PDF, CDF and their inverses (except F). The inverse usually requires just a few iterations since the initial guesses are close to the true result. Switching between CDF/PDF and inverse is done by setting Flag 1 (think of "1nverse" ;-)). There are common input routines for the degrees of freedom (with error checks), the random variable X and the probability p.

Numerically the major limits are the far Student's distribution tails, so if you need accurate CDFs below, say, 1E–10 a more sophisticated algorithm has to be used. That's why I am trying to implement the regularized Beta function. Which currently works fine with a series expansion for most practical cases of the F distribution.

Here are some results:

Code:

Normal distribution
    [XEQ] N [ENTER]     M?
Enter my
     0 [R/S]            S?
Enter sigma
     1 [R/S]            X?
Enter random variable
     1,5 [R/S]          0,129517596   PDF
                        0,933192799   CDF
Do another calculation
     [R/S]              X?
     -2 [R/S]           0,053990967   PDF
                        0,022750132   CDF
Calculate a quantile
     SF 1 [R/S]         P?
     0,95 [R/S]         1,637300567
                        1,644853622

     [R/S]              1,644853622
                        1,644853627

Student's t distribution
    [XEQ] T [ENTER]     J?
Enter degress of freedom
     10  [R/S]          P?
Determine 90% quantile
     0,9 [R/S]          1,349994793
                        1,372181030
     [R/S]
                        1,372181030
                        1,372183641
     [R/S]
                        1,372183641
                        1,372183641
Switch back to CDF mode
     CF 1
and start over
     [XEQ] T [ENTER]    J?
     [R/S]              X?
     2,5  [R/S]         0,026938728   PDF
                        0,984276578   CDF

Chi² distribution
    [XEQ] C [ENTER]     J?
Enter degrees of freedom
     5,7 [R/S]          J?
rejected, as dof must be a positive integer
     5 [R/S]            X?
Enter random variable
     9,5 [R/S]          0,033688016   PDF
                        0,909292608   CDF
Determine 90% quantile
     SF 1 [R/S]         P?
     0,9  [R/S]         9,166552063
                        9,236250476
     [R/S]
                        9,236250476
                        9,236356899
     [R/S]
                        9,236356899
                        9,236356900
Reset to CDF mode
     CF 1

Fisher's F distribution
    [XEQ] F [ENTER]     J?
nominator degrees of freedom
     20 [R/S]           K?
denominator degrees of freedom
     10 [R/S]           X?
Enter random variable
     2,8 [R/S]          0,048540378   1-CDF
                        0,951459622   CDF

(Fisher PDF and quantile are not yet implemented)

The F distribution is evaluated via an (experimental) regularized Beta function. Once this runs sufficiently fast and accurate, this can also be used for the Student case.

Dieter

[rprosperi]

(10-30-2015 07:43 PM)PedroLeiva Wrote:  My apologies, I expressed badly. The system don´t let me send the * .EP files of the emulator programs that I type while testing what you were designing. Try your self; what the system warned me is that do not accept that file format for sending.

The forum will only allow upload of file types that it recognizes. The simplest solution for you is to add the files into a .zip file, then upload them. They are also smaller this way, which helps with your quota.

[PedroLeiva]

Dieter, to contribute to your work on Statistical Distribution Functions, herewith a pdf with instructions for execution of the compiled program you are developing. I ask for corrections and comments. From already thank you very much for your effort.
NOTE: it was impossible to attach an Excel file

[Dieter]

(11-01-2015 12:44 PM)PedroLeiva Wrote:  Dieter, to contribute to your work on Statistical Distribution Functions, herewith a pdf with instructions for execution of the compiled program you are developing. I ask for corrections and comments. From already thank you very much for your effort.

Thank you – I appreciate your efforts. You did quite some work for a program that is still under development. ;-)

The instructions say that "x" is calculated for a quantile. Actually the quantile is x. The programs determine the quantile x for a given probability p. The latter is not given in % or ‰ (c.f. Student distr.). After each iteration step, Y and X show the previous resp. current approximation. This is helpful in some extreme cases where both the CDF and the quantile currently can only be calculated with limited accuracy.

The returned CDF is always the integral from –∞ to x (resp. 0 to x for Chi²), i.e. the "left side integral", which is usually called P(x). This is also true for the F distribution where HP's manual says that P(x) is the "right side integral", which usually is called Q(x). So the current F implementation returns Q(x) in Y and P(x) in X.

I know want to change the general approach in that I'd like to provide two basic routines that evaluate the reqularized incomplete Beta and Gamma function (actually four routines as both functions are evaluated either by a series expansion or continued fractions). Beta can be used to calculate the Student and Fisher distribution, and Gamma does it for the Normal and Chi² distribution. Execution times seem to be a bit longer, but on the other hand now even cases far out in the distribution tails can be handled. The plan is to get reasonably accurate results even for p=1E–499. ;-)

(11-01-2015 12:44 PM)PedroLeiva Wrote:  NOTE: it was impossible to attach an Excel file

That's easy: simply zip it.

Dieter

[Dieter]

(11-01-2015 09:17 PM)Dieter Wrote:  Thank you – I appreciate your efforts. You did quite some work for a program that is still under development. ;-)

Update: I am still working on this program package. The current versions return both the PDF and the CDF simultaneously. Now I wonder if the PDF really is that useful. Maybe it's better to return both P(x) and Q(x) = 1–P(x), especially for values far out in the distribution tails where P(x) rounds to 1. What do you think?

Dieter

[PedroLeiva]

Update: I am still working on this program package. The current versions return both the PDF and the CDF simultaneously. Now I wonder if the PDF really is that useful. Maybe it's better to return both P(x) and Q(x) = 1–P(x), especially for values far out in the distribution tails where P(x) rounds to 1. What do you think?

Dieter
[/quote]
I think you are right, at list only P(x) and Q(x) for Fisher distribution, for the rest, its ok to keep PDF calculation. Always is importan to have "x" value from probability (inverse)

[Dieter]

(11-05-2015 02:16 PM)PedroLeiva Wrote:  I think you are right, at list only P(x) and Q(x) for Fisher distribution, for the rest, its ok to keep PDF calculation. Always is importan to have "x" value from probability (inverse)

I am still working on it.
Here is a list of what the current version can do:

Normal distribution:

• PDF
• lower and upper CDFs P(x) resp. Q(x)
• two-sided symmetric CDF A(x) and its complement B(x)
• quantile function x(p)
  

Student's t distribution:

• PDF
• lower and upper CDFs P(t) resp. Q(t)
• two-sided symmetric CDF A(t) and its complement B(t)
• quantile function t(p)
  

Chi² distribution:

• PDF
• lower and upper CDFs P(c) resp. Q(c)
• quantile function c(p)
  

Fisher's F distribution:

• PDF
• lower and upper CDFs P(F) resp. Q(F)
• quantile function not yet implemented
  

In addition, some internal special functions can be used as well:

• Euler's Beta function B(a,b)
• Regularized Incomplete Beta functions Ix(a,b) and 1–Ix(a,b)
• Regularized Incomplete Gamma functions P(a,x) and Q(a,x)
  

The distribution results are stored in registers A...D, so the user can easily adjust the values that are returned in Y and X.

I just don't know when and if this will ever get finished. ;-)

Dieter

[PedroLeiva]

[I am still working on it.
Here is a list of what the current version can do:

Normal distribution:

• PDF
• lower and upper CDFs P(x) resp. Q(x)
• two-sided symmetric CDF A(x) and its complement B(x)
• quantile function x(p)
  

Student's t distribution:

• PDF
• lower and upper CDFs P(t) resp. Q(t)
• two-sided symmetric CDF A(t) and its complement B(t)
• quantile function t(p)
  

Chi² distribution:

• PDF
• lower and upper CDFs P(c) resp. Q(c)
• quantile function c(p)
  

Fisher's F distribution:

• PDF
• lower and upper CDFs P(F) resp. Q(F)
• quantile function not yet implemented
  

In addition, some internal special functions can be used as well:

• Euler's Beta function B(a,b)
• Regularized Incomplete Beta functions Ix(a,b) and 1–Ix(a,b)
• Regularized Incomplete Gamma functions P(a,x) and Q(a,x)
  

The distribution results are stored in registers A...D, so the user can easily adjust the values that are returned in Y and X.

I just don't know when and if this will ever get finished. ;-)

Dieter

Dear Dieter
Could you try to publish the F Fisher program individually, while the compile versión is finished?
With my sincere recognition of your hard work, Pedro

[Dieter]

(11-13-2015 08:04 PM)PedroLeiva Wrote:  Could you try to publish the F Fisher program individually, while the compile versión is finished?

The current program pack uses a number of common routines for the various distributions. For instance, the Fisher module requires the regularized Beta function, Euler's Beta and two common input routines. This approach also removes the contraints of the HP67/97 stat pac's Fisher routine where at least one of the df has to be even.

All four distributions are evaluated using either the regularized Gamma or Beta function. These two functions even allow an easy implementation of the Binomial and Poisson distribution.

In other words: extracting only the Fisher code does not make much sense since several other routines are required as well.

I recently posted a short 35s program for Student's t distribution. It runs quite fast and works well for most practical cases, i.e. for probabilities not too close to 0 or 1. That's what I like about it. The current version is based on the regularized Beta function, which does not run that fast but it allows some nice things like evaluating the quantile for, say, 10 df and P(x)=1E–400 (result = –2,56452571895 E+40).

The complete code now has about 2600 Bytes (without the old Student version). Especially the quantile functions include some know-how from the respective WP34s routines. IIRC the whole pack consists of eleven routines in total. It's getting bigger than I initially thought. ;-)

Dieter

[Dieter]

(11-13-2015 06:59 PM)I Wrote:  I am still working on it.
Here is a list of what the current version can do:

(...)

Fisher's F distribution:

• PDF
• lower and upper CDFs P(F) resp. Q(F)
• quantile function not yet implemented
  

Update: I know implemented a (not too sophisticated) function for the Fisher quantile. It works quite well for both df > 2, but there are some issues for df=1 or 2, at least if at the same time the probability is close to 0 or 1. Otherwise it seems OK. At least the results are good enough to show the one or other inaccuracy in the Fisher quantile tables at nist.gov. And they are better than just three decimals. ;-)

Let me know if you want to play around with an emulator .ep file.

Dieter

[PedroLeiva]

[/quote]
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. From already thank you very much. Pedro

[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.

[PedroLeiva]

Dear Dieter: I have used your “Distribution compile functions file”. It works perfect for the four statistical functions. Very nice and detailed instructions & examples for Student “t” distribution. In spite the four distributions are handled much the same, I wish to have instructions for the other three. Maybe I am asking for too much…….I will cooperate trying this task myself. Pedro

[Dieter]

(11-16-2015 01:13 PM)PedroLeiva Wrote:  Dear Dieter: I have used your “Distribution compile functions file”. It works perfect for the four statistical functions.

I just found an error in the Normal quantile code. You do not have to replace the whole .ep file. Just go to step N078 XEQN024 and insert RCL S  STOx D behind it.
The attachment has been replaced by a corrected version (with an additional RCL/ S, which is one step shorter).

(11-16-2015 01:13 PM)PedroLeiva Wrote:  Very nice and detailed instructions & examples for Student “t” distribution. In spite the four distributions are handled much the same, I wish to have instructions for the other three.

I'll see what I can do.

Dieter

[Dieter]

(11-16-2015 07:08 PM)Dieter Wrote:  I'll see what I can do.

OK, here is a zip file with the PDFs for the four distributions.
Maybe I'll add another version for the additional functions later (Gamma, Beta etc.).

Dieter

  Stat_distributions_instructions.zip (Size: 301.7 KB / Downloads: 16)

[PedroLeiva]

(11-16-2015 10:47 PM)Dieter Wrote:  

(11-16-2015 07:08 PM)Dieter Wrote:  I'll see what I can do.

OK, here is a zip file with the PDFs for the four distributions.
Maybe I'll add another version for the additional functions later (Gamma, Beta etc.).

Dieter

Thank you very much Dieter. It is a great job. Pedro

[Dieter]

(11-17-2015 12:44 PM)PedroLeiva Wrote:  Thank you very much Dieter. It is a great job. Pedro

I just found another error in the incomplete Gamma function. It only shows up when the function is called independently from the distributions. Line G010 must be GTO G013 (not G014).

Here is a corrected version along with a complete and reworked set of instructions in one single PDF file. Please check this and report any errors you find.

Edit: the second example for the F-distribution has an error. The results for F(5) with 10 and 15 d.o.f. of course are 0,00275 and 0,99725. Which the program correctly returns. Maybe I'll post a revised PDF together with a slightly tweaked program version that can do one or two more tricks.

Dieter

  Distributions_35s_3.zip (Size: 139.25 KB / Downloads: 20)

[Gerald H]

Could you, Dieter, please include a pdf file of the latest programmes in human readable form to enter on the calculator?