RandomVariate cannot sample from a custom multivariate ProbabilityDistribution

Question

Let us try to define 2D distribution through PDF:

f[x_, y_] := 1/((0.05 + (-0.79 + x)^2) (0.05 + (-0.89 + y)^2)) + 1/((0.05 + (-0.40 + x)^2) (0.05 + (-0.88 + y)^2)) + 1/((0.05 + (-0.66 + x)^2) (0.05 + (-0.43 + y)^2))

The plot of PDF looks good:

Plot3D[f[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> Full, ImageSize -> 400]

Now define a distribution (with "Normalize" method, as f(x,y) is not a normalized PDF)

ProDis = ProbabilityDistribution[f[x, y], {x, 0, 1}, {y, 0, 1}, Method -> "Normalize"]

and try to generate a sample

RandomVariate[ProDis, 10]

Mathematica returns the error:

RandomVariate::noimp: Sampling from ProbabilityDistribution[0.00387923 (1/((0.05 +Plus[<<2>>]^2) (0.05 +Plus[<<2>>]^2))+1/((0.05 +Plus[<<2>>]^2) (0.05 +Plus[<<2>>]^2))+1/((0.05 +Plus[<<2>>]^2) (0.05 +Plus[<<2>>]^2))),{[FormalX]1,0,1},{[FormalX]2,0,1}] is not implemented.

How can I generate samples from my distribution?

Likely related: Hot to generate a RandomVariate of a custom distribution. Unfortunately, the answer might be that you are out of luck with the built in methods that may not be capable of handling your custom distribution. Could yours be expressed as a mixture, truncation, or other combination of existing distributions? You may have more luck then. — MarcoB, May 31 '20 at 18:08
Yes, you may have to roll a rejection sampling method yourself; for more than one dimension, RandomVariate[] has poor support for arbitrary ProbabilityDistribution[] objects. — J. M.'s missing motivation, May 31 '20 at 18:40

score 10 · Accepted Answer · answered May 31 '20 at 21:48

Given your specific example, nothing is "wrong" with ProbabilityDistribution. @MarcoB and @J.M.'stechnicaldifficulties have given the issue: ProbabilityDistributon is more limited with multivariate distributions it knows and your bivariate distribution is not a standard or common distribution.

One can take random samples by finding the marginal distribution for one of the random variables and then finding the condition distribution for the other.

f[x_, y_] := (1/((0.05 + (-0.79 + x)^2) (0.05 + (-0.89 + y)^2)) + 
   1/((0.05 + (-0.40 + x)^2) (0.05 + (-0.88 + y)^2)) + 
   1/((0.05 + (-0.66 + x)^2) (0.05 + (-0.43 + y)^2)))

(* Constant of integration *)
cxy = Integrate[f[x, y], {x, 0, 1}, {y, 0, 1}]
(* 257.783 *)

(* Marginal density for x *)
fx[x_] := Integrate[f[x, y], {y, 0, 1}]/cxy

(* (Conditional) distribution of y given x *)
yGivenx[x_] := ProbabilityDistribution[(f[x, y]/cxy)/fx[x], {y, 0, 1}]

(* Take random sample from the bivariate distribution *)
n = 20;
SeedRandom[12345];
xx = RandomVariate[ProbabilityDistribution[fx[x], {x, 0, 1}], n];
yy = RandomVariate[yGivenx[#], 1][[1]] & /@ xx;
data = Transpose[{xx, yy}]
(* {{0.784308, 0.811649}, {0.544349, 0.482533}, {0.374593, 0.668343}, 
    {0.56832, 0.588582}, {0.752874, 0.660399}, {0.732304, 0.382063},
    {0.79024, 0.856047}, {0.239503, 0.778425}, {0.506675, 0.365711}, 
    {0.48217, 0.801524}, {0.625505, 0.36676}, {0.747901, 0.551866}, 
    {0.736869, 0.330566}, {0.626223, 0.267467}, {0.383792, 0.684872},
    {0.655435, 0.579795}, {0.292907, 0.573538}, {0.683136, 0.682127}, 
    {0.75874, 0.345907}, {0.22253, 0.71164}}  *)

Thanks, JibB. Surely, I'm aware of this method which is given in each primary course of stochastic modeling. However, I supposed that RandomVariate function includes this method and it is enough formally substitute any distribution in order to get a sample. But, obviously, RandomVariate is limited in its capability. — Konstantin, Jun 01 '20 at 06:26
You might want to send your request for more general multivariate distributions to Wolfram, Inc. Especially if Maple or MATLAB offers this capability. — JimB, Jun 01 '20 at 15:31

RandomVariate cannot sample from a custom multivariate ProbabilityDistribution

1 Answers1

Linked