33
$\begingroup$

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

$\endgroup$
5
  • 8
    $\begingroup$ 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? $\endgroup$
    – Szabolcs
    Oct 20, 2015 at 12:53
  • $\begingroup$ @Szabolcs mathematica.stackexchange.com/help/badges/57/reversal :D $\endgroup$ Oct 20, 2015 at 13:14
  • 2
    $\begingroup$ To add to the comment by @Szabolcs, this question is far from trivial. $\endgroup$ Oct 20, 2015 at 14:43
  • 2
    $\begingroup$ This is relevant here too: en.wikipedia.org/wiki/Bertrand_paradox_(probability) $\endgroup$
    – Szabolcs
    Oct 20, 2015 at 19:15
  • $\begingroup$ 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. $\endgroup$ Aug 11, 2017 at 22:21

6 Answers 6

17
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Oct 20, 2015 at 12:01
  • 2
    $\begingroup$ in what way is it not working? $\endgroup$
    – Jason B.
    Oct 20, 2015 at 12:04
  • $\begingroup$ @ChandanDatta qs = RandomReal[{0.5, 1}, 100];ps = RandomReal[{0.5, 1/(2 #)}] & /@ qs; $\endgroup$ Oct 20, 2015 at 12:10
  • 1
    $\begingroup$ And then just to verify the results: (0.5 < #1 < 1.0 && 0.5 < #2 < 1.0 && #2 < 1/(2 #1)) & @@@ Transpose[{qs, ps}] $\endgroup$
    – Jason B.
    Oct 20, 2015 at 12:13
  • $\begingroup$ Thanks to both of you. Now Its working fine. $\endgroup$ Oct 20, 2015 at 12:17
35
$\begingroup$

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

Mathematica graphics


Documentation

We are using ImplicitRegion

Mathematica graphics

And RandomPoint

Mathematica graphics

$\endgroup$
4
  • 2
    $\begingroup$ @rhermans neat and thank you for introducing me to RandomPoint +1 :) $\endgroup$
    – ubpdqn
    Oct 20, 2015 at 12:47
  • 1
    $\begingroup$ Is this unreasonably slow for anyone else? When I evaluate list = Table[RandomPoint[region], {1000}]; // AbsoluteTiming it takes about 100 seconds $\endgroup$
    – Jason B.
    Oct 20, 2015 at 12:51
  • 1
    $\begingroup$ That does seem like a bug, I'm using linux version 10.2 $\endgroup$
    – Jason B.
    Oct 20, 2015 at 13:06
  • 1
    $\begingroup$ The overhead is mainly the symbolic computations, numerical approximation reduces it significantly. See here $\endgroup$
    – rhermans
    Oct 20, 2015 at 15:32
28
$\begingroup$

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]

bivariate sample

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

$\endgroup$
6
  • 4
    $\begingroup$ +1; people should really get into the habit of thinking about distributions when dealing with random number/point generation problems. $\endgroup$ Oct 20, 2015 at 15:54
  • 1
    $\begingroup$ 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. $\endgroup$ Oct 20, 2015 at 18:55
  • $\begingroup$ Its not that unreasonable to assume uniform distribution as default if no other information is given. Arguably it is the simplest distribution. $\endgroup$
    – rhermans
    Oct 22, 2015 at 14:49
  • 3
    $\begingroup$ I agree. But I think any answers should make that assumption explicit. $\endgroup$
    – JimB
    Oct 22, 2015 at 17:36
  • 1
    $\begingroup$ @JimBaldwin thank you for this instructive answer +1. :) $\endgroup$
    – ubpdqn
    Oct 26, 2015 at 5:19
19
$\begingroup$

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]

enter image description here

compared with:

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

enter image description here

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

$\endgroup$
9
  • $\begingroup$ Also campare with ListPlot[{RandomReal[{0.5, 1/(2 #)}], #} & /@ RandomReal[{0.5, 1}, 10000], PlotStyle -> Red, Frame -> True]. $\endgroup$
    – Karsten7
    Oct 20, 2015 at 12:38
  • $\begingroup$ @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... $\endgroup$
    – ubpdqn
    Oct 20, 2015 at 12:43
  • $\begingroup$ @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' $\endgroup$
    – ubpdqn
    Oct 20, 2015 at 12:55
  • 1
    $\begingroup$ You had my +1 before I added my comment. In my opinion pointing out the differences between these two approaches is essential. $\endgroup$
    – Karsten7
    Oct 20, 2015 at 12:56
  • 2
    $\begingroup$ 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. $\endgroup$
    – march
    Oct 20, 2015 at 22:20
12
$\begingroup$

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

region

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

pdf

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

points

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.

$\endgroup$
0
$\begingroup$

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]]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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