An elliptically symmetric angular Gaussian distribution

We define a distribution on the unit sphere Sd-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {S}^{d-1}$$\end{document} called the elliptically symmetric angular Gaussian distribution. This distribution, which to our knowledge has not been studied before, is a subfamily of the angular Gaussian distribution closely analogous to the Kent subfamily of the general Fisher–Bingham distribution. Like the Kent distribution, it has ellipse-like contours, enabling modelling of rotational asymmetry about the mean direction, but it has the additional advantages of being simple and fast to simulate from, and having a density and hence likelihood that is easy and very quick to compute exactly. These advantages are especially beneficial for computationally intensive statistical methods, one example of which is a parametric bootstrap procedure for inference for the directional mean that we describe.


Introduction
A natural way to define a distribution on the unit sphere S d−1 is to embed S d−1 in R d , specify a distribution for a random variable z ∈ R d , then consider the distribution of z either conditioned to lie on, or projected onto, S d−1 . The general Fisher-Bingham and angular Gaussian distributions, defined respectively in Mardia (1975) and Mardia and Jupp (2000) can both be constructed this way by taking z to be multivariate Gaussian in R d . Then the Fisher-Bingham distribution is the conditional distribution of z conditioned on z = 1, and the angular Gaussian is the distribution of the projection z/ z . The choice of the mean, μ, and covariance matrix, V , of z controls the concentration and the shape of the contours of the induced probability density on S d−1 .
It is usually not practical to work with the general Fisher-Bingham or angular Gaussian distributions, however, because they have too many free parameters to be identified well by data. This motivates working instead with subfamilies that have fewer free parameters and stronger symmetries.
In the spherical case, d = 3, the general distributions have 8 free parameters. Respective subfamilies with 3 free parameters are the Fisher and the isotropic angular Gaussian (IAG) distributions. Both are "isotropic" in the sense that they are rotationally symmetric about the mean direction, i.e., contours on the sphere are small circles centred on the mean direction. Respective subfamilies with 5 free parameters are the Bingham and the central angular Gaussian distributions, both of which are antipodally symmetric.
An important member of this Fisher-Bingham family is the Kent distribution (Kent, 1982). For d = 3, it has 5 free parameters, and it has ellipse-like contours on the sphere. This offers a level of complexity well suited to many applications, since the distribution is flexible enough to model anisotropic data yet its parameters can usually be estimated well from data. To our knowledge, nobody to date has considered its analogue in the angular Gaussian family. The purpose of this paper is to introduce such an analogue, which we call the elliptically symmetric angular Gaussian (ESAG), and establish some of its basic properties.
The motivation for doing so is that in some ways the angular Gaussian family (and hence ESAG) is much easier to work with than the Fisher-Bingham family (and hence the Kent distribution). In particular, simulation is easy and fast, not requiring rejection methods (which are needed for the Fisher-Bingham family Kent et al. 2013), and the density is free of awkward normalising constants, so the likelihood can be computed quickly and exactly. Hence in many modern statistical settings the angular Gaussian family is the more natural choice; see for example Presnell et al. (1998) who use it in a frequentist approach for circular data, and Wang and Gelfand (2013) and Hernandez-Stumpfhauser et al. (2017) who use it in Bayesian approaches for circular and spherical data, respectively.
In the following section, we introduce ESAG, first for general d before specialising to the case d = 3.

The general angular Gaussian distribution
The angular Gaussian distribution is the marginal distribution of the directional component of the multivariate normal distribution. Let denote the multivariate normal density in R d with d ×1 mean vector μ, d × d covariance matrix V , assumed non-singular, and where |V | denotes the determinant of V . Then, writing z = r y, where r = z = (z z) 1/2 and y = z/ z ∈ S d−1 , and using dz = r d−1 dr dy, where dy denotes Lebesgue, or geometric, measure on the unit sphere S d−1 , and integrating r over r > 0, leads to where C d = 1/(2π) (d−1)/2 and, for real α, where φ(·) and (·) are the standard normal probability density function and cumulative density function, respectively; more details about the M d (α) are given in Sect. A.1.

An elliptically symmetric subfamily
The subfamily of (2) that we shall call the elliptically symmetric angular Gaussian distribution (ESAG) is defined by the two conditions under which the angular Gaussian density (2) simplifies to From (5), the positive definite matrix V has a unit eigenvalue. If the other eigenvalues are then the inverse of V can be written where ξ 1 , …, ξ d−1 and ξ d = μ/ μ is a set of mutually orthogonal unit vectors. Moreover, constraint (6) means that Once the d parameters in μ are fixed, then from (5) and (6) there are d − 2 remaining degrees of freedom for the eigenvalues of V , and d(d−1)/2−(d−1) degrees of freedom for its unit eigenvectors. The total number of free parameters is thus (d − 1)(d + 2)/2, the same as for the multivariate normal in a tangent space R d−1 to the sphere. Condition (5) imposes symmetry about the eigenvectors of V . Without loss of generality, suppose that the eigenvectors are parallel to the coordinates axes; that is, each element of the vector ξ j equals 0 except the jth which equals 1. Then if y = (y 1 , . . . , y d ) , In this case, the density (7) depends only on y j through y 2 j for j = 1, . . . , d − 1. Consequently, the density is invariant with respect to sign changes of the y 1 , . . . , y d−1 , that is, which implies reflective symmetry about 0 along the axes defined by ξ 1 , …, ξ d−1 . This type of symmetry is implied by ellipse-like contours of constant density inscribed on the sphere, and such contours arise when the density (7) is unimodal. Whether the density is unimodal depends on the nature of the stationary point at y = μ/ μ , which is characterised by the following proposition.
Proposition 1 Write α = μ and assume without loss of generality that (8) The proof of Proposition 1 is given in Appendix A. 2. We conjecture that if the stationary point at y = ξ d is a local maximum, then it is also a global maximum, and that in this case the distribution is unimodal; and in the case that the stationary point is a local minimum, then the distribution is bimodal. A rigorous proof appears difficult, but the conjecture is strongly supported by some extensive numerical investigations that we have performed and describe as follows. For each of a wide variety of combinations of ESAG parameters with d = 3 -in particular, choosing for (α, γ 1 , γ 2 ) values on a 9 × 9 × 9 rectangular lattice on [0.2, 20] × [−5, 5] × [−5, 5] -we performed numerical maximisation to find Using the Manopt implementation (Boumal et al., 2014) of the trust region approach of Absil et al. (2007), for each (α, γ 1 , γ 2 ) , we performed the optimisation multiple times from distinct starting points with the rationale that if the distribution is indeed multimodal, then optimisations from different starting points will converge to different local optima. We chose to use 42 different starting points since this enabled the points to be exactly equi-spaced on S 2 using the method of Teanby (2006). For each (α, γ 1 , γ 2 ) , we hence computed y max 1 , . . . , y max 42 . Instances that converged to the same mode had values of y max that were not quite identical, owing to the finite tolerance of the numerical optimisation. To account for this, we identified the number of modes according to clustering of the {y max i }, designating the distribution unimodal if and only if y max In cases identified as non-unimodal by this criterion, we used k-means clustering to identify k = 2 clusters; in each such case, every y max i was within a distance 10 −6 of its cluster centre indicating bimodality. In agreement with the conjecture, amongst the 9 3 = 729 parameter cases we considered, in every 553 cases with ρ d−1 ≤ H d (α), the foregoing procedure identified the distribution to be unimodal, and in every 176 cases with The next proposition concerns the limiting distribution of a sequence of unimodal ESAG distributions as the sequence becomes more highly concentrated. Without loss of generality, we fix ξ d = (0, . . . , 0, 1) , take ξ 1 , . . . , ξ d−1 to be the other coordinate axes and define Neglecting the hemisphere defined by y d negative is no drawback when considering α = μ → ∞ as follows, because in this limit the distribution becomes increasingly concentrated about y = ξ d .
Proposition 2 Assume the conditions (5) and (6), and therefore (10), hold, where the ρ j are assumed to be fixed, and suppose that α = μ → ∞. Then Remark 1 In the general case, we replace the coordinate vectors ξ 1 , . . . , ξ d by an arbitrary orthonormal basis, and then the limit distribution lies in the vector subspace spanned by ξ 1 , . . . , ξ d−1 .
Remark 2 Proposition 2 is noteworthy because it is atypical for high-concentration limits within the angular Gaussian family to be Gaussian.

A parameterisation of ESAG for d = 3
An important practical question is how to specify a convenient parameterisation for the matrix V so that it satisfies the constraints (5) and (6). With d = 3, such a V has two free parameters.
Henceforth, we will use this parameterisation. For a random variable y with ESAG distribution, we will write y ∼  Remark 4 (A test for rotational symmetry.) For a sample of observations y 1 , . . . , y n assumed independent and identically distributed, a standard large-sample likelihood ratio test can be used to test H 0 : y i ∼ IAG(μ) vs H 1 : y i ∼ ESAG(μ, γ ). Letl 0 andl 1 be the values of the maximised log-likelihoods under H 0 and H 1 , respectively. The models are nested and differ by two degrees of freedom, and by Wilks' theorem, when n is large the statistic T = 2(l 1 −l 0 ) has approximately a χ 2 2 distribution if H 0 is true, and H 0 is rejected for large values of T . The null distribution can alternatively be approximated using simulation by the parametric bootstrap, that is, by simulating a sample of size n from the null model H 0 at the maximum likelihood estimate of the parameters, computing the test statistic, T , and then repeating this a large number, say B, times. The empirical distribution of the resulting bootstrapped statistics T * 1 , . . . , T * B approximates the null distribution of T . where I μ,γ = −E μ,γ ∂ 2 log f 1 /∂μ∂γ . Moreover, the directional and magnitudinal components μ and γ , that is, μ , μ/ μ , γ and γ / γ , are all mutually orthogonal. The proofs follow easily from the symmetries of f ESAG and are omitted.
Often in applications, there is particular interest in the directional mean, m = μ/ μ . A parametric bootstrap procedure to construct confidence regions for m, which exploits both the ease of simulation and the parameter orthogonality, is as follows.

Parametric bootstrap confidence regions for m = μ/ μ
Let l(μ, γ ) denote the log-likelihood for a sample y 1 , . . . , y n each of which is assumed to be an independent ESAG(μ, γ ) random variable and define Denote the maximum likelihood estimate of (μ, γ ) by (μ,γ ), letm =μ/||μ||,ˆ = (μ,γ ) and letξ 1 andξ 2 denote the maximum likelihood estimates of the axes of symmetry (15), such thatm,ξ 1 andξ 2 are mutually orthogonal unit vectors. Define the matrixξ = (ξ 1 ,ξ 2 ) , then test statistic is suitable for defining confidence regions for m as m ∈ S 2 : T (m) ≤ c ; see Fisher et al. (1996) for discussion of test statistics of this form in the context of non-parametric bootstrap procedures. For a given significance level, α, the constant c can be determined as follows: simulate B bootstrap samples each of size n from ESAG(μ,γ ), and hence with m =m, and for each sample compute the test statis-tic (19), withξ,μ,ˆ replaced by corresponding quantities calculated from the bootstrap sample; then c is the (1 − α) quantile of the resulting statistics T * 1 (m), . . . , T * B (m). Examples of confidence regions calculated by this algorithm are shown in Fig. 2 (right).

An example: estimates of the historic position of Earth's magnetic pole
This dataset contains estimates of the position of the Earth's historic magnetic pole calculated from 33 different sites in Tasmania (Schmidt 1976). The data are shown in Fig. 2  Twice the difference in maximised log-likelihoods equals 14.12, which when referred to a χ 2 2 distribution (see Remark 4) corresponds to a p-value of less than 10 −3 , indicating strong evidence in favour of ESAG over the IAG. Figure 3 shows contours and transects of the densities of ESAG and Kent distributions. The parameter values for each are computed by fitting the two models to a large sample of independent and identically distributed data from ESAG(μ,γ ), with μ = (0, 0, 2.5) and γ = (0.75, 0) . For the inner contours ESAG is more anisotropic than the matched Kent distribution and appears slightly more peaked at the mean. Besides these small differences, the figure shows that ESAG and Kent distributions are very similar distributions in this example, as we have found them to be more generally. Indeed, preliminary results, not presented here, suggest that for typical sample sizes it is usually very difficult to distinguish between them using a statistical criterion. This warrants making the modelling choice between using the Kent distribution or the ESAG on grounds of practical convenience. The Kent distribution is a member of the exponential family, but its density involves a non-closed-form normalising constant, and simulation requires a rejection algorithm (Kent et al. 2013). The ESAG distribution has a density that is less tidy than the Kent density, hence less suited to computing moment estimators, etc., but this is not much of a drawback given that its density can be computed exactly so that the exact likelihood can be easily maximised. Moreover, The results are for three different cases, involving pairs of ESAG and Kent distributions with parameters "matched" in the sense that the were fitted to an initial common set of data. For each combination of model and parameter set, we simulated b = 500 Monte Carlo samples of n = 100 observations. The measures of error ofm,ξ 1 andξ 2 are defined in the text. Simulation times, given in seconds, are cumulative over the b = 500 Monte Carlo runs simulating from ESAG is particularly quick and easy (see Remark 3). Table 1 shows the results of a simulation study, including computational timings, for fitting ESAG and Kent densities to ESAG and Kent simulated data. To approximate the Kent normalising constant when fitting the Kent distribution, we use a saddlepoint approximation method (see Kume et al. 2013), and for simulating from the Kent distribution, we use the rejection method of Kent et al. (2013). A notion of accuracy of the fitted model is how well the mean direction of the fitted model,m, corresponds with the population mean direction, m. A measure we use for this is

A comparison of ESAG with the Kent distribution
with the expectation approximated by Monte Carlo; hence in the tables, we report wherem (i) is the mean direction of the fitted model for the ith run out of b Monte Carlo runs. We also consider accuracy of the major and minor axes of the fitted model. Since the signs ofξ 1 andξ 2 are arbitrary, in this case we define and similar forξ 2 . Note that in interpreting the results in Table 1, the different simulation times of ESAG and Kent should be compared across columns, whereas the fitting times and accuracies should be compared across rows.
The results show, as expected, that the accuracy ofm,ξ 1 , andξ 2 is typically better when the data-generating model is fitted. However, the accuracy is not dramatically worse when the non-data-generating model is fitted, i.e., when ESAG is fitted to Kent data, or the Kent distribution is fitted to ESAG data. There is a very notable difference in computation times between ESAG and Kent: for both simulation and fitting, ESAG is typically more than an order of magnitude faster than Kent.

Conclusion
In the pre-computer days of statistical modelling, the Fisher-Bingham family was perhaps favoured over the angular Gaussian family on account of having a simpler density, which makes it more amenable to constructing classical estimators such as moment estimators. However, in the era of computational statistics, the less simple form of the angular Gaussian density is hardly a barrier and is more than compensated by having a normalising constant that is trivial to evaluate. The likelihood can consequently be computed quickly and exactly, and maximised directly. Wang and Gelfand (2013) have recently argued in favour of the general angular Gaussian distribution as a model for Bayesian analysis of circular data. For spherical data, a major obstacle to using the general angular Gaussian distribution is that its parameters are poorly identified by the data. The ESAG subfamily overcomes this problem, and is a direct analogy of the Kent subfamily of the general Fisher-Bingham distribution. Besides having a tractable likelihood, the ease and speed with which ESAG can be simulated makes it especially well suited to methods of simulation-based inference. Natural wider applications of ESAG include using it as an error distribution for spherical regression models with anisotropic errors; for classification on the sphere (as a model for class-conditional densities); and for clustering spherical data (based on ESAG mixture models). Code written in MATLAB for performing calculations in this paper is available at the second author's web page.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A.1 Properties of the function M d (α)
Integrating (3) by parts, we obtain Moreover, differentiating (3) with respect to α, and exchanging the order of differentiation and integration on the RHS, we obtain where a prime denotes differentiation; and, using (21) and (22) together, it is found that Further points to note, which are easily established by induction, are that we can write where P d (α) is a polynomial of order d and Q d is a polynomial of order d − 1; and P d (α) inherits the properties (21) and (23) (21) but not (23); instead of (23), we have Q d (α) = Q d+1 (α) − P d (α). One relevant property of P d (α) is that the coefficient of its leading term α d is 1; this follows from either (21) or (23) applied repeatedly to P d (α). As an easy consequence, we obtain a result that is used below: for any fixed integer d ≥ 0,
As a consequence, we can restrict attention to y of the form y = (0, . . . , u, (1 − u 2 ) 1/2 ) where u is in a neighbourhood of 0. Substituting this form of y into (7) and differentiating, it is seen that ∂ log f 1 ∂u (u) u=0 = 0, and, using (22), it is found that where δ = (ρ d−1 − 1)/ρ d−1 . The point u = 0, which corresponds to y = ξ d , is a maximum or minimum depending on whether (25) is negative or positive and, using this fact, Proposition 2.1 follows after some elementary further manipulations.
Consequently as the ρ j are remaining fixed, it is easy to see that y V −1 y → 1.