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 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?
(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.
(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.
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
(10-26-2015 08:39 PM)Thomas Klemm Wrote: [ -> ]Here's a program for the Student's t-distribution:
(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
(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.
(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.
(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
(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
(10-27-2015 02:29 PM)PedroLeiva Wrote: [ -> ]What is the sequence of keys to enter STEP C006.
C006 ∫FN d T ; integrate pdf
Quote:To run the label I, the question is A ?, rather X?, right?
(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 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!
(10-27-2015 07:32 PM)Dieter Wrote: [ -> ]BTW, at the moment I am looking at an 35s implementation of the HP67 routines.
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
[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)
(10-27-2015 05:22 PM)PedroLeiva Wrote: [ -> ]the program do not stop integration!
(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.
(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-27-2015 11:41 PM)PedroLeiva Wrote: [ -> ]2° time: x=0.5, y=3.850 E-7; variable X= 0.001
I018 GTO I007 ; not yet