Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

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?

share|improve this question
    
Are you working with reg polygons? –  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[[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

4 Answers 4

up vote 16 down vote accepted

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

enter image description here

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

enter image description here

Also some fun with nontrivial 2D regions.

enter image description here

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]

enter image description here

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. enter image description here

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

enter image description here

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}

share|improve this answer
    
+1 but you have a leak at the north west :) –  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 –  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
2  
I have to unupvote because your distribution is not uniform :( It is clearly seen for 10000 point. –  ybeltukov Sep 17 '13 at 20:08

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]
    

    enter image description here

    Let's check the distribution

    Histogram3D[RandomVariate[Metropolis[pdf, {2, 0}], 20000000], 25]
    

    enter image description here

    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]
    

    enter image description here

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

    enter image description here


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];
share|improve this answer
    
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
1  
@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}]

enter image description here

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]

enter image description here

This implementation is about 3x faster than the original solution of PlatoManiac and produce uniform distribution in the whole range of interest.

share|improve this answer

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

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.