# How to generate random points in a region?

The Mathematica 10 documentation was updated for FindInstance adding support for regions.

In my use case, I'm trying to sample points in a set of disks:

region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]
FindInstance[{x, y} ∈ region, {x, y}, Reals, 2] // N


However the above code fails and generates the following error:

FindInstance::elemc: "Unable to resolve the domain or region membership condition {x,y}∈"


What's going wrong here?

-
Thanks @RunnyKine how can we escalate it to get a response? – user5601 Aug 22 '14 at 19:39
It would be nice if more people reported it to them and linked to these questions here like I did. – RunnyKine Aug 22 '14 at 19:40
FindInstance is not for generating random numbers ... If you already have the disks, it's trivial to find an instance of a point which lies within them. Take the centre of any disk. It doesn't even matter that you have several disks, you can just use a single one. If you have any expectation that these numbers will be uncorrelated and uniformly distributed (within the disks) then do not use FindInstance. – Szabolcs Aug 22 '14 at 20:09
Some of the answers here should point the way for WRI to include this functionality in Mma in the future. – Michael E2 Apr 23 at 19:21
@MichaelE2 your wish has been granted for Version 10.2! reference.wolfram.com/language/ref/RandomPoint.html – blochwave Jul 10 at 7:00

There are already good answers, but I'm going to improve the performance, generalize to any region in any dimensions and make the function more convenient. The main idea is to use DirichletDistribution (the uniform distribution on a simplex, e.g. triangle or tetrahedron). This idea was implemented by PlatoManiac and me in the related question obtaining random element of a set given by multiple inequalities (there is also Metropolis algorithm, but it is not suitable here).

The code is relatively short:

RegionDistribution /:
RandomDistributionVector[RegionDistribution[reg_MeshRegion], n_Integer, prec_?Positive] :=
Module[{d = RegionDimension@reg, cells, measures, s, m},
cells = DeveloperToPackedArray@MeshPrimitives[reg, d][[All, 1]];
s = RandomVariate[DirichletDistribution@ConstantArray[1, d + 1], n];
measures = PropertyValue[{reg, d}, MeshCellMeasure];
m = RandomVariate[#, n] &@EmpiricalDistribution[measures -> Range@Length@cells];
#[[All, 1]] (1 - Total[s, {2}]) + Total[#[[All, 2 ;;]] s, {2}] &@
cells[[m]]]


## Examples

### Random disks (2D in 2D)

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]


{0.004473, Null}

Precise test

pts = RandomVariate[RegionDistribution[region], 200000000]; // AbsoluteTiming


{85.835022, Null}

Histogram3D[pts, 50, "PDF", BoxRatios -> {Automatic, Automatic, 1.5}]


It is fast for $2\cdot10^8$ points and the distribution is really flat!

### Intervals (1D in 1D)

region = DiscretizeRegion[Interval[{0, 1}, {2, 4}]];
pts = RandomVariate[RegionDistribution[region], 100000]; // AbsoluteTiming
Histogram[Flatten@pts]


{0.062430, Null}

### Random circles (1D in 2D)

region = DiscretizeRegion@RegionUnion[Circle /@ RandomReal[10, {100, 2}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPlot[pts, AspectRatio -> Automatic]


{0.006216, Null}

### Balls (3D in 3D)

region = DiscretizeRegion@RegionUnion[Ball[{0, 0, 0}], Ball[{1.5, 0, 0}], Ball[{3, 0, 0}]];
pts = RandomVariate[RegionDistribution[region], 10000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]


{0.082202, Null}

### Surface cow disctribution (2D in 3D)

region = DiscretizeGraphics@ExampleData[{"Geometry3D", "Cow"}];
pts = RandomVariate[RegionDistribution[region], 2000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]


{0.026357, Null}

### Line in space (1D in 3D)

region = DiscretizeGraphics@ParametricPlot3D[{Sin[2 t], Cos[3 t], Cos[5 t]}, {t, 0, 2 π}];
pts = RandomVariate[RegionDistribution[region], 1000]; // AbsoluteTiming
ListPointPlot3D[pts, BoxRatios -> Automatic]


{0.005056, Null}

-
Very nice (+1). – kale Oct 9 '14 at 13:53
@ybeltukov Very nice answer. Could you just explain briefly why you use DirichletDistribution? – s.s.o Oct 9 '14 at 14:33
@s.s.o It is the uniform distribution on a simplex (e.g. triangle, tetrahedron). – ybeltukov Oct 9 '14 at 15:42
Ah, I should have known my randomBarycentric function was (the special case of) a well-known probability distribution. – Rahul Oct 9 '14 at 19:14
I love the cow distribution, +1. – rcollyer Apr 23 at 20:17

It would be nice if UniformDistribution worked on arbitrary regions, then we could simply do RandomVariate[UniformDistribution[region]]. Someone at Wolfram should get on that.

In the meantime, it seems we have to write our own sampling routines. @m_goldberg's answer is very nice (vote it up!) and uses rejection sampling, which works for arbitrary regions. However it will become slow if the measure of the region is small compared to that of its bounding box, for example if the disks are small and far apart. On the other hand, since you have a MeshRegion, it is actually possible to generate uniformly distributed samples without any rejection. We will pick a mesh cell randomly with probability proportional to its measure, then pick a point uniformly within that cell.

Sorry, I haven't bothered to wrap the code up into a neat module.

SeedRandom[0];
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]


d = RegionDimension[region];
coords = MeshCoordinates[region];
cells = MeshCells[region, d];
cellMeasures = PropertyValue[{region, d}, MeshCellMeasure];

randomCell[] := RandomChoice[cellMeasures -> cells]
randomBarycentric[] := Differences@Flatten@{0, Sort@RandomReal[1, d], 1}
randomRegionPoint[] := randomBarycentric[].coords[[First@randomCell[]]]

ListPlot[Table[randomRegionPoint[], {10000}], AspectRatio -> Automatic]


-
This is very nice upgrade version of goldberg's answer – Junho Lee Aug 23 '14 at 23:02
You can collect the already available cellMeasures like this: PropertyValue[{region, 2}, MeshCellMeasure] – RunnyKine Aug 23 '14 at 23:26
@RunnyKine: Thanks, added! – Rahul Aug 23 '14 at 23:59
@RahulNarain very nice...learned a lot...+1 – ubpdqn Aug 24 '14 at 8:34

Good news! Version 10.2 of Mathematica has this built-in with the function RandomPoint[]. From the documentation:

• RandomPoint can generate random points for any RegionQ region that is also ConstantRegionQ.
• RandomPoint will generate points uniformly in the region reg.

The first example given is a simple disk, but there are a whole host of neat applications given in the documentation!

pts = RandomPoint[Disk[], 5000];
Graphics[{PointSize[Tiny], Point[pts]}]


-

Here is a work-around that will let you extract as many randoms points from a region as you wish.

SeedRandom[0];
region =
DiscretizeRegion @ RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}]


randomFromRegion[region_, trials_] :=
Module[{bounds, randPts},
bounds = RegionBounds@region;
randPts =
Transpose @ {RandomReal[bounds[[1]], trials], RandomReal[bounds[[2]], trials]};
Select[randPts, RegionMember[region]]]


Let's see what we get from 5000 trials.

pts = randomFromRegion[region, 5000];
Length @ pts

2550


We extracted 2550 points that lie within the region. Here is how they are distributed.

Graphics @ Point @ pts


-
RegionBounds make rectangle region from circle. From the beginning, You ought to have made rectangle region. – Junho Lee Aug 22 '14 at 21:20
Sorry, but I don't understand your objection. – m_goldberg Aug 22 '14 at 21:30
@JunhoLee. If you were referring to my assertion concerning approximating area, I agree it was erroneous and have removed it. – m_goldberg Aug 22 '14 at 22:42
SeedRandom[0];

region = RegionUnion @ Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];

DiscretizeRegion @ region


points = RandomReal[{-1, 5}, {10000, 2}];

circles = List @@ (region /. Disk -> Circle);

Graphics[{AbsolutePointSize[2],
Transpose[{RegionMember[region, points] /.
{False -> White, True -> Gray}, Point /@ points}], circles}, Frame -> True]


-
thanks, but (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal[1]], {10000}]; RegionMember[region, RandomReal[{0, 100}, {100000, 2}]]) crashes the kernel – user5601 Aug 22 '14 at 20:46
@user5601 Your question is "How to generate random points in a region". It's not about FindInstance (see several comments), and it's not about the endless possibilities of "crashing the kernel". I answered your question using your own definition of region. – eldo Aug 22 '14 at 21:05
Right! I should have specified that I need many points! – user5601 Aug 25 '14 at 15:08
region = RegionUnion @ Table[Disk[RandomReal[4, {2}], DiscretizeRegion is very slow, so if the number of disks exceeds 10000 then it freezes... – user5601 Aug 25 '14 at 15:09

In version 10.1 the undocumented function RandomRandomPointVector is useful:

region = DiscretizeRegion@RegionUnion@
Table[Disk[RandomReal[4, {2}], RandomReal[1]], {10}];

Graphics@Point@
RandomRandomPointVector[region, 1000, Automatic, Automatic]


The two Automatic arguments appear to be working precision and a method option - other allowed values are "Mesh" and "Rejection"

-
I think this may have developed into 10.2's RandomPoint[] – blochwave Jul 10 at 8:26

FindInstance have to apply to a expression like an equation or an inequality. DiscretizeRegion does not make an inequality. You should apply like the following example.

region = RegionUnion@Table[Disk[RandomReal[4, 2], RandomReal[1]],
{10}
];
DiscretizeRegion[region]


RegionQ[region]


True

FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N


{{x -> 2.4597, y -> 4.09129}, {x -> 1.31552, y -> 2.00384}}

-
But DiscretizeRegion is a region. Because RegionQ[] returns True – RunnyKine Aug 22 '14 at 20:48
Seems like it is too slow: (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal[1]], {10000}]; FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N) – user5601 Aug 22 '14 at 20:48
@user5601 Please check technique applying method of FindInstance. I might guess you find the solution like this posting example. – Junho Lee Aug 22 '14 at 20:50
@RunnyKine you right. DiscretizeRegion is not a quantified system of equations and inequalities. – Junho Lee Aug 22 '14 at 20:55
@user5601 RegionMember[region, {x, y}] make inequalities. too slow is natural.. – Junho Lee Aug 22 '14 at 20:56

you can look at this question here

This method also may work fine:

p1 = {x, y} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
p2 = {y, x} /. FindInstance[x^2 + y^2 < 1, {x, y}, Reals, 1000];
point = RandomChoice[Join[p1, p2], 1000];
Graphics[{Circle[], {Red, Point[point]}}]

-
FindInstance is not a random number generator. There's no guarantee that it will give you uncorrelated random numbers with a uniform distribution. What the OP is asking for a misuse of FindInstance that is going to lead to wrong results. – Szabolcs Aug 22 '14 at 20:12
I know FindInstance is not a random number generator but I think they will be uniformly distributed even if the pointes are correlated. then I used RandomChoice to get random pointes. I don't know if this way is correct or not. – Algohi Aug 23 '14 at 17:02