The 8-parameter Fisher-Bingham distribution on the sphere

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 ($\mathrm{FB}_8$) is an eight-parameter family of probability density functions (PDF) on $S^2$ 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 $\mathrm{FB}_8$, such as the Kent ($\mathrm{FB}_5$) 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 $\mathrm{FB}_8$ 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 $\mathrm{FB}_8$ can be expressed as an infinite sum consisting of hypergeometric functions, similar to that of the $\mathrm{FB}_5$. This allows the normalization to be computed under summation with adequate stopping conditions. I then fit the $\mathrm{FB}_8$ to a synthetic dataset using a maximum-likelihood approach and show its improvements over a fit with the more restrictive $\mathrm{FB}_5$ distribution.

Equations

$$f_5(\vec{x})=c_5(\kappa ,\beta {)}^{-1}\exp \left\{\kappa {\vec{\gamma }}_1\cdot \vec{x}+\beta [({\vec{\gamma }}_2\cdot \vec{x}{)}^2-({\vec{\gamma }}_3\cdot \vec{x}{)}^2]\right\},$$

$$f_4(\vec{x})$$

$$=c_4(\kappa ,\beta {)}^{-1}\exp \left\{\kappa {\vec{\gamma }}_1\cdot \vec{x}+\beta [1-({\vec{\gamma }}_1\cdot \vec{x}{)}^2]\right\}.$$

$$f_6(\vec{x})$$

$$f_8(\vec{x})$$

$$=c_8(\kappa ,\beta ,\eta ,\vec{\nu }{)}^{-1}\exp \left\{\kappa \vec{\nu }\cdot {\mathbf{\Gamma }}^T\vec{x}+\beta [({\vec{\gamma }}_2\cdot \vec{x}{)}^2-\eta ({\vec{\gamma }}_3\cdot \vec{x}{)}^2]\right\},$$

$$f_8({\vec{x}}^{\ast })=c_8(\kappa ,\beta ,\eta ,\vec{\nu }{)}^{-1}\exp \left\{\kappa \vec{\nu }\cdot {\vec{x}}^{\ast }+\beta (x_2^{\ast 2}-\eta x_3^{\ast 2})\right\}.$$

$${\vec{x}}^{\ast }=(\cos \theta ,\sin \theta \cos \varphi ,\sin \theta \sin \varphi )$$

$$f_8({\vec{x}}^{\ast })=\frac{\exp \left\{\kappa ({\nu }_1\cos \theta +{\nu }_2\sin \theta \cos \varphi +{\nu }_3\sin \theta \sin \varphi )+\beta {\sin }^2\theta ({\cos }^2\varphi -\eta {\sin }^2\varphi )\right\}}{c_8(\kappa ,\beta ,\eta ,\vec{\nu })},$$

$$c_8(\kappa ,\beta ,\eta ,\vec{\nu })$$

$$\equiv {\int }_0^{\pi }{\int }_0^{2\pi }\mathcal{I}\mathrm{d}\varphi \mathrm{d}\theta .$$

$$\mathcal{I}$$

$$\displaystyle=e^{\kappa\nu_{1}\cos\theta}\sum_{l,k,j=0}^{\infty}\sum_{i=0}^{j}\bigg{\{}\frac{\kappa^{l+k}\beta^{j}\nu_{2}^{l}\nu_{3}^{k}(-\eta)^{i}}{l!k!j!}\binom{j}{i}\sin^{2j+l+k+1}\theta\cos^{2(j-i)+l}\phi\sin^{2i+k}\phi\bigg{\}}.$$

$$c_8(\kappa ,\beta ,\eta ,\vec{\nu })$$

$$\displaystyle\quad\times\sum_{i=0}^{j}(-\eta)^{i}\binom{j}{i}B\left(j-i+l+\frac{1}{2},i+k+\frac{1}{2}\right)\bigg{\}}\mathrm{d}\theta$$

$$\displaystyle=2\sqrt{\pi}\sum_{l,k,j=0}^{\infty}\bigg{\{}\frac{\kappa^{2(l+k)}\beta^{j}\nu_{2}^{2l}\nu_{3}^{2k}}{(2l)!(2k)!j!}\left|\frac{\kappa\nu_{1}}{2}\right|^{-j-l-k-\frac{1}{2}}I_{j+l+k+\frac{1}{2}}(|\kappa\nu_{1}|)$$

$$\displaystyle\quad\times\sum_{i=0}^{j}(-\eta)^{i}\binom{j}{i}\Gamma\left(j-i+l+\frac{1}{2}\right)\Gamma\left(i+k+\frac{1}{2}\right)\bigg{\}}$$

$$\displaystyle=2\sqrt{\pi}\sum_{l,k,j=0}^{\infty}\bigg{\{}\frac{\kappa^{2(l+k)}\beta^{j}\nu_{2}^{2l}\nu_{3}^{2k}}{(2l)!(2k)!j!}\frac{\Gamma\left(k+1/2\right)\Gamma\left(j+l+1/2\right)}{\Gamma\left(j+l+k+3/2\right)}$$

$$\displaystyle\quad\times{}_{0}F_{1}\left(;j+l+k+\frac{3}{2};\frac{\kappa^{2}\nu_{1}^{2}}{4}\right){}_{2}F_{1}\left(-j,k+\frac{1}{2};\frac{1}{2}-j-l;-\eta\right)\bigg{\}}$$

$$c_{6}(\kappa,\beta,\eta)=2\pi\sum_{j=0}^{\infty}\bigg{\{}\frac{\beta^{j}\Gamma\left(j+1/2\right)}{j!\Gamma\left(j+3/2\right)}{}_{0}F_{1}\left(;j+\frac{3}{2};\frac{\kappa^{2}}{4}\right){}_{2}F_{1}\left(-j,\frac{1}{2};\frac{1}{2}-j;-\eta\right)\bigg{\}},$$

$$c_6(\kappa ,\beta ,\eta )$$

$$=e^{\beta \left(1+\frac{{\kappa }^2}{4{\beta }^2}\right)}{\int }_0^{\pi }{\int }_0^{2\pi }{e}^{-\beta {\left(\cos \theta -\frac{\kappa }{2\beta }\right)}^2}{e}^{-\beta (1+\eta ){\sin }^2\theta {\sin }^2\varphi }\mathrm{d}\varphi \sin \theta \mathrm{d}\theta$$

$$=e^{\beta \left(1+\frac{{\kappa }^2}{4{\beta }^2}\right)}{\int }_{-1}^1{\int }_0^{2\pi }{e}^{-\beta {\left(z-\frac{\kappa }{2\beta }\right)}^2}{e}^{-\beta (1+\eta )(1-z^2){\sin }^2\varphi }\mathrm{d}\varphi \mathrm{d}z,$$

$$c_6(\kappa ,\beta ,\eta )$$

$$\approx 2\pi e^{\beta \left(1+\frac{{\kappa }^2}{4{\beta }^2}\right)}{}_1{F}_1\left(\frac{1}{2};1;\beta (1+\eta )\left(\frac{{\kappa }^2}{4{\beta }^2}-1\right)\right){\int }_{-1}^1{e}^{-\beta {\left(z-\frac{\kappa }{2\beta }\right)}^2}\mathrm{d}z$$

$$\approx 2\pi e^{\beta \left(1+\frac{{\kappa }^2}{4{\beta }^2}\right)}{}_1{F}_1\left(\frac{1}{2};1;\beta (1+\eta )\left(\frac{{\kappa }^2}{4{\beta }^2}-1\right)\right)\sqrt{\frac{\pi }{\beta }}(\kappa <2\beta ).$$

$$c_6(\kappa ,\beta ,\eta )\approx 2\pi e^{\kappa }[(\kappa -2\beta )(\kappa +2\beta \eta ){]}^{-\frac{1}{2}}(\kappa >2\beta ),$$

$$A_{L,K,J}$$

$$B_{L,K,J}$$

Full text

K****eywords Directional statistics  $\cdot$ Fisher-Bingham distribution  $\cdot$ Kent distribution  $\cdot$ Von Mises-Fisher distribution

1 Introduction

Directional statistics involves the study of probability density functions (PDF) with supports on $S^{n}\subset\mathbb{R}^{n+1}$. This paper will focus on distributions on the sphere, $S^{2}$. 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 $S^{2}$ 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 ($\mathrm{FB}_{5}$) distribution [4],

[TABLE]

where $c_{5}(\kappa,\beta)$ is the normalization constant, $\vec{x}$ is a unit vector on $S^{2}$, $\kappa$ and $\beta$ are non-negative parameters, and $\vec{\gamma}_{i}$ are unit vectors that correspond to the columns of a $3\times 3$ orthogonal matrix, $\mathbf{\Gamma}$, which determines the orientation of the PDF. An additional constraint, $\kappa>2\beta$, is required to interpret the $\mathrm{FB}_{5}$ 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 $\beta=0$. Both the vMF and $\mathrm{FB}_{5}$ distributions have been well-studied, and have found use in many applications. However, $\mathrm{FB}_{5}$ suffers from the restriction that its PDFs must be symmetric across two great circles intersecting at $90^{\circ}$ on the sphere. Data that clusters along small circles, for example the non-equatorial lines of latitude, are ill-described by $\mathrm{FB}_{5}$.

An alternative to $\mathrm{FB}_{5}$ is the small-circle distribution proposed in [5]. This is a four-parameter subclass ($\mathrm{FB}_{4}$) of the Fisher-Bingham family and can be written as,

[TABLE]

With $\kappa<2\beta$, $f_{4}$ describes small-circle distributions on the sphere as shown in the left panel of Fig. 1. Removing the constraint on $\kappa$ and $\beta$, the only difference between $\mathrm{FB}_{5}$ and $\mathrm{FB}_{4}$ is the sign of the term in square-brackets. The $\mathrm{FB}_{4}$ distribution is completely specified by four parameters, since $\vec{\gamma}_{1}$ 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 $\mathrm{FB}_{4}$ and $\mathrm{FB}_{5}$.

In order to perform a maximum likelihood fit of directional data, the PDFs need to be normalized. It was shown in [5] that $c_{4}(\kappa,\beta)$ can be written in terms of the confluent hypergeometric function. It was shown in [4] that $c_{5}(\kappa,\beta)$ can be written as an infinite sum consisting of modified Bessel functions of the first kind. These approaches motivated the calculation of the $\mathrm{FB}_{8}$ normalization discussed in Section 3, but first a natural generalization of $\mathrm{FB}_{4}$ and $\mathrm{FB}_{5}$ is given in Section 2.

2 The Fisher-Bingham distribution

It is simple to construct a 6-parameter PDF ($\mathrm{FB}_{6}$) that generalizes $\mathrm{FB}_{4}$ and $\mathrm{FB}_{5}$ [6]. This is given here as

[TABLE]

with $\kappa\geq 0$, $\beta\geq 0$ and $|\eta|\leq 1$. Clearly, $\eta=1$ corresponds to Eq. 1 and $\eta=-1$ to Eq. 2. For $\kappa>2\beta$, $f_{6}$ has a single maximum, corresponding to where $\vec{x}$ is aligned with $\vec{\gamma}_{1}$. For $\kappa<2\beta$, $f_{6}$ can describe either the small-circle distribution of [5] or a bimodal distribution where the modes are $180^{\circ}$ degrees apart on a small circle as shown in the left panel of Fig. 2. As is the case for $\mathrm{FB}_{5}$, the $\mathrm{FB}_{6}$ 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 $\mathrm{FB}_{4}$ distributions. This can be generalized further by combining the vMF and $\mathrm{FB}_{6}$ distributions, which results in the Fisher-Bingham distribution ($\mathrm{FB}_{8}$). It is parametrized here as,

[TABLE]

where $\vec{\nu}$ 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 $\mathbf{\Gamma}\vec{\nu}-\vec{\gamma}_{1}$, and the second term is proportional to Eq. 4. Thus, $\mathrm{FB}_{8}$ does not necessarily have to be symmetric about a great circle. When $\vec{\nu}=(1,0,0)$ Eq. 6 reduces to Eq. 4.

3 Calculating the $\mathbf{\mathrm{\textbf{FB}}_{6}}$ and $\mathbf{\mathrm{\textbf{FB}}_{8}}$ normalizations

3.1 Exact series solution

The $\mathrm{FB}_{8}$ distribution was proposed in [8, 9], though it has not been widely applied due to difficulties in computing $c_{8}(\kappa,\beta,\eta,\vec{\nu})$ [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 $c_{8}(\kappa,\beta,\eta,\vec{\nu})$ is desirable to perform maximum likelihood inference using Eq. 6. This is the subject of this Section and Section 4.

Since $\mathbf{\Gamma}$ simply enacts a rotation of the sphere, the normalization constants above are independent of $\mathbf{\Gamma}$ and it is simpler to work in the standard frame, where the coordinate axes are defined by the columns of $\mathbf{\Gamma}$ with the $z$-axis corresponding to $\vec{\gamma}_{1}$ [11]. The coordinate transformation $\vec{x}^{*}=\mathbf{\Gamma}^{T}\vec{x}$ allows us to write Eq. 6 as,

[TABLE]

In spherical coordinates

[TABLE]

and

[TABLE]

where

[TABLE]

Taylor expanding $\mathcal{I}$ 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 $\phi$ vanishes unless $k$ and $l$ are both even,

[TABLE]

where $\Gamma$ denotes the gamma function, $B$ the beta function, $I_{v}$ the modified Bessel function of the first kind, ${}_{0}F_{1}$ the confluent hypergeometric limit function, and ${}_{2}F_{1}$ 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 $c_{8}(\kappa,\beta,\eta,\vec{\nu})$ can be computed to good approximation with adequate stopping conditions on $j$, $k$, and $l$. By setting $\vec{\nu}=(1,0,0)$, the only nonzero terms occur when $k=l=0$ and Eq. 16 simplifies to

[TABLE]

the normalization for $\mathrm{FB}_{6}$. Further setting $\eta=1$ and applying Eq. (15.1.21) in [12] recovers $c_{5}(\kappa,\beta)$ as computed in [4].

3.2 Closed-form approximation for $\mathbf{c_{6}(\kappa,\beta,\eta)}$

If $\kappa$ or $\beta$ is large, $c_{6}(\kappa,\beta,\eta)$ can be approximated piecewise in two separate regimes: $\kappa<2\beta$ and $\kappa>2\beta$. For $\kappa<2\beta$, note that we can rewrite

[TABLE]

where $z=\cos\theta$. In order to factorize the integration over $\phi$ and $z$, we make the assumption that the term modulated by $\sin^{2}\phi$ can be fixed to $z=\kappa/(2\beta)$, where the maximum of $f_{6}$ occurs in latitude. Then, using Eq. (13.1.27) in [12],

[TABLE]

The last line uses an approximation for large $\beta$ [5].

In the case of $\kappa>2\beta$, $f_{6}$ is maximal at $\theta=0$ and for large $\kappa$ becomes approximately a bivariate normal distribution. By setting $\phi=0$ and $\phi=\pi/2$ and Taylor expanding in $\theta$, 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 $\kappa=2\beta$. The saddlepoint approximation of [11] is also shown (dotted orange).

These approximations are useful when working with large $\kappa$ or $\beta$, 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 $\ln c_{6}(\kappa,\beta,\eta)$, which can be approximated with Eqs. 23 and 24 without running into computational overflows. For $c_{8}(\kappa,\beta,\eta,\vec{\nu})$, 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 $c_{6}(\kappa,\beta,\eta)$ but not for $c_{8}(\kappa,\beta,\eta,\vec{\nu})$.

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 $a_{l,k,j}$ be the summand of Eq. 16 such that $c_{8}(\kappa,\beta,\eta,\vec{\nu})=\sum_{l,k,j=0}^{\infty}a_{l,k,j}$. For $\eta\leq 0$, $a_{l,k,j}$ is non-negative for all $l,k,j$ while for $\eta>0$, $a_{l,k,j}$ is guaranteed to be non-negative only for even $j$. As such, $a_{l,k,j}$ may be an alternating sequence in $j$. Furthermore, $|a_{l,k,j}|$ 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 $c_{8}(\kappa,\beta,\eta,\vec{\nu})$ would be to reindex and order $b_{n}=a_{l,k,j}$ such that $|b_{n+1}|\geq|b_{n}|$ for all $n\in\mathbb{W}$, and perform the summation starting from $n={0,1,2\ldots}$ up to certain tolerance. However, this ordering is difficult to evaluate. In practice, a robust calculation can be obtained by setting a step size $s\in 2\mathbb{N}$ and defining

[TABLE]

This ensures that an even number of $a_{l,k,j}$ terms are summed at each step and $c_{8}(\kappa,\beta,\eta,\vec{\nu})=\sum_{L,K,J=0}^{\infty}A_{L,K,J}$. Then for some tolerance $\epsilon>0$, the stopping algorithm is

⬇

$L$, $P_{L}$, $\tilde{c}_{8}$ = 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 $\tilde{c}_{8}$ denotes the series calculation of $c_{8}(\kappa,\beta,\eta,\vec{\nu})$. This nested summation routine loops through $J$, $K$, and $L$ in that order and breaks once the contribution to the partial sum of $B_{L,K,J}$ for the current index is within the tolerance and is less than the previous term. In tests, $\epsilon=10^{-12}$ and $s=13$ seemed to perform well. A comparison of $\ln c_{8}$ 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 $J$ while setting $l=k=0$. This amounts to just running the innermost loop. One final thing to note is that $a_{l,k,j}$ 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 $\exp(\ln a_{l,k,j})$, explicitly taking the logarithm before exponentiation.

5 Example application

To illustrate the performance of the $\mathrm{FB}_{8}$ distribution in modeling directional data, a synthetic dataset was randomly sampled from an $\mathrm{FB}_{8}$ distribution that peaked along a small circle on the sphere. An unbinned maximum likelihood fit was then performed using the $\mathrm{FB}_{5}$ and $\mathrm{FB}_{8}$ distributions [4, 1]. The SLSQP routine was used to perform a constrained fit of $\mathrm{FB}_{5}$ and the L-BFGS-B routine was used for a bounded fit of $\mathrm{FB}_{8}$ [13].

The results are shown in Fig. 4. The synthetic data (black points) is poorly described by the $\mathrm{FB}_{5}$ distribution, but well described by $\mathrm{FB}_{8}$. The best-fit negative log-likelihoods are $31.7$ for the $\mathrm{FB}_{5}$ and $-12.4$ for the $\mathrm{FB}_{8}$, indicating much better data agreement using $\mathrm{FB}_{8}$.

6 Conclusion

In this paper, I have calculated the normalization of the 8-parameter Fisher-Bingham distribution on $S^{2}$ using its series expansion. This is given in Eq. 16. By construction, the normalization for the $\mathrm{FB}_{6}$ distribution [6] is a simplification and given in Eq. 17. Further, a piecewise approximation of $c_{6}(\kappa,\beta,\eta)$ was derived in closed form and seems to perform well for large $\kappa$ or $\beta$, away from the region where $\kappa=2\beta$.

An algorithm for computing Eq. 16 numerically was described in Section 4. As the sequence of $a_{l,k,j}$ 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 $a_{l,k,j}$ into $A_{L,K,J}$ and $B_{L,K,J}$ as defined in Eq. 25. The stopping condition is then described using partial sums of $B_{L,K,J}$. 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 $c_{8}(\kappa,\beta,\eta,\vec{\nu})$ 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 $\mathrm{FB}_{5}$ ill-described the data, while fits using $\mathrm{FB}_{8}$ exhibited better agreement. As $\mathrm{FB}_{8}$ is a superset of $\mathrm{FB}_{5}$, it should allow for more flexible descriptions of directional data. A Python package, extended from the implementation of $\mathrm{FB}_{5}$ 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.