# Marginal density of uniform spherical distribution

Suppose that $$X$$ is distributed uniformly in the scaled $$n$$-sphere $$\sqrt{n} \mathbf{S}^{n-1} \subset \mathbf{R}^n$$. Then apparently the distribution of $$(X_1, \dots, X_k)$$, the first $$k < n$$ coordinates of $$X$$ has density $$p(x_1, \dots, x_k)$$ with respect to Lebesgue measure in $$\mathbf{R}^k$$, moreover if $$r^2 = x_1^2 + \cdots + x_k^2$$, then it is proportional to $$\left(1 - \frac{r^2}{n}\right)^{(n-k)/2 - 1}, \quad \text{if}~0 \leq r^2 \leq n,$$ and otherwise is 0. I tried to compute this using the fact that $$(X_1, \dots, X_k) \stackrel{\rm d}{=} \sqrt{n} (g_1, \dots, g_k)/\sqrt{g_1^2 + \cdots + g_n^2}$$, when $$g_i$$ are iid standard normal variables, but it was somewhat unclear to me how to use even this representation to compute the density. Can anyone sketch the details for me?

$$\newcommand{\R}{\mathbb{R}} \newcommand{\x}{\mathbf{x}} \newcommand{\X}{\mathbf{X}}$$ This is to present a formalization of the answer by Carlo Beenakker, without explicit use of the delta function.

We are going to assume that $$\X=(X_1,\dots,X_n)$$ is uniformly distributed on the unit sphere $$\mathbb S^{n-1}$$, rather than on $$\sqrt n\,\mathbb S^{n-1}$$.

For each real $$t\in(0,1)$$, define the measures $$\mu_t$$ and $$\nu_t$$ over $$\R^n$$ by the conditions $$\begin{equation*} \int f\,d\mu_t=\int_{\R^n}d\x f(\x)1_{1-t<|\x|^2\le1} \end{equation*}$$ and $$\begin{equation*} \int f\,d\nu_t=\frac{\int f\,d\mu_t}{\int d\mu_t} \end{equation*}$$ for all (say) nonnegative continuous functions $$f\colon\R^n\to\R$$, where $$|\cdot|$$ denotes the Euclidean norm. Then $$\nu_t$$ is a probability measure converging (as $$t\downarrow0$$) to the Haar measure $$h$$ on the unit sphere in $$\R^n$$, in the sense that (say) $$\begin{equation*} \int f\,d\nu_t\to\int f\,dh \end{equation*}$$ for all nonnegative continuous functions $$f\colon\R^n\to\R$$.

Take now any function $$f\colon\R^n\to\R$$ such that $$\begin{equation*} f(\x)=g(\x_{n-1}) \end{equation*}$$ for some nonnegative continuous function $$g\colon\R^{n-1}\to\R$$ and all $$\x\in\R^n$$, where $$\x_j:=(x_1,\dots,x_j)$$ for $$\x=(x_1,\dots,x_n)\in\R^n$$ and $$j=1,\dots,n-1$$. Then
\begin{align*} \int f\,d\mu_t&=\int_{\R^{n-1}}d\x_{n-1}\,g(\x_{n-1})\int_\R du\, 1_{1-t<|\x_{n-1}|^2+u^2\le1} \\ &=\int_{\R^{n-1}}d\x_{n-1}\,g(\x_{n-1}) (1+o(1))t\,(1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{align*} Also, $$\begin{equation*} \int d\mu_t=\int_{\R^n}d\x\, 1_{1-t<|\x|^2\le1}\propto(1+o(1))t, \end{equation*}$$ where $$\propto$$ means an equality up to a constant factor, depending only on $$n,k$$. So, $$\begin{equation*} \int f\,dh=\lim_{t\downarrow0} \int f\,d\nu_t\propto\int_{\R^{n-1}}d\x_{n-1}\,g(\x_{n-1}) (1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{equation*}$$ Thus, the joint pdf of $$\X_{n-1}=(X_1,\dots,X_{n-1})$$ is given by $$\begin{equation*} p_{n-1}(\x_{n-1})\propto(1-|\x_{n-1}|^2)^{-1/2}\,1_{|\x_{n-1}|<1}. \end{equation*}$$ Now successively integrating $$p_{n-1}(\x_{n-1})$$ ($$n-1-k$$ times) in $$x_{n-1},\dots,x_{k+1}$$ and each time using the formula $$\begin{equation*} \int_0^{b^{1/2}}(b-u^2)^p du=c_p b^{p+1/2} \end{equation*}$$ for real $$b>0$$ and $$p>-1$$ with $$c_p:=\int_0^1(1-u^2)^pdu\in(0,\infty)$$ (so that $$1/2$$ is added to the exponent $$p$$ after such an integration), we see that the joint pdf of $$\X_k=(X_1,\dots,X_k)$$ is given by $$\begin{equation*} p_k(\x_k)\propto(1-|\x_k|^2)^{(n-k)/2-1}\,1_{|\x_k|<1}, \end{equation*}$$ as desired.

• This is a wonderfully simple calculation and was illuminating. Thank you for the follow up questions, your persistence, and both answers. May 8 '20 at 20:35
• @DrewBrady : Thank you for your kind words. However, I feel grateful to Carlo for teaching me something of value, with patience. May 8 '20 at 21:03

With $$X$$ uniformly distributed over the unit $$n$$-sphere, the joint probability distribution of all $$n$$ elements of $$X$$ is a Dirac delta function, $$P(X_1,X_2,\ldots X_n)\propto\delta\left(1-\sum_{j=1}^n X_j^2\right).\qquad\qquad(1)$$ Now you integrate out elements one by one, to obtain the marginal distribution $$P_k$$ of $$k elements. The first integration gives $$P_{n-1}(X_2,X_3,\ldots X_n)\propto\left(1-\sum_{j=2}^n X_j^2\right)^{-1/2}\theta\left(1-\sum_{j=2}^n X_j^2\right),\qquad(2)$$ with $$\theta$$ the unit step function. The second integration gives $$P_{n-2}(X_3,\ldots X_n)\propto\theta\left(1-\sum_{j=3}^n X_j^2\right),$$ the third integration $$P_{n-3}(X_4,\ldots X_n)\propto\left(1-\sum_{j=2}^n X_j^2\right)^{1/2}\theta\left(1-\sum_{j=2}^n X_j^2\right),$$ and so on. Each additional integration increases the power by 1/2, $$P_{k}(X_{n-k+1},\ldots X_n)\propto\left(1-\sum_{j=n-k+1}^n X_j^2\right)^{(n-k)/2-1}\theta\left(1-\sum_{j=n-k+1}^n X_j^2\right).$$ This is the answer in the OP (without rescaling the radius of the $$n$$-sphere, so $$r^2/n\mapsto r^2$$).

As requested in the comments, a more detailed exposition of the various steps.

First step: the delta function. Denote the surface measure on the unit $$n$$-sphere as $$d\Omega$$, and $$\int d\Omega=A_n$$ the surface area. Uniformity of a distribution on the unit $$n$$-sphere means uniformity with measure $$d\Omega$$. I maintain that the joint probability distribution of the components of the vector $${\mathbf X}=(X_1,X_2,\ldots X_n)$$, uniformly distributed on the unit $$n$$-sphere, is given by Eq. (1) with normalization constant $$2/A_n$$. Let us check this by calculating the expectation value of an arbitrary function $$f$$ of $${\mathbf X}$$. For that purpose I transform to hyperspherical coordinates $$r,\phi_1,\phi_2,\ldots\phi_{n-1}$$, $$\mathbb{E}[f(\mathbf{X})]=\int dX_1\int dX_2\cdots\int dX_n \,f(X_1,X_2,\ldots X_n)P(X_1,X_2,\ldots X_n)$$ $$\qquad=\int_0^\infty r^{n-1} dr \int d\Omega\, f(r,\phi_1,\phi_2,\ldots\phi_{n-1})\frac{2}{A_n}\delta(1-r^2)$$ $$\qquad=\frac{1}{A_n}\int d\Omega\, f(r=1,\phi_1,\phi_2,\ldots\phi_{n-1}).$$ In the last step I used that $$\int_0^\infty r^{n-1}\delta(1-r^2)\,dr=1/2$$ for $$n\geq 2$$.

Second step: integration of the delta function, to arrive at Eq. (2). From now on I will ignore the normalization constants, these can easily be recovered at the end. Let me abbreviate $$\sum_{j=2}^n X_j^2=s_2$$. The marginal distribution $$P_1(X_2,X_3,\ldots X_n)$$ is obtained by definition upon integration of $$P(X_1,X_2,X_3,\ldots X_n)$$ over $$X_1$$. I carry out this integration in cartesian coordinates, changing variables to $$q=X_1^2$$, $$P_1(X_2,X_3,\ldots X_n)\propto \int_{-\infty}^\infty dX_1\delta(1-s_2-X_1^2),$$ $$\qquad=\int_0^\infty\delta(1-s_2-q)\frac{dq}{\sqrt q}=(1-s_2)^{-1/2}\theta(1-s_2).$$

Third and following steps: The following steps, subsequent integrations of $$X_2,X_3,\ldots$$ are now immediate consequences of the integral $$\int_0^a(a^2-x^2)^p\,dx=c_p a^{1+2p}.$$

• Hm. What does it mean that this distribution "is a delta function"? May 7 '20 at 16:24
• @ Carlo Beenakker - are you sure you understand what a joint distribution is?
– R W
May 7 '20 at 18:51
• certainly, let's work this out for $n=2$; we then have in polar coordinates $x_1=r\cos\phi$, $x_2=r\sin\phi$ and we would expect the uniform distribution on the unit circle to be $P(\phi)=\text{constant}$. So let's check that the delta function distribution indeed gives this: take $P(x_1,x_2)\propto\delta(1-x_1^2-x_2^2)$ and integrate to obtain $P(\phi)=\int_0^\infty P(r\cos\phi,r\sin\phi) rdr\propto\int_0^\infty \delta(1-r^2)rdr=\text{constant}$. The multivariate Gaussian is a convenient way to generate the random coordinates, but for this calculation the delta function is easier. May 7 '20 at 21:02
• @ Carlo Beenakker - how does the RHS of your formula know that you are dealing with the uniform distribution on the sphere?
– R W
May 7 '20 at 21:19
• Carlo, I'd really like to understand your solution, especially because it promises an improvement over the Gaussian approach. Alas, I don't understand almost any part of your answer. Especially hard for me to understand your step-by-step integration. Can you write down in detail the integrals for at least the first two steps and how you take them? May 7 '20 at 23:29

Let $$G_1,\dots,G_n$$ be iid standard normal random variables. Then the random vector $$\begin{equation*} (Y_1,\dots,Y_n):=\Big(\frac{G_1}{\sqrt{\sum_1^n G_j^2}},\dots, \frac{G_n}{\sqrt{\sum_1^n G_j^2}}\Big) \end{equation*}$$ is uniformly distributed on the unit sphere $$\mathbb S^{n-1}$$. Let \begin{equation*} \begin{aligned} Z_i&:=Y_i=\frac{G_i}{\sqrt{\sum_1^n G_j^2}}&\text{ if }i\le k,\\ Z_i&:=G_i&\text{ if }i> k. \end{aligned} \end{equation*} We want to find the joint pdf of $$(Z_1,\dots,Z_k)$$, which is the same as the joint pdf of $$(X_1,\dots,X_k)/\sqrt n$$.

The vector $$(Z_1,\dots,Z_n)$$ is obtained from $$(G_1,\dots,G_n)$$ by the transformation given by \begin{equation*} \begin{aligned} z_i&:=\frac{g_i}{\sqrt{\sum_1^n g_j^2}}&\text{ if }i\le k,\\ z_i&:=g_i&\text{ if }i> k. \end{aligned} \end{equation*} The transformation inverse to this is given by \begin{equation*} \begin{aligned} g_i&:=\sqrt{s_2}\frac{z_i}{\sqrt{1-s_1}}&\text{ if }i\le k,\\ g_i&:=z_i&\text{ if }i> k, \end{aligned} \tag{1} \end{equation*} where $$\begin{equation*} s_1:=\sum_1^k z_j^2,\quad s_2:=\sum_{k+1}^n z_j^2. \end{equation*}$$ The Jacobian determinant of the inverse transformation is $$\begin{equation*} J=\det(cM)=c^k\det M, \end{equation*}$$ where $$\begin{equation*} c:=s_2^{1/2}(1-s_1)^{-3/2},\quad M:=(1-s_1)I_k+UU^T, \end{equation*}$$ $$I_k$$ is the $$k\times k$$ identity matrix, and $$U:=[z_1,\dots,z_k]^T$$.

Write $$U=|U|Qe_1$$, where $$|U|=\sqrt{s_1}$$ is the Euclidean norm of $$U$$, $$Q$$ is some orthogonal matrix, and $$e_1:=(1,0,\dots,0)$$. Then it is clear that the matrix $$M$$ is similar to $$N:=(1-s_1)I_k+|U|^2e_1e_1^T=(1-s_1)I_k+s_1e_1e_1^T$$, whence $$\det M=\det N=(1-s_1)^{k-1}$$. So, $$\begin{equation*} J=s_2^{k/2}(1-s_1)^{-k/2-1}. \tag{2} \end{equation*}$$ Also, the joint pdf of $$(G_1,\dots,G_n)$$ is given by $$\begin{equation*} (2\pi)^{-n/2}\exp\Big\{-\frac12\sum_1^n g_j^2\Big\}. \end{equation*}$$ So, in view of (1) and (2), the joint pdf of $$(Z_1,\dots,Z_n)$$ is given by $$\begin{equation*} f_n(z_1,\dots,z_n) =(2\pi)^{-n/2}\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\}s_2^{k/2}(1-s_1)^{-k/2-1}. \end{equation*}$$ So, the joint pdf of $$(Z_1,\dots,Z_k)$$ is given by \begin{align*} f_k(z_1,\dots,z_k)&=\int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,f_n(z_1,\dots,z_n) \\ &=(2\pi)^{-n/2}(1-s_1)^{-k/2-1} \\ &\times \int_{\mathbb R^{n-k}}dz_{k+1}\dots dz_n\,s_2^{k/2}\,\exp\Big\{-\frac12\frac{s_2}{1-s_1}\Big\} \\ &\propto(1-s_1)^{(n-k)/2-1}, \end{align*} because $$s_2=\sum_{k+1}^n z_j^2$$. So, we have the desired result.

More specifically, \begin{align*} f_k(z_1,\dots,z_k) =\frac{2^{(k-n)/2}\Gamma(n/2)}{\pi^{n/2}\Gamma((n-k)/2)}(1-s_1)^{(n-k)/2-1} \end{align*} for $$s_1=\sum_1^k z_j^2\in(0,1)$$.

$$\newcommand{\R}{\mathbb{R}} \newcommand{\x}{\mathbf{x}} \newcommand{\X}{\mathbf{X}}$$ Here is yet another solution, which is partly informal but I think not hard to completely formalize. Its advantage is a strong and hopefully convincing appeal to geometric intuition.

Again, we are going to assume that $$(X_1,\dots,X_n)$$ is uniformly distributed on the unit sphere $$\mathbb S^{n-1}$$, rather than on $$\sqrt n\,\mathbb S^{n-1}$$. Let $$\X_j:=(X_1,\dots,X_j)$$ and $$\x_j:=(x_1,\dots,x_j)$$ for $$\x=(x_1,\dots,x_n)\in\R^n$$ and $$j=1,\dots,n-1$$. Let $$|\cdot|$$ denotes the Euclidean norm.

The main point is that the probability density $$p_{n-1}(\x_{n-1})$$ of $$\X_{n-1}$$ at a point $$\x_{n-1}\in\R^{n-1}$$ with $$|\x_{n-1}|<1$$ is proportional to the ratio $$r_{n-1}(\x_{n-1}):=vol_{n-1}(dS)/vol_{n-1}(dA)$$, where $$vol_{n-1}$$ is of course the $$(n-1)$$-volume, $$dA$$ is an infinitesimal neighborhood of the point $$\x_{n-1}$$ in $$\R^{n-1}$$ and $$dS$$ is the preimage of $$dA$$ under the projection of the upper hemisphere $$\mathbb S^{n-1}_+:=\{\x\in\mathbb S^{n-1}\colon\x\cdot e_n\ge0\}$$ onto the closed unit ball in $$\R^{n-1}$$; this projection is given by $$\mathbb S^{n-1}_+\ni(\x_{n-1},u)\mapsto\x_{n-1}$$; here $$e_n:=(0,\dots,0,1)$$ and $$\cdot$$ denotes the dot product. But $$\begin{equation} r_{n-1}(\x_{n-1})=\frac{vol_{n-1}(dS)}{vol_{n-1}(dA)}=\frac1{\cos\phi}, \end{equation}$$ where $$\phi$$ is the angle between the hyperplane $$\R^{n-1}\times\{0\}$$ of $$\R^n$$ and the tangent hyperplane to $$\mathbb S^{n-1}$$ at the point $$(\x_{n-1},\sqrt{1-|\x_{n-1}|^2})\in\mathbb S^{n-1}_+$$; that is, $$\phi$$ is the angle between the corresponding normal vectors $$e_n$$ and $$(\x_{n-1},\sqrt{1-|\x_{n-1}|^2})$$ of these two hyperplanes.

Thus, $$\begin{equation} p_{n-1}(\x_{n-1})\propto r_{n-1}(\x_{n-1})\propto\frac1{\cos\phi}=(1-|\x_{n-1}|^2)^{-1/2}, \end{equation}$$ where $$\propto$$ means an equality up to a constant factor, depending only on $$n,k$$.

Now successively integrating $$p_{n-1}(\x_{n-1})$$ ($$n-1-k$$ times) in $$x_{n-1},\dots,x_{k+1}$$ and each time using the formula $$\begin{equation*} \int_0^{b^{1/2}}(b-u^2)^p du=c_p b^{p+1/2} \end{equation*}$$ for real $$b>0$$ and $$p>-1$$ with $$c_p:=\int_0^1(1-u^2)^pdu\in(0,\infty)$$ (so that $$1/2$$ is added to the exponent $$p$$ after such an integration), we see that the joint pdf of $$\X_k=(X_1,\dots,X_k)$$ is given by $$\begin{equation*} p_k(\x_k)\propto(1-|\x_k|^2)^{(n-k)/2-1}\,1_{|\x_k|<1}, \end{equation*}$$ as desired.

Here is a picture, for $$n=3$$, showing the upper hemisphere $$\mathbb S^{n-1}_+$$; a small neighborhood of a point $$\x_{n-1}$$ in the projection of $$\mathbb S^{n-1}_+$$ to the horizontal plane $$\R^{n-1}\times\{0\}$$ of $$\R^n$$; the preimage of that neighborhood under that projection; and the normal vectors of the horizontal plane and the tangent plane to the sphere -- with $$\phi$$ being the angle between these two normal vectors. I enjoyed thinking about these answers and this is my attempt to put them into (nonrigorous) geometrical terms. Writing the joint density compositionally as

$$p(\mathbf{x}_k \mid |\mathbf{x}| = 1)p(\mathbf{x}_{n-k} \mid \mathbf{x}_k, |\mathbf{x}| = 1) = p(\mathbf{x} \mid |\mathbf{x}| = 1) \propto 1,$$

we want to solve for the first term on the left. But since our density is proportional to a constant, this is just

$$p(\mathbf{x}_k \mid |\mathbf{x}| = 1) \propto \frac{1}{p(\mathbf{x}_{n-k} \mid \mathbf{x}_k, |\mathbf{x}| = 1)}.$$

Accordingly, instead of performing our calculation by integrating out $$X_{k+1} \dots X_{n}$$, we can think about it in terms of the conditional density for sampling $$\mathbf{X}_{n-k}$$, given $$\mathbf{X}_k$$ and the norm constraint $$|\mathbf{X}| = 1$$, denoted in the denominator above.

I propose a two step procedure. First, draw a point uniformly from within the $$n-k-1$$ dimensional ball with radius $$r_{n-k} = \sqrt{1 - |\mathbf{x}_k|^2}.$$ Each such vector has density proportional to $$\left(1 - |\mathbf{x}_k|^2\right)^{-(n-k-1)/2}.$$

This corresponds to the second and higher integrations in the previous answers, whereas here we directly use the formula for the volume of a ball.

Next, the $$(n-k)$$th coordinate must satisfy $$|\mathbf{x}| = 1$$, which is achieved by any point on the circle with radius $$r_{n-k}$$, a set with measure proportional to $$r_{n-k}$$. Proving this is the first integration in the previous answers.

Putting these two steps together and taking the reciprocal gives

$$p(\mathbf{x}_k \mid |\mathbf{x}| = 1) \propto \left(1 - |\mathbf{x}_k|^2\right)^{(n-k)/2 - 1}.$$