Likelihood Inference for Skew-Normal and Skew-T Models

Likelihoods for Skew-Normal and Skew-t Samples

Likelihood inference turns a model density into a function of its unknown parameters, evaluated at the observed data. In skew models, this involves both an elliptical density and a distribution-function factor. The latter is what allows asymmetry, and it also makes the likelihood less straightforward to optimize than its symmetric counterpart.

Use the following parameterization. For a \(p\)-dimensional observation \(x\), let \(\xi\in\mathbb R^p\), let \(\Omega\) be a positive-definite scale matrix, and set \(D=\operatorname{diag}(\sqrt{\Omega_{11}},\ldots,\sqrt{\Omega_{pp}})\). Write \(z=D^{-1}(x-\xi)\) and \(q=(x-\xi)^\mathsf{T}\Omega^{-1}(x-\xi)\). The shape parameter is \(\alpha\in\mathbb R^p\), and \(\eta=\alpha^\mathsf{T}z\). The matrix \(\Omega\) is a scale parameter; for skew distributions it is not generally the covariance matrix.

For independent observations \(x_1,\ldots,x_n\), the likelihood is the product of the individual densities and the log-likelihood is their sum. For example, the skew-normal log-likelihood is

The skew-\(t\) log-likelihood is

These expressions presume \(\Omega\) is positive definite and, for the skew-\(t\), \(\nu>0\). They specify the likelihood, but do not guarantee a finite interior maximizer or that numerical optimization will find one.

Location and Shape Score Components

A score is a derivative of the log-likelihood with respect to a parameter. For compactness, define the inverse Mills ratio \(r(u)=\phi(u)/\Phi(u)\), where \(\phi\) is the standard normal density. For the skew-\(t\), define \(h_m(u)=t_m(u)/T_m(u)\), with \(t_m\) the univariate \(t\) density. The following derivatives hold for one observation with \(\Omega\) fixed; sample scores are obtained by summing the corresponding expressions.

For the skew-normal, differentiating the Gaussian log-density with respect to \(\xi\) gives \(\Omega^{-1}(x-\xi)\). Also, \(\nabla_\xi\eta=-D^{-1}\alpha\), and the derivative of \(\log\Phi(\eta)\) with respect to \(\eta\) is \(r(\eta)\). Combining the two terms gives its location score. Since \(\nabla_\alpha\eta=z\), differentiation with respect to \(\alpha\) gives the stated shape score.

For the skew-\(t\), the derivative of its symmetric log-density with respect to \(\xi\) is \((\nu+p)\Omega^{-1}(x-\xi)/(\nu+q)\). For the skewing argument \(s=\eta k\), the chain rule gives \(\nabla_\xi s=-kD^{-1}\alpha+s\Omega^{-1}(x-\xi)/(\nu+q)\), because \(\nabla_\xi q=-2\Omega^{-1}(x-\xi)\). Multiplying this derivative by \(h_{\nu+p}(s)\), the derivative of \(\log T_{\nu+p}(s)\) with respect to \(s\), yields the skew-\(t\) location score. Finally, \(\nabla_\alpha s=kz\), giving the shape score. These are complete score components for \(\xi\) and \(\alpha\) with \(\Omega,\nu\) held fixed; scatter derivatives must also account for the dependence of \(D\), \(q\), and \(\eta\) on \(\Omega\).

Degrees of Freedom, Optimization, and Uncertainty

The skew-\(t\) degrees-of-freedom parameter affects both its symmetric density and its skewing factor. Put \(m=\nu+p\), and let \(H_m(s)=\partial\log T_m(s)/\partial m\), with \(s\) held fixed in this partial derivative. Differentiating the one-observation log-density gives

where \(\psi\) is the digamma function. The last two terms account for the dependence of the \(t\) distribution function on its degrees of freedom and the change in \(s\) with \(\nu\), respectively. In practice, \(H_m(s)\) can be evaluated using a reliable derivative routine or finite differences checked for numerical stability.

Optimization should respect the parameter constraints. A Cholesky factorization \(\Omega=LL^\mathsf{T}\), with positive diagonal entries of \(L\) represented on a log scale, maintains positive definiteness. Optimizing \(\log\nu\) maintains \(\nu>0\). The shape vector is unconstrained in this parameterization. Stable evaluation of \(\log\Phi(\eta)\) and \(\log T_m(s)\) is important in the tails; direct calculation of a tiny distribution-function value followed by a logarithm can underflow. Multiple starting values and checks of the final gradient, Hessian, and likelihood are prudent because skewness, scale, and tail thickness can be difficult to distinguish in finite samples.

At an interior maximum \(\hat\theta\), the observed information is \(J(\hat\theta)=-\nabla^2\ell(\hat\theta)\). When the model is correctly specified, the true parameter is interior, the information is nonsingular, and usual differentiability and identifiability conditions hold, likelihood theory motivates the approximate covariance \(J(\hat\theta)^{-1}\) and a local normal approximation for \(\hat\theta\). A likelihood-ratio statistic compares \(2\{\ell(\hat\theta)-\ell(\theta_0)\}\) to a chi-square reference distribution under regularity conditions for the tested parameter. Such approximations should not be used automatically for boundary or nonidentified cases; an ordinary interior chi-square approximation may fail there.

Common Pitfall and Practical Interpretation

A common mistake is to interpret a nonzero shape score at \(\alpha=0\) as proof that the model is skewed, or to assume a numerical optimizer must have found a finite maximum because it returned parameter values. The score is a local derivative for the observed sample, not a diagnosis of population asymmetry. Likewise, a large finite estimate can indicate a likelihood that continues to improve toward a boundary, poor numerical conditioning, or weak information in the data. Inspect likelihood profiles and gradients, compare multiple starts, and report boundary behavior rather than relying on a single optimizer status code.

Likelihood-based uncertainty summaries are most informative when the fitted solution is interior and the local curvature is well behaved. If the profile likelihood is markedly asymmetric or flat, a symmetric standard-error interval based on the inverse Hessian may be misleading. A profile likelihood keeps the parameter of interest fixed while reoptimizing the others, making departures from a quadratic likelihood more visible. Boundary fits and weak identification call for particular caution in interpreting Wald and likelihood-ratio approximations.