Beefy Boxes and Bandwidth Generously Provided by pair Networks
good chemistry is complicated,
and a little bit messy -LW
 
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 taking refuge in the Monastery: (14)
As of 2014-07-10 22:33 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    When choosing user names for websites, I prefer to use:








    Results (217 votes), past polls