# Distribution of Function of Random Sum of Random Variables

Let $S = \sum_{i=1}^n X_i$ where:

• Each $X_i$ is independently 3 or 9 (with equal probability), and
• The sample size $n$ is itself an independent random variable where $N \sim \text{NegativeBinomial}(r,p)$ e.g. $r = 5$ and $p = \frac34$

Let $W = \begin{cases}S-10 & S > 10 \\ 0 & S \leq 10 \end{cases} . \quad$

Find: (i) the pmf of $S\quad$ (ii) the pmf of $W \quad$ (iii) If not the former, can one at least find $\mathbb{E}[W]$ ?

How can I obtain the distribution of $W$ using mathematica? Although it may not be able to give me a reasonable PDF, can I at least find the expected value of $W$ with mathematica? Also, can I also draw random samples from the distribution of $W$? If so, how would I do this? Thank you.

• The closest built-in is CompoundPoissonDistribution[], but that assumes $N$ is Poisson-distributed, and not negative binomial-distributed as in your desired application. – J. M.'s ennui Mar 7 '18 at 7:39
• @J.M., the point is that I can define a function of several random variables and then draw empirical samples from them. Can I define my own random variables, define some function of them, and then treat that as a new random variable in Mathematica? That's the goal. – jippyjoe4 Mar 7 '18 at 7:49
• Yes, I understood what you wanted; I was only saying that my perusal of the docs did not show me a more general version of CompoundPoissonDistribution[] which would do what you want. – J. M.'s ennui Mar 7 '18 at 7:51
• To the OP: Before posting a solution, there is an aspect to your question that causes some 'complexity' which you may not be intending (or realising), and which I seek clarification on. You define $S$ as the sum from 1 to $N$ random variables, where $N$ is a NegBinomial random variable. A NegBinomial random variable includes 0 in its domain of support, so some of your runs will have 0 values in them. Such runs will return neither a 3 nor a 9 ... they will have length of zero, and the sum of the run will be 0. Is that actually what you are intending? – wolfies Mar 8 '18 at 16:46
• Yes, in the case that $N$ is zero, then the whole sum should be zero as well. – jippyjoe4 Mar 9 '18 at 18:36

Here is a basic principles approach to obtaining the probability mass function for both $S$ and $W$ (which uses the clarification asked for by @wolfies):

(* Set parameters *)
r = 5;
p = 3/4;
pBin = 1/2;

(* Function to generate the combinations of n (sample size) and
k (number of successes) that result in S: {n,k} as S = 9n - 6k *)
f[s_] := Table[{(s - 2 i)/3, s/3 - i}, {i, 0, s/3, 3}]

(* Individual probability for a combination of n and k *)
pr[n_, k_] := PDF[BinomialDistribution[n, pBin], k]*
PDF[NegativeBinomialDistribution[r, p], n]

(* Probability mass function for S *)
(* Note that only non-negative multiples of 3 have a positive probability *)
prS[s_] := If[IntegerQ[s/3], Total[pr[#[[1]], #[[2]]] & /@ f[s]], 0]

(* Probability mass function for W *)
prW[w_] := If[w == 0, prS[0] + prS[3] + prS[6] + prS[9],
If[IntegerQ[(w + 1)/3], prS[w + 10], 0]]

(* Approximate expectation of W *)
N[Sum[w prW[w], {w, 0, 500}]]
(* 3.79827 *)


Plots of the probability mass functions are as follows:

ListPlot[Table[prS[s], {s, 0, 50}], PlotRange -> All,
Frame -> True, FrameLabel -> {"S", "Probability"}, Filling -> 0]
ListPlot[Table[prW[w], {w, 0, 50}], PlotRange -> All,
Frame -> True, FrameLabel -> {"W", "Probability"}, Filling -> 0]


Simulating random samples is straight-forward:

r = 5;
p = 0.75;
m = 10000000;
Nlist = RandomVariate[NegativeBinomialDistribution[r, p], m];
W = Ramp[Subtract[Total[RandomChoice[{3, 9}, #] & /@ Nlist, {2}], 10]];


You could obtain the empiric probability density function by

distro = EmpiricalDistribution[W];
ρ = PDF[distro, #] &;


However, using BinCounts is much faster:

density = N[BinCounts[W, {0, Max[W], 1}]/m];
ListPlot[density, PlotRange -> All, Filling -> 0]