Beefy Boxes and Bandwidth Generously Provided by pair Networks
"be consistent"

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

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% --
      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?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://441814]
and the radiator hisses contentedly...

How do I use this? | Other CB clients
Other Users?
Others browsing the Monastery: (2)
As of 2017-04-28 09:33 GMT
Find Nodes?
    Voting Booth?
    I'm a fool:

    Results (519 votes). Check out past polls.