I'm interested in finding the CDF and PDF of $U_i$ defined as follows, $$U_i=\frac g{d^{\alpha}}$$ where $g$ is a gamma distributed random variable with shape $k$ and scale $\theta$, and $d$ is a random variable with distribution $f_d(x)=\frac {2x}{R^2}, 0\le x \le R$.
To find the CDF of $U_i$ $$F_{U_i}(y)=P(U_i<y)=\int_0^R F_g(yx^\alpha)f_d(x)dx$$ where $F_g(y)$ is the CDF of gamma distributed random variable.
After performing the integration, which requires the use of the expansion of incomplete gamma function, I came to the following CDF: $$F_{U_i}(y)=\sum_{n=0}^\infty \frac {2(-1)^nR^{k\alpha+n\alpha}y^{k+n}}{\Gamma(k)n!\theta^{k+n}(k+n)(k\alpha+n\alpha+2)}$$ I want the CDF in this format because my ultimate goal is to find the CDF of $U=\sum_{i=1}^NU_i$ which requires me to find the Laplace transform of $F_{U_i}(y)$ and raise it to power $N$.
Now it's easy to find the PDF of $U_i$, i.e. $f_{U_i}(y)=\frac {d}{dy}F_{U_i}(y)$.
I'm looking for someone to clarify these points to me:

  • My derived $f_{U_i}(y)$ doesn't make sense to me, because it's based on a divergent series. How can it's integration from zero to infinity equals one?
  • In Matlab I took the $K^{th}$ elements of the series, and plotted it along with the histogram of 1000 $U_i$'s, the two plots are not related.
  • If you look at the definition of $U_i$, $U_i$ can take values from 0 to infinity, high values take less probability. This fact is not clear in the derived distribution $f_{U_i}(y)$. What is the domain of the derived PDF? What is the integration limit if I wannna find $F_\epsilon(\epsilon)=\int f(y)F_{U_i}(\epsilon(y+\sigma^2))dy$? where $f(y)$ is another density function.
  • How can I simulate the derived PDF in Matlab to check if it's really correlated with $U_i$
  • $\endgroup$
    • $\begingroup$ Please explain why you think this series diverges. As far as I can tell, for most values of $\alpha$ (including all positive values) it is an entire function: it converges everywhere in the complex plane. $\endgroup$
      – whuber
      Jan 18, 2020 at 15:12
    • $\begingroup$ @whuber if it converges how do you proof $\int_0^\infty f_{U_i}(y) dy=1$? $\endgroup$
      – William
      Jan 18, 2020 at 20:04
    • $\begingroup$ That's irrelevant, because a more basic problem is that $F_{U_i}$ obviously is not a CDF, since it is directly proportional to $\theta^{-k}.$ Thus, there is at most one $\theta$ that could possibly make this a valid CDF and for arbitrary $\theta$ and $k$ it cannot be a CDF. $\endgroup$
      – whuber
      Jan 18, 2020 at 20:18
    • $\begingroup$ If it's not a CDF then was my derivation wrong? $\endgroup$
      – William
      Jan 18, 2020 at 20:29
    • $\begingroup$ It must have been. There are some simple checks. One of them is that because $\theta$ is a scale parameter for $g,$ it must be a scale parameter for $U_i,$ whence the CDF must be a function of $y/\theta$--but it is not. Another is to plug in simple values of parameters. E.g., $\alpha=1,$ $k=1,$ $R=1$ gives a sum that evaluates to $(y^2 + 2e^{-y}(y+1)-2)/y^2,$ which works, so you might be close. $\endgroup$
      – whuber
      Jan 18, 2020 at 20:38

    1 Answer 1


    The series converges to a distribution function. It can be evaluated in closed form.

    Upon identifying the terms varying with $n,$ write your function in a simpler form as

    $$\begin{equation*} F_{U_i}(y)=\frac{2(R^{\alpha}y/\theta)^k}{\Gamma(k)}\sum_{n=0}^\infty \frac {(-R^\alpha y / \theta)^n}{n!(k+n)((k+n)\alpha+2)} \\ = \frac{2x^k}{\Gamma(k)}\sum_{n=0}^\infty \frac{(-x)^n}{n!\,(k+n)((k+n)\alpha+2)} \end{equation*}$$

    for $x = R^\alpha y /\theta.$ In other words, $\sigma=\theta/R^\alpha$ is a scale parameter.

    Assuming (from the form of the expression $g/d^\alpha$) that $\alpha \gt 0,$ every one of the terms in the sum is bounded above in size by $x^n / n!$ showing the sum is dominated by the absolutely convergent series for $\exp(x),$ whence the sum converges absolutely for all $x.$

    To evaluate such a sum we will use partial fractions. Consider the simpler function

    $$h(x, a) = \sum_{n=0}^\infty \frac{(-1)^n x^n}{n! (a+n)}$$

    with $a \ge 0.$ Similar considerations demonstrate absolute convergence so we may differentiate term by term to obtain

    $$\frac{\mathrm{d}}{\mathrm{d}x} \left(x^a h(x,a)\right) = \sum_{n=0}^\infty \frac{(-1)^n x^{n+a-1}}{n!} = x^{a-1}\sum_{n=0}^\infty \frac{(-1)^n x^{n}}{n!} = x^{a-1} e^{-x}.$$


    $$h(x,a) = x^{-a} \int^x t^{a-1} e^{-t}\mathrm{d}t = C + x^{-a}\,\gamma(a, x),$$

    where $\gamma$ is the lower Incomplete Gamma Function and $C$ is a constant of integration. It can be found by applying L'Hopital's Rule to $x^{-a}\int_0^x t^{a-1} e^{-t}\mathrm{d}t$ as $x\to 0^+$ and comparing the resulting limit of $C+1/a$ to $h(0,a) = (-0)^0/(0!(a+0)) = 1/a$ to conclude $C=0.$

    Now because

    $$\frac{1}{(k+n)((k+n)\alpha+2)} = \frac{1}{\alpha}\frac{1}{(k+n)(k+n+2/\alpha)} = \frac{1}{2}\left(\frac{1}{k+n} - \frac{1}{k+n+2/\alpha}\right),$$

    we again exploit the absolute convergence of the series for $F_{U_i}$ and set $a=k,$ $b=k+2/\alpha$ to express it as

    $$F_{U_i}(y) = \frac{2x^k}{\Gamma(k)}\sum_{n=0}^\infty \frac{(-x)^n}{n!\,(k+n)((k+n)\alpha+2)} = \frac{\gamma(k,x) - x^{-\frac{2}{\alpha}}\gamma(k+\frac{2}{\alpha},x)}{\Gamma(k)}$$ where $$x = y/\sigma = R^\alpha y / \theta \text{ and }\sigma=\theta/R^\alpha.$$

    Differentiating this to obtain the density is straightforward, yielding

    $$\frac{\mathrm{d}}{\mathrm{d}x} F^\prime_{U_i}(y) = \frac{2\gamma(k+2/\alpha, x)}{\alpha \Gamma(k) x^{1 + 2/\alpha}} \ge 0$$

    and it is also elementary to establish that

    $$\lim_{y\to\infty} F_{U_i}(y) = 1 \text{ and } \lim_{y\to 0^+} F_{U_i}(y) = 0.$$

    Therefore $F_{U_i}$ is the CDF of a continuous random variable.

    Almost any statistical computing platform will compute the CDF of Gamma variables (or, equivalently, of Chi-squared variables). This is the normalized version of $\gamma.$ For instance, here are R implementations of $F_{U_i}$ and its derivative:

    # CDF
    pFU <- function(x, k, alpha, scale=1) {
      x <- x / scale
      h <- function(x, a) 
        ifelse(x > 0, exp(pgamma(x, a, log.p=TRUE) + lgamma(a)), 1/a)
      (h(x, k) - x^(-2/alpha) * h(x, k + 2/alpha)) / gamma(k)
    # PDF
    dFU <- function(x, k, alpha, scale=1) {
      x <- x / scale
      ifelse(x <= 0, 0, 
             2/alpha * exp(-(1 + 2/alpha)*log(x) + pgamma(x, k+2/alpha, log.p=TRUE) + 
                             lgamma(k+2/alpha) - lgamma(k))) / scale

    Here are plots of them using the curve function:

    Figure showing graphs of the CDF and PDF

    The underlying red curves use these functions. As a check, overplotted in black are direct implementations of the original series (for the CDF) and a numerical derivative of pFU (for the PDF).

    • $\begingroup$ Great!!! Thank you very much. $\endgroup$
      – William
      Jan 21, 2020 at 2:50
    • $\begingroup$ How do you integrate the derived pdf from 0 to infinity and you get 1? How can I calculate the expected value? $\endgroup$
      – William
      Feb 2, 2020 at 20:31
    • $\begingroup$ (1) $F_{U_i}$ already is the integral of the PDF, so just take its limit as the argument goes to $\infty.$ (2) Compute the expected value by integrating $1-F_{U_i}(x)$ from $0$ to $\infty.$ Because $F_{U_i}$ has been expressed in terms of $\gamma,$ this looks is fairly simple. I would integrate by parts. $\endgroup$
      – whuber
      Feb 2, 2020 at 21:21
    • $\begingroup$ please have a look at my new question. I need $E[x^2] $also $\endgroup$
      – William
      Feb 2, 2020 at 21:26

    Your Answer

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

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