2
$\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$
2
  • $\begingroup$ Related: (13038), (32501), (33652) $\endgroup$
    – Mr.Wizard
    Aug 15, 2014 at 1:23
  • $\begingroup$ @Mr.Wizard It seems pretty close to a duplicate of 33652 -- see my answer. $\endgroup$
    – Michael E2
    Aug 15, 2014 at 3:12

3 Answers 3

5
$\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$
6
  • $\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, 2014 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, 2014 at 1:16
  • $\begingroup$ See links above for advanced methods for this class of problem. $\endgroup$
    – Mr.Wizard
    Aug 15, 2014 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, 2014 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, 2014 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$
1
  • $\begingroup$ I think this is a good extension of the linked question rather than a duplicate. +1 $\endgroup$
    – Mr.Wizard
    Aug 15, 2014 at 3:47
5
$\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.