TL;DR
This paper introduces a new method to compute the normalization constant of the 8-parameter Fisher-Bingham distribution on the sphere, enabling more accurate modeling of directional data beyond existing subclasses.
Contribution
The paper derives an infinite sum representation for the Fisher-Bingham distribution's normalization constant, facilitating efficient computation and improved data fitting.
Findings
Normalization constant expressed as an infinite sum of hypergeometric functions
Enhanced fit of directional data using FB8 over FB5
Demonstrated computational feasibility of the new method
Abstract
The Fisher-Bingham distribution () is an eight-parameter family of probability density functions (PDF) on that, under certain conditions, reduce to spherical analogues of bivariate normal PDFs. Due to difficulties in computing its overall normalization constant, applications have been mainly restricted to subclasses of , such as the Kent () or von Mises-Fisher (vMF) distributions. However, these subclasses often do not adequately describe directional data that are not symmetric along great circles. The normalizing constant of can be numerically integrated, and recently Kume and Sei showed that it can be computed using an adjusted holonomic gradient method. Both approaches, however, can be computationally expensive. In this paper, I show that the normalization of can be expressed as an infinite sum…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
The 8-parameter Fisher-Bingham distribution on the sphere
Tianlu Yuan
Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center
University of Wisconsin
Madison, WI 53706
Abstract
The Fisher-Bingham distribution () is an eight-parameter family of probability density functions (PDF) on that, under certain conditions, reduce to spherical analogues of bivariate normal PDFs. Due to difficulties in computing its overall normalization constant, applications have been mainly restricted to subclasses of , such as the Kent () or von Mises-Fisher (vMF) distributions. However, these subclasses often do not adequately describe directional data that are not symmetric along great circles. The normalizing constant of can be numerically integrated, and recently Kume and Sei showed that it can be computed using an adjusted holonomic gradient method. Both approaches, however, can be computationally expensive. In this paper, I show that the normalization of can be expressed as an infinite sum consisting of hypergeometric functions, similar to that of the . This allows the normalization to be computed under summation with adequate stopping conditions. I then fit the to a synthetic dataset using a maximum-likelihood approach and show its improvements over a fit with the more restrictive distribution.
K****eywords Directional statistics Fisher-Bingham distribution Kent distribution Von Mises-Fisher distribution
1 Introduction
Directional statistics involves the study of probability density functions (PDF) with supports on . This paper will focus on distributions on the sphere, . Such distributions have found applications in fields as varied as earthquake modeling to paleomagnetism of lava flows to reconstruction of radio pulses [1, 2]. A simple and commonly used PDF on that is the analogue to an isotropically distributed, bivariate normal distribution is the von Mises-Fisher (vMF) distribution [3]. A more general distribution that is the analogue to a general bivariate normal distribution is the Kent () distribution [4],
[TABLE]
where is the normalization constant, is a unit vector on , and are non-negative parameters, and are unit vectors that correspond to the columns of a orthogonal matrix, , which determines the orientation of the PDF. An additional constraint, , is required to interpret the distribution as an analogue of the general bivariate normal distribution [4], as shown in the right panel of Fig. 1. The vMF distribution corresponds to the trivial case of . Both the vMF and distributions have been well-studied, and have found use in many applications. However, suffers from the restriction that its PDFs must be symmetric across two great circles intersecting at on the sphere. Data that clusters along small circles, for example the non-equatorial lines of latitude, are ill-described by .
An alternative to is the small-circle distribution proposed in [5]. This is a four-parameter subclass () of the Fisher-Bingham family and can be written as,
[TABLE]
With , describes small-circle distributions on the sphere as shown in the left panel of Fig. 1. Removing the constraint on and , the only difference between and is the sign of the term in square-brackets. The distribution is completely specified by four parameters, since is a unit vector. However, it is only a good description of data that is evenly distributed along a small circle. Generalizations are needed in order to model data that falls between the extremes described by and .
In order to perform a maximum likelihood fit of directional data, the PDFs need to be normalized. It was shown in [5] that can be written in terms of the confluent hypergeometric function. It was shown in [4] that can be written as an infinite sum consisting of modified Bessel functions of the first kind. These approaches motivated the calculation of the normalization discussed in Section 3, but first a natural generalization of and is given in Section 2.
2 The Fisher-Bingham distribution
It is simple to construct a 6-parameter PDF () that generalizes and [6]. This is given here as
[TABLE]
with , and . Clearly, corresponds to Eq. 1 and to Eq. 2. For , has a single maximum, corresponding to where is aligned with . For , can describe either the small-circle distribution of [5] or a bimodal distribution where the modes are degrees apart on a small circle as shown in the left panel of Fig. 2. As is the case for , the distribution is symmetric, which means that it cannot describe distributions that lie along small circles with a unique mode.
In order to describe unimodal distributions that lie along small circles, [7] proposed a distribution that is a natural combination of the vMF and distributions. This can be generalized further by combining the vMF and distributions, which results in the Fisher-Bingham distribution (). It is parametrized here as,
[TABLE]
where is a unit vector on the sphere, and an example is shown in the right panel of Fig. 2. The first term in Eq. 5 is proportional to a vMF distribution with mean direction aligned along , and the second term is proportional to Eq. 4. Thus, does not necessarily have to be symmetric about a great circle. When Eq. 6 reduces to Eq. 4.
3 Calculating the and normalizations
3.1 Exact series solution
The distribution was proposed in [8, 9], though it has not been widely applied due to difficulties in computing [1]. An exact calculation involving holonomic functions was given in [10], which requires solving ordinary differential equations with the Runge-Kutta method. It can also be estimated using numerical integration. Both of these methods, however, can be computationally expensive. A faster approximation given in [11] relies on the saddlepoint method, but this is known to be inexact [10]. A fast and accurate calculation of is desirable to perform maximum likelihood inference using Eq. 6. This is the subject of this Section and Section 4.
Since simply enacts a rotation of the sphere, the normalization constants above are independent of and it is simpler to work in the standard frame, where the coordinate axes are defined by the columns of with the -axis corresponding to [11]. The coordinate transformation allows us to write Eq. 6 as,
[TABLE]
In spherical coordinates
[TABLE]
and
[TABLE]
where
[TABLE]
Taylor expanding gives,
[TABLE]
The integration proceeds as in [4] by applying Eq. (6.2.1) and Eq. (9.6.18) in [12]. Noting that the integral over vanishes unless and are both even,
[TABLE]
where denotes the gamma function, the beta function, the modified Bessel function of the first kind, the confluent hypergeometric limit function, and the Gaussian hypergeometric function. The last equality can be derived from Eq. (9.6.47) and Eq. (15.4.1) in [12]. These special functions can be evaluated numerically, and can be computed to good approximation with adequate stopping conditions on , , and . By setting , the only nonzero terms occur when and Eq. 16 simplifies to
[TABLE]
the normalization for . Further setting and applying Eq. (15.1.21) in [12] recovers as computed in [4].
3.2 Closed-form approximation for
If or is large, can be approximated piecewise in two separate regimes: and . For , note that we can rewrite
[TABLE]
where . In order to factorize the integration over and , we make the assumption that the term modulated by can be fixed to , where the maximum of occurs in latitude. Then, using Eq. (13.1.27) in [12],
[TABLE]
The last line uses an approximation for large [5].
In the case of , is maximal at and for large becomes approximately a bivariate normal distribution. By setting and and Taylor expanding in , we see that
[TABLE]
which is similar to Eq. (3.5) in [4]. A comparison of Eq. 17 to Eqs. 23 and 24 is shown in the right panel of Fig. 3 (dashed blue). This approximation is accurate away from . The saddlepoint approximation of [11] is also shown (dotted orange).
These approximations are useful when working with large or , where it may not be possible to numerically compute Eq. 17 due to extremely large terms in the summand. With maximum likelihood estimation, for example, it is often simpler to work with , which can be approximated with Eqs. 23 and 24 without running into computational overflows. For , unfortunately, no closed-form approximation is known, and the options are to perform numerical integration, use the method proposed in [10], or use Eq. 16. The saddlepoint method was tested to be accurate for but not for .
4 Numerical computation
The infinite series in Eqs. 16 and 17 can be evaluated by truncation under an appropriate stopping condition. To simplify the notation, let be the summand of Eq. 16 such that . For , is non-negative for all while for , is guaranteed to be non-negative only for even . As such, may be an alternating sequence in . Furthermore, does not in general decrease monotonically in any of the indices, but does so only after a certain point.
The most efficient algorithm to estimate would be to reindex and order such that for all , and perform the summation starting from up to certain tolerance. However, this ordering is difficult to evaluate. In practice, a robust calculation can be obtained by setting a step size and defining
[TABLE]
This ensures that an even number of terms are summed at each step and . Then for some tolerance , the stopping algorithm is
, , = 0
while True:
$K$, $P_{L,K}$, $S_{L}$ = 0
**while** True:
$J$, $P_{L,K,J}$, $S_{L,K}$ = 0
**while** True:
$\tilde{c}_{8}$ += $A_{L,K,J}$
$S_{L}$ += $B_{L,K,J}$
$S_{L,K}$ += $B_{L,K,J}$
**if** $B_{L,K,J}<|\tilde{c}_{8}|\epsilon$ **and** $B_{L,K,J}\leq P_{L,K,J}$:
**break**
$P_{L,K,J}$ = $B_{L,K,J}$
$J$ += 1
**if** $S_{L,K}<|\tilde{c}_{8}|\epsilon$ **and** $S_{L,K}\leq P_{L,K}$:
**break**
$P_{L,K}$ = $S_{L,K}$
$K$ += 1
**if** $S_{L}<|\tilde{c}_{8}|\epsilon$ **and** $S_{L}\leq P_{L}$:
**break**
$P_{L}$ = $S_{L}$
$L$ += 1
where denotes the series calculation of . This nested summation routine loops through , , and in that order and breaks once the contribution to the partial sum of for the current index is within the tolerance and is less than the previous term. In tests, and seemed to perform well. A comparison of computed using the series summation and a numerical integration routine (QUADPACK) is shown in the left panel of Fig. 3.
The evaluation of Eq. 17 follows the same procedure as above, but with a single summation over while setting . This amounts to just running the innermost loop. One final thing to note is that often contains large terms in its numerator and denominator that may often cancel each other. Computationally it may be difficult to evaluate them separately, and a better approach is to work with , explicitly taking the logarithm before exponentiation.
5 Example application
To illustrate the performance of the distribution in modeling directional data, a synthetic dataset was randomly sampled from an distribution that peaked along a small circle on the sphere. An unbinned maximum likelihood fit was then performed using the and distributions [4, 1]. The SLSQP routine was used to perform a constrained fit of and the L-BFGS-B routine was used for a bounded fit of [13].
The results are shown in Fig. 4. The synthetic data (black points) is poorly described by the distribution, but well described by . The best-fit negative log-likelihoods are for the and for the , indicating much better data agreement using .
6 Conclusion
In this paper, I have calculated the normalization of the 8-parameter Fisher-Bingham distribution on using its series expansion. This is given in Eq. 16. By construction, the normalization for the distribution [6] is a simplification and given in Eq. 17. Further, a piecewise approximation of was derived in closed form and seems to perform well for large or , away from the region where .
An algorithm for computing Eq. 16 numerically was described in Section 4. As the sequence of is not, in general, non-negative and only decreases in absolute value to zero after a certain point, a truncation tolerance based on successive partial sums is not sufficient to robustly calculate the normalization. Instead, the proposed technique groups contiguous into and as defined in Eq. 25. The stopping condition is then described using partial sums of . The series calculation of the normalization is shown to be robust and matches that obtained from numerical integration. The summation can be computationally much faster than numerical integration, although this depends on their respective tolerance settings.
With in hand, exact maximum likelihood fits can be performed using Eq. 6. As an example, a synthetic dataset was generated along a small circle on the sphere. Maximum likelihood fits performed using ill-described the data, while fits using exhibited better agreement. As is a superset of , it should allow for more flexible descriptions of directional data. A Python package, extended from the implementation of in [2], is available and contains all the distributions described in this paper111https://github.com/tianluyuan/sphere.
Acknowledgements
This work stemmed out of discussions with Dmitry Chirkin and attempts to fit directional data from [14]. I would also like to thank Kareem Farrag and Austin Schneider for discussions on other potential approaches to calculate the normalization. The author is supported in part by NSF grant PHY-1607644 and by the University of Wisconsin Research Committee with funds granted by the Wisconsin Alumni Research Foundation.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] K. V. Mardia and P. E. Jupp. Directional statistics . Wiley Series in Probability and Statistics. John Wiley & Sons, Ltd., Chichester, 2000.
- 2[2] Eric Daniël Fraenkel. From Radio Pulse to Elusive Particle . Ph D thesis, University of Groningen, 2014.
- 3[3] Ronald Fisher. Dispersion on a sphere. Proc. Roy. Soc. London. Ser. A. , 217:295–305, 1953.
- 4[4] John T. Kent. The Fisher-Bingham distribution on the sphere. J. Roy. Statist. Soc. Ser. B , 44(1):71–80, 1982.
- 5[5] C. Bingham and K. V. Mardia. A small circle distribution on the sphere. Biometrika , 65(2):379–389, 1978.
- 6[6] Louis-Paul Rivest. On the information matrix for symmetric distributions on the hypersphere. Ann. Statist. , 12(3):1085–1089, 09 1984.
- 7[7] Byungwon Kim, Stephan Huckemann, Jörn Schulz, and Sungkyu Jung. Small-sphere distributions for directional data with application to medical imaging. Scandinavian Journal of Statistics , 0(0), 2019.
- 8[8] K. V. Mardia. Statistics of directional data . Academic Press, London-New York, 1972. Probability and Mathematical Statistics, No. 13.
