Density Contours and Modes of Skew-Normal Models

Why the Mode and Contours Need Separate Attention

For a multivariate normal distribution, the mode is its location parameter, and its density contours are centered ellipsoids. A skew-normal density behaves differently: its contours are asymmetric, and the mode generally moves away from the location parameter. The shape matrix still informs the geometry, but it does not by itself describe either the contour displacement or the covariance of the distribution.

A useful way to analyze this geometry is to transform the density to whitened coordinates. There, the skewness acts in one distinguished direction, and the mode can be found by solving a one-dimensional equation. Mapping the result back to the original coordinates clarifies why the mode displacement need not point along a principal axis of the shape matrix.

Density and Whitened Coordinates

Let $X$ be a $p$-dimensional skew-normal random vector with location $\xi\in\mathbb R^p$, positive-definite shape matrix $\Omega$, and shape vector $\alpha\in\mathbb R^p$. Write $\omega=\operatorname{diag}(\sqrt{\Omega_{11}},\ldots,\sqrt{\Omega_{pp}})$. The Azzalini skew-normal density is

f_X(x)=2\phi_p(x;\xi,\Omega)\, \Phi\!\left(\alpha'\omega^{-1}(x-\xi)\right),

where $\phi_p(x;\xi,\Omega)$ is the $N_p(\xi,\Omega)$ density and $\Phi$ is the standard normal distribution function. Here $\Omega$ is the matrix in the Gaussian kernel; when $\alpha\ne 0$, it is generally not the covariance matrix of $X$.

Use the symmetric positive-definite square root of $\Omega$ to define $y=\Omega^{-1/2}(x-\xi)$, $c=\omega^{-1}\alpha$, and $b=\Omega^{1/2}c$. Put $\lambda=\|b\|$. The standardized density, up to the change-of-variables factor, is

f_Y(y)=2\phi_p(y)\Phi(b'y).

Thus the skewing factor depends on $y$ only through its projection onto $b$. When $\lambda>0$, let $e=b/\lambda$. Any orthogonal transformation of $y$ that fixes $e$ leaves the density unchanged. In the original coordinates this becomes an affine symmetry around the line through $\xi$ in direction $\Omega^{1/2}e=\Omega c/\lambda$. It need not be ordinary Euclidean rotational symmetry.

Definition: A density contour is a level set on which the density has a fixed value. For the skew-normal density in whitened coordinates, these contours are exactly the level sets of $H(y)=\|y\|^2-2\log\Phi(b'y)$. When $b\ne0$, their distinguished symmetry axis maps to the line through $\xi$ in direction $\Omega c$.

Indeed, the logarithm of the whitened density is a constant minus $H(y)/2$. If $t=e'y$ and $r=y-te$, then $b'y=\lambda t$, $\|y\|^2=t^2+\|r\|^2$, and

H(y)=t^2+\|r\|^2-2\log\Phi(\lambda t).

At a fixed $t$, the contour's cross-section in the directions perpendicular to $e$ depends only on $\|r\|$. But the term involving $\Phi(\lambda t)$ treats positive and negative $t$ differently. Consequently, the contours are generally asymmetric along the axis even though their cross-sections have a simple radial form.

Worked Example: Comparing Densities on Opposite Sides of the Axis

Take $p=2$, $\xi=0$, $\Omega=I_2$, and $\alpha=(1,0)'$. Then $b=(1,0)'$, so the symmetry axis is the first coordinate axis. Compare points on the axis at equal distances from the origin.

At $y=(1,0)'$, the density is $2\phi_2(y)\Phi(1)$. At $y=(-1,0)'$, it is $2\phi_2(y)\Phi(-1)$. The Gaussian factors are equal because both points have squared norm 1. Therefore the density ratio is

\frac{f_Y(1,0)}{f_Y(-1,0)} =\frac{\Phi(1)}{\Phi(-1)} \approx\frac{0.8413}{0.1587} \approx 5.30.

The equal-distance points do not have equal density: the contour shape is asymmetric along the skewness direction. In contrast, the points $(0,1)'$ and $(0,-1)'$ have the same density, since both have $b'y=0$ and the same squared norm. This checks the symmetry in directions perpendicular to the axis.

Locating the Mode

For fixed $t=e'y$, the log density is largest when $r=0$, because its only dependence on $r$ is the term $-\|r\|^2/2$. The mode must therefore lie on the symmetry axis. Along that axis, maximizing the log density reduces to maximizing $h(t)=-t^2/2+\log\Phi(\lambda t)$. Define the inverse Mills ratio $m(z)=\phi(z)/\Phi(z)$, which is well-defined for every finite real $z$.

Theorem: If $b\ne0$, the skew-normal density has a unique mode. Its whitened coordinate is $y_{\text{mode}}=s_*b/\lambda$, where $s_*$ is the unique solution of $$ s_*=\lambda\frac{\phi(\lambda s_*)}{\Phi(\lambda s_*)}. $$ In the original coordinates, $$ x_{\text{mode}}=\xi+\frac{s_*}{\lambda}\Omega c. $$ If $b=0$, the density is normal and its mode is $y=0$, or $x=\xi$.

Proof. Differentiating the axis log density gives $h'(t)=-t+\lambda m(\lambda t)$. Any stationary point therefore solves $t=\lambda m(\lambda t)$, which is the stated equation for $s_*=t$. The right side is positive, so any solution has $t>0$.

To prove uniqueness, differentiate $m$: $m'(z)=-m(z)\{z+m(z)\}$. For $z\ge0$, $z+m(z)>0$. For $z<0$, use $\phi'(u)=-u\phi(u)$ to obtain $\phi(z)=\int_{-\infty}^z(-u)\phi(u)\,du>(-z)\Phi(z)$; hence $m(z)>-z$, again giving $z+m(z)>0$. Thus $m'(z)<0$ for every $z$. It follows that $h''(t)=-1+\lambda^2m'(\lambda t)<0$, so $h'$ is strictly decreasing.

At $t=0$, $h'(0)=\lambda m(0)>0$. As $t$ tends to infinity, $m(\lambda t)$ tends to zero, so $h'(t)$ tends to minus infinity. Continuity and strict decrease imply that $h'$ has exactly one zero. Since the density tends to zero as $\|y\|$ tends to infinity and, for each fixed $t$, is maximized at $r=0$, this stationary point gives the unique global mode. Finally $y_{\text{mode}}=(s_*/\lambda)b$, so $x_{\text{mode}}=\xi+\Omega^{1/2}y_{\text{mode}} =\xi+(s_*/\lambda)\Omega c$. $\square$

Worked Example: Computing a Mode When the Skewness Is Stronger

Let $\Omega=I_2$, $\xi=0$, and $\alpha=(2,0)'$. Then $b=(2,0)'$, $\lambda=2$, and the mode has the form $(s_*,0)'$. Its coordinate solves

s_*=2\frac{\phi(2s_*)}{\Phi(2s_*)}.

At $s=0.52$, the right side is $2\phi(1.04)/\Phi(1.04)\approx0.546$, which exceeds $0.52$. At $s=0.54$, it is $2\phi(1.08)/\Phi(1.08)\approx0.518$, which is less than $0.54$. The unique root is therefore between these values; numerical solution gives $s_*\approx0.531$. Hence the mode is approximately $(0.531,0)'$, displaced from $\xi=0$ in the positive skewness direction. Strict concavity in the proof ensures that this numerical root identifies the unique mode.

Contour Displacement Versus Covariance Elongation

The axis and mode describe density geometry; they should not be mistaken for covariance eigenvectors. For this parameterization, the moments of the multivariate skew-normal give $\operatorname{Cov}(X)=\Omega-(2/\pi)\delta\delta'$, where $\delta=\Omega c/\sqrt{1+\lambda^2}$. Thus skewness changes covariance as well as shifting the mode, but these are distinct effects. The contour axis is determined by $\Omega c$, while covariance principal directions are the eigenvectors of the resulting covariance matrix.

Worked Example: A Correlated Bivariate Model

Take $\Omega=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}$, $\xi=0$, and $c=(1,0)'$. This $\Omega$ is positive definite because its leading diagonal entry is 4 and its determinant is $4-1.2^2=2.56>0$. Since $\omega=\operatorname{diag}(2,1)$, this choice corresponds to $\alpha=\omega c=(2,0)'$.

The skewness strength is $\lambda^2=c'\Omega c=4$, so $\lambda=2$, and the scalar mode coordinate is the same $s_*\approx0.531$ as in the preceding example. The symmetry-axis direction in the original coordinates is $\Omega c=(4,1.2)'$. The mode is therefore

x_{\text{mode}} =\frac{0.531}{2}(4,1.2)' \approx(1.062,0.319)'.

This is a displacement along the affine symmetry axis, whose slope is $1.2/4=0.30$. It is not the principal direction of the Gaussian kernel: the largest-eigenvalue direction of $\Omega$ has slope approximately $0.351$.

The covariance is also different from $\Omega$. Here $\delta=(4,1.2)'/\sqrt5$, so

\operatorname{Cov}(X) =\Omega-\frac{2}{5\pi} \begin{pmatrix}16&4.8\\4.8&1.44\end{pmatrix} \approx \begin{pmatrix} 1.963&0.589\\ 0.589&0.817 \end{pmatrix}.

Its largest-eigenvalue direction has slope approximately $0.423$, different from both the symmetry-axis slope $0.30$ and the Gaussian-kernel principal-direction slope $0.351$. The mode displacement, kernel contour orientation, and covariance elongation are related through the model parameters, but none should be substituted for another.

Common Pitfall

A frequent error is to assume that skew-normal contours are simply ellipses centered at the mean, or that the mode must lie along a principal eigenvector of $\Omega$. The density formula contradicts both assumptions: the factor $\Phi(b'y)$ changes the density along the $b$-direction, and whitening maps that direction back to $\Omega c$, which need not be an eigenvector. In addition, $\Omega$ is the Gaussian-kernel shape matrix, not generally the covariance.

The reliable procedure is to inspect the density in whitened coordinates, identify the skewness direction $b$, solve the one-dimensional mode equation, and then transform back. For questions about variability or elongation measured by second moments, use the covariance matrix and its eigenvectors instead.

Key takeaway: Skew-normal contours are symmetric in cross-sections perpendicular to a distinguished affine axis but asymmetric along it. The unique mode lies on that axis, generally away from the location parameter; covariance-driven elongation is a separate geometric feature.

Check Your Understanding

Use the density and mode results to answer the following.

In whitened coordinates, which scalar projection of $y$ determines the skewing factor?
For $b\ne0$, why must the mode have zero component perpendicular to $b$?
What equation determines the mode coordinate $s_*$, and why does it have only one solution?
In the correlated bivariate example, why is the original-coordinate symmetry-axis direction $\Omega c$ rather than $c$?
Explain why the largest-eigenvalue direction of the covariance matrix need not coincide with the contour symmetry axis.