Convolve discrete random variables efficiently

Question

Assume $X_i$ are random variables, which are identically and independently distributed and obtain values in $\{1,..,n\}$. We know their distribution:

$P[X_i=k]=p_k$

The question is how to efficiently compute the probability that the sum of $m$ such random variables (the convolution) equals a given integer $r$, that is

$P[\sum_i^m X_i =r]$

Is there an efficient way for this computation?

An example: compute $P[\sum_i^m X_i =r]$ with $r=100$, $m=40$, $n=10$ and the distribution

$p_k = \left( \Pi_{j=1}^{k-1} \frac{e^{-j}}{1 + e^{-j}} \right) * \begin{cases} \frac{1}{1 + e^{-k}} &\text{if}\ 1 \leq k < n \\ 1 &\text{if}\ k=n \end{cases} $

I tried a "brute force" approach, which is not very efficient with large numbers. With the code below, the example takes about 15 seconds on my computer. The idea was simply to calculate all sets which $m$ elements sum up to $r$, and compute the probabilities that the single random variables obtain the values of the set.

(*the probability that a single random variable has value k (between 1 and n)*)
p[k_, n_] := 
  p[k, n] = 
   N@Product[E^(-j)/(1 + E^(-j)), {j, 1, k - 1}]*
    Piecewise[{{1/(1+ E^(-k)), k < n}, {1, k == n}}];

(*the probability that the m-fold convolution results in the value r*)
pconvolute[r_, m_, n_] :=
  Module[
   {prob, permuts, partitions},

   (*partitions of m integers summing up to r*)
   partitions = IntegerPartitions[r, {m}, Range[n]];

   (*number of permutations of a partition*)
   permuts[partition_] := 
    Factorial[Length[partition]]/
     Apply[Times, Factorial /@ Tally[partition][[All, 2]]];

   (*probability for all permutations of a partition*)
   prob[partition_] := 
    permuts[partition] Apply[Times, p[#, n] & /@ partition];

   (*sum up the probabilities*)
   Total[prob /@ partitions]
   ];

(*the example*)
AbsoluteTiming[pconvolute[100, 40, 10]]

Answer 1

You are given a discrete random variable with pmf $f(x)$, defined on domain of support $x = 1, \dots, n$. I will work here with $n = 5$ to keep the output neat, but will conduct timings with $n = 10$ (as per your example). For a given value of $n$, we can enter your pmf $f(x)$ in a List Form as:

f = With[{n = 5}, 
  Table[ Product[E^-j/(1 + E^-j), {j, 1, k-1}] If[k<n, 1/(1+E^-k), 1], 
{k, 1, n}] ] // Simplify

and domain of support:

domain[f] = {x, Range[5]}  && {Discrete};

Here is a plot of your pmf:

Check that your pmf is well-defined and sums to 1:

Plus @@ f // Simplify

1

All good.

Solution

The pgf (probability generating function) is $E[t^X]$, which can be evaluated in Mma as say:

 pgf = (t^Range[5]) . f 

Let $S = \sum_{i=1}^m X_i$ denote the sum of $m$ such independent random variables. Then, the pgf for the sum $S$ of $m$ such variables is:

$$\begin{align*}\displaystyle \text{pgfS} &= E[t^S] \\&= E[t^{X_1}] E[t^{X_2}] \dots E[t^{X_{m}}] \text{ (... by independence)} \\ &= \text{pgf }^m \end{align*}$$

So, for the sum of $m = 40$ variables, the pgf of the sum is:

  pgfS = pgf^40;

Finally, given pgfS, we can find $P(S = s)$, for any desired value $s$. For example, $P(S = 55)$ is:

Coefficient[pgfS, t, 55] // N // AbsoluteTiming

{0.002342, 0.0756939}

i.e. this takes 0.002342 seconds to evaluate.

Timing comparison for the OP's example

With $n = 10$, $m = 40$ and $P(S=100)$:

 AbsoluteTiming[pconvolute[100, 40, 10]]   (* OP code *)

{6.31383, 1.21213*10^-26}

Above approach:

Coefficient[pgfS, t, 100] // N // AbsoluteTiming

{1.27931, 1.21213*10^-26}

The two key advantages of the pgf approach are:

• Faster (about 5 times faster here)
• Exact: if you remove the //N, you will get back the exact symbolic solution for $P(S=100)$.

Answer 2

Here's a purely numerical approach, related to @wolfies more abstract method. I use the convolution theorem and the FFT:

convolveNtimes[list_, n_] := 
  Join[ConstantArray[0.0, n - 1], 
    Normalize[
      Re[InverseFourier[
        Fourier[Join[list, 
          ConstantArray[0.0, (Length[list] - 1)*(n - 1)]]]^n]], Total]]

(All the fungus around the FFT stuff is getting padding and normalization right, and killing the imaginary errors)

Now, calculate the distribution. I memoize it to speed subsequent calculations.

pdist[m_, n_] := pdist[m, n] = convolveNtimes[Table[p[k, n], {k, 1, n}], m]

The final function just selects the element from the distribution:

pffconvolve[r_, m_, n_] := pdist[m, n][[r]]

This is very fast:

AbsoluteTiming[pffconvolve[80, 40, 10]]
(* {0.00034, 1.04703*10^-11} *)

Compare to your original:

AbsoluteTiming[pconvolute[80, 40, 10]]
(* {0.533918, 1.04703*10^-11} *)

Memoization makes subsequent accesses to the same distribution even faster:

AbsoluteTiming[pffconvolve[81, 40, 10]]
(* {7.*10^-6, 2.5095*10^-12} *)

This method has the disadvantage of absolute errors of order 10^-17:

pffconvolve[100, 40, 10]
(* -2.60968*10^-17 *)

If you don't need accurate estimates of extremely small probabilities, I think this is a good, fast method.

Answer 3

The approach with the moment generating functions is definitely the way to go, but out of curiosity as a programming problem, here is an alternate. The IntegerPartition function and the Tally are both fast, a lot faster than I was expecting, so I was curious to see if the final operation on it could be sped up. First precompute a vector of the probabilities.

pkn = (n = 10; Table[p[k, n], {k, 1, n}])

Use this in the function.

pconvolute4[r_, m_, n_] := 
  Module[{prob, permuts partitions},

   (*partitions of m integers summing up to r*)
   partitions = Tally /@ IntegerPartitions[r, {m}, Range[n]];

   prob[partition_] := ((Multinomial @@ (partition[[All, 2]]))) * 
         Times @@ Map[pkn[[(# // First)]]^(# // Last) &, partition];

   (*sum up the probabilities*)
   Total[prob /@ partitions]
   ];

This routine

AbsoluteTiming[pconvolute4[100, 40, 10]]
(*    {2.58898, 1.21213*10^-26}    *)

Your original routine

(*    {5.578, 1.21213*10^-26}    *)

EDIT

As an OBTW on some lessons learned in fiddling with this.

• Prerunning the probability function and turning it into a vector, vice using memoization (which should have preset all of the function values in the first few calls), resulted in a significant overall saving in time (5-10%)
• Letting Mathematica calculate the permutations via the Multinomial function was 20% saving.
• Tried using the MultinomialDistribution and PDF functions to let Mathematica do all the work, but this was costly by an order of magnitude or more.
• As a preparation for the MultinomialDistributioncall, needed to run BinCounts over the partitions, which proved to be painfully slow. Proved better to home-grow one using Tally first, although ultimately the MultinomialDistribution and PDF call was too expensive in the first place.