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?
 A: 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<n$ 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=4}^n X_j^2\right)^{1/2}\theta\left(1-\sum_{j=4}^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}.$$

A: $\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. 
A: 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{\Gamma(n/2)}{\pi^{k/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)$.
A: $\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. 

A: Although great answers have already been provided, the one provided below is perhaps be the shortest possible answer to the question.

Let $x_1,\ldots,x_k \in \mathbb R$ such that $r_k^2:=\sum_{i=1}^k x_i^2 < 1$, and note that the marginal density of $(X_1,\ldots,X_k)$ at $(x_1,\ldots,x_k)$ equals $p_n(r_k^2)$, where
$$
p_n(t) \propto \int_0^\infty\ldots \int_0^\infty dx_{k+1}\ldots dx_n\delta(\sum_{i=1}^n x_i^2-t) \,\forall t \in \mathbb R.
$$
The Laplace transform w.r.t $t$ is given by
$$
\begin{split}
\hat{p}_n(s) = \int_0^\infty e^{-ts}p_n(t)dt &\propto \int_0^\infty e^{-s\sum_{i=1}^n x_i^2}dx_{k+1}\ldots dx_n \\
&= e^{-s\sum_{i=1}^k x_i^2}\int_0^\infty e^{-s\sum_{i=k+1}^n x_i^2}dx_{k+1}\ldots dx_n\\
&= e^{-sr_k^2}\left(\int_0^\infty e^{-sz^2}dz\right)^{(n-k)} \propto e^{-sr_k^2}s^{-(n-k)/2}.
\end{split}
$$
Evaluating (e.g via mathematica, etc.) the inverse Laplace transform of the last term at $t=r_k^2$, we deduce that
$p_n(r_k^2) \propto (1-t)^{(n-k)/2-1}\delta(t-r_k^2)\bigg\rvert_{t=r_k^2} = (1-r_k^2)^{(n-k)/2-1}$.
A: 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}.$$
