How to generate random points in a region?

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 /: 
 Random`DistributionVector[RegionDistribution[reg_MeshRegion], n_Integer, prec_?Positive] :=
  Module[{d = RegionDimension@reg, cells, measures, s, m},
   cells = Developer`ToPackedArray@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}

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]

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]}]