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?

Answer 1

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}}  *)