# obtaining random element of a set given by multiple inequalities

I try to obtain random elements from a set given by multiple inequalities. This is of general interest to me, but let me present one example here:

I have multiple inequalities specifying half spaces in 2D given by:

Eta[a_] := {Cos[a], Sin[a]}
NI[a_] := {Cos[a], Sin[a]}
Table[Dot[{Phi1, Phi2}, Eta[b]]] <= Norm[NI[Pi] - Eta[b]]^2 + 2,
{b, 0, 2 Pi, 2 Pi/10}]


The set I am interested in is the intersection of all these inequalities (half spaces). Here is a RegionPlot of the set:

RegionPlot[And @@ Table[Dot[{Phi1, Phi2}, Eta[b]]] <= Norm[NI[Pi] - Eta[b]]^2 + 2,
{b, 0, 2 Pi, 2 Pi/10}], {Phi1, -7, 7}, {Phi2, -7, 7}] Now I want to pick random vectors (elements) from this set.

My current approach is to define an indicator function for the set like this:

set[x1_,x2_] = And @@ Table[Dot[{Phi1, Phi2}, Eta[b]]] <=
Norm[NI[Pi] - Eta[b]]^2 + 2,
{b, 0, 2 Pi, 2 Pi/10}] /. {Phi1 -> x1, Phi2 -> x2}


Then use a random number generator for x1 and x2 with appropriate bounds, that I get from inspecting the RegionPlot of the set to obtain a function that gives me n "random" elements from my set:

randomVector[n_] := Module[{list = {}, vector, i = 0},
While[i < n,
vector = {RandomReal[{-3, 6}], RandomReal[{-5, 5}]};
If[set[Sequence @@ vector], i++; AppendTo[list, vector];];
];
Return[list];]


Do you know of a more elegant way to do this?

• Are you working with reg polygons? – Dr. belisarius Sep 17 '13 at 13:36
• @belisarius: In the above case, yes, but I am more interested in a general solution that handles sets defined by (preferably) nonlinear inequalities. – Wizard Sep 17 '13 at 14:19
• One approach would be to use your method above to generate a collection of random points. Then use EmpiricalDistribution or SmoothKernelDistribution on these points to define a distribution. Using this, it should be quicker and easier to generate lots more data that lies in the shape you want. – bill s Sep 17 '13 at 15:06
• Fast code for the special case of generating n points in a regular something-gon, with given scale factor: randomPoints = Compile[{{n, _Integer}, {scale, _Real}, {sectors, _Integer}}, Module[ {v, r, rot, pt, c, s}, r = scale*Sec[Pi/sectors]; v = r*Sin[Pi/N[sectors]]; Table[rot = RandomInteger[{0, sectors - 1}]; pt = {RandomReal[{0, scale}], RandomReal[{0, v}]}; c = Cos[2.*Pi*rot/sectors]; s = Sin[2.*Pi*rot/sectors]; If[pt[] > pt[]*v/scale, pt = {scale - pt[], pt[] - v }]; {{c, s}, {-s, c}}.pt , {n}]]]; – Daniel Lichtblau Sep 23 '13 at 2:55

## 5 Answers

Lets call your plot res.

res = RegionPlot[And @@ Table[
Dot[{Phi1, Phi2}, Eta[b]] <= Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0,
2 Pi, 2 Pi/10}], {Phi1, -7, 7}, {Phi2, -7, 7}];


Lets extract the mesh Mathematica is generating by default. Use more PlotPoints to get more triangular mesh of your 2D region.

pts = res[[1, 1]]; (* Vertices *)
{triangles, qd} = Cases[res[], Polygon[{a___}] -> {a}, Infinity]; (* Triangle and Quads*)
quadTotri = Flatten[{Drop[#, {2}], Take[#, 3]} & /@ qd, 1];(* Quad to Triangle*)
Graphics[{FaceForm[], EdgeForm[Red],GraphicsComplex[pts, Polygon@triangles],
FaceForm[],EdgeForm[{Blue, Dashed}], GraphicsComplex[pts, Polygon@quadTotri]}] Now create a distribution based on the area of those triangles.

allTrig = triangles~Join~quadTotri;
vetices = (Extract[pts, Transpose@{#}] & /@ allTrig);
area = .5 Det@((Append[#, 1]) & /@ #) &; (* Calculate area of triangle *)
dat = area /@ (Extract[pts, Transpose@{#}] & /@ (allTrig));
d = EmpiricalDistribution[dat -> Range[Length@allTrig]];


Use random barycentric coordinates to choose the random points from those triangles.

Randpts = (Dot[#/Total[#] &@RandomReal[1, 3], #] & /@
Transpose@vetices[[#]]) & /@ RandomVariate[d, 1000];
Show[res, Graphics[{Red, PointSize[Small], Point /@ Randpts}]]


Tuning:

Using suggestions from @ybeltukov the above can be tied in a function. It takes as argument the RegionPlot graphics output of your region and the number $n$ of random points you want to pick/sample from the region.

Clear[RegionRandom];
RegionRandom[plot_Graphics, n_] := Block[{pts, triangles, qd, quadTotri,
allTrig, vetices, areas,
empdist, u0, u1, u2, CustomDistribution, rp},
pts = plot[[1, 1]];
{triangles, qd} =Cases[plot[], Polygon[{a___}] -> {a}, Infinity];
allTrig =
triangles~Join~Flatten[{Drop[#, {2}], Take[#, 3]} & /@ qd, 1];
vetices = (Extract[pts, Transpose@{#}] & /@ allTrig);
areas =
0.5 Abs[#1[[All, 2]] #2[[All, 1]] - #1[[All, 1]] #2[[All,
2]]] &[#[[All, 2]] - #[[All, 1]], #[[All, 3]] - #[[All,
1]]] &@vetices;
empdist = EmpiricalDistribution[areas -> Range@Length@vetices];
u0 = vetices[[All, 1]];
u1 = vetices[[All, 2]] - u0;
u2 = vetices[[All, 3]] - u0;
CustomDistribution /:
RandomDistributionVector[CustomDistribution[], p_Integer,
prec_?Positive] :=
Module[{s =
RandomVariate[DirichletDistribution[{1, 1, 1}], p,
WorkingPrecision -> prec],
m = RandomVariate[empdist, p, WorkingPrecision -> prec]},
u0[[m]] + s[[All, 1]] u1[[m]] + s[[All, 2]] u2[[m]]
];
rp = RandomVariate[CustomDistribution[], n]
];


Lets test it with above RegionPlot named res for $40000$ random points.

pt = RegionRandom[res, 40000];
Show[res, Graphics[{Red, PointSize[Tiny], Point /@ pt}]] Also some fun with nontrivial 2D regions. Result is not too bad! You can use this answer to create nicer 2D mesh of your region than the default one. I may update that later!

Comparison:

Histogram: First lets see the histogram for $20000000$ sample!

Histogram3D[RegionRandom[res, 20000000], 25] Entropy of generated data: We know entropy is an information theoretic metric to measure the intrinsic randomness of sample. We use Entropy function to test which algorithm provides more randomness. Mahalanobis distance: We also compare the Mahalanobis distance pdf of both the sample generators as it gives a visual presentation of their mutual disagreement.

n=50000; (* reduce n to as the following is very memory intensive ~ I had 64 Gb RAM *)
{TrigData, MetData} = {RegionRandom[res, n],
RandomVariate[Metropolis[pdf, {2, 0}], n]};
MahalanobisDistance[data_] :=
Block[{m = Mean[data], cov = Inverse[Covariance[data]], temp},
temp = (#1 - m &) /@ data; Diagonal[temp.cov.Transpose[temp]]]
dist = SmoothKernelDistribution[MahalanobisDistance[#]] & /@ {TrigData, MetData}; Timing: The Metropolis algorithm performs faster after compilation by @ybeltukov

pt1 = RandomVariate[Metropolis[pdf, {2, 0}],20000000]; // AbsoluteTiming


{9.653125, Null}

pt = RegionRandom[res, 20000000]; // AbsoluteTiming


{14.207429, Null}

• +1 but you have a leak at the north west :) – Dr. belisarius Sep 17 '13 at 17:16
• @belisarius I guess it is some Graphics issue. If you make the pic bigger they come inside ;).... – PlatoManiac Sep 17 '13 at 17:32
• barrycentric – Dr. belisarius Sep 17 '13 at 17:38
• +1 for great answer! As an alternative you can use Dirichlet distribution to choose the random points in triangles. For example, RandomVariate[ TransformedDistribution[{1, 0} s1 + {0, 1} s2, {s1, s2} \[Distributed] DirichletDistribution[{1, 1, 1}]], 1500] – ybeltukov Sep 17 '13 at 19:09
• @ybeltukov That's why there was a "Result is not too bad!" In the concluding text. Hope you noticed! – PlatoManiac Sep 17 '13 at 20:13

# Metropolis algorithm

Update: ~15x speedup with Compile!

I propose an original solution, which consists in using the Metropolis algorithm. It is a very general approach, which is applicable for any probability density function in any dimensions.

Metropolis /:
RandomDistributionVector[
Metropolis[pdf_, u0_, s_: 1, n0_: 100, chains_: 200], n_Integer,
prec_?Positive] :=
Module[{u, du, p, p1, accept, cpdf},
cpdf = Compile @@ {{#, _Real} & /@ #, pdf @@ #, RuntimeAttributes -> {Listable},
RuntimeOptions -> "Speed"} &[Unique["x", Temporary] & /@ u0];
u = ConstantArray[u0, chains];
p = cpdf @@ Transpose[u];
(Join @@ Table[
du = RandomVariate[NormalDistribution[0, s], {chains, Length[u0]}];
p1 = cpdf @@ Transpose[u + du];
accept = UnitStep[p1/p - RandomReal[{0, 1}, chains]];
p += (p1 - p) accept;
u += du accept
, {Ceiling[(n0 + n)/chains]}])[[n0 + 1 ;; n0 + n]]
];


Here pdf is a custom probability density function, u0 is an initial point, s is a step size, n0 is a number of step to forgot initial state, chains is the number of simultaneous Markov chains.

Examples

1. Multiple inequalities

Eta[a_] := {Cos[a], Sin[a]};
NI[a_] := {Cos[a], Sin[a]};
pdf = Evaluate@Boole[And @@
Table[N@Dot[{#1,#2}, Eta[b]] <= N@Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0, 2Pi, 2Pi/10}]] &;

p = RandomVariate[Metropolis[pdf, {0, 0}], 30000];
ListPlot[p, AspectRatio -> Automatic] Let's check the distribution

Histogram3D[RandomVariate[Metropolis[pdf, {2, 0}], 20000000], 25] It is really uniform!

2. Custom probability density function

pdf = Cos[Pi Sqrt[#1^2 + #2^2]]^2 Exp[-#1^2 - #2^2] &;

p = RandomVariate[Metropolis[pdf, {0, 0}], 30000];
ListPlot[p, AspectRatio -> Automatic] Again, the distribution checking

nrm = NIntegrate[pdf[x, y], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}];
Show[Histogram3D[
RandomVariate[Metropolis[pdf, {2, 0}],
2000000], {{-2, 2, 0.1}, {0, 2, 0.1}}, "PDF",
BoxRatios -> {2, 1, 1}],
ParametricPlot3D[{x, 0, 2 pdf[x, 0]/nrm}, {x, -5, 5},
PlotStyle -> {Red, Thick}]] For users of Mathematica 7 and earlier:

Metropolis /:
RandomDistributionVector[
Metropolis[pdf_, u0_, s_: 1, n0_: 100, chains_: 200], n_Integer,
prec_?Positive] :=
Module[{u, du, p, p1, accept}, u = ConstantArray[u0, chains];
p = pdf @@@ u;
(Join @@
Table[du = RandomReal[NormalDistribution[0, s], {chains, Length[u0]}];
p1 = pdf @@@ (u + du);
accept = UnitStep[p1/p - RandomReal[{0, 1}, chains]];
p += (p1 - p) accept;
u += du accept, {Ceiling[(n0 + n)/chains]}])[[n0 + 1 ;; n0 + n]]];

p = RandomReal[Metropolis[pdf, {0, 0}], 10000];

• Ah! This is really the time-tested, practical solution that's both robust and stands on formalized basis. Lots of scientific computing over half a century has relied on it. I'd wish I could give multiple upvotes! – kirma Sep 18 '13 at 5:13
• This is nice indeed! This is pretty general..Thx for finding it. +1 – PlatoManiac Sep 18 '13 at 7:35
• By the way you may also supply a measure using some statistic (PearsonChiSquareTest or KolmogorovSmirnovTest) to show how uniformity distributed your result really is? – PlatoManiac Sep 18 '13 at 9:08
• @PlatoManiac I upgrade the algorithm and provide direct distribution comparison. Do you know how to use these tests with custom distribution? – ybeltukov Sep 18 '13 at 16:12
• @ybeltukov Very cool update! I will check about the tests. – PlatoManiac Sep 18 '13 at 19:41

I enjoy the solution of PlatoManiac and I want to improve it

Eta[a_] := {Cos[a], Sin[a]}
NI[a_] := {Cos[a], Sin[a]}

res = RegionPlot[
And @@ Table[
Dot[{Phi1, Phi2}, Eta[b]] <= Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0,
2 Pi, 2 Pi/10}], {Phi1, -7, 7}, {Phi2, -7, 7}];

points = res[[1, 1]];
{id3, id4} = Cases[res[], Polygon[{a___}] -> {a}, Infinity];
id3 = Join[id3, id4[[All, {1, 2, 3}]], id4[[All, {3, 4, 1}]]];
triangles = Partition[points[[Flatten[id3]]], 3];

Graphics[{Lighter[Red], EdgeForm[Black], Polygon@triangles}] areas = 0.5 Abs[#1[[All, 2]] #2[[All, 1]] - #1[[All, 1]] #2[[All, 2]]] &
[#[[All, 2]] - #[[All, 1]], #[[All, 3]] - #[[All, 1]]] &@triangles;
empdist = EmpiricalDistribution[areas -> Range@Length@triangles];
u0 = triangles[[All, 1]];
u1 = triangles[[All, 2]] - u0;
u2 = triangles[[All, 3]] - u0;


Now we define own generator

CustomDistribution /:
RandomDistributionVector[CustomDistribution[], n_Integer, prec_?Positive] :=
Module[{s = RandomVariate[DirichletDistribution[{1, 1, 1}], n, WorkingPrecision -> prec],
m = RandomVariate[empdist, n, WorkingPrecision -> prec]},
u0[[m]] + s[[All, 1]] u1[[m]] + s[[All, 2]] u2[[m]]]


It can be used as a standard generator

p = RandomVariate[CustomDistribution[], 30000];
ListPlot[p, AspectRatio -> Automatic] This implementation is about 3x faster than the original solution of PlatoManiac and produce uniform distribution in the whole range of interest.

My submission is some sort of exercise in bull-headed brute force with some fragile hacks sprinkled in. Apart from the commented parts, I hope it's reasonably clear what it does. Run-time and memory use of this implementation rises rapidly with shape complexity; thus here I use only a pentagon. It is likely that this method doesn't work for more complicated regions for which Integrate[Boole[...], ...] construct fails to produce meaningful symbolic results.

Eta[a_] := {Cos[a], Sin[a]}
NI[a_] := {Cos[a], Sin[a]}

region[u_, v_] := And @@ Table[
Dot[{u, v}, Eta[b]] <= Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0, 2 Pi, 2 Pi/5}]

generateRegionMapping[region_] :=
Module[{vp, ux, ua, up},

vp = Solve[
Integrate[
region[u, v], {v, -\[Infinity], x}, {u, -\[Infinity], \[Infinity]},
Assumptions -> x \[Element] Reals] /
Integrate[
region[u, v],
{v, -\[Infinity], \[Infinity]}, {u, -\[Infinity], \[Infinity]}] == p //
FullSimplify, x, Reals, Method -> Reduce] // FullSimplify;

ux = Integrate[region[u, v], {u, -\[Infinity], x},
Assumptions -> {v, x} \[Element] Reals] // FullSimplify;

ua = Integrate[region[u, v], {u, -\[Infinity], \[Infinity]},
Assumptions -> x \[Element] Reals] // FullSimplify;

(* Trickery starts here: put divisions inside and simplify. *)
ux = ua /. {val_, ineq_Inequality} -> {ux/val, ineq} // FullSimplify;

(* Remove term that would return 1 - it's unnecessary for us. *)
ux = Piecewise@Cases[ux // First, Except[{1, _}]];

(* Remove comparisons that involve x.
Those only complicate solving in our case. *)
ux = ux //. ((v_Greater | v_Less | v_GreaterEqual | v_LessEqual | v_Inequality) /;
MemberQ[v, x, Infinity]) :> True // FullSimplify;

up = Solve[ux == p, x, Reals, Method -> Reduce] // FullSimplify;

Function[{p1, p2},
With[{vv = Select[x /. (vp /. p -> p1), NumericQ] // First},
{Select[x /. (up /. {v -> vv, p -> p2}), NumericQ] // First, vv}]]]

mapper = generateRegionMapping[Boole[region@##] &]; // AbsoluteTiming


{126.051116, Null}

Show[
RegionPlot[region[Phi1, Phi2], {Phi1, -7, 7}, {Phi2, -7, 7}],
ListPlot[mapper @@@ RandomReal[{0, 1}, {3000, 2}], AspectRatio -> 1]] Since v10.2 RandomPoint has provided a way to pick uniform samples from geometric regions (which you can trivially derive from your specification using ImplicitRegion):

Eta[a_] := {Cos[a], Sin[a]};
NI[a_] := {Cos[a], Sin[a]};
reg = ImplicitRegion[
And @@ Table[
Dot[{x, y}, Eta[b]] <= Norm[NI[Pi] - Eta[b]]^2 + 2, {b, 0, 2 Pi,
2 Pi/10}], {x, y}];
Show[{RegionPlot[reg],
ListPlot[Quiet@RandomPoint[reg, 10000], PlotStyle -> Red]}] RandomPoint is pretty nifty - when it works efficiently. It is not particularly hard to construct regions where that doesn't apply. Thankfully it does work on most of simple cases one might explore...