For instance, it's easy to see that if the numbers are evenly distributed between 0 and 1, for N random numbers N / 2 would be have the highest frequency, while 0 and N would have the lowest.

Sorry for derailing this thread if you've got good reason to need the training set - just a thought.

Hi Joost, The training set isn't a problem to compute, and necessary to know what is the probability of obtaining any given sum. The problem is more on how to scale up the recursive convolution. FFRANK

0 + 0 = 0
0 + 1 = 1
0 + 2 = 2
1 + 0 = 1
1 + 1 = 2
1 + 2 = 3
2 + 0 = 2
2 + 1 = 3
2 + 2 = 4
[download]

I guess I am wondering if the sum to be gotten was 2, for example, then 3 cases would qualify with your algorithm.

 
0 + 2 = 2
1 + 1 = 2
2 + 0 = 2
[download]

The numbers you are generating now are also known as the cross product of n sets. Perhaps you just want to test all the possible combinations instead?

0 + 1 = 1
0 + 2 = 2
1 + 2 = 3
[download]

As far as your code scaling up, it depends on how large the number of trials. If you use the cross product, with 10,000 numbers there will be 10,000 * 10,000 trials. On the other hand, if you tested the combinations there would be 10,000 * 9,999 / 2 trials - half as many.

Chris

Hi Cristoforo, you getting the point right. If, for example, only three "object sums" are possible (0, 1 and 2).

$nElements = 1:
Sum = 0 -> 0
Sum = 1 -> 1
Sum = 2 -> 2

$nElements = 2:
Sum = 0 -> 0+0
Sum = 1 -> 0+1 1+0
Sum = 2 -> 0+2 1+1 2+0
Sum = 3 -> 1+2 2+1
Sum = 4 -> 2+2

$nElements = 3:
Sum = 0 -> 0+0+0
Sum = 1 -> 0+0+1 0+1+0 1+0+0
Sum = 2 -> 0+0+2 0+1+1 0+2+0 1+0+1 1+1+0 2+0+0
Sum = 3 -> 0+1+2 0+2+1 1+0+2 1+1+1 1+2+0 2+0+1 2+1+0
Sum = 4 -> 0+2+2 1+1+2 1+2+1 2+0+2 2+1+1 2+2+0
Sum = 5 -> 1+2+2 2+1+2 2+2+1
Sum = 6 -> 2+2+2

And so on...
[download]

Of course, there is some redundancy in there, for ex. Element = 3, sum = 1 could be treated as (3*((p0)**2)*p1). But, how to generate this automatically and fast ?

Thanks, FFRANK

It's not in Perl, but some of the best sources for practical algorithms are the Numerical Recipes books, which have online editions. You might also look at Mastering Algorithms With Perl

Hi David, I believe the multinomial formula does not apply here, as we deal with the probability of a sum given n trials (over the combinations of outcomes that can generate the sum over n trial).

Given: number of trials = t, sum of random variables = Z, the discrete convolution is:

Pr[Yt = Z] = ∑ (e = 0 .. Z) Pr[Yt-1 = Z-e] * Pr[e]

FFRANK

#!/usr/bin/perl
use strict;
use warnings;
use Set::CrossProduct;
use List::Util 'sum';

my $targetSum = 4;
my @freq = qw (0.1 0.2 0.7);
my $nElements = 2;
my @combo = ([0..$#freq]) x $nElements;

my $cross = Set::CrossProduct->new( \@combo );

my @ways;
while( my $array_ref = $cross->get ) {
    push @ways, $array_ref if sum(@$array_ref) == $targetSum;
}

#Now calculate P-value
my $pVal;
for my $aref (@ways) {
    my $p = 1;
    for my $operand (@$aref) {
        $p *= $freq[$operand];
        }
    $pVal += $p;
    }
print "P-value(sum of $nElements = $targetSum) = $pVal\n";
[download]

use List::Util qw/sum/;
my @a; # holds the known distribution
sub conv {
 my ($n) = @_;
 return sum map {
  ($_ < @a) ? ((($n - $_) < @a) ? $a[$_] * $a[$n - $_]
                                : 0)
            : 0
 } 0 .. $n;
}
[download]

BTW P-value has a different meaning. You really just want to call it probability.

Enjoy!

#!/usr/bin/perl -w
use strict;

#Set no. of elements, frequencies and the sum value
my @freq = qw (0.1 0.2 0.7);

my $nElements = 2;
my $targetSum = 4;

#
# If you are only interested in the first few $targetSums
# and @freq is a large array, you can speed up the computation
# and reduce space requirements by setting $maxSum to the
# largest conceivable $targetSum
#
my $maxSum=400;


sub conv
{   my ($p, $q) = @_;
    my ($np, $nq) = (scalar @$p, scalar @$q);
    my @ans = ();
    push @ans, 0 for (1..$np+$nq);
    my $ip = 0;
    for my $vp (@$p)
    {    my $iq = 0;
        for my $vq (@$q)
        {   $ans[$ip+$iq] += $vp*$vq;
            $iq++;
            last if $ip+$iq > $maxSum;
        }
        $ip++;
    }
    return \@ans;
}

# see http://en.wikipedia.org/wiki/Exponentiation_by_squaring
sub convn
{   my ($n, $p) = @_;
    my $ans = [1];
    while ($n)
    {    $ans = conv($ans, $p) if $n & 1;    # if n is odd
        $p = conv($p, $p);
        use integer;
        $n = $n >> 1;
    }
    return $ans;
}

my $ans = convn($nElements, \@freq);

print "Prob(sum of $nElements = $targetSum) = $ans->[$targetSum]\n";

#Prints: Prob(sum of 2 = 4) = 0.49

print convn(200, [0.1, 0.2, 0.3, 0.2, 0.1, 0.1])->[400];
#Prints: 0.000210606620621573
[download]

Hi b4,

Your solution is great and fast, Thanks very much !

FFRANK