PDF of volume of tetrahedron with random coordinates

Question

Question

What is the probability distribution function (PDF) of the absolute volume of a tetrahedron with random coordinates?

The 4 random tetrahedron vertices in $\mathbb{R}^3$ are $$ \mathbf{\mathrm{X}_1} =(x_1^1,x_1^2,x_1^3),\;\; \mathbf{\mathrm{X}_2}=(x_2^1,x_2^2,x_2^3),\;\; \mathbf{\mathrm{X}_3}=(x_3^1,x_3^2,x_3^3),\;\; \mathbf{\mathrm{X}_4} =(x_4^1,x_4^2,x_4^3)$$ where $x_i^j$ are independent standard normal distributed variables $$x_i^j\sim\mathcal{N}(0,1)$$ The non-oriented volume of a random tetrahedron instance is $$V=\frac{1}{6}\left| ( \mathbf{\mathrm{X}_1}- \mathbf{\mathrm{X}_4})\cdot \left(( \mathbf{\mathrm{X}_2}- \mathbf{\mathrm{X}_4} ) \times ( \mathbf{\mathrm{X}_3}- \mathbf{\mathrm{X}_4} )\right) \right| \tag{1} $$ $$=\frac{1}{6}\left|x_1^2 x_2^3 x_3^1 - x_1^1 x_2^3 x_3^2 + x_1^3 x_2^1 x_3^2- x_1^2 x_2^1 x_3^3 + x_1^1 x_2^2 x_3^3- x_1^3 x_2^2 x_3^1 + x_1^3 x_2^2 x_4^1- x_1^1 x_2^2 x_4^3 + x_1^1 x_2^3 x_4^2- x_1^2 x_2^3 x_4^1 + x_1^2 x_2^1 x_4^3- x_1^3 x_2^1 x_4^2 + x_2^3 x_3^2 x_4^1 -x_2^1 x_3^2 x_4^3 + x_2^1 x_3^3 x_4^2- x_2^2 x_3^3 x_4^1 + x_2^2 x_3^1 x_4^3- x_2^3 x_3^1 x_4^2 + x_1^2 x_3^3 x_4^1- x_1^1 x_3^3 x_4^2 + x_1^3 x_3^1 x_4^2- x_1^2 x_3^1 x_4^3 + x_1^1 x_3^2 x_4^3- x_1^3 x_3^2 x_4^1 \right|$$

Known relations

The expectation value of $V$ is $$\mathbb{E}[V]=\frac{2}{3}\sqrt{\frac{2}{\pi}}\tag{2}$$ A proof can be found in a Math Stack Exchange post.

The variance of $V$ is $$\mathbb{Var}[V]=\mathbb{E}[V^2]-(\mathbb{E}[V])^2=\frac{2}{3}-\frac{8}{9\pi}\tag{3}$$ where $\mathbb{E}[V^2]$ can be calculated by multiple integration.

Approximate relations based on empirical data

The remaining part contains only unproven statements that could give indications of the true solution.

The probability distribution of empirical data of $V$ can be fitted quite well with a function of the form $$f(V)=\text{exp}\left(-\left(\frac{V}{c_2}\right)^{c_1}\right)c_3\tag{4}$$ where $c_1,c_2,c_3$ are fit parameters. As a PDF must fulfill the conditions $$\int_0^\infty f(V) \mathrm{d}V=1\ \,\, \text{and}\ \int_0^\infty V f(V)\mathrm{d}V=\mathbb{E}[V]$$ the fit parameters $c_2$ and $c_3$ in eq.(4) can be expressed in dependence of $c_1$ $$c_2=\mathbb{E}[V]\frac{\Gamma(1/c_1)}{\Gamma(2/c_1)}\ ,\;\; c_3=\frac{c_1}{c_2\Gamma(1/c_1)}$$ with $\Gamma$ being the Gamma function. Only $c_1$ remains to be fitted. The best fit is for $c_1\approx\pi/4$, i.e. $c_1\approx 0.7854, c_2\approx 0.3491, c_3\approx 2.4944$. However it is not known whether eq.(4) is the true form of the PDF at all. It just models well experimental data.

Answer 1

Since the coordinates of the vertices are i.i.d. as standard normal variables, then the coordinates of the difference of two vertices are i.i.d normal variables with $0$ mean and variance $2$.

But the differences taken wrt the same vertex are no longer uncorrelated and this fact complicates the treatment.

Herewith I am providing, for the moment, the PDF-CDF of the volume of a random tetrahedron with a vertex at the origin, and the other three vertices having independent standard normal coordinates.

This Mathworld article indicates that the PDF of the product of three normal independent variables can be expressed through a Meijer G-Function. To find by this way the distribution of the sum of the triplets appearing in the algebraic definition of the triple product is therefore unviable.

However the distribution of each vertex has the nice property of spherical symmetry and we are going to exploit that via the trigonometric version of the triple product.

a) Polar distribution

Using Cartesian and Spherical (geographic convention) coordinates in parallel $$ \left\{ \matrix{ x = r\cos \phi \cos \theta \hfill \cr y = r\cos \phi \sin \theta \hfill \cr z = r\sin \phi \hfill \cr} \right.\quad \left| \matrix{ \; - \pi /2 \le \phi \le \pi /2 \hfill \cr \; - \pi < \theta \le \pi \hfill \cr} \right.\quad dx\,dy\,dz \leftrightarrow r^{\,2} \cos \phi \,dr\,d\theta \,d\phi $$ we can write the spatial distribution of a vertex as $$ \eqalign{ & p_p ({\bf V}) dV = \,{\cal N}_{\sigma ^{\,2} } \,(x){\cal N}_{\sigma ^{\,2} } \,(y){\cal N}_{\sigma ^{\,2} } \,(z)\,dx\,dy\,dz = \cr & = \left( {{1 \over {\sigma \sqrt {2\pi } }}} \right)^{\,3} e^{\, - \,{1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } r^{\,2} \cos \phi \,dr\,d\theta \,d\phi = \cr & = \left( {{1 \over {\sqrt {2\pi } }}} \right)^{\,3} e^{\, - \,{1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } \left( {{r \over \sigma }} \right)^{\,2} \,d\left( {{r \over \sigma }} \right)\,\cos \phi \,d\theta \,d\phi = \cr & = \left( {{1 \over {\sigma \sqrt {2\pi } }}} \right)^{\,2} {\cal N}_{\sigma ^{\,2} } \,(r)\,dA_{\,r} \,dr = {1 \over {2\pi }}{\cal N}_1 \,(r/\sigma )\,dA_{\,r/\sigma } \,d\left( {{r \over \sigma }} \right)\, = \cr & = \left( {{1 \over {2^{3/2} \Gamma \left( {3/2} \right)}}} \right)e^{\, - \,{1 \over 2} \left( {{r \over \sigma }} \right)^{\,2} } \left( {{r \over \sigma }} \right)^{\,2} \,d\left( {{r \over \sigma }} \right)\,\cos \phi \,d\theta \,d\phi = \cr & = {1 \over {4\pi }}\,\chi _3 \,(r/\sigma )\,d\left( {{r \over \sigma }} \right)\,\,\cos \phi \,d\theta \,d\phi = \cr & = {1 \over {2\pi }} \,\,((r/\sigma )^{\,2} )\,\left( {{r \over \sigma }} \right)d\left( {{r \over \sigma }} \right) \,\,\cos \phi \,d\theta \,d\phi = \cr & = \,{1 \over {4\pi }}\chi _3 \,(r/\sigma )\,{{dA_{\,r/\sigma } } \over {(r/\sigma )^{\,2} }}\,d\left( {{r \over \sigma }} \right) = {1 \over {4\pi }}\,\chi _3 \,(r/\sigma )\,d\Omega \,d\left( {{r \over \sigma }} \right) \cr} $$ where:

• we generalize to the case of zero mean and generic variance $\sigma ^{\,2} $;
• $dA$ is the surface area element;
• $d\Omega$ the solid angle in steradians;
• $\chi _3$ and $\chi _3^2$ are respectively the chi and chi-square distributions.

The radial distribution is instead $$ \eqalign{ & p_r (r)dr = \int_{\phi = - \pi /2}^{\pi /2} {\int_{\theta = - \pi }^\pi {\left( {{1 \over {\sigma \sqrt {2\pi } }}} \right)^{\,3} e^{\, - {1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } r^{\,2} \cos \phi \,dr\,d\theta \,d\phi } } = \cr & = 4\pi r^{\,2} \left( {{1 \over {\sigma \sqrt {2\pi } }}} \right)^{\,3} e^{\, - {1 \over 2} \left( {{r \over \sigma }} \right)^{\,2} } dr = 4\pi \left( {{r \over \sigma }} \right)^{\,2} \left( {{1 \over {\sqrt {2\pi } }}} \right)^{\,3} e^{\, - {1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } d\left( {{r \over \sigma }} \right) = \cr & = \left( {{{\left( {r/\sigma } \right)^{\,2} } \over {2^{\,1/2} \Gamma \left( {3/2} \right)}}} \right) e^{\, - {1 \over 2}\left( {{r \over \sigma }} \right)^{\,2} } d\left( {{r \over \sigma }} \right) = 2\left( {{r \over \sigma }} \right)^{\,2} {\cal N}_1 \,(r/\sigma )\,\,d\left( {{r \over \sigma }} \right) = \cr & = \chi _3 \,(r/\sigma )\,d\left( {r/\sigma } \right) = 2\left( {r/\sigma } \right)\chi _3^2 \,\,((r/\sigma )^{\,2} )d\left( {r/\sigma } \right) \cr} $$

b) Cross product

Passing to unitary variance for simplicity, we can compute the distribution of the modulus ($c$) of the cross product of two vectors (two vertices) by fixing one vector $\bf v$ of modulus $v$ and then the set of vectors having a component $c /v$ normal to the first, and thus lying over the cylindrical shell of radius $c/v, (c+dc)/v$ around $\bf v$, and integrate over $v$. That is $$ \eqalign{ & p_c (c)dc = \int\limits_v {p_r (v)\,dv{\cal N}_1 \,(c/v){{dc} \over v}{\cal N}_1 \,(0)2\pi {c \over v}} = \cr & = \int\limits_v {2v^{\,2} {\cal N}_1 \,(v)\,\,dv{\cal N}_1 \,(c/v){{dc} \over v} {\cal N}_1 \,(0)2\pi {c \over v}} = \cr & = {{4\pi } \over {\sqrt {2\pi } }}cdc\int\limits_v {{\cal N}_1 \,(v)\,\,{\cal N}_1 \,(c/v)} \;dv = \cr & = {2 \over {\sqrt {2\pi } }}cdc\int_{v = 0}^\infty {e^{\, - \,{1 \over 2}\left( {v^{\,2} + c^{\,2} /v^{\,2} } \right)} \,dv} \cr} $$ and the corresponding CDF being $$ \eqalign{ & P_c (c) = {2 \over {\sqrt {2\pi } }}\int_{t\, = \,0}^c {\int_{v\, = \,0}^\infty {t\,e^{\, - \,{1 \over 2}\left( {v^{\,2} + t^{\,2} /v^{\,2} } \right)} dt\,dv} } = \cr & = {2 \over {\sqrt {2\pi } }}\int_{v\, = \,0}^\infty {e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv \int_{t\, = \,0}^c {t\,e^{\, - \,{1 \over 2}\left( {t^{\,2} /v^{\,2} } \right)} dt\,} } = \cr & = {1 \over {\sqrt {2\pi } }}\int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv \int_{t\, = \,0}^c {2\,\left( {{t \over v}} \right)\,e^{\, - \, {1 \over 2}\left( {t^{\,2} /v^{\,2} } \right)} d\left( {{t \over v}} \right)\,} } = \cr & = {2 \over {\sqrt {2\pi } }}\int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv \int_{u/2\, = \,0}^{\left( {c^{\,2} /v^{\,2} } \right)/2} {\,e^{\, - \,{1 \over 2}u} d\left( {{u \over 2}} \right)\,} } = \cr & = \sqrt {{2 \over \pi }} \int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv \left( {1 - e^{\, - {1 \over 2}\,\left( {c^{\,2} /v^{\,2} } \right)} } \right)} = \cr & = \sqrt {{2 \over \pi }} \left( {\int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv} - \int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} + c^{\,2} /v^{\,2} } \right)} dv} } \right) \cr & = 1 - \sqrt {{2 \over \pi }} \int_{v\, = \,0}^\infty {v^{\,2} e^{\, - \,{1 \over 2}\left( {v^{\,2} + c^{\,2} /v^{\,2} } \right)} dv} \cr} $$

c) Dot product

For the modulus (absolute value) $q$ of the dot product of two vectors , we fix again a vector $\bf v$ and integrate over a plane normal to it at distance $q/v$, and double the result to include the symmetric plane, which means $$ \eqalign{ & p_d (q)dq = 2\int\limits_v {2v^{\,2} {\cal N}_1 \,(v)\,dv{\cal N}_1 \,(q/v){{dq} \over v}} = \cr & = 4dq\int\limits_v {v\,{\cal N}_1 \,(v)\,\,{\cal N}_1 \,(q/v)dv} = \cr & = {2 \over \pi }dq\int_{v\, = \,0}^\infty {v\,e^{\, - \,{1 \over 2}\left( {v^{\,2} + q^{\,2} /v^{\,2} } \right)} dv} \cr} $$ and $$ \eqalign{ & P_d \left( q \right) = \int_{t\, = \,0}^q {{2 \over \pi }dt\int_{v\, = \,0}^\infty {v\,e^{\, - \,{1 \over 2}\left( {v^{\,2} + t^{\,2} /v^{\,2} } \right)} dv} } = \cr & = {2 \over \pi }\int_{v\, = \,0}^\infty {v\,e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} dv \int_{t\, = \,0}^q {e^{\, - \,{1 \over 2}\left( {t^{\,2} /v^{\,2} } \right)} dt} } = \cr & = \sqrt {{2 \over \pi }} \int_{v\, = \,0}^\infty {v^{\,2} \,e^{\, - \,{1 \over 2}\left( {v^{\,2} } \right)} \, {\rm erf}\left( {{q \over {\sqrt 2 v}}} \right)dv} \cr} $$

d) Triple product

Now it is easy to combine the results above and obtain for the volume $v$ of a parallelepiped , so leaving apart the factor of $1/6$ $$ \eqalign{ & p_t (v)dv = 2\int\limits_c {p_c (c)dc{\cal N}_1 \,(v/c){{dv} \over c}} = \cr & = 2\int_{c\, = \,0}^\infty {{2 \over {\sqrt {2\pi } }}c\,dc{\cal N}_1 \,(v/c){{dv} \over c} \int_{t = 0}^\infty {e^{\, - \,{1 \over 2}\left( {t^{\,2} + c^{\,2} /t^{\,2} } \right)} \,dt} } = \cr & = {2 \over \pi }dv\int_{c\, = \,0}^\infty {e^{\, - \,{1 \over 2}\left( {v^{\,2} /c^{\,2} } \right)} dc \int_{t = 0}^\infty {e^{\, - \,{1 \over 2}\left( {t^{\,2} + c^{\,2} /t^{\,2} } \right)} \,dt} } \cr} $$ and the CDF $$ \eqalign{ & P_{\,t} (v) = {2 \over \pi }\int_{t = 0}^v {dt\int_{c = 0}^\infty {e^{\, - \,\,{1 \over 2}\left( {{{t^{\,2} } \over {c^{\,2} }}} \right)} dc\int_{u = 0}^\infty {e^{\, - \,{1 \over 2}\left( {u^{\,2} + c^{\,2} /u^{\,2} } \right)} \,du} } } = \cr & = {2 \over \pi }\int_{c = 0}^\infty {dc\int_{t = 0}^v {e^{\, - \,\,{1 \over 2}\left( {{{t^{\,2} } \over {c^{\,2} }}} \right)} dt \int_{u = 0}^\infty {e^{\, - \,{1 \over 2}\left( {u^{\,2} + c^{\,2} /u^{\,2} } \right)} \,du} } } = \cr & = \sqrt {{2 \over \pi }} \int_{c = 0}^\infty {\,\,c\;{\rm erf}\left( {{v \over {\sqrt 2 \,c}}} \right)dc \int_{u = 0}^\infty {e^{\, - \,{1 \over 2}\left( {u^{\,2} + c^{\,2} /u^{\,2} } \right)} \,du} } \cr} $$

which are definitely less complicated than expected, and allow in case to work out some approximations.

Note that all the CDF formulas above have been checked by numerical simulation, and besides that all of them correctly evaluate to 1 at $\infty$. Since they are obtained by exact integration of the PDF's, also the latter should be correct.

-- update --

I just discovered thanks to A0, that the inner integral above is discussed in many other posts, e.g. in this, and that it is simply $$ \int_{x = 0}^\infty {e^{\, - \,{1 \over 2}\left( {x^{\,2} + c^{\,2} /x^{\,2} } \right)} \,dx} = \sqrt {{\pi \over 2}} \,e^{\, - \,\,\sqrt {c^{\,2} } } $$ By that, for the triple product (volume of a parallelepiped) we get

$$ \eqalign{ & p_t (v)\,dv = dv\sqrt {{2 \over \pi }} \int_{t\, = \,0}^\infty {e^{\, - \,{1 \over 2}\left( {v^{\,2} /t^{\,2} + 2t} \right)} dt\,} \cr & P_{\,t} (v) = \int_{t = 0}^\infty {\,\,t\;e^{\, - \,\,t} \, {\rm erf}\left( {{v \over {\sqrt 2 \,t}}} \right)dt\,} \cr} $$

e) General distribution

Concerning a Tetrahedron (parallelepiped) in a general position, i.e. the scheme considered above plus a translation of the origin, the numerical simulation suggests that the $P_{\,t} (v)$ above converts to $P_{\,t} (2 v)$.
I am thriving to find a justification of that.

Answer 2

We'll try to solve the following similar problem. It's easy to see how the solution of such problem answers to your question.

Given $$ \mathbf{\mathrm{X}_1} =(x_1^1,x_1^2,x_1^3),\;\; \mathbf{\mathrm{X}_2}=(x_2^1,x_2^2,x_2^3),\;\; \mathbf{\mathrm{X}_3}=(x_3^1,x_3^2,x_3^3)$$ where $x_i^j$ are independent standard normal distributed variables $$x_i^j\sim\mathcal{N}(0,1)$$

What is the PDF of the volume of a parallelepiped whith vertices on $(0,0,0), \mathbf{\mathrm{X}_1}, \mathbf{\mathrm{X}_2}, \mathbf{\mathrm{X}_3}, \mathbf{\mathrm{X}_1}+\mathbf{\mathrm{X}_2}, \mathbf{\mathrm{X}_1}+\mathbf{\mathrm{X}_3}, \mathbf{\mathrm{X}_2}+\mathbf{\mathrm{X}_3}, \mathbf{\mathrm{X}_1}+\mathbf{\mathrm{X}_2}+\mathbf{\mathrm{X}_3}$?

The idea to solve this problem is to decompose each random vector in two random variables that correspond to their unit vector and their (eucidean) norm. Luckily, the variables obtained can be shown to be independent.

So, we define $$\mathbf{\mathrm{U}_i}=\mathbf{\mathrm{\hat{X}}_i}$$ $$\mathbf{\mathrm{N}_i}=||\mathbf{\mathrm{X}_i}||_2$$

We observe that the volume (with sign) $V$ of the parallelepiped can be expressed as $$V=N_1N_2N_3(U_3\cdot(U_1\times U_2))$$ The distribution of the $N_i$'s is known to be the chi distribution with parameter $k=3$.

Let's then focus on $U_3\cdot(U_1\times U_2)$. It can be shown that:

1).the angle $\theta$ between $U_1$ and $U_2$ has a uniform probability distribution on the interval $[-\pi,\pi]$.

2).the angle $\phi$ between $U_1\times U_2$ and $U_3$ has a uniform probability distribution on the interval $[-\pi,\pi]$.

3).The random variables $\phi$ and $\theta$ are independent.

Moreover, $U_3\cdot(U_1\times U_2)=\text{sin}(\theta)\text{cos}(\phi)$. It's not too hard to find that the probability distribution of $\text{sin}(\theta)$ (which is equal to that of $\text{cos}(\phi)$) is $\frac{1}{\pi \sqrt{1-x^2}}$.

I don't think we can go further than that. If we want the PDF explicity we should apply some convolution on all those density functions, and the result (if you can get one) wouldn't be pretty at all.