5

I can generate a 2D random number within a square region, say $\{x,y\} \in [p_i,p_f]$, using

rand = RandomReal[{pi, pf}, 2];

I now would like to get a random number in this region but outside the sphere centered around a and with radius r with a+r<pf and a-r>pi.

In this thread various ways of generating random numbers inside a region are discussed. However, I couldn't figure out how to modify the responses to solve my problem.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Shasa
  • 1,043
  • 5
  • 12
  • 3
    Use RandomPoint with an appropriate region (look up RegionDifference). Obviously, the region you use must be of finite size. If this did not work, be specific about where you got stuck. – Szabolcs Jun 09 '22 at 17:45
  • @Szabolcs I can use RandomPoint[Disk[], 1] to generate a random number inside a Disk but how can I generate a number outside this Disk? – Shasa Jun 09 '22 at 17:50
  • The outside of a disk is infinite. It makes no sense to ask for a random point within an infinite region, at least not with uniform distribution. – Szabolcs Jun 09 '22 at 17:51
  • I already mentioned in my question that I would like to have my numbers to be inside a square! Anyway, I got an answer. – Shasa Jun 09 '22 at 17:55
  • 1
    That is not at all clear from the phrasing. – Szabolcs Jun 09 '22 at 17:57
  • 2
    @Szabolcs After double-reading my question can understand your point. I have modified the question for future readers. – Shasa Jun 09 '22 at 18:16

3 Answers3

13

You can use RandomPoint with RegionDifference:

reg = RegionDifference[Rectangle[{-1, -1}, {1, 1}], Disk[{0,0}, 1]];
pts = RandomPoint[reg, 10]

{{-0.865556, -0.700941}, {0.913053, -0.737421}, {0.288003, 0.962705}, {0.772499, 0.856204}, {0.956926, -0.539727}, {0.29716, -0.99774}, {-0.544408, 0.968495}, {-0.78035, 0.657943}, {-0.600353, 0.836259}, {0.705759, -0.835619}}

Visualization:

Graphics[{LightBlue, DiscretizeRegion@reg, PointSize[Large], Red, Point[pts]}]

enter image description here

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
2

Here is a slightly faster method:

samplePointC = Compile[{},
  Module[{pt = {0., 0.}},
   While[
    pt = RandomReal[{-1, 1}, 2];
    pt.pt < 1];
   pt
   ]
  ]

Which gives the following timing:

AbsoluteTiming[
 pts = Table[
    samplePointC[],
    {i, 1, 100000}];
 ]

>>> {0.0693943, Null}

The RandomPoint method gave ~ 0.986045 seconds. For those who are interested, I also tried reflecting points that are inside the diamond given by $|x|+|y|=1$ to reject about half as many points but this only caused a minor slow down. I used the following code before pt.pt<1 to do this:

u = pt[[1]] + pt[[2]];
v = pt[[1]] - pt[[2]];
If[Abs[u] > Abs[v],
  If[Abs[u] < 1, u = Sign[u] 2 - u],
  If[Abs[v] < 1, v = Sign[v] 2 - v]
  ];
pt[[1]] = 1/2 (u + v);
pt[[2]] = 1/2 (u - v);
AccidentalTaylorExpansion
  • 1,267
  • 6
  • 10
2

Here is a brute-force rejection method. The speed depends mainly on the desired sample size and somewhat on the relative size of the circle within the square.

(* Set parameters *)
pi = 1;
pf = 4;
r = 0.66;
x0 = (pf + pi)/2 + r;
y0 = (pf + pi)/2 - r;

(* Relative area outside of circle *)
p = 1 - π (r/(pf - pi))^2;


(* Desired sample size *)
m = 100000;


(* Minimum sample size needed to be relatively certain to have at least m observations outside of the circle *)
n = (18 + m - 18 p)/p + 20 Sqrt[(9 + m - 18 p - m p + 9 p^2)/p^2] // Round;


(* Generate random sample excluding area in circle )
SeedRandom[12345];
AbsoluteTiming[xy = RandomVariate[UniformDistribution[{{pi, pf}, {pi, pf}}], n];
  xy = Select[xy, (#[[1]] - x0)^2 + (#[[2]] - y0)^2 > r^2 &][[1 ;; m]];]
( {0.394952, Null} *)


ListPlot[xy, AspectRatio -> 1, PlotRange -> {{pi, pf}, {pi, pf}}]

Uniform random sample excluding area in circle

JimB
  • 41,653
  • 3
  • 48
  • 106