I have the following joint distribution:

$$f(x, y) = 3x^2y^xe^{-x^3}(1 + x),\quad x \gt 0,\ y \in (0,1).$$

I want to simulate a sample of this distribution through the inverse method but I don't know if my procedure is right.

What I've tried

In order to simulate the random variable above I'll have to compute the inverse marginal cdf for $x, y$ , $F^{-1}_X(x), F^{-1}_Y(x)$ and then for $U, V \sim U_{[0, 1]}$ I'll be able to simulate the random variable. But I don't know how to compute the marginal cdf's in this case. Any hints or advice are more than welcome!

  • $\begingroup$ Hint: you will likely have to integrate out either $x$ or $y$ to find a marginal distribution. Which one leads to a (much) easier integral? Here's another hint: $$f(x,y)=\left(\frac{\mathrm{d}}{\mathrm{d}x}\left(-e^{-x^3}\right)\right)\,\left(\frac{\mathrm{d}}{\mathrm{d}y}y^{1+x}\right).$$ $\endgroup$
    – whuber
    May 14, 2022 at 21:17
  • 1
    $\begingroup$ Once again, @whuber, you've helped me! When i find one of the marginal distribuitions I'll be able to find the other one because $f_Y(y | x) = \frac{f(x,y)}{f_X(x)}$. And then, with both marginals I can get both CDF's, and I'll just have to invert then in order to properly simulate this r.v, right? $\endgroup$
    – Occhima
    May 14, 2022 at 22:06
  • 3
    $\begingroup$ You can't generate from the two marginals to get a value from the joint; the variables are not independent. Instead you use the general decomposition into marginal and conditional ... $f_{X,Y}(x,y) = f_X(x) . f_{Y|X}(y|x)$ $\endgroup$
    – Glen_b
    May 15, 2022 at 1:04
  • $\begingroup$ Note also that a naïve implementation of the inverse cdf method in dimensions larger than one does not work. $\endgroup$
    – Xi'an
    May 15, 2022 at 8:25

1 Answer 1


For the record, integrating out $y$ gives the marginal density

$$f_X(x) = \int_0^1 f(x,y)\,\mathrm{d}y = 3x^2e^{-x^3}.$$

By inspection (or integration via a substitution for $x^3$) this has distribution function

$$F_X(x) = 1 - e^{-x^3}\quad (x \ge 0).$$

Conditional on $x,$ the density of $y$ is proportional to $y^x,$ whose distribution function is

$$F_{Y\mid X}(y\mid x) = y^{1+x}\quad (0\le y \le 1).$$

Thus, if we draw a uniform variable $U$ in the range $[0,1],$ $x = F_X^{-1}(U)$ will have the same distribution as $X;$ and then if we independently draw another uniform variable $V,$ solving the equation

$$V = F_{Y\mid X}(y\mid x)$$

for $y$ will give us one realization $(x,y)$ from the distribution with density $f.$

A practical example is this R implementation, which draws n iid values from this distribution:

rjoint <- function(n) {
  x <- (-log(runif(n)))^(1/3)
  y <- runif(n)^(1/(1+x))
  cbind(x, y)

(It saves a little time by applying $F_X^{-1}$ to $1-U,$ which also has a uniform distribution.)

The return value is an $n\times 2$ array with $x$ in the first column and $y$ in the second.

The result of drawing ten million values with x <- rjoint(1e7) (which takes two seconds on this machine) looks like this figure, where contours of the resulting 2D density are drawn over a colored contour plot of $f$ (both use the same contour intervals ranging from $0.2$ near the bottom to $2.2$ at the top).


The agreement looks good.


This problem had two obvious approaches: integrate out $y$ or integrate out $x.$ Both work, but the latter is so nasty I don't know how to do it by hand. The Wolfram Language 13 engine found an expression:

Figure 2

However, it doesn't know how to simplify or invert this function ;-).

  • 1
    $\begingroup$ I just got a free dev' license this year for the Wolfram Language. I should use it more, especially for problems where Wolfram alpha times out. Anyway, wikipedia gives some elementary definitions for the Airy-related functions in the Wolfram language output. $\endgroup$
    – Galen
    May 14, 2022 at 23:36

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.