(25) Binomial Probability Distribution

[Dave Britten]

This program calculates the cumulative binomial probability distribution between a given lower and upper value for r.

B(n, r, P) = Binomial probability mass of r successes in n independent trials, each with a chance of success P.

The program is built around these two formulas:

ln(B(n, 0, P)) = n * ln(1-P)

ln(B(n, r+1, P)) = ln(B(n, r, P)) + ln(n-r) - ln(r+1) + ln(P) - ln(1-P)

Logarithms are used to avoid over/underflows in intermediate calculations.

Usage:

Enter the program, and store your four inputs, n, P, MIN, and MAX, into the designated registers.

R1 = n
R2 = P
R3 = MIN (i.e. lower value of r)
R4 = MAX (i.e. upper value of r)

Press f PRGM, R/S, and wait for the cumulative probability to be returned.

A lower tail can be obtained by setting R3/MIN to 0, and an upper tail by setting R4/MAX equal to n. Set R3=R4 to compute probability mass for a single number of successes.

Program code:

Code:

01    0
02    STO 0
03    STO 5
04    RCL 1
05    1
06    RCL 2
07    -
08    ln
09    *
10    STO 7
11    LASTx
12    CHS
13    RCL 2
14    ln
15    +
16    STO 6
17    RCL 3
18    RCL 0
19    x<y
20    GTO 25
21    RCL 7
22    e^x
23    STO+ 5
24    RDown
25    RCL 4
26    x=y
27    GTO 44
28    RDown
29    CHS
30    RCL 1
31    +
32    ln
33    RCL 0
34    1
35    +
36    ln
37    -
38    RCL 6
39    +
40    STO+ 7
41    1
42    STO+ 0
43    GTO 17
44    RCL 5
45    GTO 00

Example:

100 trials, P = 40%, what is the probability that 0-50 trials succeed?

100 STO 1
.4 STO 2
0 STO 3
50 STO 4
f PRGM R/S

...wait approx. 2 minutes...

0.9832 (fix 4), or 98.32%

[Gene]

Here's the version from the PPC Journal HP-25 Library of programs.

#41

[Dave Britten]

(11-17-2019 09:48 PM)Gene Wrote:  Here's the version from the PPC Journal HP-25 Library of programs.

#41

That one is conceptually similar, but doesn't use logs (so p^n could quickly underflow if n is very large). Also, it only does a lower tail, rather than a cumulative range between any two min and max values. So it's quicker to type in, and probably runs a little faster, but it's a bit more limited. Use mine if you need to do an extremely large number of trials. It'll take a good amount of time to run, but it's still faster than submitting your job to the computing center and picking up the results in the morning, right?

That PPC Journal HP 25 Library is a real treasure trove. The prime factor finder (#33, I think it was) is extremely clever and efficient, even without enough memory/features to do a full mod 30 sieve.

[PedroLeiva]

Another example
EXAMPLE2:
20 trials, P= 20%, what is the probability that 0-6 trials succeed? and 1 succeed?
n= 20
P= 0,2
r MIN= 0
r MAX= 6
P(0 to 6 successes) = .0115 + .0576 + .1369 + .2054 + .2182 + .1746 + .1091 = .9133
==> so cumulative probability= 0,9133 or 91.33%
P(0 success) = .0115
==> P(1 success) = .0576 or 5.76%
also…..
P(6 successes) = .1091
(TYVM to Albert Chan for the corrections !!)

I try with RPN GO-25 and it only takes a few seconds

Pedro

[Gene]

(11-17-2019 09:57 PM)Dave Britten Wrote:  That PPC Journal HP 25 Library is a real treasure trove. The prime factor finder (#33, I think it was) is extremely clever and efficient, even without enough memory/features to do a full mod 30 sieve.

Gene: Agreed. Thanks for your version!

[Albert Chan]

It might also be of interest: An approximation of the Cumulative Binomial Probability Distribution
Note: if p > 0.5, then B(x; n, p) = 1 - B(n+1-x; n, 1-p)

(11-17-2019 09:14 PM)Dave Britten Wrote:  100 trials, P = 40%, what is the probability that 0-50 trials succeed?

For above example, x=50, n=100, p=0.4:

z = (x+½ - n*p) / √(n(p-p²)) = (50.5 - 40) / √24 ≈ 2.1433
k = (1-2p) / (6 √(n(p-p²))) ≈ 0.0068041
a = (1/(p-p²) - 6)/(24n) * (z³ - z) ≈ -0.0058838
b = k²/2 * (z^5 - 10 z^3 + 15 z) ≈ -0.00048795

B(x; n, p) ≈ CDF(z) - (k(z²-1) + a + b) PDF(z) ≈ 0.98396 - 0.01808 * 0.040123 ≈ 0.98323

Update: Camp-Paulson normal approximation is simpler and more accurate for small n:
\(\large a = {1 \over 9(n-x)} \quad,\quad b = {1\over 9(x+1)} \quad,\quad q = \sqrt[3]{{1-p \over p}{a \over b} } \)

\(\large B(x; n,p) ≈ CDF( z = { (1-b)q\,-\, (1-a) \over \sqrt{b\,q^2 + a}} ) \)

→ B(50; 100, 0.4) ≈ CDF(2.12540) ≈ 0.9832
→ B( 6 ; 20 , 0.2 ) ≈ CDF(1.36167) ≈ 0.9133