Hi,
I am attempting to compute a binomial distribution based on formulation here.

Basically binomial distributions return the probability of an event occuring $k times. In $n attempts, where the probability of it occuring in a single attempt is $p.

I tried two methods that should give the identical result. One subroutine using Math::Pari via its 'binomial' function plus logarithmic function and the other one using brute combinatorial method. In the end it's the log-Math::Pari method which I intended to use, since it is able to handle large number.

But however the result given by log-Math::Pari function is different from the correct combinatorial Method. 
Sub Binom Comb        = 0.3125    #This is correct
Sub Binom Log         = 0.06868298952623189587  #wrong

[download]
What's wrong with my log-Math::Pari subroutine? It seems to me I have constructed log version mathematically in a right way. Or have I used the Math::Pari function wrongly?
#!/usr/bin/perl -w
use strict;
use Math::Pari qw(binomial);

my $n = 6;
my $k = 3; 
my $p = 0.5;

my $sub_binom_comb = binomial_comb($k,$n,$p);
my $sub_binom_log = binomial_log($k,$n,$p);

print "Sub Binom Comb       = $sub_binom_comb\n";
print "Sub Binom Log         = $sub_binom_log\n";


#----My Subroutines ------------
sub binomial_log{

    #Find binomial distributions using Log 
    #and Math::Pari

    my ($k, $n, $p) = @_;

    my $first = log(binomial($n,$k));  
    #The above binomial function from Math::Pari

    my $second = $k * log($p);
    my $third = ($n-$k) * log (1-$p);
    
    my $log_Prob = $first + $second + $third;
      
    my $Prob = 10 ** ($log_Prob);
    return $Prob;

}


sub binomial_comb {
    #With combinatorial method
    #using brute factorial

    my ($x, $n, $p) = @_;
    return unless $x >= 0 && $x == int $x && $n > 0 &&
    $n == int $n && $p > 0 && $p <1;
    return choose($n,$k) * ($p ** $x) * ((1-$p) ** ($n-$x));

}

sub choose {
    my ($n,$k) = @_;
    my ($result,$j) = (1,1);
    
    return 0 if $k>$n||$k<0;
    $k = ($n - $k) if ($n - $k) <$k;
    
    while ($j <= $k ) {
        $result *= $n--;
        $result /= $j++;
    }
    return $result;
}

[download]
Regards,
Edward

In reply to Getting Binomial Distribution under Math::Pari (log) and Combinatorial Method by monkfan

Title:
Use:  <p> text here (a paragraph) </p>
and:  <code> code here </code>
to format your post, it's "PerlMonks-approved HTML":


  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.