
==== Front
Adv Data Anal Classif
Adv Data Anal Classif
Advances in Data Analysis and Classification
1862-5347
1862-5355
Springer Berlin Heidelberg Berlin/Heidelberg

558
10.1007/s11634-023-00558-2
Regular Article
Parsimony and parameter estimation for mixtures of multivariate leptokurtic-normal distributions
http://orcid.org/0000-0003-4543-0218
Browne Ryan P. 1
http://orcid.org/0000-0003-3129-5385
Bagnato Luca luca.bagnato@unicatt.it

2
http://orcid.org/0000-0001-7742-1821
Punzo Antonio 3
1 https://ror.org/01aff2v68 grid.46078.3d 0000 0000 8644 1405 Department of Statistics and Actuarial Science, University of Waterloo, Waterloo, ON Canada
2 https://ror.org/03h7r5v07 grid.8142.f 0000 0001 0941 3192 Department of Economic and Social Sciences, Catholic University of the Sacred Heart, Milano, Italy
3 https://ror.org/03a64bh57 grid.8158.4 0000 0004 1757 1969 Department of Economics and Business, University of Catania, Catania, Italy
27 9 2023
27 9 2023
2024
18 3 597625
30 11 2022
2 7 2023
© The Author(s) 2023
2023
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Mixtures of multivariate leptokurtic-normal distributions have been recently introduced in the clustering literature based on mixtures of elliptical heavy-tailed distributions. They have the advantage of having parameters directly related to the moments of practical interest. We derive two estimation procedures for these mixtures. The first one is based on the majorization-minimization algorithm, while the second is based on a fixed point approximation. Moreover, we introduce parsimonious forms of the considered mixtures and we use the illustrated estimation procedures to fit them. We use simulated and real data sets to investigate various aspects of the proposed models and algorithms.

Supplementary Information

The online version contains supplementary material available at 10.1007/s11634-023-00558-2.

Keywords

Leptokurtic-normal distribution
Majorization–minimization algorithm
Mixture models
Parsimony
Mathematics Subject Classification

62Hxx Multivariate analysis
62H30 Classification and discrimination; cluster analysis (statistical aspects); mixture models
62H12 Estimation in multivariate analysis
62G07 Density estimation
MUR2022XRHT8R Browne Ryan P. Next Generation EUPE0000018 Bagnato Luca http://dx.doi.org/10.13039/501100000038 Natural Sciences and Engineering Research Council of Canada RGPIN-2018-04444 Punzo Antonio Università Cattolica del Sacro CuoreOpen access funding provided by Università Cattolica del Sacro Cuore within the CRUI-CARE Agreement.

issue-copyright-statement© The Managing Editorial Board 2024
==== Body
pmcIntroduction

Finite mixtures of distributions are a powerful clustering tool (Titterington et al. 1985; Bouveyron et al. 2019; McNicholas 2016a). In this family, for a d-variate continuous random vector X, mixtures of multivariate normal distributions have received considerable attention because of their computational and theoretical convenience by assuming, in most cases, that each mixture component represents a cluster (or group) within the original data (see, e.g., Fraley and Raftery 2002 and McNicholas 2016b). The popularity of these mixtures has been further increased since the works of Banfield and Raftery (1993) and Celeux and Govaert (1995), that introduce parsimony via the eigendecomposition of the component covariance matrices. This produces a set of 14 parsimonious models that can handle various clustering scenarios and are easy to interpret.

However, in many real data applications, the tails of the component normal distributions are shorter than required. This can have a negative effect on the underlying clustering structure and also on the estimates of the component parameters (Peel and McLachlan 2000). For this reason, several mixtures of elliptical heavy-tailed distributions have been proposed in the literature. Examples are mixtures of the following multivariate distributions: t (Peel and McLachlan 2000; Andrews and McNicholas 2012), contaminated normal (Punzo and McNicholas 2016; Punzo et al. 2018, 2020), leptokurtic-normal (Bagnato et al. 2017, 2022), power exponential (Zhang and Liang 2010 and Dang et al. 2015), tail-inflated normal (Punzo and Bagnato 2021; Tomarchio et al. 2022) and shifted exponential normal (Punzo and Bagnato 2020 and Tomarchio et al. 2022). Among them, mixtures of multivariate leptokurtic-normal (MLN) distributions (refer to Sect. 2) have the advantage to be characterized by parameters directly related to the moments of practical interest, especially in each cluster.

For maximum likelihood (ML) estimation of the parameters of MLN mixtures, Bagnato et al. (2017) used the EM algorithm; in particular, they adopted the general-purpose optimizer optim(), provided in R, to update the component-wise parameters in each M-step of the algorithm. Browne (2022) introduced a fixed point (FP) algorithm and showed improvements in terms of computational time. However, the FP algorithm is not guaranteed to be monotonic although Browne (2022) provided some evidence of this through simulations.

In Sect. 4, we explore parsimonious forms of MLN mixtures by setting some parameters to be equal across components. For ML estimation we derive two algorithms, one based on the MM principle (Lange et al. 2000; Hunter and Lange 2004), in Sect. 5, and the second one based on the FP algorithm from Browne (2022), in Sect. 6. The MM algorithm is based on a quadratic function obtained by using a bound on the second derivative of the log-density. We derive this bound and show connection between the MM and FP procedures in Sect. 3.

The remaining part of the paper is structured as follows. In Sect. 7, we first use simulated data to compare the two proposed algorithms in terms of computational burden and solutions obtained (Sect. 7.1). Then, in Sect. 7.2, we use real data with the same aim as for Sect. 7.1, but also to show the existence of real data situations where a parsimonious variant of the general model can provide a better description of the data. In Sect. 8, we describe the results of a further simulation study aiming to investigate aspects such as: the evaluation of the asymptotic properties of the ML estimator (Section E.1 in the supplementary material), model selection via the BIC (Sect. 8.1), and the behavior of the parsimonious MLN mixtures under different data generating models (Sect. 8.2). We give a final overview of the results obtained in this paper in Sect. 9.

A supplementary material is also available online where we provide: an overview about the notation and the main functions used herein (Section A), useful details for the implementation of the proposed approaches (Sections B and C), and some tables, figures, and details related to the simulation studies discussed in Sects. 7.1 (Section D) and 8 (Section E). Section F of the supplementary material provides an additional simulation study to investigate the viability of the estimation procedures on large data sets and sampling variability of the estimation procedure.

Finite mixtures of MLN distributions

The MLN distribution—with mean vector μ, covariance matrix Σ, and excess kurtosis β—was introduced by Bagnato et al. (2017) and has the following joint probability density function (pdf)1 fx∣μ,Σ,β=1+βgx-μ⊤Σ-1x-μϕx∣μ,Σ,

where ϕx∣μ,Σ denotes joint pdf of the multivariate normal (MN) distribution with mean μ and covariance matrix Σ,2 gr=r2-2(d+2)r+d(d+2)8d(d+2),

andβ∈0,4d(d+2)d+4.

The MLN distribution is obtained by applying a multivariate Gram-Charlier expansion to the MN distribution and has one additional parameter, β, that, as already said above, controls the excess kurtosis.

For model-based clustering purposes, Bagnato et al. (2017) also proposed and illustrated finite mixtures of MLN distributions with joint pdf3 px∣θ=∑j=1kπjfx∣μj,Σj,βj,

where πj>0, ∑j=1kπj=1, and θ contains all parameters.

Bagnato et al. (2017) proposed (1) to model data with excess kurtosis and (3) as a robust model based clustering approach. Browne (2022) showed that finite mixtures of MLN distributions are identifiability when (μj,Σj) are different; he also suggested that, if μj,Σj=μt,Σt for some j≠t, then these two components are to be combined with β⋆=πjβj+πtβt.

ML estimation: surrogate functions for MM and FP procedures

Herein, we derive the key functions for the MM and FP procedures used to ML estimates of the parameters of the MLN distribution. First, we find a bound for the Hessian. To derive a bound we begin by rewriting the pdf in (1) as a function of both the standardized vector z=Σ-1/2x-μ and the squared Mahalanobis distance r=z⊤z=x-μ⊤Σ-1x-μ as4 ur:=logfr=logfz⊤z=-d2log2π-logΣ+log1+βgr-r2.

Browne (2022) used this form to show that the log-density is unimodal by finding a bound on the first derivative,u′r:=ddrlogfr=βg′r1+βgr-12.

Here, we expand on this by finding a bound on the Hessian. Using the chain rule we can obtain the Hessian and we findd2logf(z)dzdz⊤=d2logf(r)dr4zz⊤+2dlogf(r)drId⪯d2logf(r)dr4r+2dlogf(r)drId,

where Id denotes the d×d identity matrix and A⪯B means that B-A is positive definite. This bound depends on the following univariate function5 hr:=d2logfrdr4r+2dlogfrdr=βg′′r1+βgr4r-βg′r1+βgr24r+2βg′r1+βgr.

We can then find the maximum and minimum of h(·) by taking the derivative with respect to r. Setting the derivative equal to zero we find that the stationary points are roots to the following depressed quartic polynomial-3d(d+2)2(β+8)(βd+4β-8d)+8β(d+2)2(βd-β+12d)r+6β(d+2)(β-4d)r2+β2r4=0.

Then, maximum and minimum are obtained by evaluating h(·) at the four roots. These roots and bounds are illustrated in Fig. 1 for two different values of d (2 and 5).

Of particular interest is the lower bound M, as it can be used to construct a quadratic minorizer function. We can obtain it simply evaluating h(·) at the four roots, say r1, r2, r3, and r4, i.e.6 M=|minh(r1),h(r2),h(r3),h(r4)|.

Note, we apply the absolute value so that the quadratic minorizer has a particular form. Also, the value M only depends on β and d.Fig. 1 Examples of hr function along with the stationary points and bounds with β=4d(d+2)/(d+4). Left panel has d=2 and right panel has d=5

Applying this bound, we can construct the following quadratic minorizer, using the expansion point z0,7 logfz⊤z≥logfz0⊤z0+∇logfz⊤zz-z0-M2z-z0⊤z-z0≥C0+2u′z0⊤z0+Mz0⊤z-M2z⊤z,

where C0 depends on z0. This quadratic minorizer is used as a starting point for the MM algorithm presented in Sect. 6.

Alternatively, the following quadratic surrogate function can be constructed8 logfz⊤z≈logfz0⊤z0+ddrlogfrr=z0⊤z0z⊤z-z0⊤z0≈logfz0⊤z0+u′z0⊤z0z⊤z-z0⊤z0.

This approximation was considered by Browne (2022) and is used as a starting point for the FP algorithm presented in Sect. 6. Browne (2022) showed that even though this surrogate function is not minorizer, an algorithm can be constructed to yield monotonic log-likelihood values. Monotonicity of the log-likelihood values was shown through simulations.

Figure 2 illustrates the functions given in (7) and (8), when d=1 and β=4d(d+2)/(d+4) when using expansion points equal to z0=1,2,3.Fig. 2 Left panel: The log density, logf(z2), along with three quadratic minorizers based on (7) using the expansion points z0=1,2,3. Right panel: The log density, logf(z2), along with three quadratic surrogate functions based on (8) using the expansion points z0=1,2,3. Both figures are based on d=1 and β=4d(d+2)/(d+4)

Parsimonious finite mixtures

Here, we introduce parsimonious forms of the finite mixture model in (3). Similar to Fraley and Raftery (2002) and Celeux and Govaert (1995), we consider the following 3-term eigendecomposition of the covariance matrix9 Σj=λjΓjΨjΓj⊤,

where λj=|Ψj|1/d, Γj is a d×d orthogonal matrix of the eigenvectors of Σj, Ψj is the d×d diagonal matrix with the scaled eigenvalues of Σj (such that Ψj=1) located on the main diagonal. The decomposition in (9) has some useful practical interpretations. From a geometric point of view, λj determines the volume, Γj governs the orientation, and Ψj denotes the shape of the jth cluster. From a statistical point of view, as well-documented in Fraley and Raftery (2002), Celeux and Govaert (1995), Greselin and Punzo (2013), Bagnato and Punzo (2021), and Punzo and Bagnato (2022), the columns of Γj govern the orientation of the principal components (PCs) of the jth cluster, the diagonal elements in Ψj are the normalized variances of these PCs, and λj can be meant as the overall volume of the scatter in the space spanned by the PCs of the jth cluster.

By setting the three parameters be equal, vary or fixed across groups (See Celeux and Govaert 1995), the decomposition in (9) yields fourteen different models. These model are denoted by three letters (Fraley et al. 2015; Scrucca et al. 2016; Pocuca et al. 2022). The MLN model has an additional parameter βj which can be equal or vary across groups. This results in twenty eight possible models. This family of models is represented using a 4-letter identifier of the type X1X1X3X4, where: X1=E if λj=λ or X1=V if λj is varying across groups;

X2=I if Ψj=Id, X2=E if Ψj=Ψ or X2=V if Ψj is varying across groups;

X3=I if Γj=Id, X3=E if Γj=Γ or X3=V if Γj is varying across groups;

X4=E if βj=β or X4=V if βj is varying across groups.

Just as an example, VEIE indicates Σj=λjΨ and βj=β. However, we will also use a 3-letter identifier X1X1X3 when βj, j=1,…,d, will be considered as fixed.

ML estimation

Similar to Bagnato et al. (2017) and Browne (2022), we apply the expectation-maximization (EM) algorithm (Dempster et al. 1977) to obtain ML estimates. In the context of finite mixture models, we introduce the component memberships as missing data for each observation xi in the sample. We represent the missing data using an indicator variable, denoted by wij, such that wij=1 if observation i was generated by component j and wij=0 otherwise. At iteration q, the expected value of the missing data given the observed data xi and the current parameter values (at iteration q) iswij(q)=πj(q)fxi∣μj(q),Σj(q),βj(q)∑h=1kπh(q)fxi∣μh(q),Σh(q),βh(q).

Using these expected values, we can form the expected complete log-likelihood10 C+∑i=1n∑j=1kwij(q)logπj-rij2+log1+βjgrij,

where rij=xi-μj⊤Σj-1xi-μj and C does not depend on the parameters.

The M-step involves maximizing the expected complete log-likelihood in (10). The update for πj is πj(q+1)=nj(q)/n, where nj(q)=∑i=1nwij(q). The update for βj, when equal and varying across components, is given in Sect. 4.2. Then, for (μ1,Σ1,…,μk,Σk) we consider two alternative methods for estimation; a MM algorithm and a FP algorithm.

The MM algorithm is based the minorizer derived in (7), while the FP algorithm is based on the surrogate function in (8).

Section 5 gives an overview of the MM algorithm but all the technical derivations for each model are quite long, so we provide the details in Section B of the supplementary material. The FP algorithm updates are shown in Sect. 6. The FP algorithm have form related to Celeux and Govaert (1995) and as such we can use formulae similar to Celeux and Govaert (1995) and Browne and McNicholas (2014a) to construct updates for all the eigendecomposed models.

Updating βj or β

Varying βj

For the parameter βj, Bagnato et al. (2017) maximizes the weighted log-likelihood with respect to a transformed variable ηj=log[βj/(βmax-βj)], where βmax is the largest value βj can be. However, Browne (2022) showed that the log-likelihood (10), as a function of βj, is self concordant. Self concordant functions yield a bound on the number of Newton steps required to minimize a function (for more details see Boyd and Vandenberghe 2006).

An unconstrained Newton step for βj isβj(q+1)=βj(q)+∑i=1nwij(q)γiμj(q),Σj(q),βj(q)∑i=1nwij(q)γiμj(q),Σj(q),βj(q)2,

whereγiμ,Σ,β=βgxi-μ⊤Σ-1xi-μ1+βgxi-μ⊤Σ-1xi-μ.

We apply this update while maintaining that βj stays within the range of unimodality. Alternatively, one could apply an univariate numerical optimization or incorporate backtracking with Newton’s method.

Equal βj=β

Let βj be equal across components, i.e., βj=β. Since self concordance is preserved under summation and affine transformation, the log-likelihood is again concordant. Update for β, based on an unconstrained Newton step, is given byβ(q+1)=β(q)+∑i=1n∑j=1kwij(q)γiμj(q),Σj(q),β(q)∑i=1n∑j=1kwij(q)γiμj(q),Σj(q),β(q)2.

MM algorithm for updating μj and Σj

Applying βj to (6) we obtain Mj and then using (7), we can obtain the following minorizer functionC-nj(q)2logΣj+TrΣj-1/2∑i=1nwij(q)2u′rij(q)+Mjxi-μjzij(q)⊤-Mj2TrΣj-1∑i=1nwij(q)xi-μjxi-μj⊤,

where C does not depend on the parameters, rij(q)=zij(q)⊤zij(q) and zij(q)=Σj(q)-1/2xi-μj(q). Taking the derivative, and solving, yields the following update for μjμj(q+1)=x¯j(q)-1Mjnj(q)∑i=1nwij(q)2u′rij(q)+MjΣj(q)1/2zij(q)⊤,

wherex¯j(q)=1nj(q)∑i=1nwij(q)xi.

The minorizer, in terms of the covariance parameters, with updated μj(q+1), is11 ∑j=1k-nj(q)2logΣj+TrΣj-1/2Aj-12TrΣj-1Bj,

whereAj=∑i=1nwij(q)2u′rij(q)+Mjxi-μj(q+1)zij(q)⊤

andBj=Mj∑i=1nwij(q)xi-μj(q+1)xi-μj(q+1)⊤.

Equation (11) is the starting point for all the eigendecomposed models.

The technical derivations for each model are quite long, so we provide the details for each model in Section B of the supplementary material. When considering a model, we repeatedly use reparameterization and variety of methods to obtain updates for the parameters. For example, for the VVV model we let Ξj=Σj-1/2 then use the algebraic Riccati equation (for review, see Wonham 1968) to obtain an unique symmetric positive definite update based method from Laub (1979). Estimation of Γ with the constraint Γ⊤Γ=Id has been considered by Browne and McNicholas (2014b) and Absil et al. (2008), but here we will follow Browne and McNicholas (2014a) and Kiers (2002).

FP algorithm for updating μj and Σj

When we are in the generic jth component, using (8), an alternative surrogate function for the expected complete log-likelihood (10) isC-nj(q)2logΣj-M2TrΣj-1∑i=1nwij(q)κiμj(q),Σj(q),βj(q)xi-μjxi-μj⊤,

whereκiμ,Σ,β=1-2βg′xi-μ⊤Σ-1xi-μ1+βgxi-μ⊤Σ-1xi-μ.

Note that κi(·) is bounded to the interval [0, 1]. Taking the derivative, and solving, yields the following update for μj12 μj(q+1)=∑i=1nwij(q)κiμj(q),Σj(q),βj(q)xi∑i=1nwij(q)κiμj(q),Σj(q),βj(q).

Then, using the updated μj, the update for Σj is based on the following quantityRj(q)=1nj(q)∑i=1nwij(q)κiμj(q+1),Σj(q),βj(q)xi-μj(q+1)xi-μj(q+1)⊤.

Browne (2022) noticed using Rj(q) to update Σj can cause fluctuations in the log-likelihood values; so, instead, the author suggests using a weighted average13 Sj(q)=1-a(q)Rj(q)+a(q)Rj(q)wherea(q)=11+β(q)d+44d(d+2),

so that the weight is 1 when β(q)=0 and the weight is 1/2 when β(q) is equal to the largest possible value of 4d(d+2)/(d+4). In Section C of the supplementary material we explore different values for the weight.

We use S1(q),…,Sk(q) and follow Celeux and Govaert (1995) to update all the eigendecomposed models, except for VVE and EVE where we follow Browne and McNicholas (2014a). Just as an example, for the VVV (or varying Σj) model, the update is Σj(q+1)=Sj(q+1), and for the EEE (Σj=Σ) model the update isΣ(q+1)=∑j=1kπj(q+1)Sj(q).

Using the updates (12) and (13) does not guarantee an increase in the log-likelihood. However, Browne (2022) noted that the sequence of log-likelihoods generated from the EM exhibited monotonicity and no fluctuations in the log-likelihood values. Section C of the supplementary material illustrates fluctuations in the sequence of likelihood values when a=1 and β=4d(d+2)/(d+4 and investigates a range of weights. Similar to Browne (2022), when using (12) and (13), the log-likelihood sequences obtained from the simulation studies in Sects. 7 and 8 exhibited monotonicity. This provides additional evidence towards a monotonicity property for this FP procedure.

Since the updating procedure for the weighted log-likelihood can be viewed as a generalized EM algorithm, convergence criteria, such as the Aitken acceleration (see Böhning et al. 1994 and Lindsay 1995), can be applied.

Analyses on simulated and real data

Using the R package Browne et al. (2023), we simulated and real data to compare the MM and FP algorithms in terms of computational time and solutions obtained when using: (1) the random starting values and (2) the Aitken acceleration (Aitken 1926, Böhning et al. 1994 and Lindsay 1995) as stopping criterion, using a tolerance equal to 10-6. The computational time is calculated using the wall-clock time on an iMac with 4 GHz Quad-Core Intel Core i7. Real data are also used to show the importance of considering parsimonious variants of the general, unconstrained, MLN mixture model (VVVV; Sect. 7.2). Always in this section, we adopt the Bayesian information criterion (BIC; Schwarz 1978) for model selection and the adjusted Rand index (ARI; Hubert and Arabie 1985) for clustering/classification assessment. Note that, in the formulation used herein, the BIC values are multiplied by -1; this means that, in analogy with the log-likelihood, this criterion needs to be maximized.

Simulated data: comparing MM and FP algorithms

Here, we compare the MM and FP algorithms on simulated data sets while varying the number of variables and covariance model. We hold βj fixed and equal to 4d(d+2)/(d+4) since both procedures use the same updates for this parameter. So, we only focus on the updates for μj and Σj, j=1,…,k. We assess the average computational time for a fixed number of iterations, the total time to convergence, and the number of iterations at convergence. We generate n=500 observations from a two component (k=2) MLN mixture model with parameters equal to14 (π1,π2)=12,12,μ1,μ2=-51d,+51d,Σ1=Σ2=Id,andβ1=β2=4d(d+2)(d+4).

Both algorithms use the same initial value, which is generated randomly. We varied the number of variables d to be 2, 8, or 16. We compare the computational time, and the log-likelihood values at convergence, from fitting a k=2 component MLN mixture model. We used 100 replications to reduce simulation error, i.e., differences due to a particular data set.

As an initial comparison, we evaluate the computational time for each procedure when performing 100 EM iterations. Table 1 reports the ratio of computational times (MM over FP) and indicates the FP algorithm is faster than the numerical optimization procedures except for the models with covariance equal to VVE and EVE. This is due to the FP algorithm requiring the eigenvalues from a positive definite matrix whereas the MM algorithm does not.Table 1 Ratio of computational times (MM over FP), for 100 iterations, while varying the number of variables d and the covariance structure

d	EEE	VVV	VEE	EVV	EEV	VEV	EVI	EII	VEI	VVI	EEI	VII	VVE	EVE	
2	1.6	1.6	1.5	1.2	1.4	1.2	1.4	1.3	1.3	1.3	1.2	1.2	0.7	0.8	
8	1.5	1.5	1.4	1.3	1.4	1.3	1.5	1.4	1.4	1.4	1.4	1.4	0.8	0.8	
16	2.0	1.9	1.8	1.8	1.7	1.6	1.5	1.4	1.4	1.4	1.4	1.4	1.0	0.8	

The second comparison quantifies the computational time, the total number of iterations and the log-likelihood values at convergence. Table 2 displays the average computational time and iterations from the FP and MM algorithms when d=2 or 8. The results in the case d=16 are given in Section D of the supplementary material. The tables also show the mean difference in the log-likelihoods, say l^max-l^FP and l^max-l^MM, where l^max=maxl^FP,l^MM and l^FP and l^MM are the log-likelihood values obtained at convergence for the FP and MM algorithms.

When the procedures converge, they tend to have roughly the same number of iterations and log-likelihood values. However, for some models, the MM algorithm reaches the maximum number of iterations. This indicates we might want to use a stricter convergence criterion for the MM procedure. Interestingly, EVE and VVE models had slower computational time for 100 iterations, but here we see that requiring fewer iterations yields a faster time to convergence.Table 2 Average computational time and number of iterations until convergence (with the maximum number of iterations fixed to 1000), for each estimation algorithm over 100 replications

	d=2, Mean from (FP, MM) of	d=8, Mean from (FP, MM) of	
Time	Iterations	l^max-l^	Time	Iterations	l^max-l^	
EVV	0.012	1.085	16	980	0.0	1.7	0.051	1.437	55	1000	0.0	28.5	
VEV	0.015	1.029	19	950	0.0	1.5	0.072	1.427	74	1000	0.0	27.6	
EEV	0.012	0.969	16	930	0.0	0.1	0.074	1.378	80	1000	0.0	3.6	
EVE	0.034	1.002	26	928	0.0	1.7	0.908	1.381	505	1000	0.0	26.2	
VVE	0.036	0.894	27	874	0.0	0.0	0.872	1.337	487	1000	0.0	1.7	
VVV	0.015	0.015	20	15	0.0	0.0	0.039	0.061	44	50	0.0	0.0	
EEE	0.010	0.011	14	12	0.0	0.0	0.019	0.019	22	17	0.0	0.0	
VVI	0.011	0.010	17	12	0.0	0.0	0.011	0.015	15	15	0.0	0.0	
EEI	0.009	0.009	13	11	0.0	0.0	0.009	0.011	11	11	0.0	0.0	
VEE	0.012	0.013	17	11	0.0	1.1	0.025	0.013	28	9	0.0	34.5	
VII	0.010	0.010	14	11	0.0	0.0	0.006	0.008	8	8	0.0	0.0	
VEI	0.011	0.010	16	11	0.0	0.5	0.010	0.008	12	8	0.0	3.6	
EVI	0.009	0.009	13	10	0.0	1.0	0.011	0.007	14	6	0.0	7.0	
EII	0.007	0.008	10	10	0.0	0.0	0.006	0.007	7	6	0.0	0.0	
The columns with header l^max-l^ show the average difference between the log-likelihood at convergence and the maximum log-likelihood for that data set. The output is displayed in pairs (FP, MM), which is the result from the FP and MM algorithms

Real data: parsimony and computational efficiency

Here, we use real data with a twofold aim. Firstly, we show the existence of real-data situations where a parsimonious MLN mixture can provide better performance than its unconstrained version. Secondly, we highlight the gain of the two fitting algorithms here introduced, in terms of computation time, with respect to the one illustrated in Bagnato et al. (2017). We consider the ais (Australian Institute of Sports) data set (Cook and Weisberg 1994) already considered by Bagnato et al. (2017) as an example of a real data situation where the MLN mixture is to be preferred to well-known alternatives in the mixture of elliptical heavy-tailed distributions world; for other applications of this benchmark data set in the model-based clustering literature see, e.g., Azzalini et al. (2016), Dang et al. (2017), Galimberti and Soffritti (2014), and Morris et al. (2019).

The data set contains measurements on n=202 athletes, subdivided in k=2 groups (100 female and 102 male), and is available in the R-packages alr3 (Weisberg 2011) and sn (Azzalini 2015). As in Bagnato et al. (2017), we analyze a subset of d=6 variables: height in centimeters (Ht), weight in kilograms (Wt), lean body mass (LBM), red cell count (RCC), white cell count (WCC), and Hematocrit (Hc). We show the scatter plot of the data, with labelling based on the gender, in Fig. 3.Fig. 3 ais data. Scatter plot (M denotes male and F female)

On these data, we fit two-component MLN mixtures with four different approaches:EM algorithm where a BFGS algorithm, with exact derivatives, is used in each M-step to update all the component parameters;

EM algorithm where a BFGS algorithm, with numerical derivatives, is used in each M-step to update all the component parameters;

MM algorithm;

FP algorithm.

The first two approaches are those considered by Bagnato et al. (2017) for the VVVV model only, while the remaining ones are those proposed herein to fit all the 28 parsimonious models.

For the VVVV model only, Table 3 compares the first two fitting algorithms in terms of Log-likelihood, BIC, ARI, and ratio of the elapsed computational time (over the FP) at convergence. As we can see, they provide the same results, with the exception of the computational time which is much higher (about 690 times higher than the time required by the FP algorithm to fit the same model) when numerical derivatives are used. Anyway, even for the EM with exact derivatives, the computational time required is about 77 times higher than the one required by the FP algorithm.Table 3 Log-likelihood, BIC, ARI, and ratio of the elapsed computational time (over the FP) at convergence of the algorithms adopted by Bagnato et al. (2017) to fit the VVVV MLN mixture

Fitting algorithm	Log-lik	BIC	ARI	Time (over FP)	
EM with exact derivatives	-2777.461	-5857.493	0.811	76.732	
EM with numerical derivatives	-2777.461	-5857.493	0.811	690.049	

For all the 28 parsimonious MLN mixtures, Table 4 compares the MM and FP algorithms in terms of log-likelihood values and computational time required. As for the first comparison we note that, for some models, the two algorithms provide different results. A background gray color is used to highlight the algorithm working better. Here we realize that the maximum log-likelihood value provided by the FP algorithm is always greater than, or equal to, the one from to the MM algorithm. This is in line with the results obtained in Sect. 7.1. In terms of computational times, instead, the competition is more heated; however, 15 out of 28 times, the FP algorithm is faster in providing the final solution, and this further corroborates the results in Sect. 7.1. Because of the results on the log-likelihood, Table 3 also reports BIC and ARI values for the solution from the FP algorithm only. For the sake of interpretation, we also report the ranking, among parsimonious model, induced by these criteria.Table 4 Log-likelihood values and ratio of elapsed computational times (MM over FP) at convergence of the FP and MM algorithms for all the parsimonious MLN mixtures

A background gray color for the best log-likelihood values separately for each model. For the models estimated by the FP algorithm, BIC and ARI values are reported along with the corresponding rankings

The first aspect we note comparing Table 3 with Table 4 is that, for the VVVV model, all the competing fitting algorithms provide the same results, but the MM and FP algorithms, which behave similarly in terms of computational time (ratio of 1.286), obtain these results in a much more fast way. As for the parsimony, the model selected by the BIC (EEVV) provides exactly the same partition of the data given by the VVVV model, but this partition in now obtained via a simpler model with a gain both in terms of estimated parameters (6 parameters less) and interpretation: thanks to the selected model, we realize the two inferred groups are the same in terms volume and shape, but differ in terms of orientation (take a look at the scatterplot between Wt and LBM in Fig. 3) and different excess kurtosis. Finally, by looking at the ARI values in Table 4, a better classification of the data could be obtained by considering the still more parsimonious VEEE and EEEE models, which make the ARI to increase to 0.884.

To find out more about the EEVV model selected by the BIC, we used the parametric bootstrap technique with 10000 bootstrapped samples. The estimates of the parameters along with standard errors (in round brackets) computed by using a parametric bootstrap technique with 10000 bootstrapped samples are showed in Table G.8 shown in Section G of the supplementary material. To tackle the well-known label switching issue in each bootstrap replication, we assigned the labels by minimizing the distance between the estimated and true means, as in Stephens (2000), Farcomeni and Punzo (2020), Gallaugher et al. (2022), and Tomarchio et al. (2022), just to cite a few. The excess kurtosis estimates highlight a further interesting aspect; while one group needs a heavy-tailed distribution to be described in a good way (with a model-based estimated excess kurtosis of 17.083), the other one is well-fitted by a classical normal distribution.

Investigating some aspects of the MLN mixture

We further investigate the properties of the MLN mixture model when using only the FP method of estimation (Browne et al. 2023). We consider three simulation studies to quantify the sampling variability of the estimation procedure (given in Section E.1 of the supplementary material), the ability of the BIC to pick the correct model (Sect. 8.1), and the effect when data is generated from the t or the generalized hyperbolic distributions (in Sect. 8.2).

For these simulations we vary the number of components (k), the number of variables (d), the distance between the components (δ), the covariance structure (governed by ρ and λ, with ρ∈[0,1) and λ>0), and the sample size (n). We fix the β parameter to its maximum value, which depends on d, in all the simulations. We provide the model parameters, θkδ,λ,ρ,d, in a compact form below; for an expanded form, see Section E of the supplementary material. The covariance or scale structure is defined by15 Ξdρ=1-ρId+ρJd,

where Jd denotes a d×d matrix of ones. The two-component setup is16 θ2δ,λ,ρ,d=πj=12,μj=(3-2j)δ1d,Σj=λ3-2jΞd(ρ),βj=4d(d+2)(d+4),

Table 5 Connection between the covariance model structure and the parameters λ and ρ

Parameters	EIIE	VIIE	EEEE	VEEE	
(λ,ρ)	(1, 0)	(2, 0)	(1, 0.5)	(2, 0.5)	

where j=1,2. The three-component setup is17 θ3δ,λ,ρ,d=πj=13,μj=j-2δ1d,Σj=λ(j-2)Ξdρ,βj=4d(d+2)(d+4),

where j=1,2,3. The four-component setup is18 θ4δ,λ,ρ,d=πj=14,μj=δ(-1)j1d(-1)j+121dΣj=λj/2-1Ξd(ρ),βj=4d(d+2)(d+4),

where j=1,2,3,4 and ⌊·⌋ is the floor function, i.e. ⌊1.6⌋=1 and ⌊-0.6⌋=-1.

This setup allows to use λ and ρ to move through the eigendecomposed or covariance model space. The particular parameter values are given within the forthcoming subsections but of special interest are models given in Table 5. The table shows the connection between some of the covariance models and the parameters.

Figure 4 shows three examples of 500 observations generated from θkδ=5,λ=2,ρ=0.5,d=2 when varying k∈{2,3,4}. This is an example of the VEEE model since λ=2 and ρ=0.5.Fig. 4 Examples of scatterplots, with n=500 observations, from a data generating MLN mixture model with θkδ=5,λ=2,ρ=0.5,d=2, with k=2 (left), k=3 (middle), and k=4 (right)

Model selection

In this subsection, we investigate how well the BIC picks the correct number of components while varying the distance between clusters (δ), the dimension (d), and the sample size (n). Three simulation studies were performed to investigate properties of using the BIC: (1) to select the model while fixing k, (2) to select k while fixing the model, and (3) to jointly select k and the model. The parameters explored were generally the same for each of the three studies but some care was needed when setting δ.

First, we examined the ability of the BIC to pick the model while fixing k. We considered 72 scenarios arising from the combination of the factors n∈{50,100}, ρ∈{0,0.5}, λ∈{1,2}, d∈{2,8,16}, k∈{2,3,4}, and δ=3. For each scenario, 100 data sets were generated and all the models were fitted while using the correct k; then, the BIC was used to select the model. Table 6 gives the number of times each model (on the row) was selected by varying the model type, described by the pair λ,ρ, and the sample size. The values we provide are aggregated over k and d.Table 6 Number of times each parsimonious model was selected by the BIC

	Data generated from model with λ,ρ	
EIIE (1, 0)	VIIE (2, 0)	EEEE (1, 0.5)	VEEE (2, 0.5)	
n	50	100	50	100	50	100	50	100	
EIIE	816	859	86	37	137	10	28	5	
VIIE	41	20	720	830	72	1	192	6	
EEEE	3	2	6	2	371	569	93	47	
VEEE	1		4	1	27	12	254	517	
VVEE	2		2		36	61	140	231	
EVEE	2		3		125	200	26	6	
VVEV					8	6	42	62	
VIIV			37	10	6		34		
EVEV	1				44	24	4	1	
EEEV					44	3	8	2	
VEEV					2		40	13	
EEVE	3		2	2	16	14	11		
EEIE	18	16	8	3			2		
VEVE			4	2	5		19	9	
VEIE	1		11	6	1		2		
EVIE	6	1	6	2			1		
EIIV	6	1	4		5				
VVIE			4	2					
EVVE			1		1		2		
VVVE			1				1	1	
EVVV		1	1						
VEIV				1					
VVIV				1					
EEVV							1		
VVVV				1					
EEIV									
EVIV									
VEVV									
Total	900	900	900	900	900	900	900	900	
Each column represents a summary of 900 data sets generated from a MLN mixture model with θkδ,λ,ρ,d using k∈{2,3,4}, d∈{2,8,16}, and δ=3

Second, we examined the ability of the BIC to pick k while fixing the model type. We considered 72 scenarios arising from the combination of the factors n∈{100,200,300}, ρ∈{0,0.5}, λ∈{1,2}, d∈{2,16}, and k∈{2,3,4}. The parameter δ was chosen so that when n varied across its levels, {100,200,300}, the frequency of the number of components selected, say k~, would be not a degenerate distribution (see Table 7), i.e., if k=2 and k~=2 across n would be not illustrate the effect of increasing the sample size. For each scenario, 100 data sets were generated. For each of them, we fitted the MLN mixture model for k∈{1,…,8} while using the correct model; then, the BIC was used to select k. Table 7 gives the number of times each value of k was selected by the BIC, for each scenario, when d=2. For each model, as we increase n, the BIC tends to increase the frequency of selecting the correct model. Moving from EIIE to VIIE, or from EEE to VEEE, tends to increases the frequency of k~=k. When k=4, there is more variety in the selected number of components. For the results related to the case d=16, see Table E.6 in the supplementary material.Table 7 Number of times each value of k is selected by the BIC when fixing the model structure and varying the number of components

Data generated from θkδ,λ,ρ,d=2	Selected k (k~)	
k	Model	ρ	λ	δ	n	1	2	3	4	5	6	7	8	
2	EIIE	0.0	1	0.50	100	70	30							
EIIE	0.0	1	0.50	200	23	77							
EIIE	0.0	1	0.50	300	6	94							
EEEE	0.5	1	0.75	100	81	19							
EEEE	0.5	1	0.75	200	54	46							
EEEE	0.5	1	0.75	300	34	66							
VIIE	0.0	2	0.50	100	13	86	1						
VIIE	0.0	2	0.50	200	1	97	2						
VIIE	0.0	2	0.50	300		99					1		
VEEE	0.5	2	0.75	100	17	82	1						
VEEE	0.5	2	0.75	200		98	1					1	
VEEE	0.5	2	0.75	300		100							
3	EIIE	0.0	1	1.50	100		35	64	1					
EIIE	0.0	1	1.50	200		6	91	3					
EIIE	0.0	1	1.50	300		2	96	2					
EEEE	0.5	1	2.00	100	23	41	36						
EEEE	0.5	1	2.00	200		15	85						
EEEE	0.5	1	2.00	300		2	98						
VIIE	0.0	2	1.50	100		50	43	6			1		
VIIE	0.0	2	1.50	200		12	85	3					
VIIE	0.0	2	1.50	300		5	93	2					
VEEE	0.5	2	2.00	100	5	73	21	1					
VEEE	0.5	2	2.00	200		37	63						
VEEE	0.5	2	2.00	300		16	82	1	1				
4	EIIE	0.0	1	1.25	100	69		2	24	5				
EIIE	0.0	1	1.25	200	3			94	3				
EIIE	0.0	1	1.25	300				99	1				
EEEE	0.5	1	1.25	100	20	33	1	40	6				
EEEE	0.5	1	1.25	200		1		95	4				
EEEE	0.5	1	1.25	300				100					
VIIE	0.0	2	1.25	100	39	25	17	18	1				
VIIE	0.0	2	1.25	200	7	12	16	62	3				
VIIE	0.0	2	1.25	300			2	95	3				
VEEE	0.5	2	1.25	100	18	20	27	33	1		1		
VEEE	0.5	2	1.25	200		1	7	90	1	1			
VEEE	0.5	2	1.25	300			2	96	2				
Each row refers to 100 data sets generated from a MLN mixture with θkδ,λ,ρ,d=2

Third, we examined the ability of the BIC to pick both the model and k. As for the previous case, we considered 72 scenarios arising from the combination of the factors n∈{100,200,300}, ρ∈{0,0.5}, λ∈{1,2}, d∈{2,16}, and k∈{2,3,4}. Again, the parameter δ is chosen so that when n varied across the values {100,200,300}, the frequency of k~ would not be a degenerate distribution (see Table 8), i.e., δ was chosen so that k~≠k for each replication with n∈{100,200,300}. For each scenario, 100 data sets were generated. For each of them, we found the pair of model and k selected by the BIC when the search space was k∈{1,…,8} and all the possible model structures. Table 8 gives the number of times each pair model,k was selected by the BIC, for each scenario, when d=2. The last two columns give the frequency of picking the number of components and the correct model (second last column) and picking the number of components and one of models in the subset shown in third column group of Table 8. This means that, by definition, the values on the second last column are lower than, or equal to, the values on the last column. For each combination, there is an improvement in the model and component selection when increasing n both marginally and jointly. The case k=3 has the mostly variety in the chosen model and k=4 has the mostly variety in k~. For the results when d=16, see Table E.7 in the supplementary material.Table 8 Counts of the BIC choices, in the case d=2, when the search is over k=1,…,8 and all or some of the parsimonious structures

Data generated from θkδ,λ,ρ,d=2	Selected k (k~)	Selected fitted model	When k=k~ and	
k	Model	λ	ρ	δ	n	1	2	3	4	≥5	EIIE	VIIE	VEEE	EEEE	VVEE	Model	Subset	
2	EIIE	1	0.0	0.5	100	70	29			1	73	5	4	10		22	27	
2	EIIE	1	0.0	0.5	200	48	52				76	2	6	9		50	52	
2	EIIE	1	0.0	0.5	300	28	71		1		66	6	6	9	3	64	71	
2	EEEE	1	0.5	0.8	100	75	22	3			12		22	21	6	8	20	
2	EEEE	1	0.5	0.8	200	41	57			2			55	13	4	48	51	
2	EEEE	1	0.5	0.8	300	23	77						79	6	1	75	75	
2	VIIE	2	0.0	0.5	100	33	60	6		1	43	37	2	3		35	48	
2	VIIE	2	0.0	0.5	200	21	67	9	3		21	59	4	2		55	57	
2	VIIE	2	0.0	0.5	300	13	68	16	2	1	10	65	1	2	2	58	58	
2	VEEE	2	0.5	0.8	100	17	77	5	1		4	9	9	6	54	54	70	
2	VEEE	2	0.5	0.8	200		98			2				2	91	90	92	
2	VEEE	2	0.5	0.8	300		98			2				2	92	92	94	
3	EIIE	1	0.0	1.5	100	18	50	30	2		44	7	13	8	2	24	24	
3	EIIE	1	0.0	1.5	200	3	24	70	2	1	69	1	16	3		65	65	
3	EIIE	1	0.0	1.5	300		12	86	1	1	84	2	11			83	85	
3	EEEE	1	0.5	2.0	100	23	40	35	1	1	10		51	7	10	23	33	
3	EEEE	1	0.5	2.0	200	1	22	74	1	2	1		90		6	73	74	
3	EEEE	1	0.5	2.0	300		6	92	1	1			92	1	2	87	89	
3	VIIE	2	0.0	1.5	100	5	70	23	2		8	56	1	5	13	10	17	
3	VIIE	2	0.0	1.5	200		36	62	1	1	4	65		3	15	51	57	
3	VIIE	2	0.0	1.5	300		8	91		1		89			4	87	87	
3	VEEE	2	0.5	2.0	100	8	70	17	5		7	2	7	18	51	7	14	
3	VEEE	2	0.5	2.0	200		52	45	1	2		1	4	21	68	38	44	
3	VEEE	2	0.5	2.0	300		24	72	4					11	86	71	71	
4	EIIE	1	0.0	1.2	100	37	39	1	22	1	60	8	5			21	22	
4	EIIE	1	0.0	1.2	200	4	6	2	88		91					87	87	
4	EIIE	1	0.0	1.2	300			1	99		96					96	96	
4	EEEE	1	0.5	1.2	100	16	26	5	51	2	31	7	47		3	34	50	
4	EEEE	1	0.5	1.2	200	2	3	2	88	5	10	2	84			79	86	
4	EEEE	1	0.5	1.2	300				98	2			96		1	95	96	
4	VIIE	2	0.0	1.2	100	33	32	23	9	3	42	31		3	1	3	7	
4	VIIE	2	0.0	1.2	200	3	9	15	69	4	32	55				39	66	
4	VIIE	2	0.0	1.2	300		1	4	93	2	16	82				77	92	
4	VEEE	2	0.5	1.2	100	20	37	17	23	3	29	22	20		7	3	21	
4	VEEE	2	0.5	1.2	200		5	12	76	7	9	8	29		47	43	75	
4	VEEE	2	0.5	1.2	300			4	93	3	2		24		70	67	92	
There are four blocks of columns. In the first one, there are the parameters/quantities used to generated the 100 data sets. From the second block onwards, there is the number of times: each value of k is picked (2nd block), each model is selected (3rd block), and each pair model,k is selected (4nd block), where “model” belongs to the whole family, in the 2nd last column, and to the subset of models of the 3rd block in the last column

Data generated from other distributions

In this subsection, we examine the effect on the fitted MLN mixture when the data is generated from another mixture model which exhibits component-wise heavy tails and/or skewness. We consider generating data from finite mixtures of multivariate generalized hyperbolic (GHD) (Browne and McNicholas 2015; Barndorff-Nielsen and Halgreen 1977) and t (Peel and McLachlan 2000; Andrews and McNicholas 2012) distributions. We consider a parameter setup similar to the one given in (16) for the two-component MLN mixture, but we replace β with the corresponding parameter(s). In particular, for the t mixture we have19 θtν=π1=π2=12,(μ1,μ2)=(-31d,+31d),Σj=Id,ν1=ν2=ν,

where the degrees of freedom ν affects the tails, while for the GHD mixture we have20 θGHDλ,α=πj=12,(μ1,μ2)=(-31d,+31d),Σj=Id,ωj=1,λj=λ,(α1,α2)=α(-1d,0d),

where λ affects the tails and α affects the skewness for one of the components.

Figure 5 shows two examples of the MLN mixture model fitted to generated data. The data in the left panel are generated from a two-component t mixture with θt(4). The BIC selected the VIIE MLN mixture model with k=3 componnets. Two components match the location of the two multivariate-t to the centers, and then the third component is used to model the observations that deviate from normality. The data in the right panel are generated from a GHD mixture with θGHD(1,-1); they display skewness in the bottom-left component. In this case, the BIC still selected the VIIE MLN mixture model, but now with k=4 components. One component is used to model the group on the top-right with no skewness, but three components are needed to model the skewness of the other group. In both examples, the ARI between fitted and generated labels was 0.921 (t mixture) and 0.729 (GHD mixture).Fig. 5 Examples of MLN mixtures fitted to data generated from a 2-component t mixture with ν=4 (left panel), and from a 2-component GHD mixture with (λ,ω,α)=(1,1,-1) (right panel). The BIC selected the VIIE model with k=3, in the left panel, and with k=4 in the right panel

To investigate the effect of a misspecified model, we generate samples of size n=500 from either the t or GHD mixture and we fit the MLN mixture with k=1,…,8 and all the 28 possible parsimonious structures. The BIC is used to pick k and the model; then, we calculate the ARI between the fitted and generated labels. For the t-mixture we use ν={2,10} and d={2,8} while for the GHD mixture we use (λ,α)=(1,0),(1,-1),(-1,0),(-1,-1). For both mixtures, we consider d={2,8} and, for each combination of parameters, we generate 100 data sets.

Table 9 shows the frequency of the number of components and the model picked by the BIC when fitting the MLN mixture to data generated from two-component t or GHD mixtures. We display all possible k=1,…,8 whereas we only display the four selected models with the highest frequency. The last two columns shows the mean and standard deviations for the ARI. This indicates, on average, that the MLN mixture does will at the group structure even thought the BIC does not pick the correct number of components; e.g., when data are generated from a GHD-mixture with α=-1, the BIC picks k=2 three times from the 400 replications.Table 9 The frequency of the number of components (k) and model picked by the BIC when fitting the MLN mixture to data generated from two-component t or GHD mixtures

θtν	d	Selected k (k~)	Selected model	ARI	
1	2	3	4	5	6	7	8	EIIE	VIIE	VIIV	VVEE	Avg	SD	
4	2		41	48	11					88	11	1		0.970	0.019	
4	8		4	29	17	19	15	9	7	77	19	4		0.967	0.047	
8	2		83	15	2					98	2			0.993	0.007	
8	8		74	23	3					92	7			0.999	0.002	
θGHD(λ,α)																
(1,0)	2		88	9	3					94	4			0.957	0.017	
(1,0)	8		71	25	4					75	24	1		0.985	0.050	
(-1,0)	2		74	19	6	1				94	6			0.994	0.015	
(-1,0)	8		33	31	16	20				88	10	1		0.994	0.010	
(1,-1)	2		2	46	34	10	5	3		13	36	15	22	0.753	0.060	
(1,-1)	8			27	48	20	4	1			16	40	16	0.691	0.078	
(-1,-1)	2		1	51	29	16	2	1		33	47	9	2	0.863	0.071	
(-1,-1)	8			61	19	6	4	4	6	22	33	40		0.825	0.059	

Conclusions

In this paper, we have improved the mixture of multivariate leptokurtic-normal distributions from two points of view. Firstly, we have defined parsimonious variants of the model by putting convenient constraints on the component covariance matrices and excess kurtoses. The use of parsimonious models can lead to more precise parameter estimates and data classifications (Tomarchio et al. 2022). These better performances occur despite the reduction in the number of estimated parameters, allowing for a more parsimonious data modelization. Furthermore, as well documented in the literature (Flury 1988; Greselin et al. 2011; Greselin and Punzo 2013), it should be remarked that parsimonious approaches can disclose richer information on the underlying data structure than the classical unconstrained method. Secondly, we have made maximum likelihood estimation more computationally efficient by introducing two algorithms based on the majorization-minimization algorithm and the fixed point approximation. We highlighted the advantages of our methodological and computational proposals using simulated and real data.

Supplementary Information

Below is the link to the electronic supplementary material.Supplementary file 1 (pdf 1497 KB)

Acknowledgements

Antonio Punzo has been partially supported by MUR, Grant Number 2022XRHT8R)—“The SMILE project: Statistical Modelling and Inference to Live the Environment”. This paper was also developed within the project funded by Next Generation EU)—“GRINS–Growing Resilient, INclusive and Sustainable” project (PE0000018), National Recovery and Resilience Plan (NRRP)—PE9–Mission 4, C2, Intervention 1.3. Ryan Browne has been partially supported by the Natural Sciences and Engineering Research Council of Canada through their Discovery Grants program (RGPIN-2018-04444).

Funding

Open access funding provided by Università Cattolica del Sacro Cuore within the CRUI-CARE Agreement.

Declarations

Conflict of interest

Ryan P. Browne declares that he has no conflict of interest. Luca Bagnato declares that he has no conflict of interest. Antonio Punzo declares that he has no conflict of interest.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

Absil P-A Mahony R Sepulchre R Optimization algorithms on matrix manifolds 2008 Princeton, NJ Princeton University Press
Absil P-A, Mahony R, Sepulchre R (2008) Optimization algorithms on matrix manifolds. Princeton University Press, Princeton, NJ
Aitken AC On Bernoulli’s numerical solution of algebraic equations Proc R Soc Edinburgh 1926 46 289 305
Aitken AC (1926) On Bernoulli’s numerical solution of algebraic equations. Proc R Soc Edinburgh 46:289–305
Andrews JL McNicholas PD Model-based clustering, classification, and discriminant analysis via mixtures of multivariate t-distributions Stat Comput 2012 22 5 1021 1029
Andrews JL, McNicholas PD (2012) Model-based clustering, classification, and discriminant analysis via mixtures of multivariate t-distributions. Stat Comput 22(5):1021–1029
Azzalini A (2015) R package sn: The skew-normal and skew- distributions. Version 1.3-0
Azzalini A Browne RP Genton MG McNicholas PD On nomenclature for, and the relative merits of, two formulations of skew distributions Stat Probab Lett 2016 110 201 206
Azzalini A, Browne RP, Genton MG, McNicholas PD (2016) On nomenclature for, and the relative merits of, two formulations of skew distributions. Stat Probab Lett 110:201–206
Bagnato L Punzo A Unconstrained representation of orthogonal matrices with application to common principal components Comput Stat 2021 36 2 1177 1195
Bagnato L, Punzo A (2021) Unconstrained representation of orthogonal matrices with application to common principal components. Comput Stat 36(2):1177–1195
Bagnato L Punzo A Zoia MG The multivariate leptokurtic-normal distribution and its application in model-based clustering Can J Stat 2017 45 1 95 119
Bagnato L, Punzo A, Zoia MG (2017) The multivariate leptokurtic-normal distribution and its application in model-based clustering. Can J Stat 45(1):95–119
Bagnato L Punzo A Zoia MG Leptokurtic moment-parameterized elliptically contoured distributions with application to financial stock returns Commun Stat Theory Methods 2022 51 2 486 500
Bagnato L, Punzo A, Zoia MG (2022) Leptokurtic moment-parameterized elliptically contoured distributions with application to financial stock returns. Commun Stat Theory Methods 51(2):486–500
Banfield JD Raftery AE Model-based Gaussian and non-Gaussian clustering Biometrics 1993 49 3 803 821
Banfield JD, Raftery AE (1993) Model-based Gaussian and non-Gaussian clustering. Biometrics 49(3):803–821
Barndorff-Nielsen O Halgreen C Infinite divisibility of the hyperbolic and generalized inverse Gaussian distributions Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 1977 38 309 311
Barndorff-Nielsen O, Halgreen C (1977) Infinite divisibility of the hyperbolic and generalized inverse Gaussian distributions. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 38:309–311
Böhning D Dietz E Schaub R Schlattmann P Lindsay B The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family Ann Inst Stat Math 1994 46 373 388
Böhning D, Dietz E, Schaub R, Schlattmann P, Lindsay B (1994) The distribution of the likelihood ratio for mixtures of densities from the one-parameter exponential family. Ann Inst Stat Math 46:373–388
Bouveyron C, Celeux G, Murphy TB, Raftery AE (2019) Model-based clustering and classification for data science: with applications in R. Cambridge series in statistical and probabilistic mathematics. Cambridge, United Kingdom: Cambridge University Press
Boyd SP Vandenberghe L Convex optimization 2006 corrections Cambridge Cambridge University Press
Boyd SP, Vandenberghe L (2006) Convex optimization, corrections. Cambridge University Press, Cambridge
Browne R McNicholas P Estimating common principal components in high dimensions Adv Data Anal Classif 2014 8 2 217 226
Browne R, McNicholas P (2014) Estimating common principal components in high dimensions. Adv Data Anal Classif 8(2):217–226
Browne R McNicholas P Orthogonal Stiefel manifold optimization for eigen-decomposed covariance parameter estimation in mixture models Stat Comput 2014 24 2 203 210
Browne R, McNicholas P (2014) Orthogonal Stiefel manifold optimization for eigen-decomposed covariance parameter estimation in mixture models. Stat Comput 24(2):203–210
Browne RP (2022) Revitalizing the multivariate elliptical leptokurtic-normal distribution and its application in model-based clustering. Stat Probab Lett
Browne RP McNicholas PD A mixture of generalized hyperbolic distributions Can J Stat 2015 43 2 176 198
Browne RP, McNicholas PD (2015) A mixture of generalized hyperbolic distributions. Can J Stat 43(2):176–198
Browne RP, Punzo A, Bagnato L (2023) LeptokurticMixture: implements parsimonious finite mixtures of multivariate elliptical leptokurtic-normals. R package version 1.0
Celeux G Govaert G Gaussian parsimonious clustering models Pattern Recogn 1995 28 5 781 793
Celeux G, Govaert G (1995) Gaussian parsimonious clustering models. Pattern Recogn 28(5):781–793
Cook RD Weisberg S An introduction to regression graphics 1994 New York Wiley
Cook RD, Weisberg S (1994) An introduction to regression graphics. Wiley, New York
Dang UJ Browne RP McNicholas PD Mixtures of multivariate power exponential distributions Biometrics 2015 71 4 1081 1089 26134429
Dang UJ, Browne RP, McNicholas PD (2015) Mixtures of multivariate power exponential distributions. Biometrics 71(4):1081–108926134429
Dang UJ Punzo A McNicholas PD Ingrassia S Browne RP Multivariate response and parsimony for Gaussian cluster-weighted models J Classif 2017 34 1 4 34
Dang UJ, Punzo A, McNicholas PD, Ingrassia S, Browne RP (2017) Multivariate response and parsimony for Gaussian cluster-weighted models. J Classif 34(1):4–34
Dempster AP Laird NM Rubin DB Maximum likelihood from incomplete data via the EM algorithm J Roy Stat Soc B 1977 39 1 1 38
Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood from incomplete data via the EM algorithm. J Roy Stat Soc B 39(1):1–38
Farcomeni A Punzo A Robust model-based clustering with mild and gross outliers TEST 2020 29 4 989 1007
Farcomeni A, Punzo A (2020) Robust model-based clustering with mild and gross outliers. TEST 29(4):989–1007
Flury BN Common principal components and related multivariate models 1988 New York Wiley
Flury BN (1988) Common principal components and related multivariate models. Wiley, New York
Fraley C Raftery AE Model-based clustering, discriminant analysis, and density estimation J Am Stat Assoc 2002 97 611 631
Fraley C, Raftery AE (2002) Model-based clustering, discriminant analysis, and density estimation. J Am Stat Assoc 97:611–631
Fraley C Raftery AE Scrucca L Murphy TB Fop M mclust: normal mixture modeling for model-based clustering, classification, and density estimation R Package Version 2015 5 1
Fraley C, Raftery AE, Scrucca L, Murphy TB, Fop M (2015) mclust: normal mixture modeling for model-based clustering, classification, and density estimation. R Package Version 5:1
Galimberti G Soffritti G A multivariate linear regression analysis using finite mixtures of t distributions Comput Stat Data Anal 2014 71 138 150
Galimberti G, Soffritti G (2014) A multivariate linear regression analysis using finite mixtures of distributions. Comput Stat Data Anal 71:138–150
Gallaugher MP Tomarchio SD McNicholas PD Punzo A Model-based clustering via skewed matrix-variate cluster-weighted models J Stat Comput Simul 2022 92 13 2645 2666
Gallaugher MP, Tomarchio SD, McNicholas PD, Punzo A (2022) Model-based clustering via skewed matrix-variate cluster-weighted models. J Stat Comput Simul 92(13):2645–2666
Greselin F Ingrassia S Punzo A Assessing the pattern of covariance matrices via an augmentation multiple testing procedure Stat Methods Appl 2011 20 2 141 170
Greselin F, Ingrassia S, Punzo A (2011) Assessing the pattern of covariance matrices via an augmentation multiple testing procedure. Stat Methods Appl 20(2):141–170
Greselin F Punzo A Closed likelihood ratio testing procedures to assess similarity of covariance matrices Am Stat 2013 67 3 117 128
Greselin F, Punzo A (2013) Closed likelihood ratio testing procedures to assess similarity of covariance matrices. Am Stat 67(3):117–128
Hubert L Arabie P Comparing partitions J Classif 1985 2 1 193 218
Hubert L, Arabie P (1985) Comparing partitions. J Classif 2(1):193–218
Hunter DL Lange K A tutorial on MM algorithms Am Stat 2004 58 1 30 37
Hunter DL, Lange K (2004) A tutorial on MM algorithms. Am Stat 58(1):30–37
Kiers H Setting up alternating least squares and iterative majorization algorithms for solving various matrix optimization problems Comput Stat Data Anal 2002 41 157 170
Kiers H (2002) Setting up alternating least squares and iterative majorization algorithms for solving various matrix optimization problems. Comput Stat Data Anal 41:157–170
Lange K Hunter DR Yang I Optimization transfer using surrogate objective functions J Comput Graph Stat 2000 9 1 1 20
Lange K, Hunter DR, Yang I (2000) Optimization transfer using surrogate objective functions. J Comput Graph Stat 9(1):1–20
Laub A A Schur method for solving algebraic Riccati equations IEEE Trans Autom Control 1979 24 6 913 921
Laub A (1979) A Schur method for solving algebraic Riccati equations. IEEE Trans Autom Control 24(6):913–921
Lindsay BG (1995) Mixture models: theory, geometry and applications. In NSF-CBMS regional conference series in probability and statistics, Volume 5, pp 63–65. California: Institute of Mathematical Statistics: Hayward
McNicholas PD Mixture model-based classification 2016 Boca Raton CRC Press
McNicholas PD (2016) Mixture model-based classification. CRC Press, Boca Raton
McNicholas PD Model-based clustering J Classif 2016 33 3 331 373
McNicholas PD (2016) Model-based clustering. J Classif 33(3):331–373
Morris K Punzo A McNicholas PD Browne RP Asymmetric clusters and outliers: Mixtures of multivariate contaminated shifted asymmetric Laplace distributions Comput Stat Data Anal 2019 132 145 166
Morris K, Punzo A, McNicholas PD, Browne RP (2019) Asymmetric clusters and outliers: Mixtures of multivariate contaminated shifted asymmetric Laplace distributions. Comput Stat Data Anal 132:145–166
Peel D McLachlan GJ Robust mixture modelling using the t distribution Stat Comput 2000 10 4 339 348
Peel D, McLachlan GJ (2000) Robust mixture modelling using the distribution. Stat Comput 10(4):339–348
Pocuca N Browne RP McNicholas PD Mixture: mixture models for clustering and classification R Package Version 2022 2 5
Pocuca N, Browne RP, McNicholas PD (2022) Mixture: mixture models for clustering and classification. R Package Version 2:5
Punzo A Bagnato L Allometric analysis using the multivariate shifted exponential normal distribution Biom J 2020 62 6 1525 1543 32240556
Punzo A, Bagnato L (2020) Allometric analysis using the multivariate shifted exponential normal distribution. Biom J 62(6):1525–154332240556
Punzo A Bagnato L The multivariate tail-inflated normal distribution and its application in finance J Stat Comput Simul 2021 91 1 1 36
Punzo A, Bagnato L (2021) The multivariate tail-inflated normal distribution and its application in finance. J Stat Comput Simul 91(1):1–36
Punzo A Bagnato L Multiple scaled symmetric distributions in allometric studies Int J Biostat 2022 18 1 219 242
Punzo A, Bagnato L (2022) Multiple scaled symmetric distributions in allometric studies. Int J Biostat 18(1):219–242
Punzo A Blostein M McNicholas PD High-dimensional unsupervised classification via parsimonious contaminated mixtures Pattern Recogn 2020 98 107031
Punzo A, Blostein M, McNicholas PD (2020) High-dimensional unsupervised classification via parsimonious contaminated mixtures. Pattern Recogn 98:107031
Punzo A Mazza A McNicholas PD ContaminatedMixt: an R package for fitting parsimonious mixtures of multivariate contaminated normal distributions J Stat Softw 2018 85 10 1 25 30505247
Punzo A, Mazza A, McNicholas PD (2018) ContaminatedMixt: an R package for fitting parsimonious mixtures of multivariate contaminated normal distributions. J Stat Softw 85(10):1–2530505247
Punzo A McNicholas PD Parsimonious mixtures of multivariate contaminated normal distributions Biom J 2016 58 6 1506 1537 27510372
Punzo A, McNicholas PD (2016) Parsimonious mixtures of multivariate contaminated normal distributions. Biom J 58(6):1506–153727510372
Schwarz G Estimating the dimenson of a model Ann Stat 1978 6 461 464
Schwarz G (1978) Estimating the dimenson of a model. Ann Stat 6:461–464
Scrucca L Fop M Murphy TB Raftery AE mclust 5: clustering, classification and density estimation using Gaussian finite mixture models R J 2016 8 1 289 317 27818791
Scrucca L, Fop M, Murphy TB, Raftery AE (2016) mclust 5: clustering, classification and density estimation using Gaussian finite mixture models. R J 8(1):289–31727818791
Stephens M Dealing with label switching in mixture models J R Stat Soc Ser B 2000 62 795 809
Stephens M (2000) Dealing with label switching in mixture models. J R Stat Soc Ser B 62:795–809
Titterington DM Smith AFM Makov UE Statistical analysis of finite mixture distributions 1985 New York Wiley
Titterington DM, Smith AFM, Makov UE (1985) Statistical analysis of finite mixture distributions. Wiley, New York
Tomarchio SD Bagnato L Punzo A Model-based clustering via new parsimonious mixtures of heavy-tailed distributions AStA Adv Stat Anal 2022 106 2 315 347
Tomarchio SD, Bagnato L, Punzo A (2022) Model-based clustering via new parsimonious mixtures of heavy-tailed distributions. AStA Adv Stat Anal 106(2):315–347
Tomarchio SD Punzo A Maruotti A Parsimonious hidden Markov models for matrix-variate longitudinal data Stat Comput 2022 32 3 53 35730052
Tomarchio SD, Punzo A, Maruotti A (2022) Parsimonious hidden Markov models for matrix-variate longitudinal data. Stat Comput 32(3):5335730052
Weisberg S (2011) alr3: data to accompany applied linear regression 3rd edn. http://CRAN.R-project.org/package=alr3
Wonham WM On a matrix Riccati equation of stochastic control SIAM J Control 1968 6 4 681 697
Wonham WM (1968) On a matrix Riccati equation of stochastic control. SIAM J Control 6(4):681–697
Zhang J Liang F Robust clustering using exponential power mixtures Biometrics 2010 66 4 1078 1086 20163406
Zhang J, Liang F (2010) Robust clustering using exponential power mixtures. Biometrics 66(4):1078–108620163406
