# How to compute the marginal CDF of a joint density?

I am trying to compute the cumulative distribution function of a random variable $$u$$ that has the following density:

$$f(u) = \int_{1}^\infty \frac{e^{-4uv}}{v^5}dv$$

for $$u \gt 0$$.

### What I've tried

I've tried computing this integral, giving a function $$f(u)$$, and then calculating $$\int_{-\infty}^x f(u) du$$, obtaining the CDF of the density. But I always get some weird results from this integral. I know that there must be a better way to solve this using the fact that this is a joint density, but I don't know how to do it. Can someone please give me any hint or advice?

• What are your "weird results" ? I think your joint density is missing a multiplicative factor of 20 so that it integrates to 1.
– JimB
Commented May 11, 2022 at 20:43
• dlmf.nist.gov/8.19 gives you a good set of tools to evaluate this.
– whuber
Commented May 11, 2022 at 21:13
• @DilipSarwate Not sure who you're asking but it would appear to be $\frac{20 \exp (-4 u v)}{v^5}$ (once the correct multiplier is applied) and if $u$ ranges from $0$ to $\infty$.
– JimB
Commented May 12, 2022 at 2:38
• Note that $u$ must not be negative, for then the integral no longer exist. Moreover, you can simply excahnge the order of integrations, which leads to a different power in the denominator of the integrand. The resulting integral can be looked up on the page to which @whuber linked. Commented May 12, 2022 at 8:01

More generally, consider the distribution whose density is

$$f_p(u, \theta) = \theta C_p\int_1^\infty v^{-p} e^{-\theta u v}\,\mathrm{d}v = \theta C_pE_p(\theta u)$$

for $$u \gt 0$$ and parameters $$p\gt 0,$$ $$\theta \gt 0.$$ $$E_p$$ is the Exponential Integral function and $$C_p$$ is a normalizing constant (which we will find at the end). Notice, for future reference, that

$$E_{p+1}(0) = \int_1^\infty v^{-(p+1)}\,\mathrm{d}v = \frac{1}{p}.$$

Moreover, since for any $$z\gt 0$$

$$0\le E_p(z) = \int_1^\infty v^{-p} e^{-zv}\,\mathrm{d}v\le \int_1^\infty e^{-zv}\,\mathrm{d}v = \frac{1}{z},$$

it follows from the Squeeze Theorem (for limits) that

$$\lim_{z\to\infty} E_p(z) = 0.$$

The distribution function (cdf) is, by definition,

$$F_p(u,\theta) = \int_0^u f_p(z)\,\mathrm{d}z = C_p \int_0^u E_p(\theta z)\,\theta\mathrm{d}z = C_p \int_0^{\theta u} E_p(z)\,\mathrm{d}z.$$

The integrand defining $$E_p$$ is smooth and decreasing so quickly at its upper limit that we may interchange the operations of differentiating and integrating, showing

\begin{aligned} \frac{\mathrm{d}}{\mathrm{d}z}E_{p+1}(z) &= \frac{\mathrm{d}}{\mathrm{d}z}\int_1^\infty v^{-(p+1)} e^{-zv}\,\mathrm{d}v\\&=\int_1^\infty \frac{\mathrm{d}}{\mathrm{d}z}v^{-(p+1)} e^{-zv}\,\mathrm{d}v\\&=-\int_1^\infty v^{-p} e^{-zv}\,\mathrm{d}v\\&=-E_p(z). \end{aligned}

It is immediate (by the Fundamental Theorem of Calculus) that

$$F_p(u,\theta) = C_p \left[E_{p+1}(0) - E_{p+1}(\theta u)\right] = C_p \left[\frac{1}{p} - E_{p+1}(\theta u)\right].$$

Finally, since $$\lim_{u\to\infty}F_p(u,\theta)=1$$ (by the Law of Total Probability), we find

$$1 = \lim_{u\to\infty}C_p \left[\frac{1}{p} - E_{p+1}(\theta u)\right] = \frac{C_p}{p},$$

showing $$C_p = p$$ and giving

$$F_p(u,\theta) = 1 - pE_{p+1}(\theta u).$$

Set $$p=5$$ and $$\theta=4$$ for the answer to the question. Here are graphs of the density and its cdf:

These were drawn in gray using the analytical solution. As a check, over them are plotted, in red, the values obtained from numerically computing the single integral for $$f$$ and the double integral for $$F.$$ The R code that generated these graphs follows.

#
# Solution.
#
library(expint)
f <- function(u, p, theta) theta * p * expint_En(theta * u, p)
F <- function(u, p, theta) 1 - p * expint_En(theta * u, p+1)
#
# Brute force verification using numerical integration.
#
f. <- Vectorize(function(u, p, theta, ...) {
theta * p * integrate(function(v) exp(-theta * u * v) / v^p, 1, Inf, ...)$$value }, "u") F. <- Vectorize(function(u, p, theta, ...) { integrate(function(z) f.(z, p, theta), 0, u, ...)$$value
}, "u")
#
# Example.
#
p <- 5
theta <- 4
par(mfrow=c(1,2))
curve(f(x, p, theta), 0, 1.5, n=501, lwd=2, col=gray(0.25),
xlab="u", ylab="Density",
main=bquote(f[.(p)](u, .(theta))))
curve(f.(x, p, theta), add = TRUE, col="Red", lwd=2, lty=2)

curve(F(x, p, theta), 0, 1.5, n=501, lwd=2, col=gray(0.25),
yaxp=c(0,1,1),
xlab="u", ylab="Probability",
main=bquote(F[.(p)](u, .(theta))))
curve(F.(x, p, theta), add = TRUE, col="Red", lwd=2, lty=2)
par(mfrow=c(1,1))

• many thanks, so great when your answer is accompanied with fully replicable r code! +1 Commented May 12, 2022 at 18:12
• @Max I just hope my redefinition of F doesn't cause any problems ;-).
– whuber
Commented May 12, 2022 at 19:00

Warning: this answer/approach lacks mathematical rigor and relies on Mathematica to do the heavy lifting.

From looking at what is given for the marginal density for $$U$$

$$f(u)=\int_1^\infty \frac{e^{-4 u v}}{v^5} dv$$

one can surmise that the joint pdf for $$U$$ and $$V$$ is proportional to $$\frac{e^{-4 u v}}{v^5}$$. I say "proportional" because the integral over the values of $$U$$ and $$V$$ integrates to 1/20 and not to 1. Therefore a multiplier of 20 is needed.

Integrate[Exp[-4 u v]/v^5, {v, 1, \[Infinity]}, {u, 0, \[Infinity]}]
(* 1/20 *)


So the joint pdf is $$\frac{20 e^{-4 u v}}{v^5}$$. Integrating over $$v$$ gets the pdf for $$U$$:

Integrate[20 Exp[-4 u v]/v^5, {v, 1, \[Infinity]}, Assumptions -> u > 0]


$$\frac{5}{3} \left(e^{-4 u} \left(3-4 u \left(8 u^2-2 u+1\right)\right)-128 u^4 \text{Ei}(-4 u)\right)$$

where $$\text{Ei}$$ is the exponential integral function. The cdf is then found with

Integrate[5/3 (E^(-4 u) (3 - 4 u (1 - 2 u + 8 u^2)) - 128 u^4 ExpIntegralEi[-4 u]), {u, 0, z},
Assumptions -> z > 0] /. z -> u


$$\frac{1}{3} \left(-128 u^5 \text{Ei}(-4 u)+e^{-4 u} \left(u \left(3-4 u \left(8 u^2-2 u+1\right)\right)-3\right)+3\right)$$

Here are plots of the pdf and cdf:

Update

I followed @whuber 's suggestion below to treat the pdf symbolically using the general form in his answer to see if Mathematica would then produce a simpler form for the pdf and cdf. And it did.

pdf = \[Theta] p Integrate[v^(-p) Exp[-\[Theta] u v], {v, 1, \[Infinity]},
Assumptions -> u > 0 && \[Theta] > 0 && p > 0] // FunctionExpand


$$\theta p (\theta u)^{p-1} \Gamma (1-p,u \theta )$$

cdf = Simplify[Integrate[pdf, {u, 0, w},
Assumptions -> \[Theta] > 0 && p > 0 && w > 0] /. w -> u // FunctionExpand,
Assumptions -> u > 0 && \[Theta] > 0 && p > 0]


$$(\theta u)^p \Gamma (1-p,u \theta )-e^{\theta (-u)}+1$$

where $$\Gamma (1-p,u \theta )$$ is the incomplete gamma function.

• These formulas are correct and agree with my solution. They arise because it appears Mathematica doesn't know how to simplify the general Exponential Integral function. Perhaps applying FullSimplify would help clean up the formula.
– whuber
Commented May 12, 2022 at 17:02
• @whuber. Thanks. I did try FullSimplify (with Assumptions -> u>0) and FunctionExpand and got no simplification.
– JimB
Commented May 12, 2022 at 17:04
• Did you try this using the constants $5$ and $4$ as stated in the question or by treating them as symbolic variables? In the latter case, if Mathematica can evaluate the integral it almost certainly will give it in a simpler form (perhaps as a confluent hypergeometric function, incomplete Gamma function, or Exponential Integral). (I no longer have access to Mathematica or I would try this myself...)
– whuber
Commented May 12, 2022 at 17:10
• @whuber I have not tried that but I will shortly. In the meantime there is the "free" version of Mathematica called "Wolfram Engine for Developers". Essentially it is "free" and "legal" as long as you don't make money using it. (Of course, the license is more specific than that but I think it's essentially a similar license as "SAS on demand for academics".)
– JimB
Commented May 12, 2022 at 17:20
• Thank you for that tip! FullSimplify isn't actually needed. wolframscript -code "Assuming[u > 0, Integrate[Exp[-4 u x]/x^p, {x,1,Infinity}]]" returns ExpIntegralE[p, 4*u], as one would hope. Following that up with an evaluation of "Integrate[ExpIntegralE[5, 4*u], {u, 0, y}]" gives (1 - 5*ExpIntegralE[6, 4*y])/20.
– whuber
Commented May 12, 2022 at 18:58