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
[download]