Bayesian Model Choice for Directional Data

Abstract

This article is concerned with the problem of choosing between competing models for directional data. In particular, we consider the question of whether or not two independent samples of axial data come from the same Bingham distribution. This is not a straightforward question to answer, due to the intractable nature of the parameter-dependent normalizing constant of the Bingham distribution. We propose three different methods to perform this task within a Bayesian framework, and apply the methodology to a real dataset on earthquakes in New Zealand. R code to run our methods is available in online supplementary materials.

KEYWORDS:

• Bingham distribution
• Doubly intractable distributions
• Model choice
• Reversible jump MCMC

1 Introduction

The Bingham distribution, introduced by Bingham (1974), is one of the most important and useful models for data which arise as unsigned directions (axes). Bingham (1974) derived its main properties, and since then it has been studied further; see for example Kent (1987), Mardia and Jupp (2000), Kume and Walker (2014), and Dryden (2005), who studied the properties of the Bingham distribution in very high dimensions.

The Bingham distribution can be constructed by conditioning a zero-mean multivariate Normal (MVN) distribution in ${\mathbb{R}}^q$ to lie on the unit sphere ${\mathcal{S}}^{q-1}$. For $\mathit{x}\in {\mathcal{S}}^{q-1}$, the density of the Bingham distribution with respect to the uniform measure on ${\mathcal{S}}^{q-1}$ is given by (1) $$f(\mathit{x};A)=\frac{\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}(-{\mathit{x}}^{T}A\mathit{x})}{c(A)},\hspace{1em}{\mathit{x}}^T\mathit{x}=1\hspace{1em}\text{and}\hspace{1em}\mathit{x}\in {\mathbb{R}}^q,$$(1) where the parameter A is a real symmetric q × q matrix and c(A) is the corresponding normalizing constant.

Let $A=V\Lambda V^T$ be the Singular Value Decomposition of A. Then maximum likelihood estimates of V (whose columns are the eigenvectors of A) are easily obtained. However, inferring the eigenvalues of A (the elements of the diagonal matrix Λ) is a rather challenging task, since evaluation of the parameter-dependent normalizing constant c(A) is required, which, in the general case, is not available in closed form. This poses significant challenges when fitting the Bingham distribution to some data in either a frequentist or Bayesian setting. There have been several approaches proposed in the literature to approximate/estimate the normalizing constant, and hence facilitate maximum likelihood estimation, based on asymptotic expansions (Kent 1987), saddlepoint approximations (Kume and Wood 2005) and more recently, holonomic gradient methods (Sei and Kume 2015). However, these are nontrivial numerical methods to implement, which in any case involve some level of approximation of c(A). Likewise, standard approaches to model selection involving Bingham distributions would require evaluation of c(A).

In the Bayesian setting, MCMC is also complicated by c(A), since standard Metropolis-Hastings updates would require evaluation of this constant at current and proposed configurations. This is known as a doubly intractable problem, for which various algorithms have been proposed (e.g., Møller et al. 2006; Murray, Ghahramani, and MacKay 2006; Walker 2011). The specific case of the Bingham distribution was considered by Walker (2014) and Fallaize and Kypraios (2016).

The above discussion relates to inference given data from a single Bingham distribution. In this article, we consider, within a Bayesian framework, the problem of model selection in the context of Bingham distributions. That is, suppose we have two independent datasets consisting of some axial data, considered to be realizations from Bingham distributions. A natural question that often arises is the following: do these two datasets come from the same Bingham distribution? A common approach is to compute model probabilities/Bayes factors for the two competing scenarios. Again, standard techniques are complicated by the intractability of the constant c(A). Recall that, generally speaking, there are two approaches for calculating Bayes factors/model probabilities using MCMC:

• within-model estimation. This involves estimating the model evidence (marginal likelihood) for each model separately, for example, using MCMC samples and Chib’s method. In the doubly intractable case, although samples can be obtained from each posterior distribution using one of the aforementioned algorithms, one would also have to plug in an approximation to the unknown normalizing constant (c(A) for the Bingham distribution). So in practice, the model evidence has to be approximated.

• across-model estimation. Reversible jump can be used to run one sampler which jumps between competing models and thus allows direct (Monte Carlo) estimation of the model probabilities. In the doubly intractable case, the Reversible Jump Exchange algorithm (Caimo and Friel 2013) is an across-model approach that enables exact (in a Monte Carlo sense) calculation of posterior model probabilities without the need to approximate the intractable constant, but this is more computationally-intensive and requires careful design for constructing an efficient sampler.

Although these general approaches are available, implementations and performance are context-specific. The main purpose of this article is to show how methods in each of the above categories can be successfully applied to perform model selection for Bingham distributions, and to investigate performance in terms of ability to discriminate between competing models depending on sample size and parameter configurations. Specifically, we investigate three different approaches: the Reversible Jump Exchange algorithm of Caimo and Friel (2013); a within-model approach in which the model probabilities are approximated; an alternative to the calculation of Bayes Factors, based on calculating (multivariate) highest posterior density regions for a certain parameter of interest, one value of which corresponds to distributional equality. The latter two approaches are computationally fast and easy to implement. They only require the practitioner being able to fit a Bingham distribution to directional data within a Bayesian framework, something that can efficiently be done using the algorithm of Fallaize and Kypraios (2016). They also provide a useful cross-check for the Reversible Jump Exchange algorithm, which is appealing in the sense that it provides exact (Monte Carlo) calculation of posterior model probabilities, but requires more careful sampler design and is more computationally intensive. Hence, all three approaches are useful in practice.

The article is structured as follows. In Section 2 we introduce a motivating example of data from earthquakes in New Zealand. In Section 3, we review how Bayesian inference can proceed for the Bingham distribution, and in Section 4 we describe our three strategies for performing model choice for the Bingham distribution in the Bayesian framework. In Section 5, we illustrate the methods on simulated data, and perform a simulation study to investigate the extent to which it is possible to discriminate between models as parameters and sample sizes vary. In Section 6, we apply our methods to the New Zealand earthquake data, before concluding with a discussion.

2 Motivating Application

Here, we introduce a motivating example, which relates to earthquake data in New Zealand. An earthquake gives rise to three orthogonal axes, and geophysicists are interested in analyzing such data in order to compare earthquakes at different locations and/or at different times. An earthquake gives rise to a pair of orthogonal axes, known as the compressional (P) and tensional (T) axes, from which a third axis, known as the null (B) axis is obtained via $B=P\times T$. Each of these quantities are determined only up to sign, and so models for axial data are appropriate. Arnold and Jupp (2013) treated these data as orthogonal axial 3-frames in ${\mathbb{R}}^3$, but, following Fallaize and Kypraios (2016) we will illustrate our methods using the B axes. In general, an orthogonal axial r-frame in ${\mathbb{R}}^q,r\le q$, is an ordered set of r axes, $\{\pm u_1,\dots ,\pm u_r\}$, where $u_1,\dots ,u_r$ are orthonormal vectors in ${\mathbb{R}}^q$ (Arnold and Jupp 2013). The familiar case of data on the sphere ${\mathcal{S}}^2$ is the special case corresponding to $q=3,r=1$, which is the case we consider here.

The data consist of three clusters of observations derived from earthquakes in New Zealand, where each observation represents a B axis obtained from a particular earthquake. The first two clusters each consist of 50 observations from earthquakes near Christchurch which took place before and after a large earthquake on the 22nd of February 2011. For these two clusters, the P and T axes are quite highly concentrated in the horizontal plane, and as a result the majority of the B axes are concentrated about the vertical axis (see Figure 1). From an applied viewpoint, it is of interest to geophysicists to establish whether there is a difference between the pattern of earthquakes before and after the large earthquake. The third cluster is a more diverse set of 32 observations obtained from earthquakes in the north of New Zealand’s South Island. Visually, this cluster appears to exhibit more dispersion than the other two samples. Hence, one approach is to model the B axes using the Bingham distribution, and the principal task is to assess the evidence that the different samples are from populations with the same Bingham distribution or not.

Fig. 1 The B axes for the earthquake data. Top left: Christchurch before the large earthquake. Top right: Christchurch after the large earthquake. Bottom: South Island. The data have been rotated to principal axes, such that the principal axis is the North–South axis. Relatively, the two Christchurch clusters are more concentrated about the vertical axis, and the South Island cluster more dispersed. All points have been plotted in the same hemisphere to enable easier comparison of the variability in each sample. In some cases, the antipodal point was the observation—as these are axial data, both are treated as the same observation and it does not matter which is observed.

3 Bayesian Inference

Consider the probability density function of the Bingham distribution as given in (1). If $A=V\Lambda V^T$ is the Singular Value Decomposition of A where V is orthogonal and $\Lambda =\text{diag}({\lambda }_1,\dots ,{\lambda }_q)$, then it can be shown that if $\mathit{x}$ is drawn from a distribution with probability density function $f(·;A)$, the random vector $V^T\mathit{x}$ is drawn from a distribution with density $f(·;\Lambda )$ (see, e.g., Kume and Walker 2006; Kume and Wood 2007). Therefore, for ease of exposition and without loss of generality, we assume that $A=\Lambda =\text{diag}({\lambda }_1,\dots ,{\lambda }_q)$. Moreover, to ensure identifiability, we assume that ${\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_q=0$ (Kent 1987) and denote $\mathit{\lambda }:=({\lambda }_1,\dots ,{\lambda }_{q-1})$, the ordered vector of parameters to be inferred. Hence, the probability density function becomes (2) $$f(\mathit{x};\mathit{\lambda })=\frac{\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{-\sum _{i=1}^{q-1}{\lambda }_i{x}_i^2\right\}}{c(\mathit{\lambda })}$$(2)

with respect to the uniform measure on the sphere and $$c(\mathit{\lambda })={\int }_{\mathit{x}\in S^{q-1}}\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{-\sum _{i=1}^{q-1}{\lambda }_i{x}_i^2\right\}\text{d}{S}^{q-1}(\mathit{x}).$$

(Throughout, we write $\mathit{\lambda }$ in place of A or Λ, since it is to be understood that A is of the form $A=\Lambda =\text{diag}(\mathit{\lambda },0)$.) This integral can be expressed as a hypergeometric function of matrix argument (see, e.g., Mardia and Jupp 2000), which is an infinite series involving zonal polynomials and hence not a tractable function to work with in practice. We note that in the case when there are only two distinct eigenvalues of A, the distribution reduces to a Watson distribution, in which case the normalizing constant can be easily evaluated. The Watson distribution is appropriate when there is rotational symmetry about some axis, which is not the case for the Bingham distribution in general. Here, we assume that this possibility has been discounted, and that the general Bingham distribution is appropriate.

Let ${\mathit{x}}_j^T=(x_j^1,x_j^2,\dots ,x_j^q)\in {\mathcal{S}}^{q-1}$ be a unit vector, and suppose $\tilde{\mathit{x}}=({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n)$ is a random sample of n unit vectors from the Bingham distribution with density (2). Then the likelihood function is given by (3) $$\begin{array}{c}\pi (\tilde{\mathit{x}}|\mathit{\lambda })=\frac{1}{c{(\mathit{\lambda })}^n}\mathrm{\text{exp}}\hspace{0.17em}\left\{-\sum _{i=1}^{q-1}{\lambda }_i\sum _{j=1}^n{\left(x_j^i\right)}^2\right\}\\ =\frac{1}{c{(\mathit{\lambda })}^n}\mathrm{\text{exp}}\hspace{0.17em}\left\{-n\sum _{i=1}^{q-1}{\lambda }_i{\tau }_i\right\}\\ =\frac{1}{c{(\mathit{\lambda })}^n}{\pi }_u(\tilde{\mathit{x}}|\mathit{\lambda })\end{array}$$(3) where ${\tau }_i=\frac{1}{n}{\sum }_{j=1}^n{(x_j^i)}^2$ and ${\pi }_u(\tilde{\mathit{x}}|\mathit{\lambda })=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{-n{\sum }_{i=1}^{q-1}{\lambda }_i{\tau }_i\right\}$. The data can therefore be summarized by $(n,{\tau }_1,\dots ,{\tau }_{q-1})$, and $({\tau }_1,\dots ,{\tau }_{q-1})$ are sufficient statistics for $({\lambda }_1,\dots ,{\lambda }_{q-1})$.

The likelihood contains the normalizing constant $c(\mathit{\lambda })$ and the fact that there does not exist a closed form expression for it makes Bayesian inference for $\mathit{\lambda }$ challenging. Denote by $\pi (\mathit{\lambda })$a prior distribution that is placed on ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{q-1}$, with support ${\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_{q-1}\ge 0$. Hence, the posterior distribution of $\mathit{\lambda }$ is given by (4) $$\pi (\mathit{\lambda }|\tilde{\mathit{x}})\propto \frac{{\pi }_u(\tilde{\mathit{x}}|\mathit{\lambda })}{c{(\mathit{\lambda })}^n}\pi (\mathit{\lambda }),{\lambda }_1\ge {\lambda }_2\ge \cdots {\lambda }_{q-1}\ge 0,$$(4) which also involves the intractable normalizing constant $c(\mathit{\lambda })$. As explained in Fallaize and Kypraios (2016), in principle one can employ one of the proposed methods in the literature to approximate the normalizing constant (see, e.g., Kent 1987; Kume and Wood 2005; Sei and Kume 2015) and then use, for example, a Metropolis-Hastings algorithm to sample from $\pi (\mathit{\lambda }|\tilde{\mathit{x}})$. However, despite such an approach being feasible, it is not entirely satisfactory. First, the approximation would have to be computed at each iteration of the Markov chain and second, and most importantly, the stationary distribution of such an MCMC algorithm will not be the distribution of interest, but rather an approximation to it (despite how good this might be). Fallaize and Kypraios (2016) showed how one can bypass the need of calculating the normalizing constant by making use of recent developments in computational methods for doubly intractable distributions.

Of course, in general we will not have V = I. In practice, we do the following. Let X be the n × q data matrix, that is, row j of X is ${\mathit{x}}_j^T$, and $\overline{X}=\frac{1}{n}(X^{T}X)$. Then let $Q\Phi Q^T$ be the Singular Value Decomposition of $\overline{X}$, whereby the diagonal elements of $\Phi$ are the τ_i above. These are unchanged under any rotation of X. Hence, in practice, we take our data ${\tau }_i,i=1,\dots ,q$, to be the eigenvalues of $\overline{X}$. Due to the identifiability constraints on $\mathit{\lambda }$, we work with the ordered eigenvalues (with τ₁ being the smallest and so on), and τ_q is redundant since the constraint ${\lambda }_q=0$ is imposed. This procedure is equivalent to performing inference for $\mathit{\lambda }$ as above, with V = I.

This can be viewed as (and is equivalent to) the following hybrid algorithm for performing inference for V and $\mathit{\lambda }$. Q is in fact the maximum likelihood estimate (mle) of V. Given the data X in the original coordinates, compute Q and set Y = XQ. Then the eigenvalues of $\frac{1}{n}(Y^{T}Y)$ are the same as those of $\overline{X}$, which will also be equal to $\frac{1}{n}{\sum }_{j=1}^n{(y_j^i)}^2$. So working with the eigenvalues of $\overline{X}$ computed using the original data is equivalent to taking the mle of V, using this to rotate the data, and performing inference for $\mathit{\lambda }$ using the posterior distribution (4) with the rotated data.

4 Bayesian Model Choice

4.1 Preliminaries

Given a set of K competing models indexed by a model indicator $k\in \{1,\dots ,K\}$, the principal interest in Bayesian model selection lies in quantifying the extent to which the observed data support each model k. The Bayesian treatment for model choice involves joint inference for k and a parameter vector ${\mathit{\lambda }}_k$ where the model indicator determines the dimension of the parameter, and this may vary from model to model.

Suppose data $\tilde{\mathit{x}}$ is assumed to have been generated by model k; the posterior distribution is $$\pi ({\mathit{\lambda }}_k|\tilde{\mathit{x}},k)=\frac{\pi (\tilde{\mathit{x}}|{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_k|k)}{\pi (\tilde{\mathit{x}}|k)},$$where $\pi (\tilde{\mathit{x}}|{\mathit{\lambda }}_k,k)$ represents the likelihood function and $\pi ({\mathit{\lambda }}_k|k)$ the prior distribution of the parameters of model k. The model evidence for model k is the integral $$\pi (\tilde{\mathit{x}}|k)={\int }_{{\mathit{\lambda }}_k}\pi (\tilde{\mathit{x}}|{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_k|k)d{\mathit{\lambda }}_k.$$

It represents the probability of the data $\tilde{\mathit{x}}$ given a particular model k and is vital in Bayesian model choice. The reason is that it allows us to make statements about posterior model probabilities which can be written as $$\pi (k|\tilde{\mathit{x}})=\frac{\pi (\tilde{\mathit{x}}|k)\pi (k)}{\sum _{k=1}^K\pi (\tilde{\mathit{x}}|k)\pi (k)},$$where $\pi (k)$ is the prior probability of model k. (In our examples we have two competing models and take $\pi (k)=0.5$ as a default choice in the absence of any prior information to the contrary.) These posterior model probabilities give rise to the Bayes Factors which give the relative posterior likelihoods of pairs of models, for example $$\frac{\pi (k_1|\tilde{\mathit{x}})}{\pi (k_2|\tilde{\mathit{x}})}=\frac{\pi (k_1)}{\pi (k_2)}\times \frac{\pi (\tilde{\mathit{x}}|k_1)}{\pi (\tilde{\mathit{x}}|k_2)}=\frac{\pi (k_1)}{\pi (k_2)}\times B_{12}$$ where B₁₂ is the Bayes factor of model k₁ relative to k₂, given the observed data $\tilde{\mathit{x}}$, assuming that the prior model probabilities, $\pi (k_1)$ and $\pi (k_2)$ are independent of $\tilde{\mathit{x}}$. The larger the Bayes factor is, the greater the evidence in favor of model k₁ compared to model k₂. Jeffreys (1998) presents a scale in which one can interpret the strength of evidence of one model against another.

4.2 Computing Bayes Factors

In practice, the posterior model probabilities, which are essential to conduct Bayesian model choice, are very rarely analytically available. Therefore, very often, we rely on computational methods. One of the most favored computational methods is Markov chain Monte Carlo (MCMC). The main idea behind MCMC is to construct a Markov chain which has the posterior distribution of interest as its invariant distribution.

In general, there are two main approaches for computing Bayes factors: (i) across-model and (ii) within-model estimation. The within-model approaches aim to approximate/estimate the model evidence $\pi (k|\tilde{\mathit{x}})$, for each model $k=1,\dots ,K$. There are a number of approaches ranging from Laplace approximations to path sampling; see Friel and Wyse (2012) for a comprehensive review. The across-model strategy relies on constructing a Markov chain that crosses the joint model and parameter space and samples from $$\pi ({\mathit{\lambda }}_k,k|\tilde{\mathit{x}})\propto \pi (\tilde{\mathit{x}}|{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_k|k)\pi (k),$$and one of the most popular approaches used in this context is the reversible jump MCMC (RJMCMC) algorithm of Green (1995).

In the context of the Bingham distribution, both within-model strategies and across-model strategies (such as RJMCMC) cannot be used straightforwardly because the normalizing constant $c(\mathit{\lambda })$ in the Bingham likelihood is not available analytically. Here, we present an implementation of the Reversible jump (RJ) exchange algorithm of Caimo and Friel (2013), which is a modification of the reversible jump algorithm of Green (1995) based on the exchange algorithm of Murray, Ghahramani, and MacKay (2006). Essentially, the RJ exchange algorithm enables us to implement a RJ algorithm without having to compute/approximate the unknown normalizing constant $c(\mathit{\lambda })$. We will also present a within-model approach for estimating the model evidence in the context of the Bingham distribution. The motivation for this is to provide a useful method for low-dimensional models to use as an alternative reference to compare with the RJ exchange algorithm.

4.3 Reversible Jump Exchange Algorithm

As described in Section 4.2, the evidence in favor of competing models can be obtained in an across-model fashion. In this case, the parameter space is $({\mathit{\lambda }}_k,k)$, and a Markov chain which samples from the joint posterior distribution $\pi ({\mathit{\lambda }}_k,k|\tilde{\mathit{x}})$ can in principle be constructed using RJMCMC. Given such a chain, inference for the model index k is readily available, since it is directly incorporated as a parameter. Given a current state $({\mathit{\lambda }}_k,k)$, a RJMCMC algorithm essentially consists of two steps: (i) within-model moves which perform updates of the form ${\mathit{\lambda }}_k\to {\mathit{\lambda }}_k^{\prime }$ for the current model k, and (ii) across-model moves which perform updates of the form $({\mathit{\lambda }}_k,k)\to ({\mathit{\lambda }}_{k^{\prime }},k^{\prime })$, whereby a move from model k to model $k^{\prime }$ is proposed, along with a parameter vector ${\mathit{\lambda }}_{k^{\prime }}$ for that model. In our setting, both types of move are complicated by the fact that the normalizing constant from the Bingham likelihood would need to be evaluated when computing the acceptance probabilities for both moves.

For the within-model updates we use the exchange algorithm of Fallaize and Kypraios (2016); see Algorithm 1. In the algorithm, $q(·|\mathit{\lambda })$ is a user-defined proposal distribution and $f(·|\mathit{\lambda })$ is the Bingham distribution (2) with parameter $\mathit{\lambda }$. In our applications, $\mathit{\lambda }=({\lambda }_1,{\lambda }_2)$, and we find setting $q(·|\mathit{\lambda })=N_2(\mathit{\lambda },I)$ works perfectly well.

Algorithm 1

Within-model moveinput: current state $({\mathit{\lambda }}_k,k)$ output: new state $({\mathit{\lambda }}_k^{\prime },k)$

1. Draw ${\mathit{\lambda }}_k^{\prime }\sim q(·|{\mathit{\lambda }}_k)$

2. Draw $\tilde{\mathit{y}}\sim f(·|{\mathit{\lambda }}_k^{\prime })$

3. Accept ${\mathit{\lambda }}_k^{\prime }$ with probability $$\mathrm{\text{min}}\left(1,\frac{{\pi }_u(\tilde{\mathit{x}}|{\mathit{\lambda }}_k^{\prime })\pi ({\mathit{\lambda }}_k^{\prime })q({\mathit{\lambda }}_k|{\mathit{\lambda }}_k^{\prime }){\pi }_u(\tilde{\mathit{y}}|{\mathit{\lambda }}_k)}{{\pi }_u(\tilde{\mathit{x}}|{\mathit{\lambda }}_k)\pi ({\mathit{\lambda }}_k)q({\mathit{\lambda }}_k^{\prime }|{\mathit{\lambda }}_k){\pi }_u(\tilde{\mathit{y}}|{\mathit{\lambda }}_k^{\prime })}\right).$$

We now show how across-model updates can be performed in our setting using an exchange version of RJMCMC proposed by Caimo and Friel (2013), in order to provide model evidence in favour of competing Bingham models (Algorithm 2). We denote model k by an indicator m_k. Within the algorithm, a move from model k to $k^{\prime }$ is proposed from a proposal distribution $h(k^{\prime }|k)$, and given the proposed value $k^{\prime }$, a candidate value for the model-specific parameter ${\mathit{\lambda }}_{k^{\prime }}$ is then drawn from a proposal density $w_{k^{\prime }}({\mathit{\lambda }}_{k^{\prime }}|{\mathit{\varphi }}_{k^{\prime }},m_{k^{\prime }})$. The latter is a model-specific proposal density with parameters ${\mathit{\varphi }}_{k^{\prime }}$ from which ${\mathit{\lambda }}_{k^{\prime }}$ is drawn independently of the current value of ${\mathit{\lambda }}_k$. Hence, this is a form of independence sampler, which obviates the need to tune parameters for the proposal densities for trans-dimensional moves in RJMCMC, which is often a very challenging task.

Algorithm 2

Across-model moveinput: current state $({\mathit{\lambda }}_k,k)$ output: new state $({\mathit{\lambda }}_{k^{\prime }},k^{\prime })$

1. Draw $k^{\prime }\sim h(k^{\prime }|k)$;

2. Given $k^{\prime }$ draw ${\mathit{\lambda }}_{k^{\prime }}\sim w_{k^{\prime }}({\mathit{\lambda }}_{k^{\prime }}|{\mathit{\varphi }}_{k^{\prime }},m_{k^{\prime }})$

3. Given ${\mathit{\lambda }}_{k^{\prime }}$ draw $\tilde{\mathit{y}}\sim f(·;{\mathit{\lambda }}_{k^{\prime }},m_{k^{\prime }})$

4. Accept the move $({\mathit{\lambda }}_k,k)\to {\mathit{\lambda }}_{k^{\prime }},k^{\prime })$ with probability: $\alpha =\mathrm{\text{min}}\{1,r\}$ where r is given (5).

Values for the parameters ${\mathit{\varphi }}_k$ which parameterize the densities $w_k({\mathit{\lambda }}_k|{\mathit{\varphi }}_k,m_k)$ can be obtained by fitting the Bingham distribution to the individual models $k=1,\dots ,K$ using pilot MCMC runs on each model separately (Fallaize and Kypraios 2016), prior to running the main RJMCMC sampler. In our examples, we find that a multivariate normal distribution (MVN) can approximate sufficiently $\pi ({\mathit{\lambda }}_k|\tilde{\mathit{x}})$, and hence ${\mathit{\varphi }}_k$ consists of a mean vector and covariance matrix for the parameters $\{{\lambda }_1^{(k)},\dots ,{\lambda }_{q-1}^{(k)}\}$ corresponding to model k, estimated from the pilot runs. (The conclusion that the MVN approximation for use as a proposal distribution is adequate was determined by monitoring performance of the MCMC, where good mixing and acceptance rates were observed, in addition to visually inspecting the fitted MVN contours with those from a kernel density estimate.) Proposals which do not respect the constraints on the parameters are automatically rejected.

Given a proposal ${\mathit{\lambda }}_{k^{\prime }}$, auxilary data $\tilde{\mathit{y}}$ are drawn from $f(·;{\mathit{\lambda }}_{k^{\prime }},m_{k^{\prime }})$ (Step 3), the sampling mechanism of the data under the model $k^{\prime }$ with parameter ${\mathit{\lambda }}_{k^{\prime }}$. In our case, this will be a form of Bingham distribution, for which exact samples can be drawn using the algorithm of Kent, Ganeiber, and Mardia (2018). The overall acceptance probability is given by $\alpha (({\mathit{\lambda }}_k,k);({\mathit{\lambda }}_{k^{\prime }},k^{\prime }))=\mathrm{\text{min}}\{1,r\},$ where (5) $$\begin{array}{c}r=\frac{f(\tilde{\mathit{x}};{\mathit{\lambda }}_{k^{\prime }},k^{\prime })f(\tilde{\mathit{y}};{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_{k^{\prime }}|k^{\prime })p(k^{\prime })w_k({\mathit{\lambda }}_k|{\mathit{\varphi }}_k,k)h(k|k^{\prime })}{f(\tilde{\mathit{y}};{\mathit{\lambda }}_{k^{\prime }},k^{\prime })f(\tilde{\mathit{x}};{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_k|k)p(k)w_{k^{\prime }}({\mathit{\lambda }}_{k^{\prime }}|{\mathit{\varphi }}_{k^{\prime }},k^{\prime })h(k^{\prime }|k)}\\ =\frac{{\pi }_u(\tilde{\mathit{x}};{\mathit{\lambda }}_{k^{\prime }},k^{\prime }){\pi }_u(\tilde{\mathit{y}};{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_{k^{\prime }}|k^{\prime })p(k^{\prime })w_k({\mathit{\lambda }}_k|{\mathit{\varphi }}_k,k)h(k|k^{\prime })}{{\pi }_u(\tilde{\mathit{y}};{\mathit{\lambda }}_{k^{\prime }},k^{\prime }){\pi }_u(\tilde{\mathit{x}};{\mathit{\lambda }}_k,k)\pi ({\mathit{\lambda }}_k|k)p(k)w_{k^{\prime }}({\mathit{\lambda }}_{k^{\prime }}|{\mathit{\varphi }}_{k^{\prime }},k^{\prime })h(k^{\prime }|k)}.\end{array}$$(5)

This is the Metropolis-Hastings ratio of the move on the augmented state space $\{\tilde{\mathit{y}},{\mathit{\lambda }}_k,k\}$, and the choice of proposal for the auxiliary data ($\tilde{\mathit{y}}\sim f(·;{\mathit{\lambda }}_{k^{\prime }},k^{\prime })$) ensures the ratio $c({\mathit{\lambda }}_{k^{\prime }})/c({\mathit{\lambda }}_k)$ appears in both numerator and denominator and hence cancels. This a RJMCMC version of the exchange algorithm (Murray, Ghahramani, and MacKay 2006), as described in Caimo and Friel (2013). The auxiliary data $\tilde{\mathit{y}}$ are drawn from the same sampling mechanism as the observed data, with parameters set at the candidate state $({\mathit{\lambda }}_{k^{\prime }},k^{\prime })$. Intuitively, the roles of $\tilde{\mathit{x}}$ and $\tilde{\mathit{y}}$ are reversed (or “exchanged”)—the auxiliary data $\tilde{\mathit{y}}$ are offered to the current state $({\mathit{\lambda }}_k,k)$ in exchange for the observed data $\tilde{\mathit{x}}$, in order to persuade a move to the candidate state $({\mathit{\lambda }}_{k^{\prime }},k^{\prime })$.

Therefore, this is a form of reversible jump MCMC algorithm which avoids the need to evaluate the awkward normalizing constants in the same spirit as the exchange algorithm of Murray, Ghahramani, and MacKay (2006), which Caimo and Friel (2013) call the Auto-RJ exchange algorithm.

Algorithm 3

Reversible Jump exchange algorithm

1. for $k=1,\dots ,K$ do

2. Draw samples $\left\{{\mathit{\lambda }}_k\right\}$ from $\pi ({\mathit{\lambda }}_k|\tilde{\mathit{x}},k)$

3. Set ${\mathit{\mu }}_k$ and Σ_k as sample mean and sample covariance of $\left\{{\mathit{\lambda }}_k\right\}$

4. Assign $w_k({\mathit{\lambda }}_k|{\mathit{\varphi }}_k,m_k)\equiv \mathit{\text{MVN}}({\mathit{\mu }}_k,{\Sigma }_k)$ in Step 2 of Algorithm 2

5. end for

6. repeat

7. perform a within-model update of ${\mathit{\lambda }}_k$ using Algorithm 1

8. perform an across-model move using Algorithm 2

9. until desired posterior samples required

The process is summarized in Algorithm 3. In our applications, there are K = 2 competing models, and we set $h(k^{\prime }|k)=1,k^{\prime }\ne k$, so that a move to the other model is proposed with certainty.

4.4 Estimating Model Evidence

In this section we present a within-model approach for estimating the model evidence based on the work of Chib (1995). Suppose that data $\tilde{\mathit{x}}=({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n)$ are assumed to have been generated from a Bingham distribution with parameter $\mathit{\lambda }$. The approach follows from noticing that for any ${\mathit{\lambda }}^{\star }$, (4) implies that $$\pi (\tilde{\mathit{x}})=\frac{\pi (\tilde{\mathit{x}}|{\mathit{\lambda }}^{\star })\pi ({\mathit{\lambda }}^{\star })}{\pi ({\mathit{\lambda }}^{\star }|\tilde{\mathit{x}})}=\frac{{\pi }_u(\tilde{\mathit{x}}|{\mathit{\lambda }}^{\star })}{c{({\mathit{\lambda }}^{\star })}^n}\frac{\pi ({\mathit{\lambda }}^{\star })}{\pi ({\mathit{\lambda }}^{\star }|\tilde{\mathit{x}})}.$$

Hence, to estimate the model evidence $\pi (\tilde{\mathit{x}})$ we require both an estimate of the (normalized) posterior density at ${\mathit{\lambda }}^{\star },\pi ({\mathit{\lambda }}^{\star }|\tilde{\mathit{x}})$, as well as an estimate of the normalizing constant $c({\mathit{\lambda }}^{\star })$. For the latter, we used the method of (Kent 1987), based on asymptotic expansions. Another possibility is the holonomic gradient method of (Sei and Kume 2015), as implemented in the R package hgm. We also used this method, but occasionally encountered some numerical instability in computing the normalizing constant for some of the datasets we simulated. Avoiding such potential problems is one argument in favour of methods which do not require computation of the constant. With regards to estimating the posterior density at ${\mathit{\lambda }}^{\star }$, once we have samples from $\pi (\mathit{\lambda }|\tilde{\mathit{x}})$ using the exchange algorithm of Fallaize and Kypraios (2016), then we can produce a kernel density estimate of the posterior density at ${\mathit{\lambda }}^{\star }$. Of course, in practice, because of the curse of dimensionality this approach can only be used in the case where the number of parameters is fairly small. However, this is not uncommon in the applications where the Bingham distribution is applied to, where observed data are typically points on ${\mathcal{S}}^2$.

4.5 Model Selection via Contour Probabilities

We now consider an alternative approach, which uses posterior contour probabilities (Held 2004) of a suitable parameter in order to provide a measure of the evidence against the simpler model. The probability can be interpreted in a similar fashion to a p-value, and henceforth any reference to p-values will mean the posterior contour probability described below.

Suppose that we are given two independent samples of unit vectors in ${\mathcal{S}}^{q-1},\tilde{\mathit{x}}=({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n)$ and $\tilde{\mathit{y}}=({\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_m)$, respectively. Let us assume that $\tilde{\mathit{x}}$ and $\tilde{\mathit{y}}$ are independent random samples from some Bingham distributions with probability density function $f(\mathit{x};{\lambda }_{\mathit{x}})$ and $f(\mathit{y};{\lambda }_{\mathit{y}})$, respectively, where ${\mathit{\lambda }}_{\mathit{x}}=({\lambda }_1^{(x)},\dots ,{\lambda }_{q-1}^{(x)})$ and ${\mathit{\lambda }}_{\mathit{y}}=({\lambda }_1^{(y)},\dots ,{\lambda }_{q-1}^{(y)})$.

Assessing whether or not the samples $\tilde{\mathit{x}}$ and $\tilde{\mathit{y}}$ come from the same Bingham distribution is equivalent to testing the hypothesis that ${\mathit{\lambda }}_{\mathit{x}}={\mathit{\lambda }}_{\mathit{y}}$ versus ${\mathit{\lambda }}_{\mathit{x}}\ne {\mathit{\lambda }}_{\mathit{y}}$. We propose fitting a Bingham distribution to the combined dataset $\tilde{\mathit{\omega }}=(\tilde{\mathit{x}},\tilde{\mathit{y}})=({\mathit{x}}_1,{\mathit{x}}_2,\dots ,{\mathit{x}}_n,{\mathit{y}}_1,{\mathit{y}}_2,\dots ,{\mathit{y}}_m)$. By independence, the joint density of $\tilde{\mathit{\omega }}$ is ${\prod }_{i=1}^{n+m}f({\tilde{\mathit{\omega }}}_i;{\mathit{\lambda }}^{(i)})$. Assuming a particular structure for the $\left\{{\mathit{\lambda }}^{(i)}\right\}$ allows us to assess the evidence against ${\mathit{\lambda }}_{\mathit{x}}={\mathit{\lambda }}_{\mathit{y}}$. Specifically, let $z_i=\mathbb{I}({\tilde{\mathit{\omega }}}_i\text{is a}\mathit{y}\text{sample}),i=1,\dots ,n+m$, where $\mathbb{I}(·)$ is the indicator function. Then, we endow the $\left\{{\mathit{\lambda }}^{(i)}\right\}$ with the structure (6) $${\mathit{\lambda }}^{(i)}={\mathit{\lambda }}_{\mathit{x}}+z_i\mathit{\alpha },$$(6) where $\mathit{\alpha }={({\alpha }_1,\dots ,{\alpha }_{q-1})}^T\in {\mathbb{R}}^{q-1}$. In other words, the combined data constitute n + m independent observations from the Bingham distribution, where the $\tilde{\mathit{x}}$ samples take parameter value ${\mathit{\lambda }}_{\mathit{x}}=({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_{q-1})$ and the $\tilde{\mathit{y}}$ samples take parameter value ${\mathit{\lambda }}_{\mathit{y}}=({\lambda }_1+{\alpha }_1,{\lambda }_2+{\alpha }_2+\cdots ,{\lambda }_{q-1}+{\alpha }_{q-1})$, with ${\alpha }_1,\dots ,{\alpha }_{q-1}\in \mathbb{R}.$

Testing ${\mathit{\lambda }}_{\mathit{x}}={\mathit{\lambda }}_{\mathit{y}}$ versus ${\mathit{\lambda }}_{\mathit{x}}\ne {\mathit{\lambda }}_{\mathit{y}}$ is equivalent to testing $H_0:{\alpha }_1={\alpha }_2=\cdots ,{\alpha }_{q-1}=0$ versus $H_1:\text{at least one}\alpha \text{is not zero}$. Therefore, we can fit the Bingham distribution to the combined dataset $\tilde{\mathit{\omega }}$ within a Bayesian framework, where each observation ${\tilde{\mathit{\omega }}}_i$ has parameter ${\mathit{\lambda }}^{(i)}$ of the form (6), and then determine if the posterior distribution of the coefficients ${\alpha }_i,i=1,\dots ,q-1$ is “significantly” different from the zero vector. At first glance, one might investigate the corresponding univariate credible intervals of some predefined level separately. However, it is rather unclear how this will translate to a simultaneous probability statement since the coefficients are not independent. One option is to calculate simultaneous credible bands on various levels and determine the smallest level such that the point of interest (i.e., the zero vector in our case) is contained in the simultaneous credible band. Held (2004) proposed instead to use (multivariate) highest posterior density (HPD) regions to check the support of the posterior distribution for a certain parameter vector of interest. Held (2004) defined the (posterior) contour probability $\mathit{\text{Pr}}({\mathit{\alpha }}_0|\tilde{\mathit{\omega }})$ as 1 minus the content of the HPD region of the posterior distribution $\pi (\mathit{\alpha }|\tilde{\mathit{\omega }})$ which just covers ${\mathit{\alpha }}_0$: $$\mathrm{\text{Pr}}({\mathit{\alpha }}_0|\tilde{\mathit{\omega }})=\mathit{\text{Pr}}(\pi (\mathit{\alpha }|\tilde{\mathit{\omega }})\le \pi ({\mathit{\alpha }}_0|\tilde{\mathit{\omega }})|\tilde{\mathit{\omega }}).$$

If the functional form of the posterior density of $\mathit{\alpha }$ is available analytically, then we can obtain an unbiased estimate of the contour probability by Monte Carlo; in other words, by computing the proportion of posterior samples with density value smaller or equal than $\pi ({\mathit{\alpha }}_0|\tilde{\mathit{\omega }})$. In our context, the functional form of the posterior density is unavailable, and therefore we use the modified Rao-Blackwell estimate of the simple Monte Carlo estimator proposed by Held (2004) and implemented in the bayesSurv R package. Using the algorithm that was proposed in Fallaize and Kypraios (2016), one can sample from the posterior distribution of the parameters (${\mathit{\lambda }}_{\mathit{x}},\mathit{\alpha }$), and then determine the posterior simultaneous credible region of the coefficients ${\alpha }_i,i=1,\dots ,q-1$.

4.6 Computational Intensity

In applications of practical interest, observations are unit vectors in ${\mathbb{R}}^3$, which is the situation for our numerical results in Sections 5 and 6. The most intensive case we have considered in our simulations is the reversible jump exchange algorithm with a total of 1000 observations (which is larger than typical axial datasets). Running the sampler for 100,000 sweeps took 4 min on 3.4-GHz cpu using a single core. The other two approaches require only runs of the within-model sampler of Fallaize and Kypraios (2016) which is less intensive. The additional computation of the HPD region for the contour method is minimal, and ran in a matter of seconds using the simult.pvalue function from the bayesSurv R package.

For higher dimensions, that is, data in ${\mathbb{R}}^q$ for larger q than 3, there will of course be increased computational cost. However, our focus in this article is on investigating the use of the computational methods considered, in the cases of most practical importance for real applications of axial data.

5 Simulation Study

The aims of this section are to (i) investigate the performance of the proposed methods and the extent to which they agree and (ii) assess the extent to which it is possible to discriminate between models as parameters and sample sizes vary. We address these aims by performing simulation studies in which datasets of varying size using different parameter combinations are generated, to which we apply our model selection methods.

In this section, and also the real data application in Section 6, the observations are unit vectors in ${\mathbb{R}}^3$, so the Bingham distribution has 2 (eigenvalue) parameters ${\lambda }_1\ge {\lambda }_2\ge 0$ to be inferred. For a prior distribution, we use a product of exponential distributions subject to the above constraint, that is, $$\pi ({\lambda }_1,{\lambda }_2)\propto \hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left\{-{\lambda }_1{\mu }_1-{\lambda }_2{\mu }_2\right\},{\lambda }_1\ge {\lambda }_2\ge 0,$$where μ_i are rate parameters. In all of our numerical results we have used ${\mu }_i=0.01$.

5.1 Comparison

In order to test the performance of the different methods we have proposed, we performed simulations for two different parameter cases. In the first case, we set ${\mathit{\lambda }}_x={\mathit{\lambda }}_y=(20,10)$ and in the second case we set ${\mathit{\lambda }}_x=(20,10)$ and ${\mathit{\lambda }}_y=(20,6)$. In both cases we simulated data for two different sample sizes, $n=m=20$ and $n=m=200$. For the four different combinations of parameter and sample size, we simulated 200 sets of data and performed inference using each of the three approaches we have introduced. The results are presented in Figure 2, where m₁ denotes the event that both populations have the same Bingham distribution. The first row shows the results for ${\mathit{\lambda }}_x={\mathit{\lambda }}_y=(20,10)$ when n = 20 (left) and n = 200 (right). For both sample sizes, there is very good agreement between the RJMCMC and Chib’s method. As expected, for the larger sample size, for most of the datasets the (true) simpler model is favored, whereas there is more variability in $p(m_1)$ for the smaller sample size of 20. For the contour probabilities, with a sample size of 20, the smallest p-value was 0.046, and all other p-values were greater than 0.1 (75% were greater than 0.5). As expected, when the sample size was increased to 200, there were some smaller p-values, but all but 9 were greater than 0.05, as we would expect. In summary, the methods agree well and correctly favor the simpler model.

Fig. 2 Results from the simulation experiments described in Section 5.1. The plots show the model evidence for the simpler model, $p(m_1)$, using Chib’s method and RJMCMC. The bottom row plots also encode information about the contour probabilities. Top row: ${\mathit{\lambda }}_x={\mathit{\lambda }}_y=(20,10)$. Bottom row: ${\mathit{\lambda }}_x=(20,10)$ and ${\mathit{\lambda }}_y=(20,6)$. Left column: $n=m=20$. Right column: $n=m=200$.

The second row shows the results when ${\mathit{\lambda }}_x=(20,10)$ and ${\mathit{\lambda }}_y=(20,6)$. Again, we see good agreement between the two approaches for calculating $p(m_1)$. In the case when n = 20 (left), the majority of the points are concentrated in the upper right of the plot, with $p(m_1)>0.8$, despite the samples coming from different Bingham distributions. In fact, in the next section we will investigate further the ability to discriminate between models for various parameter settings, and see that it is quite hard to distinguish between different Bingham distributions for the parameters used here when sample sizes are small. When n = 200 (right), there is a dense cluster of points in the bottom left of the plot ($p(m_1)<0.2$), but still a fairly uniform spread of values across the remaining values, indicating that it is still reasonably challenging to discriminate between the two different Binghams when the sample size is 200 per group. We also note that, for the points corresponding to larger estimates of $p(m_1)$, the p-values from the contour method are larger too, for example, the points in the upper right of the plot have p-values greater than 0.01. So the three approaches are again in good agreement.

5.2 Model Discrimination

Now that we have a reliable method for determining model evidence, we attempt to answer the following question: how well can we discriminate between models in different parameter and sample size settings? Here, we present some further simulation experiments to investigate this. In particular, we wish to know how well we can discriminate between models when the parameters in each population are close/far apart, for different sample sizes. Specifically, we set ${\mathit{\lambda }}_{\mathit{x}}=(20,10)$ and consider ${\mathit{\lambda }}_{\mathit{y}}$ of the form $(20,{\lambda }_2^{(y)})$, with ${\lambda }_2^{(y)}\in \{2,4,6,8,10\}$. Hence, the cases range from equal parameters through to quite different parameters, in terms of the second eigenvalue. In each case, we also consider varying the sample size, with n = m and $n\in \{20,100,200,500\}$. For each scenario, we simulated 200 datasets, and ran the reversible jump method. Let p₁ be the posterior probability of the simpler model, that is, the model where both populations have the same Bingham distribution. We estimate $P(p_1>0.5)$ via the proportion of the 200 repetitions for which our estimate of p₁ was > 0.5; the results are plotted in Figure 3. We see that, as expected, it becomes easier to discriminate between models as n increases and as the parameters are either very similar or very different. More specifically, it seems that once we have n = 100, there is a high chance of choosing the correct model when ${\lambda }_2^{(y)}=2,4,10$. When n = 200 or 500, most cases are clear, apart from when ${\lambda }_2^{(y)}=8$ in which case there is still a good chance that the simpler model will be selected, even though the parameters are different. For the case n = 200, ${\lambda }_2^{(y)}=6$, there is also a probability of around 0.4 that the simpler model will be preferred, which agrees with our observations in Section 5.1.

Fig. 3 Results from the simulation experiments described in Section 5.2. For each sample size/parameter combination, we simulated 200 datasets and computed the proportion of times that the simpler model was favored, as determined by the rule $p_1>0.5$.

6 Earthquake Data Results

We now turn our attention to the motivating application described in Section 2. Recall that the data consist of three orthogonal axes arising from earthquakes in three situations: 50 observations near Christchurch which took place before a large earthquake on the 22nd of February 2011; 50 observations near Christchurch which took place after this large earthquake; and 32 observations obtained from earthquakes in the north of New Zealand’s South Island. We will label the three samples as CCA (Christchurch before the large earthquake), CCB (Christchurch after the large earthquake) and SI (South Island). We will model the B axes from each earthquake, and assess the evidence that the samples come from the same Bingham distribution.

Fallaize and Kypraios (2016) analyzed these data in a Bayesian framework by fitting Bingham models to the B axes from each of the individual clusters and considering the posterior distributions of the parameters of the diagonal component of the Bingham parameter matrix. Let us denote these parameters from the CCA, CCB, and SI models as ${\lambda }_i^A,{\lambda }_i^B$, and ${\lambda }_i^S$, respectively, i = 1, 2. The authors investigated informally if there is any evidence of a difference between the two Christchurch clusters by considering the bivariate quantity $({\lambda }_1^A-{\lambda }_1^B,{\lambda }_2^A-{\lambda }_2^B)$ and approximating its posterior distribution with a bivariate Normal distribution to visually assess if the origin $(0,0$) was contained comfortably within a 95% probability region. Here, we perform formal model selection to assess the evidence that the samples of B axes come from the same, or different, Bingham distributions. We will perform two pairwise comparisons: CCA-CCB and CCA-SI.

Let the data for the two samples under consideration be stored in matrices $\tilde{\mathit{x}}$ and $\tilde{\mathit{y}}$, where the columns of $\tilde{\mathit{x}}$ are the observed axes ${\mathit{x}}_i\in {\mathcal{S}}^2,i=1,\dots ,n$ and similarly for $\tilde{\mathit{y}}$, which contains m observed axes. The combined data are contained in the $3\times (n+m)$ matrix $\tilde{\mathit{\omega }}$. In practice, we take our data for $\tilde{\mathit{x}}$ to be the eigenvalues ${\tilde{\mathit{\tau }}}^x=\{{\tau }_j^x,j=1,\dots ,3|{\tau }_1^x<{\tau }_2^x<{\tau }_3^x,\Sigma {\tau }_j^x=1\}$ of the sample moment of inertia matrix $\tilde{\mathit{x}}{\tilde{\mathit{x}}}^T/n$, and similarly for $\tilde{\mathit{y}}$. This is equivalent to first rotating to principal axes and taking V to be the identity matrix. (Recall that $A=V\Lambda V^T$ is the Singular Value Decomposition of the Bingham parameter matrix A.) For the combined data (i.e., when fitting one Bingham model to the combined samples) we can then take the data to be the weighted average of ${\tilde{\mathit{\tau }}}^x$ and ${\tilde{\mathit{\tau }}}^y$.

The results are summarized in Table 1. For the CCA-CCB comparison, the model evidence strongly favors the simpler model, whereas for the CCA-SI comparison, the model evidence supports the case that the data come from two different Bingham populations. In terms of the contour probabilities, the same conclusions are reached. This suggests that there is no difference in the pattern of Christchurch earthquakes before and after the large event, but that Christchurch and South Island earthquakes follow different patterns. Regarding a test of equality of two populations based on full axial frame data, (Arnold and Jupp 2013) obtained a p-value of 0.890 for the CCA-CCB comparison and a p-value of less than 0.001 for the CCA-SI comparison, so our results agree with these conclusions.

Table 1 Results from earthquake data.

Data	$p(m_1)$	p-value
CCA-CCB	0.9988	>0.5
CCA-SI	0.0036	0.004

The plots from the contour method are shown in Figure 4. For the CCA-CCB comparison, the point (0, 0) is comfortably contained within the main point cloud, and vice-versa for the CCA-SI comparison, agreeing with the conclusions from the numerical results in Table 1.

Fig. 4 Samples of $({\alpha }_1,{\alpha }_2)$ from the contour method for the CCA-CCB data (left) and CCA-SI data (right). The cross marks the point (0, 0) and the triangle denotes the maximum likelihood estimate.

7 Discussion

In this article, we have introduced three strategies for performing model choice for the Bingham distribution in the Bayesian setting, which is not a straightforward task due to the intractable nature of the normalizing constant. We have illustrated three approaches, two of which avoid the need to approximate this constant, and one which uses an approximation. Application to simulated data show that the methods can successfully discriminate between models, but if the parameters are not very different, relatively large sample sizes are needed to detect a true difference. Analysis of a real dataset shows how formal model choice can be performed in a situation of practical interest to geophysicists.

Although we have focused on the Bingham distribution, the methods could be useful for other distributions in use in directional statistics which have intractable normalizing constants. Also, although the RJMCMC and Chib’s methods are applicable to any number of models in theory, we have focused on choosing between two models in our applications. It will be interesting to investigate performance for other distributions and as the number of competing models increases.

Supplementary Materials

R codeThe R code to run our methods is available in the online supplementary materials. This includes a “readme” file explaining how to run the code for each of our methods. (BinghamCode.tar.gz, GNU zipped tar file).

Disclosure Statement

The authors report there are no competing interests to declare.