# Simulating a joint distribution with the inverse method

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!

• 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).$$
– whuber
May 14 at 21:17
• 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? May 14 at 22:06
• 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)$ May 15 at 1:04
• Note also that a naïve implementation of the inverse cdf method in dimensions larger than one does not work. May 15 at 8:25

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.

### 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:

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

• 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. May 14 at 23:36