Generate random points in a region with higher probability in a sub-region

I am trying to generate random samples of points in the unit square in such a way that points inside a given circle are twice as likely to be selected as points outside of that circle.

I tried to approach the problem by creating a 2D distribution. To do that I need to transform the function

f[x_,y_]:= If[ (x-.5)^2 + (y-.5)^2 < .1, 1, .5]

into a probability distribution for both x and y between 0 and 1. Then I would use the PDF to generate points and eventually plot them. My main struggle comes from turning the If function into a two-variable PDF.

• What have you tried so far? Are you asking to plot this? If yes, what plotting function did you try and what difficulties did you have? – Szabolcs Jan 16 '14 at 21:53
• Listplot is the plotting function but the problem comes from making the If function be a probability distribution for both the x and the y coordinates of the points. I could only figure out how to do a PDF for each coordinate separately before. – Ryan Summers Jan 16 '14 at 21:56

You want to scale $f$ so that $\int_0^1\int_0^1\, f = 1$.

f[x_,y_]:=If[(x-.5)^2+(y-.5)^2<.1,1,.5]
k=Integrate[f[x,y],{x,0,1},{y,0,1}] (* pdf is f/k *)
Plot3D[f[x,y]/k,{x,0,1},{y,0,1},
PlotPoints->100,PlotRange->{0,Full}] Edit

Alright, what you want can be stated in a simpler way : random sample of points on the unit square with double probability in the circle. No need for pdf's. Just generate random points on the square. If a point lies in the circle, add it. If not, toss a coin, add the point only if you get Head. Here is a way to do that.

n = 0;
list = {};
r := RandomReal[];
coin := RandomChoice[{True, False}];

While[n < 10000, {x, y} = {r, r};
If[((x - .5)^2 + (y - .5)^2 <= .1 || coin), AppendTo[list, {x, y}]; n++];]

(a variant of rejection sampling. If you are very unlucky it will take a very long time -- in practice no need to worry).

Now look at it:

Show[
ListPlot[list, PlotStyle -> {PointSize[.005]}, AspectRatio -> 1],
Graphics[Circle[{.5, .5}, .1^(1/2)]]] Variants

• You can vary the density ratio between the circle and outside. If you want the circle to be, say, 5 times more likely, redefine coin := RandomReal[] < 1/5;

• You can replace the circle by any region as long as you have a membership function, just use that function in place of (x - .5)^2 + (y - .5)^2 <= .1.

• I like that but what I need now is a 2D listplot where there are twice as many points inside the circle because the function is twice as much there. – Ryan Summers Jan 16 '14 at 22:27
• Your question is not clear. Would "put two points at the center and one on a corner" do? If not, why? – A.G. Jan 16 '14 at 22:46
• Sorry I've learned I'm not good at clarifying lol. Ok so say we're going to plot 1000 points. Take one of those points in particular. The probability of it having a location (x,y) in the listplot has to be given by the If function of x and y. Since the If function has twice the value inside the circle, any given point has twice the probability of falling inside the circle. – Ryan Summers Jan 16 '14 at 22:51
• So there will be a jump in point density along the boundary of the circle, an edge if you will. – Ryan Summers Jan 16 '14 at 22:52
• @Ryan You should look at the documentation how to make a random variable from a pdf and then how to generate random samples from a random variable. – A.G. Jan 16 '14 at 23:05

I'm arriving a bit late to the party, but I have an alternative approach that could be of interest to some readers.

All of the other solutions proposed (as of this writing) use variations on the acceptance/rejection technique of sampling. This works well enough for your specific problem, where the difference in the intensity of the process between the areas inside/outside the circle is comparatively small. However, as others have already pointed out, efficiency can become a problem in the limit of large differences between the "background" and the "features" in your problem -- you could have to generate a large number of "rejects" for each "good" sample.

However, if we think of your problem as an example of a (nonhomogeneous) Poisson spatial point process, we can use the additivity property of the Poisson distribution to assemble a sample from your complicated intensity function by combining realizations of simpler, easier to generate, point processes on the square and the circle in isolation, with little or no "wasted" computational effort.

Before we start, let's build a couple of helper routines to generate single, uniformly distributed, points on general rectangles and circles:

generateUniformDistributedPointOnRectangle[{{left_, bottom_}, {right_, top_}}] :=
{RandomVariate[UniformDistribution[{left, right}]],
RandomVariate[UniformDistribution[{bottom, top}]]}

generateUniformDistributedPointOnCircle[{{centreX_, centreY_}, radius_}] :=
Module[{u, theta, r},
u = Total[RandomVariate[UniformDistribution[{0, radius}], {2}]];
theta = RandomReal[{0, 2 Pi}];
r = If[u > radius, 2 radius - u, u];
{centreX, centreY} + {r Cos[theta], r Sin[theta]}
]

Our plan of attack is now to generate $N_{S}$ samples from the entire square and $N_{C}$ samples from just the circular area such that when we combine the two samples, we are left with $N$ samples from our target distribution. So, how to determine the correct values (and distributions) for $N_{S}$ and $N_{C}$?

Our intensity function $\lambda(x,y)$ for the process is given by $$\lambda(x,y)=\begin{cases}\lambda_{0} & (x,y)\in S \setminus C \\ \lambda_{1} & (x,y) \in C \end{cases},$$ where $S$ and $C$ denote the square and the circle, respectively. The expected number of points in the pattern is then given by $\mu = \lambda_{0}(1-\pi / 10) + \pi\lambda_{1}/10$; the actual number of points in a given realization will be a stochastic variable distributed as $N \sim \mathrm{Poisson}(\mu)$.

Now, additivity tells us that if we have two Poisson distributed variables $X_{1}$ and $X_{2}$ with means $\mu_{1}$ and $\mu_{2}$, respectively, then their sum $X = X_{1}+X_{2}$ is also Poisson distributed with mean $\mu_{1} + \mu_{2}$. Moreover, if we condition on a particular total $x$, then we have $$X_{1} \sim \mathrm{Binomial}\left(x,\frac{\mu_{1}}{\mu_{1}+\mu_{2}}\right),$$ and analogously for $X_{2}$. We will use this result to determine how many points to sample from the square and how many to sample from the circle.

To this end, let us decompose the intensity function into the sum of a function $\lambda_{S}(x,y)=\lambda_{0}$ defined on the entire square and a function $$\lambda_{C}(x,y)=\begin{cases}\lambda_{1}-\lambda_{0} & (x,y) \in C \\ 0 & \text{otherwise}\end{cases}.$$ Therefore, if we want a total of $N$ points in our realization, we should draw a random integer distributed as $$N_{S} \sim \mathrm{Binomial}\left(N,\frac{10}{10+\pi(\lambda_{1}/\lambda_{0} -1)}\right),$$ for the number of points to draw from the uniform distribution over the entire square, and then augment that with $N-N_{s}$ points from the uniform distribution over the circle. Note that the probability only depends on the ratio of the intensities between the "background" and the "feature", and not their absolute magnitudes. (A consequence of conditioning on a particular number of points in the pattern.)

So what to use for the ratio? Here things get a little fuzzy, as it is not entirely clear what you mean by points are "twice as likely" to be drawn from the circle as from the rest of the area. @A.G. and @Rojo have interpreted this to be a direct statement on the ratio of the intensities of the point process within these two areas, i.e., $\lambda_{1}/\lambda_{0} = 2$. On the other hand, @ubpdqn has interpreted this as meaning that 2/3 of the points fall within the circle, and your later post seems to favour this interpretation. These are not equivalent statements, as in the latter case, the binomial probability becomes $10/(30-3\pi)$.

So, let's code up both options (in the same function, no less):

generatePointSetConditionedOnNumberOfPoints[num_Integer,
opts : OptionsPattern[{"Interpretation" -> "TwoThirds"}]] :=
Module[{numBackground},
Which[OptionValue["Interpretation"] == "TwiceRate",
numBackground = RandomVariate[BinomialDistribution[num, 10/ (10 + Pi)]]
,
True,
numBackground = RandomVariate[BinomialDistribution[num, 10/(3 (10 - Pi))]]
];
Join[Table[generateUniformDistributedPointOnRectangle[{{0, 0}, {1, 1}}], {numBackground}],
Table[generateUniformDistributedPointOnCircle[{{1/2, 1/2}, Sqrt[1/10]}],
{num - numBackground}]]
]

In the above function, specifying the option "Interpretation"->"TwiceRate" gives the pattern corresponding to the interpretation used by @A.G. and @Rojo, while the (default option) "Interpretation"->"TwoThirds" corresponds to the choice by @ubpdqn.

Let's take a look at the difference on realizations consisting of 1000 points generated using each method:

pointSet1 =
generatePointSetConditionedOnNumberOfPoints[1000, "Interpretation" -> "TwiceRate"];
pointSet2 =
generatePointSetConditionedOnNumberOfPoints[1000, "Interpretation" -> "TwoThirds"];

GraphicsRow[
Map[ListPlot[#, Axes -> False, AspectRatio -> 1,
Epilog -> {EdgeForm[Black], FaceForm[None], Rectangle[],
Directive[GrayLevel[0.2], Dashed],
Circle[{0.5, 0.5}, Sqrt[1/10]]
},
PlotRange -> {{-0.025, 1.025}, {-0.025, 1.025}}] &,
{pointSet1, pointSet2}], Scaled[0.02], ImageSize -> 600] The "twice the rate" interpretation is on the left, while the "two thirds of the points fall in the circle" interpretation is on the right. There is a fairly substantial difference between the two; in particular, the edge between the two regions is much harder to discern on the left than on the right. (Grey dashed circle added to help in identification.)

I apologize for the length; hopefully those brave souls who read the whole post found it useful and/or learned something along the way. And now ... off to bed.

@Rojo has provided the backbone for this answer (and I have voted 1+ for his answer). I post this as my interpretation of the question is the generating a plot of a square region (0,1)x (0,1) in which there is sampling from a mixed distribution such that being within a given circle was twice as likely as not within this region. The annular region of the randomBoth leads to overinclusion of points inside the circle. This can be tested by making an emprical distribution of the sample and testing the probability.

Modifying the code of @rojo for what I interpret as the request:

randomInCircle[r_] :=
NestWhile[RandomReal[r {-1, 1}, 2] &, Null, ! TrueQ[Norm[#] < r] &];
rout[r_] :=
NestWhile[RandomReal[0.5 {-1, 1}, 2] &,
Null, ! TrueQ[Norm[#] >= r] &];
rnd[r_] :=
If[RandomChoice[{w1, w2} -> {1, 2}] == 1, randomInCircle[r],
rout[r]] + {c1, c2}

Then plottting:

ListPlot[emp = Table[rnd[Sqrt[0.1]], {5000}],
AspectRatio -> Automatic]

yields: Using the sample as the basis for an empirical distribution:

empd = EmpiricalDistribution[emp];
NProbability[(x - .5)^2 + (y - .5)^2 <= 0.1, {x, y} \[Distributed]
empd]

yields:

0.6682, i.e. approximately 2/3 as desired.

I note that if the aim was to do probability calculations you could generate a mixed probability distribution as follows. It is, however, cumbersome.

p = ProbabilityDistribution[
Boole[(x - .5)^2 + (y - .5)^2 < 0.1]/(Pi 0.1), {x, 0, 1}, {y, 0, 1}]
q = ProbabilityDistribution[
Boole[(x - .5)^2 + (y - .5)^2 >= 0.1]/(1 - (Pi 0.1)), {x, 0, 1}, {y,
0, 1}]
mix = MixtureDistribution[{2, 1}, {p, q}]

Visualizing the PDF:

Plot3D[Evaluate[PDF[mix, {x, y}]], {x, 0, 1}, {y, 0, 1},
PlotRange -> {0, 3}, Exclusions -> None, PlotPoints -> 200] The CDF is a little more problematic but NProbability is useful:

ListPlot3D[
Table[NProbability[i < x && j < y, {i, j} \[Distributed] mix], {x, 0,
1, 0.1}, {y, 0, 1, 0.1}],
Ticks -> {Thread[{Range[0, 10], Range[0, 1, 0.1]}],
Thread[{Range[0, 10, 2], Range[0, 1, 0.2]}], Automatic},
PlotRange -> {{0, 10}, {0, 10}, {0, 1}}, Mesh -> False] Unfortunately, you cannot sample (RandomVariate) for the mixture distribution but you can use NProbability for questions, e.g.

NProbability[(x - .5)^2 + (y - .5)^2 <= 0.1, {x, y} \[Distributed]
mix]

yields: 0.666667

Again, I have voted for @Rojo but posted this as my interpretation of the desired outcome modifying (for reasons expressed above).

Straightforward method

One can do it very efficient without rejections and loops (100-500 times faster than posted methods)

n = 10000000;
r2 = 0.1;

AbsoluteTiming[
choise = RandomChoice[{π r2, 1} -> {0, 1}, n];
box = RandomReal[1, {n, 2}];
circle = Transpose@{0.5 + # Cos@#2, 0.5 + # Sin@#2} &[Sqrt@RandomReal[r2, n],
RandomReal[2 π, n]];
pts = box choise + circle (1 - choise);]
(* {1.972023, Null} *)

ListPlot[pts[[;; ;; 300]], AspectRatio -> Automatic] Histogram3D[pts, 25, "PDF"] Metropolis algorithm

Besides the above algorithm you can use the Matropolis algorithm. It is a very general algorithm and it cat take any probability density functions. See here for details and the definition of Metropolis

pdf = Function[{x, y}, UnitStep[x] UnitStep[y] UnitStep[1 - x] UnitStep[
1 - y] (1 + UnitStep[0.1 - (x - .5)^2 - (y - .5)^2])];
p = RandomVariate[Metropolis[pdf, {0.5, 0.5}, 0.1], 10000000]; // AbsoluteTiming
(* {11.082236, Null} *)

ListPlot[p[[;; ;; 300]], AspectRatio -> Automatic] Histogram3D[p, 20, "PDF"] There is only one limitation: pdf should be Listable.

You can also use RandomPoint with regions. First we construct a 3D region where the volume of the cylinder is twice the total volume of the cuboid minus the cylinder. Since the volume of the cylinder is $$\pi r^2 h_{cyl}$$ and the volume of the cuboid minus the cylinder is $$(1-\pi r^2) h_{cub}$$, we can find the values of $$h_{cyl}$$ and $$h_{cub}$$ when the volume of the cylinder is twice as large:$$\unicode{960}$$

Solve[π r^2 hc + (1 - π r^2) hs == 1 && π r^2 hc == 2 (1 - π r^2) hs, {hc, hs}]

{{hc -> 2/(3 π r^2), hs -> -(1/(3 (-1 + π r^2)))}}

Then, the region we want to construct is:

reg[x_, y_, r_] := RegionUnion[
Cuboid[{0, 0, 0}, {1, 1, 1/(3(1-π r^2))}],
Cylinder[{{x, y, 0}, {x, y, 2/(3 π r^2)}}, r]
]

For the OP example, then, we have:

r = reg[.5, .5, Sqrt[.1]]

BooleanRegion[#1 || #2 &, {Cuboid[{0, 0, 0}, {1, 1, 0.486021}], Cylinder[{{0.5, 0.5, 0}, {0.5, 0.5, 2.12207}}, 0.316228]}]

Let's check that the total volume is 1:

Volume[r]

1.

and that the volume of the Cylinder portion is 2/3:

Volume[Cylinder[{{0.5,0.5,0},{0.5,0.5,2.1220659078919377}},0.31622776601683794]]

0.666667

Now, we use RandomPoint to sample from the region, and keep only the first 2 coordinates:

rand = RandomPoint[r, 10^4]; //AbsoluteTiming

{0.055251, Null}

Visualization:

Graphics[{Blue, Point @ rand[[All, ;;2]]}] Let's check that the number of points in the cylinder (disk) are about 2/3 of the total (10000):

Count[RegionMember[Disk[{.5,.5}, Sqrt[.1]], rand[[All, ;;2]]], True]

6640