1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Related: (13038), (32501), (33652) $\endgroup$ – Mr.Wizard Aug 15 '14 at 1:23
  • $\begingroup$ @Mr.Wizard It seems pretty close to a duplicate of 33652 -- see my answer. $\endgroup$ – Michael E2 Aug 15 '14 at 3:12
4
$\begingroup$

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

$\endgroup$
  • $\begingroup$ I believe that generation and filtering is a fine approach when the selection percentage is relatively high as is the case here. +1 $\endgroup$ – Mr.Wizard Aug 15 '14 at 1:13
  • $\begingroup$ @Mr.Wizard Thanks. I thought about a different approach but this seems to do the job in this case. $\endgroup$ – RunnyKine Aug 15 '14 at 1:16
  • $\begingroup$ See links above for advanced methods for this class of problem. $\endgroup$ – Mr.Wizard Aug 15 '14 at 1:26
  • $\begingroup$ @Mr.Wizard. Thanks, I actually voted in 2 out of the 3 links you provided. Just forgot about those. $\endgroup$ – RunnyKine Aug 15 '14 at 1:42
  • 3
    $\begingroup$ @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%. $\endgroup$ – Michael E2 Aug 15 '14 at 3:30
5
$\begingroup$

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

$\endgroup$
  • $\begingroup$ I think this is a good extension of the linked question rather than a duplicate. +1 $\endgroup$ – Mr.Wizard Aug 15 '14 at 3:47
4
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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