Multimodal Distributions for Circular Axial Data

The density function of a circular random variable $\phi\in(0,2\pi]$, based on non-negative trigonometric sums (NNTS) with a complex parameter vector $\underline{c}$, is defined as

where $i=\sqrt{-1}$, and $c_{k}=c_{re,k}+ic_{im,k}$ are complex parameters for $k=0,1,\dots,M$, with $c_{re,k}$ and $c_{im,k}$ being the real and imaginary parts, respectively. Here, $e^{ik\phi}=\cos(k\phi)+i\sin(k\phi)$. The NNTS density is thus the squared norm of a complex sum. To define a valid density function, the parameters must satisfy

where $||c_{k}||^{2}=c_{re,k}^{2}+c_{im,k}^{2}$. This implies that $c_{im,0}=0$, making $c_{0}$ a positive real number.

A density function for an axial random variable $\theta\in[0,\pi)$ must satisfy the condition $f(\theta)=f(\theta\pm\pi)$ since an axis represents an undirected unit vector. If $\phi$ is the angle of a unit vector on a circle, then the density should have identical values at $\phi$ and $\phi\pm\pi$, representing opposite directions. Thus, the condition for the NNTS circular density becomes

Evaluating this condition:

Since $e^{i(\pm k\pi)}=\cos(k\pi)$, and $\cos(k\pi)=1$ for even $k$ and $\cos(k\pi)=-1$ for odd $k$, the condition $f_{NNTS}(\phi)=f_{NNTS}(\phi\pm\pi)$ implies:

To satisfy these conditions, we must have $c_{k}=0$ for all odd $k$. Thus, the NNTS density for an axial random variable $\theta\in[0,\pi)$, with a complex parameter vector $\underline{v}=(v_{0},v_{1},\dots,v_{M})^{T}$, is defined by considering only even indexes in the summation of Equation (1):

This corresponds to doubling the original circular angle $\phi$, a common technique for analyzing axial data. For this to be a valid density, the parameters must satisfy

implying $v_{0}$ must be a non-negative real number. Alternatively, the NNTS axial density can be expressed as:

with $\sum_{k=0}^{M}||v_{k}||^{2}=1$, placing parameters on the unit complex hypersphere $CS^{M}=\{\underline{v}:\sum_{k=0}^{M}||v_{k}||^{2}=1\}$.

The number of free parameters in $\underline{v}$ equals $2M$, and $M$ determines the maximum number of modes for the NNTS axial density. The support of the axial random variable can be $(0,\pi]$ due to $f_{NNTSaxial}(\theta\mid\underline{v})=f_{NNTSaxial}(\theta\pm\pi\mid\underline{v})$, although the interval $[-\frac{\pi}{2},\frac{\pi}{2})$ is also common in literature. Note that NNTS axial models are nested; models with smaller $M^{*}$ are special cases of those with larger $M^{**}$, facilitating likelihood ratio tests. Specifically, the case $M=0$ corresponds to the uniform axial density.

The trigonometric moments, $E(e^{ir\theta})$, for $r$ a non-negative integer, of the NNTS axial distribution in Equation (4), are given by:

The integral $\int_{0}^{\pi}e^{i(2(k-m)+r)\theta}\,d\theta$ equals zero when $2(k-m)+r\neq 0$, and equals $\pi$ when $2(k-m)+r=0$. Since $k$ and $m$ are integers, the trigonometric moment is zero for odd $r$ and simplifies to:

for $r=2,4,\ldots,M$ when $M$ is even, and $r=2,4,\ldots,M-1$ when $M$ is odd. From these trigonometric moments, one can derive other characteristics of the NNTS axial distribution, such as the mean direction, circular variance, asymmetry coefficient, and kurtosis (see Mardia and Jupp, 2000).

In many practical applications, it is necessary to model axial data that are symmetric with respect to an angle $\mu\in(0,\pi]$, known as the axis of symmetry. A NNTS axial model is symmetric if the complex parameter vector $\underline{v}$ is restricted to be real (see Fernández-Durán and Gregorio-Domínguez, 2025). The corresponding density function is defined as:

Here, $\underline{v}_{R}$ is a vector of real coefficients. Alternatively, the complex vector $\underline{v}^{*}$ can be defined with $v_{k}^{*}=v_{Rk}e^{-i2k\mu}$. The total number of free parameters in a symmetric NNTS axial model is $M+1$.

Since the proposed NNTS axial density can be expressed as a quadratic form,

with $\underline{v}=(v_{0},v_{1},\ldots,v_{M})^{T}$ as the parameter vector and $\underline{e}=(1,e^{i2\theta},e^{i4\theta},\ldots,e^{i2M\theta})^{T}$ as the vector of trigonometric statistics, the estimation algorithm described by Fernández-Durán and Gregorio-Domínguez (2010) can be used. This involves a modified Newton algorithm on the complex unit hypersphere (manifold) $CS^{M}=\{\underline{v}:\sum_{k=0}^{M}||v_{k}||^{2}=1\}$. It is recommended to have at least $7M^{*}$ observations to reliably fit an NNTS axial model with $M=M^{*}$.

Since the model parameters are estimated via maximum likelihood, likelihood ratio tests can be defined as

where $\underline{\theta}$ is the vector of observed axial data, $\hat{\underline{v}}_{G}$ is the MLE under the general model, $\hat{\underline{v}}_{R}$ is the MLE under the restricted model, and $\hat{\ell}(\cdot)$ denotes the corresponding maximized log-likelihood values. Under standard regularity conditions, $LLR_{RG}$ asymptotically follows a chi-squared distribution with degrees of freedom equal to the number of constraints imposed to derive $\underline{v}_{R}$ from $\underline{v}_{G}$.

For instance, a test for symmetry can be performed by comparing the general NNTS axial model to a symmetric version. In this case, the test statistic follows a chi-squared distribution with $M-1$ degrees of freedom. To test for homogeneity across $P$ different populations with independent samples, one compares the sum of the maximized log-likelihoods from each population to the log-likelihood obtained from pooling all data into a single model. The degrees of freedom for the test are given by:

where $2M_{k}$ is the number of parameters in the fitted model for population $k$, and $2M$ is the number of parameters in the pooled model.

Feldspar crystals within rocks often align in parallel or subparallel patterns, providing valuable insights into the processes of rock formation and historical geological events, such as volcanic eruptions. Dataset B2 from Fisher (1993) contains measurements of the long-axis orientations of 133 feldspar laths in basalt. Table 1 (columns 2 to 4) presents maximized log-likelihood, Bayesian Information Criterion (BIC), and Akaike Information Criterion (AIC) values for fitted NNTS axial models with parameter $M$ ranging from 0 to 13. Figure 1 (left plot) shows a circular dot plot and histogram of the dataset along with the best-fitted NNTS axial densities selected by BIC ($M=0$) and AIC ($M=4$). A uniformity test by Fernández-Durán and Gregorio-Domínguez (2024), applied by doubling the axial angle, rejects the null hypothesis of uniformity implied by the BIC-selected model ($M=0$) in favor of the multimodal NNTS axial model chosen by AIC ($M=4$), yielding a p-value of less than 0.01.

Dataset B5 from Fisher (1993) contains orientations of an additional 60 feldspar laths in basalt. Fisher emphasized the importance of testing uniformity and determining the number of modes in this dataset. Table 1 (columns 5 to 7) presents maximized log-likelihood, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) values for NNTS axial models with parameter $M$ ranging from 0 to 6. Figure 1 (center plot) shows the histogram along with the best-fitting BIC and AIC models, both of which correspond to a unimodal and symmetric density with $M$=1. We conducted a likelihood ratio test comparing a symmetric NNTS model with $M=1$ against an asymmetric NNTS model with $M=2$, depicted as a dashed line in Figure 1 (center plot). The resulting p-value was 0.6077, indicating insufficient evidence to reject the simpler symmetric model ($M=1$) in favor of the asymmetric alternative ($M=2$). When comparing the symmetric and asymmetric NNTS axial models with ($M=2$) using a likelihood ratio test for symmetry, the null hypothesis of symmetry is rejected at a 5% significance level with a p-value equal to 0.0459.

Additionally, the uniformity test by Fernández-Durán and Gregorio-Domínguez (2024) produced a p-value between 0.01 and 0.05, leading to the rejection of the null hypothesis of uniformity at a 5% significance level when comparing the uniform density model ($M=0$) with the best AIC- and BIC-selected NNTS axial models with $M=1$. Our results are consistent with those of Fisher (1993), supporting both the non-uniformity and unimodality of the dataset.

Face cleats are longitudinal fractures that represent significant geological features within coal seams-structures composed of coal deposits embedded within layers of rock. Variations in the orientation of these fractures can indicate hazardous mining conditions, including the potential presence of gas. Dataset B22 from Fisher (1993) includes 63 median directions of face cleats measured at the Wallsend Borehole Colliery in New South Wales, Australia. The median directions were measured at intervals of 20 meters along the coal seam. Consequently, the dataset exhibits sequential ordering, implying that the observations may not be independent. Despite the presence of potential dependence structures within the dataset, we analyzed the data as if it were a random sample to facilitate a direct comparison with results presented by Arnold and SenGupta (2006).

Table 1 (columns 8 to 10) provides maximized log-likelihood, Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) values for fitted NNTS axial models with parameter $M$ ranging from 0 to 7. Figure 1 (right plot) shows the histogram along with the best-fitting AIC and BIC models, both of which select $M=6$. The selected NNTS model with $M=6$ displays two prominent modes within the observed range of data. This result contrasts with findings by Arnold and SenGupta (2006), who considered only unimodal densities for the same dataset.

Dataset B13 from Fisher (1993) contains measurements of termite mound orientations for the species Amitermes laurensis, collected from 14 sites on the Cape York Peninsula, North Queensland, Australia. Observed orientations were standardized to the interval between 0 and $\pi$ by subtracting $\pi$ from any measurement greater than $\pi$. Figure 2 displays histograms for the orientation data at each of the 14 sites, along with sample sizes and geographical coordinates (latitude and longitude in degrees). The figure also shows the best-fitting NNTS axial densities according to the Bayesian Information Criterion (BIC) for each site. The final plot in Figure 2 presents the circular dot plot, histogram and the best BIC-fitted NNTS axial density for the combined data from all 14 sites, which serves as the basis for a likelihood ratio test of homogeneity, with the null hypothesis stating that the NNTS axial density is the same across all 14 sites. Table 2 provides maximized log-likelihood, AIC, and BIC values for NNTS axial densities with $M=0,1,\ldots,5$.

Applying a likelihood ratio test for homogeneity using the best-fitting BIC-selected NNTS models for the 14 individual sites compared to the best-fitting BIC NNTS model for the combined dataset yields a p-value of $2.3648\times 10^{-8}$, assuming an asymptotic chi-squared distribution with 140 degrees of freedom. An alternative test, comparing the fitted NNTS models with $M=5$ at each site against a joint model, produces a p-value of 0.0062 with 130 degrees of freedom. In both tests, the null hypothesis of homogeneity across the 14 sites is rejected, indicating that variations exist in the orientation densities of termite mounds depending on their geographic location (latitude and longitude). These findings confirm the results obtained by Fisher (1993), who performed separate tests for the equality of mean directions and dispersions using a von Mises density model by doubling the axial angles to analyze them as circular data. However, it should be noted that the von Mises model is inherently unimodal and therefore unsuitable for modeling multimodal datasets, such as several of those shown in Figure 2, particularly evident in the histograms and fitted NNTS axial densities for site 8.

Fisher (1993) conducted tests of homogeneity specifically for pairs of sites 5 and 14, and sites 6 and 8. Using a von Mises density model by doubling the axial angles to analyze them as circular data, Fisher rejected the null hypothesis of equal mean directions for both pairs. However, for the pair of sites 6 and 8, Fisher concluded that the hypothesis of equal dispersions could not be rejected. To comprehensively compare our results with Fisher’s (1993) analysis, we performed likelihood ratio homogeneity tests for the same pairs (sites 5 and 14, and sites 6 and 8) using an NNTS axial model with parameter $M=4$. The resulting p-values were 0.0051 for sites 5 and 14, and 0.0010 for sites 6 and 8, leading us to reject the null hypothesis of homogeneity for both pairs. These outcomes align with Fisher’s results, although Fisher’s tests separately evaluated mean directions and dispersions, while our approach employs a single integrated likelihood ratio test of overall homogeneity.

Using satellite imagery, deer bedding patterns in snow, and field observations, Begall et al. (2008) analyzed the body orientations of domestic cattle, red deer, and roe deer during grazing and resting, with respect to true (geographic) North, at various locations worldwide with differing magnetic declinations. Their analysis suggested that the primary factor influencing the animals’ orientations is alignment with Earth’s magnetic field, even after accounting for factors such as site-specific declination, wind direction, solar position, and temperature. For each of the three samples, Begall et al. (2008) clearly rejected the null hypothesis of uniformity using the Rayleigh uniformity test. We utilized the dataset from Begall et al. (2008), consisting of mean axial vector directions for herds observed at 308 locations for domestic cattle, 40 locations for red deer, and 201 locations for roe deer. Using these data, we fitted NNTS axial models and conducted homogeneity tests on the orientations across the three populations. Particular attention was paid to assessing homogeneity between red and roe deer populations.

Table 3 shows maximized log-likelihood, AIC, and BIC values for fitted NNTS axial models with different values of $M$ (from 0 to 5 for red deer, and from 0 to 10 for cattle and roe deer). When testing homogeneity across the three populations, we clearly rejected the null hypothesis, obtaining p-values of $4.85\times 10^{-22}$ (best BIC model: $M=1$ for cattle, $M=3$ for red deer, $M=6$ for roe deer, and $M=4$ for all three populations combined) and $8.31\times 10^{-22}$ (best AIC model: $M=2$ for cattle, $M=6$ for red deer, $M=5$ for roe deer, and $M=5$ for all three populations combined).

When testing homogeneity specifically between red deer and roe deer, we obtained p-values of 0.0511 (best BIC model with $M=6$) and 0.9999 (best AIC model with $M=6$). Figure 3 displays circular dot plots, histograms, and fitted density curves based on the best BIC and AIC NNTS axial models for cattle, red deer, and roe deer. For each population, the modes correspond to the average magnetic declination of the various sampling locations worldwide, as reported by Begall et al. (2008).

The spatial orientation of plant leaves, particularly leaf inclination angles, provides essential information regarding key biological processes such as photosynthesis efficiency, reflectance, temperature regulation, and overall resource utilization. Leaf inclination angle distributions vary according to plant species and geographic location. Recent advancements in digital imaging technology have significantly increased the availability of leaf inclination angle measurements, as algorithms can now efficiently determine these angles directly from digital images (see Pisek and Adamson, 2020).

In botanical studies, the leaf inclination angle $\theta_{leaf}$ typically ranges from 0 to $\frac{\pi}{2}$. It is common practice in the literature to model the transformed linear variable $X=\frac{2\theta_{leaf}}{\pi}$ using a Beta distribution with parameters $\alpha$ and $\beta$, and a probability density proportional to $x^{\alpha-1}(1-x)^{\beta-1}$. Based on parameter estimates for different plant species, de Wit (1965) proposed six theoretical distributions for leaf orientation angles: planophile (leaves predominantly horizontal), erectophile (leaves predominantly vertical), plagiophile (leaves predominantly inclined), extremophile (leaves predominantly horizontal and vertical), spherical (leaves oriented isotropically like a sphere), and uniform. The parameters of the Beta distribution for the plagiophile case satisfy $\alpha=\beta$ with $\alpha,\beta>1$, whereas for the extremophile case, they satisfy $\alpha=\beta$ with $\alpha,\beta<1$, indicating that the Beta distribution takes equal values at the extremes 0 and 1 (see Chianucci et al., 2018). For the uniform distribution, the parameters satisfy $\alpha=\beta=1$. Therefore, the angle $\theta=2\theta_{leaf}$ can be treated as an axial angle, and the distributions of $\theta$ for plants exhibiting plagiophile, extremophile, or uniform leaf distributions can be modeled using NNTS axial models.

The dataset from Pisek and Adamson (2020) contains leaf inclination angle measurements for 71 species of gum trees (genus Eucalyptus). From these, we selected 11 species observed at the Huntington Library, Art Collections, and Botanical Gardens in Pasadena, California (latitude 34.125, longitude -118.114, altitude 207 m.a.s.l.), classified by Pisek and Adamson (2020) as either uniform (9 species) or plagiophile (2 species). Table 4 provides species names, sample sizes, Beta distribution parameter estimates, the optimal $M$ values for the best-fitting NNTS axial models based on the Bayesian Information Criterion (BIC), and the corresponding BIC values. Figures 4 and 5 display circular dot plots, histograms, and the best BIC-fitted NNTS axial densities for species classified as uniform (U) and plagiophile (PG), respectively.

Notably, for species classified as uniform (cases 1 balladoniensis and 10 shirleyi), neither the best BIC- nor AIC-selected NNTS axial densities corresponded to a uniform distribution (Figure 4). This observation is further confirmed by the uniformity test of Fernández-Durán and Gregorio-Domínguez (2024) shown in Table 4, where the null hypothesis of uniformity is rejected for cases 1 and 10 when comparing against either a NNTS axial model with $M=1$ or the best BIC-selected NNTS axial model. Moreover, when using the best AIC-selected NNTS axial model as an alternative, the hypothesis of uniformity is additionally rejected for cases 2, 6, and 8 at a 5% significance level.