Characteristic functions can make computations involving the sums and differences of random variables really easy. Mathematica has lots of functions to work with statistical distributions, including a builtin to transform a distribution into its characteristic function.
I'd like to illustrate this with two concrete examples: (1) Suppose you wanted to determine the results of rolling a collection of dice with differing numbers of sides, e.g., roll two six-sided dice plus one eight-sided die (i.e., 2d6+d8)? Or (2) suppose you wanted to find the difference of two dice rolls (e.g., d6-d6)?
An easy way to do this would be to use the characteristic functions of the underlying discrete uniform distributions. If a random variable $X$ has a probability mass function $f$, then its characteristic function $\varphi_X(t)$ is just the discrete Fourier Transform of $f$, i.e., $\varphi_X(t) = \mathcal{F}\{f\}(t) = E[e^{i t X}]$. A theorem tells us:
If the independent random variables $X$ and $Y$ have corresponding probability mass functions $f$ and $g$, then the pmf $h$ of the sum $X + Y$ of these RVs is the convolution of their pmfs $h(n) = (f \ast g)(n) = \sum_{m=-\infty}^\infty f(m) g(n-m)$.
We can use the convolution property of Fourier Transforms to restate this more simply in terms of characteristic functions:
The characteristic function $\varphi_{X+Y}(t)$ of the sum of independent random variables $X$ and $Y$ equals the product of their characteristic functions $\varphi_{X}(t) \varphi_{Y}(t)$.
This Mathematica function will make the characteristic function for an s-sided die:
MakeCf[s_] :=
Module[{Cf},
Cf := CharacteristicFunction[DiscreteUniformDistribution[{1, s}],
t];
Cf]
The pmf of a distribution can be recovered from its characteristic function, because Fourier Transforms are invertible. Here is the Mathematica code to do it:
RecoverPmf[Cf_] :=
Module[{F},
F[y_] := SeriesCoefficient[Cf /. t -> -I*Log[x], {x, 0, y}];
F]
Continuing our example, let F be the pmf that results from 2d6+d8.
F := RecoverPmf[MakeCf[6]^2 MakeCf[8]]
There are $6^2 \cdot 8 = 288$ outcomes. The domain of support of F is $S=\{3,\ldots,20\}$. Three is the min because you're rolling three dice. And twenty is the max because $20 = 2 \cdot 6 + 8$. If you want to see the image of F, compute
In:= F /@ Range[3, 20]
Out= {1/288, 1/96, 1/48, 5/144, 5/96, 7/96, 13/144, 5/48, 1/9, 1/9, \
5/48, 13/144, 7/96, 5/96, 5/144, 1/48, 1/96, 1/288}
If you want to know the number of outcomes that sum to 10, compute
In:= 6^2 8 F[10]
Out= 30
If the independent random variables $X$ and $Y$ have corresponding probability mass functions $f$ and $g$, then the pmf $h$ of the difference $X - Y$ of these RVs is the cross-correlation of their pmfs $h(n) = (f \star g)(n) = \sum_{m=-\infty}^\infty f(m) g(n+m)$.
We can use the cross-correlation property of Fourier Transforms to restate this more simply in terms of characteristic functions:
The characteristic function $\varphi_{X-Y}(t)$ of the difference of two independent random variables ${X,Y}$ equals the product of the characteristic function $\varphi_{X}(t)$ and $\varphi_{Y}(-t)$ (N.B. the negative sign in front of the variable t in the second characteristic function).
So, using Mathematica to find the pmf G of d6-d6:
G := RecoverPmf[MakeCf[6] (MakeCf[6] /. t -> -t)]
There are $6^2 = 36$ outcomes. The domain of support of G is $S=\{-5,\ldots,5\}$. -5 is the min because $-5=1-6$. And 5 is the max because $6-1=5$. If you want to see the image of G, compute
In:= G /@ Range[-5, 5]
Out= {1/36, 1/18, 1/12, 1/9, 5/36, 1/6, 5/36, 1/9, 1/12, 1/18, 1/36}