Generate two random numbers with constraints

Question

I want to generate two random numbers, $p$ and $q$, between $0.5$ and $1$.

They are connected by the constraint $1/(2q) > p$.

How do I generate $p$ and $q$?

Answer 1

Assuming you just ommitted the parentheses and meant $1/(2q) > p$, then this works

q = RandomReal[{0.5, 1}]
p = RandomReal[{0.5, 1/(2 q)}]

Edit: The above code will generate 2 random numbers that satisfy the given criteria, they won't be sampling the same distribution.

qlist = RandomReal[{0.5, 1.0}, 10000];
plist = RandomReal[{0.5, 1/(2 #)}] & /@ qlist; 
ListPlot[ Transpose[{qlist, plist}]]

You see that the first variable samples {0.5,1.0} uniformly while the other does not. But the inequality is symmetric, and can be written $1/(2p)>q$ so that isn't right. rhermans's answer fixes this, but it runs unreasonably slow on my machine. For example, generating a list of 1000 pairs takes about 100 seconds for me using ImplicitRegion and RandomPoint.

But I can do it in about a tenth of a second using this inelegant code

list2 = Reap[
     i = 0;
     While[
      i < 10001,
      test = RandomReal[{0.5, 1.0}, 2];
      If[test[[1]] < 1/(test[[2]] 2), i++; Sow[test]];
      ]
     ][[2, 1]]; // AbsoluteTiming
ListPlot@list2

Answer 2

This assumes uniform distribution. See answer by @JimBaldwin for discussion on limitations (implicit assumptions) of my answer.

Answer

region = ImplicitRegion[0.5 < q < 1. && 0.5 < p < 0.5/q, {p, q}];
RandomPoint[region]
(* {0.793318, 0.550934} *)

---

Visual

Show[
 RegionPlot[region]
 , ListPlot[RandomPoint[region, 1000], PlotStyle -> Red]
 ]

---

Documentation

We are using ImplicitRegion

And RandomPoint

Answer 3

The question does not state an essential piece of information which is the joint distribution of $p$ and $q$. All of the previous answers (so far) jump to a solution without making the joint distribution explicit (at least prior to what one sees in the code and the resulting figures).

The answers using regions assume that $p$ and $q$ have uniform distributions, $U(0.5,1)$, but they are restricted to the region $0.5 \le p \le 1$, $0.5 \le q \le 1$, and ${1\over{2q}}>p$.

The other answer given is that $p\sim U(0.5,1)$ and $q|p\sim U(0.5,{1\over{2p}}))$.

There is no reason why one couldn't have $p\sim 0.5+0.5\,\mathrm{Beta}(\alpha,\beta)$ with $q|p\sim U(0.5,{1\over{2p}})$. A random sample from a linear function of a beta random variable is also a legitimate random sample.

The question most likely is about restricting two independent random variables from uniform distributions to a region, but that should be made explicit.

Addition: Doing it the hard way

Using regions and RandomPoint is the way to go as @rhermans describes (especially if the region of interest is not simple). But if you want to go about it in a brute force way, here is an option.

First we assume that without the additional restrictions that $p$ and $q$ are independently distributed on $U(0.5,1)$. The joint probability density function in the square of interest is

f[p_, q_] := 4

Now determine the joint probability density when $0.5\le p<1/(2q)\le 1$:

c = Integrate[f[p, q], {q, 1/2, 1}, {p, 1/2, 1/(2 q)}];
g[p_, q_] := f[p, q]/c
(* 4/(-1 + Log[4]) *)

Find the marginal probability density function for $p$:

gp[p_] := FullSimplify[Integrate[g[p, q], {q, 1/2, 1/(2 p)}]]
(* 2(1-p)/(p*(-1+Log[4])) *)

Now find the conditional distribution of $q$ given $p$ (which is just a uniform distribution on $0.5$ to $1/(2p)$):

gqGivenp[q_, p_] := g[p, q]/gp[p]
(* 2p/(1-p) *)

Define a ProbabilityDistribution for $p$:

dp = ProbabilityDistribution[gp[p], {p, 1/2, 1}];

Finally, generate a set of bivariate random samples:

n = 5000;  (* Number of samples *)
rp = RandomVariate[dp, n];
rq = 0.5 + RandomReal[1, n]*(1/(2 rp) - 0.5);

Plotting the results we have

ListPlot[Transpose[{rp, rq}], PlotRange -> {{0.5, 1}, {0.5, 1}}, AspectRatio -> 1]

This approach might also be successful when the original joint probability function is a bit more complex.

Answer 4

Perhaps more a comment: note differences:

r = RandomReal[{0.5, 1}, {10000, 2}];
Show[ListPlot[Sort@GatherBy[r, Times @@ # < 0.5 &], 
  PlotStyle -> {{Red}, Blue}], 
 Plot[1/(2 x), {x, 0.5, 1}, PlotStyle -> Green], Frame -> True]

compared with:

ListPlot[{#, RandomReal[{0.5, 1/(2 #)}]} & /@ 
  RandomReal[{0.5, 1}, 10000], PlotStyle -> Red, Frame -> True]

Desired outcome depends on aim. Latter understandable given narrowing of uniform distribution of "q" as "p"->1

Answer 5

Here's another approach with RandomVariate and ProbabilityDistribution. First we find the area that can be occupied by points satisfying the criteria:

norm = Integrate[Boole[1/(2 y) > x], {x, 1/2, 1}, {y, 1/2, 1}]
(* (Log[4] - 1)/4 *)

Let p run along the horizontal axis and q along the vertical axis.

RegionPlot[ImplicitRegion[1 > 1/(2 p) > q > 1/2, {p, q}], PlotRange -> {0, 1}]

Here's the region of allowed values for {p,q}.

equal small elements of area from this region must have equal likelihood of being occupied. The area corresponding to a range $(p, p + dp)$ is then simply $\frac{dp}{2 p} - \frac{dp}{2}$. Dividing by $dp$ we'll get the probability distribution of p. We define a distribution:

dist = ProbabilityDistribution[(1/(2 x) - 1/2)/norm, {x, 1/2, 1}]

which must be normalized (/norm).

Plot[PDF[dist, x], {x, 2/5, 6/5}]

Now we generate a list of random ps with the above distribution:

(p = RandomVariate[dist, 100000]) // AbsoluteTiming // First
(* 0.211 seconds *)

and find a random q corresponding to the constraint of each of the random p:

(q = RandomReal[{.5, 1/(2 #)}] & /@ p) // AbsoluteTiming // First
(* ~ 9 msec *)

Then plot every 100th point (not all of them, otherwise it's difficult to assess if it's uniform or not):

ListPlot[Transpose[{p, q}][[;; ;; 100]]]

This is about 5 times faster than JasonsReap-Sow` approach.

Most likely, better results can be achieved by using a proper combination of built-in distribution functions. Unfortunately, RandomVariate refuses to accept multivariate distributions from ProbabilityDistribution, which is why there's a generation of first a list of ps and then a list of qs, instead of directly generating pairs.

Answer 6

Probably not the most efficient way but

FindInstance[1/(2 q) > p && 0.5 < p < 1 && 0.5 < q < 1, {p, q}, 100, 
RandomSeed -> RandomInteger[100]]