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?

HP Forums > HP Software Libraries > General Software Library > (35S) Statistical Distributions Functions

You're currently viewing a stripped down version of our content. View the full version with proper formatting.

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?

10-25-2015, 01:29 PM

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

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

- 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

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

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.

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

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

10-26-2015, 08:39 PM

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

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

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

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

10-26-2015, 09:28 PM

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

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:

with parameters \(d_1\) and \(d_2\).

To calculate the beta function we can use:

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

Cheers

Thomas

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:

with parameters \(d_1\) and \(d_2\).

To calculate the beta function we can use:

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

Cheers

Thomas

10-26-2015, 10:29 PM

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

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.

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

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

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:

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:

with parameters \(d_1\) and \(d_2\).

To calculate the beta function we can use:

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

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?

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

10-27-2015, 05:22 PM

(10-27-2015 04:48 PM)Thomas Klemm 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!(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

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

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

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.

10-27-2015, 11:41 PM

(10-27-2015 10:29 PM)Thomas Klemm Wrote: [ -> ]For routing T the value after R/S is 0.385(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 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

10-27-2015, 11:48 PM

(10-27-2015 07:32 PM)Dieter Wrote: [ -> ]Thank you for your interest Dieter. I am ussing FIX 3, but the problem persist(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

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.

HP Forums > HP Software Libraries > General Software Library > (35S) Statistical Distributions Functions

**HP Forums:**https://www.hpmuseum.org/forum/index.php**:**

Powered By MyBB, © 2002-2019 MyBB Group