# Generate two random numbers with constraints

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$?

• Come on people, if an answer was worth 6 upvotes, surely the question must be worth a few too, no matter how trivial. Many did learn right here about RandomPoint, no? – Szabolcs Oct 20 '15 at 12:53
• – Dr. belisarius Oct 20 '15 at 13:14
• To add to the comment by @Szabolcs, this question is far from trivial. – Daniel Lichtblau Oct 20 '15 at 14:43
• This is relevant here too: en.wikipedia.org/wiki/Bertrand_paradox_(probability) – Szabolcs Oct 20 '15 at 19:15
• I think a much more interesting question would be to find points $(p,q)$ obeying the constraint but that $P(p) \sim {\cal U}[.5, 1]$ and $P(q) \sim {\cal U}[.5, 1]$, which none of these solutions ensure. – David G. Stork Aug 11 '17 at 22:21

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/(test[] 2), i++; Sow[test]];
]
][[2, 1]]; // AbsoluteTiming
ListPlot@list2 • I have tried that one. Its working for one single random number. But If I want to produce say 100 numbers then its not working. – Chandan Datta Oct 20 '15 at 12:01
• in what way is it not working? – Jason B. Oct 20 '15 at 12:04
• @ChandanDatta qs = RandomReal[{0.5, 1}, 100];ps = RandomReal[{0.5, 1/(2 #)}] & /@ qs; – Dr. belisarius Oct 20 '15 at 12:10
• And then just to verify the results: (0.5 < #1 < 1.0 && 0.5 < #2 < 1.0 && #2 < 1/(2 #1)) & @@@ Transpose[{qs, ps}] – Jason B. Oct 20 '15 at 12:13
• Thanks to both of you. Now Its working fine. – Chandan Datta Oct 20 '15 at 12:17

This assumes uniform distribution. See answer by @JimBaldwin for discussion on limitations (implicit assumptions) of my 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  • @rhermans neat and thank you for introducing me to RandomPoint +1 :) – ubpdqn Oct 20 '15 at 12:47
• Is this unreasonably slow for anyone else? When I evaluate list = Table[RandomPoint[region], {1000}]; // AbsoluteTiming it takes about 100 seconds – Jason B. Oct 20 '15 at 12:51
• That does seem like a bug, I'm using linux version 10.2 – Jason B. Oct 20 '15 at 13:06
• The overhead is mainly the symbolic computations, numerical approximation reduces it significantly. See here – rhermans Oct 20 '15 at 15:32

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


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


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.

• +1; people should really get into the habit of thinking about distributions when dealing with random number/point generation problems. – J. M. will be back soon Oct 20 '15 at 15:54
• I agree with you and on this particular subject I have seen quite some ignorance, possibly resulting from a lack of formal language. On the other hand, it appears some people will use the word "random" in isolation to mean uniformly distributed quite consistently. – Jacob Akkerboom Oct 20 '15 at 18:55
• Its not that unreasonable to assume uniform distribution as default if no other information is given. Arguably it is the simplest distribution. – rhermans Oct 22 '15 at 14:49
• I agree. But I think any answers should make that assumption explicit. – JimB Oct 22 '15 at 17:36
• @JimBaldwin thank you for this instructive answer +1. :) – ubpdqn Oct 26 '15 at 5:19

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

• Also campare with ListPlot[{RandomReal[{0.5, 1/(2 #)}], #} & /@ RandomReal[{0.5, 1}, 10000], PlotStyle -> Red, Frame -> True]. – Karsten 7. Oct 20 '15 at 12:38
• @Karsten7. Yes, if you think I am being unhelpful I am happy to delete...I note 2 down votes so I may have misunderstood something... – ubpdqn Oct 20 '15 at 12:43
• @Karsten7. By yes I mean appreciate the symmetry and perhaps I should have used aspect ratio 1 and amplified point showing dense points on vertical axis...to illustrate what was my point which is care about 'random' – ubpdqn Oct 20 '15 at 12:55
• You had my +1 before I added my comment. In my opinion pointing out the differences between these two approaches is essential. – Karsten 7. Oct 20 '15 at 12:56
• One of my favorite "teaching" moments when teaching computational physics is to ask the students to generate random points in a circle, using polar coordinates; invariably, they uniformly choose the radius of the circle and are confused when the density of points decreases with the radius. – march Oct 20 '15 at 22:20

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 - 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.

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]
`