(11C) Poisson distribution

[Gamo]

The Poisson distribution is popular for modelling the number of times an event occurs in an interval of time or space.

Formula: P(k events in interval) = (e^-λ)(λ^k) / k!

where:
λ (lambda) is the average number of events per interval
e is the number 2.71828... (Euler's number) the base of the natural logarithms
k takes values 0, 1, 2, …
k! = k × (k − 1) × (k − 2) × … × 2 × 1 is the factorial of k.

Example Problem:
Ugarte and colleagues report that the average number of goals in a World Cup soccer match is approximately 2.5
Because the average event rate is 2.5 goals per match, λ = 2.5
What is the probability of gold of P(k) = 0, 1, 2, 3, 4, 5, 6, 7

Program:

Code:


LBL A (λ)
STO 1
RTN
LBL B (k)
STO 2
RTN
LBL C (P)
RCL 1
CHS
e^x
RCL 1
RCL 2
Y^x
x
RCL 2
X!
/
RTN

Run Program:
2.5 A
0 B
C 0.082
1 B
C 0.205
2 B
C 0.257
3 B
C 0.213
.
.
.
.
7 B
C 0.010

The table below gives the probability for 0 to 7 goals in a match.

k P(k goals in a World Cup soccer match)
0 0.082
1 0.205
2 0.257
3 0.213
4 0.133
5 0.067
6 0.028
7 0.010

Credit to Wikipedia for information and example problem.

Gamo

[Dieter]

(12-18-2017 09:38 AM)Gamo Wrote:  What is the chance of gold of P(k) = 0, 1, 2, 3, 4, 5, 6, 7

Gold ?-)  I assume this is supposed to mean
"What is the probability P(k) for k = 0, 1, 2, 3, 4, 5, 6 or 7 goals".

But why do you use two separate labels for k and P(k)? This way calculating the PDF always requires pressing two keys, B and C.

Here is another version that also calculates the CDF, i.e. P(k₁ ≤ k ≤ k₂).

Code:

LBL A
STO 0
RTN
LBL B
RCL 0
x<>y
y^x
LastX
x!
/
RCL 0
CHS
e^x
*
RTN
LBL C
ABS
INT
EEX
3
/
x<>y
ABS
INT
+
STO I
LastX
GSB B
x=0?
RTN
ENTER
ENTER
LBL 1
ISG      // (ISG I on the 15C)
GTO 2
x<>y
RTN
LBL 2
RCL 0
*
RCL I
INT
/
+
LastX
GTO 1

Example for λ = 2,5:

Code:

Enter λ

 2,5        [A]    2,5000

Calculate the probability for 1, 2 or 3 goals

  1         [B]    0,2052   P(1)
  2         [B]    0,2565   P(2)
  3         [B]    0,2138   P(3)

Calculate the probability of a match with 1 to 4 goals

1 [ENTER] 4 [C]    0,8091   P(1 ≤ k ≤ 4)

If desired: [R↓]   0,1336   P(4)
            [R↓]   0,2052   P(1)

Direct evaluation of the Poisson PDF often leads to overflow errors. Even cases where k>69 can not be handled this way. Too bad there is no lnΓ function available, this could provide an easy fix.

But there are two workarounds:

1. The recursive method of the CDF routine significantly extends the useable range for λ and k, and this can also be used for calculating the PDF:
Simply enter 0 [ENTER] k [C], and when the result is displayed press [R↓] or [x<>y] to get P(k).

Example:
Evaluate P(80) for λ=90.

Code:

   90        [A]  90,0000
   80        [B]   8,1940 -40  flashing
                   Overflow error while trying to evaluate 90^80 and 80!

0 [ENTER] 80 [C]   0,1582
            [x<>y] 0,0250   P(80)

However, the iteration with k required loops may take some time on a hardware 11C.

Also please note that for λ > 227,9559242 the expression e^–λ will underflow to zero. In this case 0 is returned.

2. For large λ and k the following code may be used. It implements a Stirling-based approximation, and since there is no iteration the result is returned immediately.

Code:

LBL E
STO I
RCL 0
RCL I
/
LN
1
+
*
RCL 0
-
e^x
RCL I
6
1/x
+
Pi
*
2
*
sqrt
/
RTN

Example:
Evaluate P(80) for λ=90.

Code:

   90        [A]  90,0000
   80        [E]   0,0250   P(80)

Dieter

[Gamo]

Thanks Dieter

Very nice program detail with in-depth information.

Gamo

[SlideRule]

Gamo

Check the HP-25 post for Poisson Distribution from 1977.

BEST!
SlideRule

[Dieter]

(12-23-2017 06:42 PM)SlideRule Wrote:  Check the HP-25 post for Poisson Distribution from 1977.

That's three short programs in this thread.
I just posted an improved version. Which also includes a correct initialization of the summation registers to avoid erroneous results.
And there's even some information on the author of the original programs.

Dieter