12

This is inspired by this answer which features code that I cannot get to work at all, and I tried to figure out why.

It boils down to using the function RandomPoint, which is new and a little bit tricky to use (see here and here), in combination with ImplicitRegion.

Here are two and three dimensional regions,

twoDRegion = 
  ImplicitRegion[{0 <= a <= 1, 0 <= b <= 1, a + b == 1}, {a, b}];
threeDRegion = 
  ImplicitRegion[{0 <= a <= 1, 0 <= b <= 1, 0 <= c <= 1, 
    a + b + c == 1}, {a, b, c}];

These regions are the set of all positive tuples that sum to 1,

RandomPoint works on the 2D version, but not the 3D version,

RandomPoint /@ {twoDRegion, threeDRegion}
(* {{0.146978, 0.853022}, 
 RandomPoint[
  ImplicitRegion[
   0 <= a <= 1 && 0 <= b <= 1 && 0 <= c <= 1 && a + b + c == 1, {a, b,
     c}]]} *)

Interestingly, RegionPlot also has trouble with the 3D version,

RegionPlot /@ {twoDRegion, threeDRegion}

enter image description here

but DiscretizeRegion can handle both

DiscretizeRegion /@ {twoDRegion, threeDRegion}

enter image description here

This led me to think that I just need to discretize before selecting the random point, and this in fact works for the 3D case

RandomPoint@DiscretizeRegion@# & /@ {twoDRegion, threeDRegion}
(* {{0.579605, 0.420395}, {0.154819, 0.0404491, 0.804732}} *)

But it won't work for higher dimensions since DiscretizeRegion is limited to 3 dimensions or lower.

Is this a bug?

Community
  • 1
Jason B.
  • 68,381
  • 3
  • 139
  • 286
  • I have run into the problem you describe here myself and have accepted the behavior as either an incomplete implementation or a design choice by the developers, not a bug. That is only a speculation on my part, of course. I hope you get an informed answer, perhaps from someone connected to WRI. – m_goldberg Feb 08 '16 at 12:03
  • Dear OP, did you find any solution for this issue? I also have similar requirement to use RandomPoint in higher dimension region. Thanks – Prashanth Nov 07 '16 at 02:57
  • Any word on a fix for this? I can't sample from BooleanRegions of with n>3. – Edmund May 18 '18 at 23:27

1 Answers1

3

In 12.1.1 this is still a problem. For a workaround I've noticed that making the a + b + c == 1 plane 'thicker' using 1 - 10^-6 <= a + b + c <= 1 allows you to generate points, and you can correct the small error by using #/Total[#] & so that a + b + c is exactly 1 without introducing any significant bias into the sampling:

#/Total[#] &@
 RandomPoint@
  ImplicitRegion[{0 <= a <= 1, 0 <= b <= 1, 0 <= c <= 1, 
    1 - 10^-6 < a + b + c <= 1}, {a, b, c}]

Interestingly, the thinner the final inequality gets the harder Mathematica has to work, and when close to the $MachineEpsilon it hangs indefinitely on my machine. There is definitely something Mathematica doesn't like about sampling in 3D with a strict equality constraint.

flinty
  • 25,147
  • 2
  • 20
  • 86