Beefy Boxes and Bandwidth Generously Provided by pair Networks
good chemistry is complicated,
and a little bit messy -LW

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]
[Eily]: I still don't understand how the Turkish AA fit into the German+Czech joke though :P
[LanX]: new Firefox + cb sidebar do random auto expand on submit
[LanX]: probably need to start pm discussion
[LanX]: they have a constitutional referendum in turkey, kind of "do you want a dictator" and everybody opting no gets problems ...
[Corion]: LanX: Random Auto Expand?
[Corion]: LanX: Well, everybody opting "yes" will also get problems. The question is more like "Do you want problems now or problems later?"
[LanX]: and the AA had posters with a big say NO with a small "to alcohol"
[ambrus]: LanX: is it the kind of free and secret vote where there's only one box you can check?

How do I use this? | Other CB clients
Other Users?
Others imbibing at the Monastery: (9)
As of 2017-03-27 12:14 GMT
Find Nodes?
    Voting Booth?
    Should Pluto Get Its Planethood Back?

    Results (320 votes). Check out past polls.