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

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;
}
[download]

Here are the different results I get:

Sub Binom Comb       = 0.3125
Sub Binom Log        = 0.3125000000000001128
Sub Binom Pari       = 0.3125000000000000000
[download]

Comment on Re: Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method Select or Download Code

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% --` [download] Regards, Edward	[reply] [d/l]
Re^3: Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method by tilly (Archbishop) on Mar 24, 2005 at 03:08 UTC
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(2PIn)) = nlog(n) - n + (log(2 + log(PI) + log(n))/2` [download] 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(2nPI + 1/3)n*n/exp(n)`. Algebraically this is less convenient, but the improved accuracy may matter for whatever you're trying to do.	[reply] [d/l] [select]


P is for Practical
	PerlMonks