Perl Monk, Perl Meditation PerlMonks

### Re: Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method

by tall_man (Parson)
 on Mar 23, 2005 at 16:22 UTC ( #441814=note: print w/replies, xml ) Need Help??

You could get more accuracy from the Math::Pari version if you did the powers and multiplies on the Pari side and didn't use logs:
```sub binomial_pari {
my (\$k, \$n, \$p) = @_;
my \$first = binomial(\$n,\$k);
# gpui is imported from Math::Pari.  It is x ** y.
my \$second = gpui(\$p,\$k);
my \$third = gpui(1.0 - \$p, \$n - \$k);
my \$Prob = \$first * \$second * \$third;
return \$Prob;
}
Here are the different results I get:
```Sub Binom Comb       = 0.3125
Sub Binom Log        = 0.3125000000000001128
Sub Binom Pari       = 0.3125000000000000000

Replies are listed 'Best First'.
Re^2: Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method
by monkfan (Curate) on Mar 24, 2005 at 01:35 UTC
When speed matters (apart from ability to deal with *BIG* number), the performance are:
```                 Rate  binomial_log binomial_pari binomial_comb
binomial_log   6164/s            --          -61%          -87%
binomial_pari 15709/s          155%            --          -66%
binomial_comb 45985/s          646%          193%            --
Regards,
Edward
If you just want an approximation, then performance may be better still if you use the fact that from Stirling's formula the log of n! is approximately
```log(n**n * exp(-n) * sqrt(2*PI*n))
= n*log(n) - n + (log(2 + log(PI) + log(n))/2
Plug that into the fact that n choose m is n!/(m!*(n-m)!) and you can come up with a good approximation that uses Perl's native arithmetic.

This approximation may well turn out to be the fastest approach for large numbers.

Update: From Mathworld I just found out that a significantly better approximation of n! is sqrt(2*n*PI + 1/3)*n**n/exp(n). Algebraically this is less convenient, but the improved accuracy may matter for whatever you're trying to do.

Create A New User
Node Status?
node history
Node Type: note [id://441814]
help
Chatterbox?
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others examining the Monastery: (6)
As of 2018-04-19 16:14 GMT
Sections?
Information?
Find Nodes?
Leftovers?
Voting Booth?
My travels bear the most uncanny semblance to ...

Results (74 votes). Check out past polls.

Notices?