Transforming from a simplex in 4 dimensions to a 3D tetrahedron

Question

I have a collection of points that lie in the unit simplex in 4 dimensions. That is, each point is a vector of four non-negative real numbers that sum to 1. (They are probability distributions.)

This simplex in 4 dimensions is equivalent to a three dimensional tetrahedron. I would like to make a 3D plot of my points, each one plotted inside a regular tetrahedron, but I can't seem to work out the formula to transform my four-dimensional points into three-dimensional space in the correct way.

That is, I'm looking for a set of (presumably linear) equations that map the vectors $(1,0,0,0)$, $(0,1,0,0)$, $(0,0,1,0)$ and $(0,0,0,1)$ to the four corners of a regular tetrahedron, and the vector $(\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4})$ to its centre. (The position, scale and rotation of the tetrahedron don't matter.) Any insights into how to construct this mapping would be most appreciated.

Answer 1

This is a self-answer for the sake of others, generated from the helpful hints in the comments.

Ethan Bolker wrote:

Wikipedia gives several sets of coordinates for the vertices of a regular tetrahedron in $\mathbb{R}^3$. Just map your four vertices to one of those sets of four. en.wikipedia.org/wiki/Tetrahedron

and Rahul wrote:

Once you've chosen your vertices (say $\mathbf p_1,\dots,\mathbf p_4$), the mapping is just $(x_1,\dots,x_4)\mapsto x_1\mathbf p_1+\dots+x_4\mathbf p_4$.

This is straightforward once you realise it. For future visitors who just want a quick solution, here's my first attempt in Python, which is much simpler than what I was trying before, but puts the tetrahedron in a strange orientation, with an edge facing downwards:

def p_to_xyz(p):
    x = p[0] + p[1] - p[2] - p[3]
    y = p[0] - p[1] + p[2] - p[3]
    z = p[0] - p[1] - p[2] + p[3]
    return x, y, z

I prefer to draw it as a triangular pyramid, with a face pointing downward. After some further hints from this page, here's some code that does that:

def p_to_xyz(p):
    s3 = 1/math.sqrt(3.0)
    s6 = 1/math.sqrt(6.0)
    x =  -1*p[0] +  1*p[1] +    0*p[2] +    0*p[3]
    y = -s3*p[0] - s3*p[1] + 2*s3*p[2] +    0*p[3]
    z = -s3*p[0] - s3*p[1] -   s3*p[2] + 3*s6*p[3]
    return x, y, z

(If you're worried about efficiency you can remove the terms where something gets multiplied by zero, as I just put those in to make the structure clear.)