Probability distribution function expressed in terms of a divergent series

Question

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$

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

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:

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