4
$\begingroup$

Let x[1], x[2], ..., x[n] be independent random variables each of which is distributed exponentially, i.e. with the PDF

f[x_,a_] := a Exp[-a x], a>0, x>=0

and let

t = Product[x[i],{i,1,n}]

Calculate the PDF and the CDF of t for arbitrary n = 1, 2, 3, ...

Remark: As I have indicated in my answer to Probability: Calculating a multiple integral this is an interesting problem which I have solved in Mathematica. I propose it here for others to find possible different paths to the solution.

$\endgroup$
5
$\begingroup$

The Problem

Let $(X_1, \dots, X_n)$ denote independent and identically distributed variables, each with common Exponential pdf $f(x)$:

f = a Exp[-a x];       domain[f] = {x, 0, Infinity}  &&  {a > 0}; 

We seek the pdf of $\prod_{i=1}^n X_i$, for $n = 2, 3, \dots$

Solution

The pdf of the product of two Exponentials is simply:


(source: tri.org.au)

where I am using the TransformProduct function from the mathStatica package for Mathematica.

The domain of support is:

 domain[g] = {x, 0, Infinity}  &&  {a > 0};

The product of 3, 4 and 5 Exponentials is obtained by iteratively applying the same function (here 3 times):


(source: tri.org.au)

By induction, the pdf of the product of $n$ iid Exponentials is thus:

$$a^n \text{MeijerG}[\{ \{ \}, \{ \} \}, \{ \{0_1, \dots, 0_n \}, \{ \} \}, a^n x] \quad \quad \text{ for } x > 0 $$

Quick Monte Carlo check: Here is a quick check comparing:

  • the theoretical pdf just obtained (when $n = 3$ and $a=2$): RED DASHED curve,
  • to the empirical Monte Carlo pdf based: squiggly BLUE curve


(source: tri.org.au)

And all is well.

CDF

In the same way, we can find the cdf at any iteration $n$. For the current value $n=5$, the cdf $P_g(X<x)$ is:


(source: tri.org.au)

... using mathStatica's Prob function. This takes the same form for each value of $n$. By induction, the general form for the cdf is:

MeijerG[{{1}, {}}, {lis, {0}}, a^n x] where lis is an $n$-length vector of 1's.

$\endgroup$
  • $\begingroup$ @ wolfies: very elegant solution in accordance with mine (+1). As I have version 8 I could not use TransformProduct. $\endgroup$ – Dr. Wolfgang Hintze Feb 9 '15 at 22:36
  • $\begingroup$ Hi Wolfgang. TransformProduct is a function from the mathStatica package for Mathematica ... i.e. a stats add-on to mma. According to the compatibility table at mathstatica.com/faq/compatibility.html, Mma 8 is compatible with mathStatica 2.5 ... which was also the first version of mathStatica to include the TransformProduct function. So it should work fine with Mma 8. It would also be interesting to see your solution using standard Mathematica functions, if you have a chance to post it. $\endgroup$ – wolfies Feb 10 '15 at 4:10
3
$\begingroup$

Although we have now the elegant solution by wolfies I shall post my solution which uses only standard functions of MMA version 8 and which might be of interest because of a strategy to overcome difficulties with integration.

We calculate the probability distribution functions PDF ($fp$) and CDF ($fc$) of a product of $n$ independent random variables $x_1, x_2, \text{...}, x_n$ distributed according to exponential distributions $a_i e^{-a_i x_i}$ for arbitrary integer $n > 0$.

First we note that it is sufficient to study random variables distributed with the same PDF and with the Parameter set to $a = 1$. In fact, the reader can easily show that for the CDF $fc$ with a set ${a_1, a_2, ..., a_n}$ of parameters we have

$fc\left(n,t,\{a_1,a_2, ..., a_n\}\right)=fc\left(n,t*{\prod _{i=1}^n a_i},\{1_1,1_2, ..., 1_n\}\right)$

Now let

f[x_]:=Exp[-x]

The CDF for n = 2 is defined by

fci[2, s_] = 
 Integrate[f[x] f[y] UnitStep[s - x y], {x, 0, \[Infinity]}, {y, 0, \[Infinity]}, Assumptions -> s > 0]

(* 1 - 2 Sqrt[s] BesselK[1, 2 Sqrt[s]] *)

For n = 3 we have

fci[3, s_] = 
 Integrate[f[x] f[y] f[z] UnitStep[s - x y z], {x, 0, \[Infinity]}, {y,0, \[Infinity]}, {z, 0, \[Infinity]}, Assumptions -> s > 0]

(* 
Integrate[
 E^-x (1 - 2 Sqrt[s/x] BesselK[1, 2 Sqrt[s/x]]), {x, 0, \[Infinity]}, 
 Assumptions -> s > 0] 
*)

Unfortunately, already for n = 3 Mathematica does not perform all the necessary integrations.

In order to make progress let's perform some "heuristic experimental mathematics" with Mathematica. To this end we shall first study fci[3,s] for some small integer values of s and look for patterns.

Table[{s, fci[3, s]}, {s, 1, 5}]

$\left\{\left\{1,1-\text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},1\right]\right\},\left\{2,1-\sqrt{2} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},2\right]\right\},\left\{3,1-\sqrt{3} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},3\right]\right\},\left\{4,1-2 \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},4\right]\right\},\left\{5,1-\sqrt{5} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},5\right]\right\}\right\}$

Here the special function MeijerG appears which is, fortunately, contained as a standard function in Mathematica.

The pattern is easily seen, and we "conclude" that for n = 3 we have

fcx[3, t_] :=  1 - Sqrt[t] MeijerG[{{}, {}}, {{-(1/2), 1/2, 1/2}, {}}, t]

Now let's increase n

fci[4, s_] := 
 Integrate[f[x] f[y] f[z] f[u] UnitStep[s - x y z u], {x, 0, \[Infinity]}, {y,
    0, \[Infinity]}, {z, 0, \[Infinity]}, {u, 0, \[Infinity]}, 
  Assumptions -> s > 0]

Again for some integer arguments

Table[fci[4, s], {s, 1, 5}]

$\left\{1-\text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},1\right],1-\sqrt{2} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},2\right],1-\sqrt{3} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},3\right],1-2 \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},4\right],1-\sqrt{5} \text{MeijerG}\left[\{\{\},\{\}\},\left\{\left\{-\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right\},\{\}\right\},5\right]\right\}$

We find very similar expressions to those for n = 3, only the list in the MeijerG-function is longer by one

fcx[4, s_] = 1 - Sqrt[s ] MeijerG[{{}, {}}, {{-(1/2), 1/2, 1/2, 1/2}, {}}, s]

We "conclude" that generally for integer n > 0 and real s > 0 the CDF is given by

fcxx[n_, s_] := 
 1 - Sqrt[s ] MeijerG[{{}, {}}, {Join[{-(1/2)}, Array[1/2 &, n - 1]], {}}, s]

Calling it leads to an automatic simplification:

fcxx[4, s]

(*
1 - MeijerG[{{}, {}}, {{0, 1, 1, 1}, {}}, s]
*)

Hence we have found good reasons to assume that the general CDF is given by

fc[n_, s_] := 1 - MeijerG[{{}, {}}, {Join[{0}, Array[1 &, n - 1]], {}}, s]

This is in fact the same as the results of wolfies

fc1[n_, s_] := MeijerG[{{1}, {}}, {Array[1 &, n], {0}}, s] (* wolfies, 07.02.15 *)

as can be shown by plotting both functions for some values of n.

The PDF is given by the derivative with respect to s

fpx[n_, s_] := D[fc[n, s], s]

The first few are

Table[fpx[n, s], {n, 1, 5}]

(*
{E^-s, BesselK[0, 2 Sqrt[s]] - BesselK[1, 2 Sqrt[s]]/Sqrt[s] + 
  BesselK[2, 2 Sqrt[s]], MeijerG[{{}, {}}, {{0, 0, 0}, {}}, s], 
 MeijerG[{{}, {}}, {{0, 0, 0, 0}, {}}, s], 
 MeijerG[{{}, {}}, {{0, 0, 0, 0, 0}, {}}, s]}
*)

We conclude for the PDF the expression

fp[n_, t_] := MeijerG[{{}, {}}, {Array[0 &, n], {}}, t]

This is exactly the same expression as that obtained by wolfies.

EDIT #1 Moments, Asymptotics, Complex Integral

Moments

Because the variables are independent we have for the k-th moment of t

$m\left(t^k\right) = m\left((x_1 x_2 ... x_n)^k\right) = m\left(x_1^k \right) m\left(x_2^k\right) ... m\left(x_n^k \right) = m\left(x^k \right)^n$

and because

$m\left(x^k \right) = \int_0^{\infty } \exp (-x) x^k \, dx = \Gamma (k+1) = k!$

the k-th moment of t is just $(k!)^n$.

This result is confirmed by direct calculation

Table[1/(k!)^2 Integrate[t^k fp[2, t], {t, 0, \[Infinity]}], {k, 1, 10}]

(*
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
*)

Table[1/(k!)^3 Integrate[t^k fp[3, t], {t, 0, \[Infinity]}], {k, 1, 10}]

(*
{1, 1, 1, 1, 1, 1, 1, 1, 1, 1}
*)

and so on for

n = 4, 5, ...

Asymptotics at infinity

For the first few n we have

Series[fp[1, t], {t, \[Infinity], 1}] // Normal

$e^{-t}$

Series[fp[2, t], {t, \[Infinity], 1}] // Normal

$e^{-2 \sqrt{t}} \left(\sqrt{\pi } \left(\frac{1}{t}\right)^{1/4}-\frac{1}{16} \sqrt{\pi } \left(\frac{1}{t}\right)^{3/4}\right)$

Series[fp[3, t], {t, \[Infinity], 1}] // Normal

$e^{-3 t^{1/3}} \left(\frac{4 \pi \left(\frac{1}{t}\right)^{1/3}}{\sqrt{3}}-\frac{4 \pi \left(\frac{1}{t}\right)^{2/3}}{9 \sqrt{3}}\right)$

Looking at the expressions up to, say n = 10, we can guess for the leading term the following expression

$\text{fp}(n,t\to \infty ) = \sqrt{\frac{2}{\pi n t}}\left(\sqrt{2\pi t}\right)^n \exp \left(-n t^{\frac{1}{n}}\right)$

We can see that for large n the exponential decay becomes very slow.

Asymptotics at zero

For the first few n we get

Series[fp[1, t], {t, 0, 1}] // Normal // TraditionalForm

$1-t$

Series[fp[2, t], {t, 0, 1}] // Normal // TraditionalForm

$t (-\log (t)-2 \gamma +2)-\log (t)-2 \gamma$

where $\gamma$ is the constant EulerGamma.

Series[fp[3, t], {t, 0, 1}] // Normal // TraditionalForm

$\frac{1}{4} t \left(-2 \log ^2(t)-12 \gamma \log (t)+12 \log (t)-\pi ^2-18 \gamma ^2+36 \gamma -24\right)+\frac{1}{4} \left(2 \log ^2(t)+12 \gamma \log (t)+\pi ^2+18 \gamma ^2\right)$

...

Hence the general leading term for n>1 is a power of $log(1/t)$:

$\text{fp}(n,t\to 0 ) = \frac{\log ^{n-1}\left(\frac{1}{t}\right)}{(n-1)!}$

MeijerG as a complex integral

Looking at the General deifition of this function in terms of a complex integral (from http://functions.wolfram.com/HypergeometricFunctions/MeijerG/02/0001/)

$\text{MeijerG}\left[\left\{\left\{a_1,a_n\right\},\left\{a_{1+n},a_p\right\}\right\},\left\{\left\{b_1,b_m\right\},\left\{b_{1+m},b_q\right\}\right\},z\right]==\frac{1}{2\pi i } \text{ContourIntegrate}\left[\frac{z^{-s} \left(\prod _{k=1}^n \text{Gamma}\left[1-s-a_k\right]\right) \prod _{k=1}^m \text{Gamma}\left[s+b_k\right]}{\left(\prod _{k=1+n}^p \text{Gamma}\left[s+a_k\right]\right) \prod _{k=1+m}^q \text{Gamma}\left[1-s-b_k\right]},\{s,\mathcal{L}\}\right]$

We find that all $a_i$ are missing, ${b_1, ... , b_m} = {0_1, ... , 0_n}$ and all other b's are missing so that

$fp(n,t)$ has the representation

$fp(n,t) = \frac{1}{2 \pi i}\int _{c - i \infty }^{c + i \infty }t^{-s} \Gamma (s)^nds$

where c is a real number >0. Changing the integration variable to s = 1 + I u we find

$fp(n,t) = \frac{1}{2 \pi }\int _{-\infty }^{+ \infty }t^{-(1+ i u)} \Gamma (1+ i u)^ndu$

The identity of this expression to the original one can be shown with examples using NIntegrate.

Observation: the expression for $fp(n,t)$ is nothing but the inverse Mellin transformation of the n-th power of the Mellin transform of the exponential Distribution, which is igiven by

$M( f(t), s) = \int_0^{\infty } f(t) t^{s-1} \, dt$

specifically

$\int_0^{\infty } \text{Exp}(-t) t^{s-1} \, dt = \Gamma (s)$

$\endgroup$
  • $\begingroup$ @ sysadmin: I don't see the Latex typewriting in the editor anymore. It seems that it is transformed only after saving the file. For example the input "$m\left(t^k\right)$" is left unchanged. How can I recover this feature? $\endgroup$ – Dr. Wolfgang Hintze Feb 12 '15 at 9:48
  • $\begingroup$ @Mr. Wizard. I'm ready. Please also have a look at my previous comment. $\endgroup$ – Dr. Wolfgang Hintze Feb 15 '15 at 21:51
  • $\begingroup$ Sorry, I am not sure what you mean. Are you saying that the preview area below the edit box does not show live LaTeX, but that it previously did? $\endgroup$ – Mr.Wizard Feb 15 '15 at 22:18
  • $\begingroup$ @Mr. Wizard, yes, exactly. $\endgroup$ – Dr. Wolfgang Hintze Feb 15 '15 at 23:16
  • $\begingroup$ Strange; I still get live preview. I suggest you post a report on Meta Stack Exchange as this is beyond my control. $\endgroup$ – Mr.Wizard Feb 15 '15 at 23:26

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.