Doubly intractable models are encountered in a number of fields, e.g. social networks, ecology and epidemiology. Inference for such models requires the evaluation of a likelihood function, whose normalising factor depends on the model parameters and is assumed to be computationally intractable. The normalising constant of the posterior distribution and the additional normalising factor of the likelihood function result in a so-called doubly intractable posterior, for which it is difficult to directly apply Markov chain Monte Carlo methods (MCMC). We propose a signed pseudo-marginal Metropolis-Hastings algorithm with an unbiased block-Poisson estimator to sample from the posterior distribution of doubly intractable models. As the estimator can be negative, the algorithm targets the absolute value of the estimated posterior and uses an importance sampling estimator to ensure simulation-consistent estimates of the posterior mean of a function of the parameters. The importance sampling estimator can perform poorly when its denominator is close to zero. We derive a finite-sample concentration inequality that ensures, with high probability, that this pathological case does not occur. Our estimator for doubly intractable problems has three advantages over existing estimators. First, the estimator is well-suited for efficient parallelisation and vectorisation. Second, its structure is ideal for correlated pseudo-marginal methods, which are well known to dramatically increase sampling efficiency. Third, the estimator enables the derivation of heuristic guidelines for tuning its hyperparameters under simplifying assumptions. We demonstrate the superior performance of our method in the standard benchmark example that models correlated spatial data using the Ising model, as well as the Kent distribution model for spherical data.

Abstract Inclination is the angle of a magnetization vector from horizontal. Clastic sedimentary rocks often experience inclination shallowing whereby syn‐ to post‐depositional processes result in flattened detrital remanent magnetizations relative to local geomagnetic field inclinations. The deviation of recorded inclinations from true values presents challenges for reconstructing paleolatitudes. A widespread approach for estimating flattening factors (f) compares the shape of an assemblage of magnetization vectors to that derived from a paleosecular variation model (the elongation/inclination [E/I] method). Few studies exist that compare the results of this statistical approach with empirically determined flattening factors and none in the Proterozoic Eon. In this study, we evaluate inclination shallowing within 1.1 billion‐year‐old, hematite‐bearing red beds of the Cut Face Creek Sandstone that is bounded by lava flows of known inclination. Taking this inclination from the volcanics as the expected direction, we found that detrital hematite remanence is flattened with whereas the pigmentary hematite magnetization shares a common mean with the volcanics. Using the pigmentary hematite direction as the expected inclination results in . These flattening factors are consistent with those estimated through the E/I method supporting its application in deep time. However, all methods have significant uncertainty associated with determining the flattening factor. This uncertainty can be incorporated into paleomagnetic poles with the resulting ellipse approximated with a Kent distribution. Rather than seeking to find “the flattening factor,” or assuming a single value, the inherent uncertainty in flattening factors should be recognized and incorporated into paleomagnetic syntheses.

This paper presents a Gaussian tracking algorithm with direction-of-arrival (DOA) measurements modelled via the Kent distribution. The key aspect of the algorithm is that the Kent distribution directly models the specific characteristics of DOA measurements in the 3-D space, and can account for different uncertainties in azimuth and elevation, which the von Mises-Fisher distribution cannot. At each update step, the algorithm performs iterated statistical linear regressions. We provide two implementations of the algorithms, one based on sigma-points and the other on analytical linearisation. The effectiveness of the approach is evaluated via numerical simulations.

This paper proposes an efficient numerical integration formula to compute the normalizing constant of Fisher–Bingham distributions. This formula uses a numerical integration formula with the continuous Euler transform to a Fourier-type integral representation of the normalizing constant. As this method is fast and accurate, it can be applied to the calculation of the normalizing constant of high-dimensional Fisher–Bingham distributions. More precisely, the error decays exponentially with an increase in the integration points, and the computation cost increases linearly with the dimensions. In addition, this formula is useful for calculating the gradient and Hessian matrix of the normalizing constant. Therefore, we apply this formula to efficiently calculate the maximum likelihood estimation (MLE) of high-dimensional data. Finally, we apply the MLE to the hyperspherical variational auto-encoder (S-VAE), a deep-learning-based generative model that restricts the latent space to a unit hypersphere. We use the S-VAE trained with images of handwritten numbers to estimate the distributions of each label. This application is useful for adding new labels to the models.

In the original publication of the article, the corrections in Eq. (13) were missed, in which 2v − 1 was changed to 2v in the exponent.

The Fisher–Bingham distribution ( $$\mathrm {FB}_8$$ ) is an eight-parameter family of probability density functions (PDF) on the unit sphere that, under certain conditions, reduce to spherical analogues of bivariate normal PDFs. Due to difficulties in its interpretation and estimation, 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 two datasets using a maximum-likelihood approach and show its improvements over a fit with the more restrictive $$\mathrm {FB}_5$$ distribution.

Robust Feature-Based Point Registration Using Directional Mixture Model

We propose a new distribution for analyzing paleomagnetic directional data, that is, a novel transformation of the von Mises–Fisher distribution. The new distribution has ellipse-like symmetry, as does the Kent distribution; however, unlike the Kent distribution the normalizing constant in the new density is easy to compute and estimation of the shape parameters is straightforward. To accommodate outliers, the model also incorporates an additional shape parameter, which controls the tail-weight of the distribution. We also develop a general regression model framework that allows both the mean direction and the shape parameters of the error distribution to depend on covariates. The proposed regression procedure is shown to be equivariant with respect to the choice of coordinate system for the directional response. To illustrate, we analyses paleomagnetic directional data from the GEOMAGIA50.v3 database. We predict the mean direction at various geological time points and show that there is significant heteroscedasticity present. It is envisaged that the regression structures and error distribution proposed here will also prove useful when covariate information is available with (i) other types of directional response data; and (ii) square-root transformed compositional data of general dimension. Supplementary materials for this article are available online. Code submitted with this article was checked by an Associate Editor for Reproducibility and is available as an online supplement.

Existing parametric regression models in the literature for response data on the unit sphere assume that the covariates have particularly simple structure, for example that they are either scalar or are themselves on the unit sphere, and/or that the error distribution is isotropic. In many practical situations, such models are too inflexible. Here, we develop richer parametric spherical regression models in which the covariates can have quite general structure (for example, they may be on the unit sphere, in Euclidean space, categorical or some combination of these) and in which the errors are anisotropic. We consider two anisotropic error distributions—the Kent distribution and the elliptically symmetric angular Gaussian distribution—and two parametrisations of each which enable distinct ways to model how the response depends on the covariates. Various hypotheses of interest, such as the significance of particular covariates, or anisotropy of the errors, are easy to test, for example by classical likelihood ratio tests. We also introduce new model-based residuals for evaluating the fitted models. In the examples we consider, the hypothesis tests indicate strong evidence to favour the novel models over simpler existing ones.

Abstract We propose novel parametric concentric multi‐unimodal small‐subsphere families of densities forp − 1 ≥ 2‐dimensional spherical data. Their parameters describe a common axis forKsmall hypersubspheres, an array ofKdirectional modes, one mode for each subsphere, andKpairs of concentrations parameters, each pair governing horizontal (within the subsphere) and vertical (orthogonal to the subsphere) concentrations. We introduce two kinds of distributions. In its one‐subsphere version, the first kind coincides with a special case of the Fisher–Bingham distribution, and the second kind is a novel adaption that models independent horizontal and vertical variations. In its multisubsphere version, the second kind allows for a correlation of horizontal variation over different subspheres. In medical imaging, the situation ofp − 1 = 2 occurs precisely in modeling the variation of a skeletally represented organ shape due to rotation, twisting, and bending. For both kinds, we provide new computationally feasible algorithms for simulation and estimation and propose several tests. To the best knowledge of the authors, our proposed models are the first to treat the variation of directional data along several concentric small hypersubspheres, concentrated near modes on each subsphere, let alone horizontal dependence. Using several simulations, we show that our methods are more powerful than a recent nonparametric method and ad hoc methods. Using data from medical imaging, we demonstrate the advantage of our method and infer on the dominating axis of rotation of the human knee joint at different walking phases.

Fracture orientation is a prominent factor in determining the reservoir fluid flow direction in a formation because fractures are the major paths through which fluid flow occurs. Hence, a true modeling of orientation leads to a reliable prediction of fluid flow. Traditionally, various distributions are used for orientation modeling in fracture networks. Although they offer a fairly suitable estimation of fracture orientation, they would not consider any spatial structure for simulated fracture orientations, and would not able to properly reproduce the histograms, and the stereogram of dip and azimuth. To respect this geostatistical and statistical parameters, in this paper a new approach has been presented in which the observed fractures on the image log are firstly clustered, and the major facture families are categorically simulated over the area of study. Afterwards, azimuths are simulated using the probability field obtained from categorical simulation, and dips are conditionally simulated to azimuths. The method is illustrated through a case and the results show that the histograms and stereograms are completely reproduced. In addition, the connectivity of modeled fracture network using the presented method is surveyed in comparison with modeled fracture network using Kent distribution.

We define a distribution on the unit sphere mathbb {S}^{d-1} 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.

ABSTRACTCompositional data are vectors of proportions defined on the unit simplex and this type of constrained data occur frequently in Government surveys. It is also possible for the compositional data to be correlated due to the clustering or grouping of the observations within small domains or areas. We propose a new class of the mixed model for compositional data based on the Kent distribution for directional data, where the random effects also have Kent distributions. One useful property of the new directional mixed model is that the marginal mean direction has a closed form and is interpretable. The random effects enter the model in a multiplicative way via the product of a set of rotation matrices and the conditional mean direction is a random rotation of the marginal mean direction. In small area estimation settings, the mean proportions are usually of primary interest and these are shown to be simple functions of the marginal mean direction. For estimation, we apply a quasi-likelihood method which results in solving a new set of generalized estimating equations and these are shown to have low bias in typical situations. For inference, we use a nonparametric bootstrap method for clustered data which does not rely on estimates of the shape parameters (shape parameters are difficult to estimate in Kent models). We analyze data from the 2009–2010 Australian Household Expenditure Survey CURF (confidentialized unit record file). We predict the proportions of total weekly expenditure on food and housing costs for households in a chosen set of domains. The new approach is shown to be more tractable than the traditional approach based on the logratio transformation.

This paper considers the 3-D spatial fading correlation (SFC) resulting from an angle-of-arrival (AoA) distribution that can be modeled by a mixture of Fisher–Bingham distributions (FB-distributions) on the sphere. By deriving a closed-form expression for the spherical harmonic transform for the component FB-distributions, with arbitrary parameter values, we obtain a closed-form expression of the 3-D SFC for the mixture case. The 3-D SFC expression is general and can be used in arbitrary multiantenna array geometries and is demonstrated for the cases of a 2-D uniform circular array (UCA) and a 3-D regular dodecahedron array (RDA). In computational aspects, we use recursions to compute the spherical harmonic coefficients and give pragmatic guidelines on the truncation size in the series representations to yield machine precision accuracy results. The results are further corroborated through numerical experiments to demonstrate that the closed-form expressions yield the same results as significantly more computationally expensive numerical integration methods.

We apply the holonomic gradient method to compute the distribution function of a weighted sum of independent noncentral chi-square random variables. It is the distribution function of the squared length of a multivariate normal random vector. We treat this distribution as an integral of the normalizing constant of the Fisher–Bingham distribution on the unit sphere and make use of the partial differential equations for the Fisher–Bingham distribution.

Geochemical surveys collect sediment or rock samples, measure the concentration of chemical elements, and report these typically either in weight percent or in parts per million (ppm). There are usually a large number of elements measured and the distributions are often skewed, containing many potential outliers. We present a new robust principal component analysis (PCA) method for geochemical survey data, that involves first transforming the compositional data onto a manifold using a relative power transformation. A flexible set of moment assumptions are made which take the special geometry of the manifold into account. The Kent distribution moment structure arises as a special case when the chosen manifold is the hypersphere. We derive simple moment and robust estimators (RO) of the parameters which are also applicable in high-dimensional settings. The resulting PCA based on these estimators is done in the tangent space and is related to the power transformation method used in correspondence analysis. To illustrate, we analyze major oxide data from the National Geochemical Survey of Australia. When compared with the traditional approach in the literature based on the centered log-ratio transformation, the new PCA method is shown to be more successful at dimension reduction and gives interpretable results.

Remotely sensed data products are now routinely used to study various aspects of the Earth's atmosphere. These remote sensing datasets are typically very high dimensional, have near global coverage and exhibit nonstationary spatial correlation structures. Proper statistical analysis of these datasets should be sufficiently flexible to account for all these aspects. To this end, we develop a kernel convolution construction of spatial processes on a sphere. As is the case with kernel convolution constructions on the plane, we establish a link between stationary kernels and a stationary covariance function on the sphere via the spherical harmonic decomposition of the kernel. We also introduce the Kent distribution as an appropriate kernel with interpretable parameters to be used in the kernel convolution construction. We demonstrate the discrete kernel convolution model using a dataset of remotely sensed CO2concentrations over the globe. Copyright © 2014 John Wiley & Sons, Ltd.

We propose an accelerated version of the holonomic gradient descent and apply it to calculating the maximum likelihood estimate (MLE) of the Fisher---Bingham distribution on a $$d$$ d -dimensional sphere. We derive a Pfaffian system (an integrable connection) and a series expansion associated with the normalizing constant with an error estimation. These enable us to solve some MLE problems up to dimension $$d=7$$ d = 7 with a specified accuracy.

In an earlier paper Kume & Wood (2005) showed how the normalizing constant of the Fisher--Bingham distribution on a sphere can be approximated with high accuracy using a univariate saddlepoint density approximation. In this sequel, we extend the approach to a more general setting and derive saddlepoint approximations for the normalizing constants of multicomponent Fisher--Bingham distributions on Cartesian products of spheres, and Fisher--Bingham distributions on Stiefel manifolds. In each case, the approximation for the normalizing constant is essentially a multivariate saddlepoint density approximation for the joint distribution of a set of quadratic forms in normal variables. Both first-order and second-order saddlepoint approximations are considered. Computational algorithms, numerical results and theoretical properties of the approximations are presented. In the challenging high-dimensional settings considered in this paper the saddlepoint approximations perform very well in all examples considered. Copyright 2013, Oxford University Press.

Compositional data can be transformed to directional data by the square root transformation and then modelled by using the Kent distribution. The current approach for estimating the parameters in the Kent model for compositional data relies on a large concentration assumption which assumes that the majority of the transformed data is not distributed too close to the boundaries of the positive orthant. When the data is distributed close to the boundaries with large variance significant folding may result. To treat this case we propose new estimators of the parameters derived based on the actual folded Kent distribution which are obtained via the EM algorithm. We show that these new estimators significantly reduce the bias in the current estimators when both the sample size and amount of folding is moderately large. We also propose using a saddlepoint density approximation for the Kent distribution normalising constant in order to more accurately estimate the shape parameters when the concentration is small or only moderately large.