(HP-65) Binomial distribution with large input support

[Dave Britten]

This program calculates the binomial probability mass function, or cumulative binomial probabilities, with support for large (>69) numbers of trials. To do this, it uses Nemes' formula to approximate ln(x!) rather than calculating x! directly and overflowing.

Recommended Key Labels

A: CALC
E: ln(x!)

Usage

There weren't enough steps left for a fancy user interface, so inputs must first be stored directly into R1-R4 before running the program. Store these values:

R1: n - number of trials
R2: p - percentage of success on a single trial
R3: xmin - minimum number of successes
R4: xmax - maximum number of successes

Then press A to calculate the probability mass/cumulative probability. If xmin=xmax, the program calculates the probability mass for that value of x. If xmax>xmin, then the program calculates cumulative probability for the specified range of x values. If xmin>xmax, or xmax>n, you'll get an error.

Note that Nemes' formula is an approximation, which means this program produces approximate results. Accuracy improves with larger inputs, but the worst-case scenario appears to be about 4 digits of accuracy (e.g. n=2, p=0.5, xmin=0, xmax=2 yields about 100.037%). So it's no worse than a slide rule!

Probability mass is calculated in constant time - about 10 seconds - and a cumulative probability with xmax-xmin=10 takes about 30 seconds. Execution time in seconds can be estimated using 10 + 2(xmax - xmin).

Also, you can press E to calculate ln(x!).

Example

If you roll a fair six-sided die 4000 times, what are the odds that you will get any given result between 662 and 672 times? (The mean of this distribution is 666.67.)

4000 STO 1
6 g 1/x STO 2
662 STO 3
672 STO 4
A: 0.1844

There is approximately an 18.44% chance that any specific result will be rolled between 662 and 672 times.

Code:

RCL 1   34 01
E       15
RCL 3   34 03
STO 8   33 08
E       15
-       51
RCL 1   34 01
RCL 3   34 03
-       51
E       15
-       51
RCL 3   34 03
RCL 2   34 02
f       31
LN      07
STO 5   33 05
*       71
+       61
RCL 1   34 01
RCL 3   34 03
-       51
1       01
RCL 2   34 02
-       51
f       31
LN      07
STO 6   33 06
*       71
+       61
f-1     32
LN      07
STO 7   33 07
LSTx    35 00
LBL     23
1       01
RCL 8   34 08
RCL 4   34 04
g x=y   35 23
RCL 7   34 07
RTN     24
g RDown 35 08
RCL 1   34 01
g x><y  35 07
-       51
RCL 8   34 08
1       01
+       61
STO 8   33 08
/       81
f       31
LN      07
+       61
RCL 5   34 05
+       61
RCL 6   34 06
-       51
f-1     32
LN      07
STO     33
+       61
7       07
g LSTx  35 00
GTO     22
1       01
LBL     23
E       15
1       01
+       61
STO 7   33 07
1       01
2       02
*       71
.       83
1       01
RCL 7   34 07
/       81
-       51
g       35
1/x     04
RCL 7   34 07
+       61
f       31
LN      07
1       01
-       51
RCL 7   34 07
*       71
2       02
g       35
PI      02
*       71
RCL 7   34 07
/       81
f       31
LN      07
2       02
/       81
+       61
RTN     24

[Thomas Klemm]

As a rule of thumb the binomial distribution can be approximated by the normal distribution if both \(np > 5\) and \(nq > 5\).
This is given in the example since:

\(
\begin{align}
4000 \cdot \frac{1}{6} &\approx 666.6667 > 5 \\
4000 \cdot \frac{5}{6} &\approx 3333.3333 > 5 \\
\end{align}
\)

We can calculate mean \(\mu\) and standard deviation \(\sigma\) with:

\(
\begin{align}
\mu &= n \cdot p \\
\sigma &= \sqrt{n \cdot p \cdot q} \\
&= \sqrt{n \cdot p \cdot (1 - p)} \\
\end{align}
\)

4000
ENTER
6
1/x
×
STO 0

666.66667

1
LSTx
-
×
\(\sqrt{x}\)

STO 1

23.570226

The following program for the HP-15C can be used to integrate the pdf of the normal distribution:

Code:

   001 {    42 21 11 } f LBL A
   002 {       43 11 } g x²
   003 {           2 } 2
   004 {          10 } ÷
   005 {          16 } CHS
   006 {          12 } eˣ

In order to get the best approximation we subtract 0.5 from \(x\) or add 0.5 to \(x\).

FIX 5

661.5
RCL - 0
RCL ÷ 1

-0.21920

672.5
RCL - 0
RCL ÷ 1

0.24749

\(\int_y^x\) A

0.46244


Now we just have to compensate for the missing factor \(\frac{1}{\sqrt{2 \pi}}\) in the pdf:

2
\(\pi\)
×
\(\sqrt{x}\)
÷

0.18449

---

I haven't checked but I assume that there are already programs for the HP-65 that deal with the normal distribution.

[Thomas Klemm]

From HP-65 STAT PAC 1, pp. 26-27:

Quote:NORMAL DISTRIBUTION

For a standard normal distribution

\(
\begin{align}
f(x) &= \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} \\
\\
Q(x) &= \frac{1}{\sqrt{2 \pi}} \int_{x}^{\infty}e^{-\frac{t^2}{2}}\mathrm{d}t \\
\end{align}
\)

For \(x \geqslant 0\), polynomial approximation is used to compute \(Q(x)\):

\(
\begin{align}
Q(x) =f(x) (b_1 t + b_2 t^2 + b_3 t^3 + b_4 t^4 + b_5 t^5) + \epsilon(x)
\end{align}
\)

where \(|\epsilon(x)| < 7.5 \times 10^{-8}\)

\(
\begin{align}
t &= \frac{1}{1 + rx} \\
r &= 0.2316419 \\
\end{align}
\)

\(
\begin{align}
b_1 &= .31938153 \\
b_2 &= -.356563782 \\
b_3 &= 1.781477937 \\
b_4 &= -1.821255978 \\
b_5 &= 1.330274429 \\
\end{align}
\)

Note: \(f(-x)=f(x), Q(-x) = 1 - Q(x)\)

Reference: Handbook of Mathematical Functions, Abramowitz and
Stegun, National Bureau of Standards, 1968

Examples:

\(
\begin{matrix}
1. & f(1.18) &= 0.20 \\
& Q(1.18) &= 0.12 \\
2. & f(2.28) &= 0.03 \\
& Q(2.28) &= 0.01 \\
\end{matrix}
\)

Calculations

Upper Limit

672
ENTER
.5
+

4000
ENTER
6
÷
-

LST x
5
×
6
÷
\(\sqrt{x}\)
÷

2.47487372-01

A
3.869098599-01

B
4.022655818-01

Lower Limit

662
ENTER
.5
-

4000
ENTER
6
÷
-

LST x
5
×
6
÷
\(\sqrt{x}\)
÷

-2.192031036-01

A
3.894719103-01

B
5.867540454-01

Difference

From this value we can subtract the result of the upper limit:

.4022655818
-

1.844884636-01


Comparison

We can compare this with the result from: Sum[PDF[BinomialDistribution[4000, 1/6], k],{k,662,672}]

0.1844441329465701303644228283076550814475904333448317000634794690...


Programs

These programs can be loaded directly into this HP-65 Microcode Emulator.
The first program has to be loaded and started once with A to load the constants.
The second program can then be used to calculate \(f(x)\) and \(Q(x)\).

STAT 1-10A 1

Code:

PROG
100
001: 23    : LBL
002: 11    : A
003: 83    : .
004: 02    : 2
005: 03    : 3
006: 01    : 1
007: 06    : 6
008: 04    : 4
009: 01    : 1
010: 09    : 9
011: 33 03 : STO 3
012: 01    : 1
013: 83    : .
014: 03    : 3
015: 03    : 3
016: 00    : 0
017: 02    : 2
018: 07    : 7
019: 04    : 4
020: 04    : 4
021: 02    : 2
022: 09    : 9
023: 33 04 : STO 4
024: 01    : 1
025: 83    : .
026: 08    : 8
027: 02    : 2
028: 01    : 1
029: 02    : 2
030: 05    : 5
031: 05    : 5
032: 09    : 9
033: 07    : 7
034: 08    : 8
035: 42    : CHS
036: 33 05 : STO 5
037: 01    : 1
038: 83    : .
039: 07    : 7
040: 08    : 8
041: 01    : 1
042: 04    : 4
043: 07    : 7
044: 07    : 7
045: 09    : 9
046: 03    : 3
047: 07    : 7
048: 33 06 : STO 6
049: 83    : .
050: 03    : 3
051: 05    : 5
052: 06    : 6
053: 05    : 5
054: 06    : 6
055: 03    : 3
056: 07    : 7
057: 08    : 8
058: 02    : 2
059: 42    : CHS
060: 33 07 : STO 7
061: 83    : .
062: 03    : 3
063: 01    : 1
064: 09    : 9
065: 03    : 3
066: 08    : 8
067: 01    : 1
068: 05    : 5
069: 03    : 3
070: 33 08 : STO 8
071: 24    : RTN
072: 35 01 : g NOP
073: 35 01 : g NOP
074: 35 01 : g NOP
075: 35 01 : g NOP
076: 35 01 : g NOP
077: 35 01 : g NOP
078: 35 01 : g NOP
079: 35 01 : g NOP
080: 35 01 : g NOP
081: 35 01 : g NOP
082: 35 01 : g NOP
083: 35 01 : g NOP
084: 35 01 : g NOP
085: 35 01 : g NOP
086: 35 01 : g NOP
087: 35 01 : g NOP
088: 35 01 : g NOP
089: 35 01 : g NOP
090: 35 01 : g NOP
091: 35 01 : g NOP
092: 35 01 : g NOP
093: 35 01 : g NOP
094: 35 01 : g NOP
095: 35 01 : g NOP
096: 35 01 : g NOP
097: 35 01 : g NOP
098: 35 01 : g NOP
099: 35 01 : g NOP
100: 35 01 : g NOP
CARD
6
Title: STAT 1-10A 1
A: 
B: 
C: 
D: 
E: 
HELP
1

END

STAT 1-10A 2

Code:

PROG
100
001: 23    : LBL
002: 11    : A
003: 33 01 : STO 1
004: 41    : ENTER
005: 71    : x
006: 02    : 2
007: 81    : /
008: 42    : CHS
009: 32    : f-1
010: 07    : LN
011: 35    : g
012: 02    : pi
013: 02    : 2
014: 71    : x
015: 31    : f
016: 09    : SQRT
017: 81    : /
018: 33 02 : STO 2
019: 24    : RTN
020: 23    : LBL
021: 12    : B
022: 34 01 : RCL 1
023: 00    : 0
024: 35 24 : g x>y
025: 22    : GTO
026: 01    : 1
027: 23    : LBL
028: 13    : C
029: 01    : 1
030: 34 01 : RCL 1
031: 34 03 : RCL 3
032: 71    : x
033: 61    : +
034: 35    : g
035: 04    : 1/x
036: 41    : ENTER
037: 41    : ENTER
038: 41    : ENTER
039: 34 04 : RCL 4
040: 71    : x
041: 34 05 : RCL 5
042: 61    : +
043: 71    : x
044: 34 06 : RCL 6
045: 61    : +
046: 71    : x
047: 34 07 : RCL 7
048: 61    : +
049: 71    : x
050: 34 08 : RCL 8
051: 61    : +
052: 71    : x
053: 34 02 : RCL 2
054: 71    : x
055: 24    : RTN
056: 23    : LBL
057: 01    : 1
058: 34 01 : RCL 1
059: 42    : CHS
060: 33 01 : STO 1
061: 13    : C
062: 01    : 1
063: 35 07 : g x<>y
064: 51    : -
065: 24    : RTN
066: 35 01 : g NOP
067: 35 01 : g NOP
068: 35 01 : g NOP
069: 35 01 : g NOP
070: 35 01 : g NOP
071: 35 01 : g NOP
072: 35 01 : g NOP
073: 35 01 : g NOP
074: 35 01 : g NOP
075: 35 01 : g NOP
076: 35 01 : g NOP
077: 35 01 : g NOP
078: 35 01 : g NOP
079: 35 01 : g NOP
080: 35 01 : g NOP
081: 35 01 : g NOP
082: 35 01 : g NOP
083: 35 01 : g NOP
084: 35 01 : g NOP
085: 35 01 : g NOP
086: 35 01 : g NOP
087: 35 01 : g NOP
088: 35 01 : g NOP
089: 35 01 : g NOP
090: 35 01 : g NOP
091: 35 01 : g NOP
092: 35 01 : g NOP
093: 35 01 : g NOP
094: 35 01 : g NOP
095: 35 01 : g NOP
096: 35 01 : g NOP
097: 35 01 : g NOP
098: 35 01 : g NOP
099: 35 01 : g NOP
100: 35 01 : g NOP
CARD
6
Title: STAT 1-10A 2
A: f(x)
B: Q(x)
C: 
D: 
E: 
HELP
1

END