Quantile for CompoundPoissonDistribution

Question

I need Quantile for CompoundPoissonDistribution, for example

Quantile[CompoundPoissonDistribution[1, GammaDistribution[100, 200]], 0.95]

I have Mathematica 9. Any idea how to get it?

Thx for answers.

Answer 1

You could find it via simulation:

In[187]:= 
sample = RandomVariate[
   CompoundPoissonDistribution[1, GammaDistribution[100, 200]], 10^7];

In[188]:= Quantile[sample, 0.95`]

Out[188]= 59877.1

Answer 2

You can approximate the Quantile using data generated with RandomVariate

dist = CompoundPoissonDistribution[1, GammaDistribution[100, 200]];

While you can calculate the Mean and StandardDeviation, PDF or CDF and by extension Quantile don't evaluate with this CompoundPoissonDistribution

{μ, σ} = {Mean[dist], StandardDeviation[dist] // N}

(*  {20000, 20099.8}  *)

PDF and CDF return unevaluated

{PDF[dist, 20000.], CDF[dist, 20000.]}

(*  {PDF[CompoundPoissonDistribution[1, GammaDistribution[100, 200]], 20000.], 
 CDF[CompoundPoissonDistribution[1, GammaDistribution[100, 200]], 20000.]}  *)

Approximating

SeedRandom[1];

data = RandomVariate[dist, 10000];

{μest, σest} = {Mean[data], StandardDeviation[data]}

(*  {19998.9, 19910.6}  *)

Comparing with the theoretical values

({μest, σest} - {μ, σ})/{μ, σ}

(*  {-0.0000548537, -0.00940959}  *)

The approximate Quantile is

q95est = Quantile[data, 0.95]

(*  59850.  *)

Answer 3

The OP seeks a symbolic/theoretical solution. To proceed symbolically, consider first what the OP is actually asking. The question is this:

The Question

Let $X \sim \text{Gamma}(a,b)$, and let $\{X_1, X_2,\dots, X_m\}$ denote an iid sample of size $m$, where the sample size $m$ (instead of being fixed) is itself a Poisson random variable $M=m$. The OP seeks the distribution of the sample sum:

$$Y = X_1 + X_2 + \dots + X_m \quad \quad \text{where} \quad M \sim \text{Poisson}(1)$$

As $M$ is Poisson, and the domain of support of a Poisson includes 0, it follows that the sample size $M$ can be 0, in which case $Y = 0$. This is important, because it means that $P(Y = 0)$ will have discrete mass.

The OP has sought to implement this model by using the syntax:

    CompoundPoissonDistribution[1, GammaDistribution[a, b]]

but unfortunately it does not seem to work with PDF or CDF etc

Solution

To proceed, first note that the sum of $m$ independent identical $\text{Gamma}(a,b)$ variables has a $\text{Gamma}(m a,b)$ distribution i.e. $Y$ has pdf $f(y \; \big| \; M = m)$:

_{(source: tri.org.au)}

where parameter $M \sim \text{Poisson}(1)$ with pmf $g(m)$:

_{(source: tri.org.au)}

We seek the parameter mixture distribution of $Y$ and $M$, which in standard Mathematica syntax, is:

PDF[ParameterMixtureDistribution[GammaDistribution[m a, b], 
  m \[Distributed] PoissonDistribution[1]], y]

Unfortunately, this too does not work, irrespective of whether one enters numbers for the parameters $a$ and $b$, or symbols, or even numerical values for $y$ -- Mma just whirrs or returns Undefined. So, let us try a different approach ...

Unconditional pdf of $Y$

• Discrete Part: $Y = 0$

$Y = 0$ iff $m = 0$. This occurs with probability $P(M=0)$:

_{(source: tri.org.au)}

• Continuous Part: The parameter-mix distribution, for $Y>0$, is given by:

_{(source: tri.org.au)}

where:

• I am using the Expect function from the mathStatica package for Mathematica

• The OP has specified some numerical values for $a$ and $b$ (e.g. $a = 100$ and $b = 200$). Mma can find a solution as a function of arbitrary $b$, so long as a numerical value for $a$ is set. Here, without loss of generality, we have set $a = 1$. It works just as well with $a=100$ ... the answer will just take longer to produce, and be more messy.

In summary, the unconditional pdf of $Y$ is:

$$\text{pdf}(Y) = \left\{ \begin{array}{cc} \frac{1}{e} & \text{ if } y = 0 \\ \text{sol} & \text{ if } y > 0 \\ \end{array}\right.$$

which is a mixed discrete-continuous distribution.

Quantiles

Mathematica cannot integrate the Hypergeometric0F1Regularized above, but we can still make our own numerical probability or cdf function (here, with say $b =10$):

NProb[w_?NumericQ] := 1/E + NIntegrate[sol /. b -> 10, {y, 0, w}]

For example, when $a = 1$ and $b=10$, the $P(Y\leq 20)$ is:

NProb[20]

0.817415

Then, the 0.95 quantile is found with:

FindRoot[NProb[y] == 0.95, {y, 10, 100}]

{y -> 39.1804}

In the OP's case, with $a = 100$ and $b=200$, the same approach yields:

FindRoot[NProb[y] == 0.95, {y, 20000, 100000}]

{y -> 59883.1}

Verify that it is correct:

NProb[59883.1]

0.95

And all is good.

Checking and Testing Simulation Answers posted by others
Just for fun, the two simulation answers produced (see other answers) return:

NProb[59850]

0.949766

NProb[59877.1]

0.949958