2

I'd like to create a list of triplets $(x,y,z)$ which satisfy the following properties: $$0<x<1/2 \\ 0<y<1/2-x \\ -1<z<-2(x+y)$$ Basically I would like to uniformly generate random points inside a pyramid which is what these inequalities represent. I'm not sure how to place these restrictions on the RandomReal command or how to tell Mathematica to exclude points that do not satisfy these conditions. Any help?

Thank you.

Sheheryar Zaidi
  • 209
  • 1
  • 4

3 Answers3

5

You can use RandomPoint with ImplicitRegion:

RandomPoint[
    ImplicitRegion[0<x<1/2 && 0<y<1/2 && -1<z<-2(x+y), {x,y,z}],
    10
]

{{0.0644104, 0.289938, -0.812196}, {0.162676, 0.0447583, -0.805135}, {0.29663, 0.0801703, -0.955995}, {0.0845065, 0.0982267, -0.542481}, {0.00591773, 0.0586124, -0.222642}, {0.395283, 0.0611798, -0.965819}, {0.281329, 0.0247549, -0.985743}, {0.0997621, 0.156397, -0.969405}, {0.136775, 0.176456, -0.934246}, {0.165527, 0.234188, -0.836819}}

Or, if you know the corners of the tetrahedron, you could use:

RandomPoint[
    Tetrahedron[{{0,0,-1},{0,1/2,-1},{1/2,0,-1},{0,0,0}}],
    10
]

{{0.103469, 0.00460676, -0.804016}, {0.155514, 0.266517, -0.934628}, {0.0258173, 0.0425049, -0.344449}, {0.161183, 0.0999294, -0.639421}, {0.124769, 0.194545, -0.914526}, {0.223348, 0.08202, -0.620869}, {0.0379937, 0.0701453, -0.401479}, {0.0608969, 0.0313861, -0.458453}, {0.00207286, 0.136834, -0.477351}, {0.246812, 0.0634065, -0.89431}}

Carl Woll
  • 130,679
  • 6
  • 243
  • 355
5

The following seems to work:

Cases[RandomReal[{-1, 1}, {1500000, 3}], {x_, y_, z_} /; 
  0 < x < 0.5 && 0 < y < 0.5 - x && -1 < z < -2 (x + y)]

OR

Select[RandomReal[{-1, 1}, {1500000, 3}], 
 0 < #[[1]] < 0.5 && 0 < #[[2]] < 0.5 - #[[1]] && -1 < #[[3]] < -2 (#[[1]] + #[[2]]) &]

Here is what it looks like:

Mathematica graphics

RunnyKine
  • 33,088
  • 3
  • 109
  • 176
  • I believe that generation and filtering is a fine approach when the selection percentage is relatively high as is the case here. +1 – Mr.Wizard Aug 15 '14 at 01:13
  • @Mr.Wizard Thanks. I thought about a different approach but this seems to do the job in this case. – RunnyKine Aug 15 '14 at 01:16
  • See links above for advanced methods for this class of problem. – Mr.Wizard Aug 15 '14 at 01:26
  • @Mr.Wizard. Thanks, I actually voted in 2 out of the 3 links you provided. Just forgot about those. – RunnyKine Aug 15 '14 at 01:42
  • 3
    @Mr.Wizard I think it should be pointed out that about 0.5% of the points will be selected. RunnyKine, it would more efficient to use separate RandomReal calls for each coordinate, transpose, and select. The expected yield would be 1/6 or about 17%. – Michael E2 Aug 15 '14 at 03:30
  • @MichaelE2 I missed that and assumed the larger figure. Good catch. – Mr.Wizard Aug 15 '14 at 03:46
5

Ray Koopman's answer to Uniformly distributed n-dimensional probability vectors over a simplex applied to a simplex with given vertices leads to the following way to get (exactly) n uniformly distributed random points in the simplex (and do it quickly):

pts = #.vertices/Total[#, {2}] &@ Log @ RandomReal[1, {n, 4}];

Here are 2000 points in the OP's simplex:

vertices = Append[
  {x, y, z} /. 
     Last@Maximize[{#, 
         0 < x < 1/2 && 0 < y < 1/2 (1 - 2 x) && -1 < z < -2 x - 2 y} /. 
        Less -> LessEqual, {x, y, z}] & /@ {x, y, -z},
  {0, 0, 0}
  ]
(* {{1/2, 0, -1}, {0, 1/2, -1}, {0, 0, -1}, {0, 0, 0}} *)

pts = #.vertices/Total[#, {2}] &@ Log @ RandomReal[1, {2000, 4}];

Graphics3D[Point[pts], Axes -> True, 
 PlotRange -> {{0, 1/2}, {0, 1/2}, {-1, 0}}]

Mathematica graphics

Community
  • 1
Michael E2
  • 235,386
  • 17
  • 334
  • 747