4
$\begingroup$

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 at 15:12
    • $\begingroup$ @whuber if it converges how do you proof $\int_0^\infty f_{U_i}(y) dy=1$? $\endgroup$ – Aymen Kareem Jan 18 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 at 20:18
    • $\begingroup$ If it's not a CDF then was my derivation wrong? $\endgroup$ – Aymen Kareem Jan 18 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 at 20:38
    5
    +50
    $\begingroup$

    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

    $$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)}$$

    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}.$$

    Therefore

    $$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).

    $\endgroup$
    • $\begingroup$ Great!!! Thank you very much. $\endgroup$ – Aymen Kareem Jan 21 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$ – Aymen Kareem Feb 2 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 at 21:21
    • $\begingroup$ please have a look at my new question. I need $E[x^2] $also $\endgroup$ – Aymen Kareem Feb 2 at 21: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.