{'fraction_non_alphanumeric': 0.034107760751359364, 'fraction_numerical': 0.018801638302379776, 'mean_word_length': 4.545276554087127, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Minimum Message Length based Mixture Modelling using Bivariate von Mises Distributions with Applications to Bioinformatics

The modelling of empirically observed data is commonly done using mixtures of probability distributions. In order to model angular data, directional probability distributions such as the bivariate von Mises (BVM) is typically used. The critical task involved in mixture modelling is to determine the optimal number of component probability distributions. We employ the Bayesian information-theoretic principle of minimum message length (MML) to distinguish mixture models by balancing the trade-off between the model's complexity and its goodness-of-fit to the data. We consider the problem of modelling angular data resulting from the spatial arrangement of protein structures using BVM distributions. The main contributions of the paper include the development of the mixture modelling apparatus along with the MML estimation of the parameters of the BVM distribution. We demonstrate that statistical inference using the MML framework supersedes the traditional methods and offers a mechanism to objectively determine models that are of practical significance.

The strategies based on Bayesian inference, and more specifically, using the minimum message length (MML) framework have increasingly found support in mixture modelling tasks (Wallace and Boulton, 1968;Wallace, 1986;Roberts et al, 1998;Figueiredo and Jain, 2002;Kasarapu and Allison, 2015). The MML-based inference framework decomposes a modelling problem into two parts: the first part determines the model complexity by encoding all the parameters of the model, and the second part corresponds to encoding of the observed data using the chosen parameters. Thus, a two-part message length is obtained for a model under consideration. A model that results in the least total message length is then determined to be the optimal model under this framework (Oliver and Baxter, 1994).
The MML-based search method developed by Kasarapu and Allison (2015) is demonstrated to outperform the traditionally used approaches and is the current state-of-the art. Kasarapu and Allison (2015) have designed the mixture modelling apparatus to include Gaussian distributions to model data in the Euclidean space and von Mises-Fisher (vMF) distributions to model directional data distributed on the surface of a sphere. While Gaussian mixtures are ubiquitously used because of their computational tractability (McLachlan and Peel, 2000), they are ineffective to model directional data. In this context, analogues of the Gaussian distribution defined on the surfaces of the appropriate Reimannian manifolds are typically considered. The vMF is the most fundamental directional probability distribution defined on the spherical surface. It is the spherical analogue of a symmetrical Gaussian wrapped around a unit hypersphere (Fisher, 1953;Watson and Williams, 1956) and is demonstrated to be useful in large-scale text clustering (Banerjee et al, 2003;Kasarapu and Allison, 2015) and gene expression analyses (Banerjee et al, 2005). A general form of the vMF distribution is the Fisher-Bingham (FB 5 ) distribution which is used to model asymmetrically distributed data on the spherical surface (Kent, 1982). Mixtures of FB 5 distributions have been employed by Peel et al (2001) to identify joint sets in rock masses, and by Hamelryck et al (2006) to sample random protein conformations. The FB 5 distribution has increasingly found support in machine learning tasks for structural bioinformatics (Kent and Hamelryck, 2005;Boomsma et al, 2006;Hamelryck, 2009).
A two-dimensional version of the vMF distributions called the von Mises circular is used to model data distributed on the boundary of a circle. Each data point on the circle has a domain [−π, π). If such data occur as pairs, then the resulting manifold in three dimensions would be a torus. The bivariate von Mises (BVM) distributions are used to model such data distributed on the toroidal surface and serve as the Gaussian analogue. Motivated by its practical applications in bioinformatics, the BVM distributions are widely studied. The mixtures of BVM distributions have been previously used in modelling protein dihedral angles (Dowe et al, 1996;Mardia et al, 2007Mardia et al, , 2008. However, these approaches have some limitations. Dowe et al (1996) treat the pairs of angles to be independent of each other and do not account for their correlation. This is akin to conflating two von Mises circular distributions together. As explained in Section 4, such an approximation leads to inefficient mixtures. Although Mardia et al (2007) use BVM distributions that account for the correlation between the angular pairs, they do not have a rigorous search method to determine the optimal mixtures. This limits their ability to correctly distinguish among models that, while being of different type, have the same number of model parameters.
This paper develops the mixture modelling apparatus to address these limitations using the MML framework. Further, different variants of the BVM distribution obtained by constraining some of its characterizing parameters (see Section 2). Mardia et al (2007) have evaluated the utility of these variants in the context of modelling the protein dihedral angles. We adopt the MML principle in objectively assessing the mixture distributions of these variants (see Section 4). We have developed a search method to determine the optimal number of mixture components and their parameters that describe the given data in a completely unsupervised setting. The use of the MML modelling paradigm and our proposed search method is explored on real-world data corresponding to the dihedral, that is, torsion angles of protein structures. We demonstrate that mixtures of BVM distributions facilitate the design of reliable computational models for protein structural data.
In addition to determining the optimal number of mixture components, the parameters of the individual component distributions need to be estimated. Traditionally, the optimum parameters are obtained by maximum likelihood (ML) or Bayesian maximum a posteriori probability (MAP) estimation. For a mixture distribution, the parameters are estimated by maximizing the likelihood of the data by employing an expectationmaximization (EM) algorithm that iteratively updates the mixture parameters (Dempster et al, 1977). The key differences between ML, MAP and MML-based estimation is: (1) unlike ML, MML uses a prior over the parameters and considers their precision while encoding; (2) unlike MAP, MML estimators are invariant under non-linear transformations of the parameters (Oliver and Baxter, 1994). The estimation of parameters using ML ignores the cost of stating the parameters, and MAP based estimation uses the probability density of parameters instead of their probability measure. In contrast, the MML inference process takes into account the optimal precision to which parameters should be stated and uses it to determine a corresponding probability value. Parameter estimation using the MML framework has been carried out on various probability distributions (Wallace, 2005). Kasarapu and Allison (2015) have demonstrated that the MML estimators outperform the traditionally used estimators in the case of Gaussian and vMF distributions. Furthermore, for a FB 5 distribution, Kasarapu (2015) have shown that the MML estimators have lower bias and error as compared to the ML and MAP estimators. Contributions: The main contributions of this paper are as follows: -We derive the MML-based estimates of the parameters of a BVM distribution. The MML estimators are demonstrated to have lower bias and mean squared error when compared to their traditional counterparts. We consider two variants of the BVM distribution, namely the Independent (Dowe et al, 1996) and the Sine variant (Singh et al, 2002). -We design a search method to infer the optimal number of BVM mixture components that best describe the angular data distributed on the toroidal surface. -We demonstrate the utility of the MML framework in determining the suitability of the two variants of the BVM distribution in modelling the protein dihedral angle data. We show that the Sine variant that includes the correlation term explains the data much more effectively than the Independent version. -We demonstrate the effectiveness of the mixture modelling method by applying it to cluster protein dihedral angles. We demonstrate that the resulting mixtures closely correspond to the commonly observed secondary structural regions in protein structures.
The rest of the chapter is organized as follows: Section 2 describes the BVM distribution, the Independent and the Sine variants, and their relevance in modelling data distributed on the toroidal surface. Section 3 describes the MML framework and outlines the differences between parameter estimation using ML, MAP and MML methods. It includes the derivation of the MML estimators of the parameters of the BVM distribution. We empirically demonstrate that the MML estimators outperform the traditionally used ML and MAP estimators by having lower bias and mean squared error. Section 4 discusses the search and inference of mixtures of bivariate von Mises (BVM) distributions using the MML framework. As a specific application, we employ the mixtures to model protein dihedral angles. We demonstrate that our search method is able to infer meaningful clusters that directly correspond to frequently occuring conformations in protein structures.

Bivariate von Mises probability distribution
The class of bivariate von Mises (BVM) distributions was introduced by Mardia (1975a,b) to model data distributed on the surface of a 3D torus. The study of these distributions has been partly motivated by biological research, where it is required to model the protein dihedral angles (see Section 4.2). The probability density function of the BVM distribution has the general form f (x; Θ) ∝ exp{κ 1 cos(θ 1 − µ 1 ) + κ 2 cos(θ 2 − µ 2 ) + (cos θ 1 , sin θ 1 )A(cos θ 2 , sin θ 2 ) T } (1) where x = (θ 1 , θ 2 ), such that θ 1 , θ 2 ∈ [−π, π) and the parameter vector Θ = (µ 1 , µ 2 , κ 1 , κ 2 , A), such that µ 1 , µ 2 ∈ [−π, π) are the mean angles, κ 1 ≥ 0 and κ 2 ≥ 0 are the concentration parameters, and A is a 2 × 2 real-valued matrix. The term exp{κ 1 cos(θ 1 − µ 1 )} corresponds to a von Mises distribution on a circle characterized by the parameters µ 1 and κ 1 , Hence, the BVM distribution (Equation 1) can be explained as a product of two von Mises circular distributions, with an additional exponential term involving A, that accounts for the correlation. The general form of the BVM distribution has 8 free parameters. In order to draw an analogy to the bivariate Gaussian distribution (with 5 free parameters), sub-models of the BVM distribution have been proposed by restricting the values that A can take (Jupp and Mardia, 1980). A 6-parameter version was explored by Rivest (1988) and has the form In particular, when α = 0 and β = λ, the above density reduces to the following 5-parameter version, which is called the BVM Sine model (Singh et al, 2002).
where c(κ 1 , κ 2 , λ) is the normalization constant of the distribution defined as and I v is the modified Bessel function of first kind and order v. The 5-parameter vector will be Θ = (µ 1 , µ 2 , κ 1 , κ 2 , λ) where λ is a real number. If λ = 0, the probability density function (Equation 3) will just be the product of two independent von Mises circular distributions, and corresponds to the case when there is no correlation between the two variables θ 1 and θ 2 . The probability density function in such a case is given as where c(κ 1 , κ 2 ) is the normalization constant defined as c(κ 1 , κ 2 ) = 1 2πI 0 (κ 1 ) . and corresponds to the product of the normalization constants for the respective von Mises circular distributions. Alternatively, when α = −β, the form of Equation 2 results in a different reduced form called the BVM Cosine model (Mardia et al, 2007). The Sine and the Cosine models serve as natural analogues of the bivariate Gaussian distribution on the 3D torus. In fact, for huge concentrations, Singh et al (2002) approximate the Sine model to a bivariate Gaussian distribution with the 2 × 2 covariance matrix C = [c ij ], i, j ∈ {1, 2}, whose elements are given by The limiting case approximation is valid when κ 1 κ 2 > λ 2 . Also, from the covariance matrix, the correlation coefficient ρ can be determined as (Pearson, 1895): In order to better understand the interaction of κ 1 , κ 2 , and the correlation coefficient ρ, we provide an example in Figure 1, where the distribution is shown for values of ρ = 0.1 (low correlation), ρ = 0.5 (moderate correlation), and ρ = 0.9 (high correlation). Note that ρ can take negative values, in which case the resultant distribution will just be a reflection in some axis (Mardia et al, 2007). The modelling of directional data using the BVM Sine and Cosine models has been previously explored by Mardia et al (2007). For estimating the parameters of the distribution, maximum likelihood based optimization is used. We discuss ML and MAP based estimators, which are the traditionally used methods of parameter estimation.

Maximum likelihood parameter estimation
In applications involving modelling directional data using the BVM Sine distributions, the maximum likelihood (ML) estimates are typically used (Boomsma et al, 2006;Mardia et al, 2007Mardia et al, , 2008. For BVM Sine distributions, the moment and ML estimates are the same, as the BVM Sine distribution belongs to the exponential family of distributions (Mardia et al, 2008).
, the ML estimates of the parameter vector Θ = (µ 1 , µ 2 , κ 1 , κ 2 , λ) are obtained by minimizing the negative log-likelihood expression of the data given by The ML estimates satisfy ∂L ∂Θ = 0. However, as no closed form soultions exist because of the complicated form of c(κ 1 , κ 2 , λ), an optimization library is used. We use NLopt 1 , a non-linear optimization library, to compute the parameter estimates.

Maximum a posteriori probability (MAP) estimation
For an independent and identically distributed sample D, the MAP estimates are obtained by maximizing the posterior density Pr(Θ|D). This requires the definition of a reasonable prior Pr(Θ) on the parameter space. The MAP estimates are sensitive to the nature of parameterization of the probability distribution and this limitation is discussed here. We demonstrate that the MAP estimators are inconsistent and are subjective to the parameterization. We consider two alternative parameterizations in the case of the BVM Sine distribution.

Non-linear transformations of the parameter space
We consider non-linear transformations of the parameter space, in order to demonstrate that the MAP estimates are not invariant in different parameterizations of the probability distribution. We discuss a simple non-linear transformation of the parameter space involving the correlation parameter λ. Additionally, we also describe a parameterization that transforms all the five parameters.
An alternative parameterization involving λ: The BVM Sine probability density function (Equation 3) can be reparameterized in terms of the correlation coefficient ρ, instead of λ, by using the relationship λ = ρ √ κ 1 κ 2 (as per Equation 6). If Θ = (µ 1 , µ 2 , κ 1 , κ 2 , ρ) denotes the modified vector of parameters, the modified prior density Pr(Θ ) is obtained by dividing Pr(Θ) with the Jacobian of the transformation J = ∂ρ ∂λ With this transformation, the posterior density Pr(Θ |D) can be computed, and subsequently used to determine the MAP estimates.
An alternative parameterization involving Θ: In addition to the transformation of the correlation parameter λ, we study another transformation that was proposed by Rosenblatt (1952). The method transforms a given continuous k-variate probability distribution into the uniform distribution on the k-dimensional unit hypercube. Such a transformation applied on the prior density of the parameter vector Θ results in the prior transforming to a uniform distribution. Hence, estimation in this transformed parameter space is equivalent to the corresponding maximum likelihood estimation. For the 5-parameter vector Θ = (µ 1 , µ 2 , κ 1 , κ 2 , λ), the Rosenblatt (1952) As the cumulative densities are bounded by 1, the above transformation results in 0 ≤ z i ≤ 1, i = {1, . . . , 5}. Further, Rosenblatt (1952) argue that each z i is uniformly and independently distributed on [0, 1], so that the prior density in this transformed parameter space is Pr(Θ ) = Pr(z 1 , z 2 , z 3 , z 4 , z 5 ) = 1 (10) In order to achieve such a transformation, we need to express z i in terms of the original parameters. Based on the assumptions made in the formulation of the prior Pr(Θ), we derive the following relationships: Based on the independence assumption in the formulation of priors of angular and scale parameters, we have z 3 = F 3 (κ 1 |µ 2 , µ 1 ) = F 3 (κ 1 ), and therefore we have Hence, κ 1 = tan(arccos(1 − z 3 )) and κ 2 = tan(arccos(1 − z 4 )) Further, F 5 (λ|κ 2 , κ 1 , µ 2 , µ 1 ) = F 5 (λ|κ 2 , κ 1 ), as λ is independent of µ 1 and µ 2 . Hence, the invertible transformation corresponding to λ is as follows so that λ can be expressed as a function of z 3 , z 4 , and z 5 . The transformed BVM Sine probability density function f (x, Θ ) is obtained by substituting the expressions of Θ in terms of In summary, we considered two additional parameterizations of the BVM Sine probability density. For statistical invariance, the estimates of the parameters should also be affected by the same transformation in alternative parameterizations. The MAP estimation does not satisy this property, as illustrated by the following example.

An example demonstrating the effects of alternative parameterizations
An example of estimating parameters using the posterior distributions resulting from the various prior densities (Equations 8 -10) is described here. A random sample of size N = 10 is generated from a BVM Sine distribution (Singh et al, 2002). The true parameters of the distribution are µ 1 = µ 2 = π/2, κ 1 = κ 2 = 10, and λ = 9 (corresponding to a correlation coefficient of ρ = 0.9).
The MAP estimators are obtained by maximizing the posterior densities using the non-linear optimization library NLopt (Johnson, 2014) in conjunction with derivative-free optimization (Powell, 1994). The differences in the estimates are explained below.

Minimum Message Length (MML) Inference
In this section, we describe the model selection paradigm using the Minimum Message Length criterion and proceed to give an overview of MML-based parameter estimation for any distribution.
3.1 Model selection using minimum message length criterion Wallace and Boulton (1968) developed the first practical criterion for model selection based on information theory. As per Bayes's theorem: where D denotes observed data, and H some hypothesis about that data. Further, Pr(H & D) is the joint probability of data D and hypothesis H, Pr(H) and Pr(D) are the prior probabilities of hypothesis H and data D respectively, Pr(H|D) is the posterior probability, and Pr(D|H) is the likelihood.
As per Shannon (1948), given an event E with probability Pr(E), the length of the optimal lossless code to represent that event requires I(E) = − log 2 (Pr(E)) bits. Applying Shannon's insight to Bayes's theorem, Wallace and Boulton (1968) got the following relationship between conditional probabilities in terms of optimal message lengths: The above equation can be intrepreted as the total cost to encode a message comprising of the following two parts: 1. First part: the hypothesis H, which takes I(H) bits, 2. Second part: the observed data D using knowledge of H, which takes I(D|H) bits.
As a result, given two competing hypotheses H and H , Hence, Pr(H |D) = 2 ∆I Pr(H|D) gives the log-odds posterior ratio between the two hypotheses. The framework provides a rigorous means to objectively compare two competing hypotheses. Clearly, the message length can vary depending on the complexity of H and how well it can explain D. A more complex H may explain D better but takes more bits to be stated itself. The trade-off comes from the fact that (hypothetically) transmitting the message requires the encoding of both the hypothesis and the data given the hypothesis, that is, the model complexity I(H) and the goodness of fit I(D|H).

MML-based parameter estimation
Wallace and Freeman (1987) introduced a generalized framework to estimate a set of parameters Θ given data D. The method requires a reasonable prior h(Θ) on the hypothesis and evaluating the determinant of the Fisher information matrix |F(Θ)| of the expected second-order partial derivatives of the negative loglikelihood function, L(D|Θ). The parameter vector Θ that minimizes the message length expression (given by Equation 11) is the MML estimate according to Wallace and Freeman (1987).
where d is the number of free parameters in the model, and q d is the d-dimensional lattice quantization constant (Conway and Sloane, 1984). The total message length I(Θ, D), therefore, comprises of two parts: (1) the cost of encoding the parameters, I(Θ), and (2) the cost of encoding the data given the parameters, I(D|Θ). A concise description of the MML method is presented in Oliver and Baxter (1994).
The key differences between ML, MAP, and MML estimation techniques are as follows: in ML estimation, the encoding cost of parameters is, in effect, considered constant, and minimizing the message length corresponds to minimizing the negative log-likelihood of the data (the second part). In MAP based estimation, a probability density rather than the probability is used. It is self evident that continuous parameter values can only be stated to some finite precision; MML incorporates this in the framework by determining the region of uncertainty in which the parameter is located.
gives a measure of the volume of the region of uncertainty in which the parameter Θ is centered. This multiplied by the probability density h(Θ) gives the probability of a particular Θ as Pr(Θ) = h(Θ)V . This probability is used to compute the message length associated with encoding the continuous valued parameters (to a finite precision).

MML estimation of the parameters of the BVM distribution
In this section, we outline the derivation of the MML-based parameter estimates of a BVM Sine distribution. As explained in Section 3.2, the derivation of the MML estimates requires the formulation of the message length expression (Equation 11) for encoding some observed data using the BVM Sine distribution.
The formulation requires the use of a suitable prior density on the parameters. We use the parameterization Θ and the corresponding prior Pr(Θ) that was formulated in the MAP analyses in Section 2.2. It is to be noted that the MML estimation is invariant to the parameterization used (Oliver and Baxter, 1994).
Notations: Before describing the MML approach, the following notations are defined as these are used in the following discussion. The partial derivatives of the normalization constant c(κ 1 , κ 2 , λ) of the BVM Sine distribution would be required later on. The following are the notations adopted to represent them.
We also require the determinant of the Fisher information for the MML estimation of parameters. We use the above notations in the following computation of the Fisher information. The computation of these partial derivatives is explained in Section 3.4.

Computation of Expectations
In order to proceed with the derivation of the Fisher information, we first outline the derivation of some of the required expectation quantities. For random variables θ 1 , θ 2 sampled from the BVM Sine distribution (Equation 3), we compute the following quantities: Singh et al (2002) derived the normalization constant as an infinite series expansion given by Equation 4. We use the following integral form of the normalization constant to derive the above mentioned expectations, as a function of κ 1 , κ 2 , and λ.
On differentiating the above integral with respect to κ 1 , we get Hence, the expectation can be represented using the above defined notation as On differentiating twice the integral form of c(κ 1 , κ 2 , λ) with respect to κ 1 , κ 2 , and λ, we get the following relationships

Computation of the Fisher information
As described in Section 3.2, the computation of the determinant of the Fisher information matrix requires the evaluation of the second order partial derivatives of the negative log-likelihood function with respect to the parameters of the distribution. As per the density function (Equation 3), the negative log-likelihood of a datum x = (θ 1 , θ 2 ) is given by where Θ = (µ 1 , µ 2 , κ 1 , κ 2 , λ) as indicated before. Let F 1 (Θ) denote the Fisher information for a single observation. the Fisher information matrix F 1 (Θ) in the case of an FB 5 distribution is a 5 × 5 symmetric matrix. Further, the determinant |F 1 (Θ)| is decomposed as a product of |F A | and |F S |, where F A is the Fisher matrix associated with the angular parameters µ 1 and µ 2 , and F S is the Fisher matrix associated with the scale parameters κ 1 , κ 2 , and λ.
Fisher matrix (F A ) associated with µ 1 , µ 2 : F A is a 2 × 2 symmetric matrix whose elements are the expected values of the second order partial derivatives of L with respect to µ 1 and µ 2 . On differentiating Equation 14 with respect to µ 1 , we get and On taking the derivative of Equation 15 with respect to µ 2 , we get Fisher matrix (F S ) associated with κ 1 , κ 2 , λ: F S is a 3 × 3 symmetric matrix whose elements are the expected values of the second order partial derivatives of L with respect to κ 1 , κ 2 , and λ. On differentiating Equation 14 with respect to κ 1 , κ 2 , and λ, we get Fisher matrix F(Θ) associated with the 5-parameter vector Θ: On differentiating Equation 15 with respect to κ 1 and computing the expectation of the differential, we get Hence, E ∂ 2 L ∂κ 1 ∂µ 1 = 0 = F κ 1 µ 1 and E ∂ 2 L ∂λ∂µ 1 = 0 = F λµ 1 This allows for the computation of |F 1 (Θ)| as the product of |F A | and |F S |, that is, Then, the Fisher information for some observed data D = {x 1 , . . . , x N } is given by as each element in |F 1 (Θ)| is multiplied by the sample size N .

Message length formulation
The message length to encode some observed data D can now be formulated by substituting the prior density Pr(Θ) (Equation 8), the Fisher information |F(Θ)| and the negative log-likelihood of the data (Equation 7) in the message length expression (Equation 11). The MML parameter estimates are the ones that minimize the total message length. As there is no analytical form of the MML estimates, the solution is obtained, as for the maximum likelihood and MAP cases, by using the NLopt optimization library (Johnson, 2014). At each stage of the optimization routine, the Fisher information needs to be calculated. However, this involves the computation of complex entities such as the normalization constant c(κ, β) and its partial derivatives. The computation of these intricate mathematical forms using numerical methods is discussed next in Section 3.4.

Computation of the normalization constant and its derivatives
The computation of the negative log-likelihood and the message length requires the normalization constant and its associated derivatives. In this section, the description of the methods that can be employed to efficiently compute these complex functions is explored. We will utilize the properties of Bessel functions to implement the normalization constant and the necessary partial derivatives as limiting order summations for the BVM Sine distribution.
Computation of the series S (m,n) 1 : We first establish that f j+1 < f j ∀j ≥ 0 and show that S (m,n) 1 converges to a finite sum as j → ∞. Consider the logarithm of the ratio of consecutive terms f j and f j+1 in S for p, v > 0, I p+v < I p , and the ratio I p+v I p → 0 as p → ∞ (Amos, 1974). Further, e < 1 implies the above equation is the sum of negative terms. Hence, log f j+1 is a convergent series.
For a practical implementation of the sum, we need to express S (m,n) 1 as the modified summation S (m,n) 1 where each f j is divided by the maximum term f 0 . For each j > 0, log f j is calculated using the previous term log f j−1 (Equation 21). The new term t j = f j /f 0 is then computed 2 as exp(log f j − log f 0 ). This is because computing the difference with the maximum value and then taking the exponent ensures numerical stability.
The summation is terminated when the ratio t j j k=0 t k < (a small threshold ∼ 10 −6 ). -Let the j th term dependent on κ 1 in Equation 4 be represented as g j (κ 1 ) = I j /κ j 1 , where I j implicitly refers to I j (κ 1 ). Based on the relationship between the Bessel functions I j , I j+1 , and the derivative I j in Equation 23 (Abramowitz and Stegun, 1965), the expressions for the first and second derivatives of g j (κ 1 ) (Equation 24) are derived as -Let S(c κ 1 ) = log c κ 1 : Because of the similar forms of g j (κ 1 ) and g j (κ 1 ), the expression for S(c κ 1 ) will be similar to that of S(c) with a change in order of the Bessel functions from m = 0 in Equation 20 to m = 1. Hence, S(c κ 1 ) = S where m, n ∈ {0, 1}, δ 2 = 8π 2 λ , and e = λ 2 4κ 1 κ 2 . Note that S (m,n) 2 is a convergent series (the proof is based on the same reasoning as for S (m,n) 1 ).
Let the j th term dependent on λ, κ 1 in Equation 4 be represented as g j (λ, κ 1 ) = λ 2j I j κ j 1 . Its partial derivatives are given below. These derivatives are the terms in the series S (m,n) 2 (after factoring out the common elements as δ 2 ). -The expression to compute S(c λλ ) = log c λλ is given by The practical implementation of S (m,n) 2 and S(c λλ ) is similar to that of S (m,n) 1 given by Equation 22. However, in these cases, the expressions of f j and consequently t j , are modified depending on their specific forms. Also, the series begin from j = 1 and, hence, the respective maximum terms will correspond to f 1 .

Evaluation of the MML estimates
For a given BVM Sine distribution characterized by concentration parameters κ 1 , κ 2 and correlation coefficient ρ, a random sample of size N is generated using the method proposed by Mardia et al (2007). The angular parameters of the true distribution are set to {µ 1 , µ 2 } = π/2. The scale parameters κ 1 , κ 2 , and ρ are varied to obtain different BVM Sine distributions and corresponding random samples. The parameters are estimated using the sampled data and the different estimation methods. The procedure is repeated 1000 times for each combination of N, κ 1 , κ 2 , and ρ.

Methods of comparison
For every randomly generated sample from a BVM Sine distribution, we compute the the ML, MAP, and MML estimators of the parameters, and these are compared with each other across all the simulations. The results include the three versions of MAP estimates resulting from the three forms of the posterior distributions (Equations 8-10): MAP1 corresponds to the posterior with parameterization Θ = (µ 1 , µ 2 , κ 1 , κ 2 , λ), MAP2 corresponds to the posterior with parameterization Θ = (µ 1 , µ 2 , κ 1 , κ 2 , ρ), and MAP3 corresponds to the posterior with parameterization Θ = (z 1 , z 2 , z 3 , z 4 , z 5 ). As noted in Section 2.2, the MAP3 estimator will be the same as the ML estimator due to the Rosenblatt (1952) transformation of Θ to Θ .
In order to compare the various estimators, we use the mean squared error (MSE) and Kullback-Leibler (KL) distance as the objective evaluation metrics. The estimates are also compared using statistical hypothesis testing. For a parameter vector Θ characterizing a true BVM Sine distribution, and its estimate Θ, we analyze the MSE and KL distance of Θ with respect to the true parameter vector Θ. The analytical form of the KL distance between two BVM distributions is derived in Appendix A. We analyze the percentage of times (wins) the KL distance of a particular estimator is smaller than that of others. When the KL distance of different estimates is compared, because of three different versions of MAP estimation, three separate frequency plots are presented. corresponding to the MAP1, MAP2, and MAP3 estimators.
With respect to statistical hypothesis testing, the likelihood ratio test statistic is asymptotically approximated as an χ 2 distribution with five degrees of freedom For the various parameter estimates compared here, it is expected that at especially large sample sizes, the estimates are close to the ML estimate. In other words, the empirically determined test statistic is expected to be lower than the critical value τ = 13.086, corresponding to a p-value greater than 0.01.

Empirical analyses
As per the experimental setup, we present the results for when the original distribution from which the data is sampled has κ 1 = 1 and κ 2 = 10. The correlation coefficient ρ is varied between 0 and 1, so that we obtain different values for the correlation parameter λ (Equation 6). We discuss the results for varying values of sample sizes N , and ρ = 0.1, 0.5, 0.9, corresponding to a low, moderate, and high correlation, respectively.
For ρ = 0.1 : ρ = 0.1 : ρ = 0.1 : The results are presented in Figure 2. Compared to the ML estimators, the MAP and MML estimators result in lower bias and MSE for all values of N . Both the bias and MSE continue to decrease as the sample size increases, as the estimation improves with more evidence for all methods. When compared with MAP1 and MAP2, the MML estimators have greater bias and greater MSE. As with the FB 5 distribution, we observe that that MAP1 and MAP2 result in different estimators, and therefore, result in different bias and MSE values.
The KL distance with respect to MAP1 is in favour of the MAP1 estimators. The MAP1 estimates result in lower KL distance as compared to the other estimators almost 50% of the 1000 simulations for each N (Figure 2c). However, the MML estimators win when the MAP2 and MAP3 versions are used. When MAP3 is used, the MML estimators have a smaller KL distance in close to 70% of the simulations (Figure 2e). Further analysis using statistical hypothesis testing illustrates that the null hypotheses corresponding to the MAP and MML estimators are accepted (p-values greater than 0.01 in Figure 2f). at the 1% significance level.
For ρ = 0.5 : ρ = 0.5 : ρ = 0.5 : Similar to when ρ = 0.1, we observe that the bias and MSE of the MAP and MML estimators are lower than the ML estimators for different values of N . In contrast to ρ = 0.1, the bias of the MML estimator is lower than the MAP1 estimator but higher than the MAP2 estimator ( Figure 3a). As with the previous case, MAP-based estimation result in different estimators. Further analysis of the estimators using KL distance and statistical hypothesis testing follow the same pattern as when ρ = 0.1.
For ρ = 0.9 : ρ = 0.9 : ρ = 0.9 : The results are presented in Figure 4. As with the previous two cases, we observe that the ML estimators have the greatest bias and MSE for all values of N . The bias of the MML estimators is lower than all the MAP estimators. However, the MSE of the MML estimators is greater compared to the MAP1 or MAP2 estimators. Contrary to the previous two cases, we observe that the frequency of wins of KL distance for the MML estimators is lower when compared to MAP2 estimation (Figure 4e). Further, the results following the statistical hypothesis testing follow the same trend as the previous two cases. As the same size increases, the different estimators converge to the ML estimators as seen from the high p-values (Figure 4g).
The empirical analyses of the controlled experiments discussed above indicate that the ML estimators of the parameters of a BVM distribution are biased. The same was observed with other directional probability distributions such as the vMF (Kasarapu and Allison, 2015) and FB 5 (Kasarapu, 2015). Also, we observe that the MAP estimation method result in different estimators depending on how the distribution is parameterized. We have shown that the MAP estimators are not invariant under non-linear transformations of the parameter space. In this context, the MML estimators are empirically demonstrated to have lower bias than the traditional ML estimators and are invariant to alternative parameterizations unlike the MAP estimators.

Mixtures of bivariate von Mises distributions
We consider two kinds of bivariate von Mises (BVM) distributions in mixture modelling. In addition to the Sine variant (Equation 3) that has the correlation parameter λ, we also consider the independent variant obtained when λ = 0. The independent version assumes zero correlation between the data distributed on the torus (see Equation 5). We provide a comparison for the mixture models obtained using both versions of the BVM distributions.
Previous work on MML-based modelling of protein dihedral angles used independent BVM distributions (Dowe et al, 1996). Their work used the Snob mixture modelling software (Wallace and Dowe, 1994b). As pointed out by Dowe et al (1996), Snob does not have the functionality to account for the correlation between the data. We therefore study the BVM Sine distributions and demonstrate how they can be integrated with our generalized MML-based mixture modelling method.

Approach for BVM distributions
We extend the search method described in Kasarapu and Allison (2015) to infer mixtures of BVM distributions.
To infer the optimal number of mixture components, the mixture modelling apparatus is now modified to handle the directional data distributed on the surface of a torus. As in the case of the vMF and FB 5 distributions, the split operation detailed in Kasarapu (2015) is tailored for the BVM mixtures. The basic idea behind splitting a parent component is to identify the means of the child components so that they are on either side of the parent mean and are reasonably apart from each other. Recall that for a Gaussian parent component, we computed the direction of maximum variance and selected the initial means, along this direction, that are one standard deviation away on either side of the parent mean. We employ the same strategy for BVM distributions. For , we compute the direction of maximum variance in the (φ, ψ)-space. This allows us to compute the initial means of the child components.
The delete and merge operations are carried out in the same spirit. During merging BVM components, the KL distance is evaluated to determine the closest pair. We derive the KL distance for BVM Sine and BVM Independent distributions as shown in Appendix A. Further, in all the operations, the MML estimators of the BVM Sine distribution, derived in Section 3.3 are used in the update step of the EM algorithm.

Mixture modelling of protein main chain dihedral angles
We consider the spatial orientations resulting from the interactions of the main chain atoms in protein structures. A protein main chain is comprised of a chain of amino acids, each os which is characterized by a central carbon C α . The angular data corresponds to the spatial orientations of the the planes containing the atoms from successive amino acids. A protein main chain is characterized by a sequence of φ, ψ, and ω angles. These angles uniquely determine the geometry of the protein backbone structure (Richardson, 1981). However, in a majority of protein structures, ω = 180 • and, hence the sequence of C α -C-N-C α atoms lie in a plane (see the dotted planar representation in Figure 5a). As a result, the angles φ and ψ are typically analyzed (Ramachandran et al, 1963). The angles φ and ψ are called the dihedral angle pair corresponding to an amino acid residue with a central carbon atom C α along the protein main chain. Geometrically, a dihedral angle is the angle between any two planes defined using four non-collinear points. In Figure 5(a), φ is the angle between the two planes formed by C-N-C α and N-C α -C. Similarly, ψ is the angle between the two planes formed by N-C α -C and C α -C-N. The dihedral angles φ and ψ are measured in a consistent manner. For example, in order to measure φ, the four atoms C-N-C α -C are arranged such that φ is calculated as the deviation between N-C and C α -C when viewed in some consistent orientation. As an illustration, in Figure 5(b), view the arrangement of the four atoms through the N-C α bond such that C α is behind the plane of the paper and N eclipses the C α atom. Also, the C atom directly attached to N is at the 12 o' clock position. In this orientation, φ is given as the angle of rotation required to align the N-C bond with the C α -C bond in the plane of the paper. Further, if it is a clockwise rotation, it is considered a positive value. This ensures that φ ∈ [−π, π). The dihedral angle ψ is measured by following the same convention with the four atoms being N-C α -C-N. The (φ, ψ) pair measured in this way can be plotted on the surface of a 3D torus. Each (φ, ψ) pair corresponds to a point on the toroidal surface. The angle φ is used to identify a particular cross-section (circle) of a torus, while ψ locates a point on this circle (see Figure 6).
We generate the entire set of dihedral angle data from the 1802 experimentally determined protein structures in the ASTRAL SCOP-40 (version 1.75) database (Murzin et al, 1995) representing the "β class" proteins. The number of (φ, ψ) dihedral angle pairs resulting from this data set is 253,165. We model this generated set of dihedral angles using BVM Sine distributions.
A random sample from this empirical distribution consisting of 10,000 points is shown in Figure 7. The plot is a heat map showing the density of the data distribution on the toroidal surface. Note that there are regions on the torus which are highly concentrated (yellow), corresponding to the helical regions in the protein. The ellipse-like patches (mostly in blue) roughly correspond to the β strands in proteins. Furthermore, the data is multimodal which motivates its modelling using mixtures of BVM distributions. We consider the effects of using the BVM Sine distribution as compared to the BVM Independent variant in this context. Fig. 7 A sample of 10,000 points randomly generated from the empirical distribution of (φ, ψ) pairs. The figure shows the random sample from different viewpoints.

Search of BVM Independent and BVM Sine mixtures
The search method inferred a 32-component BVM Independent mixture and terminated after 42 iterations involving split, delete, and merge operations. In the case of modelling using BVM Sine distributions, our search method inferred 21 components and terminated after 29 iterations. In each of these iterations, for every intermediate K-component mixture, each constituent component is split, deleted, and merged (with an appropriate component) to generate improved mixtures.
The progression of the search method for the optimal BVM Independent mixture begins with a single component. The search method results in continuous split operations until the 17 th iteration when a 17component mixture is inferred (see Figure 8a). This corresponds to a progressive increase in the first part of the message (red curve). Between the 17 th and the 21 st iterations, we observe a series of delete/merge and split operations leading to a stable 19-component mixture. The search method again continues to favour the split operations until the 28 th iteration when a 26-component mixture is inferred. Thereafter, a series of deletions and splits yield a stable 29-component mixture at the end of the 35 th iteration. The search method eventually terminates when a 32-component mixture is inferred with a characteristic step-like behaviour towards the end indicating perturbations involving split and delete/merge operations (see Figure 8a).
In the case of searching for the optimal BVM Sine mixture, our proposed search method continues to split the components thereby increasing the mixture size. This occurs until 21 iterations. At this stage, there are 21 mixture components. This can be observed in Figure 8(b), when the first part of the message (red curve) continually increases until the 21 st iteration. During this period, observe that the second part (blue) and the total message length (green) continually decrease signifying an improvement to the mixtures.
After the 21 st iteration, we observe a step-like behaviour as in the case of mixture modelling using the BVM Independent distributions. The behaviour characterizes the reduction or increase in the number of mixture components corresponding to a decrease or increase to the first part of the message. After the 24 th iteration, we observe that the mixture has 22 components. However, the final mixture stabilizes in the subsequent iterations to a 21-component mixture. After the 29 th iteration, there is no further improvement to the total message length and the search method terminates. We observe a characteristic increase in the mixture size initially followed by some perturbations that stabilize the intermediate mixture (step-like behaviour), eventually resulting in an optimal mixture (see Figure 8). There is an initial sharp decrease in the total message length until about 7 iterations for BVM mixtures. Because of the multimodal nature of the directional data (see Figure 7), the initial increase in the number of components would explain the data distribution corresponding to those modes that are clearly distinguishable. This leads to a substantial improvement to the total message length as the minimal increase in the first part is dominated by the gain in the second part. However, towards the end of the search, when the increase in first part dominates the reduction in second part, the method stops. Thus, we see the trade-off of model complexity (as a function of the number of components and their parameters), and the goodness-of-fit being balanced using the search based on the MML inference framework.

Comparison of BVM mixture models of protein data
The existing work of MML-based mixture modelling of protein dihedral angles by Dowe et al (1996) inferred 27 clusters using the BVM Independent distributions. In contrast, our search method inferred 32 clusters. However, their data consists of only 41,731 (φ, ψ) pairs generated from the protein structures known at that time. In contrast, we have used 253,165 pairs of dihedral angles along with a different search method as explained previously (see Section 4.2). So, there is some consensus on the rough number of component distributions if the protein dihedral angles were modelled using BVM distributions assuming no correlation between φ and ψ.
The visualization of the dihedral angles is commonly done by the Ramachandran plot (Ramachandran et al, 1963) who first analyzed the various possible protein configurations and represented them as a two-dimensional plot. An example of one such plot is provided in Lovell et al (2003) and reproduced here (Figure 9). Such a plot is indicative of the allowed conformations that protein structures can adopt. There are vast spaces in the dihedral angle space where few data are present. The conformations corresponding to those regions are not possible. We consider the plot to explain the similarities between our inferred mixture models and the one that is traditionally used.
Our resulting mixtures of BVM Independent and the Sine variants are shown in Figure 10. The contours of the constituent components plotted in the (φ, ψ)-space can be seen in the diagram. For visualization purposes, we display the contour of each component that corresponds to 80% of the data distribution. The data in Figure 10 corresponds to a random sample drawn from the empirical distribution (same as in Figure 7) visualized in the (φ, ψ)-space.
In Figure 9, we observe that the top-left region corresponds to the β strands in protein structures. The empirical distribution of dihedral angles we generated also has this characterstic. We observe a concentrated mass in the top-left in Figure 10. Furthermore, our inferred mixtures are able to model this region using the appropriate components. Note that smaller or highly compact contours correspond to BVM distributions that have greater concentration parameters (κ 1 and κ 2 in Equation 3).
We note that components numbered 1-11 ( Figure 10a) and components 1-8 ( Figure 10b) are used to describe this region. These correspond to the components of the BVM Independent and BVM Sine models respectively. Clearly, more number of components are required to model roughly the same amount of data (corresponding to the β strands) using the BVM Independent mixture.
Similarly, in Figure 9, we observe another concentration of mass in the middle-left portion of the figure. This corresponds mainly to right-handed α-helices, which are very frequent in protein structures. In Figure 10, we have the corresponding mass and also note the dense region (bright yellow). As per our inferred mixtures, component 17 ( Figure 10a) and component 12 (Figure 10b) are used to predominantly describe this dense region. The other surrounding regions in the dihedral angle space of the right-handed helices are described by components 12-19 ( Figure 10a) and by components 9-13 ( Figure 10b). Again, we observe that the similar data is described using 8 components by the BVM Independent mixture as opposed to 5 components by the BVM Sine mixture. Lovell et al (2003) display another region of concentrated mass in the middle-right of Figure 9. This region corresponds to the infrequent left-handed helices in protein structures. We see a corresponding mass in the empirical distribution in Figure 10. The components 20-25 of our inferred BVM Indenpendent mixture describe this region (Figure 10a). The same region is described by components 14-16 of our inferred BVM Sine mixture (Figure 10b). Notice how this region is described by components 15 and 16. These two components describe the dense mass within this region while component 14 is responsible for mainly modelling the data that is further away from this clustered mass. We again observe that the same region is modelled by greater number of components when using BVM Independent distributions.
The remaining mixture components describe the insignificant mass present in other regions of the dihedral angle space. The ability of our inferred mixtures to identify and describe specific regions of the protein conformational space in a completely unsupervised setting is remarkable. Further, we have qualified the effects of using the BVM distributions which do not account for the correlation between the dihedral angle pairs. In this regard, the BVM Sine mixtures fare better when compared to mixtures of BVM Independent distributions. We now quantify these effects in terms of the total message length.
Our proposed search method to infer an optimal mixture involves evaluating the encoding cost of the mixture parameters or the first part (model complexity), and encoding the data using those parameters or the second part (goodness-of-fit). The progression of the search method continues until there is no improvement to the total message length. We observe that the resulting 21-component BVM Sine mixture has a first part of 966 bits and a corresponding second part of 5.735 million bits (see Table 1). A BVM Independent mixture with the same number of components has a first part of 872 bits and a corresponding second part of 5.751 million bits. Although the model complexity is lower for the BVM Independent mixture (difference of ∼ 94 bits), the BVM Sine mixture has an additional compression of ∼ 16, 000 bits in its goodness-of-fit. Thus, the significant gain in the second part dominates the minimal increase in the first part of the BVM Sine mixture.
Further, if we compare the 21-component BVM Independent mixture with the inferred 32-component BVM Independent mixture, we observe that the first part is more in the 32-component case. This is expected because there are more number of mixture parameters to encode in the 32-component mixture. There is a difference of 1292 − 872 = 420 bits (see Table 1). However, the 32-component mixture results in an extra compression of ∼ 15, 000 bits. So, the total message length is lower for the 32-component mixture, and is therefore, preferred to the 21-component BVM Independent mixture. When the inferred 32-component BVM Independent and the 21-component BVM Sine mixtures are compared, we observe that the total message length is lower for the BVM Sine mixture. In this case, both the first and second parts are lower for the Sine mixture leading to an overall gain of about ∼ 1000 bits. Thus, the BVM Sine mixture is more appropriate as compared to the BVM Independent mixture in describing the protein dihedral angles. This exercise shows how an optimal mixture model is selected by achieving a balance between the trade-off due to the complexity and the goodness-of-fit to the data. Furthermore, as in the case of the vMF and FB 5 distributions, we can devise null model descriptions of protein dihedral angles based on the BVM mixtures. For comparison, we consider a uniform distribution on the torus, which is referred to as the uniform null model in the equation below.
where R and r are the radii that define the size of the torus (see Figure 6). When R = r = 1, the surface area of the torus is 1/4π 2 . The null models based on the BVM mixtures have the same form as the vMF and FB 5 mixtures given as mixture distributions (Kasarapu, 2015) with the number of respective components being K = 32 and K = 21 corresponding to the Independent and the Sine variants respectively. Compared to the uniform model, both the BVM mixtures result in additional compression (see Table 2). The message length to encode the entire collection of 253,165 dihedral angle pairs using the uniform null model is 6.388 million bits which amounts to 25.234 bits per residue. In comparison, the BVM Independent mixture results in a compression of 5.735 million bits which amounts to 22.656 bits per residue. The additional compression is therefore, close to 2.58 bits per residue (on average). The BVM Sine mixture further leads to an additional compression of 323 bits over the BVM Independent mixture. This is equivalent to an additional saving of 0.0013 bits per residue (on average). These results indicate that the BVM mixtures are superior compared to the uniform model. This can be argued from the fact that the empirical distribution (see Figure 7) has empty regions in the dihedral angle space. This is also confirmed from the Ramachandran plot (Figure 9). However, the BVM Independent and the BVM Sine variants are in close competition with each other. Noting that we need more mixture components in the Independent case and because the Sine mixture can describe the data more effectively, we conclude that the BVM Sine mixture supersedes the BVM Independent mixture. The ability of the BVM Sine mixture to model correlated data leads to improved description of the protein dihedral angles.

Conclusion
We have considered the problem of modelling directional data using the bivariate von Mises distributions. We have demonstrated that the MML-based estimation results in parameters that have a lower bias and MSE compared to the traditional ML estimators, and contrast to MAP estimators, they are invariant to transformations of the parameter space. To model empirically distributed data with multiple modes, we have used mixtures of BVM distributions. We have addressed the important problems of selecting optimal number of mixture components along with their parameters using the MML inference framework. We employed the designed framework to model protein dihedral angles using mixtures of BVM distributions. The empirical distribution of the pairs of dihedral angles represented on a toroidal surface clearly suggests correlation between the angle pairs. As such, the BVM Sine mixtures are shown to be appropriate. Both the BVM Independent and the Sine mixtures effectively model the dihedral angle space. The ability of the search method to correctly identify components corresponding to the regions of critical protein configurations is remarkable. This is more so because our search method does not rely on any prior information and infers the mixtures in a completely unsupervised setting.