(35S) Statistical Distributions Functions

(35S) Statistical Distributions Functions - Printable Version

+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (35S) Statistical Distributions Functions (/thread-5005.html)

Pages: 1 2 3 4

(35S) Statistical Distributions Functions - PedroLeiva - 10-25-2015 12:19 AM

I need programs of Statistical Distributions: t Student, F Snedecor and Chi Square. The Hp 67 Statistical Pakage was excellent, but I am not able to translate them to HP 35s. Does somebody could help me in this task?

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-25-2015 01:29 PM

The HP 35s has two features that are missing in the HP-67:

integration
gamma function by using $\Gamma(n)=(n-1)!$

Thus I wouldn't recommend to translate these HP-67 programs directly to the HP 35s. I just had a quick look at Chi-Square Distribution for the HP-67 and it uses a series approximation to evaluate the cumulative distribution. With the HP 35s you can use numeric integration instead.

Just have a look at the program Normal and Inverse–Normal Distributions in the user's guide (pp 16-11). You can replace the subroutine F that calculates the integrand for the normal distribution. It calculates the "upper tail" area but that is easy to change. Since integration and solver can't be mixed a Newton-iteration is used to calculate the inverse. The calculation of the derivative is just evaluating the probability density function due to the fundamental theorem of calculus.

Let us know if you struggle to translate the formulas of the probability density functions to programs. Using an algebraic equation might be easier but will be probably slower.

HTH
Thomas

RE: HP 35s Statistical Distributions Functions - SlideRule - 10-26-2015 03:03 PM

(10-25-2015 12:19 AM)PedroLeiva Wrote: I need programs of Statistical Distributions: t Student, F Snedecor and Chi Square. The Hp 67 Statistical Pakage was excellent, but I am not able to translate them to HP 35s. Does somebody could help me in this task?

Might I suggest...
NCS35TB Stat Pac
For the HP-35s Calculator
... by Namir Shamas
at NCSSTATPAC

BEST!
SlideRule

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-26-2015 04:41 PM

Thank you very much for your recommendation, but Namir programs are not what I need. In Appendix A, there are detailed calculations based on estimates of the inverse probability of each 4 statistical functions. I also need to calculate probability according to a value of each function.

HP 67 Statistical Distributions programs (Normal, Chi Square, "F" and "t" Student) are perfect for this purpose.

RE: HP 35s Statistical Distributions Functions - Dieter - 10-26-2015 08:30 PM

(10-26-2015 04:41 PM)PedroLeiva Wrote: Thank you very much for your recommendation, but Namir programs are not what I need. In Appendix A, there are detailed calculations based on estimates of the inverse probability of each 4 statistical functions. I also need to calculate probability according to a value of each function.

OK – what exactly do you need? The CDFs for Normal, Student, Chi² and Fisher's F distribution? Also their quantile (inverse) functions? What accuracy is required?

BTW, the 34s does all this. I provided some code for these four quantile functions which are evaluated with mostly 30+ digit accuracy. This is possible due to CDFs that are (mostly) equally accurate. These CDFs in turn are calculated using incomplete Beta and Gamma functions, coded by Pauli.

(10-26-2015 04:41 PM)PedroLeiva Wrote: HP 67 Statistical Distributions programs (Normal, Chi Square, "F" and "t" Student) a perfect for this purpose.

Could you provide a pointer to this package? How accurate are the results? Do these programs also calculate the inverse (quantile) functions?
If you really appreciate these HP67 programs so much it should be easy to translate them for the 35s.

Dieter

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-26-2015 08:39 PM

Here's a program for the Student's t-distribution:

Code:

T001 LBL T       ; Student's t

T002 INPUT K     ; degrees of freedom

T003 1

T004 -

T005 2

T006 /

T007 !           ; Γ((k+1)/2) = (k-1)/2 !

T008 RCL K

T009 2

T010 /

T011 1

T012 -

T013 !           ; Γ(k/2) = (k/2 - 1)!

T014 /

T015 PI

T016 RCL* K

T017 SQRT

T018 /

T019 STO G

T020 RTN

I001 LBL I       ; inverse CDF ( CDF -- X )

I002 INPUT A     ; area

I003 0.5

I004 -

I005 RCL/ G

I006 STO X       ; initial guess

I007 XEQ C003    ; calulate cdf(X)

I008 RCL- A      ; we want this difference to become 0

I009 RCL X

I010 STO T

I011 R↓

I012 XEQ D001    ; derivative at X = pdf(X)

I013 /

I014 STO- X      ; improved guess

I015 ABS

I016 RCL E       ; small like 1e-6

I017 x<y?        ; close enough?

I018 GTO I007    ; not yet

I019 RCL X

I020 RTN

C001 LBL C       ; CDF

C002 INPUT X

C003 0           ; lower limit

C004 RCL X       ; upper limit

C005 FN= D       ; the distribution

C006 ∫FN d T     ; integrate pdf

C007 0.5         ; we add 0.5 since

C008 +           ; the lower limit was 0

C009 RTN

P001 LBL P       ; PDF

P002 INPUT X

P003 STO T

D001 LBL D       ; Distribution: ( T -- PDF )

D002 RCL T

D003 x^2

D004 RCL/ K

D005 1

D006 +           ; 1 + t^2/k

D007 LASTx

D008 RCL+ K

D009 2

D010 /           ; (1 + k)/2

D011 y^x

D012 RCL G

D013 x<>y

D014 /

D015 RTN

You have to initialize the program with the degrees of freedom:

XEQ T ENTER
K?
7
R/S

To calculate the PDF use the program P:

XEQ P ENTER
X?
1.6
R/S
0.110665

PDF[StudentTDistribution[7],1.6]

To calculate the CDF use the program C:

XEQ C ENTER
X?
1.5
R/S
0.911351

CDF[StudentTDistribution[7],1.5]

For the next program set the accuracy to a reasonable value:

1E-6
STO E

To calculate the Inverse CDF use the program I:

XEQ I ENTER
X?
0.99
R/S
2.997951

HTH
Thomas

RE: HP 35s Statistical Distributions Functions - Dieter - 10-26-2015 09:25 PM

(10-26-2015 08:39 PM)Thomas Klemm Wrote: Here's a program for the Student's t-distribution:

Hint: the Normal, Student, Chi² and Fisher quantiles will converge faster if the Halley method (instead of simple Newton) is used. This is easy here since the second derivative simply is the PDF times a factor. Another tweak that can provide better convergence (especially in the Chi² and F case) is solving ln(CDF/p)=0 instead of CDF–p=0.

Both methods are used in the WP34s quantile functions.

Dieter

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-26-2015 09:28 PM

For the Chi-squared distribution you just have to calculate a different PDF:

$[Image: a107033ede4e5de93c3b7817fc5aeb33.png]$

The factor $\frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)}$ has to be calculated only once as long as the degrees of freedom k doesn't change. This can be done similar to program T.
The value of the distribution has to be calculated in program D. The other programs don't change.

For the F-distribution you have to calculate:

$[Image: b8d9434be6497007b94e7ae5f22c34a3.png]$

with parameters $d_1$ and $d_2$.

To calculate the beta function we can use:

$[Image: fcd016ccf3e5aca1621c9b41f077a339.png]$

Together with the program for the Student's t-distribution it should be possible to create the missing programs.

Cheers
Thomas

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-26-2015 10:29 PM

[attachment=2697][attachment=2697][attachment=2697][attachment=2697]

(10-26-2015 08:30 PM)Dieter Wrote:
(10-26-2015 04:41 PM)PedroLeiva Wrote: Thank you very much for your recommendation, but Namir programs are not what I need. In Appendix A, there are detailed calculations based on estimates of the inverse probability of each 4 statistical functions. I also need to calculate probability according to a value of each function.

OK – what exactly do you need? The CDFs for Normal, Student, Chi² and Fisher's F distribution? Also their quantile (inverse) functions? What accuracy is required?

BTW, the 34s does all this. I provided some code for these four quantile functions which are evaluated with mostly 30+ digit accuracy. This is possible due to CDFs that are (mostly) equally accurate. These CDFs in turn are calculated using incomplete Beta and Gamma functions, coded by Pauli.

(10-26-2015 04:41 PM)PedroLeiva Wrote: HP 67 Statistical Distributions programs (Normal, Chi Square, "F" and "t" Student) a perfect for this purpose.

Could you provide a pointer to this package? How accurate are the results? Do these programs also calculate the inverse (quantile) functions?
If you really appreciate these HP67 programs so much it should be easy to translate them for the 35s.

Dieter

Attached are three HP 67 programs, with descriptions and formulas, some examples or test cases, and the program itself (keystrokes sequence). I need exactly what they allow to calculate, and with the same precision.

I am not a specialist in mathematics; I simply use a lot of statistic and want to become independent of the function tables.

RE: HP 35s Statistical Distributions Functions - Dieter - 10-26-2015 11:04 PM

(10-26-2015 10:29 PM)PedroLeiva Wrote: Attached are three HP 67 programs, with descriptions and formulas, some examples or test cases, and the program itself (keystrokes sequence). I need exactly what they allow to calculate, and with the same precision.

So you only need the PDFs and CDFs – that's simple and straightforwand on the 35s. Actually the implementation is easier as the 35s features a Gamma function.

(10-26-2015 10:29 PM)PedroLeiva Wrote: I am not a specialist in mathematics; I simply use a lot of statistic and want to become independent of the function tables.

Get a 34s. ;-)

Dieter

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-27-2015 11:59 AM

(10-26-2015 09:28 PM)Thomas Klemm Wrote: For the Chi-squared distribution you just have to calculate a different PDF:

$[Image: a107033ede4e5de93c3b7817fc5aeb33.png]$

The factor $\frac{1}{2^{\frac{k}{2}}\Gamma\left(\frac{k}{2}\right)}$ has to be calculated only once as long as the degrees of freedom k doesn't change. This can be done similar to program T.
The value of the distribution has to be calculated in program D. The other programs don't change.

For the F-distribution you have to calculate:

$[Image: b8d9434be6497007b94e7ae5f22c34a3.png]$

with parameters $d_1$ and $d_2$.

To calculate the beta function we can use:

$[Image: fcd016ccf3e5aca1621c9b41f077a339.png]$

Together with the program for the Student's t-distribution it should be possible to create the missing programs.

Cheers
Thomas

Thomas, with time and patience I will try to program my HP 35s with the suggestions that you just gave me. Thank you

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-27-2015 02:29 PM

(10-26-2015 08:39 PM)Thomas Klemm Wrote: Here's a program for the Student's t-distribution:

Code:

T001 LBL T ; Student's t T002 INPUT K ; degrees of freedom T003 1 T004 - T005 2 T006 / T007 ! ; Γ((k+1)/2) = (k-1)/2 ! T008 RCL K T009 2 T010 / T011 1 T012 - T013 ! ; Γ(k/2) = (k/2 - 1)! T014 / T015 PI T016 RCL* K T017 SQRT T018 / T019 STO G T020 RTN I001 LBL I ; inverse CDF ( CDF -- X ) I002 INPUT A ; area I003 0.5 I004 - I005 RCL/ G I006 STO X ; initial guess I007 XEQ C003 ; calulate cdf(X) I008 RCL- A ; we want this difference to become 0 I009 RCL X I010 STO T I011 R↓ I012 XEQ D001 ; derivative at X = pdf(X) I013 / I014 STO- X ; improved guess I015 ABS I016 RCL E ; small like 1e-6 I017 x<y? ; close enough? I018 GTO I006 ; not yet I019 RCL X I020 RTN C001 LBL C ; CDF C002 INPUT X C003 0 ; lower limit C004 RCL X ; upper limit C005 FN= D ; the distribution C006 ∫FN d T ; integrate pdf C007 0.5 ; we add 0.5 since C008 + ; the lower limit was 0 C009 RTN P001 LBL P ; PDF P002 INPUT X P003 STO T D001 LBL D ; Distribution: ( T -- PDF ) D002 RCL T D003 x^2 D004 RCL/ K D005 1 D006 + ; 1 + t^2/k D007 LASTx D008 RCL+ K D009 2 D010 / ; (1 + k)/2 D011 y^x D012 RCL G D013 x<>y D014 / D015 RTN

You have to initialize the program with the degrees of freedom:

XEQ T ENTER
K?
7
R/S

To calculate the PDF use the program P:

XEQ P ENTER
X?
1.6
R/S
0.110665

PDF[StudentTDistribution[7],1.6]

To calculate the CDF use the program C:

XEQ C ENTER
X?
1.5
R/S
0.911351

CDF[StudentTDistribution[7],1.5]

For the next program set the accuracy to a reasonable value:

1E-6
STO E

To calculate the Inverse CDF use the program I:

XEQ I ENTER
X?
0.99
R/S
2.997951

HTH
Thomas

Thomas, two questions. What is the sequence of keys to enter STEP C006. To run the label I, the question is A ?, rather X?, right?

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-27-2015 04:48 PM

(10-27-2015 02:29 PM)PedroLeiva Wrote: What is the sequence of keys to enter STEP C006.

Code:

C006 ∫FN d T     ; integrate pdf

That's ∫ followed by T.

[ ↰ ] [EQN] [ 9 ]

Quote:To run the label I, the question is A ?, rather X?, right?

Input is area (that's why A) and the result is $x$ so that $\int_{-\infty}^x pdf(t)\,dt = A$.
You may rename it to something different of course. Just don't use variables that are already used in the program. Thus X might be confusing.

Cheers
Thomas

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-27-2015 05:22 PM

(10-27-2015 04:48 PM)Thomas Klemm Wrote:
(10-27-2015 02:29 PM)PedroLeiva Wrote: What is the sequence of keys to enter STEP C006.

Code:

C006 ∫FN d T ; integrate pdf

That's ∫ followed by T.

[ ↰ ] [EQN] [ 9 ]

Quote:To run the label I, the question is A ?, rather X?, right?

Input is area (that's why A) and the result is $x$ so that $\int_{-\infty}^x pdf(t)\,dt = A$.
You may rename it to something different of course. Just don't use variables that are already used in the program. Thus X might be confusing.

Cheers
Thomas

Something is not working well. When I run the routine I (XEQ I ENTER), with a 1E-6 E value in register E, the program do not stop integration!

RE: HP 35s Statistical Distributions Functions - Dieter - 10-27-2015 07:32 PM

(10-27-2015 05:22 PM)PedroLeiva Wrote: Something is not working well. When I run the routine I (XEQ I ENTER), with a 1E-6 E value in register E, the program do not stop integration!

The 1E–6 error threshold in E is irrelevant for the integration routine. The 35s "Integrate" function stops if the result is correct to display precision. You probably have set your 35s to ALL display mode, or maybe FIX 9 or SCI 9, so the 35s calculates the integral until is thinks it has all (!) 12 digits (resp. 9) correct. This will take some time. Solution: simply set something like FIX 4 or FIX 6 (you may add this step to the program right before the Integrate command), so you get 4 resp. 6 correct decimals.

BTW, at the moment I am looking at an 35s implementation of the HP67 routines.

Dieter

RE: HP 35s Statistical Distributions Functions - Dieter - 10-27-2015 10:25 PM

(10-27-2015 07:32 PM)Dieter Wrote: BTW, at the moment I am looking at an 35s implementation of the HP67 routines.

OK, here is a first experimental version of the Chi² algorithm used in the HP67 Stat Pac. PDF and CDF are calculated simultaneously, and both are returned in Y and X. Here is the code:

Code:

C001  LBL C

C002  INPUT J

C003  e

C004  IP

C005  STO T

C006  ÷

C007  !

C008  STO G

C009  INPUT X

C010  XEQ C013

C011  STOP

C012  GTO C009

C013  RCL X

C014  RCL÷ T

C015  +/-

C016  e^x

C017  LASTx

C018  +/-

C019  RCL J

C020  RCL÷ T

C021  y^x

C022  x

C023  RCL÷ G

C024  STO A

C025  RCL J

C026  STO B

C027  SGN

C028  ENTER

C029  -

C030  ENTER

C031  ENTER

C032  LASTx

C033  RCLx X

C034  RCL B

C035  RCL+ T

C036  STO B

C037  ÷

C038  +

C039  x≠y?

C040  GTO C030

C041  RCL A

C042  RCL X

C043  x=0?

C044  e^x

C045  ÷

C046  RCLx J

C047  RCL÷ T

C048  x<>y

C049  RCLx A

C050  RCL+ A

C051  RTN

Usage (results as shown in FIX 9 mode):

Code:

[XEQ] C [ENTER]             J=?

Enter degrees of freedom

    10  [R/S]               X=?

Enter random variable X

    6,3 [R/S]               0,087896860   // =PDF(6,3|10)

                            0,210539719   // =CDF(6,3|10)

Do another calculation

        [R/S]               X=?

    20  [R/S]               0,009458319   // =PDF(20|10)

                            0,970747312   // =CDF(20|10)

As usual, comments, suggestions and error reports are welcome. ;-)

Since the main computation routine (C013 ff.) returns both the PDF and the CDF, adding a function that determines the inverse (quantile) is rather simple. I tried an approach using a second-order (Halley) method and it worked well for me. This required about 30 additional lines of code.

Dieter

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-27-2015 10:29 PM

(10-27-2015 05:22 PM)PedroLeiva Wrote: the program do not stop integration!

Let's try to pin that down. Dieter already gave a hint: the accuracy in register E and the display precision should match.
But since the program I uses program C I first want to know if you get the same result for this example:

XEQ T ENTER
K?
7
R/S

XEQ C ENTER
X?
1.5
R/S
0.911351

This is just to make sure that program C works fine.

If that's okay please make sure to set FIX 3 and E = 1E-3 before running the next example.

XEQ I ENTER
X?
0.99
R/S
2.997

It should flicker about four times INTEGRATING and then display the result.

If the integration doesn't stop then please add an R/S command between lines I007 and I008 and report both the value displayed and the contents of register X (use VIEW X). Just use R/S to run the next round of the loop for a couple of times.

Hopefully this helps to find the problem.

Cheers
Thomas

PS: Thanks for posting the HP-67 programs. It's always interesting to learn something new. For instance I didn't know about the formula used to calculate the CDF of the Student's t-distribution. But then I must admit that I'm not familiar with statistics. Still I'm curious how they came up with that formula.

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-27-2015 11:41 PM

(10-27-2015 10:29 PM)Thomas Klemm Wrote:
(10-27-2015 05:22 PM)PedroLeiva Wrote: the program do not stop integration!

Let's try to pin that down. Dieter already gave a hint: the accuracy in register E and the display precision should match.
But since the program I uses program C I first want to know if you get the same result for this example:

XEQ T ENTER
K?
7
R/S

XEQ C ENTER
X?
1.5
R/S
0.911351

This is just to make sure that program C works fine.

If that's okay please make sure to set FIX 3 and E = 1E-3 before running the next example.

XEQ I ENTER
X?
0.99
R/S
2.997

It should flicker about four times INTEGRATING and then display the result.

If the integration doesn't stop then please add an R/S command between lines I007 and I008 and report both the value displayed and the contents of register X (use VIEW X). Just use R/S to run the next round of the loop for a couple of times.

Hopefully this helps to find the problem.

Cheers
Thomas

PS: Thanks for posting the HP-67 programs. It's always interesting to learn something new. For instance I didn't know about the formula used to calculate the CDF of the Student's t-distribution. But then I must admit that I'm not familiar with statistics. Still I'm curious how they came up with that formula.

For routing T the value after R/S is 0.385
For routing C the value is the same as yours, 0.911351
The program did not stop integration again. After R/S adition the figures are:
1° time: x=0.878, y=3.782 E-2; variable X= 1.273
2° time: x=0.5, y=3.850 E-7; variable X= 0.001
# times: the same as 2° time

RE: HP 35s Statistical Distributions Functions - PedroLeiva - 10-27-2015 11:48 PM

(10-27-2015 07:32 PM)Dieter Wrote:
(10-27-2015 05:22 PM)PedroLeiva Wrote: Something is not working well. When I run the routine I (XEQ I ENTER), with a 1E-6 E value in register E, the program do not stop integration!

The 1E–6 error threshold in E is irrelevant for the integration routine. The 35s "Integrate" function stops if the result is correct to display precision. You probably have set your 35s to ALL display mode, or maybe FIX 9 or SCI 9, so the 35s calculates the integral until is thinks it has all (!) 12 digits (resp. 9) correct. This will take some time. Solution: simply set something like FIX 4 or FIX 6 (you may add this step to the program right before the Integrate command), so you get 4 resp. 6 correct decimals.

BTW, at the moment I am looking at an 35s implementation of the HP67 routines.

Dieter

Thank you for your interest Dieter. I am ussing FIX 3, but the problem persist

RE: HP 35s Statistical Distributions Functions - Thomas Klemm - 10-28-2015 12:11 AM

(10-27-2015 11:41 PM)PedroLeiva Wrote: 2° time: x=0.5, y=3.850 E-7; variable X= 0.001

That was helpful. The value of E got copied to X.
Turns out there was a bug in this line:

Code:

I018 GTO I007    ; not yet

Should go to line I007 instead of I006.

Cheers
Thomas

PS: Fixed it in the listing.