Calculate propabilities without recursion? Coin toss.

[Microcebus]

Dear Monks, I want to calculate the propability to get side 1 at least k times, tossing a coin n times.

The following code works fine with numbers up to about 150 but when I use numbers above, I get the result -1.#IND

Is there a better way to solve this problem?

#!/usr/bin/perl -w

$coin_toss=1000; # =n
$toss_side1=60; # =k

$total_p=0;
foreach($toss_side1..$coin_toss)
    {
    # calculate factorials
    $n_fac=fac($coin_toss);
    $n_minus_k_fac=fac($coin_toss-$_);
    $k_fac=fac($_);
    
    # calculate propability for getting side 1 exactly $_ times
    $partial_p=($n_fac/($n_minus_k_fac*$k_fac))*(0.5**$_)*(0.5**($coin
+_toss-$_));
    
    # add propability to total propapility
    $total_p=$total_p+$partial_p;
    }
print"Propability to get side 1 at least $toss_side1 times tossing a c
+oin $coin_toss times: $total_p\n";
system("pause");
exit;

# subroutine for calculating factorials
sub fac
    {
    $_[0]>1?$_[0]*fac($_[0]-1):1;
    }

[roboticus]

Microcebus:

I'd suggest computing the binomial expansion directly without factorials, as the factorials can get rather slow. I recently did so (with repellent's help in Re^5: Monte Carlo - Coin Toss in James_H's thread Monte Carlo - Coin Toss. Let me know if you have any trouble with it.

...roboticus

When your only tool is a hammer, all problems look like your thumb.

[Microcebus]

Could you show me how the script would look like? Unfortunately I'm just a biologist and a very poor programmer. Many Thanks...

[roboticus]

Microcebus:

As others have said, if you're a scientist, you should definitely find a programming language to become proficient at!

Anyway, here's a start. I don't know exactly what you're trying to do, so I just hacked together a couple things from the posts I cited in my previous response. Here ya go:

use strict;
use warnings;
use List::Util qw(reduce);
use bignum;

# Number of coins to toss
my $coin_toss=1000;

# We cumulative probability of side1 appearing anywhere between 1 and 
+60 times
my $side1=60;

# Compute the binomial coefficients for row 1000 of pascal's triangle
my @triangle = robo_2(1000);

print <<EOHDR;
Tails      Count     %
-----   ----------  ------
EOHDR
my $total_prob=0;
my $runs = reduce { $a + $b } @triangle;
for (my $i = 0; $i < @triangle; $i++) {
        my $prob = $triangle[$i]/$runs;
        $total_prob += $prob unless $i > $side1;

        print "$i\t$triangle[$i]\t",100*$prob, "\n";
}

print "\n\nProbability of 0..60 tails is: ", $total_prob*100;

sub robo_2 {
        my $row = shift;
        my @cols = (1);
        ++$row;
        $cols[$_] = $cols[$_-1] * ($row-$_)/$_ for 1 .. $row-1;
        return @cols;
}

Note: The output gets pretty darned wide with 1000 tosses and bignum in the picture...

...roboticus

When your only tool is a hammer, all problems look like your thumb.

[broomduster]

see also Stirling's Approximation

[educated_foo]

Ratazong answered your recursion question above. You're probably getting -1 because numbers in intermediate results are too big. If you're trying to calculate something like

  N! / (K! * (N-K)!)
= (N * (N-1) * ... * (K+1)) / ((N-K) * (N-K-1) * ... * 1)

you need to either mix up the multiplications and divisions so the intermediate results stay near the magnitude of the final result, or operate in log-space:

= exp((log N + log (N-1) + ... + log (K+1)) - (log (N-K) + log (N-K-1)
+ + ... + 0))

[Ratazong]

Your function fac results in n * (n-1) * (n-2) ... * 1. This can be easily calculated using a loop:

my $res = 1;
for(my $i = 2; $i <= $_[0]; $i++) { $res *= $i; }
return $res;

update: However that is probably not solving your issue with -1.#IND, which indicates that your numbers are getting too big; others propose to add

use bignum;

to the beginning of your script...

[Microcebus]

Thanks allot, using bignum works. But makes the whole very very slowly. I want to run this script up to 10000 times with n=up to 5000. This will take weeks... Any suggestions?

[happy.barney]

Memoize (Attribute::Memoize) fac.
Btw: n!/k! = List::Util::reduce { $a * $b } $k+1 .. $n

[ikegami]

It wouldn't surprise me if 5000! is larger than the number of atoms in the universe. Would you care to reconsider that?

[Anonymous Monk]

This is the factorial function, which is a special case of the gamma function. If you are going to use this number in ratios, especially in combination with other gamma functions, what you might be better off calculating is the log-gamma function.

[ikegami]

That's a very naïve implementation. You calculate n! a total of (n-k+1) times. You calculate 0.5^x * 0.5^n-x just as many times even though it's always equal to 0.5ⁿ. At the very least, doing some simple math legwork first would go a long way.

[TomDLux]

A scientist who doesn't use and create software is a rare breed and likely headed for extinction.

I think you need to consider learning essential concepts in software design, software engineering and related topics. One site with articles specifically designed for scientists is http://software-carpentry.org/ .. Unfortunately, they deal with Python .. not that Python is bad for you, but it doesn't advance Perl.

[dwm042]

A scientist who doesn't use and create software is a rare breed and likely headed for extinction.

In fact, it's entirely the opposite in the experimental sciences. No one cares how well you handle data. Good scientists can hire programmers or postdocs to do that. But if you're magic with the test tubes, you have a job for life.