# 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? – belisarius has settled 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[[2]] > pt[[1]]*v/scale, pt = {scale - pt[[1]], pt[[2]] - v }]; {{c, s}, {-s, c}}.pt , {n}]]]; – Daniel Lichtblau Sep 23 '13 at 2:55

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[[1]], 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],


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[[1]], 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 :) – belisarius has settled 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 – belisarius has settled 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[[1]], 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...

-