HP Forums
(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.