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?
$\begingroup$
$\endgroup$
2
1 Answer
$\begingroup$
$\endgroup$
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).
RandomVariate[]still doesn't play nice withProbabilityDistribution[]in the multivariate case, so you'll still need to put in some effort. $\endgroup$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. $\endgroup$