Estimating Degrees of Freedom in Multivariate T Models

Why Estimate the Degrees of Freedom?

The degrees of freedom $\nu$ control how quickly a multivariate Student t distribution places probability in its tails relative to a Gaussian distribution. Smaller $\nu$ means heavier tails; as $\nu$ grows, the t density approaches a Gaussian density. Estimating $\nu$ therefore gives a data-driven way to describe tail weight, but it requires care: the likelihood equation is nonlinear, and its inputs depend on how location and scatter are fitted.

Let $x_1,\ldots,x_n$ be independent observations in $\mathbb{R}^p$. Use $\mu$ for the location vector and $\Sigma$ for the positive-definite scatter matrix. In the t model, $\Sigma$ is not generally the covariance matrix. When $\nu>2$, the covariance is $\nu\Sigma/(\nu-2)$; when $\nu\leq 2$, the covariance does not exist.

Definition: The $p$-variate Student t model with degrees of freedom $\nu>0$, location $\mu$, and positive-definite scatter $\Sigma$ has density $$ f(x)=\frac{\Gamma((\nu+p)/2)}{\Gamma(\nu/2)(\nu\pi)^{p/2}|\Sigma|^{1/2}} \left(1+\frac{\delta(x)}{\nu}\right)^{-(\nu+p)/2}, \qquad \delta(x)=(x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu). $$ The quantity $\delta(x)$ is the squared Mahalanobis distance measured using the scatter matrix.

The Likelihood Score When Location and Scatter Are Fixed

First hold $\mu$ and $\Sigma$ fixed. Write $\delta_i=(x_i-\mu)^\mathsf{T}\Sigma^{-1}(x_i-\mu)$. Up to terms that do not depend on $\nu$, the log-likelihood is

\ell(\nu)=n\left\{\log\Gamma\left(\frac{\nu+p}{2}\right)-\log\Gamma\left(\frac{\nu}{2}\right)-\frac{p}{2}\log\nu\right\} -\frac{\nu+p}{2}\sum_{i=1}^n\log\left(1+\frac{\delta_i}{\nu}\right).

The gamma-function derivative is the digamma function, $\psi(a)=d\log\Gamma(a)/da$. Differentiating the first term gives $\frac{n}{2}\{\psi((\nu+p)/2)-\psi(\nu/2)-p/\nu\}$. For each observation,

\frac{d}{d\nu}\left[-\frac{\nu+p}{2}\log\left(1+\frac{\delta_i}{\nu}\right)\right] =-\frac12\log\left(1+\frac{\delta_i}{\nu}\right) +\frac{\nu+p}{2\nu}\frac{\delta_i}{\nu+\delta_i}.

Adding these derivatives gives the score. A finite interior maximum with fixed location and scatter must make this score zero.

Result: For fixed $\mu$ and $\Sigma$, the likelihood score for $\nu$ is $$ U(\nu)=\frac{n}{2}\left\{\psi\left(\frac{\nu+p}{2}\right)-\psi\left(\frac{\nu}{2}\right)-\frac{p}{\nu}\right\} -\frac12\sum_{i=1}^n\log\left(1+\frac{\delta_i}{\nu}\right) +\frac{\nu+p}{2\nu}\sum_{i=1}^n\frac{\delta_i}{\nu+\delta_i}. $$ Thus a finite interior likelihood estimate solves $U(\nu)=0$. A root is a stationary point, not automatically the global maximum.

Worked Example: Evaluating the Score at a Candidate Value

Take $p=1$, fix $\mu=0$ and $\Sigma=1$, and suppose the three observations are $0,1,2$. Their squared distances are $\delta_1=0,\delta_2=1,\delta_3=4$. Evaluate the score at $\nu=4$.

The digamma values are $\psi(5/2)\approx0.70316$ and $\psi(2)\approx0.42278$. The first part of the score is \[ \frac{3}{2}\left\{0.70316-0.42278-\frac14\right\}\approx0.04557. \] The log-distance contribution is \[ -\frac12\left\{\log(1)+\log(1.25)+\log(2)\right\}\approx-0.45815. \] The final contribution is \[ \frac{5}{8}\left\{0+\frac15+\frac12\right\}=0.4375. \] Consequently, $U(4)\approx0.0249$. The likelihood is increasing at this candidate value, so $\nu=4$ is not a stationary point. This calculation alone does not locate a maximum; the score must be evaluated or optimized over the admissible range.

Profiling Out Location and Scatter

In a fitted model, the distances $\delta_i$ usually change as $\nu$ changes because location and scatter are estimated from the same observations. Define the profile log-likelihood by $\ell_{\mathrm{prof}}(\nu)=\ell(\nu,\hat\mu_\nu,\hat\Sigma_\nu)$, where, for each fixed $\nu$, $(\hat\mu_\nu,\hat\Sigma_\nu)$ maximizes the likelihood over location and positive-definite scatter. This definition presumes that such an interior maximizer exists for the values being considered.

At a differentiable interior profile solution, the derivative still has the score form above, evaluated at the fitted distances. The reason is the chain rule: derivatives through $\hat\mu_\nu$ and $\hat\Sigma_\nu$ are multiplied by the location and scatter likelihood scores, which are zero at an interior maximum.

Theorem: Suppose $\ell(\nu,\mu,\Sigma)$ is differentiable in its arguments on an open parameter region, and for each $\nu$ in an interval it has a differentiable interior maximizer $(\hat\mu_\nu,\hat\Sigma_\nu)$. Then $$ \frac{d}{d\nu}\ell_{\mathrm{prof}}(\nu) = \left.\frac{\partial}{\partial\nu}\ell(\nu,\mu,\Sigma)\right|_{\mu=\hat\mu_\nu,\,\Sigma=\hat\Sigma_\nu}. $$ For the multivariate t likelihood, this is the score $U(\nu)$ with $\delta_i=(x_i-\hat\mu_\nu)^\mathsf{T}\hat\Sigma_\nu^{-1}(x_i-\hat\mu_\nu)$.

Proof. Apply the chain rule to the profile:

\frac{d}{d\nu}\ell(\nu,\hat\mu_\nu,\hat\Sigma_\nu) = \frac{\partial\ell}{\partial\nu} + \left(\nabla_\mu\ell\right)^\mathsf{T}\frac{d\hat\mu_\nu}{d\nu} + \left\langle\nabla_\Sigma\ell,\frac{d\hat\Sigma_\nu}{d\nu}\right\rangle.

The first-order conditions for an interior maximizer give $\nabla_\mu\ell=0$ and $\nabla_\Sigma\ell=0$. The last two terms therefore vanish, leaving just the partial derivative with $\mu$ and $\Sigma$ held fixed. Substituting the t log-likelihood derivative derived above proves the stated score formula. $\square$

This result does not say that one can calculate the profile score once and keep the same distances for every $\nu$. The distances must be recalculated at the fitted location and scatter for each candidate value. Nor does the result apply without qualification at a boundary or a point where the nuisance-parameter maximizer is not differentiable.

Numerical Solution in Practice

The profile equation generally has no closed-form solution. A reliable computation treats $\nu$ as a scalar outer parameter and fits location and scatter at each trial value. For fixed $\nu$, the t likelihood’s location and scatter score equations can be solved iteratively; the weights in those equations depend on the current distances, with $w_i=(\nu+p)/(\nu+\delta_i)$. After obtaining the nuisance estimates, evaluate the profile likelihood and its score.

A practical approach is to search on $\eta=\log\nu$, which enforces $\nu>0$, and use a bracketing or safeguarded method rather than an unguarded Newton iteration. A root-finding method locates stationary points; compare their profile likelihood values and inspect the limiting Gaussian model as $\nu\to\infty$. If the profile keeps increasing toward that limit, the likelihood estimate may be reported as a Gaussian-boundary solution rather than as a finite root. Multiple roots are possible, and a root by itself does not establish a maximum.

Choose a trial $\nu$.
Use a positive value, often represented numerically through $\eta=\log\nu$.

Fit location and scatter.
Maximize the likelihood over $\mu$ and positive-definite $\Sigma$ at this fixed $\nu$.

Evaluate the profile.
Compute the fitted Mahalanobis distances, then evaluate the profile likelihood and score.

Check the optimum.
Use a safeguarded search, compare candidate maxima, and assess the Gaussian limit rather than accepting any score root automatically.

A Moment-Based Comparison

A moment estimate provides a useful contrast to likelihood estimation, but it requires a finite fourth moment. For a multivariate t vector with $\nu>4$, let $C=\operatorname{Cov}(X)=\nu\Sigma/(\nu-2)$ and define the covariance-standardized squared radius $Q_C=(X-\mu)^\mathsf{T}C^{-1}(X-\mu)$. The population radial kurtosis is $E[Q_C^2]$; for a Gaussian vector it equals $p(p+2)$.

Result: For a $p$-variate t vector with $\nu>4$, $$ E[Q_C^2]=p(p+2)\frac{\nu-2}{\nu-4}. $$ Therefore, if $\beta=E[Q_C^2]$, the corresponding moment equation gives $$ \nu=4+\frac{2}{\beta/[p(p+2)]-1}, $$ provided $\beta>p(p+2)$.

Derivation. The t scale-mixture representation gives $X-\mu=\Sigma^{1/2}Z/\sqrt{V/\nu}$, where $Z$ is standard normal in $\mathbb{R}^p$, $V$ has a chi-squared distribution with $\nu$ degrees of freedom, and $Z$ and $V$ are independent. Thus the scatter-standardized squared radius is $Q=(X-\mu)^\mathsf{T}\Sigma^{-1}(X-\mu)=\nu Z^\mathsf{T}Z/V$. Since $C^{-1}=(\nu-2)\Sigma^{-1}/\nu$, we have $Q_C=(\nu-2)Q/\nu$. For $\nu>4$, the chi-squared moments give $E[(Z^\mathsf{T}Z)^2]=p(p+2)$ and $E[V^{-2}]=1/\{(\nu-2)(\nu-4)\}$. Independence then yields

E[Q_C^2] =\left(\frac{\nu-2}{\nu}\right)^2\nu^2 \,p(p+2)\frac{1}{(\nu-2)(\nu-4)} =p(p+2)\frac{\nu-2}{\nu-4}.

Dividing by $p(p+2)$ and solving $\beta/[p(p+2)]=(\nu-2)/(\nu-4)$ gives the stated formula. The condition $\nu>4$ is essential: without a finite fourth moment, this population kurtosis is not defined. $\square$

Worked Example: Inverting a Radial-Kurtosis Moment

Suppose $p=2$ and a covariance-standardized radial-kurtosis calculation gives $\hat\beta=12$. The Gaussian reference value is $p(p+2)=2(4)=8$, so the normalized excess is $12/8-1=0.5$. The moment estimate is \[ \hat\nu_{\mathrm{mom}}=4+\frac{2}{0.5}=8. \] Checking against the population formula, $\nu=8$ predicts $E[Q_C^2]=8(8-2)/(8-4)=12$, matching the supplied moment. This is an inversion of a moment equation, not a likelihood calculation; with observed data, the sample radial moment is an estimate and may be noisy.

Worked Example: Converting Scatter to Covariance

Take a bivariate t model with $\nu=8$ and $\Sigma=\begin{pmatrix}1&0\\0&4\end{pmatrix}$. Because $\nu>2$, its covariance is \[ C=\frac{8}{8-2}\Sigma=\frac43\Sigma =\begin{pmatrix}4/3&0\\0&16/3\end{pmatrix}. \] For the displacement $x-\mu=(2,0)^\mathsf{T}$, the squared radius using scatter is $Q=2^2/1=4$. The covariance-standardized squared radius is $Q_C=2^2/(4/3)=3$. Equivalently, $Q_C=(\nu-2)Q/\nu=(6/8)4=3$. Moment-based radial kurtosis uses covariance-standardized radii; confusing them with scatter-standardized distances changes the quantity being estimated.

What the Two Estimates Do—and Do Not—Share

Likelihood estimation uses the full shape of the t density, including the gamma-function normalizer, and can be applied for any $\nu>0$. The radial-kurtosis inversion uses only a fourth-moment relationship and requires $\nu>4$. In a sample, a plug-in moment estimate typically uses an estimated center and covariance to compute radial distances; it is not automatically unbiased, and its value can be unstable when the estimated excess kurtosis is close to zero. If that excess is nonpositive, the displayed inversion does not yield a finite $\nu>4$ estimate.

The likelihood and moment estimates need not agree. They weight information differently, and the moment estimate can be especially sensitive to a few extreme radii. Conversely, a likelihood fit can be sensitive to model misspecification and to optimization failures. A sensible analysis checks the fitted likelihood, the optimization, and the moment estimate’s assumptions rather than treating either number as definitive.

Common pitfall: Setting a score root equal to the likelihood estimate without checking it is a maximum is incorrect. The correction is to evaluate the profile likelihood at candidate stationary points and compare with relevant boundaries, including the Gaussian limit. A second common error is to use $E[Q_C^2]$ for a t model with $\nu\leq4$; that moment does not exist, so its inversion is not justified.

Check Your Understanding

Answer these questions using the definitions and equations above.

Why must the Mahalanobis distances in a profile score generally be recomputed when the candidate degrees of freedom changes?
What does the envelope argument remove from the derivative of the profile log-likelihood, and what conditions does that argument require?
Why is a root of the degrees-of-freedom score not, by itself, proof that the likelihood has a finite maximum there?
For what values of $\nu$ does the covariance-standardized radial-kurtosis moment exist?
How does the covariance-standardized squared radius relate to the scatter-standardized squared radius when $\nu>2$?