Likelihood Contributions in an Elliptical Model
In an elliptical model, observations far from the proposed center need not contribute to the likelihood in the same way as observations close to it. Their effect depends on the radial generator: the function governing how density changes with squared Mahalanobis distance. This tutorial derives how that dependence enters the likelihood scores for location and scatter.
Suppose \(X_1,\ldots,X_n\) are independent observations in \(\mathbb{R}^p\), each with the same density \(f(x;\mu,\Sigma)=|\Sigma|^{-1/2}g((x-\mu)^\mathsf{T}\Sigma^{-1}(x-\mu))\). Here \(\mu\in\mathbb{R}^p\) is the location parameter, \(\Sigma\) is a symmetric positive-definite scatter matrix, and the known generator \(g\) is positive and differentiable at the distances considered. Assume \(g\) is normalized so this is a probability density. The scatter parameter is not automatically the covariance matrix: their relationship depends on the radial distribution.
The phrase “squared Mahalanobis distance” matters: \(d_i\) depends on both the proposed center and the scatter matrix. A point can be far in ordinary Euclidean distance but not far relative to a scatter matrix that assigns substantial scale in that direction.
Deriving the Two Score Equations
Ignoring no parameter-dependent terms, the log-likelihood is \(\ell(\mu,\Sigma)=-\frac{n}{2}\log|\Sigma|+\sum_{i=1}^n\log g(d_i)\). The location and scatter scores are coupled because each \(d_i\), and therefore each weight, depends on both parameters.
Derivation of the location score. Write \(r_i=x_i-\mu\). For a location perturbation \(d\mu\), the scatter matrix is fixed and \(d d_i=-2r_i^\mathsf{T}\Sigma^{-1}d\mu\). Consequently, \(d\log g(d_i)=(g'(d_i)/g(d_i))d d_i =w_i r_i^\mathsf{T}\Sigma^{-1}d\mu\). Summing over observations identifies the gradient with respect to \(\mu\) as \(\sum_i w_i\Sigma^{-1}r_i\). Setting this vector to zero gives the first equation.
Derivation of the scatter score. Perturb \(\Sigma\) by a symmetric matrix \(d\Sigma\), holding \(\mu\) fixed. The determinant differential is \(d\log|\Sigma|=\operatorname{tr}(\Sigma^{-1}d\Sigma)\). Also, \(d d_i=-\operatorname{tr}(\Sigma^{-1}r_i r_i^\mathsf{T}\Sigma^{-1}d\Sigma)\). Thus the differential of the log-likelihood is \[ d\ell=\operatorname{tr}\left[\left(-\frac{n}{2}\Sigma^{-1} +\frac{1}{2}\sum_{i=1}^n w_i\Sigma^{-1}r_i r_i^\mathsf{T}\Sigma^{-1}\right)d\Sigma\right]. \] The matrix in brackets is symmetric. Because the differential must vanish for every symmetric perturbation at an interior stationary point, that matrix must be zero. Multiplying the resulting equation on the left and right by \(\Sigma\) gives \(\sum_i w_i r_i r_i^\mathsf{T}=n\Sigma\). This proves both stated score equations. \(\square\)
When \(\sum_i w_i\ne 0\), the location equation can also be written as a weighted mean: \(\mu=\sum_i w_i x_i/\sum_i w_i\). The weights are not fixed in advance: they must be evaluated at the distances formed using the parameters being estimated. The location and scatter equations therefore generally have to be solved together.
Worked Examples
Worked Example: The Gaussian Score Equations
For a \(p\)-dimensional Gaussian model, \(g(t)=(2\pi)^{-p/2}\exp(-t/2)\). Its logarithmic derivative is \(g'(t)/g(t)=-1/2\), so \(w(t)=1\) at every distance. The score equations reduce to \(\sum_i(x_i-\mu)=0\) and \(\sum_i(x_i-\mu)(x_i-\mu)^\mathsf{T}=n\Sigma\).
Take the three two-dimensional observations \(x_1=(0,0)^\mathsf{T}\), \(x_2=(2,0)^\mathsf{T}\), and \(x_3=(0,2)^\mathsf{T}\). The location equation gives \(\hat\mu=(x_1+x_2+x_3)/3=(2/3,2/3)^\mathsf{T}\). Their residuals are \((-2/3,-2/3)^\mathsf{T}\), \((4/3,-2/3)^\mathsf{T}\), and \((-2/3,4/3)^\mathsf{T}\).
The sum of their outer products is \[ \begin{pmatrix}4/9&4/9\\4/9&4/9\end{pmatrix} +\begin{pmatrix}16/9&-8/9\\-8/9&4/9\end{pmatrix} +\begin{pmatrix}4/9&-8/9\\-8/9&16/9\end{pmatrix} = \begin{pmatrix}8/3&-4/3\\-4/3&8/3\end{pmatrix}. \] Dividing by \(n=3\), the scatter score gives \(\hat\Sigma=\begin{pmatrix}8/9&-4/9\\-4/9&8/9\end{pmatrix}\). Its eigenvalues are \(4/9\) and \(4/3\), so it is positive definite. These are the Gaussian maximum-likelihood estimates for this full-rank sample; the covariance estimate uses divisor \(n\), not \(n-1\).
Worked Example: Student t Weights at a Candidate Parameter
For a \(p\)-dimensional Student t model with known degrees of freedom \(\nu>0\), the generator is proportional to \(g(t)\propto(1+t/\nu)^{-(\nu+p)/2}\). Differentiating its logarithm gives \(g'(t)/g(t)=-(\nu+p)/(2(\nu+t))\), hence \(w(t)=(\nu+p)/(\nu+t)\).
Set \(p=2\), \(\nu=4\), and consider \(x_1=(1,0)^\mathsf{T}\), \(x_2=(0,2)^\mathsf{T}\) at the candidate \(\mu=(0,0)^\mathsf{T}\), \(\Sigma=I\). The squared distances are \(d_1=1\) and \(d_2=4\). Their weights are \(w_1=6/(4+1)=6/5\) and \(w_2=6/(4+4)=3/4\).
The weighted residual sum for the location equation is \((6/5)(1,0)^\mathsf{T}+(3/4)(0,2)^\mathsf{T}=(6/5,3/2)^\mathsf{T}\), which is not zero. The weighted outer-product sum for the scatter equation is \((6/5)\begin{pmatrix}1&0\\0&0\end{pmatrix} +(3/4)\begin{pmatrix}0&0\\0&4\end{pmatrix} =\begin{pmatrix}6/5&0\\0&3\end{pmatrix}\). But \(n\Sigma=2I\), so the scatter equation also fails. Thus this candidate is not stationary. The calculation illustrates how each observation’s distance sets its weight, but the weights alone do not establish a fitted solution.
Worked Example: Exponential-Power Radial Weighting
Consider the known generator \(g(t)=C\exp(-t^2/2)\), where \(C\) is its normalizing constant for the chosen dimension. This generator is integrable in every finite dimension. Since \(g'(t)/g(t)=-t\), its radial weight is \(w(t)=2t\).
For \(x_1=(1,0)^\mathsf{T}\) and \(x_2=(0,2)^\mathsf{T}\), evaluate the scores at \(\mu=0\) and \(\Sigma=I\). The distances are \(1\) and \(4\), giving weights \(2\) and \(8\). The weighted location sum is \(2(1,0)^\mathsf{T}+8(0,2)^\mathsf{T}=(2,16)^\mathsf{T}\). The weighted outer-product sum is \(2\begin{pmatrix}1&0\\0&0\end{pmatrix} +8\begin{pmatrix}0&0\\0&4\end{pmatrix} =\begin{pmatrix}2&0\\0&32\end{pmatrix}\). The location sum is nonzero, and the outer-product sum differs from \(n\Sigma=2I\); this candidate is not stationary. Compared with the Student t example, the more distant point receives a larger weight here. Radial weighting is determined by the chosen generator, not by a universal rule that distant points must be downweighted.
A Useful Identity and What the Equations Mean
At a stationary point, the scatter equation implies a scalar identity that can be used to check calculations. Premultiply its matrix equation by \(\Sigma^{-1}\) and take the trace. Since \(\operatorname{tr}(\Sigma^{-1}r_i r_i^\mathsf{T})=r_i^\mathsf{T}\Sigma^{-1}r_i=d_i\), the left side becomes \(\sum_i w_i d_i\). The right side is \(n\operatorname{tr}(I_p)=np\). Therefore \[ \sum_{i=1}^n w(d_i)d_i=np. \] This identity is a necessary consequence of the scatter score equation, not an additional independent equation.
In the Gaussian case, \(w(d_i)=1\), so the identity says that the sum of fitted squared Mahalanobis distances is \(np\). For other generators, the distances are modified by the radial weights. This relation helps reveal a common error: using the correct-looking weighted outer products but omitting the factor \(n\) on the right would imply the wrong trace identity.
Common Pitfall and Practical Interpretation
- Treating the weights as constants. The weights depend on \(d_i\), and \(d_i\) depends on both \(\mu\) and \(\Sigma\). A proposed solution must evaluate the distances and weights at that same proposed parameter pair.
- Calling every score root a maximum. The equations identify interior stationary points. They do not prove that a root exists, is unique, or maximizes the likelihood; those questions require additional analysis.
- Confusing scatter with covariance. The matrix in the density is the scatter parameter. It equals the covariance only under an appropriate radial moment condition and scaling; do not substitute covariance notation without checking that relationship.
- Assuming all non-Gaussian generators downweight distant observations. The Student t generator in the example assigns smaller weights at larger distances, whereas the exponential-power generator shown assigns larger weights. The derivative of the specified generator determines the score.
The practical value of these score equations is that they separate the geometry of the data from the radial behavior of the model. The residual outer product \(r_i r_i^\mathsf{T}\) records direction and magnitude, while \(w(d_i)\) determines how strongly that observation contributes under the chosen generator. This is the core calculation behind likelihood-based fitting of elliptical location and scatter models.
Check Your Understanding
Use the definitions and score equations above to answer the following questions.
- For a generator \(g(t)\), what is the radial score weight \(w(t)\), and which quantity supplies its argument for observation \(i\)?
- Why do the location and scatter equations generally need to be solved together?
- For the Student t generator, explain how \(w(t)=(\nu+p)/(\nu+t)\) changes as \(t\) increases.
- What scalar identity follows by taking the trace of the stationary scatter equation?
- Why does satisfying both score equations not, by itself, establish that a parameter pair is a likelihood maximum?