# Random number generation of a Multivariate PDF

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?

• Last I checked RandomVariate[] still doesn't play nice with ProbabilityDistribution[] in the multivariate case, so you'll still need to put in some effort. – J. M.'s technical difficulties Oct 3 '17 at 5:30
• Your comment on the 5-year old post is essentially "Does it still work?" I'd say "Yes". If there's something that doesn't work, then you should post a question with a specific joint density which would allow us to pinpoint the issue. My guess is that the problem is because the joint density above isn't normalized and needs to be multiplied by $\frac{4}{2 \sqrt{\pi }+\pi }$ or the Method="Normalize" option is necessary. I would also characterize the approaches in the old post as "direct" rather than "crudely coded". Having secretive black boxes is not always a good thing. – JimB Oct 4 '17 at 16:19

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][] & /@ ySample;
(* Random sample from bivariate distribution *)
data = Transpose[{xSample, ySample}]


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