Binomial Probability Distribution

05242015, 05:58 PM
(This post was last modified: 05262015 03:46 PM by Dave Britten.)
Post: #1




Binomial Probability Distribution
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 6sided 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.419629920E6 in just over a minute. Compare to the version on my 48 using a direct summation of terms and repeated COMB calls returning 6.4196298769E6. That's an error of 6.71E7%, 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 

05252015, 09:51 PM
Post: #2




RE: Binomial Probability Distribution
(05242015 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 

05252015, 10:11 PM
Post: #3




RE: Binomial Probability Distribution
(05252015 09:51 PM)Dieter Wrote:(05242015 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. 

05252015, 10:38 PM
(This post was last modified: 05252015 10:40 PM by Dieter.)
Post: #4




RE: Binomial Probability Distribution
(05252015 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:
Dieter 

05252015, 10:51 PM
Post: #5




RE: Binomial Probability Distribution
(05252015 10:38 PM)Dieter Wrote:(05252015 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)*(nr). 

05262015, 03:47 PM
Post: #6




RE: Binomial Probability Distribution
First post updated with a new listing. Thanks for the algebra optimization tips.


05262015, 10:38 PM
(This post was last modified: 05262015 10:58 PM by Dieter.)
Post: #7




RE: Binomial Probability Distribution
(05242015 05:58 PM)Dave Britten Wrote: Running it with inputs 120, 1/6, 40, 120 yields 6.419629920E6 in just over a minute. Compare to the version on my 48 using a direct summation of terms and repeated COMB calls returning 6.4196298769E6. That's an error of 6.71E7%, 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:
Dieter 

05272015, 12:51 PM
(This post was last modified: 05272015 01:13 PM by Dieter.)
Post: #8




RE: Binomial Probability Distribution
(05262015 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 

05282015, 01:44 AM
Post: #9




RE: Binomial Probability Distribution
Good idea with the early exit. I'll have to play around with that a bit.


05282015, 11:11 AM
Post: #10




RE: Binomial Probability Distribution
(05282015 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 tendigit value of 1/6. :) Dieter 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)