Binomial Probability Distribution

[Dieter]

(05-26-2015 03:47 PM)Dave Britten Wrote:  First post updated with a new listing. Thanks for the algebra optimization tips.

Great. Here's my solution. It has 74 steps including its own COMB routine. If you want to use your own external function, simply replace steps 16...32 with your COMB command. The program then is down to 58 steps.

The code below implements another feature that I consider very useful especially for large n resp. wide ranges between k1 and k2.
Take for instance this example:

n=100, p=0,1, k=10...90
The result is 0,5487 and requires the summation of 80 single PDFs. The previous version of my program required 43,5 seconds for this.

However, there is not much sense in adding the PDFs for k=50, 60 or even 90. They are so small that they do not change the final result. At k=40 the PDF already has dropped to 2,47 E–15, so the addition loop can stop here, or even earlier.

This is checked by the program version below: if an additional PDF does not change the CDF sum, the program exits. In the example above the last added term is the PDF for k=36, which was 1,16 E–11 so that it did not change the final result of 0,5487. Which of course is returned much faster now: instead of 43,5 s only 18,5 s are required – less than half the time.

Finally, here is the code:

Code:

01 LBL"BINCDF"
02 x<y?
03 x<>y
04 INT
05 STO 04
06 RDN
07 INT
08 STO 03
09 RDN
10 STO 02
11 RDN
12 INT
13 STO 01
14 RCL 01
15 RCL 03
16 -
17 RCL 03
18 x<y?
19 x<>y     // leave min(n, n-k1) in Y
20 SIGN     // set x:=1
21 x>y?     // min(n, n-k1) = 0 ?
22 GTO 03   // then Comb=1, skip Comb loop
23 LBL 01
24 RCL Z    // calculate Comb(n, k1) iteratively
25 *
26 RCL Y
27 /
28 DSE Z
29 LBL 00   // dummy label = NOP
30 DSE Y
31 GTO 01
32 LBL 03
33 RCL 02
34 RCL 03
35 y^x
36 *
37 1
38 RCL 02
39 -
40 ST/ 02   // prestore p/(1-p)
41 RCL 01
42 RCL 03
43 -
44 STO 01
45 y^x
46 *
47 RCL 04
48 RCL 03
49 -
50 x=0?     // # PDF loops = 0 ?
51 GTO 04   // then quit
52 STO 04
53 x<>y
54 ENTER
55 LBL 02   // stack: PDF(k) CDF(k) z  t
56 RCL 01
57 RCL 02
58 *
59 *        // this way CDF(k) is pushed to T
60 SIGN
61 ST- 01
62 ST+ 03
63 x<> L
64 RCL 03
65 /        // stack: X=PDF(k), Y=CDF(k), Z=CDF(k-1), T=CDF(k-1)
66 +        // PDF in L, current and previous CDF in X and Y
67 x=y?     // result same as before?
68 GTO 04   // then quit
69 LastX    // else recall last PDF
70 DSE 04   // decrement counter
71 GTO 02   // and continue with next k
72 LBL 04
73 RDN
74 END

Now try your example:

Code:

120 [ENTER] 6 [1/x] 40 [ENTER] 120
XEQ"BINCDF"   =>   6,419629920 E–6

The result now is returned in 25 seconds. That's more than twice as fast as the standard approach since only PDFs between k=41 and k=59 had to be summed up. That's merely 19 loops instead of 80(!):

Code:

LastX  => 1,7254 E–16  (last summed = first negligible PDF(k))
RCL 03 => 59  (last k)
RCL 04 => 62  (number of saved loops + 1)

Dieter