(35S) Statistical Distributions Functions

[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