BAMBI: An R package for Fitting Bivariate Angular Mixture Models
Saptarshi Chakraborty Samuel W. K. Wong






§ INTRODUCTION
 

Statistical analyses of angular or directional data have found applications in a variety of fields, such as geology (earth's magnetic poles), meteorology (wind directions) and bioinformatics (backbone structures of proteins). Directional data can be univariate or multivariate, and one way of representing such data is via angles measured on a circle [0, 2π) (element-wise when multivariate), and hence the name angular. Angular methods are also applicable to any interval that wraps around (e.g., [0, L) or [-L/2, L/2) for some L > 0) when transformed to the circle [0, 2π).  The wraparound condition on the support invalidates direct applicability of many standard statistical methods.  There is substantial literature devoted to the development of descriptive and inferential techniques for directional data (see, e.g.,  <cit.>), with the traditional univariate case as the primary focus, although the bivariate case is gaining increasing interest <cit.> along with the emergence of new applications. Bivariate angular data can now be found in  a variety of modern scientific problems, with many notable applications arising from the field of computational biology <cit.>. A major area of research in protein bioinformatics involves modeling and predicting protein 3-D structures, which requires proper handling of the paired backbone torsion angles. Formal analyses of these bivariate angle pairs thus require rigorous statistical techniques and models.

A unique feature in the modeling of  directional data is the use of angular probability distributions, or  mixtures thereof (see Section <ref>), which are inherently different from their linear (Euclidean) counterparts because of the wraparound nature of their supports.  Bayesian methods provide flexible tools for analyzing and modeling  such data. First, one may incorporate prior information, if available, into modeling. Second, one may use powerful computational methods, i.e., Markov chain Monte Carlo (MCMC, see Section <ref>) for sampling from the posterior, to fit such models and assess the fitted models.  Third, one may readily compute posterior quantities of interest while coherently accounting for uncertainty in the model parameters.  Within this context, this package was developed for fitting Bivariate Angular Mixtures using Bayesian Inference (hence, the package name BAMBI).  In we implement the two most popular angular distributions, namely the wrapped normal (or Gaussian) and the von Mises distributions, and consider both univariate and bivariate versions of these.  provides functionality for modeling univariate and bivariate angular data using these distributions, and for fitting finite mixture models of these distributions.   We first introduce the basics of these distributions and mixture models. It should be noted that the bivariate distributions considered in this paper have support [0, 2π)^2 (i.e., on a torus), which are distinct from those defined on the surface of the unit sphere, such as the von Mises - Fisher distribution.



 §.§ Wrapped Normal Distributions
 
For univariate continuous data, the angular analogue of the normal distribution on the real line is the wrapped normal distribution obtained by wrapping a normal random variable around the unit circle  (see, e.g., <cit.>). Formally, let X be a normal random variable with mean μ and variance σ^2 > 0. Then the distribution of ψ = X  2π is called the wrapped normal distribution with mean μ and variance σ^2 and is denoted by (μ, σ^2). The density of ψ∼(μ, σ^2) is given by:

    f_(ψ|μ, σ) = 1/σ√(2π)∑_ω∈exp[-1/2 σ^2(ψ - μ - 2πω)^2 ];   ψ∈ [0, 2π)

where  denotes the set of all integers. Since the density contains a summation over entire , without loss of generality, we let μ∈ [0, 2π) to ensure identifiability. Figure <ref> displays the univariate wrapped density with μ = π and κ = 0.01, 1 and 10, which shows that the density is symmetric around μ and becomes more concentrated as κ increases.


The multivariate generalization of the above distribution is straightforward <cit.>.  The distribution of a random vector = (ψ_1, ⋯, ψ_p)^⊤ with probability density

    1/√(|Σ|(2π)^p)∑_∈^pexp[-1/2( -  - 2π)^⊤Σ^-1( -  - 2π) ];   ∈ [0, 2π)^p

with ∈ [0, 2π)^p and Σ  positive definite, is called the p-variate wrapped normal distribution with mean vector  and variance matrix Σ, denoted by _p(, Σ).  Although (<ref>) and (<ref>) are the most common parameterizations of the wrapped normal distributions found in the literature, to facilitate comparability with the von Mises distribution (defined in Section <ref>), we shall use the equivalent representation(s) obtained through the re-parameterization(s) κ = 1/σ^2 and Δ = Σ^-1.  BAMBI handles the univariate and bivariate cases, namely p = 1 and p = 2. Thus, the form of the univariate  wrapped normal density we use is

    f_(ψ|μ, κ) = √(κ/2π)∑_ω∈exp[-κ/2(ψ - μ - 2πω)^2 ];   ψ∈ [0, 2π)

with μ∈ [0, 2π) and κ > 0; and that of the bivariate density is

    f__2(ψ_1, ψ_2 | μ_1, μ_2, κ_1, κ_2, κ_3) 
       = √(κ_1 κ_2 - κ_3^2)/2π∑_(ω_1, ω_2) ∈^2exp[-1/2{κ_1 (ψ_1 - μ_1 - 2πω_1)^2 + κ_2 (ψ_2 - μ_2 - 2πω_2)^2  . .
                                     . . + 2 κ_3 (ψ_1 - μ_1 - 2πω_1) (ψ_2 - μ_2 - 2πω_2) }]

where ψ_1, ψ_2, μ_1, μ_2 ∈ [0, 2π), κ_1, κ_2 > 0 and κ_3^2 ≤κ_1κ_2, obtained by letting = (μ_1, μ_2)^⊤ and

    Δ = [ κ_1 κ_3; κ_3 κ_2 ].

Figure <ref> plots the univariate wrapped normal densities with μ = π and κ = 0.01, 1 and 10, which shows that the density is symmetric around μ and becomes more concentrated as κ increases. 


Similarly, the bivariate wrapped normal  density is also symmetric around (μ_1, μ_2) and becomes more concentrated as κ_1 and/or κ_2 increases, while the parameter κ_3 regulates the association between the random coordinates.
This can be visualized from Figure <ref> displaying the surfaces of the density created via BAMBI function surface_model,  for different parameter combinations. (Codes for generating these plots can be found in the replication R script for this paper.) The upper panels of Figure <ref> show how the density becomes more concentrated when κ_1 and κ_2 are increased (while keeping κ_3 fixed). In contrast,  the lower panels of Figure <ref> display density surfaces showing how the association between the random coordinates changes (from positive to negative), when κ_3 is changed (from negative to positive, since κ_3 is the diagonal element of the inverse covariance matrix) while keeping κ_1 and κ_2 are fixed. 



Note that when κ→ 0 (or Δ→ 0_2 × 2) then the distribution of ψ = X  2π converges to the uniform distribution over [0, 2π) (or [0, 2π)^2). Hence, we shall include the cases κ = 0 and κ_1 = κ_2 = κ_3 = 0 in the support of these parameters, and define the associated densities by their limits.





The precision parameter κ (κ_1, κ_2 in the bivariate case) is (are) conceptually similar to the concentration parameters in the von Mises distribution (see Section <ref>). Therefore to aid comparability, we shall call κ (κ_1 and κ_2) the concentration parameter(s) of the univariate (bivariate) wrapped normal model. In BAMBI, evaluation of univariate and bivariate wrapped normal densities are implemented through the function dwnorm and dwnorm2 respectively. Random data  from these models can be generated using rwnorm and rwnorm2 respectively. 




 §.§ von Mises distributions
 
Wrapped normal models have a high computational cost in practice. Although the sum over  in the expression for the density can be well-approximated by a sum over the set A = {-3,-2,-1,0,1,2,3} (i.e., 3 integer displacements, covering ± 3 standard deviations from the mean), it can be seen that the number of terms in the sum grows exponentially as the dimension increases. For instance, in the bivariate case, even if  is approximated  by set A, the (double) sum in the density  consists of 49 terms.

Because of this difficulty, the von Mises distribution is an alternative that is widely used; it is able to approximate the wrapped normal while being less computationally intensive <cit.>. A random variable ψ is said to follow the von Mises distribution (also called  the circular normal distribution, <cit.>)  with mean parameter μ and concentration parameter κ, denoted ψ∼(μ, κ), if ψ has the density

    f_(ψ|μ, σ) = 1/2π I_0(κ)exp(κcos(ψ-μ));   ψ∈ [0, 2π)

where μ∈ [0, 2π), κ≥ 0 and I_r(·) denotes the modified Bessel function of the first kind and order r.  Letting κ = 0 makes (<ref>) the uniform density over [0, 2π), and when κ→∞, (<ref>) converges to a normal density. An intuitive explanation of the latter result follows from the fact that when  the concentration parameter κ is large, ψ - μ≈ 0, so that cos (ψ - μ) ≈ 1-(ψ-μ)^2/2, which makes the exponent in the density (<ref>) approximately proportional to the  (μ, (1/√(κ))^2) density. A formal proof can be found in  <cit.>. 

Figure <ref> plots the von Mises densities with μ = π and κ = 0.01, 1 and 10, which shows that the density is symmetric around μ and becomes more concentrated as κ increases, and that the density is broadly similar to the associated univariate wrapped normal density.

A multivariate generalization for the univariate von Mises distribution is however not as straightforward as the wrapped normal distribution, as there is not a unique way of defining a  multivariate distribution with univariate von Mises-like marginals.  In the bivariate case, two versions of the bivariate von Mises distribution have been suggested for practical use, namely the sine model <cit.> and the cosine model <cit.>. They are comparable to the bivariate normal model both in terms of number of parameters (five), and the interpretability of those parameters.  Other generalizations with more parameters have been studied theoretically <cit.>.

Let = (ψ_1, ψ_2)^⊤ be a random vector on ^2 with support [0, 2π)^2. Then  is said to follow the (bivariate) von Mises sine distribution with mean parameters μ_1, μ_2, concentration parameters κ_1, κ_2, and association parameter κ_3, denoted ∼(μ_1, μ_2, κ_1, κ_2, κ_3)), if  has the probability density

    f_(ψ_1, ψ_2 |μ_1, μ_2, κ_1, κ_2, κ_3) 
       = C_s(κ_1, κ_2, κ_3)exp [ κ_1 cos(ψ_1 - μ_1) + κ_2 cos(ψ_2 - μ_2) + κ_3 sin(ψ_1 - μ_1) sin(ψ_2 - μ_2) ]

where κ_1, κ_2 ≥ 0, -∞ < κ_3 < ∞, μ_1, μ_2 ∈ [0, 2π) and the normalizing constant is given by

    C_s(κ_1, κ_2, κ_3)^-1 = 4 π^2 ∑_m = 0 ^∞2mm(κ_3^2/4 κ_1 κ_2)^m I_m(κ_1) I_m(κ_2).


In contrast,  is said to follow the (bivariate) von Mises cosine distribution with mean parameters μ_1, μ_2, concentration parameters κ_1, κ_2, and association parameter κ_3, denoted ∼(μ_1, μ_2, κ_1, κ_2, κ_3), if  has the probability density [<cit.> define the density with -κ_3 instead of κ_3 in the exponent. However, that makes the normalizing constant equal to C_c(κ_1, κ_2, -κ_3) in our current notation (i.e., in the form shown in (<ref>)) and not C_c(κ_1, κ_2, κ_3) as given in the paper. See Appendix <ref> for a proof.]

    f_(ψ_1, ψ_2 |μ_1, μ_2, κ_1, κ_2, κ_3) 
       =  C_c(κ_1, κ_2, κ_3) exp [ κ_1 cos(ψ_1 - μ_1) + κ_2 cos(ψ_2 - μ_2) + κ_3 cos(ψ_1 - μ_1 - ψ_2 + μ_2) ].

Here, similar to the sine model, κ_1, κ_2 ≥ 0, -∞ < κ_3 < ∞, μ_1, μ_2 ∈ [0, 2π) and the normalizing constant is given by

    C_c(κ_1, κ_2, κ_3)^-1 = 4 π^2 { I_0(κ_1) I_0(κ_2) I_0(κ_3) + 2 ∑_m = 0 ^∞ I_m(κ_1) I_m(κ_2) I_m(κ_3)  }.


From (<ref>) and (<ref>) it is easy to see that when κ_3 = 0, both the von Mises sine and cosine densities become products of univariate von Mises densities, implying independence between the two random coordinates. In addition, when κ_1 and κ_2 are also zero, both densities become uniform over [0, 2π)^2. <cit.> and <cit.> provide explicit forms for the marginal and conditional distributions in the sine and cosine models; the conditional distributions in both sine and cosine models are univariate von Mises, whereas the marginal distributions, although not von Mises, are symmetric around μ_1 and μ_2. 

One key difference between the bivariate wrapped normal model and the bivariate von Mises models is that  κ_3^2 is not required to be bounded above by κ_1 κ_2 in the latter, and thus can take any value in (-∞, ∞). Consequently, the densities can be bimodal; <cit.> show that the sine (cosine) joint density is unimodal if  κ_3^2 < κ_1 κ_2 (κ_3 ≥ - κ_1κ_2/(κ_1+κ_2)), and bimodal otherwise. This flexibility gives the two bivariate von Mises distributions  richer sets of possible contour plots and the ability to model a larger class of angular data.







Figures <ref> and <ref> display the surfaces of the von Mises sine and von Mises cosine densities respectively with μ_1 = μ_2 = π, κ_1 = κ_2 = 1 and various κ_3's. From Figure <ref>, it can be seen that the density is bimodal when κ_3 = ± 2 (or more generally for |κ_3| ≥ 1 when κ_1 = κ_2 = 1), and unimodal when |κ_3| < 1. It can also be seen that the density surface (or the contours) of a sine model with κ_3 = ξ is essentially a mirror image of that with κ_3 = -ξ, for any ξ∈ (-∞, ∞); see, e.g., the upper-left and the lower-right panels of Figure <ref>. Such is however, not the case for the cosine density, as depicted in Figure <ref>. The cosine density is bimodal when κ_3 is very negative (κ_3 ≤ -0.5 when  κ_1 = κ_2 = 1, see, e.g., the upper-left and upper-middle panels of Figure <ref>), and is unimodal otherwise. Moreover, flipping the sign of κ_3 does not yield density surfaces (or contours) that are mirror images of each other. 

An interesting feature of both sine and cosine densities is that they both approximate the regular bivariate normal density (on ^2) when the concentration parameters κ_1 and κ_2 are large, and the densities are unimodal (<cit.>, <cit.>). This property is analogous to the univariate von Mises distribution. A heuristic explanation of this result again follows from the fact that when the distributions are unimodal and κ_1, κ_2 are large, then ϕ_1 and ϕ_2 are highly concentrated around μ_1 and μ_2. This means ϕ_i - μ_i ≈ 0 so that sin(ϕ_i - μ_i) ≈ (ϕ_i - μ_i) and cos(ϕ_i - μ_i) ≈ 1 - (ϕ_i - μ_i)^2/2 for i=1,2.  




 §.§ Summary Measures for Univariate and Bivariate Angular Distributions
 
Circular summary measures are useful for describing various aspects of angular distributions. The circular mean or mean direction (see <cit.>) of an angular random variable ψ is defined as

    E_c (ψ) = arctan[ E(sinψ)/E(cosψ)]

and the circular variance of ψ is given by

    _c(ψ) = 1 - E[cos (ψ - E_c(ψ))].
  
Note that 0 ≤_c(ψ) ≤ 1.

When considering the joint distribution of paired angular random variables (ϕ, ψ), their association can be measured using circular correlation. Multiple parametric circular correlation coefficients have been proposed in the literature, and here we consider two of them.  Let μ_1 and μ_2 be the circular means of ψ_1 and ψ_2 respectively. Then the Jammalamadaka-Sarma (JS) circular correlation coefficient <cit.> is defined as

    ρ_(ψ_1, ψ_2) = E[sin (ψ_1 - μ_1) sin (ψ_2 - μ_2) ]/√(E[sin^2 (ψ_1 - μ_1) ] E[sin^2  (ψ_2 - μ_2)]).
 
Now let (ψ_1^(1), ψ_2^(1)) and (ψ_1^(2), ψ_2^(2)) be independent and identically distributed (IID) copies of (ψ_1, ψ_2). Then the Fisher-Lee (FL) circular correlation coefficient <cit.> is defined by

    ρ_(ψ_1, ψ_2) = E[sin(ψ_1^(1) - ψ_1^(2)) sin(ψ_2^(1) - ψ_2^(2)) ]/√(E[sin^2 (ψ_1^(1) - ψ_1^(2)) ] E[sin^2  (ψ_2^(1) - ψ_2^(2)) ]).
 
Both ρ_ and ρ_ have properties  similar to the ordinary correlation coefficient. In particular, ρ_, ρ_∈ [-1, 1] and they are equal to 1 (-1) under perfect positive (negative) toroidal-linear (T-linear) relationship <cit.>.

Note that all distributions considered in have circular mean(s) equal to the respective mean parameter(s). For the univariate models, the circular variances are just functions of the associated concentration parameter (see <cit.>). In particular, if ψ∼(μ, κ) then _c(ψ) = 1 - exp(-σ^2/2) with σ^2 = 1/κ, and for ψ∼(μ, κ), _c(ψ) = 1 - I_1(κ)/I_0(κ). For a bivariate wrapped normal model with  Σ = (Σ_ij) = Δ^-1, the marginal circular variance of the first coordinate is 1 - exp(-Σ_11/2), ρ_ = sinh(2Σ_12) /√(sinh(2Σ_11) sinh(2Σ_22)) and ρ_ = sinh(Σ_12) /
√(sinh(Σ_11) sinh(Σ_22)) <cit.>, where sinh denotes the hyperbolic sine function. For bivariate von Mises models (both sine and cosine forms), these expressions, provided in Appendix <ref>, are much more complicated, and involve infinite series of product of modified Bessel functions (see <cit.>). In we implement circular variances and correlation coefficients  for all the three bivariate models considered in this article. In addition, a function for calculating sample circular correlation coefficients is also provided, where the sample analogs of ρ_ and ρ_, along with two other non-parametric circular correlation coefficients are considered (see Section <ref>). 



 §.§ Mixture Models
 

Mixture models are convex combinations (mixtures) of two or more probability distributions, and provide a semi-parametric approach to modeling complex datasets with multiple noticeably distinct clusters. Mixture models of both univariate and multivariate (non-wrapped) normal distributions are well studied in the literature (e.g., see <cit.>), and implemented in many statistical packages, such as the R <cit.> packages mixtools <cit.>, mclust <cit.>, and Rmixmod <cit.>.  However, these are not applicable to mixture models for angular data.  This is a key motivation for our creation of BAMBI, which considers finite mixture models of univariate and bivariate angular distributions (the single function fit_angmix handles the fitting of all such models; see Section <ref> Category <ref>). 


Let K denote the number of components (where K is finite), { f(·|_j): j = 1, ⋯, K } denote the component densities (f can be univariate or bivariate) with _j denoting the parameter vector associated with the j-th component, and let = (p_1, ⋯, p_K)^T  denote the vector of mixing proportions (or weights) with p_j ≥ 0 and  ∑_j=1^K p_j = 1. Then the mixture density is defined as

    (·|; _1, ⋯, _K) = ∑_j=1^K p_j f(·|_j)

In practice, the number of components K necessary to fit the data is usually unknown,  and thus should be estimated on the basis of the data itself.  (See Section <ref> for a discussion on number of components estimation.) 


An important special case of the general mixture model (<ref>) is the mixture of product components, also called a conditional independence model. Here, one assumes each multivariate component density f(·|_j) to be a product of univariate densities; specifically for the bivariate angular models considered in BAMBI, this is achieved by letting κ_3=0 in each component. Note that a mixture of product components does not imply independence in the final mixture density. In fact, such a model can reasonably approximate a wide class of more complicated models, while being computationally less involved (see <cit.>);  however, one often needs a larger K compared to a general (non-product) mixture model to achieve similar results, thus offsetting some of the potential computational gains. In a  product component mixture can be fitted via fit_angmix by setting the argument cov.restrict = "ZERO" (see Section <ref> under Category <ref>).

It is also noteworthy to mention the aspect of bimodality of bivariate von Mises distributions in the context of mixture modeling. In practice, often each  component of a mixture model is used to represent one single (unimodal) cluster in data. However, as discussed in Section <ref>, both von Mises sine and cosine models can be bimodal depending on the values of the concentration and association parameters. When bimodality is present in some of the component specific densities, the final mixture model can be harder to interpret. To avoid this issue, it is possible to restrict the parameter spaces associated with the concentration and association parameters (by letting κ_3^2 < κ_1κ_2 in the sine model, and κ_3 ≥ -κ_1κ_2 / (κ_1+κ_2) in the cosine model) in these angular models to force unimodality in each component specific density.  Consequently, a larger K may be needed to achieve similar results, which increases model complexity. In  we provide an option of having only unimodal von Mises component densities. This is achieved by setting the logical argument unimodal.component = TRUE in fit_angmix (defaults to FALSE). See the  discussion in Section <ref>, Category <ref>.   



 §.§ Related Work and Motivation for BAMBI





  §.§.§ Literature


Several papers have addressed inferential problems relating to mixtures of bivariate angular distributions.  <cit.> consider the mixture of bivariate von Mises cosine distributions, and suggest an EM algorithm for frequentist estimation of the associated parameters. Their approach is used in <cit.> in the context of modeling protein backbone angles. In other work,  <cit.> consider a Bayesian non-parametric model involving an infinite mixture of von Mises sine distributions. In we focus on classical finite mixtures, providing a unifying framework for Bayesian estimation of all three bivariate angular models presented earlier.








  §.§.§ Software



To the best of our knowledge, no previous packages or libraries handle finite mixture modeling for univariate or bivariate angular data, whether in R or otherwise. In fact, the only available software (as of the time of writing this manuscript) that  has functionality for bivariate von Mises models  is the C++ library mocapy++ in the context of Dynamic Bayesian Networks  <cit.>. However, mocapy++ does not implement bivariate wrapped normal models.  




The overarching goal of BAMBI is to create a unified platform that implements descriptive and  inferential statistical tools required to analyze bivariate and univariate angular data. First, provides functions for  density evaluation, computation of various summary measures (such as circular mean, variance and correlation coefficient), and random data generation from bivariate and univariate angular models and their mixtures. Second, it has functions for fitting these models to real angular data using Bayesian methods. Third, it implements a number of post-processing steps required in any Bayesian statistical analysis. For example, visual and numerical assessment of the goodness of fits can be done using a number of native functions, as well as coda package functions, which are applicable on outputs (angmcmc objects) through a convenient  as.mcmc.list method.  Furthermore, has functions for model selection as well as random data generation and density evaluation from fitted models, which are useful in posterior predictive analyses.   









It is to be noted that while it is possible to use general-purpose MCMC samplers such as stan <cit.>, JAGS <cit.> and WinBUGS <cit.> for fitting the angular mixture models considered in BAMBI, there are important motivations for developing specialized implementations for these models.
First, special care needs to be taken while handling the normalizing constants in the von Mises sine and cosine densities, which contain infinite series of product of Bessel functions that can be numerically unstable for some ranges of parameter values; such cases are handled in BAMBI via (quasi) Monte Carlo approximations. Second, computations for Bayesian mixture modeling benefit from using a latent allocation structure, as done in BAMBI (see Section <ref>), which allows independent sampling of  the component specific parameters.  Such an approach cannot be used in stan due to the discreteness of the allocation (p. 79, Section 6.2 of the reference manual v2.18.0); consequently, the dimensionality of the parameter vector can hinder convergence of MCMC sampling for mixtures with many components. 
In contrast, JAGS/WinBUGS allows incorporation of discrete latent allocation; however, 


their sampling techniques 
do not make use of the gradient of the target (log) posterior density. As discussed in Section <ref>, Hamiltonian Monte Carlo uses the gradient  and hence is typically more efficient for sampling from intractable distributions.


Finally, the  analytic gradients necessary for efficient MCMC sampling in these models are built into BAMBI. 













 



























 §.§ Organization of the Paper


The remainder of this article is organized as follows. In Section <ref>, we  review Bayesian methods for fitting angular mixture models to data. In Section <ref> we describe the capabilities of BAMBI, by describing all functions and datasets available in BAMBI, and providing brief overviews on their usage. Following, in Section <ref> we illustrate angular mixture modeling on datasets included in BAMBI.  The paper concludes with a brief summary and possible directions for future development in Section 5. A derivation for the von Mises cosine model normalizing constant, formulas for circular variances and correlation coefficients in the von Mises sine and cosine models, analytic forms of gradients needed for efficient MCMC sampling (discussed in Section <ref>), and MCMC parameter traceplots associated with one of the examples considered in Section <ref> are  provided in the Appendices.




§ METHODS
 




 §.§ Overview

We adopt a Bayesian approach for fitting angular mixture models to data. Let ^⊤ = (_1, ⋯, _n) be the  data matrix (or  data vector in the univariate case)  with each _i being a bivariate vector of angles (or a univariate angle) [0, 2π)^2 (or in  [0, 2π)).  We are interested in fitting a mixture density of the form (<ref>) for a given number of components K. For example, in bivariate wrapped normal mixtures, the density for the j-th component is given by f_j  ≡ f_ (·|_j) =: f_, j, where _j^⊤ = (κ_1j, κ_2j, κ_3j, μ_1j ,μ_2j) denotes the vector of (model) parameters for the jth component, j = 1, ⋯, K, and the mixture density is given by _ = ∑_j=1^K p_j f_, j. For a specified K, our objective is to estimate the parameter vector ^⊤ = (^⊤, ^⊤),  which consists of the model parameters ^⊤ = (^⊤_1, ⋯, ^⊤_K) and the mixing proportions ^⊤ = (p_1, ⋯, p_K), based on . Often, K itself will also need to be estimated. In the following, we review some commonly used techniques in Bayesian mixture model fitting.




 §.§ Bayesian mixture modeling
 

Under a Bayesian framework a prior distribution must be specified for the parameter vector, which can be non-informative (or diffuse) if a priori information is unavailable. Let π(, ) denote the joint prior density for . Often the prior distributions of  and  are assumed to be independent so that (with a slight abuse of notation; here π(y) stands for the appropriate prior density of the random variable y) π(, ) = π() π(). Moreover, parameters from different components are often assumed to be independent, so that π() = ∏_j=1^Kπ(_j). Let L(|, ) = ∏_i = 1^n (_i |, ) denote the likelihood function of the data. Then the posterior density of  given the data is

    π(, |) ∝ L(|, )  π() ∏_j=1^Kπ(_j),

which is the basis for Bayesian inference on . It is to be noted that the prior densities π(_j)'s all need to be proper in order to ensure that the posterior density π(, |) is proper <cit.>.   Specific comments about the choice of priors used in the current setting are provided in Section <ref>. Note that the associated posterior mean, median or mode, commonly used as point estimates of the parameters, are not available  in closed form for our distributions of interest.  Additionally, π(, |) is intractable for directly simulating IID samples, and thus some kind of Markov Chain Monte Carlo (MCMC) technique is used in practice as an alternative.  Starting from some initial point, an MCMC algorithm generates a Markov chain which has the target posterior density π(, |) as the invariant distribution. Various summary measures of the posterior distributions – such as mean, mode (known as the maximum a posteriori or MAP parameter value), and quantiles –  can then be approximated based on the MCMC realizations. In practice, the MCMC algorithm must be run long enough for the Markov chain to converge, so that the realizations approximately follow the target posterior distribution.  For this purpose the chain is given a burn-in period, where the initial iterations are discarded.





 

In the function fit_angmix fits a Bayesian angular mixture model with a specified number of components, and the function fit_incremental_angmix fits angular mixtures with incremental number of components to determine an optimum number of components. In the following we briefly review the MCMC generation techniques Gibbs sampler (GS), Metropolis Hastings and Hamiltonian Monte Carlo (HMC), and describe how they are used for sampling from the posterior distributions of model parameters and mixing proportions in these two functions.




 §.§ Gibbs sampler (GS)
 

The Gibbs sampler (GS) <cit.> breaks the Markov chain updates for the parameter vector into blocks.  For example, when = (_1, _2) the GS generates the N-th state of the Markov Chain (_1^(N), _2^(N)) from the previous state (_1^(N-1), _2^(N-1)) with the steps



  * Generate _1^(N) from π(_1 |_2^(N-1),).

  * Generate _2^(N) from π(_2 |_1^(N),).

The GS is most effective when it is easy to sample from the (full) conditional posterior densities π(_1 |_2,) and π(_2 |_1,).  Note that when _1 and _2 are vectors, this is sometimes called the blocked Gibbs sampler.

For mixture models, an efficient Gibbs sampling step for the mixing proportions  (when K > 1) can be obtained by adopting a so-called Data Augmentation scheme, where one introduces (“augments”) unobserved data to make the conditional distributions simpler <cit.>. Here, we introduce (hidden) component indicators _i^⊤ = (ζ_i1, ⋯, ζ_iK) corresponding to each observation _i where ζ_ij is 1 if the ith observation comes from the jth component, and 0 otherwise, for i = 1, ⋯,n and j=1,⋯,K. Thus, given ζ_ij = 1, the density of _i is simply f(_i |_j), i.e., the density in the j-th component evaluated at _i. Moreover,  n_j := ∑_i=1^n ζ_ij is the total number of observations coming  from this density. It is customary to assume a Dirichlet() prior for , where ^⊤ = (α_1, ⋯, α_K) with α_j > 0 for all j, so that π() ∝∏_j=1^K p_j^α_j-1. Note that α_j = 1 for all j represents the uniform prior. Let Z^⊤ = (_1^⊤, ⋯, _n^⊤) and let ^(N-1), ^(N-1) and Z^(N-1) be the (N-1)th MCMC realizations of ,  and Z respectively. Then the Nth realization of  (and Z) are obtained as follows:


  * For i = 1, ⋯,n, generate _i^(N) from Multinomial(1; _i1^(N-1), ⋯, _iK^(N-1)) independently, and define n_j^(N) := ∑_i=1^n ζ_ij ^(N), where

    _ij^(N-1) = p_j^(N-1) f(_i |_j^(N-1))/∑_h=1^K p_h^(N-1) f(_i |_h^(N-1))

are the posterior membership probabilities.


  * Generate ^(N) from Dirichlet(α_1+n_1^(N), ⋯, α_K+n_K^(N)).



Thus when K > 1, the latent allocation _i's generated during the Gibbs sampling step for  leads to simplifications that reduce the computational burden substantially.  Note that, conditional on _i's, all _i's have independent single component densities f(·|_j_i), with j_i being the non-zero position of _i. Thus, given _i's, all _j's are independent with only data points coming from component j contributing to the respective likelihoods. Consequently _j's can be sampled independently from their (component-specific) conditional posterior densities.

To complete the GS scheme for the mixture model, it remains to sample _j's from π(_j | Z, , ). As these distributions are still intractable for direct IID simulation, we use a Markov chain simulation technique for sampling, and then combine this step with the GS updates for  and Z. In the following we describe two such Markov chain simulation techniques, and how they are used in BAMBI.




 §.§ Metropolis-Hastings Algorithm
 

The Metropolis-Hastings algorithm  <cit.>  is simple and widely-used for Markov chain simulation. 
Formally, let x be the current state of a Markov chain Φ with stationary density q. Let (·| x) be a proposal density defined on the state space of Φ that is easy to sample from.  Then the next state x' of the Markov chain Φ is obtained as follows:

	
  * Generate x^* from _x.
	
  * Define r(x^*, x) = min{1, q(x^*)/q(x) (x | x^*)/(x^* | x)}, and define the next state x' equal to x^* with probability r(x^*, x) and equal to x with probability 1-r(x^*, x).


The random walk variant of Metropolis-Hastings (RWMH) uses a proposal density (·| x) that is symmetric about x; e.g., by taking (·| x) to be the density of Y_x = x + Y_0,  where Y_0 is a normal random variable with mean zero. Under RWMH, (x | x^*) = (x^* | x), and hence r(x^*, x) = min{1, q(x^*)/q(x)}, thus simplifying computations. In BAMBI, RWMH is implemented with independent normal proposals.

Note that the variance of the density (·| x) strongly affects the acceptance probabilities r(x^*, x). Convergence of the Markov chain will be slow if the variance of (·| x) is too large or too small. <cit.> suggest maintaining an acceptance rate of 20-30% as a general rule-of-thumb. In we provide an auto-tuning feature that implements adaptive tuning during the burn-in period.  Briefly, the acceptance rate and scale of the sampled parameters  are monitored at regular intervals, and the proposal variances are adjusted accordingly (see the documentation of fit_angmix for details).  We limit adaptation to the burn-in period, so that the desired properties of the final MCMC samples are retained.






 §.§ Hamiltonian/Hybrid Monte Carlo (HMC)
 
Simple RWMH can become quite inefficient in multi-dimensional problems.  A powerful alternative to RWMH when the gradient of the posterior density has an analytical form is  Hamiltonian (also called Hybrid) Monte Carlo (HMC) <cit.>.  HMC makes use of the gradient of the log posterior density and an auxiliary random variable, and incorporates tools from molecular dynamics to furnish proposal states coming from high posterior density regions. This allows a much faster exploration of the state space than a RWMH scheme. A gentle and detailed introduction to HMC with applications to statistical problems can be found in <cit.>.  Briefly, in HMC first an auxiliary random variable  called momentum is considered along with the variable of interest (vector of model parameters  in our case), which is classically called the position in physical problems, denoted by .[In classical HMC literature, the auxiliary variable is denoted by ; however, we will keep that notation for mixing proportions.] Furthermore, two energy functions U()  and K() are introduced, followed by a Hamiltonian function H(, ) which is usually the sum of those two energies, i.e., H(, ) = U() + K(). U(), called the  potential energy, is defined as the negative log posterior density of  (plus any fixed constant) in MCMC applications, and  K(), called the kinetic energy, is usually defined as K() = ^⊤ M^-1 for some fixed positive definite matrix M.  This form for K() corresponds to the negative log  density (plus a constant) of the zero-mean normal distribution with variance matrix M. In practice, M is typically taken to be diagonal, often the identity matrix (as used in BAMBI), or a scalar multiple of the identity matrix. Let ∇ U() denote the gradient vector of U() with respect to . Further, let ϵ > 0 be a small real number, called the step-size, and L≥ 2, a positive integer, called the number of  leapfrog steps.  Then one step of HMC that updates (via leapfrog method) the current state  to the next state ' can be described as follows:

	
  * Generate  from N(0, M) and let ^(0) = and ^(0) =  - (ϵ/2) ∇ U(^(0)).
	
  * For t = 1, ⋯, L define ^(t) = ^(t-1) + ϵ ^(t-1) and ^(t) = ^(t-1) - (ϵ/γ_l) ∇ U(^(t)), where γ_l = 1 for l = 1, ⋯, L-1 and γ_L = 2.
	
  * Let ^* = ^(L) and ^* = -^(L), and define β(^*,^*; ,) = min{ 1,  exp[ H(^*, ^*) - H(, ) ] }.
	
  * Finally, define the new state ' equal to ^* with probability β(^*,^*; ,), and equal to  with probability 1-β(^*,^*; ,).

Special care needs to be taken for the cases where the variables being sampled are constrained: for our angular models, μ_i's are angles in [0,2π), and the (raw) concentration parameters are positive. See  section 5.5.1.5 of <cit.>  for more details.

Since HMC approximates the dynamics by discretization, the step-size ϵ needs to be sufficiently small for the proposals to have a high acceptance rate. However,  if ϵ is too small, convergence of the Markov chain will be slow. Thus, ϵ requires tuning to obtain a reasonable acceptance rate (∼40-90%, with 65% being optimal, as suggested by <cit.>).   In we provide an auto-tune feature for ϵ similar to the one for the proposal standard deviation in RWMH (see Section <ref>), which adaptively tunes ϵ during burn-in to ensure a reasonable acceptance rate (60-90% by default).

Care is required for choosing the number of leapfrog steps L, since a L that is too large or too small can lead to poor convergence.  While setting an appropriate L can be challenging for high dimensional parameter vectors, here the independence of components π(_j | Z, , ) means that only two (for univariate models) or five (for bivariate models) parameters need to be sampled at a time.  Thus, the default L=10 used in BAMBI, which works well empirically, suffices for mixtures with any number of components. As suggested in <cit.>, randomly choosing ϵ and L from some fairly small interval at the beginning of every HMC step may improve convergence of the chain. In ϵ is by default randomly chosen at each iteration from an interval of the form (ϵ_0(1-δ), ϵ_0(1+δ)) for a fixed ϵ_0 > 0 (can be auto-tuned in BAMBI) and a given δ∈ (0, 1), while L is kept fixed. However, these settings can be changed; in particular, L can also be randomly chosen  from the set of integers contained in an interval (L_0/exp(A), L_0 exp(A)) for some given L_0 > 0 and A > 0, or both ϵ and L can be specified to be non-random. See the documentation of fit_angmix for more details. 

When properly tuned, HMC can achieve faster convergence and better exploration of the target density than RWMH, for a similar computational cost. Note that the computational cost for each HMC iteration is higher due to L gradient evaluations, however, HMC usually requires fewer iterations to reach stationarity and successive samples have lower autocorrelation.  Hence, HMC is our recommended sampling approach in BAMBI.  HMC, while powerful, does not solve all the challenges associated with MCMC sampling algorithms; in particular, both RWMH and HMC can get trapped in local modes. One possible remedy is to use multiple independent chains, see Section <ref>.

By default, uses HMC to sample .  All angular densities considered here,  both univariate and bivariate, admit analytic gradients for efficient programming implementation. Expressions for the conditional log posterior density and its gradients are provided in the following section.




 §.§ Using RMWH or HMC for angular mixture models

Consider the mixture model (<ref>) with density f(·|_j) for the j-th component, j=1,⋯,K. It follows that given the component indicators Z, information on  is superfluous, and the complete-data (i.e., given  and Z) likelihood for = (_1, ⋯, _K) is given by:

    likelihood(| Z, ) ∝∏_i = 1^n ∏_j=1^K  f(_i  |_j)^ζ_ij.

Recall that the joint prior density of  is ∏_j=1^K π(_j). Hence, the complete-data posterior density of  is given by:

    π(| Z, ) ∝{∏_i = 1^n ∏_j=1^K  f(_i  |_j)^ζ_ij}∏_j=1^K π(_j).

Therefore, by taking the logarithm, the complete-data log posterior density (LPD) for = (_1, ⋯, _K) given the component indicators Z is obtained as

    _complete-data() := logπ(| Z, )    = C + ∑_i=1^n ∑_j=1^Kζ_ijlog f(_i  |_j) + ∑_j=1^Klogπ(_j) 
       = C + ∑_j=1^K{∑_i=1^n ζ_ijlog f(_i  |_j) + logπ(_j) }
       = C + ∑_j=1^K{∑_i ζ_ij = 1log f(_i  |_j) + logπ(_j) }

where C is a constant free of . The above expression shows that conditional on Z, _j's are independent, and that the complete-data log posterior density of _j is of the form

    _j(_j) = C_j + ∑_i ζ_ij = 1log f(_i  |_j) + logπ(_j)
  
where C_j's are  constants (free of ).  Given the current GS draw of Z, samples from the conditional posterior density _j in (<ref>) can therefore be drawn independently for all j=1,⋯, K. For each j ≥ 1,  we let _j play the role of the target density q (see Section <ref>) in RWMH, or let -_j play the role of  the potential energy U (see Section <ref>) in HMC; the gradient of U with respect to _j, ∇ U, is therefore the negative of the gradient ∇_j. From (<ref>), it follows that


    ∇_j()    =  ( ∑_i ζ_ij = 1∂log f(_i  |_j)/∂_j) + ∂logπ(_j)/∂_j
       =  (∑_i ζ_ij = 11/f(_i  |_j)·∂ f(_i  |_j)/∂_j) + ∂logπ(_j)/∂_j
	 

For the von Mises distributions (both univariate and bivariate), form (<ref>) is easier to work with, whereas form (<ref>) is more useful for the wrapped normal distributions. Full analytic expressions for all model specific gradients are provided in Appendix <ref>.

Note that parameters with a non-negative support are often sampled more efficiently on the log scale; we use this strategy for sampling the concentration parameters κ (in univariate models) and κ_1, κ_2 (in bivariate models). 










 §.§ Choice of priors
 
Selection of prior constitutes an important step in Bayesian analyses, as they play a key role in the final inference. This is comparatively more standard for the component-specific model parameters . As discussed, proper prior distributions for the model parameters are required to ensure posterior propriety.  For the mean parameters μ (in univariate models) and μ_1, μ_2 (in bivariate models), their prior distributions can be taken to be a member of the same family of the distribution which are being used in the mixture model (e.g., von Mises sine prior for (μ_1, μ_2) in a von Mises sine mixture model) to aid conjugacy. <cit.> use this conjugate prior for the mean parameters in their von Mises sine (infinite) mixture model. Note that conjugacy for the mean parameter is not achievable except in trivial cases in the wrapped normal distributions (both univariate and bivariate).  In we set a uniform prior over [0, 2π) (if univariate) or [0, 2π)^2 (if bivariate)  for the mean parameter(s), which can be viewed as a special case of the von Mises and wrapped normal distributions (see Sections <ref> and <ref>). Conjugacy is also possible for the concentration and association parameters, e.g., <cit.> consider such a family for von Mises sine model. However, that approach does not aid sampling, as the resulting unnormalized densities involve infinite sums of products of modified Bessel functions.  
As a simple alternative, we suggest using independent normal distributions with zero mean as the prior for the association parameter κ_3, as well as for the log of the concentration parameters κ, κ_1, and  κ_2 (i.e.,  the prior for concentration parameters are log normal). These prior distributions can be made informative or diffuse through appropriate choices of the variance hyper-parameter. Priors are assigned independently to each parameter, and truncation is performed to reflect any specified constraints in the model (such as κ_3^2 < κ_1 κ_2 in a bivariate wrapped normal model, and a von Mises sine model with unimodal density). 

Care is required in the selection of prior for the mixing proportions , as an ill-chosen prior may result in very poor fits. This is particularly true when K is too large (i.e., the mixture is overfitted). Note that overfitting is a necessary step when the true number of components is unknown and needs to be estimated, see Section <ref> for more details.  
It is customary to assume a Dirichlet() prior for , where ^⊤ = (α_1, ⋯, α_K) with α_j > 0, often with the special case  α_j = α_0 for all j. When the mixture is overfitted, the asymptotic results in <cit.> show that α_j's strongly influence how the spurious mixture components are handled by the limiting posterior density. In particular, if max_j=1,⋯, Kα_j < d /2, where d = _j, then the spurious components vanish asymptotically. On the other hand, if min_j=1,⋯, Kα_j > d /2, then the spurious components asymptotically get  superimposed on some of the existing components with positive mixing proportions <cit.>. The subsequent estimation of K depends on which way the overfitting is handled by the posterior density (see Section <ref>); thus α_j's all need to be appropriately either small or large <cit.>.  A uniform prior with α_j=α_0 = 1 is a rather poor choice in this regard. In estimation of K is done assuming the use of α_j > d/2 for all j in conjunction with a model selection criterion;  our default is α_j = α_0 =  (r +r(r+1)/2)/2 + 3 as used in <cit.>, where r denotes the dimension of the data, i.e., r is 1 or 2 according as whether the model is univariate or bivariate (and consequently, all α_j's are either 4 or 5.5).






 §.§ Estimating the number of components K from data
 


Suppose the data were generated from a mixture of K_ (non-empty, non-identical) components. In practice, K_ will not be known, and therefore mixture modeling requires estimating the appropriate number of components from the data.

In the Bayesian setting, the estimation of K_ requires an overfitted mixture model, i.e., one that has spurious or superfluous components. There are two ways of introducing superfluous components to overfit a mixture model, and the subsequent estimation of K_ should reflect which way is taken. First, the superfluous components can be arbitrarily introduced at regions with no data points (“leave some groups empty”), and assigned zero mixing proportions. Then, the number of non-empty components in the fitted mixture provides a good estimate of K_. Second, the spurious components can be superimposed on some of the existing components (“let two component-specific parameters be identical”), and assigned positive mixing proportions. Here, the total number of components in the fitted mixture, after accounting for model complexity (via some model selection criterion), provides a reasonable estimate for K_.  Note that the prior distribution of the mixing proportion  affects the way overfitting is handled by the posterior, and hence the associated prior hyper-parameters need to be wisely chosen (see Section <ref>). A detailed discussion on the estimation of the number of components can be found in <cit.>. 

In we assume that the superfluous components are introduced in the second (“let  two component-specific parameters be identical”) way. Consequently, K_ is estimated by first incrementally fitting the data with one additional component (starting from K=1), until a model with K+1 component fails to improve upon the previous fit with K component (as determined through a model selection criterion); that value of K is then used as an estimate of K_. There exist multiple model selection criteria in the literature; we review six such criteria implemented in BAMBI and comment on their applicability in MCMC simulations. In the following, = (, ) denotes the entire parameter vector, and =(y_1, …, y_n) is the vector/matrix of n independent  observations.




    	
  * Watanabe-Akaike Information Criterion (WAIC) (<cit.>). Given the dataset y_1, ⋯, y_n, the Markov chain realizations {_1, ⋯, _N } of the parameter vector, and the pointwise densities {p(y_i |_s): i = 1, ⋯, n;   s = 1, ⋯, N }, define the computed log pointwise posterior predictive density
    	
    LPPD = ∑_i=1^n log(1/N∑_s=1^N p(y_i |_s) ).

    	Then WAIC is defined as
    	
    WAIC = LPPD - p_W

    	where p_W is a correction term to adjust for effective number of parameters. Two forms for the adjustment terms are proposed in the literature, both being approximations based on Bayesian cross validation. In the first approach, (computed) p_W is defined as
    	
    p_W = 2 ∑_i=1^n [ log(1/N∑_s=1^N p(y_i |_s) ) -1/N∑_s=1^N log p(y_i |_s)  ]

    	whereas, in the second approach, (computed) p_W is defined by p_W = ∑_i=1^n var log  p(y_i |), where for i = 1, ⋯, n, var log  p(y_i |) denotes the estimated variance of p(y_i |) based on the realizations _1, ⋯, _N.
    	
    	
  * Leave One Out Cross Validation Information Criterion (LOOIC) <cit.>. Under the same set-up as WAIC, the LOOIC is defined as
    	
    LOOIC = ∑_i=1^n log( ∑_s=1^N w_i^s p(y_i |_s)/∑_s=1^N w_i^s)

    	where for each s = 1, ⋯, N, w^s =(w_1^s, ⋯, w_n^s) is a vector of importance sampling weights, typically calculated via the Pareto smoothed importance sampling method <cit.>. Because of the importance sampling weights,  LOOIC can be more stable in practice than WAIC. See <cit.> for a gentle and thorough introduction to both WAIC and LOOIC, including applications and case studies. 
    	
    	Because both WAIC and LOOIC are based on the mixture likelihood and do not explicitly depend on the sampled model parameters (thus, remain unaffected by the presence of multiple permutation and non-permutation modes), in BAMBI we recommend using either of these two criteria for selecting the number of mixture components.  Both WAIC and LOOIC are made available in via their implementations in the R package loo <cit.>, which also provides a compare() function for comparing WAICs/LOOICs based on estimated difference in expected log predictive density (elpd). In BAMBI, during an incremental model fitting via fit_incremental_angmix with crit = 'WAIC' or crit = 'LOOIC', a test of hypothesis H_0K: elpd for the fitted model with  K components ≥ elpd for the fitted model with K + 1 components, is performed at every K ≥ 1. The test-statistic used for the test is an approximate z-score based on the normalized estimated elpd difference between the two models (obtained from compare(), which provides estimated elpd difference along with its standard error). The incremental fitting stops if  H_0K cannot be rejected (at level alpha, defaults to 0.05) for some K ≥ 1; this K is then regarded as the optimum number of components. 

	
  * Marginal Likelihood (ML). Marginal likelihood is arguably the most natural and intuitive model selection criterion that is used in the Bayesian paradigm. As the name suggests, marginal likelihood is the likelihood obtained by integrating out the parameters from the joint density of the data and parameters, and provides a natural way of measuring the “marginal” effect of data. In the context of Bayesian model selection, marginal likelihood provides a way of selecting an optimum model in that the model with largest marginal likelihood provides the best fit.  Given the likelihood L(|) and prior density π(), marginal likelihood is the constant (dependent only on the data):
	
    m() = ∫_ℰ L(|) π()   d

	where ℰ denotes the support of the parameter vector . Note that m() is the reciprocal of the normalizing constant required to define the posterior density. Evaluation of the marginal likelihood m() in practice however is typically challenging, as it tends to be a high-dimensional intractable integral (as in our case). Multiple estimation techniques based on samples from the posterior density π(|) have been proposed in the literature; in we implement bridge sampling <cit.>.


	Briefly, the key idea is to first write m() as
	
    m() = ∫_ℰ h() L(|) π() g() d/∫_ℰ h() g()  π(|)  d
    	= E_g[h() L(|) π()] /E_π(·|)[h() g() ]

	where g is a density, called the proposal density, and h is a function, called the bridge function. Then one approximates the above ratio by
	
    m̂() = 1/n_2∑_j=1^n_2 h(_j^*) L(|_j^*) π(_j^*)/1/n_1∑_j=1^n_1h(_j^†) g(_j^†)

	where _1^†, ⋯, _n_1^† are MCMC samples from the posterior density π(·|) and _1^*, ⋯, ^*_n_2 are samples from the proposal density g. Note that h and g play crucial roles in the estimation of m(), and must be optimally chosen for accurate results. See <cit.> for a gentle and detailed tutorial on bridge sampling. In BAMBI, marginal likelihood can be used to select the optimal number of mixture component in a fit_incremental_angmix run, by specifying crit = 'LOGML'. This will ensure computation of the log marginal likelihood via bridge sampling for every mixture model during the incremental run, and the model attaining the first minimum negative log marginal likelihood will be treated as the optimum model. 
	
	It should however be noted  that for mixture models, optimal selection of h and g is typically difficult  due to the multi-modality of the posterior density; see <cit.>, for a review of some of the available methods. In  BAMBI, marginal likelihood is computed by leveraging the function bridge_sampler from the R package bridgesampling <cit.>, and the authors of bridgesampling warn against the use of  bridge_sampler in mixture models <cit.>. As such, use of this method in is not recommended, even though the functionality is available. 

	
	
	
  * Akaike Information Criterion (AIC) <cit.>.  Let L̂ be the maximum value of the likelihood function for the model and  let m be the number of estimated parameters in the model. Then AIC is defined as
	
    AIC = -2 logL̂ + 2m.

	
	
  * Bayesian Information Criterion (BIC) <cit.>. Under the same setup, if n denotes the number of data points, BIC is defined by
	
    BIC = - 2 logL̂ + m log(n).

	Observe that both AIC and BIC depend on the maximum value L̂, which, in general, is not directly available in MCMC simulations. A possibly suboptimal estimate of the global maximum is given by the maximum value of the likelihood function computed at the MCMC samples. 
	
	During model selection, the model with minimum AIC (or BIC) can be treated as the optimal model. In AIC/BIC can be used for selecting optimum number of components in an fit_incremental_angmix run by specifying crit = 'AIC' or crit = 'BIC'. This will ensure computation of AIC/BIC of every mixture model fitted during the incremental fitting; the model attaining the first minimum AIC/BIC will be treated as the optimum model.
	
	It should to be noted, however, that AIC and BIC are both based on asymptotic normality results that do not hold for  mixture models with multiple modes, and hence their use in selecting the number of mixture components may lead to inconsistent results. Thus, though implemented, using AIC or BIC is not recommended in BAMBI.  
	
	
  * Deviance Information Criterion (DIC) <cit.>. DIC is another model selection criterion, which, similar to AIC and BIC, is based on an asymptotic result for large samples. Let   D( ) = -2log p(| )  denote the deviance, where  denotes the vector of all parameters in the model and p(|) denotes the likelihood. Let {_1, ⋯, _N } denote the MCMC realizations of the parameters. Define (estimated) effective number of parameters p_D by  p_D = D̅() - D(), where D̅() = N^-1∑_s=1^N D(_s)  and = N^-1∑_s=1^N _s. Another commonly used form for p_D is given by p_D = var  D() / 2, <cit.> where var  D() denotes the estimated variance of D() based on _1, ⋯, _N. Then DIC is defined as
	
    DIC = p_D + D̅() = D () + 2 p_D.

	In AIC/BIC can be used for selecting optimum number of components in an fit_incremental_angmix run by specifying crit = 'DIC'. This will ensure computation of DIC of every mixture model fitted during the incremental fitting; the model attaining the first minimum DIC will be treated as the optimum model.
	
	It should be noted that use of DIC can be unstable in practice. For example,  if the first form of p_D, i.e., p_D = D̅() - D()) is used, DIC becomes heavily dependent on the plug-in estimator . However, in Bayesian mixture modeling the posterior mean is not always a suitable plug-in estimator for the parameter vector as it may lie between different modes of the posterior density <cit.>. The problem is exacerbated by the presence of label switching in the MCMC samples (see Section <ref>). Moreover, depending on how the information on latent component indicators are handled, multiple versions of DIC can be constructed here. <cit.> consider no less than eight variants, but are unable to recommend any of them for practical use. Likewise we caution against the use of DIC, although the functionality is available in BAMBI. 
	
	
 






 §.§ Label switching
 
Label switching is a fundamental aspect of Bayesian mixture modeling that requires proper care. Briefly, when exchangeable priors π(_j)'s are placed on the model parameters _j's, the resulting posterior distribution becomes invariant to permutation in the component label j's. As a result, the posterior density  consists of symmetric or permutation modes that are identical up to permutation of component labels. A well mixing  MCMC algorithm will explore these permutation modes,  causing the component labels to switch over the course of an MCMC simulation. This phenomenon is called label switching, and is required in MCMC-based Bayesian mixture modeling to justify convergence. A fundamental limitation of 
MCMC-based Bayesian mixture modeling is that the chains may become trapped at local modes, rather than fully exploring the symmetric modes. One possible remedy is to run multiple independent chains to improve exploration of the posterior.  Alternatively, one may embed a deliberate random relabeling step into the sampler, i.e., adopt a so-called permutation sampling <cit.> scheme: after each draw of the  random allocation, components are relabeled according to a random permutation of 1, ⋯, K. This is, in fact, a specific example of a sandwich algorithm <cit.>, where a computationally inexpensive step (of drawing a random permutation of {1,…, K}, and then relabeling the components according to that random permutation) is sandwiched in between the two steps (drawing the allocation vector, and drawing the component-specific parameters) of a Data Augmentation algorithm. Sandwich algorithms often converge faster  than the original Data Augmentation algorithm; in the case of permutation sampling, the chain is forced to visit the permutation modes (and potentially the non-permutation modes) more frequently. These potential improvements in convergence  would not be achieved by simply randomly switching labels of the MCMC samples post-hoc. In permutation sampling, inclusion of the random label switching step results in a modified MCMC algorithm that  is theoretically proven to be at least as good as the original Data Augmentation algorithm in terms of convergence rates <cit.>. 

However, care must be taken if RWMH or HMC updates for the component specific parameters are adaptively tuned according to the scales and variabilities of the sampled model parameters (to do so properly requires keeping track of each component label). In permutation sampling can be done after burn-in, by setting perm_sampling = TRUE (defaults to FALSE) in a fit_angmix call. 



Although label switching is required for MCMC convergence in Bayesian mixture modeling, its presence in MCMC samples makes inference on the different components via posterior means or quantiles challenging (note that MAP estimation is not affected). A number of techniques have been proposed to handle this problem; see, e.g., <cit.> and <cit.>. The available methods either need to be applied during MCMC sampling (on-line) or after simulating the entire chain (post processing).

Several post-processing techniques that undo label switching are implemented in the useful R package label.switching <cit.>, and in BAMBI, we provide a wrapper called fix_label for the main label.switching function from that package. All the methods available in label.switching are appropriately implemented in fix_label, which takes an angmcmc object (see Section <ref>) as input, and may require additional user inputs depending on the method used. The Kullback-Leibler divergence based method by <cit.> (method = 'STEPHENS') is used by default if the permutation sampling is performed during original MCMC run; otherwise, the default method is the data-based algorithm of <cit.> (method = 'DATA-BASED'); neither requires any additional input other than an angmcmc object.




 §.§ Initialization of the parameters, and the use of multiple chains for faster convergence



MCMC algorithms can converge faster if the initial values are chosen well.


The function fit_angmix, when called without supplying starting parameter values, will automatically initialize the latent allocation to the mixture components. The default (and recommended) option is initialization via a k-means algorithm:  toroidal angle pairs are first projected onto the surface of a unit sphere, and then Cartesian coordinates of the projected spherical points are clustered.  Random initial allocation is also provided as an option, but is not recommended as it may lead to slow convergence.  Once an initial allocation is obtained, component specific parameters are  estimated via method of moments (see <cit.> for more details on these estimators), and the mixing proportions are estimated by the sample proportions.  

When explicit starting values of the model parameters and mixing proportions are provided to fit_angmix,  no initial allocation is necessary. This is particularly useful for estimating the number of components, when mixture models are being fitted incrementally (e.g., via fit_incremental_angmix).  Under incremental model fitting the parameters of a K+1 component mixture can be  initialized directly from the parameter estimates from a K component mixture; the  extra component is simply taken as a “copy” of an existing component (preferably the one with the largest mixing proportion), and the associated mixing proportion is distributed equally between the two identical components. This method is expected to work well when the posterior density handles overfitting in the “let two component-specific parameters be identical” way  (see Sections <ref> and <ref>), which is the approach taken in BAMBI. As such, this is the default method of initializing parameters in fit_incremental_angmix  when K > 2; however k-means allocation followed by moment estimation can also be used, by setting prev_par = FALSE. 


Finally, we note that even with good initial values, MCMC samplers can still get trapped in local modes for a large number of iterations,  rather than fully exploring the posterior density. 
One possible remedy is to run multiple independent chains to improve exploration of the posterior, which is implemented in BAMBI. The argument n.chains (set to 3 by default) specifies the number of independent chains to run in fit_angmix. These chains can be run in parallel, see Section <ref> Category IV for more details.      




§ BAMBI PACKAGE
 
This section overviews the functionalities of BAMBI. At the core of the package is the angmcmc object, which is created when a model fitting function is used. In the following we first describe the angmcmc objects, then describe the datasets included in BAMBI, and finally discuss the functions available in BAMBI and comment on their usability. However, this is not an exhaustive manual; all functions in BAMBI include R documentations, which should serve as the definitive resources.



 §.§ "angmcmc" objects
 

angmcmc objects are classified lists belonging to the S3 class angmcmc that are created when the function fit_angmix is used. An angmcmc object contains a number of elements, including the dataset and its dimension (i.e., univariate or bivariate), the model being fitted, the tuning parameters used,  MCMC samples of the parameter vector, and at each iteration the (hidden) component indicators for data points,  log-likelihood and log posterior density values (up to additive constants). When printed, an angmcmc object returns a brief summary of the function arguments used and the acceptance rate of the proposal states (in HMC and RWMH). An angmcmc object can be used as an argument for the diagnostic and post-processing functions available in BAMBI for making further inferences.
















































































 §.§ Datasets


contains two illustrative datasets, namely wind (univariate) and tim8 (bivariate), each measured in the radian scale [0, 2π).


	
  wind The wind dataset consists of 239 observations on wind direction (originally measured in 10s of degrees, and then converted into radians) measured at Saturna Island, British Columbia, Canada during October 1-10, 2016 (obtained from Environment Canada website). There was a severe storm during October 4-7 in Saturna Island, which caused significant fluctuations in wind direction.
	
    
  tim8 The dataset tim8  consists of 490 pairs of backbone dihedral angles (ϕ, ψ) for Triose Phosphate Isomerase (8TIM). The three dimensional structure of  8TIM is available from the protein data bank (PDB).  The  protein is an example of a TIM barrel, a common type of protein fold exhibiting alternating
    α-helices and β-sheets. The backbone angles for this protein were obtained by using the DSSP software (<cit.>) on the PDB file for 8TIM, and then converted into radians. 			





 §.§ Functions
 
In BAMBI, all five models described in Section <ref>, namely the univariate von Mises (vm), univariate wrapped normal (wnorm), bivariate von Mises sine (vmsin), bivariate von Mises cosine (vmcos) and bivariate wrapped normal (wnorm2), and their (within same model) mixtures  are implemented. The functions in can be classified into six major categories:




  
Functions for evaluating density and generating random samples from an angular distribution.  
  The functions dvm, dwnorm, dvmsin, dvmcos and dwnorm2 evaluate the density and the functions rvm, rwnorm, rvmsin, rvmcos and rwnorm2  generate random samples from the models vm, wnorm, vmsin, vmcos and wnorm2 respectively. The parameters of the models are specified as arguments; otherwise, default values (zero means, unit concentrations, and zero association) are used.
  
Density evaluations require computation of the normalizing constants, which for the vmcos model requires proper care, especially when κ_1, κ_2 or |κ_3| is large. This is because the analytic expression involves infinite (alternating if κ_3 < 0) series of product of modified Bessel functions, which become numerically unstable when these parameters are large. As such, when κ_3 < -5 or max(κ_1, κ_2, |κ_3|) > 50, the reciprocal of the integral for the normalizing constant is evaluated numerically using a (quasi) Monte Carlo method. By default, n_qrnd = 10^4 pairs of Sobol numbers are used for this purpose; however, n_qrnd, or a two-column matrix qrnd containing random/quasi-random numbers between 0 and 1 can be supplied for this approximation[The user may perform a Monte Carlo approximation for the normalizing constant even when numerical evaluation of the analytic formula is stable, by changing force_approx_const to TRUE from its default value FALSE.]. For vmsin model, evaluation of the constant via its analytic form is much more stable, as the associated infinite series consists only of non-negative terms.  For univariate and bivariate wrapped normal models, the default absolute integer displacement for approximating the Wrapped Normal sum is 3, which can be changed to any value in {1, 2, 3, 4, 5}, through the argument int.displ. Note that int.displ regulates how many terms would be used to approximate the infinite sum present in the univariate and bivariate wrapped normal densities in (<ref>) and (<ref>). For example, int.displ = M implies that the infinite sum in the univariate wrapped normal density will be approximated by a finite sum of 2M +1 values, with the summation index ω ranging over {0, ± 1, …, ± M }. For a bivariate wrapped normal density, setting int.displ = M will ensure that the infinite double sum is approximated by a finite double sum, with the paired summation index (ω_1, ω_2) ranging over {0, ± 1, …, ± M }^2.
   

Random data generation from the von Mises models (both univariate and bivariate) is done via   rejection sampling. In the univariate case, the von Mises random deviates are efficiently generated using a rejection sampling scheme from a wrapped Cauchy distribution <cit.>. For the bivariate models, two forms of random samplings are implemented. In the first method, random deviates are generated via a naive bivariate rejection sampler with uniform proposal density (the majorization constant is numerically evaluated). In the second method (proposed in the web appendix of <cit.>), random deviates are first generated from the marginal distribution of one coordinate, then the other coordinate is drawn from the corresponding conditional distribution (which is von Mises in both models). The authors note that this latter scheme has a typical efficiency rate of at least 60%. It is to be noted that while this scheme is usually more efficient than the naive rejection sampler (especially when the concentration is high), it does have an often substantial overhead due to the numerical computations required for determining appropriate proposal density parameters. These overheads often outweigh efficiency gains, especially if the sample size and/or concentration parameters are small. In BAMBI, therefore, the naive rejection sampler is used by default when the sample size is moderate or small (< 100), or when the concentration parameters are small (< 0.1).[When the concentration parameters are large, the density becomes concentrated in (a) very narrow region(s). As such, efficiency of the naive sampler, which draws proposal random deviates from a uniform density over the entire support, can be 15-20% or less. However, even then the overall runtimes of the naive method are often still comparable to <cit.>'s method when the sample size is moderate or small.] 
For wrapped normal distributions (both univariate and bivariate) a random deviate is easily obtained by sampling from the unwrapped normal distribution (using rnorm if univariate, and rvmnorm from package mvtnorm if bivariate), and then wrapping into [0, 2π).



  
Functions for evaluating density and generating random samples from a finite mixture model with a fixed number of components.  Analogous to the functions for single component densities, the functions dvmmix, dwnormmix, dvmsinmix, dvmcosmix and dwnorm2mix evaluate the density and the functions rvmmix, rwnormmix, rvmsinmix, rvmcosmix and rwnorm2mix  generate random samples from mixtures of vm, wnorm, vmsin, vmcos and wnorm2 respectively. All model parameters and mixing proportions must be provided as input arguments.



  
Functions for visualizing and  summarizing bivariate  models.  To visualize the density for any of the three bivariate angular mixture models (with specified parameters and number of components) considered in this paper, the functions surface_model and contour_model can be used, which respectively plot the surface and the contour of a mixture density. To compute summary statistics for a single bivariate  angular distribution, the function circ_varcor_model can be used, which calculates the circular variance and  correlation coefficients (both Jammalamada-Sarma and Fisher-Lee forms, see Section <ref>).  However, summarizing angular mixture models via circular variances and correlations is not recommended, as interpretations of the results can be challenging when multiple clusters are present in the data.[To calculate the circular variances and correlations for a mixture density, one can simulate from the density first, and then approximate the population quantities by their sample analogs on the basis of the simulated data.] The function circ_cor implements the sample Jammalamada-Sarma and Fisher-Lee circular correlation coefficients, as well as two forms of Kendall's tau <cit.> as non-parametric measures. The sample circular variance can be computed using the var.circular function from R package circular <cit.>.



  
Functions for fitting a single component model or a finite mixture model with specified number of components to a given dataset using MCMC.  Given a dataset, using methods discussed in Section <ref>, the function fit_angmix  generates MCMC samples for parameters in an angular  mixture model with a specified number of components. Available models for bivariate input data (which must be supplied as a two-column matrix or data frame) are vmsin, vmcos and wnorm2, and for univariate data are vm and wnorm.  The argument ncomp specifies number of components  in the mixture model, with ncomp = 1 representing the single component case (i.e., fitting a single density). A Gibbs sampler is used to generate latent component indicators, and conditional on this allocation the model parameters are sampled either by HMC (default), or by RWMH (can be specified through the argument method). A permutation sampling step can be added after burn-in by setting the logical argument perm_sampling to TRUE. The tuning parameters epsilon and L in HMC, and propscale in RWMH have pre-specified default values, and there is an auto-tuning feature for epsilon and propscale which is used by default, but can be turned off by setting the logical argument autotune = FALSE. The burn-in proportion can be specified through the argument burnin.prop, which is set to 0.5 by default. For HMC, the option to use random epsilon and L at each iteration is specified via the logical arguments epsilon.random and L.random respectively. In practice, using multiple chains is recommended, and the argument n.chains specifies the number of chains to be used. These chains can be run in parallel, if the logical argument chains_parallel is set to TRUE. The parallelization is implemented using future_lapply from R package future.apply; an appropriate future::plan() must be set in advance to ensure that the chains run in parallel   (otherwise the chains will run sequentially), see Section <ref> for an example. To retain reproducibility while running multiple chains in parallel, the same RNG state is passed at the beginning of each chain. This is done by specifying future.seed = TRUE in future.apply::future_lapply call. Then at the beginning of the i-th chain, before drawing any parameters, i-many Uniform(0, 1) random numbers are generated using runif(i) (and then thrown away). This ensures that the RNG states across chains prior to sampling of the parameters are different (but reproducible), and hence, no two chains can become identical, even if they have the same starting and tuning parameters. This however, creates a difference between a fit_angmix call with multiple chains which is run sequentially by setting chains_parallel = FALSE, and another call which is run sequentially because of a sequential plan() (or no plan()), with chains_parallel = TRUE. In the former, base::lapply instead of future_lapply is used, which means that different RNG states are passed at the initiation of each chain.

There are options for choosing prior hyperparameters. The prior for the association parameter κ_3 in bivariate models, and the log of the concentration parameters κ (in univariate models), κ_1 and κ_2 (in bivariate models) are taken to be the normal distribution (i.e., the priors for κ, κ_1, κ_2 are log normal), all with zero mean. The default variance for the normal prior is 1000, which provide diffuse priors, although they can be set by the user via the argument norm.var.  A fixed non-informative Uniform(0, 2π) prior is used for the mean parameters. The Dirichlet prior parameters α_j's for the mixing proportions p_j's can be supplied through the argument pmix.alpha, which can either be a positive real number (same for all α_j), or a vector of the same length as pmix. It is recommended that α_j's be chosen large for proper handling of overfitted mixtures;  following <cit.>, all α_j's default to (r+r(r+ 1)/2)/2 + 3, where r denotes dimension of the data (i.e., r=1 for univariate data, and r=2 for bivariate data). See Sections <ref> and <ref> for more details.


The argument cov.restrict specifies any (additional) restriction to be imposed on the component specific association parameters while fitting the model. The available choices are "POSITIVE", "NEGATIVE", "ZERO" and "NONE". Note that when cov.restrict = "ZERO", fit_angmix fits a mixture with product components. By default, cov.restrict = "NONE", which does not impose any (additional) restriction.  When model is "vmsin" or "vmcos", the component densities can be bimodal. However, one can restrict these densities to be unimodal, by setting the logical argument unimodal.component to TRUE (defaults to FALSE). For "wnorm" and "wnorm2" models, the default absolute integer displacement for approximating the Wrapped Normal sum is 3, which can be changed to any value in {1, 2, 3, 4, 5}, through the argument int.displ. For "vmcos" model, the normalizing constant is numerically approximated using quasi Monte Carlo method when analytic evaluation suffers from numerical instability. The arguments qrnd and n_qrnd can be used to alter the default settings used for these approximations. See the documentation of fit_angmix for more details.

The function fit_angmix creates an angmcmc object, which can be used for assessing the fit, post processing, and estimating parameters. 



  
Functions for assessing the fit. Goodness of fit for MCMC-based Bayesian modeling depends on both convergence of the Markov chain and the appropriateness of the model used.  contains a number of functions that can be used to examine these two aspects. The functions paramtrace and lpdtrace respectively plot the parameter and log posterior density traces for visual assessment of convergence. These two functions are called together in the plot method for angmcmc objects. The as.mcmc.list method for angmcmc objects provides a convenient way of converting an angmcmc object to an mcmc.list object from package coda, which provides several additional functions for convergence diagnostics. 
 
Once convergence is justified, the appropriateness of the fitted model can be visually assessed by the S3 functions densityplot (from lattice)   and contour. The first function plots the density surface (for bivariate data) or density curve (for univariate data) of the fitted mixture model, and the second plots the associated contour of a bivariate model. Note that these plots provide visual assessment of the goodness of fit by assuming the Markov chains have converged and the parameters can be estimated on the basis of the MCMC samples. As such, convergence of the MCMC samples must be ensured prior to this step. Otherwise these visual diagnostics will lead to misleading conclusions. 

The comparative goodness of fit for two mixture models can be assessed on the basis of model selection criteria implemented in BAMBI, namely, marginal likelihood, AIC, BIC, DIC, WAIC and LOOIC via the functions bridge_sampler.angmcmc, AIC, BIC, DIC, waic.angmcmc and loo.angmcmc. 
As with the diagnostic plots, one general caveat for using any of these model selection criteria is that one should ensure convergence of the associated Markov chain first; otherwise, the results may be misleading.

 

  
Functions for post-processing and estimating parameters.
  BAMBI provides several post-processing functions to aid inference on the basis of the generated MCMC samples. The function add_burnin_thin adds additional burn-in and/or thinning to an angmcmc object and the function select_chains extracts a subset of chains. These two functions can be helpful if convergence diagnostics indicate that some of the chains are poor mixing and/or require additional burn-in and thinning for convergence.
  
  As described in Section <ref>, care should be taken to ensure that there is no label switching in the MCMC output if inference is being made on the basis of posterior mean/median. If present, label switching can be fixed by applying the wrapper function fix_label on the angmcmc object, which will output another angmcmc object with label switching fixed.  
  
   Point estimates of the parameters are obtained using the function pointest on a fitted angmcmc object.  The function pointest calculates point estimates by applying fn on the MCMC samples, where fn is either a function, or a character string specifying the name of the function. Default for fn is mean, which computes posterior mean. If fn is "MODE" or "MAP" then the (MCMC-based approximate) MAP estimate is returned. Posterior quantiles can be estimated by (sample) quantiles of the MCMC realizations using the S3 function quantile.angmcmc. These quantiles can be used to construct credible sets. For example, if ξ_ζ denotes the ζ-th (sample) quantile of the MCMC observations for 0 < ζ < 1, the central 95% credible interval is given by (ξ_0.025, ξ_0.975). Both of these functions can be applied on specific parameters and/or component labels by setting the arguments par.name and comp.label accordingly. The S3 function summary.angmcmc prints (estimated) posterior means and the central 95% credible intervals for all the parameters.
   
   The estimated latent allocation from an angmcmc object can be obtained using the function latent_allocation, which first estimates parameters via pointest, computes posterior membership probabilities (see (<ref>)) for each data point, and then assigns each data point the class with largest membership probability.  The (estimated) log-likelihood of an angmcmc object can be extracted as a logLik object using the S3 function logLik.angmcmc. Note that there are two methods of obtaining the log-likelihood from an angmcmc object. In the default method (method = 1), the final  log-likelihood is computed by applying a function fn (defaults to max) on the iteration wise log-likelihood values obtained during the original MCMC sampling. On the other hand, if method = 2, first the parameters are estimated (using pointest), and then the log-likelihood is computed at the estimated parameters. 

   Density evaluations and random data generation from a fitted model can be done using the functions d_fitted and r_fitted respectively. Both functions take an angmcmc object as input and apply the appropriate model specific density evaluation and random data generation functions with the estimate   of the parameter vector  obtained via pointest. The actual MCMC samples for one or more parameters in one or more components from one or more chains can be extracted via the function extractsamples on an angmcmc object for further analysis.

   

  
Functions for incremental mixture model fitting and number of components estimation.  Using the methods and model selection criteria described in Sections <ref> and <ref>, the  function fit_incremental_angmix fits mixture models with incremental number of components by calling fit_angmix at each step, and uses a Bayesian model selection criterion to determine an optimal number of components.  The arguments start_ncomp and max_ncomp provide the starting and maximum number of components to be used in the incremental fitting, which are respectively set to 1 and 10 by default. The available  model selection criterion to use (specified via the argument crit) are 'LOGML', 'AIC', 'BIC', 'WAIC', 'LOOIC' and 'DIC', which is computed for every intermediate fit. The initial values of the starting model (or a model with ≤ 2 components) are obtained by default using moment estimation on k-means clusters (they can also be directly supplied by the user). For the subsequent models (with number of components ≥ 3), the initial values are by default obtained from the MCMC-based MAP parameter estimates for the previous model with one fewer component (see Section <ref>).  This can be overridden by setting prev_par = FALSE, to use k-means clustering followed by moment estimation instead.   By default, only the “best” chain, i.e., the one  with maximum average log posterior density, is used for computation of model selection criterion, and parameter estimation (if prev_par is set to TRUE). This default helps safeguard against situations where some of the chains may get trapped at local optima. However, samples from all chains can be used for these computations by setting use_best_chain = FALSE. The function stops when crit achieves its first minimum, or when max_ncomp is reached, and returns a list with the following elements:
   
     
  * fit.best is an angmcmc object corresponding to the optimum or best fit.

     
  * crit.all provides a vector (list) of model selection criterion values for each incremental model fitted.

     
  * crit.best is the value of the model selection criterion for the model with optimal number of components.
     
     
  * maxllik.all is the maximum (obtained from MCMC iterations)  log-likelihood for all fitted models.
       
     
  * maxllik.best is maximum log-likelihood for the optimal model.

     
  * ncomp.best is the optimal number of components associated with the “best” model.
     
     
  * fit.all  is a list consisting of angmcmc objects for all number of components fitted during the model selection process. Any element of this list can be used as an argument for any function that takes angmcmc objects as input. However, this can be very memory intensive, and as such, by default not returned (can be returned by setting return_all = TRUE).

   

The angmcmc object  corresponding to the best fit and the associated value of the model selection criterion  can also be extracted from the output of fit_incremental_angmix  extracted using the convenience functions bestmodel and bestcriterion respectively.






§ ILLUSTRATIONS
 

In this section we illustrate functionalities of by fitting mixture models to the angular datasets included in the package.  The following command

R> library(BAMBI)

loads the package after it has been installed.  For reproducibility of the results presented, the same random seed 12321 is used for all of our examples.



 §.§ Fitting mixture models on the tim8 bivariate data
  
The tim8 dataset consists of 490 backbone torsion angle pairs (ϕ, ψ) for the protein 8TIM. The protein is an example of a TIM barrel, a common type of protein fold exhibiting alternating α-helices and β-sheets. Its Ramachandran plot (i.e., a scatterplot of (ϕ, ψ) pairs) is generated using the following R command and shown in Figure <ref>.

R> plot(tim8, pch = 16,
+       xlim = c(0, 2*pi), ylim = c(0, 2*pi),
+       main = "8TIM", 
+       col = scales::alpha("black", 0.6))


Note that this scatterplot projects the torus onto a 2D surface, which cannot show the wraparound nature of the angles.  This projection is not unique and affects the appearance of the scatterplot depending on how the angles are represented, e.g., in [-π, π) instead of [0, 2π). Moreover, one should be careful to note that the  top and bottom boundaries in these plots join together, as do the left and right boundaries. In Figure <ref>, about 3-6 visually distinct clusters can be seen; however, we need to note that the points around (4.5, 0) and (5.5, 6) in fact may form a single cluster.  Such features cannot be correctly modeled with statistical methods that ignore  circularity. Thus, this is a suitable example for illustrating the need for mixtures of (bivariate) angular distributions.





To fit a bivariate mixture model with a specified number of components to the dataset, we can use the BAMBI function fit_angmix by specifying a model and the number of components to be used.  For example, the following R command fits a 4 component vmsin mixture by generating 3 MCMC chains with 20,000 samples each for the mixture model parameters. HMC  is used for sampling the model parameters, with tuning parameter epsilon adaptively tuned during burn-in (which is by default constituted by the first half of all iterations, i.e., first 10,000 iterations), and L taking its default value 10.  A fit_angmix call creates an angmcmc object, which can then be used for various post-processing tasks, including convergence assessment, parameter estimation and visualizing goodness of fit. 



R> set.seed(12321)
R> fit.vmsin.4comp <- fit_angmix("vmsin", tim8, ncomp = 4, n.iter =  2e4,
+                                n.chains = 3)


Note that in order for the independent chains to be run in parallel, an appropriate plan() from R package future needs to be set first; otherwise the chains will run sequentially. For example, running the commands 

R> library(future)
R> plan(multiprocess, gc = TRUE)

before the fit_angmix call will ensure that the three chains are run in parallel, provided resources are available.  We suggest setting gc = TRUE in plan, to allow proper garbage collection from the parallel workers, even though it adds some overhead. This is because the parallel workers can end up leaving heavy memory footprints,  especially when mixture models are being fitted incrementally.




In the previous example with fit_angmix, the number of components K was specified through ncomp. However, the “true” K generally needs to be estimated from the data.  For this purpose, we use fit_incremental_angmix with start_ncomp = 2, which fits angular mixtures with incremental number of components (starting with 2 components), and uses a Bayesian model selection criterion to determine an optimal model. We use this function to fit optimal mixtures of vmsin, vmcos and wnorm2 models separately. By specifying n.chains = 3 and n.iter = 20,000 in fit_incremental_angmix, each incremental model will be fitted using three chains with 20,000 iterations each (with first 10,000 iterations treated as burn-in, where epsilon in HMC is tuned).

By default, the algorithm uses MCMC-based MAP estimates from the preceding fitted model with one fewer component (if available) as starting parameter values (see Section <ref>). We use the default leave one out cross validation information criterion ('LOOIC') as the model selection criterion for determining the optimal model. 
When the function stops, we extract the best fitted model using  the function bestmodel, and assess convergence of the associated chains. After justifying convergence, we provide point and interval estimates of the parameters and visually examine goodness of fit.  




  
Fitting the vmsin mixture model.

We start with the vmsin model. The R commands are as follows:


R> set.seed(12321)
R> fit.vmsin <- fit_incremental_angmix(model="vmsin", data = tim8,
+                                      crit = "LOOIC",
+                                      start_ncomp = 2,
+                                      max_ncomp = 10,
+                                      n.iter = 2e4,
+                                      n.chains = 3)



The algorithm stops at 5 components and determines the optimal number of components to be 4 on the basis of the 'LOOIC' values. The MCMC-based maximum log-likelihood estimates for the intermediate models are  -945.3684, -853.3546, -803.4193, and -793.8154, which are steadily increasing. This is expected, since a “smaller” mixture should be nested within a “larger” mixture when properly fitted. We extract the optimum fitted model from the output of fit_incremental_angmix via  bestmodel:

R> fit.vmsin.best <- bestmodel(fit.vmsin)	

	
Before estimating parameters, we first need to assess convergence and stationarity of the Markov chains. For this purpose, we first look at the (non-normalized) log posterior density (LPD) trace plots, which can be obtained using the function lpdtrace:

R> lpdtrace(fit.vmsin.best)






The resulting plot displayed in Figure <ref> shows that all three chains have stabilized into similar LPD ranges  after burn-in, without noticeable trends or patterns.  Next, we look at the parameter traces, plotted using function paramtrace:

R> paramtrace(fit.vmsin.best)

These trace plots are displayed in the panels of Figure <ref> through Figure <ref> in Appendix <ref>, which show adequate signs of convergence and stationarity for samples within each chain. Stationarity of a chain can be formally tested using Geweke's convergence diagnostic  <cit.>, which tests equality (of means) of the first and last part of a Markov chain using   standard z-scores. The test is implemented in R package coda, and can be applied on fit.vmsin.best, first by converting it into a  coda mcmc.list object via S3 function as.mcmc.list:

R> mcmc.vmsin.best <- coda::as.mcmc.list(fit.vmsin.best)

and then by applying the coda function geweke.diag on the output. We apply geweke.diag on mcmc.vmsin.best by setting both frac1 and frac2 equal to 0.5, to test the equality of the first and second halves. This produces a list of size 3 containing z-statistics for each chain. These values are displayed in  Figure <ref> as barplots. The R commands are as follows:

R> geweke_res <- coda::geweke.diag(mcmc.vmsin.best, frac1 = 0.5, 
+                                  frac2 = 0.5)
R> par(mfrow = c(1, 3), mar = c(5, 6, 2, 1))
R> for(j in 1:3) 
+      barplot(geweke_res[[j]]z, horiz = TRUE,
+              names = names(geweke_res[[j]]z),
+              xlim = range(geweke_res[[j]]z, -3, 3),
+              ylab="", xaxt='n', las = 2)
+      axis(1, las = 1)
+      title(main = paste("Chain", j), xlab = "geweke.diag")
+  


From these barplots, we see that all the z-scores more or less lie between±3, thus indicating adequate similarity between the first and second halves of the chains. 

Package coda provides functions for a number of additional diagnostic plots and formal tests that may be used after conversion to a mcmc.list object.  Note that not all tests are applicable in every situation.  For example, the Gelman-Rubin test <cit.>, available via coda function gelman.diag, assumes normality of the posterior density; this assumption is clearly violated when the posterior density is multimodal. In this example, as well as for mixture models with several components in general, multimodality in the posterior density can be commonly observed.

Posterior multimodality is evident in the parameter traceplots for the current fit:  the sample values of the same parameter across the independent chains exhibit noticeable differences (see, e.g., the panels of Figure <ref>). Note that this can be due to both permutation (i.e., label-switching) and non-permutation modes. These modes are from similar posterior density regions, as the  LPD traces show. This also indicates that various regions of the posterior density are being well explored by the three chains together. 

Next, we consider parameter estimation and assessing goodness of fit. Since  the combined MCMC samples are multimodal, the posterior mean point estimates from raw MCMC samples will not be meaningful, as they will lie in between modes. A simple alternative is to use the MCMC based approximate MAP estimate, which is unaffected by multimodality; proper care must be taken otherwise. The inferential difficulties associated with having permutation modes in the MCMC samples can be (potentially) solved by undoing label switching. The  function fix_label (with the default settings) can be used for this purpose:

R> fit.vmsin.best <- fix_label(fit.vmsin.best)

The parameter traceplots obtained (using paramtrace after undoing label switching) are displayed in Figures <ref> through <ref> in Appendix <ref>. Compare these traceplots to the ones displayed in Appendix <ref>. It can be seen that this procedure indeed removes some of the permutation modes, in that the parameter traces for the three independent chains now largely overlap. However, some non-unique modes are still present, as can be seen  in the traces ofκ_3andμ_2in component 4, displayed in panel (d) and (f) of Figure <ref>: e.g.,κ_3“jumps” between modes at approximately 2 and -2 over the course of the MCMC simulation. These modes might be genuine non-permutation modes, or permutation modes that fix_label is unable to resolve.


In BAMBI,  parameter estimates are computed using the function pointest, which can find point estimates of the whole parameter vector, as well as its sub-vectors. Note that the function supports multiple methods of estimation. In particular, the argument fn in pointest specifies what function to evaluate on the MCMC samples for estimation. For example, fn = mean computes the MCMC posterior mean, while fn = "MODE" returns an MCMC based approximate MAP estimate. We use pointest to find the MAP and posterior mean estimates (after applying fix_label), and then note their differences. 
The R commands are as follows.

R> round(pointest(fit.vmsin.best, fn = "MODE"), 2)

            1     2     3     4
pmix     0.43  0.16  0.35  0.06
kappa1  36.26  6.15  4.63  3.72
kappa2  27.93  2.25  7.80  0.00
kappa3 -12.03 -0.45 -0.58 -1.85
mu1      5.22  4.66  4.46  1.84
mu2      5.54  6.14  2.41  4.94

R> round(pointest(fit.vmsin.best, fn = mean), 2)

            1     2     3     4
pmix     0.43  0.17  0.34  0.06
kappa1  36.39  7.35  4.34  4.32
kappa2  29.19  2.06  8.12  0.07
kappa3 -12.74 -0.38 -1.08 -0.19
mu1      5.23  4.67  4.43  1.71
mu2      5.55  6.17  2.43  3.60


We note that both the approximate MAP estimate and the posterior mean estimate reasonably agree on the first three components. However, they disagree on the remaining fourth component regarding the value of mu2 and kappa3. This is not surprising, since the MCMC samples for these two parameters have (possibly non-permutation) multiple modes, as we saw earlier. In this case, their posterior mean estimates lie in between modes, and hence are not good point estimates. 

To visualize the differences between these estimates, we plot the  contours and  surfaces of the corresponding fitted model densities using the S3 functions contour (from graphics) and densityplot (from lattice) for angmcmc objects:

R> contour(fit.vmsin.best, fn = "MODE")
R> lattice::densityplot(fit.vmsin.best, fn = "MODE")
R> contour(fit.vmsin.best, fn = mean)
R> lattice::densityplot(fit.vmsin.best, fn = mean)




The contour plots are shown in Figures <ref> and <ref>, and the density surfaces in Figures <ref> and <ref>. As can be seen, these plots are visually quite similar, despite the differences in the point estimates. This is due to the fact the component mostly affected by the existence of multiple modes has a low mixing proportion.  Nonetheless, in the current setting, the MAP estimate is the better of the two for the reasons described.  

Finally, we compute interval estimates of the parameters. This done by the S3 function summary for angmcmc objects, which computes the MCMC posterior mean along with a 95% credible interval: 

R> summary(fit.vmsin.best)

                            1                   2
pmix        0.43 (0.38, 0.48)   0.17 (0.12, 0.22)
kappa1   36.39 (29.01, 45.41)  7.35 (4.99, 10.68)
kappa2   29.19 (21.99, 38.96)   2.06 (1.21, 3.16)
kappa3 -12.74 (-18.49, -7.50) -0.38 (-1.76, 1.16)
mu1         5.23 (5.20, 5.26)   4.67 (4.54, 4.80)
mu2         5.55 (5.51, 5.58)   6.17 (5.99, 6.28)
                          3                    4
pmix      0.34 (0.29, 0.38)   0.064 (0.043, 0.09)
kappa1    4.34 (3.45, 5.43)     4.32 (1.64, 7.76)
kappa2   8.12 (6.00, 10.61) 0.073 (0.00012, 0.55)
kappa3 -1.08 (-2.29, 0.025)   -0.19 (-3.10, 3.05)
mu1       4.43 (4.35, 4.52)     1.71 (1.44, 2.05)
mu2       2.43 (2.37, 2.49)     3.60 (0.66, 5.97)




























We also use this example to illustrate r_fitted, which generates random deviates from a fitted model, with parameters estimated using pointest. The corresponding function d_fitted evaluates the density. These can be useful for posterior predictive checks. We draw observations from the best (4 component) fitted vmsin model, construct the Ramachandran plot for the generated dataset (exhibited in Figure <ref>) and compare it with the original Ramachandran plot. The following are the R commands used for this purpose.


R> set.seed(12321)
R> vmsin.data <- r_fitted(nrow(tim8), fit.vmsin.best, fn = "MODE")
R> plot(vmsin.data, xlab = "phi", ylab = "psi",
+       xlim = c(0, 2*pi), ylim = c(0, 2*pi), 
+       pch = 16, col = scales::alpha("black", 0.6))
R> title("Data generated from best fitted vmsin")



Observe the similarity between Figures <ref> and <ref>. The different clusters of the actual data points are reproduced well in the simulated data, which corroborate the goodness of fit.







  
Fitting the vmcos mixture model.

Next, we fit vmcos mixtures to the data, by using fit_incremental_angmix with model = "vmcos". We set n_qrnd = 1e4 (the default used in dvmcos), which specifies that 10,000 pairs of quasi-random Sobol numbers would be used to approximate the vmcos normalizing constant in cases where its analytic computaion is unstable. In a small dimensional problem with finite variance (such as ours), the Sobol sequence (or low discrepancy quasi-random sequences in general) often provides a better Monte Carlo approximation than (pseudo-) random sequences. In fact, for two dimensional problems, the rate of convergence for a Sobol sequence based Monte Carlo approximation isO((logN)^2/N)as opposed toO(1/√(N))for an ordinary (pseudo-) random sequence based Monte Carlo approximation <cit.>, whereNdenotes the number of (quasi) random pairs used.  See the documentation of dvmcos for examples comparing analytic, quasi Monte Carlo, and ordinary Monte Carlo approximations  of the (normalizing constant of the) vmcos density. For fit_incremental_angmix (or more specifically in fit_angmix),  Monte Carlo approximations based on10^4pairs of Sobol numbers typically provide reasonable approximations, while keeping the computational burden moderate. 

The following are the R commands used for incrementally fitting  vmcos  mixture models.


R> set.seed(12321)
R> fit.vmcos <- fit_incremental_angmix(model="vmcos", data = tim8,
+                                      crit = "LOOIC",
+                                      start_ncomp = 2,
+                                      max_ncomp = 10,
+                                      n.iter = 2e4,
+                                      n.chains = 3,
+                                      use_best_chain = FALSE,
+                                      n_qrnd = 1e4)



Similar to the vmsin case, here also the algorithm stops at 5 components and determines the optimal number of components to be 4. We first extract the “best” fitted model, via bestmodel:

R> fit.vmcos.best <- bestmodel(fit.vmcos)

and then plot the log posterior and parameter traces. These plots show similar convergence properties, and are omitted for brevity. For parameter estimation, we compute both (approximate) MAP and posterior mean estimates (after undoing label switching), and plot the contour and surface of the associated fitted model densities. The following are the associated R commands:

R> fit.vmcos.best <- fix_label(fit.vmcos.best)
R> contour(fit.vmcos.best, fn = "MODE")
R> lattice::densityplot(fit.vmcos.best, fn = "MODE")
R> contour(fit.vmcos.best, fn = mean)
R> lattice::densityplot(fit.vmcos.best, fn = mean)




The fitted contours are displayed in Figures <ref> and <ref>, and the density surfaces are displayed in Figures <ref> and <ref>, are noticeably similar.  They  are also broadly similar to the ones obtained for the fitted vmsin models.  Estimated posterior means along with estimated 95% credible intervals are   obtained using the S3 function summary.angmcmc as follows.


















R> summary(fit.vmcos.best)

                            1                     2
pmix        0.43 (0.38, 0.49)     0.16 (0.12, 0.21)
kappa1   48.70 (37.89, 61.88)    7.97 (5.11, 12.56)
kappa2   41.49 (30.79, 55.64)    2.15 (0.002, 4.57)
kappa3 -12.61 (-18.64, -7.35) -0.0077 (-1.90, 2.12)
mu1         5.23 (5.20, 5.26)     4.67 (4.54, 4.81)
mu2         5.55 (5.51, 5.58)     6.17 (5.99, 6.28)
                         3                    4
pmix     0.34 (0.29, 0.38) 0.065 (0.043, 0.095)
kappa1   5.16 (3.80, 6.69)    3.62 (0.90, 6.77)
kappa2  8.79 (6.36, 11.81) 0.19 (0.00013, 1.33)
kappa3 -0.98 (-2.15, 0.12)  -0.20 (-1.54, 0.90)
mu1      4.43 (4.34, 4.51)    1.59 (1.31, 1.93)
mu2      2.43 (2.36, 2.49)    3.13 (0.14, 6.13)


  
Fitting the wnorm2 mixture model.

Finally, we fit wnorm2 mixtures to the data. The R commands used are as follows:


R> set.seed(12321)
R> library(future)
R> plan(multiprocess(workers = 3))
R> fit.wnorm2 <- fit_incremental_angmix(model="wnorm2", data = tim8,
+                                       crit = "LOOIC",
+                                       start_ncomp = 2,
+                                       max_ncomp = 10,
+                                       n.iter = 2e4,
+                                       n.chains = 3,
+                                       use_best_chain = FALSE)


Here also, the function stops at 5 components and determines the optimal number of components to be 4.[It should be noted that the runtime for wrapped normal fitting is considerably longer than the von Mises sine models, due to the computational burden; see Section <ref>.] As done in the previous two cases, after extracting the best model, we assess convergence via trace plots (omitted for brevity). We  find MCMC-based MAP and posterior mean estimates (after undoing label switching), and also find credible interval estimates. Finally we assess goodness of fit through fitted contour and density surfaces. The following are the R commands that perform these tasks.


R> fit.wnorm2.best <- bestmodel(fit.wnorm2)
R> lpdtrace(fit.wnorm2.best)
R> paramtrace(fit.wnorm2.best)
R> fit.wnorm2.best <- fix_label(fit.wnorm2.best)
R> contour(fit.wnorm2.best, fn = "MODE")
R> lattice::densityplot(fit.wnorm2.best, fn = "MODE")
R> contour(fit.wnorm2.best, fn = mean)
R> lattice::densityplot(fit.wnorm2.best, fn = mean)
R> summary(fit.wnorm2.best)


The contours and density surfaces  displayed in Figures <ref> and <ref>) are noticeably  similar.  They are also broadly similar to the fitted vmsin and vmcos mixture model density contours and surfaces. The estimated credible intervals along with MCMC posterior means for fitted 4 component wnorm2 mixture are obtained using the S3 function summary.angmcmc as follows.


R> summary(fit.wnorm2.best)

                          1                  2
pmix      0.44 (0.38, 0.49)  0.16 (0.12, 0.21)
kappa1 35.57 (28.05, 44.49) 7.13 (4.55, 10.90)
kappa2 28.03 (20.93, 36.91)  1.48 (0.75, 2.44)
kappa3  12.32 (7.38, 18.02) 0.20 (-0.90, 1.12)
mu1       5.23 (5.20, 5.26)  4.66 (4.54, 4.78)
mu2       5.55 (5.52, 5.58)  6.17 (5.98, 6.28)
                        3                    4
pmix    0.34 (0.29, 0.38) 0.063 (0.043, 0.087)
kappa1  3.68 (2.86, 4.65)   8.05 (2.42, 17.30)
kappa2 7.68 (5.62, 10.32)   0.46 (0.034, 0.91)
kappa3  1.00 (0.10, 2.00)   1.13 (-1.36, 3.09)
mu1     4.43 (4.35, 4.51)    1.62 (1.37, 1.94)
mu2     2.43 (2.37, 2.50)   3.06 (0.044, 6.23)


  
Comparative analysis of the three fitted models.

So far, we have considered mixtures of vmsin, vmcos and wnorm2 densities, and have fitted them to tim8 data. The associated optimal number of components were determined via leave one out information criterion (LOOIC) in incremental fitting schemes. We then plotted the fitted density contours and  surfaces to assess the goodness of fit, and noticed that these plots are broadly similar across the three optimal fitted mixture models (of vmsin, vmcos and wnorm2 densities). It is natural to then consider the question of which among these three fitted bivariate mixture models best explains the data. This can again be answered via LOOIC. For angmcmc objects LOOIC can be conveniently computed using the S3 function looic from package loo. However, here we do not need to recompute them since they were already computed during incremental model fitting, and can be extracted via the convenience function bestcriterion as follows:

R> vmsin.4.looic <- bestcriterion(fit.vmsin)
R> wnorm2.4.looic <- bestcriterion(fit.wnorm2)
R> vmcos.4.looic <- bestcriterion(fit.vmcos)


Now we compare the three models via loo::compare on the basis of their LOOIC's:

R> comp <- loo::compare(vmsin.4.looic, vmcos.4.looic, wnorm2.4.looic)
R> comp

               elpd_diff se_diff elpd_loo p_loo  looic 
vmsin.4.looic     0.0       0.0  -826.0     25.5 1652.0
wnorm2.4.looic   -0.9       3.6  -826.9     26.5 1653.8
vmcos.4.looic    -5.3       3.1  -831.3     25.2 1662.7


The documentation of loo from package loo v2.1.0 says “When comparing two fitted models, we can estimate the difference in their expected predictive accuracy by the difference in elpd_loo or elpd_waic (or multiplied by -2, if desired, to be on the deviance scale). When that difference, elpd_diff, is positive then the expected predictive accuracy for the second model is higher. A negative elpd_diff favors the first model. When using compare() with more than two models, the values in the elpd_diff and se_diff columns of the returned matrix are computed by making pairwise comparisons between each model and the model with the best ELPD (i.e., the model in the first row)”. 

Thus the above output provides a ranking among the three models based on their (estimated) expected log predictive density (elpd) values; a higher elpd indicates a better predictive accuracy and thus a better fit.  The fitted vmsin model appears to have the highest elpd (see the column elpd_diff in the above output), followed by the fitted wnorm2 model and the fitted vmcos model. However, these elpd's are estimates, and the variabilities of these estimates need to be considered when making comparisons. To address this, we make use of the standard errors of the differences provided in the column se_diff and construct approximate 95% credible interval estimates of the pairwise elpd differences (viz., elpd_diff±2 se_diff) for the fitted model pairs (wnorm2, vmsin) and  (vmcos, vmsin). A elpd difference is considered to be significant (at the 95% level) if the corresponding interval estimate does not contain zero. The R commands are as follows.

R> find_ci <- function(x, digits = 1) 
+    round(c(lower = unname(x[1] - 2*x[2]),
+            upper = unname(x[1] + 2*x[2])),
+          digits = digits)
+ 
R> t(apply(comp[-1, c("elpd_diff", "se_diff")], 1, find_ci))

               lower upper
wnorm2.4.looic  -8.1   6.3
vmcos.4.looic  -11.5   0.9
        
This shows that approximate 95% interval estimates of the elpd differences between the fitted best vmsin and the best wnorm2 model, and the best vmsin and the best vmcos model, are(-8.1, 6.3)and(-11.5, 0.9)respectively, both containing zero. It therefore follows that all three of the fitted (four component) mixture models are not significantly different in terms of their goodness of fit to these data.
































 §.§ Fitting mixture models on the wind (univariate) data


The wind data contains 239 observations on wind direction in radians  measured at Saturna Island, British Columbia, Canada, during October 1-10, 2016. As a result of a severe storm that occurred during that period, the data shows significant variability with an interesting bi- (or possibly tri-) modality. Figure <ref> shows a histogram of the data, constructed by applying the default hist function on wind[, "angle"].





Similar to the bivariate case,  we  use fit_incremental_angmix to fit mixtures of vm and wnorm separately with incremental number of components (starting at 1) and determine an optimum size in each case.  To fit  each mixture model, we first generate 20,000 MCMC samples for the parameters, with the (default) first half taken as burn-in. Except for n.iter, defaults for all other arguments are used in these examples. After generating the MCMC samples, we assess their convergence via LPD and parameter trace plots. Following, we visualize the fits via  density curves constructed using the S3 function densityplot (which requires lattice). Finally, we compute point and interval estimates for each parameter using the S3 function summary.angmcmc.



  
Fitting the vm mixture model.
We start with vm,  The R commands are as follows:

R> set.seed(12321)
R> fit.vm <- fit_incremental_angmix(model="vm", data = wind[, "angle"],
+                                   crit = "LOOIC",
+                                   start_ncomp = 1, max_ncomp = 10,
+                                   n.iter = 2e4,
+                                   n.chains = 3)


The function stops at 3 components and determines the optimal number of components to be 2. After it stops, we extract the angmcmc object corresponding to the best model from its output, and inspect its LPD and parameter traces for convergence (omitted for brevity). 

R> fit.vm.best <- bestmodel(fit.vm)
R> lpdtrace(fit.vm.best)
R> paramtrace(fit.vm.best)
























We first use fix_label to undo label switching, and then assess goodness of fit through density curves fitted using MAP and posterior mean estimation:  

R> fit.vm.best <- fix_label(fit.vm.best)
R> lattice::densityplot(fit.vm.best, fn = "MODE")
R> lattice::densityplot(fit.vm.best, fn = mean)

The plots are displayed in Figures <ref> and <ref>, which show noticeable similarity. 


Finally, we compute MCMC posterior mean and associated 95% credible interval using S3 function summary:

R> summary(fit.vm.best)

                       1                 2
pmix   0.24 (0.13, 0.43) 0.76 (0.57, 0.87)
kappa 7.81 (1.34, 22.40) 1.01 (0.59, 1.71)
mu     5.30 (5.08, 5.49) 2.75 (2.49, 3.04)


















  
Fitting the wnorm mixture model.
Next, we do similar exercises with wnorm model. The following are the R codes used.

R> set.seed(12321)
R> fit.wnorm <- fit_incremental_angmix(model="wnorm", data = wind[, "angle"],
+                                      crit = "LOOIC",
+                                      start_ncomp = 1, max_ncomp = 10,
+                                      n.iter = 2e4,
+                                      n.chains = 3)
R> fit.wnorm.best <- bestmodel(fit.wnorm)
R> lpdtrace(fit.wnorm.best)
R> paramtrace(fit.wnorm.best)
R> fit.wnorm.best <- fix_label(fit.wnorm.best)
R> lattice::densityplot(fit.wnorm.best, fn = "MODE")
R> lattice::densityplot(fit.wnorm.best, fn = mean)


Similar to the "vm" case, here also the function stops at 3 components and determines the optimal number of components to be 2.  The LPD and parameter trace plots are  omitted for brevity. The density curves fitted using MAP and poterior mean estimates are shown in Figures <ref> and <ref> respectively, which are  noticeably similar. They are also  broadly similar to the plots associated with the fitted vm mixture densities shown in Figure <ref>. 




Finally we compute the MCMC posterior mean and 95% credible interval using the S3 function summary.

R> summary(fit.wnorm.best)

                       1                 2
pmix   0.28 (0.16, 0.43) 0.72 (0.57, 0.84)
kappa 5.09 (1.16, 14.21) 0.78 (0.43, 1.35)
mu     5.35 (5.12, 5.53) 2.71 (2.45, 3.00)















  
Comparison between the two models.
Similar to the bivariate case, we  compare the fitted vm and wnorm mixture models using their LOOIC values. We first extract the LOOICs using the convenience function bestcriterion:

R> vm.2.looic <- bestcriterion(fit.vm)
R> wnorm.2.looic <- bestcriterion(fit.wnorm)

Then we compare the two models based on their estimated expected log predictive densities, by using loo::compare on vm.2.looic and wnorm.2.looic:

R> loo::compare(wnorm.2.looic, vm.2.looic)

elpd_diff        se 
      1.0       1.0 

Clearly an approximate 95% credible interval estimate for the elpd difference, obtained by elpd_diff±2 se, contains zero. This implies that the fitted vm model and the fitted wnorm model do not have a significant difference in terms of their goodness of fit to these data.




























§ CONCLUDING REMARKS AND FUTURE WORK


Angular data, both univariate and bivariate, arise naturally in a variety of modern scientific problems, and their analyses require appropriate use of rigorous statistical tools and distributions specifically developed for such data. The lack of comprehensive software implementing such methods (in R or otherwise) has hindered their applicability in practice – especially for bivariate angular models and mixtures thereof.


The package  is our contribution to this area, providing a platform that implements a set of formal statistical tools and methods for analyzing such data, and is readily accessible to  practitioners. 

There are various directions in which the software could be extended in future releases. Some possible features under consideration include the following.


  *  Implementation of additional angular distributions, such as wrapped Cauchy.

  *  Additional methods of density evaluation and random simulation from fitted models on the basis of MCMC samples. 

  *  Visualizations of bivariate angles with toroidal plots.

  *  Use of parallel tempering or related methods during MCMC simulations for faster exploration of the posterior density.

  *  Proper handling of overfitting heterogeneity that takes place in high dimensional mixture models when some of the component specific parameters in two different components are identical. <cit.> suggests the use of sparse priors for the component specific location parameters to deal with this problem. 


§ ACKNOWLEDGMENT

The authors thank the two anonymous referees who reviewed this work, and provided valuable, thorough and constructive suggestions on both the package and the manuscript. The authors are also thankful to Anahita Nodehi from Tarbiat Modares University, Iran, for notifying them of a few typos and providing helpful suggestions on the R package documentation. 









§ THE NORMALIZING CONSTANT FOR VON MISES COSINE DENSITY

The normalizing constant for the density (<ref>) is given by	

    C_c(κ_1, κ_2, κ_3) = [(2π)^2 { I_0(κ_1)I_0(κ_2)I_0(κ_3) + 2 ∑_n=1^∞  I_n(κ_1) I_n(κ_2)  I_n(κ_3)  }]^-1

Without loss of generality, we first assume that the mean parameters in the density (<ref>) are all zero, i.e., μ_1 = μ_2 = 0. Therefore, our objective boils down to evaluate the integral

    C_c(κ_1, κ_2, κ_3)^-1 =  ℐ = ∫_0^2π∫_0^2πexp(κ_1 cos x + κ_2 cos y + κ_3 cos (x-y))  dx   dy.
 Now from equation 9.6.34 of <cit.>, it follows that

    exp(κ_1 cos x)     = I_0(κ_1) + 2 ∑_l=1^∞ I_l(κ_1) cos(lx) 
    exp(κ_2 cos y)     = I_0(κ_2) + 2 ∑_m=1^∞ I_m(κ_2) cos(my) 
    and exp(κ_3 cos (x-y))     = I_0(κ_3) + 2 ∑_n=1^∞ I_n(κ_3) cos(n(x-y)).
 Therefore, the integrand in (<ref>) can be written as

    I_0(κ_1)I_0(κ_2)I_0(κ_3) + 2 {I_0(κ_2) + I_0(κ_3)}∑_l=1^∞ I_l(κ_1) cos(lx) 
          + 2 {I_0(κ_3) + I_0(κ_1)}∑_m=1^∞ I_m(κ_2) cos(my) 
          + 2 {I_0(κ_1) + I_0(κ_2)}∑_n=1^∞ I_n(κ_3) cos(n(x-y)) 
          + 8 ∑_l=1^∞∑_m=1^∞∑_n=1^∞ I_l(κ_1)  I_m(κ_2) I_n(κ_3)  cos(lx)  cos(my)  cos(n(x-y)).
 Note that for any positive integer q,

    ∫_0^2πcos(qz)   dz = ∫_0^2πsin(qz)   dz = 0

which implies, for a positive integer n,

    ∫_0^2π∫_0^2πcos(n(x-y))   dx   dy    = ∫_0^2πcos(nx)   dx ∫_0^2πcos(ny)   dy 
          + ∫_0^2πsin(nx)   dx ∫_0^2πsin(ny)   dy = 0.

(Equality of the double and the iterative integrals are ensured by the Fubini theorem, which is applicable as the integrands and the range of integrals are all finite.)

 Thus the (double) integrals of the second, third and fourth terms in (<ref>) are all zero. Hence,

    ℐ   = (2π)^2  I_0(κ_1)I_0(κ_2)I_0(κ_3)  
          + 8 ∫_0^2π∫_0^2π∑_l=1^∞∑_m=1^∞∑_n=1^∞   I_l(κ_1) I_m(κ_2)  I_n(κ_3)  cos(lx)  cos(my)  cos(n(x-y))  dx  dy.
 Now, for the second term in (<ref>), first note that

    ∫_0^2π∫_0^2π∑_l=1^∞∑_m=1^∞∑_n=1^∞|I_l(κ_1)  I_m(κ_2) I_n(κ_3)  cos(lx)  cos(my)  cos(n(x-y)) |   dx   dy 
       ≤∫_0^2π∫_0^2π∑_l=1^∞∑_m=1^∞∑_n=1^∞ I_l(κ_1)  I_m(κ_2) I_n(|κ_3|)   dx   dy 
       =  ∑_l=1^∞∑_m=1^∞∑_n=1^∞∫_0^2π∫_0^2π I_l(κ_1)  I_m(κ_2) I_n(|κ_3|)   dx   dy   (by Fubini-Tonelli)
       = (2π)^2 ( ∑_l=1^∞ I_l(κ_1)) (∑_m=1^∞  I_m(κ_2) ) ( ∑_n=1^∞   I_n(|κ_3|) ) 
       < ∞

where the equality in the third line follows from the Fubini-Tonelli theorem for non-negative integrands. Therefore, the Fubini theorem for general integrands can be applied to ensure interchangeability of the sums and the integrals in the second term in (<ref>). In particular, one can write

    ∫_0^2π∫_0^2π∑_l=1^∞∑_m=1^∞∑_n=1^∞   I_l(κ_1) I_m(κ_2)  I_n(κ_3)  cos(lx)  cos(my)  cos(n(x-y))  dx  dy 
       = ∑_l=1^∞∑_m=1^∞∑_n=1^∞   I_l(κ_1) I_m(κ_2)  I_n(κ_3) ∫_0^2π∫_0^2πcos(lx)  cos(my)  cos(n(x-y))  dx  dy.


Now, for any positive integers l, m, n,

    cos(lx)  cos(my)  cos(n(x-y))    = cos(lx) cos(nx)  cos(my)  cos(ny) 
          + cos(lx) sin(nx)  cos(my) sin(ny).
 Observe that for any two positive integers p and q,

    ∫_0^2πcos(pz) cos(qz)   dz = π_{p=q} and ∫_0^2πcos(pz) sin(qz)   dz = 0.
 Therefore, for any positive integers l, m, n,

    ∫_0^2π∫_0^2πcos(lx) cos(nx)  cos(my)  cos(ny)   dx  dy = π_{l = n}π_{m = n} = π^2 _{l = m = n}

and

    ∫_0^2π∫_0^2πcos(lx) sin(nx)  cos(my) sin(ny)   dx  dy = 0.

which implies,

    ∫_0^2π∫_0^2πcos(lx)  cos(my)  cos(n(x-y))  dx  dy  = π^2 _{l = m = n}.


Therefore, combining (<ref>), (<ref>) and (<ref>), we get

    ℐ   = (2π)^2  I_0(κ_1)I_0(κ_2)I_0(κ_3) + 8 π^2 ∑_l=1^∞∑_m=1^∞∑_n=1^∞   I_l(κ_1)  I_n(κ_3) I_m(κ_2)  _{l = m = n}
       = (2π)^2 { I_0(κ_1)I_0(κ_2)I_0(κ_3) + 2 ∑_n=1^∞  I_n(κ_1) I_n(κ_2)  I_n(κ_3)  }

This completes the proof.	


§ CIRCULAR VARIANCE AND CORRELATION COEFFICIENTS


 §.§ von Mises sine model


	Let (ψ_1, ψ_2) ∼(μ_1, μ_2, κ_1, κ_2, κ_3). Then 
		
  *  the Fisher-Lee circular correlation coefficient (<ref>) between ψ_1 and ψ_2 is given by 
	        
    ρ_ (ψ_1, ψ_2) 
    		     = (1/C̅_s∂C̅_s/∂κ_3) (1/C̅_s∂^2 C̅_s/∂κ_1 ∂κ_2)/√((1/C̅_s∂^2 C̅_s/∂κ_1^2) (1 - 1/C̅_s∂^2 C̅_s/∂κ_1^2) (1/C̅_s∂^2 C̅_s/∂κ_2^2) (1 - 1/C̅_s∂^2 C̅_s/∂κ_2^2)).

  *  the Jammalamadaka-Sarma circular correlation coefficient (<ref>) between Θ and Φ is given by 
			
    ρ_ (ψ_1, ψ_2) 
    			= 1/C̅_s∂C̅_s/∂κ_3/√((1 - 1/C̅_s∂^2 C̅_s/∂κ_1^2)  (1 - 1/C̅_s∂^2 C̅_s/∂κ_2^2)).

  *  the circular variance for ψ_i, i = 1, 2 is given by
			
    (ψ_i) = 1 - 1/C̅_s∂C̅_s/∂κ_i.
	
	
	   Here C̅_s = 1/C_s, where C_s is the normalizing constant of the von Mises sine distribution as defined in (<ref>). Infinite series expressions for partial derivatives of C̅_s constant are provided as follows.
	       
    ∂C̅_s/∂κ_1   = 4 π^2  ∑_m=0^∞2mm(κ_3^2/4κ_1 κ_2)^m I_m+1(κ_1) I_m(κ_2)  
    ∂C̅_s/∂κ_2   = 4 π^2  ∑_m=0^∞2mm(κ_3^2/4κ_1 κ_2)^m I_m(κ_1) I_m+1(κ_2)  
    ∂C̅_s/∂κ_3   =  8 π^2  ∑_m=1^∞ m 2mmκ_3^2m-1/(4κ_1 κ_2)^m I_m(κ_1) I_m(κ_2) 
    ∂^2 C̅_s/∂κ_1^2   = 4 π^2  ∑_m=0^∞2mm(κ_3^2/4κ_1 κ_2)^m     
               (I_m+1(κ_1)/κ_1 + I_m+2(κ_1)) I_m(κ_2)  
    ∂^2 C̅_s/∂κ_2^2   = 4 π^2  ∑_m=0^∞2mm(κ_3^2/4κ_1 κ_2)^m   
                 I_m(κ_1) (I_m+1(κ_2)/κ_2 + I_m+2(κ_2)) 
    ∂^2 C̅_s/∂κ_1  ∂κ_2   = 4 π^2  ∑_m=0^∞2mm(κ_3^2/4κ_1 κ_2)^m  I_m+1(κ_1) I_m+1(κ_2)


 §.§ von Mises cosine model


	Let (ψ_1, ψ_2) ∼(μ_1, μ_2, κ_1, κ_2, κ_3). Then 
		
  *  the Fisher-Lee circular correlation coefficient (<ref>) between ψ_1 and ψ_2 is given by 
	        
    ρ_ (ψ_1, ψ_2) 
    		     = (1/C̅_c{∂C̅_c/∂κ_3 -  ∂^2 C̅_c/∂κ_1 ∂κ_2}) (1/C̅_c∂^2 C̅_c/∂κ_1 ∂κ_2)/√((1/C̅_c∂^2 C̅_c/∂κ_1^2) (1 - 1/C̅_c∂^2 C̅_c/∂κ_1^2) (1/C̅_c∂^2 C̅_c/∂κ_2^2) (1 - 1/C̅_c∂^2 C̅_c/∂κ_2^2)).

  *  the Jammalamadaka-Sarma circular correlation coefficient (<ref>) between Θ and Φ is given by 
			
    ρ_ (ψ_1, ψ_2) 
    			= 1/C̅_c{∂C̅_c/∂κ_3 -  ∂^2 C̅_c/∂κ_1 ∂κ_2}/√((1 - 1/C̅_c∂^2 C̅_c/∂κ_1^2)  (1 - 1/C̅_c∂^2 C̅_c/∂κ_2^2)).

  *  the circular variance for ψ_i, i = 1, 2 is given by
			
    (ψ_i) = 1 - 1/C̅_c∂C̅_c/∂κ_i .
	
	

	   Here C̅_c = 1/C_c is the reciprocal of the von Mises cosine normalizing constant, as given in (<ref>). Infinite series expressions for partial derivatives of C̅_c are given as follows.
	   
	   		
    ∂C̅_c/∂κ_1 =  4π^2     { I_1(κ_1) I_0(κ_2) I_0(κ_3)  +  .   
           ∑_m = 1 ^∞.    I_m(κ_2) I_m(κ_3) [ I_m+1(κ_1) + I_m-1(κ_1) ] }
    ∂C̅_c/∂κ_2 = 4π^2    {  I_0(κ_1) I_1(κ_2) I_0(κ_3)  +  .   
           ∑_m = 1 ^∞.    I_m(κ_1) I_m(κ_3) [ I_m+1(κ_2) + I_m-1(κ_2) ] }
    ∂C̅_c/∂κ_3 = 4π^2    { I_0(κ_1) I_0(κ_2) I_1(κ_3)  +  .   
           ∑_m = 1 ^∞. I_m(κ_1) I_m(κ_2) [ I_m+1(κ_3) + I_m-1(κ_3) ] }. 
    ∂^2 C̅_c/∂κ_1^2 = 2π^2    {  I_0(κ_2) I_0(κ_3)[I_0(κ_1) + I_2(κ_1)] + .   
           ∑_m = 1 ^∞.   I_m(κ_2) I_m(κ_3) [I_m-2(κ_1) + 2I_m(κ_1) + I_m+2(κ_1)]   }
    ∂^2 C̅_c/∂κ_2^2 = 2π^2    {  I_0(κ_1) I_0(κ_3)[I_0(κ_2) + I_2(κ_2)] + .   
           ∑_m = 1 ^∞.   I_m(κ_1) I_m(κ_3) [I_m-2(κ_2) + 2I_m(κ_2) + I_m+2(κ_2)]   }
    ∂^2 C̅_c/∂κ_1 ∂κ_2 = 2π^2    { 2 I_1(κ_1) I_1(κ_2) I_0(κ_3) + .   
       ∑_m = 1 ^∞.  I_m(κ_3) [ I_m+1(κ_1) + I_m-1(κ_1) ] [ I_m+1(κ_2) + I_m-1(κ_2) ]}


§ GRADIENTS

For notational simplicity we shall omit the subscripts i and j. Note that, in the sequel,  stands for the parameter vector for one generic component and not the entire parameter vector of all components.



 §.§ Wrapped normal models

  * Univariate case. Here ^⊤ = (κ, μ), and

    ∂ f_(ψ| )/∂κ   = 1/2 κ^1/2√(2π)∑_ω∈exp[-κ/2(ψ - μ - 2πω)^2 ] [   1 - κ (ψ - μ - 2πω)^2] 
    ∂ f_(ψ| )/∂μ   = κ^3/2/√(2π)∑_ω∈exp[-κ/2(ψ - μ - 2πω)^2 ] (ψ - μ - 2πω).

  * Bivariate case. Here ^⊤ = (κ_1, κ_2, κ_3, μ_1, μ_2), ^⊤ = (ψ_1, ψ_2) and

    ∂ f_(| )/∂κ_1   = 1/4π√(κ_12.3) ∑_(ω_1, ω_2) ∈^2  E_ω_1, ω_2[ κ_2 - κ_12.3 (ψ_1 - μ_1 - 2πω_1)^2 ] 
    ∂ f_(| )/∂κ_2   = 1/4π√(κ_12.3) ∑_(ω_1, ω_2) ∈^2  E_ω_1, ω_2[ κ_1 - κ_12.3 (ψ_2 - μ_2 - 2πω_2)^2 ]  
    ∂ f_(| )/∂κ_3   = 1/2π√(κ_12.3) ∑_(ω_1, ω_2) ∈^2  E_ω_1, ω_2[ κ_3 - κ_12.3 (ψ_1 - μ_1 - 2πω_1)(ψ_2 - μ_2 - 2πω_2) ]  
    ∂ f_(| )/∂μ_1   = √(κ_12.3)/2π ∑_(ω_1, ω_2) ∈^2  E_ω_1, ω_2[ κ_1(ψ_1 - μ_1 - 2πω_1)  + κ_3(ψ_2 - μ_2 - 2πω_2)]   
    ∂ f_(| )/∂μ_2   = √(κ_12.3)/2π ∑_(ω_1, ω_2) ∈^2  E_ω_1, ω_2[ κ_3(ψ_1 - μ_1 - 2πω_1)  + κ_2(ψ_2 - μ_2 - 2πω_2)]

where

    E_ω_1, ω_2  =   exp   [-1/2{κ_1 (ψ_1 - μ_1 - 2πω_1)^2 + κ_2 (ψ_2 - μ_2 - 2πω_2)^2  . .
             . . + 2 κ_3 (ψ_1 - μ_1 - 2πω_1) (ψ_2 - μ_2 - 2πω_2) }]

and κ_12.3 = κ_1κ_2 - κ_3^2.


 §.§ von Mises models

  * Univariate case. Here ^⊤ = (κ, μ) and

    ∂log f_(ψ| )/∂κ   = cos(ψ - μ) - I_1(κ)/I_0(κ)
    ∂log f_(ψ| )/∂μ   = κsin(ψ - μ).

  * Bivariate sine model. Here ^⊤ = (κ_1, κ_2, κ_3, μ_1, μ_2), ^⊤ = (ψ_1, ψ_2) and

    ∂log f_(| )/∂κ_1   = cos(ψ_1 - μ_1) - ∂C̅_s(κ_1, κ_1, κ_3)/∂κ_1/C̅_s(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂κ_2   = cos(ψ_2 - μ_2) - ∂C̅_s(κ_1, κ_1, κ_3)/∂κ_2/C̅_s(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂κ_3   = sin(ψ_1 - μ_1) sin(ψ_2 - μ_2) - ∂C̅_s(κ_1, κ_1, κ_3)/∂κ_3/C̅_s(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂μ_1   = κ_1 sin(ψ_1 - μ_1) - κ_3 cos(ψ_1 - μ_1) sin(ψ_2 - μ_2) 
    ∂log f_(| )/∂μ_2   = κ_2 sin(ψ_2 - μ_2) - κ_3 sin(ψ_1 - μ_1) cos(ψ_2 - μ_2)

where C̅_s(κ_1, κ_1, κ_3) = 1/C_s(κ_1, κ_1, κ_3) and expressions for the partial derivatives are provided in Appendix <ref>.











  * Bivariate cosine model. Here ^⊤ = (κ_1, κ_2, κ_3, μ_1, μ_2), ^⊤ = (ψ_1, ψ_2) and

    ∂log f_(| )/∂κ_1   = cos(ψ_1 - μ_1) - ∂C̅_c(κ_1, κ_1, κ_3)/∂κ_1/C̅_c(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂κ_2   = cos(ψ_2 - μ_2) - ∂C̅_c(κ_1, κ_1, κ_3)/∂κ_2/C̅_c(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂κ_3   = cos(ψ_1 - μ_1 - ψ_2 + μ_2) - ∂C̅_c(κ_1, κ_1, κ_3)/∂κ_3/C̅_c(κ_1, κ_1, κ_3)
    ∂log f_(| )/∂μ_1   = κ_1 sin(ψ_1 - μ_1) + κ_3 sin(ψ_1 - μ_1 - ψ_2 + μ_2) 
    ∂log f_(| )/∂μ_2   = κ_2 sin(ψ_2 - μ_2) - κ_3 sin(ψ_1 - μ_1 - ψ_2 + μ_2)

where C̅_c(κ_1, κ_1, κ_3) = 1/C_c(κ_1, κ_1, κ_3) and infinite series expressions for the partial derviatives are provided in Appendix <ref>.









§ TRACE PLOTS FOR 4 COMPONENT VMSIN


§ TRACE PLOTS FOR 4 COMPONENT VMSIN WITH LABEL SWITCHINGS FIXED