Simulation – How to Simulate a Joint Distribution Using the Inverse Method for Better Results

cumulative distribution functionquantilessimulation

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!

Best Answer

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

Figure

The agreement looks good.

Comment

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

Related Question