3
$\begingroup$

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.

Improve this question
$\endgroup$
5
  • $\begingroup$ The closest built-in is CompoundPoissonDistribution[], but that assumes $N$ is Poisson-distributed, and not negative binomial-distributed as in your desired application. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Mar 7, 2018 at 7:39
  • $\begingroup$ @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. $\endgroup$
    – joejoejoejoe4
    Mar 7, 2018 at 7:49
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Mar 7, 2018 at 7:51
  • 2
    $\begingroup$ 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? $\endgroup$
    – wolfies
    Mar 8, 2018 at 16:46
  • $\begingroup$ Yes, in the case that $N$ is zero, then the whole sum should be zero as well. $\endgroup$
    – joejoejoejoe4
    Mar 9, 2018 at 18:36

2 Answers 2

Reset to default
1
$\begingroup$

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]

pmf of S pmf of W

Improve this answer
$\endgroup$
0
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.