(35S) Statistical Distributions Functions

[Dieter]

(11-19-2015 09:24 AM)Gerald H Wrote:  Could you, Dieter, please include a pdf file of the latest programmes in human readable form to enter on the calculator?

The current programs have been developed on HP's 35s emulator so that they can be run by simply loading the posted *.ep file into the emulator.

These programs have between 900 and 1000 lines of code. At the moment I do not know of a way to get them out of the emulator, i.e. all I could do is copy every single line by hand. If someone knows a way to transform an .ep emulator file into a readable listing, please let me know.

On the other hand you anyway would have to enter the programs into your 35s manually, line by line. Does it really make much of a difference if this is done with a PDF listing or by stepping through the programs on the emulator?

Dieter

[PedroLeiva]

[/quote]
I'm working on a list of programming keystroke steps based on the .ep file. There are 10 labels and 926 steps, a total of 11 pages 4A landscape. Once finished, I will put these at the disposal of the Forum with the authorization of Dieter, after exhaustive revision errors. Pedro

[Gerald H]

I'm amazed, Dieter, that you have developed such a large programme without having a human readable source to refer to - Congratulations!

I look forward to seeing your, PedroLeiva's, programme listing.

[Dieter]

(11-19-2015 01:17 PM)PedroLeiva Wrote:  I'm working on a list of programming keystroke steps based on the .ep file. There are 10 labels and 926 steps, a total of 11 pages 4A landscape.

I'd recommend you don't start before the final revision. ;-)

This morning the programs had 938 steps – I had added this and that, but this also allowed an error to creep in (which now has been corrected). Maybe I'll post another version in the next days.

At the moment I am thinking of adding the Binomial and Poisson distribution functions. This would require just a few lines of code (less than 50 for both) since they can be calculated by means of the regularized Beta and Gamma functions. On the other hand maybe these two should better go into a second program pack for discrete distributions. Also, implementing the respective quantile functions here are a bit trickier.

(11-19-2015 01:17 PM)PedroLeiva Wrote:  Once finished, I will put these at the disposal of the Forum with the authorization of Dieter, after exhaustive revision errors. Pedro

Yes, extensive testing is required. I suggest you use something like Excel or OpenOffice Calc for the listings. This allows simple editing as well as some other useful features.

Dieter

[PedroLeiva]

Excellent idea. The inclusion of discrete variables complement the main statistical functions. In that sense I attach the programs I wrote for the HP 35s, including the listing of programming instructions and numerical examples. I ask pardon because the instructions are written in my native tongue, Spanish.
I learned to compress files, lets see if I did it well! Pedro

[Marcio]

Excellent work!
Good thing you avoided the EOS so the codes can be easily ported into other machines.

Thanks

Marcio

[Dieter]

(11-20-2015 12:09 PM)PedroLeiva Wrote:  Excellent idea. The inclusion of discrete variables complement the main statistical functions.

Here is a short and not at all optimized listing for the Poisson and Binomial cdfs. Simply add it to the current code. It's less than 50 steps.

Code:

P001 LBL P    // Poisson CDF
P002 XEQ J012 // initialize (O=1, T=2 etc.)
P003 INPUT L  // prompt for parameter lambda
P004 x<=0?
P005 GTO P003 // accept only positive lambda
P006 INPUT X
P007 x<0?
P008 GTO P006 // reject negative x
P009 RCL+ O
P010 STO A
P011 RCL L
P012 STO Y
P013 XEQ G006  // Call regularized Gamma. Check: this points to RCL A RCL+T...
P014 x<>y      // rearrange output so that P(x,lambda) is in X
P015 STOP
P016 GTO P006

O001 LBL 0     // Binomial CDF - LBL B is already in use ;-)
O002 -14       // prompt for n and initialize
O003 XEQ J001
O004 INPUT P
O005 RCL O
O006 RCL P
O007 x<y?
O008 x<=0?
O009 GTO O004  // reject if p<=0 or p>=1
O010 XEQ A001
O011 x<0?
O012 GTO O010  // reject negative x
O013 RCL N
O014 RCL X
O015 x<y?
O016 GTO O020
O017 CLSTK     // if x>=n
O018 RCL O     // return cdf = 1
O019 GTO O029
O020 RCL+ O
O021 STO A
O022 RCL N
O023 RCL- X
O024 STO B
O025 RCL P
O026 STO Y
O027 XEQ I016   // Call regularized Beta. Check: this points to RCL A  RCL B  XEQ B010
O028 x<>y       // rearrange output as acutally 1-cdf was calculated
O029 STOP
O030 GTO O010

If this does not work you may have an older version of the other distributions. Especially if in some cases the correct Poisson result apprears in Y instead of X. ;-)

Test examples:
Poisson cdf for lambda=7,3 and x=5 is 0,26404.
Binomial cdf for n=20, p=0,2 and x=6 is 0,91331.

(11-20-2015 12:09 PM)PedroLeiva Wrote:  In that sense I attach the programs I wrote for the HP 35s, including the listing of programming instructions and numerical examples. I ask pardon because the instructions are written in my native tongue, Spanish.
I learned to compress files, lets see if I did it well! Pedro

As far as I can tell these programs calculate the density functions (probability mass functions), but not the cumulative densities.

Two questions: First, why is it required to clear all variables (CLEAR 2), i.e. have the user lose all his data, before starting the programs? It's good practice that the program itselfs clears only the registers it requires. At the moment I even do not see a reason why any register has to be cleared at all. ?!?

And second: why does the Binomial distribution calculate the combinations nCx although there is a dedicated HP35s nCr function that does this faster and with less problems (e.g. overflow)? You have to know that this function exists since you used it for the Hypergeometric distribution. ;-)

BTW, it seems that the program titled "Distribución Polinomial" actually calculates the multinomial distribution (es: distribución multinomial). Since I am not sure about the correct Spanish term: is there a difference between these two expressions or do they refer to the same distribution? If yes: the multinomial distribution can be calculated much easier and without using so many registers.

Dieter

[PedroLeiva]

Dieter, you are absolutely right, my programs calculate the probability for a variable value, not the cumulative probability. The instruction to delete all variables not is correct, I already have erased. Regarding the use of combinatorial calculation function, I did not use it because I decided to test my knowledge of statistics.

I loaded your labels for cumulative Binomial and Poisson function, they work perfectly well.
I wonder why HP says 35s has 800 programation steps, whereas now with 982 still working well ?

[Dieter]

(11-20-2015 01:19 PM)Dieter Wrote:  BTW, it seems that the program titled "Distribución Polinomial" actually calculates the multinomial distribution (es: distribución multinomial). Since I am not sure about the correct Spanish term: is there a difference between these two expressions or do they refer to the same distribution? If yes: the multinomial distribution can be calculated much easier and without using so many registers.

OK, here is a version that features some advantages:

• Usage is quite straightforward. No need to prestore values in registers or calculate sums manually.
• Uses merely two registers instead of eleven.
• Number k of events can be anything between 1 and ...infinity
• Erroneous data entry can be corrected by re-entering with a negative frequency x.
• Despite all this: shorter. Could be even shorter (not much more than 20 lines) without the error recovery function (see below)
  

Here is the listing:

Code:

M001  LBL M
M002  CLSTK
M003  STO N     n=0
M004  e^x
M005  STO P     p=1
M006  ALL
M007  CLSTK
M008  VIEW N    show accumulated n
M009  x=0?      zero or no entry?
M010  GTO M027  jump to probability calculation
M011  x<>y
M012  STO+ N    n := n + x_i
M013  SGN
M014  x<>y      save sign(x_i)
M015  LASTx
M016  ABS
M017  y^x
M018  LASTx
M019  !
M020  ÷         p_i^x_i / x_i!
M021  x<>y      if xi was negative this is -1 (instead of 1 else)
M022  y^x       which turns this into its reciprocal
M023  STOx P    update product
M024  GTO M006
M025  RCL N
M026  !
M027  STOx P    calculate probability
M028  FIX 6
M029  RCL P     return result in X
M030  VIEW P    and display labelled output
M031  GTO M001  start a new calculation with R/S

Here is an example using the data from your (Pedro's) PDF:

Code:

 x1 = 2     p1 = 0,25
 x2 = 2     p2 = 0,25
 x3 = 5     p3 = 0,25
 x4 = 3     p4 = 0,25


 [XEQ] M  [ENTER]
                       N=         0    // Data is entered by
 2 [ENTER] 0,25 [R/S]                  // x_i [ENTER] p_i  [R/S]
                       N=         2
 2 [ENTER] 0,25 [R/S]
                       N=         4
 5 [ENTER] 0,5  [R/S]
                       N=         9    // oops... wrong data
-5 [ENTER] 0,5  [R/S]                  // re-enter with negative x
                       N=         4    // and it gets removed
 5 [ENTER] 0,25 [R/S]
                       N=         9
 3 [ENTER] 0,25 [R/S]
                       N=        12
                [R/S]                  // data complete, simply press [R/S] or 0 [R/S]
                       P=  0,009913    // result

Hmmm... I just got another idea. Maybe this could be done with the Σ+ key using a sum of logs instead of a product... *8)

Edit: Yes, this will also work:

Code:

M001  LBL M
M002  CLΣ
M003  ALL
M004  Σx
M005  STO N
M006  CLSTK
M007  VIEW N
M008  x=0?
M009  GTO M030
M010  LN
M011  x<>y
M012  SGN
M013  STO N     save sign(x_i)
M014  LASTx     * x_i
M015  x         = abs(x_i)
M016  x
M017  LASTx
M018  ENTER
M019  R↓
M020  !
M021  LN
M022  -
M023  R↑
M024  RCLx N    multiply
M025  x<>y      both values
M026  RCLx N    by +1 or -1
M027  x<>y      so that this
M028  Σ+        Σ+ may become a Σ-
M029  GTO M003
M030  Σx
M031  !
M032  LN
M033  Σy
M034  +
M035  e^x
M036  STO P
M037  FIX 6
M038  VIEW P
M039  GTO M001

Program usage is the same as before. The code uses only logs of the factorial function (both "!" are followed by a "LN") so that this approach may be especially useful for calculators with lnΓ function, thus avoiding potential overflow problems.

Dieter

[Dieter]

(11-20-2015 05:45 PM)PedroLeiva Wrote:  I loaded your labels for cumulative Binomial and Poisson function, they work perfectly well.

Just to be sure: try the Poisson example with lambda=7,3 for x=0, x=7 and x=15.
You should get 0,0006755..., 0,5541... and 0,99641...

(11-20-2015 05:45 PM)PedroLeiva Wrote:  I wonder why HP says 35s has 800 programation steps, whereas now with 982 still working well ?

HP says the 35s can have up to 800+ data registers. These are the ones that are accessed indirectly by positive values in I or J. They are dynamically allocated (and de-allocated by the CLVARx command). This memory area is very handy if you want to handle large data sets, e.g. a 20x20 matrix.

As far as programs are concerned: Every labelled program can have up to 999 steps (e.g. A001...A999), which would allow for at least 26 (A...Z) x 999 = almost 26.000 steps – if there was sufficient memory, that is. So there is quite some room for more code. The whole 35s has 30 KB of memory, shared by programs, data, equations and some internal functions (e.g. Solve). Check MEM and see what's left. ;-)

Dieter

[Marcio]

(11-20-2015 08:30 PM)Dieter Wrote:  As far as programs are concerned: Every labelled program can have up to 999 steps (e.g. A001...A999), which allows for at least 26 (A...Z) x 999 = almost 26.000 steps. So there is sufficient room for more code. ;-)

Dieter

Is there any way to assess this?

I remember that I had like 1200+ lines of codes and they occupied nearly 1/3 of the device's memory, no equations anywhere.

Marcio

[Dieter]

(11-20-2015 08:37 PM)Marcio Wrote:  

(11-20-2015 08:30 PM)Dieter Wrote:  As far as programs are concerned: Every labelled program can have up to 999 steps (e.g. A001...A999), which allows for at least 26 (A...Z) x 999 = almost 26.000 steps. So there is sufficient room for more code. ;-)

Is there any way to assess this?

You can use the MEM command to see the remaining memory size in Bytes.

N.B.: I have edited my previous post before I read your reply. The 26000 steps are a theoretical limit which would only be reached if there was sufficient memory. Since every program step requires several bytes, the max. number of steps should be substantially below 10000 steps.

(11-20-2015 08:37 PM)Marcio Wrote:  I remember that I had like 1200+ lines of codes and they occupied nearly 1/3 of the device's memory, no equations anywhere.

Every single program step requires at least three bytes, some use (even substantially) more. So 1200 lines times, say on average 4 Bytes/line are already nearly 5000 Bytes or 1/6 the full capacity. Especially using numbers (even a simple "1" or "15") is very memory consuming (and slow, for that matter), so that it may well be possible that a 1200 line program can use 30% of the total memory.

I remember the days when the HP41 was new and how I tried to squeeze off one or two more bytes to make the program fit into the 224 byte capacity of a single magnetic card. ;-)

Dieter

[Dieter]

(11-20-2015 05:45 PM)PedroLeiva Wrote:  I loaded your labels for cumulative Binomial and Poisson function, they work perfectly well.

I have reworked this a bit so that now the same values as with the other distributions are returned. This includes the upper CDF (which in these cases is not just 1 minus the lower CDF) as well as the PDF (or PMF - probability mass function -, to be more precise).

Code:

P001 LBL P    // Poisson CDF
P002 XEQ J012 // initialize (O=1, T=2 etc.)
P003 INPUT L  // prompt for parameter lambda
P004 x<=0?
P005 GTO P003 // accept only positive lambda
P006 INPUT X
P007 x<0?
P008 GTO P006 // reject negative x
P009 RCL+ O
P010 STO A
P011 RCL L
P012 STO Y
P013 XEQ G006  // Call regularized Gamma. Check: this points to RCL A RCL+T...
P014 RCL L     // calculate PDF/PMF
P015 RCL X
P016 y^x
P017 RCL X
P018 !
P019 ÷
P020 RCL L
P021 +/-
P022 e^x
P023 x
P024 STO C
P025 +         // upper CDF = 1 - lower CDF + PMF
P026 x<> C
P027 STO D
P028 x<>y
P029 RCL C
P030 x<>y
P031 ENTER
P032 CLx
P033 STO A
P034 STO B
P035 R↓
P036 STOP
P037 GTO P006

O001 LBL 0     // Binomial CDF - LBL B is already in use ;-)
O002 -14       // prompt for n and initialize
O003 XEQ J001
O004 INPUT P
O005 RCL O
O006 RCL P
O007 x<y?
O008 x<=0?
O009 GTO O004  // reject if p<=0 or p>=1
O010 INPUT X
O011 x<0?
O012 GTO O010  // reject negative x
O013 RCL N
O014 RCL X
O015 x<y?
O016 GTO O023  // regular case: x < n
O017 RCL N     // handle x=n separately
O018 STO X
O019 RCL O     // return cdf = 1
O020 ENTER
O021 CLx       // and 1-cdf = 0
O022 GTO O032
O023 RCL+ O
O024 STO A
O025 RCL N
O026 RCL- X
O027 STO B
O028 RCL P
O029 STO Y
O030 XEQ I016   // Call regularized Beta. Check: this points to RCL A  RCL B  XEQ B010
O031 STO C
O032 R↓
O033 RCL N      // calculate PDF/PMF
O034 RCL X
O035 nCr
O036 RCL P
O037 RCL X
O038 y^x
O039 x
O040 RCL O
O041 RCL- P
O042 RCL N
O043 RCL- X
O044 y^x
O045 x
O046 STO D
O047 STO+ C    // upper CDF = 1 - lower CDF + PMF
O048 x<>y
O049 RCL C
O050 x<>y
O051 ENTER
O052 CLx
O053 STO A
O054 STO B
O055 R↓
O056 STOP
O057 GTO O010

These versions return:

in Z: the Poisson resp. Binomial PDF, i.e. P(X=x)
in Y: the Poisson resp. Binomial upper CDF, i.e. P(X≥x)
in X: the Poisson resp. Binomial lower CDF, i.e. P(X≤x)

This matches the returned values for the Chi² and F distribution. Like in these cases, the upper CDF and the PDF are returned in registers C and D as well.

Looks like the user instructions need an update. ;-)

At the moment I don't think I'll add the Poisson and Binomial quantile functions – that's a rather complex story.

Edit: it looks like the Poisson quantile is not that complicated. At least an experimental version seems to run without major problems. ;-)

Dieter

[Dieter]

(11-20-2015 11:06 PM)Dieter Wrote:  I have reworked this a bit so that now the same values as with the other distributions are returned. This includes the upper CDF (which in these cases is not just 1 minus the lower CDF) as well as the PDF (or PMF - probability mass function -, to be more precise).
(...)
Edit: it looks like the Poisson quantile is not that complicated. At least an experimental version seems to run without major problems. ;-)

OK, here is a first version of the complete Poisson distribution: upper and lower CDF, PMF and – the inverse cdf (quantile). Usage is the same as with other programs of this pack. Only the quantile function returns something different:

As the Poisson distribution handles a discrete random variable (x=0, 1, 2, 3, ...) the returned quantile is an integer. However, in most cases there is no exact integer for a given cumulative probability. That's why usually the quantile is defined as the smallest integer with a cdf greater than or equal to the given probability. For instance, if p=0,9 while cdf(11)=0,85 and cdf(12)=0,93 the quantile is 12.

This may lead to off-by-one errors in very close cases where the calculated cdf may be not exact in the last digits. That's why the following program returns three values. In x the calculated quantile is returned, and y at the same time holds the corresponding probability (cdf). Scrolling down the stack with R↓ displays the cdf for the next lower integer, i.e. for x–1. This way you can see the cdfs for the two adjacent integers that bracket the given probability.

Example: mean = 7,3 and p=0,9. The program returns x=11 as well as 0,9319 = cdf(11) and 0,8788 = cdf(10).

The method is quite simple and straightforward. It implements some insights of a 2013 paper by Michael Short titled "Improved Inequalities for the Poisson and Binomial Distribution and Upper Tail Quantile Functions". Actually the Poisson quantile can (almost) be calculated directly by means of Lambert's W function, but since the 35s does not offer one (unlike the 34s) I chose a different approach.

Finally, here is the code:

Code:

P001  LBL P     Poisson distribution starts here
P002  XEQ J012  initialize
P003  INPUT M   prompt for mean
P004  X<=0?
P005  GTO P003  reject non-positive entry
P006  XEQ A001  prompt for x resp. p
P007  FS? 1     quantile mode?
P008  GTO P040  go to quantile function
P009  X<0?
P010  GTO P006  reject negative x
P011  RCL+ O
P012  STO A     a = x+1
P013  RCL M
P014  STO Y     y = mean
P015  XEQ G006  call CDF = incomplete Gamma(x+1, mean)
P016  XEQ P031  call PMF(x)
P017  STO C
P018  +         add PMF to (1-CDF) to get upper CDF
P019  x<> C     store in C
P020  STO D     and move PMF to D
P021  x<>y
P022  RCL C     arrange data for output
P023  x<>y
P024  ENTER
P025  CLx
P026  STO A     clear A and B
P027  STO B
P028  R↓        and T as well
P029  STOP      return T: 0  Z: PMF  Y: upper CDF  X: lower CDF
P030  GTO P006
P031  RCL M     calculate PMF
P032  LN        use logs to avoid overflow
P033  RCLx X
P034  RCL- M
P035  e^x
P036  RCL X
P037  !         ln Gamma would be nice here
P038  ÷         = PMF
P039  RTN
P040  RCL M     start of quantile function
P041  +/-
P042  e^x
P043  RCL P
P044  x>y?      p > CDF(0) ?
P045  GTO P050  then continue
P046  CLx
P047  STO X     else return 0
P048  x<>y
P049  GTO P104  jump to exit
P050  XEQ E001  get Normal quantile estimate
P051  RCLx E
P052  3 E-4     make sure Normal quantile
P053  +         is slightly high
P054  ENTER
P055  x²
P056  6
P057  ÷
P058  x<>y
P059  RCL M
P060  √x
P061  x
P062  +
P063  RCL+ M
P064  IP
P065  x<=0?
P066  RCL T
P067  RCL+ O
P068  STO X     first guess, may be slightly high
P069  RCL+ O
P070  STO A
P071  RCL M
P072  STO Y
P073  XEQ G006  calculate CDF
P074  x>y?      CDF < 0,5?
P075  GTO P080
P076  +/-       if not, p is closer to 1
P077  RCL P     so better compare -(1-CDF)
P078  RCL- O    and -(1-p)
P079  GTO P082
P080  R↓
P081  RCL P
P082  x<>y
P083  STO W
P084  XEQ P031  calculate PMF(x)
P085  -         CDF(x) - PMF(x) = CDF(x-1)
P086  x<y?      is this less than p ?
P087  GTO P097  then exit
P088  STO W     else save CDF(x)
P089  LASTx
P090  RCLx X    calculate PMF(x) from PMF(x+1)
P091  RCL÷ M
P092  DSE X     decrement x
P093  GTO P085  and continue
P094  -         if ever executed, this should round to zero
P095  CLx       so make it exactly zero ;-)
P096  ENTER
P097  RCL W     now x and y hold the previous and current CDF
P098  x>0?      are we working with the CDFs (not with 1-CDF) ?
P099  GTO P104
P100  RCL+ O    if not, restore CDF = 1+[-(1-CDF)]
P101  x<>y
P102  RCL+ O    do so both for CDF(x) and CDF(x-1)
P103  x<>y
P104  RCL X
P105  ENTER
P106  CLx
P107  R↓        clear T
P108  STOP      return T: 0  Z: CDF(x-1)  Y:CDF(x)  X: quantile x
P109  GTO P006

Try it and see what you get. As usual, this is supposed to work with the functions J, A, E and G included in the distributions package that was posted as an .ep file for the 35s emulator. I did not do much testing, so all error reports are welcome.

Edit (2015-11-23): added line P065/P066. However, there are still a few cases with larger mean an very small p that may return a slightly low result, e.g. m=50 and p=1E-17 returns 2 instead of 4.

Edit2 (2015-11-26): this problem can be avoided by a different approach which also seems to work for the Binomial quantile.

Dieter

[PedroLeiva]

For the new example, m=50 and p=1E-17 (SF 1), I don´t get the same result as you. For me is X=0, Y=1, Z=1, T=0. What is wrong? Pedro

[Dieter]

(11-26-2015 11:07 PM)PedroLeiva Wrote:  For the new example, m=50 and p=1E-17 (SF 1), I don´t get the same result as you. For me is X=0, Y=1, Z=1, T=0. What is wrong? Pedro

Make sure you added lines P065 and P066. And check if XEQ G006 points to the correct line (RCL A). But most important, there was an error in P049 which, after adding these lines, should read GTO P104. I just corrected this.

Or, even better, forget about all this, wait a few days and use the complete recent version I just finished – including the Binomial quantile. It still requires some testing but I think I will post the .ep file this weekend.

Dieter

[Dieter]

(11-27-2015 06:41 AM)Dieter Wrote:  Or, even better, forget about all this, wait a few days and use the complete recent version I just finished – including the Binomial quantile.

Here it is. The attached file contains version 5 of the distributions package as well as a PDF with instructions for the Normal, Student, Chi² and Fisher distribution. Also included are instructions for some additional functions, while the part that handles the discrete distributions (Binomial, XEQ O and Poisson, XEQ P) is still missing.

The Binomial and Poisson quantile functions use the 35s solver, so you'll see RUNNING...SOLVING....RUNNING before the result is returned. Here X holds the quantile, Y has the corresponding CDF and Z holds the CDF of x–1. This way you can check the result and see how close a continuous (real) quantile is to the adjacent integers. Note that the lower and uppers CDFs for the discrete discrtibutions do not sum up to 1 because they are P(X≤x) and P(X≥x), which means that PMF(x) is part of both sums.

All this requires extensive testing, so load the .ep file into your 35s emulator and see what you get. There may and will be errors, so I'll appreciate all your error reports. Thank you in advance for all corrections and suggestions.

Dieter

  Distributions_35s_5.zip (Size: 156.5 KB / Downloads: 57)