6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Aug 10 '16 at 14:00
  • $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful $\endgroup$ – Michael E2 Aug 10 '16 at 14:00
  • $\begingroup$ Thx for reaction. No, it is not critical. $\endgroup$ – Jonas Vojtech Aug 10 '16 at 16:38
6
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Thx for answer. And is there any way hot to calculate exact quantiles of Poisson-Gamma distribution using Mathematic? $\endgroup$ – Jonas Vojtech Aug 10 '16 at 16:40
6
$\begingroup$

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.  *)
$\endgroup$
  • $\begingroup$ Thx for your answer. And is there any way hot to calculate exact quantiles of Poisson-Gamma distribution using Mathematic? Why I can not use Quantile[] on CompoundPoissonDistribution[], if I am able to use Mean[], Variance[] and other things... $\endgroup$ – Jonas Vojtech Aug 10 '16 at 16:42
  • $\begingroup$ @JonasVojtech - PDF, CDF, and Quantile are more complicated calculations than Mean or StandardDeviation` and Mathematica doesn't have the rules/transformations/algorithms built in to evaluate the more complicated calculations--assuming that they can even be done in closed form. $\endgroup$ – Bob Hanlon Aug 10 '16 at 17:06
6
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Thx for your answer. Is there any way how to get the mathStatica package for free or how to replace the function Expect? Thank you. $\endgroup$ – Jonas Vojtech Dec 10 '16 at 22:41
  • 1
    $\begingroup$ @JonasVojtech If you don't have a copy of mathStatica, please email me, and I'd be happy to send you a copy. $\endgroup$ – wolfies Dec 11 '16 at 6:36
  • $\begingroup$ Hello, I am sorry but I can not find your email. Could you send me your email here or info about the copy to jonas.vojtech@seznam.cz . Thank you! I tried to replace function Expect[], I used Expectation[f, m-PoissonDistribution[1]], but it is really time consuming. I hope correct. $\endgroup$ – Jonas Vojtech Dec 11 '16 at 21:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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