Post Reply 
(35S) Statistical Distributions Functions
10-28-2015, 12:24 AM
Post: #21
RE: HP 35s Statistical Distributions Functions
(10-27-2015 10:25 PM)Dieter Wrote:  
(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

Dieter, the program works perfect!. Please include the function for the inverse, no mather how long it takes
Find all posts by this user
Quote this message in a reply
10-28-2015, 12:36 AM (This post was last modified: 10-28-2015 12:38 AM by PedroLeiva.)
Post: #22
RE: HP 35s Statistical Distributions Functions
(10-28-2015 12:11 AM)Thomas Klemm Wrote:  
(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.

Thomas, I erase the STOP instruction, and corrected for "I018 GTO I007" so the program began to work perfect. Thank you
Find all posts by this user
Quote this message in a reply
10-28-2015, 01:05 AM
Post: #23
RE: HP 35s Statistical Distributions Functions
(10-27-2015 10:25 PM)Dieter Wrote:  I tried an approach using a second-order (Halley) method and it worked well for me. This required about 30 additional lines of code.

Or then you can just use the built-in solver.

Cheers
Thomas
Find all posts by this user
Quote this message in a reply
10-28-2015, 01:40 AM
Post: #24
RE: HP 35s Statistical Distributions Functions
(10-28-2015 12:36 AM)PedroLeiva Wrote:  I erase the STOP instruction, and corrected for "I018 GTO I007" so the program began to work perfect.

That was probably the most distant remote debugging session I ever had.

I've added line I005 to improve the initial guess only after writing the listing of the program. This shifted all the following commands down one line. And then I forgot to update the GTO command in the listing. We all know that the Go To Statement is Considered Harmful. Did we really need another proof?

Cheers
Thomas
Find all posts by this user
Quote this message in a reply
10-28-2015, 11:55 AM (This post was last modified: 10-28-2015 12:07 PM by Dieter.)
Post: #25
RE: HP 35s Statistical Distributions Functions
(10-28-2015 01:05 AM)Thomas Klemm Wrote:  Or then you can just use the built-in solver.

A dedicated solver is more effective and requires just one initial guess. ;-)

But maybe it is possible to give a lower and higher limit for the quantile and feed this to the solver. For all four distributions, that is. With some slight tweaks (e.g. exit if a certain accuracy level is reached since the CDF is not absolutely exact anyway) this might work as well.

EDIT: the built-in solver probably will not work here. It is known that the 35s hat problems when solving programs containing loops. So a dedicated solver seems to be the only solution.

Dieter
Find all posts by this user
Quote this message in a reply
10-28-2015, 12:34 PM (This post was last modified: 10-28-2015 12:58 PM by Dieter.)
Post: #26
RE: HP 35s Statistical Distributions Functions
(10-28-2015 12:24 AM)PedroLeiva Wrote:  Dieter, the program works perfect!. Please include the function for the inverse, no mather how long it takes

OK, here is a first version with a rough initial guess. At the X=? prompt simply enter a negative value if you want to compute the quantile for that probabily. If you already have entered the previous C program, simply add lines C010 and C011 and continue with line C054.

Code:
C001  LBL C
C002  INPUT J
C003  e
C004  IP
C005  STO T
C006  ÷
C007  !
C008  STO G
C009  INPUT X
C010  x<0?
C011  GTO C054
C012  XEQ C015
C013  STOP
C014  GTO C009
C015  RCL X
C016  RCL÷ T
C017  +/-
C018  e^x
C019  LASTx
C020  +/-
C021  RCL J
C022  RCL÷ T
C023  y^x
C024  x
C025  RCL÷ G
C026  STO A
C027  RCL J
C028  STO B
C029  SGN
C030  ENTER
C031  -
C032  ENTER
C033  ENTER
C034  LASTx
C035  RCLx X
C036  RCL B
C037  RCL+ T
C038  STO B
C039  ÷
C040  +
C041  x≠y?
C042  GTO C032
C043  RCL A
C044  RCL X
C045  x=0?
C046  e^x
C047  ÷
C048  RCLx J
C049  RCL÷ T
C050  x<>y
C051  RCLx A
C052  RCL+ A
C053  RTN
C054  +/-
C055  STO P
C056  RCLx J
C057  RCL÷ T
C058  RCLx G
C059  RCL T
C060  RCL÷ J
C061  y^x
C062  Pi
C063  x
C064  STO X
C065  XEQ C015
C066  ENTER
C067  RCL÷ P
C068  LN
C069  x
C070  ENTER
C071  ENTER
C072  RCL J
C073  RCL- X
C074  RCL- T
C075  RCL÷ X
C076  RCL÷ T
C077  RCL÷ T
C078  x
C079  R↑
C080  -
C081  ÷
C082  RCL+ X
C083  x<> X
C084  RCL X
C085  STOP
C086  GTO C064

Usage: At the X=? prompt enter the negative probability for which you want the Chi² quantile. The program will then display two estimates in Y and X (the previous and the current), so that you can see the iteration converge. Pressing R/S yields the next approximation.

Code:
[XEQ] C [ENTER]             J=?

Enter degrees of freedom
    10  [R/S]               X=?

Enter negative probability
  -0,95 [R/S]               11,176989523
                            15,310270784

        [R/S]               15,310270784
                            17,970799176

        [R/S]               17,970799176
                            18,305703292

        [R/S]               18,305703292
                            18,307038039

If 18,307 is sufficiently accurate you may stop here.
Otherwise continue:

Code:
        [R/S]               18,307038039
                            18,307038054

        [R/S]               18,307038054
                            18,307038054

When the two approximations are displayed you may also enter your own guess before pressing R/S, and the approximation will continue with this value.

NB: This program is not very sophisticated, and it does not match what's inside the 34s. For instance, the latter evaluates CDFs for large x via P(x) = 1–Q(x), and also the quantile function has a much, much better initial guess. It took me quite a while to get the guess within <10% of the true value. ;-) For instance, the first guess for the case above is 18,2..., compared to the exact value 18,307...

Some time ago I wrote 35s/41C programs for the Normal distribution. They evaluate the PDF, one- and two-sided CDFs and also one- and two-sided quantiles with several tweaks to ensure best possible accuracy and efficiency. But these are about three times the size of the small Chi² progam posted above.

Dieter
Find all posts by this user
Quote this message in a reply
10-28-2015, 06:21 PM
Post: #27
RE: HP 35s Statistical Distributions Functions
I think we have tight-fitting Student t-distribution. I attached the pdf with the program that Thomas wrote, five calculation examples (test cases), an outline of the distribición and an instruction table for input and output of data.

Comments are welcome. Pedro


Attached File(s)
.pdf  HP 35s - Student´s t - Distribution.pdf (Size: 77.62 KB / Downloads: 14)
Find all posts by this user
Quote this message in a reply
10-28-2015, 06:59 PM
Post: #28
RE: HP 35s Statistical Distributions Functions
(10-28-2015 12:34 PM)Dieter Wrote:  
(10-28-2015 12:24 AM)PedroLeiva Wrote:  Dieter, the program works perfect!. Please include the function for the inverse, no mather how long it takes

OK, here is a first version with a rough initial guess. At the X=? prompt simply enter a negative value if you want to compute the quantile for that probabily. If you already have entered the previous C program, simply add lines C010 and C011 and continue with line C054.

Code:
C001  LBL C
C002  INPUT J
C003  e
C004  IP
C005  STO T
C006  ÷
C007  !
C008  STO G
C009  INPUT X
C010  x<0?
C011  GTO C054
C012  XEQ C015
C013  STOP
C014  GTO C009
C015  RCL X
C016  RCL÷ T
C017  +/-
C018  e^x
C019  LASTx
C020  +/-
C021  RCL J
C022  RCL÷ T
C023  y^x
C024  x
C025  RCL÷ G
C026  STO A
C027  RCL J
C028  STO B
C029  SGN
C030  ENTER
C031  -
C032  ENTER
C033  ENTER
C034  LASTx
C035  RCLx X
C036  RCL B
C037  RCL+ T
C038  STO B
C039  ÷
C040  +
C041  x≠y?
C042  GTO C032
C043  RCL A
C044  RCL X
C045  x=0?
C046  e^x
C047  ÷
C048  RCLx J
C049  RCL÷ T
C050  x<>y
C051  RCLx A
C052  RCL+ A
C053  RTN
C054  +/-
C055  STO P
C056  RCLx J
C057  RCL÷ T
C058  RCLx G
C059  RCL T
C060  RCL÷ J
C061  y^x
C062  Pi
C063  x
C064  STO X
C065  XEQ C015
C066  ENTER
C067  RCL÷ P
C068  LN
C069  x
C070  ENTER
C071  ENTER
C072  RCL J
C073  RCL- X
C074  RCL- T
C075  RCL÷ X
C076  RCL÷ T
C077  RCL÷ T
C078  x
C079  R↑
C080  -
C081  ÷
C082  RCL+ X
C083  x<> X
C084  RCL X
C085  STOP
C086  GTO C064

Usage: At the X=? prompt enter the negative probability for which you want the Chi² quantile. The program will then display two estimates in Y and X (the previous and the current), so that you can see the iteration converge. Pressing R/S yields the next approximation.

Code:
[XEQ] C [ENTER]             J=?

Enter degrees of freedom
    10  [R/S]               X=?

Enter negative probability
  -0,95 [R/S]               11,176989523
                            15,310270784

        [R/S]               15,310270784
                            17,970799176

        [R/S]               17,970799176
                            18,305703292

        [R/S]               18,305703292
                            18,307038039

If 18,307 is sufficiently accurate you may stop here.
Otherwise continue:

Code:
        [R/S]               18,307038039
                            18,307038054

        [R/S]               18,307038054
                            18,307038054

When the two approximations are displayed you may also enter your own guess before pressing R/S, and the approximation will continue with this value.

NB: This program is not very sophisticated, and it does not match what's inside the 34s. For instance, the latter evaluates CDFs for large x via P(x) = 1–Q(x), and also the quantile function has a much, much better initial guess. It took me quite a while to get the guess within <10% of the true value. ;-) For instance, the first guess for the case above is 18,2..., compared to the exact value 18,307...

Some time ago I wrote 35s/41C programs for the Normal distribution. They evaluate the PDF, one- and two-sided CDFs and also one- and two-sided quantiles with several tweaks to ensure best possible accuracy and efficiency. But these are about three times the size of the small Chi² progam posted above.

Dieter

Dieter, all perfect. I got the same figures as you!. I will prepare a PDF with programme keytroke sequence, test cases and instructions for in and output data
Find all posts by this user
Quote this message in a reply
10-28-2015, 07:54 PM (This post was last modified: 10-28-2015 08:00 PM by Dieter.)
Post: #29
RE: HP 35s Statistical Distributions Functions
(10-28-2015 06:21 PM)PedroLeiva Wrote:  I think we have tight-fitting Student t-distribution. I attached the pdf with the program that Thomas wrote, five calculation examples (test cases), an outline of the distribición and an instruction table for input and output of data.

For the record, here is a version without using the Integrate function. It is based on algorithm AS 3, published as a Fortran program in Applied Statistics (1968). This version includes an improved input routine for the degrees of freedom in that it checks whether the input is a positive integer. Again, both the PDF and CDF are returned on the stack (in Y and X). Code for the inverse is not yet included.

Code:
T001  LBL T
T002  ALL
T003  INPUT J
T004  ENTER
T005  ABS
T006  IP
T007  x=y?
T008  x=0?
T009  GTO T001
T010  e
T011  IP
T012  STO T
T013  SGN
T014  -
T015  RCL÷ T
T016  !
T017  LASTx
T018  RCL T
T019  1/x
T020  -
T021  !
T022  ÷
T023  Pi
T024  RCLx J
T025  √x
T026  ÷
T027  STO G
T028  FIX 9
T029  INPUT X
T030  RCL J
T031  √x
T032  ÷
T033  STO A
T034  RCL J
T035  RCL X
T036  x²
T037  RCL+ J
T038  ÷
T039  STO B
T040  RCL J
T041  RCL T
T042  RMDR
T043  RCL+ T
T044  STO K
T045  SGN
T046  STO C
T047  ENTER
T048  ENTER
T049  RCL K
T050  RCL- J
T051  x≥0?
T052  GTO T066
T053  SGN
T054  RCL+ K
T055  RCLx B
T056  RCL÷ K
T057  RCLx C
T058  STO C
T059  +
T060  x=y?
T061  GTO T066
T062  RCL T
T063  STO+ K
T064  R↓
T065  GTO T047
T066  R↓
T067  STO C
T068  RCL J
T069  RCL T
T070  RMDR
T071  x=0?
T072  GTO T090
T073  RCL- J
T074  x=0?
T075  STO C
T076  RCL A
T077  RCLx B
T078  RCLx C
T079  Pi
T080  ÷
T081  RCL A
T082  ATAN
T083  RCL T
T084  +/-
T085  SGN
T086  ACOS
T087  ÷
T088  +
T089  GTO T095
T090  RCL B
T091  √x
T092  RCLx A
T093  RCLx C
T094  RCL÷ T
T095  RCL T
T096  1/x
T097  +
T098  RCL X
T099  x²
T100  RCL÷ J
T101  RCL T
T102  SGN
T103  +
T104  LASTx
T105  RCL+ J
T106  RCL÷ T
T107  +/-
T108  y^x
T109  RCLx G
T110  x<>y
T111  STOP
T112  GTO T028

Here is Thomas' sample case:

Code:
[XEQ] T [ENTER]             J=?

Enter degrees of freedom
     7  [R/S]               X=?

Enter random variable X
    1,5 [R/S]               0,126263061   // =PDF(1,5|7)
                            0,911350757   // =CDF(1,5|7)
Do another calculation
        [R/S]               X=?
    -2  [R/S]               0,063135337   // =PDF(-2|7)
                            0,042809664   // =CDF(-2|7)

Final remark: like the other programs discussed in this thread, this is a rather simple and straightforward implementation that works fine for most practical cases. However, far out in the distribution tails the results may lose accuracy due to rounding, and eventually values may round to 0 or 1. For instance, for 7 d.o.f. and x=–100 you will get a CDF of 2,000...E–12 while the exact result is 1,318 E–12, and for x=0 the returned result is zero while it actually is 1,69 E–17. If you want more or less exact values even for such extreme arguments, a different approach is required.

Dieter
Find all posts by this user
Quote this message in a reply
10-28-2015, 10:34 PM
Post: #30
RE: HP 35s Statistical Distributions Functions
(10-28-2015 07:54 PM)Dieter Wrote:  
(10-28-2015 06:21 PM)PedroLeiva Wrote:  I think we have tight-fitting Student t-distribution. I attached the pdf with the program that Thomas wrote, five calculation examples (test cases), an outline of the distribición and an instruction table for input and output of data.

For the record, here is a version without using the Integrate function. It is based on algorithm AS 3, published as a Fortran program in Applied Statistics (1968). This version includes an improved input routine for the degrees of freedom in that it checks whether the input is a positive integer. Again, both the PDF and CDF are returned on the stack (in Y and X). Code for the inverse is not yet included.

Code:
T001  LBL T
T002  ALL
T003  INPUT J
T004  ENTER
T005  ABS
T006  IP
T007  x=y?
T008  x=0?
T009  GTO T001
T010  e
T011  IP
T012  STO T
T013  SGN
T014  -
T015  RCL÷ T
T016  !
T017  LASTx
T018  RCL T
T019  1/x
T020  -
T021  !
T022  ÷
T023  Pi
T024  RCLx J
T025  √x
T026  ÷
T027  STO G
T028  FIX 9
T029  INPUT X
T030  RCL J
T031  √x
T032  ÷
T033  STO A
T034  RCL J
T035  RCL X
T036  x²
T037  RCL+ J
T038  ÷
T039  STO B
T040  RCL J
T041  RCL T
T042  RMDR
T043  RCL+ T
T044  STO K
T045  SGN
T046  STO C
T047  ENTER
T048  ENTER
T049  RCL K
T050  RCL- J
T051  x≥0?
T052  GTO T066
T053  SGN
T054  RCL+ K
T055  RCLx B
T056  RCL÷ K
T057  RCLx C
T058  STO C
T059  +
T060  x=y?
T061  GTO T066
T062  RCL T
T063  STO+ K
T064  R↓
T065  GTO T047
T066  R↓
T067  STO C
T068  RCL J
T069  RCL T
T070  RMDR
T071  x=0?
T072  GTO T090
T073  RCL- J
T074  x=0?
T075  STO C
T076  RCL A
T077  RCLx B
T078  RCLx C
T079  Pi
T080  ÷
T081  RCL A
T082  ATAN
T083  RCL T
T084  +/-
T085  SGN
T086  ACOS
T087  ÷
T088  +
T089  GTO T095
T090  RCL B
T091  √x
T092  RCLx A
T093  RCLx C
T094  RCL÷ T
T095  RCL T
T096  1/x
T097  +
T098  RCL X
T099  x²
T100  RCL÷ J
T101  RCL T
T102  SGN
T103  +
T104  LASTx
T105  RCL+ J
T106  RCL÷ T
T107  +/-
T108  y^x
T109  RCLx G
T110  x<>y
T111  STOP
T112  GTO T028

Here is Thomas' sample case:

Code:
[XEQ] T [ENTER]             J=?

Enter degrees of freedom
     7  [R/S]               X=?

Enter random variable X
    1,5 [R/S]               0,126263061   // =PDF(1,5|7)
                            0,911350757   // =CDF(1,5|7)
Do another calculation
        [R/S]               X=?
    -2  [R/S]               0,063135337   // =PDF(-2|7)
                            0,042809664   // =CDF(-2|7)

Final remark: like the other programs discussed in this thread, this is a rather simple and straightforward implementation that works fine for most practical cases. However, far out in the distribution tails the results may lose accuracy due to rounding, and eventually values may round to 0 or 1. For instance, for 7 d.o.f. and x=–100 you will get a CDF of 2,000...E–12 while the exact result is 1,318 E–12, and for x=0 the returned result is zero while it actually is 1,69 E–17. If you want more or less exact values even for such extreme arguments, a different approach is required.

Dieter

Dieter, perfect!. I will try this option also. Thanks
I am preparing a compiled program pdf for Chi-Square to publish in the Forum
Find all posts by this user
Quote this message in a reply
10-28-2015, 10:45 PM
Post: #31
RE: HP 35s Statistical Distributions Functions
I leave the program compilation that Dieter wrote for Chi-Square Distribution, with examples and a list of programming steps.

Much I wish you could make the attempt with F Distribution, just to complete de Statistical Distribution Functions. Remember that Normal and Inverse Distribution is already included in the HP 35s Manual!


Attached File(s)
.pdf  HP 35s - Chi-Square Distribution.pdf (Size: 95.56 KB / Downloads: 12)
Find all posts by this user
Quote this message in a reply
10-29-2015, 12:09 AM
Post: #32
RE: HP 35s Statistical Distributions Functions
(10-28-2015 10:45 PM)PedroLeiva Wrote:  I leave the program compilation that Dieter wrote for Chi-Square Distribution, with examples and a list of programming steps.

Much I wish you could make the attempt with F Distribution, just to complete de Statistical Distribution Functions. Remember that Normal and Inverse Distribution is already included in the HP 35s Manual!

Here is the program for Normal&Inverse-Normal Distributions. Just for those who will be interested in this subject


Attached File(s)
.pdf  Hp 35s - Normal&Inverse-Normal Distributions.pdf (Size: 125.71 KB / Downloads: 13)
Find all posts by this user
Quote this message in a reply
10-29-2015, 07:01 AM
Post: #33
RE: HP 35s Statistical Distributions Functions
(10-29-2015 12:09 AM)PedroLeiva Wrote:  Here is the program for Normal&Inverse-Normal Distributions. Just for those who will be interested in this subject

This is a sample program from the 35s manual. It works, but there is room for improvement. For instance the HP15C Advanced Functions Handbook has a nice 15C program that uses the Integrate function with variable transformation which works much better in the distribution tails.

Personally, I'd prefer a different approach with the classic combination of series expansion and continued fraction to evaluate the Gauss integral. For the inverse, rational approximations with various accuracy levels can be given (yes, there is much more that Hastings' three-digit method as implemented in several HP software pacs). Maybe I will eventually post an example here.

Dieter
Find all posts by this user
Quote this message in a reply
10-29-2015, 07:18 AM (This post was last modified: 10-29-2015 07:29 AM by Dieter.)
Post: #34
RE: HP 35s Statistical Distributions Functions
(10-28-2015 10:45 PM)PedroLeiva Wrote:  I leave the program compilation that Dieter wrote for Chi-Square Distribution, with examples and a list of programming steps.

Looks very nice. FTR: There are special fonts available that resemble the keys and the display of various HP calculators (simple 7 segment LCD, dot matrix, 14 segment HP41, etc.). I remember a font package by Luiz Viera and there also is a zip-file at educalc.net that inclues a 33s key font. Simply search for HP calculator fonts

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. For the time being you may include the improved J input routine from the Student program.

Dieter
Find all posts by this user
Quote this message in a reply
10-29-2015, 11:56 PM
Post: #35
RE: HP 35s Statistical Distributions Functions
(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.

I now tried this with the Chi² and Student programs. Both calculate a quite good first guess for the quantile which is based on a simple approximation of the Normal quantile. The previous Chi² example (10 dof and p=0,95) converges after just three iterations. The initial estimate was 18,24 with 18,307 being the correct result. In most cases, both for Chi² and Student, merely 3...5 iterations will do. This approach requires one or two additional subroutine(s) for the Normal estimate, and the complete Student program (including the inverse) finally has about 200 steps.

Dieter
Find all posts by this user
Quote this message in a reply
10-30-2015, 02:36 AM
Post: #36
RE: HP 35s Statistical Distributions Functions
(10-28-2015 07:54 PM)Dieter Wrote:  
(10-28-2015 06:21 PM)PedroLeiva Wrote:  I think we have tight-fitting Student t-distribution. I attached the pdf with the program that Thomas wrote, five calculation examples (test cases), an outline of the distribición and an instruction table for input and output of data.

For the record, here is a version without using the Integrate function. It is based on algorithm AS 3, published as a Fortran program in Applied Statistics (1968). This version includes an improved input routine for the degrees of freedom in that it checks whether the input is a positive integer. Again, both the PDF and CDF are returned on the stack (in Y and X). Code for the inverse is not yet included.

Code:
T001  LBL T
T002  ALL
T003  INPUT J
T004  ENTER
T005  ABS
T006  IP
T007  x=y?
T008  x=0?
T009  GTO T001
T010  e
T011  IP
T012  STO T
T013  SGN
T014  -
T015  RCL÷ T
T016  !
T017  LASTx
T018  RCL T
T019  1/x
T020  -
T021  !
T022  ÷
T023  Pi
T024  RCLx J
T025  √x
T026  ÷
T027  STO G
T028  FIX 9
T029  INPUT X
T030  RCL J
T031  √x
T032  ÷
T033  STO A
T034  RCL J
T035  RCL X
T036  x²
T037  RCL+ J
T038  ÷
T039  STO B
T040  RCL J
T041  RCL T
T042  RMDR
T043  RCL+ T
T044  STO K
T045  SGN
T046  STO C
T047  ENTER
T048  ENTER
T049  RCL K
T050  RCL- J
T051  x≥0?
T052  GTO T066
T053  SGN
T054  RCL+ K
T055  RCLx B
T056  RCL÷ K
T057  RCLx C
T058  STO C
T059  +
T060  x=y?
T061  GTO T066
T062  RCL T
T063  STO+ K
T064  R↓
T065  GTO T047
T066  R↓
T067  STO C
T068  RCL J
T069  RCL T
T070  RMDR
T071  x=0?
T072  GTO T090
T073  RCL- J
T074  x=0?
T075  STO C
T076  RCL A
T077  RCLx B
T078  RCLx C
T079  Pi
T080  ÷
T081  RCL A
T082  ATAN
T083  RCL T
T084  +/-
T085  SGN
T086  ACOS
T087  ÷
T088  +
T089  GTO T095
T090  RCL B
T091  √x
T092  RCLx A
T093  RCLx C
T094  RCL÷ T
T095  RCL T
T096  1/x
T097  +
T098  RCL X
T099  x²
T100  RCL÷ J
T101  RCL T
T102  SGN
T103  +
T104  LASTx
T105  RCL+ J
T106  RCL÷ T
T107  +/-
T108  y^x
T109  RCLx G
T110  x<>y
T111  STOP
T112  GTO T028

Here is Thomas' sample case:

Code:
[XEQ] T [ENTER]             J=?

Enter degrees of freedom
     7  [R/S]               X=?

Enter random variable X
    1,5 [R/S]               0,126263061   // =PDF(1,5|7)
                            0,911350757   // =CDF(1,5|7)
Do another calculation
        [R/S]               X=?
    -2  [R/S]               0,063135337   // =PDF(-2|7)
                            0,042809664   // =CDF(-2|7)

Final remark: like the other programs discussed in this thread, this is a rather simple and straightforward implementation that works fine for most practical cases. However, far out in the distribution tails the results may lose accuracy due to rounding, and eventually values may round to 0 or 1. For instance, for 7 d.o.f. and x=–100 you will get a CDF of 2,000...E–12 while the exact result is 1,318 E–12, and for x=0 the returned result is zero while it actually is 1,69 E–17. If you want more or less exact values even for such extreme arguments, a different approach is required.

Dieter
Dear Dieter, this is your t Student program. When I tryed with Thomas examples it works perfect, and the same when I try again with HP Stat Pack.


Attached File(s)
.pdf  HP 35s - Student´s t - Distribution_Dieter.pdf (Size: 61.87 KB / Downloads: 6)
Find all posts by this user
Quote this message in a reply
10-30-2015, 01:56 PM
Post: #37
RE: HP 35s Statistical Distributions Functions
(10-30-2015 02:36 AM)PedroLeiva Wrote:  Dear Dieter, this is your t Student program. When I tryed with Thomas examples it works perfect, and the same when I try again with HP Stat Pack.

There is room for improvement. I am working on an integrated solution that contains Normal, Chi² and Student's t distribution and calculates PDFs, CDFs and their inverses, the latter with a much improved initial guess. The F case is a bit more complicated, but I'll see.

BTW: quoting the whole previous post really is not required. If everyone would do so, the messages got longer and longer and longer as the thread continues and you would have to scroll down for a while until you find the new content. So when replying, you should simply delete everything that you do not directly refer to.

BTW2: do you have the 35s emulator? In this case I could send an .ep file that you can load into the emulator directly. This way you would not have to enter a few hundred lines of code.

BTW3: does someone know if there is a way to produce a listing or something similar from programs inside the 35s emulator?

Dieter
Find all posts by this user
Quote this message in a reply
10-30-2015, 05:20 PM
Post: #38
RE: HP 35s Statistical Distributions Functions
[/quote]
Dieter, you are absolutely right; from now on when I quote, I will erase the above text. Sorry, it would not hapend again.

I use of the Hewlett Packard HP 35s emulator, which enormously facilitates joining my work programs, and in turn I make far fewer errors. A remark about it, I wanted to send the files of the summarized previous programs in PDF, but the MoHPC system will not let me, so we should send them by email (leiva.pedro@inta.gob.ar).

It would be very useful to send a text file with instructions from the emulator program, but I do not know if that is really possible.
Find all posts by this user
Quote this message in a reply
10-30-2015, 07:16 PM
Post: #39
RE: HP 35s Statistical Distributions Functions
(10-30-2015 05:20 PM)PedroLeiva Wrote:  Dieter, you are absolutely right; from now on when I quote, I will erase the above text.

Of course not all text. Otherwise it's not a quote. ;-)

(10-30-2015 05:20 PM)PedroLeiva Wrote:  I use of the Hewlett Packard HP 35s emulator, which enormously facilitates joining my work programs, and in turn I make far fewer errors.

So you could use a .ep file and load it via the emulator menu (Calculator > Load).
I now have a collection with programs for the Normal, Student and Chi² distribution and their inverses. But there still is some work to do.

(10-30-2015 05:20 PM)PedroLeiva Wrote:  A remark about it, I wanted to send the files of the summarized previous programs in PDF, but the MoHPC system will not let me,

?!? – you have successfully attached several PDF files to your posts, e.g. the Chi² and Student programs. Where exactly is the problem?

Dieter
Find all posts by this user
Quote this message in a reply
10-30-2015, 07:43 PM (This post was last modified: 10-30-2015 07:45 PM by PedroLeiva.)
Post: #40
RE: HP 35s Statistical Distributions Functions
?!? – you have successfully attached several PDF files to your posts, e.g. the Chi² and Student programs. Where exactly is the problem?

Dieter
[/quote]
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)