6
$\begingroup$

I wish to obtain numbers from

a = 1.245710578482853*p^2*r*Sqrt[(16 + 32*p^2 - 112*p^4 + 64*p^6 - 3*r^2)/(4 - 8*p^2 + 4*p^4 - r^2)];
RandomVariate[ProbabilityDistribution[Piecewise[{{a, 0 <= r <= 1 - p^2}}, 0], {p, 0, 1}, {r, 0, 1}]]

but receive the message

Sampling from ProbabilityDistribution ... is not implemented.>>

This sort of problem was discussed in 2635 a few years ago. Do betters solutions exist now? Thanks.

$\endgroup$
2
$\begingroup$

Here is an approximate answer which takes advantage of the marginal distribution of r being very close to a Beta distribution.

First the joint density function is normalized:

a = p^2*r*Sqrt[(16 + 32*p^2 - 112*p^4 + 64*p^6 - 3*r^2)/(4 - 8*p^2 + 4*p^4 - r^2)];
const = Integrate[a, {p, 0, 1}, {r, 0, 1 - p^2}];
jointDensity = Piecewise[{{a/const, 0 <= r <= 1 - p^2}}, 0];

Now plot the joint density function:

ContourPlot[jointDensity, {p, 0, 1}, {r, 0, 1}, PlotRange -> All,
 Contours -> {0.1, 0.5, 1, 2, 3, 4, 5, 6, 7, 8},
 PlotRangePadding -> None,
 Frame -> True, FrameLabel -> (Style[#, Bold, Large] &) /@ {"p", "r"}]

Contour plot of joint density

Plot3D[jointDensity, {p, 0, 1}, {r, 0, 1}, PlotRange -> All,
 ImageSize -> Large,
 AxesLabel -> (Style[#, Bold, Larger] & ) /@ {"p", "r", "Probability density"}]

3D plot of joint density

Find the marginal distribution of r and also find a Beta distribution that has the same mean and variance:

μr = NIntegrate[r jointDensity, {r, 0, 1}, {p, 0, Sqrt[1 - r]}]
(* 0.42410218810335665 *)
σ2r = NIntegrate[(r - μr)^2 jointDensity, {r, 0, 1}, {p, 0, Sqrt[1 - r]}]
(* 0.04306834618787338 *)
solr = Solve[{μr == ar/(ar + br), σ2r == ar br/((ar + br)^2 (ar + br + 1))}, {ar, br}][[1]]
(* {ar\[Rule]1.9809707931759324,br\[Rule]2.6900043839036916} *)
Plot[{PDF[BetaDistribution[ar, br], r] /. solr, 
  NIntegrate[jointDensity, {p, 0, Sqrt[1 - r]}]}, {r, 0, 1}]

Marginal density of r and corresponding Beta distribution

One can see that there is very little difference in the two density functions. (The marginal density for p is also very similar to a Beta distribution but not as similar to a Beta distribution as r.)

So to take (approximate) random samples from the joint distribution of p and r we take a random sample from the approximate marginal distribution of r followed by the conditional distribution of p given r.

n = 1000;
rMarginal = BetaDistribution[ar, br] /. solr[[1]];
rSample = RandomVariate[rMarginal, n];
pSample = 
  Table[RandomVariate[
    ProbabilityDistribution[
       (jointDensity /. r -> rSample[[i]])/PDF[rMarginal, rSample[[i]]],
       {p, 0, Sqrt[1 - rSample[[i]]]}], 1][[1]], {i, n}];
ListPlot[Transpose[{pSample, rSample}], AspectRatio -> 1, 
 Frame -> True,
 FrameLabel -> (Style[#, Larger, Bold] &) /@ {"p", "r"},
 PlotRange -> {{0, 1}, {0, 1}}]

Display of 1000 random samples from joint distribution

$\endgroup$

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.