Binomial Probability Distribution - Printable Version +- HP Forums ( https://www.hpmuseum.org/forum)+-- Forum: HP Software Libraries ( /forum-10.html)+--- Forum: HP-41C Software Library ( /forum-11.html)+--- Thread: Binomial Probability Distribution ( /thread-3963.html) |

Binomial Probability Distribution - Dave Britten - 05-24-2015 05:58 PM
This program calculates the odds of a specified number of successes from independent trials with a known probability of success. This requires a combinations function named COMB that calculates nCr. If you need one, see here for a few options. Inputs: t: Number of trials z: Individual probability of success y: Minimum number of successes x: Maximum number of successes Probability will be returned to x. Example: What is the probability of obtaining at least two even numbers when rolling five standard 6-sided dice? 5 ENTER .5 ENTER 2 ENTER 5 XEQ BINPR Result: 0.8125 (81.25%) Edit: Rewritten with Dieter's suggested optimizations. It's a bit longer, but faster. I'm sure there's still room for more elegant stack/register use here. I'll wait until those shake out before I pencil the new listing into my Moleskine. Running it with inputs 120, 1/6, 40, 120 yields 6.419629920E-6 in just over a minute. Compare to the version on my 48 using a direct summation of terms and repeated COMB calls returning 6.4196298769E-6. That's an error of 6.71E-7%, which isn't a terrible sacrifice for having the program finish inside of a lunar cycle (or possibly solar orbit) now. Code: `t: Number of trials` RE: Binomial Probability Distribution - Dieter - 05-25-2015 09:51 PM
(05-24-2015 05:58 PM)Dave Britten Wrote: This program calculates the odds of a specified number of successes from independent trials with a known probability of success. Dave, it looks like your program evaluates COMB in every single loop so that it should not run particulary fast. But there is a simple fix: once nCr is known, the following nCr+1 can be evaluated very easily from this. Even the whole Binomial PDF for x+1 can be calculated from the PDF of the previous x, so that the slow power function in the loop can be avoided as well. Dieter RE: Binomial Probability Distribution - Dave Britten - 05-25-2015 10:11 PM
(05-25-2015 09:51 PM)Dieter Wrote:(05-24-2015 05:58 PM)Dave Britten Wrote: This program calculates the odds of a specified number of successes from independent trials with a known probability of success. Thanks for the tip. I'll see if I can incorporate that. I've just always been doing the looped nCr approach. RE: Binomial Probability Distribution - Dieter - 05-25-2015 10:38 PM
(05-25-2015 10:11 PM)Dave Britten Wrote: Thanks for the tip. I'll see if I can incorporate that. I've just always been doing the looped nCr approach. I just tried it and it works well. Let n and p be the distribution parameters, and k1 resp. k2 the lower/upper bound of the random variable. Then this is the basic idea: - Calculate Comb(n, k1).
- Use this to calculate the probability B(k1,n,p).
This is the initial value of the sum you want to get. - Repeat the following loop k2–k1 times:
- Calculate B(k+1) = B(k) * (n–k)/(k+1) * p/(1–p)
- Add this to the sum, increase k by 1.
Dieter RE: Binomial Probability Distribution - Dave Britten - 05-25-2015 10:51 PM
(05-25-2015 10:38 PM)Dieter Wrote:(05-25-2015 10:11 PM)Dave Britten Wrote: Thanks for the tip. I'll see if I can incorporate that. I've just always been doing the looped nCr approach. That appears to mesh with what I came up with: COMB(n, r+1)=COMB(n, r)/(r+1)*(n-r). RE: Binomial Probability Distribution - Dave Britten - 05-26-2015 03:47 PM
First post updated with a new listing. Thanks for the algebra optimization tips. RE: Binomial Probability Distribution - Dieter - 05-26-2015 10:38 PM
(05-24-2015 05:58 PM)Dave Britten Wrote: Running it with inputs 120, 1/6, 40, 120 yields 6.419629920E-6 in just over a minute. Compare to the version on my 48 using a direct summation of terms and repeated COMB calls returning 6.4196298769E-6. That's an error of 6.71E-7%, which isn't a terrible sacrifice for having the program finish inside of a lunar cycle (or possibly solar orbit) now. The HP41 works with 10 digits, the HP48 with 12. For the former 1/6 equals 0,166 666 666 7, for the latter it's 0,166 666 666 667. This means: - The correct result for p=1/6 is 6,4196298765918... E–6
- The correct 10-digit result for p=0,1666666667 is 6,419629908 E–6.
So the HP41 is pretty close. - The correct 12-digit result for p=0,166666666667 is 6,41962987691 E–6.
So the HP48 is very close. Since COMB is an internal function and thus most probably exact in all 12 digits, this could be expected. - Compared to the true result with p exactly 1/6, both calculators offer similar accuracy: 8 out of 10 digits resp. 10 out of 12. ;-)
Dieter RE: Binomial Probability Distribution - Dieter - 05-27-2015 12:51 PM
(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"` Now try your example: Code: `120 [ENTER] 6 [1/x] 40 [ENTER] 120` 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))` Dieter RE: Binomial Probability Distribution - Dave Britten - 05-28-2015 01:44 AM
Good idea with the early exit. I'll have to play around with that a bit. RE: Binomial Probability Distribution - Dieter - 05-28-2015 11:11 AM
(05-28-2015 01:44 AM)Dave Britten Wrote: Good idea with the early exit. I'll have to play around with that a bit. Here is an update with the Comb loop counting up the denominator: Code: `.. ...` That's three more steps, on the other hand the RCL 01 command in line 14 now can be omitted. So all in all it's two steps more. But accuracy improves – try your example once again: Code: `120 [ENTER] 6 [1/x] 40 [ENTER] 120` That's the exact result for a ten-digit value of 1/6. :-) Dieter |