I believe this to be a correct implementation.
I've tried to retain the auditability of herveus' implementation by aliasing perl's working variables to names that marry with those used in the formula, whilst avoiding copying lots of arrays.
For the small numbers in the samples it makes little difference, but by the time you get to a dozen or more, the double duplication and the spliceing in a recursive subroutine become a severe resource drain. If the numbers involved are anything like those typical for biogenetic work, they would become untenable quite quickly.
Update: Use the non-aliasing version.
The aliasing version consumes prodigous amount of memory,
This version happily processes
#! perl -slw
use strict;
use List::Util qw[ reduce ];
sub P;
sub P{ ## Non-aliasing
warn "@_\n";
return $_[ 0 ] if @_ == 1;
return reduce {
$a += ( $_[$#_ - $b + 1] - $_[$#_ - $b] )
* P2 @_[ 0 .. $b-1, $b+1 .. $#_ ];
} 0, 1 .. $#_;
}
=do not use this version
This version looks pretty and should be logically identical to the abo
+ve and it appears to work okay
for short lists.
But through what I think is a bug in multiplicity Perl
the use of c<local> cause huge memory consumption:
> 750 MB for an input list of 10 items!?
sub P { ## Pretty, aliasing, voracious memory consumer!
return $_[ 0 ] if @_ == 1;
our( @r, $i, $a, $b, $sigma );
local *r = *_;
local *i = *b;
local *sigma = *a;
my $n = $#r;
return reduce{
$sigma += ($r[$n-$i+1] - $r[$n-$i]) * P @r[0 .. $i-1, $i+1 ..
+$n];
} 0, 1 .. $n;
}
=cut
my @samples = (
## r1 r2 r3
[ qw[ 0.11 0.07 0.19 ] ],
[ qw[ 0.43 0.31 0.37 ] ],
[ qw[ 0.93 0.78 0.82 ] ],
[ qw[ 0.91 0.12 0.15 ] ],
[ qw[ 0.52 0.32 0.18 ] ],
[ qw[ 1.0 1.0 1.0 ] ],
[ qw[ 0.5 0.5 0.5 ] ],
[ qw[ 0.0 0.0 0.0 ] ],
[ qw[ 0.19 0.11 0.07 ] ],
[ qw[ 0.43 0.37 0.31 ] ],
[ qw[ 0.93 0.82 0.78 ] ],
[ qw[ 0.91 0.15 0.12 ] ],
);
print "P( @$_ ) = ", P( @$_ ) for @samples;
__END__
P( 0.11 0.07 0.19 ) = 0.001232
P( 0.43 0.31 0.37 ) = 0.004644
P( 0.93 0.78 0.82 ) = 0.016833
P( 0.91 0.12 0.15 ) = 0.547183
P( 0.52 0.32 0.18 ) = 0.045552
P( 1.0 1.0 1.0 ) = 0
P( 0.5 0.5 0.5 ) = 0
P( 0.0 0.0 0.0 ) = 0
P( 0.19 0.11 0.07 ) = 0.002128
P( 0.43 0.37 0.31 ) = 0.004644
P( 0.93 0.82 0.78 ) = 0.016833
P( 0.91 0.15 0.12 ) = 0.547183
Examine what is said, not who speaks.
"Efficiency is intelligent laziness." -David Dunham
"Think for yourself!" - Abigail
"Memory, processor, disk in that order on the hardware side. Algorithm, algorithm, algorithm on the code side." - tachyon
-
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.