Random number generation of a Multivariate PDF

Question

I want to simulate random data with a multivariate PDF. I see that five years ago this needed to be directly coded into mathematica directly, but it seems as though there is recent support for multivariate distributions. The only documentation is for distributions which are predefined. I would like to come up with my own multivariate PDF and do simulations of the data. (For example $f(x, y) = (x+y^2)e^{-(x^2+y^2)}$. Can Mathematica handle such a task more easily now?

Answer 1

If your question is about the example joint density function, then here is a direct approach.

(* Find constant of integration and define joint density function *)
c = Integrate[(x + y^2) Exp[-(x^2 + y^2)], {x, 0, ∞}, {y, -∞, ∞}]
(* 1/4 (2 Sqrt[π]+π) *)
fxy = (x + y^2) Exp[-(x^2 + y^2)]/c

Now find marginal distribution for $y$ followed by the conditional distribution of $x$ given $y$:

(* Marginal density for y *)
fy = Integrate[fxy, {x, 0, ∞}];
(* Conditional density of x given y *)
fxGiveny = fxy/fy;

(* Marginal distribution for y *)
disty = ProbabilityDistribution[fy, {y, -∞, ∞}];
(* Conditional distribution of x given y *)
distxGiveny = ProbabilityDistribution[fxy/fy, {x, 0, ∞}];

Now we can take a random sample of size $n$ from the bivariate distribution:

n = 100;
ySample = RandomVariate[disty, n];
(* Random sample of x given y *)
xSample = RandomVariate[distxGiveny /. y -> #, 1][[1]] & /@ ySample;
(* Random sample from bivariate distribution *)
data = Transpose[{xSample, ySample}]

(This is just a special case of the original post).