4
$\begingroup$

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]

Plot

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?

$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – MarcoB May 31 at 18:08
  • 2
    $\begingroup$ Yes, you may have to roll a rejection sampling method yourself; for more than one dimension, RandomVariate[] has poor support for arbitrary ProbabilityDistribution[] objects. $\endgroup$ – J. M.'s technical difficulties May 31 at 18:40
9
$\begingroup$

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}}  *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Konstantin Jun 1 at 6:26
  • 1
    $\begingroup$ You might want to send your request for more general multivariate distributions to Wolfram, Inc. Especially if Maple or MATLAB offers this capability. $\endgroup$ – JimB Jun 1 at 15:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.