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]

---

---

• Are you posting in the right place? Check out Where do I post X? to know for sure.
• Posts may use any of the Perl Monks Approved HTML tags. Currently these include the following:

  <code> <a> <b> <big> <blockquote> <br /> <dd> <dl> <dt> <em> <font> <h1> <h2> <h3> <h4> <h5> <h6> <hr /> <i> <li> <nbsp> <ol> <p> <small> <strike> <strong> <sub> <sup> <table> <td> <th> <tr> <tt> <u> <ul>

• Snippets of code should be wrapped in <code> tags not <pre> tags. In fact, <pre> tags should generally be avoided. If they must be used, extreme care should be taken to ensure that their contents do not have long lines (<70 chars), in order to prevent horizontal scrolling (and possible janitor intervention).
• Want more info? How to link or How to display code and escape characters are good places to start.