(HP-65) Binomial distribution with large input support - Dave Britten - 01-05-2020 06:37 PM
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
RE: (HP-65) Binomial distribution with large input support - Thomas Klemm - 06-25-2022 08:42 PM
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.
RE: (HP-65) Binomial distribution with large input support - Thomas Klemm - 06-26-2022 08:30 AM
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
|