Big and complex data sets are collected from various scientific fields such as atmospheric environment, social science, psychological and biomedical studies, bioinformatics, epidemiology, digital imaging information, and machine learning, to name just a few. Big Data can refer to data with big volume or velocity, high-dimensional data (Ahmed, 2017), unstructured or unusual data (Härdle et al., 2018), complex data, etc. Therefore, there is a need for developing statistical techniques other than the traditional analytical frameworks to model, interpret, and use such data in different fields of science. Data with directions are categorized as unusual data that cannot be analyzed and modeled under the Cartesian coordinate system. With the aid of directional statistics, data science meets another level of analytical methods. In this research work, we focus on the analysis of complex directional data. Bayesian methods have received extensive attention in data science because prior information can be added to enhance modeling. Therefore, here, we consider Bayesian analysis of complex directional data.

A robust roadmap with the symmetric Fisher-von Mises-Langevin (FvML) distribution as the key element from a Bayesian perspective is briefly reviewed. According to Kikuchi's collection of directional data (Kikuchi, 1982), the first attention paid to Bayesian methods for directional data was in a paper by Mardia and El-Atoum (1976). They used the Bayesian approach to estimate the location parameter of the FvML distribution when the concentration parameter was known. The author Bagchi made several contributions in this area: (i) Bagchi (1988) formulated a conjugate prior for the mean direction and a non-informative prior for the concentration parameter of the von Mises distribution; (ii) Bagchi and Guttman (1988) focused on Bayesian inference for the multi-variate FvML distribution; (iii) Bagchi and Kadane (1991) derived the Bayes estimate for the cosine of the direction parameter of the von Mises distribution when the concentration parameter is known; and (iv) Bagchi (1994) developed empirical Bayesian techniques to estimate the mean direction of the FvML distribution (see also Guttorp and Lockhart, 1988; Dowe et al., 1996). Damien and Walker (1999) presented a full Bayesian inference for the von Mises distribution implementing a Gibbs sampler while (Rodrigues et al., 2000) provided an empirical or approximate Bayesian inference for the von Mises distribution. Nuñez-Antonio and Gutiérrez-Peña (2005) presented Bayesian analysis of the FvML distribution when all of the parameters were unknown, as well as an algorithm to generate samples from the posterior distribution based on a sampling-importance-resampling method. Muralidharan and Parikh (2007) provided Bayes estimates for both location and concentration parameters of the von Mises distribution. Bhattacharya and SenGupta (2009) presented Bayesian analysis of a generalized von Mises distribution introducing a new algorithm based on importance sampling and Markov chain Monte Carlo (MCMC) to draw samples from the posterior distribution. Mardia (2010) moved the attention to the bivariate von Mises distribution (on the torus) from a Bayesian viewpoint. Infinite mixtures of FvML distributions using standard conjugate priors for the parameters and Dirichlet priors for the mixing probabilities received attention from Bangert et al. (2010). Hornik and Grün (2013) defined conjugate and Jeffreys priors for the FvML distribution while Taghia et al. (2014) worked on Bayesian inference for the FvML mixture model. From 2017 and onwards the following contributions can be highlighted: Straub (2017) presented Bayesian analysis for the FvML distribution in 3D; Røge et al. (2017) presented Bayesian inference in the case of infinite FvML mixture model assumption; Mulder et al. (2020) provided Bayesian inference for mixtures of von Mises distributions using a reversible jump MCMC sampler and focused on non-informative priors. Lastly, the interested reader is referred to Pewsey and García-Portugués (2021) for Bayesian inference of other directional distributions.

Numerous directional data sets tend to show non-trivial features such as skewness. Therefore, the underlying distribution is not always symmetric, which emphasizes the focus on skewed directional distributions. This inspired us to investigate Bayesian analysis for the general class of skew-rotationally-symmetric distributions (Ley and Verdebout, 2017b), an asymmetric extension of all rotationally symmetric distributions, when the location, concentration, and skewness parameters are unknown.

In section 2, the skew-rotationally-symmetric distribution and special cases are described. The novel contribution is given in section 3 where the posterior distributions are obtained for the skew-FvML as the likelihood model, for four different scenarios of the prior distributions for the parameters of the model. Moreover, an algorithm is provided for generating samples from these posterior distributions. The impact of the priors is explored in section 4, by implementing the Wasserstein Impact Measure (WIM). In section 5, a synthetic data analysis is conducted to show the accuracy of the Bayes estimates based on the assumptions of the skew-FvML model. We demonstrate the applicability of the Bayesian framework for well known real datasets in section 6 for dimensions p = 2 and 3.

Most of the distributions on the unit hypersphere 𝕊^p−1 = {v ∈ ℝ^p : v^⊤v = 1}, p ⩾ 2, share the common feature of being rotationally-symmetric about their location μ ∈ 𝕊^p−1. The distribution of a random variable X ∈ 𝕊^p−1 is said to be rotationally-symmetric about μ if for any orthogonal matrix O_p×p satisfying Oμ = μ it is concluded that OX is equal in distribution to X. The FvML distribution, the most common distribution in spherical studies, is obtained by conditioning on the p-variate normal distribution (see Mardia and Jupp, 2000; Ley and Verdebout, 2017a). Suppose X takes values on the non-linear manifold 𝕊^p−1 and has the FvML distribution, then its probability density function (pdf) is given by

where

τ ⩾ 0 is the concentration parameter, and I_α is the modified Bessel function of order α. For τ = 0, (1) simplifies to the uniform distribution. If p = 2 in (1), it results in the von Mises distribution and when p = 3, the Fisher distribution (Fisher, 1953) is obtained.

However, in practice, all real-life phenomena cannot be represented by symmetric models. The interested reader is referred to Downs (2003) for medical research of heart disease diagnosis, Leong and Carlile (1998) for application in neurosciences, Shearman et al. (2000) for biological research on mammalian circadian rhythms, Mardia (2013) and Ameijeiras-Alonso and Ley (2020) for application in bioinformatics especially protein structure prediction, Fisher and Lee (1994) for some studies on wind direction, Ameijeiras-Alonso et al. (2021) for biomechanics studies, and Pewsey (2002) and Ley and Verdebout (2014) for animal movement studies.

Therefore, in this paper, the focus will be on the skew-rotationally-symmetric (SRS) distributions, introduced by Ley and Verdebout (2017b) as

where f(x^⊤μ) is a rotationally-symmetric pdf about μ ∈ 𝕊^p−1, Π : ℝ → [0, 1] is a monotone increasing function satisfying Π(−t) + Π(t) = 1 for all t ∈ ℝ, and ϒμ⊤ represents the semi-orthogonal matrix such that ϒμϒμ⊤=Ip-μμ⊤ and ϒμ⊤ϒμ=Ip-1, where I_p is the p × p identity matrix. The parameter γ ∈ ℝ^p−1 is a skewness parameter vector such that γ = 0 provides the symmetric pdf f(x^⊤μ) and non-zero values of γ provide skewed pdfs. This construction allows using the full potential of existing rotationally-symmetric distributions by turning them into skewed versions.

Substituting (1) in (3) and letting ϒμ⊤x=((1-(x⊤μ)2)1/2Uμ(x)) with U_μ(x) the sign vector which is uniformly distributed on 𝕊^p−2, the skew-FvML (SFvML) distribution is obtained as

where C(τ) is defined in (2).

By using the standard cosine transformation

and choosing ϒμ⊤x=1(x⩾μ)-1(x⩽μ) in (4), for p = 2, the skew-von Mises (SvM) distribution follows as

where θ, μ ∈ (−π, π], τ > 0 and γ ∈ ℝ. Here, the scalar product x^⊤μ is cos(θ − μ). By choosing Π(x)=1+x2, x ∈ [−1, 1], in (6), the sine-skewed von Mises distribution introduced by Abe and Pewsey (2011) is obtained where γ ∈ [−1, 1].

The following lemma can be used to generate a sample from the SFvML distribution.

Lemma 1. (Ley and Verdebout, 2017b) Generate Y from the rotationally-symmetric FvML distribution in (1). Then for any uniformly distributed sign vector U_μ(Y), U_{μ; Π}(Y) is defined as

where U is uniformly distributed on (0, 1) and independent of Y. Then the vector X with pdf (4) is obtained as

In the next section, Bayesian inference with the SFvML distribution as the key element is presented with all location, concentration, and skewness parameters μ, τ, and γ unknown.

Let X = (X₁, X₂, ..., X_n) be a random sample of size n with pdf (4) where the standard normal cumulative density function (cdf) Φ replaces Π. The likelihood function is

Subsequently, four scenarios are presented to define the prior distributions for the location parameter μ, the concentration parameter τ, and the skewness parameter γ.

The prior distributions are a crucial part in Bayesian analysis. If the sample size is small, or available data provide only indirect information about the parameters of interest, the prior distribution becomes more important (Carlin and Louis, 2008). Different criteria can be used for prior selection, we refer the reader to Vehtari et al. (2017). Ghaderinezhad et al. (2022) implemented the Wasserstein Impact Measure (WIM) as a measure of the impact of the choice of the prior in a Bayesian approach. In fact it is a convenient way for quantifying prior impact which will help us to choose between two or more priors in a given situation. Suppose Ω is the vector of parameters and F₁(.) and F₂(.) are two cumulative distribution functions (cdfs) of two posterior distributions π₁(Ω|.) and π₂(Ω|.). The Wasserstein distance between these two posteriors related to two different prior sets is obtained as follows:

with D_Ω the domain of all possible values of Ω. The Wasserstein distance between two posteriors indicates, at any finite sample size n, how close the posterior distributions are and how similar the related inference will be. This is particularly interesting when considering a simple vs. a complicated, computationally intense prior; if the WIM between them is small, then one can safely use the simpler version. When n → ∞ the distance tends to 0.

In this section, a simulation study is conducted to compare the different sets of proposed priors in section 3 for p = 2 using this measure. Since the cdfs of the posteriors in (12)–(15) are not computable, Algorithm 1 and Monte Carlo integration are used to obtain the Wasserstein distance. Also, the transport package (Schuhmacher et al., 2020) in the R software offers functions for computing the Wasserstein distance between two sets of samples from different distributions. Most of the functions in this package have been designed for data with two or higher dimensions. For various combinations of the parameters we draw 200 random observations from the SvM in (6).

To compare the impact of the normal distribution and skew-normal distribution in (8) and (9) (for p = 2) as priors for the skewness parameter γ, the following steps were performed:

Figure 1 (top) shows the calculated Wasserstein distance for different values of λ and n. As expected,

To illustrate the impact of the prior selection for μ and τ the following approach was followed. Assume the normal distribution as the prior for the skewness parameter γ. The posteriors π₂(Ω|.) and π₃(Ω|.), emanating from the informative priors (10) (case 1) and (11) (case 2) were compared with the posterior resulting from the non-informative prior μ ~ Uniform(0, 2π) and π(τ) ∝ 1, denoted by π₁(Ω|.). The posteriors were sampled using Algorithm 1 for n = 10, 15, 20, 25, 30, 35, 40, 50, 100. The Wasserstein distances were calculated between them with the transport package and Monte Carlo method with 1, 000 repetitions.

Figure 1 (middle) illustrates the obtained Wasserstein distance between the posteriors π₁(Ω|.) and π₂(Ω|.) (case 1) and between π₁(Ω|.) and π₃(Ω|.) (case 2) for μ = 2, τ = 1, γ = −1 (left) and μ = 3, τ = 0.6, γ = 1 (right). Figure 1 (bottom) shows the Wasserstein distance between the posteriors π₂(Ω|.) and π₃(Ω|.) for μ = 2, τ = 1, γ = −1 (left) and μ = 3, τ = 0.6, γ = 1 (right), respectively. From Figure 1 (middle and bottom) the following observations can be obtained:

We can thus conclude that from moderate sample sizes on, both priors for all three parameters are rather similar (hence one could use the less computationally intense of both priors), but differ clearly from a non-informative one. In order to judge how large the obtained WIM values actually are, bootstrap re-sampling could be done with the original data; we leave this for future research. Our analysis here is also limited to the chosen parameter values; more simulations need to be done to get a complete picture.

A similar analysis can be performed for p = 3.

In this section, to assess the performance of the Bayesian approach for obtaining the estimates of Ω = (μ, τ, γ), a synthetic data analysis was conducted to obtain the Bayes estimates of the parameters of the SvM distribution (6). We generated samples of size N = 20, 50, 100, 500 from the posterior distributions (12)-(15) (scenarios 1–4) with a burn-in period of 5,000 and k = 500, using Algorithm 1 (the values of these parameters are written down in the respective tables). It is noteworthy that steps 2–7 in Algorithm 1 are combined for scenarios 1 and 2. Bayes estimates of the parameters μ, τ and γ were obtained under the squared error, absolute error and zero-one loss functions by calculating mean, median, and mode of the generated samples, respectively.

The results for p = 2 and p = 3 including the sample mean, standard deviation, quartiles, and mode of the posterior distribution are summarized in Table 1 for each of the scenarios. As can be seen in Table 1 the obtained Bayes estimates are close to the actual values of the parameters. In addition, for small sample sizes our proposed Bayesian approach still provides accurate estimates.

The traceplots of the generated samples from the posteriors, the compare-partial plots and the running mean plots are shown in Figure 2 (p = 2) and Figure 3 (p = 3) for each of the scenarios and p = 2 and 3 using the ggmcmc package in R (Fernández-i-Marın, 2016). A traceplot is an essential plot for evaluating convergence and diagnosing chain problems. It shows the time series of the sampling process and the expected outcome is to get a traceplot that looks completely random. A compare-partial plot provides overlapped kernel density plots that compare the last part of the chain (the last 10% of the values, in green) with the whole chain (in black). Ideally, the initial and final parts of the chain have to be sampling in the same target distribution, so the overlapped densities should be similar. In addition to the traceplots, the running mean plot of the chains is very useful to find within-chain convergence issues. A time series of the running mean of the chain allows to check whether the chain is slowly or quickly approaching its target distribution. The expected output is a line that quickly approaches the overall mean. Figures 2, 3 confirm the convergence of the chains and show that the modified Gibbs sampler recovers the values that actually come from the target posterior distributions.

Running multiple independent chains in parallel is necessary to access the representativeness of the chains. If the multiple chains are not well mixed, the convergence of the chains is suspected (Kruschke, 2014; Vehtari et al., 2021). Therefore, four independent chains were run in parallel for scenario 3 (p = 2) in Table 1 to make the inference more robust and reliable. The results are shown in Figure 4 which confirm the convergency.

To compare the efficiency of Bayes estimates with respect to the maximum likelihood estimations (MLE), the mean squared errors (MSE) of MLEs and Bayes estimates of parameters under the squared error and absolute error loss functions were obtained for scenario 2 and 3 and n = 10, 20, 30, 50, 100 using a Monte Carlo simulation with 500 replications. Then, the relative efficiency (RE) was computed as

where Ω~ is the MLE of Ω = (μ, τ, γ) and Ω^ and Ω^^ are the Bayes estimates of Ω under the squared error and absolute error loss functions, respectively. The results are shown in Figure 5 for scenario 2 (top) and scenario 3 (middle) and μ = 3, τ = 0.6, γ = 1.

From Figure 5 (top and middle) the following general conclusions can be observed:

Finally, to investigate the rule of k in Algorithm 1, the biases of the Bayes estimates of (μ, τ, γ) under the squared error loss function were obtained for scenario 3 (p = 2) in Table 1, n = 100 and different values of k = 10, 50, 100, 200, 300, 500 using a Monte Carlo simulation with 500 replications. The results are shown in Figure 5 (bottom) which demonstrate that, by increasing k, the bias tends to zero and thus, the accuracy of estimates increases.

In what follows, the proposed Bayesian approach's performance for p = 2 is demonstrated through three datasets with different sizes n = 31, 60, 725 with the skew-von Mises distribution in (6) as assumed model. The circular boxplots (Buttarazzi et al., 2018) of the datasets are shown in Figure 6 (top) and confirm the skew pattern of the datasets.

The obtained results in section 5 show all the assumed scenarios provide accurate estimates for the parameters. However, based on the obtained results in section 5 with the WIM, we propose scenario 3 or 4 for obtaining the Bayes estimates to avoid time intensive computations when the sample size is not small (see Figure 7). The justification is that scenarios 1 and 2 need the slice sampler in Algorithm 1 to generate samples from the joint prior (10). Therefore, the Bayes estimates of the parameters μ, τ, and γ were obtained based on scenario 1 for the movement of blue periwinkles dataset (n = 31); scenarios 3 and 4 for the long-axis orientations of feldspar laths (n = 60) and the thunder at Kew (n = 725) datasets, respectively. See below for the description of said datasets.

For p = 3, a dataset of size n = 40 including the expenditures of households is considered. Figure 6 (bottom) shows the scatter plot of the data. For all the datasets to obtain the Bayes estimates we generated samples of size N = 1, 000 from the posterior distributions using Algorithm 1 with a burn-in period of 10,000 and k = 500.

In what follows, we shall describe the individual datasets in detail. Since the conclusion is more or less the same for all p = 2 settings, we already write it down here. It is observed that the proposed Bayesian approach with the skew von Mises distribution as the underlying model provides a good fit to the datasets. Generally, the obtained estimates based on the squared error and absolute error loss functions are more accurate.

Since the assumption that data is rotationally-symmetric is often rejected (Pewsey, 2002; Ley and Verdebout, 2014; Ameijeiras-Alonso and Ley, 2020; Ameijeiras-Alonso et al., 2021), in this paper, we have presented a Bayesian analysis for the skew-rotationally-symmetric FvML distribution. For the first time in Bayesian analysis of directional data the impact of the proposed priors in the set up has been compared using the Wasserstein Impact Measure. Using this measure can give guidance to the practitioner to avoid computationally intensive priors if a simpler prior has similar impact. An algorithm has been used based on modified Gibbs sampling and weighted bootstrap resampling for generating samples from posterior distributions. This coming together of Bayesian methods and skew distributions in the directional domain promises new research interest.

All relevant references for data are contained within the article.

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

This work was based upon research supported in part by the Visiting Professor programme, University of Pretoria and the National Research Foundation (NRF) of South Africa, SARChI Research Chair UID: 71199; Ref.: IFR170227223754 grant No. 109214; Ref.: SRUG190308422768 grant No. 120839, and DSI-NRF Centre of Excellence in Mathematical and Statistical Sciences (CoE-MaSS), South Africa. The opinions expressed and conclusions arrived at are those of the authors and are not necessarily to be attributed to the CoE-MaSS or the NRF. Christophe Ley's research is supported by the FWO Krediet aan Navorsers grant with reference number 1510391N.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.