31
$\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$
  • 7
    $\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 '15 at 12:53
  • $\begingroup$ @Szabolcs mathematica.stackexchange.com/help/badges/57/reversal :D $\endgroup$ – Dr. belisarius Oct 20 '15 at 13:14
  • 2
    $\begingroup$ To add to the comment by @Szabolcs, this question is far from trivial. $\endgroup$ – Daniel Lichtblau Oct 20 '15 at 14:43
  • 2
    $\begingroup$ This is relevant here too: en.wikipedia.org/wiki/Bertrand_paradox_(probability) $\endgroup$ – Szabolcs Oct 20 '15 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$ – David G. Stork Aug 11 '17 at 22:21
16
$\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$
  • $\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$ – Chandan Datta Oct 20 '15 at 12:01
  • 2
    $\begingroup$ in what way is it not working? $\endgroup$ – Jason B. Oct 20 '15 at 12:04
  • $\begingroup$ @ChandanDatta qs = RandomReal[{0.5, 1}, 100];ps = RandomReal[{0.5, 1/(2 #)}] & /@ qs; $\endgroup$ – Dr. belisarius Oct 20 '15 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 '15 at 12:13
  • $\begingroup$ Thanks to both of you. Now Its working fine. $\endgroup$ – Chandan Datta Oct 20 '15 at 12:17
33
$\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$
  • 2
    $\begingroup$ @rhermans neat and thank you for introducing me to RandomPoint +1 :) $\endgroup$ – ubpdqn Oct 20 '15 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 '15 at 12:51
  • 1
    $\begingroup$ That does seem like a bug, I'm using linux version 10.2 $\endgroup$ – Jason B. Oct 20 '15 at 13:06
  • 1
    $\begingroup$ The overhead is mainly the symbolic computations, numerical approximation reduces it significantly. See here $\endgroup$ – rhermans Oct 20 '15 at 15:32
27
$\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$
  • 3
    $\begingroup$ +1; people should really get into the habit of thinking about distributions when dealing with random number/point generation problems. $\endgroup$ – J. M. is away Oct 20 '15 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$ – Jacob Akkerboom Oct 20 '15 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 '15 at 14:49
  • 3
    $\begingroup$ I agree. But I think any answers should make that assumption explicit. $\endgroup$ – JimB Oct 22 '15 at 17:36
  • 1
    $\begingroup$ @JimBaldwin thank you for this instructive answer +1. :) $\endgroup$ – ubpdqn Oct 26 '15 at 5:19
18
$\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$
  • $\begingroup$ Also campare with ListPlot[{RandomReal[{0.5, 1/(2 #)}], #} & /@ RandomReal[{0.5, 1}, 10000], PlotStyle -> Red, Frame -> True]. $\endgroup$ – Karsten 7. Oct 20 '15 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 '15 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 '15 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$ – Karsten 7. Oct 20 '15 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 '15 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, privacy policy and cookie policy

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