Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked
 
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??


in reply to 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; }
Here are the different results I get:
Sub Binom Comb = 0.3125 Sub Binom Log = 0.3125000000000001128 Sub Binom Pari = 0.3125000000000000000


Comment on Re: Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method
Select or Download Code
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.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://441814]
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others romping around the Monastery: (15)
As of 2015-07-06 13:33 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (74 votes), past polls