obtaining random element of a set given by multiple inequalities

Question

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?

Answer 1

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], 
         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[[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 /:
Random`DistributionVector[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}

Answer 2

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 /: 
  Random`DistributionVector[
   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 /: 
  Random`DistributionVector[
   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];

Answer 3

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 /: 
 Random`DistributionVector[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.

Answer 4

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

Answer 5

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...