# 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], {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? Aug 22, 2014 at 19:39
• It would be nice if more people reported it to them and linked to these questions here like I did. Aug 22, 2014 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. Aug 22, 2014 at 20:09
• Some of the answers here should point the way for WRI to include this functionality in Mma in the future. Apr 23, 2015 at 19:21
• @MichaelE2 your wish has been granted for Version 10.2! reference.wolfram.com/language/ref/RandomPoint.html Jul 10, 2015 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;
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal], {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, 2014 at 13:53
• @ybeltukov Very nice answer. Could you just explain briefly why you use DirichletDistribution? Oct 9, 2014 at 14:33
• @s.s.o It is the uniform distribution on a simplex (e.g. triangle, tetrahedron). Oct 9, 2014 at 15:42
• Ah, I should have known my randomBarycentric function was (the special case of) a well-known probability distribution.
– user484
Oct 9, 2014 at 19:14
• I love the cow distribution, +1. Apr 23, 2015 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;
region = DiscretizeRegion@RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal], {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 Aug 23, 2014 at 23:02
• You can collect the already available cellMeasures like this: PropertyValue[{region, 2}, MeshCellMeasure] Aug 23, 2014 at 23:26
– user484
Aug 23, 2014 at 23:59
• @RahulNarain very nice...learned a lot...+1 Aug 24, 2014 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;
region =
DiscretizeRegion @ RegionUnion@Table[Disk[RandomReal[4, {2}], RandomReal], {10}] randomFromRegion[region_, trials_] :=
Module[{bounds, randPts},
bounds = RegionBounds@region;
randPts =
Transpose @ {RandomReal[bounds[], trials], RandomReal[bounds[], 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. Aug 22, 2014 at 21:20
• Sorry, but I don't understand your objection. Aug 22, 2014 at 21:30
• @JunhoLee. If you were referring to my assertion concerning approximating area, I agree it was erroneous and have removed it. Aug 22, 2014 at 22:42

In version 10.1 the undocumented function RandomRandomPointVector is useful:

region = DiscretizeRegion@RegionUnion@
Table[Disk[RandomReal[4, {2}], RandomReal], {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[] Jul 10, 2015 at 8:26
SeedRandom;

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

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

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

Graphics[{AbsolutePointSize,
Transpose[{RegionMember[region, points] /.
{False -> White, True -> Gray}, Point /@ points}], circles}, Frame -> True] • thanks, but (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal], {10000}]; RegionMember[region, RandomReal[{0, 100}, {100000, 2}]]) crashes the kernel Aug 22, 2014 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, 2014 at 21:05
• Right! I should have specified that I need many points! Aug 25, 2014 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... Aug 25, 2014 at 15:09
• @eldo how to find the boudary nodes for these collocation points? Jul 13, 2021 at 16:05

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],
{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 Aug 22, 2014 at 20:48
• Seems like it is too slow: (region = RegionUnion@ Table[Disk[RandomReal[{0, 100}, {2}], RandomReal], {10000}]; FindInstance[RegionMember[region, {x, y}], {x, y}, Reals, 2] // N) Aug 22, 2014 at 20:48
• @user5601 Please check technique applying method of FindInstance. I might guess you find the solution like this posting example. Aug 22, 2014 at 20:50
• @RunnyKine you right. DiscretizeRegion is not a quantified system of equations and inequalities. Aug 22, 2014 at 20:55
• @user5601 RegionMember[region, {x, y}] make inequalities. too slow is natural.. Aug 22, 2014 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. Aug 22, 2014 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. Aug 23, 2014 at 17:02