Covariance of Student's t-distribution Another integral (this time it looks like a lot of work but maybe it can be simplified).
I have the Student's t-distribution
$$\int_{-\infty}^\infty \frac{\Gamma(D/2+\nu/2)}{\Gamma(\nu/2)}\frac{|\Lambda|^{1/2}}{(\pi \nu)^{D/2}}\left(1+\frac{(x-\mu)^{T}\Lambda(x-\mu)}{\nu}\right)^{-D/2-\nu/2}\mathrm dx$$
and I want to calculate its covariance, so I started using the definition $\text{cov}[x] = E[(x-E[x])(x-E[x])^{T}]$ which simplifies to $E[xx^{T}]-\mu\mu^{T}$
Trying to solve the expectation $E[xx^{T}]$, we get something like this
$$E[xx^{T}]=C \int \left(1+\frac{(x-\mu)^{T}\Lambda (x-\mu)}{\nu}\right)^{-\alpha} xx^{T}dx$$
with $\alpha=D/2+\nu/2$ and $C$ as a normalization constant. 
Changing variables via $y=\Lambda^{1/2}(x-\mu)$, we obtain
$$E[xx^{T}]=C \int \frac{1}{\left(1+\displaystyle \frac{y^{T}y}{\nu}\right)^\alpha} (\Lambda^{-1/2}y+\mu) (\Lambda^{-1/2}y+\mu)^{T}|\Lambda|^{-1/2}dy$$
At this point, I know the cross terms of the product $(\Lambda^{-1/2}y+\mu) (\Lambda^{-1/2}y+\mu)^{T}$ vanish and the last term is going to produce $\mu\mu^{T}$ because the distribution is normalized, but I'm having a hard time getting the first term
$$C \int_{\infty}^{\infty} \frac{\Lambda^{-1/2}y y^{T}\Lambda^{-1/2}}{\left(1+\displaystyle \frac{y^{T}y}{\nu}\right)^\alpha}|\Lambda|^{-1/2}dy$$
I tried to use spherical coordinates but since $yy^{T}$ is a $D×D$ matrix I would need to integrate every term. There should be a better procedure.
By the way, the result I should get is $\text{cov}[x]=\displaystyle \frac{\nu}{\nu-2}\Lambda^{-1}$.
Any help is appreciated.
 A: I haven't checked you algebra, but since you are interested in integration in polar coordinates, let us say that you want to compute that last integral
$$
\int_{\mathbb{R}^n} \frac{ \Lambda^{-1/2}  y y^T \Lambda^{-1/2}}{(1+\frac{y^T y}{\nu} )^{\alpha}}|\Lambda|^{-1/2} dy = |\Lambda|^{-1/2}  \Lambda^{-1/2}  \Big\{ \int_{\mathbb{R}^n} \frac{ y y^T }{(1+\frac{y^T y}{\nu} )^{\alpha}} dy \Big\} \Lambda^{-1/2}
$$
by the linearity of the integral and that $\Lambda$ is constant. Now, in general 
$$
\int_{\mathbb{R}^n} f = n v_n \int_0^\infty \int_{S^{n-1}} f(r \theta) r^{n-1} d \sigma(\theta) dr 
$$
where $\theta$ is general point on the sphere $S^{n-1}$ and $\sigma$ is uniform measure on that sphere. Here $v_n$, is the volume of the unit ball in $\mathbb{R}^n$. (See An Elementary Introduction to Modern Convex Geometry, by K. Ball, p. 5.)
Letting $y = r\theta$ and noting that $y^Ty = r^{2}$ since $\theta^T \theta = 1$, we get
$$
\int_{\mathbb{R}^n} \frac{ y y^T }{(1+\frac{y^T y}{\nu} )^{\alpha}} dy = n v_n
\int_0^\infty \frac{r^2}{(1+r^2/\nu)^\alpha} r^{n-1} dr \int_{S^{n-1}} \theta \theta^T d\sigma(\theta).
$$
Note that $\int_{S^{n-1}} \theta \theta^T d\sigma(\theta) $ is the covariance of a random variable uniformly distributed on the sphere $S^{n-1}$. By symmetry, it is equal to some constant, say $\alpha_n$, times the identity matrix $I_n$. To obtain that constant, you can verify that each coordinate of such a random variable is distributed as a Beta random variable. Thus,
$$
\int_{\mathbb{R}^n} \frac{ y y^T }{(1+\frac{y^T y}{\nu} )^{\alpha}} dy = n v_n \alpha_n I_n
\int_0^\infty \frac{r^2}{(1+r^2/\nu)^\alpha} r^{n-1} dr
$$
PS. The distribution of $\theta$ is the same as the distribution of $\frac{X}{\|X\|}$ where $X$ is a standard normal random variable. 
EDIT: If you need a rigorous proof of that integration formula in polar coordinate, one that I know of is based on the heavy machinery of geometric measure theory. In particular, if you have heard of co-area formula, then this is a special case. For example, in "Sobolev mappings, co-area formula and related topics" by Piotr Hajlasz, the following is derived (p. 232)
$$
\int_{\mathbb{R}^n} g \,dx = \int_0^\infty \Big( 
\int_{\partial B(0,r)} g \,d \mathcal{H}^{n-1} \Big) dr
$$
as an application of that general result. Here $\mathcal{H}^{n-1}$ is the $(n-1)$-dimensional Hausdorff measure and $\partial B(0,r)$ is the boundary of the ball of radius $r$ centered at origin. Note that $H^{n-1}$ restricted to $\partial B(0,r)$ is in essence the distribution of a random variable uniformly distributed on the sphere of radius $r$. I will leave it to you to show that this version and the one I mentioned earlier are the same.
