# Estimation of Parameters of Finite Mixture of Rayleigh Distribution by the Expectation-Maximization Algorithm

Authors: Noor Mohammed; Fadhaa Ali
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In the lifetime process in some systems, most data cannot belong to one single population. In fact, it can represent several subpopulations. In such a case, the known distribution cannot be used to model data. Instead, a mixture of distribution is used to modulate the data and classify them into several subgroups. The mixture of Rayleigh distribution is best to be used with the lifetime process. This paper aims to infer model parameters by the expectation-maximization (EM) algorithm through the maximum likelihood function. The technique is applied to simulated data by following several scenarios. The accuracy of estimation has been examined by the average mean square error (AMSE) and the average classification success rate (ACSR). The results showed that the method performed well in all simulation scenarios with respect to different sample sizes.


---

## Body

## 1. Introduction

Fitting an appropriate distribution of data is one of the recent issues that challenge research worldwide. One of which arises from the date when it belongs to more than one population. Therefore, finite mixture models are employed to tackle such issues. In many applications, one large population can be formed from several subpopulations that are mixed in an unknown proportion. In such a case, one can come across a population of biological or electrical elements that may be divided into two or more subpopulations based on the possible causes of failure or any other categories such as male, female, and age categories. Many researchers have introduced several inferential methods to infer some aspects of data following heavy-tailed distributions [1, 2].Many methods have been proposed to apply the mixture model to data in physical, biological, engineering sciences, and other fields as a valid technique of machine learning methods. A mixture of normal and Laplace distributions is fitted to wind shear data [3, 4]. In addition, the mixture model was used in modeling the crime and justice data [5]. A mixture of exponential distributions was considered a valid tool to infer device failure by dividing the population into several subpopulations [6, 7]. Split-and-merge operations were used to improve the likelihood via the mixture of the t-distribution model to study image compression and pattern [8]. The mixture model was used as a classifier tool to cluster individuals based on estimated density functions for each individual based on blood cell data [9].The mixture model can be structured from densities following the same families of distributions or different ones. Either way, it is useful as a confirmatory or predicting tool. Indeed, the maximum likelihood (ML) is the popular method for inferring the model [10, 11]. However, the challenge is the unknown number of allocations and observation memberships, through which the accuracy of estimation is undermined. These issues have been tackled by many proposed techniques over years. Therefore, determining the exact number of components by performing a relevant test is considered an important starting point in applying the mixture model [12, 13]. In the same manner, the most common way of handling the above issues is by using a latent variable that leads to a complete-data log likelihood rather than using the incomplete one [14], and then, the expectation-maximization (EM) algorithm is employed to estimate model parameters [15–17].Moreover, Markov hidden model can be employed by adding components-indicator labels to the mixtue model to determine observations memberships. This is an extension of Gaussian mixture model to model k-components of stationary or non-stationary auto-regressive [18], or to model time series data by non-Gaussian mixture [19]. Bayesian inference can also be used to infer the parameters of mixture models by assuming prior distribution for the parameters which in turn result in the posterior distribution of the parameters. This framework can be conducted by Markov chain Monte Carlo (MCMC) [20, 21]. This framework of estimation can avoid many difficulties in its applications that involve (ML) estimation [22–24]. The authors in [25] introduced a general form of inverse Rayleigh distribution, which is known as exponentiated inverse Rayleigh distribution (EIRD), which extends a more flexible distribution for modeling life data.In this paper, we derived the expectation-maximization (EM) algorithm to estimate parameters of the mixture of Rayleigh distributions. This can be done through the maximum likelihood of the model, and then, we can use the EM algorithm to update estimators by iterative steps until reaching convergence. Afterward, the accuracy of the model inferred is examined by simulation based on different scenarios and the number of components. The success rate of clustering is also calculated accordingly. The rest of the paper is organized to give the theoretical details of the theoretical part in Section2. In Section 3, the method is applied to the simulated data by using several scenarios. In Sections 4 and 5, the discussion and conclusion are involved to show the complete picture of the method and its results.


## 2. Materials and Methods

A random variable of Rayleigh distribution is defined as the time between two consecutive occurrences [26]. The probability density function is written as(1)fx;σ=xσ2e−x2/2σ2,wherex,σ>0.In many lifetime scenarios, data can come from more than one component, and then, known distributions become inaccurate to fit the data. In such a case, the mixture model can be used to handle the problem. When it comes to lifetime data, a finite mixture of Rayleigh distribution is used to model the data. The probability density function (PDF) of the mixture of Rayleigh distribution withk components can be written as(2)fx;θ=∑j=1kπjxσj2e−x2/2σj2,where ∑j=1kπj=1,0<πj<1,θ=π1,π2,…,πk−1,σ1,σ2,…,σkT, x,σ>0.


## 3. Inferential Method

The classical way of estimating the parameters of the mixture model is by the maximum likelihood (ML). The incomplete-data likelihood function of thek-component mixture model is defined as follows:(3)Lθ;x=∏i=1n∑j=1kπjxiσj2e−xi2/2σj2,where X=x1,x2,…,xnT.The allocations here are considered a hidden state that needs to be defined by a latent vectorZ. The variable Z is a Markovian vector, and it is defined as Z=z1,z2,…,znT, where Zi=zi1,zi2,…,zikT. zij is defined as(4)zij=1,fitemibelong to componentj,0,otherwise.Hence, the complete-data likelihood function is written as(5)Lθ;X,Z=∏i=1n∏j=1kπjxiσj2e−xi2/2σj2zij,where X=x1,x2,…,xnT.Taking the logarithm of Equation (2) results in(6)ℓθ;X,Z=∑i=1n∑j=1kzijlogπj+logxi−logσj2+−xi22σj2.Many methods have been proposed to estimate parameters of Equation (3). One of which is the expectation-maximization (EM) algorithm that is used with some mixture models [27]. In this paper, the EM algorithm is derived for the k-component mixture of Rayleigh distribution. The next section consists of the details of such a method.


## 4. EM Algorithm

Equation (6) is the complete-data log-likelihood function that needs to be maximized given the data. The EM algorithm is used here to estimate the model parameters by considering zij missing data. It consists of two iterative steps, E (expectation step) and M (maximization step). For m=1,2,…,M, we repeat the following two steps.


### 4.1.E-Step

The unobservable data (herezij) are handled by the E-step, which takes the expectation of the complete-data log likelihood conditionally on the given observed data X, using the value of θ of the current iteration. Let θ0 be the value chosen initially for θ. The first iteration of the EM algorithm computes the conditional expectation of log Lθ given X by using θ0 for θ, which can be written as(7)Qθ;θ0=Eθ0logLθ|x.As logθ is linear in the missing data zij, the E-step at the m+1th iteration requires the current expectation of the random variable Zij, which corresponds to zij, to be calculated conditionally on the given observed data X.(8)Eθ0Zij|x=PrλmZij=1|x=τijmxi;θm.Then, we get(9)τijmxi;θm=πjmfjxi;σjm∑j=1kπjmfjxi;σjm.The quantityτjmx;θm is the posterior probability, which implies that the ith element of the sample with the observed value xi belongs to the jth component of the mixture. By substituting zij into Equation (3) with τjmx;θm, the Q function becomes(10)Qθ;θm=∑i=1n∑j=1kτijmxi;θmlogπj+logxi−logσj2+−xi22σj2.


### 4.2.M-Step

To get the updated estimateθm+1, the global maximization of Qθ;θm with respect to θ over the parameter space is required by the M-step on the m+1th iteration. Given the updated τijmxi;θm at the iteration mth, the weighted parameter πjm+1 can be estimated by taking the first derivative of Qθ;θm with respect to πj and equating it to 0, which results in(11)π^jm+1=∑i=1nτijmn.The component parameterσj at the iteration m+1th can be computed by deriving Qθ;θm with respect to σj and equating it to 0, and we get(12)∂Q∂σj=∑i=1n−2σjτijσj2+τijxi2σj3⟶−2∑i=1nτijσj+∑i=1nτijxi2σj3=0⟶2∑i=1nτijσj=∑i=1nτijxi2σj32σj2∑i=1nτij=∑i=1nτijxi2.Hence, at the iterationm+1th, the component parameter can be estimated as(13)σ^jm+1=∑i=1nτijmxi22∑i=1nτijm.The algorithm is terminated as soon as we get convergence or after the prespecified number of iterations. Note that the final iterationz^ij can be estimated by generating it from multinomial distribution as(14)z^i∼Mult1,τi.The estimated value ofz^i can be used as a classifier of component members based on the definition mentioned in Equation (1).The elements of the componentj are represented by Sj, where(15)Sj=z^ij=1,i=1,2,…,n.


## 4.1.E-Step

The unobservable data (herezij) are handled by the E-step, which takes the expectation of the complete-data log likelihood conditionally on the given observed data X, using the value of θ of the current iteration. Let θ0 be the value chosen initially for θ. The first iteration of the EM algorithm computes the conditional expectation of log Lθ given X by using θ0 for θ, which can be written as(7)Qθ;θ0=Eθ0logLθ|x.As logθ is linear in the missing data zij, the E-step at the m+1th iteration requires the current expectation of the random variable Zij, which corresponds to zij, to be calculated conditionally on the given observed data X.(8)Eθ0Zij|x=PrλmZij=1|x=τijmxi;θm.Then, we get(9)τijmxi;θm=πjmfjxi;σjm∑j=1kπjmfjxi;σjm.The quantityτjmx;θm is the posterior probability, which implies that the ith element of the sample with the observed value xi belongs to the jth component of the mixture. By substituting zij into Equation (3) with τjmx;θm, the Q function becomes(10)Qθ;θm=∑i=1n∑j=1kτijmxi;θmlogπj+logxi−logσj2+−xi22σj2.


## 4.2.M-Step

To get the updated estimateθm+1, the global maximization of Qθ;θm with respect to θ over the parameter space is required by the M-step on the m+1th iteration. Given the updated τijmxi;θm at the iteration mth, the weighted parameter πjm+1 can be estimated by taking the first derivative of Qθ;θm with respect to πj and equating it to 0, which results in(11)π^jm+1=∑i=1nτijmn.The component parameterσj at the iteration m+1th can be computed by deriving Qθ;θm with respect to σj and equating it to 0, and we get(12)∂Q∂σj=∑i=1n−2σjτijσj2+τijxi2σj3⟶−2∑i=1nτijσj+∑i=1nτijxi2σj3=0⟶2∑i=1nτijσj=∑i=1nτijxi2σj32σj2∑i=1nτij=∑i=1nτijxi2.Hence, at the iterationm+1th, the component parameter can be estimated as(13)σ^jm+1=∑i=1nτijmxi22∑i=1nτijm.The algorithm is terminated as soon as we get convergence or after the prespecified number of iterations. Note that the final iterationz^ij can be estimated by generating it from multinomial distribution as(14)z^i∼Mult1,τi.The estimated value ofz^i can be used as a classifier of component members based on the definition mentioned in Equation (1).The elements of the componentj are represented by Sj, where(15)Sj=z^ij=1,i=1,2,…,n.


## 5. Assessing the Number of Components

The abovementioned EM algorithm is based on known allocationSj. Therefore, we assess the appropriate number of components by calculating the Bayesian information criteria according to the formula [28]:(16)BIC=−2logLθ+dlogn,where Lθ is the maximum likelihood function in (5) and d=2k−1. Observation memberships can be updated after getting Bayesian estimators of z^ij as previously mentioned.


## 6. Simulations

To examine the performance of the EM algorithm on the finite mixture of Rayleigh distribution, we generate data by following two-component and three-component scenarios of the mixture model. The following scenarios were used to generatexj according to predefined parameters θ=π1,θ1,θ2T, where π1,θ1,θ2=0.35,1,4,0.35,0.5,2,0.25,2,3. And θ=π1,π2,θ1,θ2,θ3T, where (π1,π2,θ1,θ2,θ3=0.3,0.6,1,4,7,0.3,0.6,0.5,2,4,0.2,0.5,1,2,3. The sample sizes n=90,120,200. A comparison was made by calculating the average mean square error (AMSE) of the estimated model and the average classification success rate (ACSR) for the number of replications equal to R by using the following formulas: (17)MSE=∑i=1nfxi;θ−fxi;θ^2n,(18)AMSE=∑r=1RMSErR.Here,MSEr is the MSE of the model on the replication nth.For calculating ACSR, we letAr be the number of items with the correctly estimated membership divided by the sample size n on the replication rth.(19)ACSR=∑r=1RArR.


## 7. Results

The results in Table1 represent the data generated from the two-component mixture model and show that the EM algorithm is good at estimating the parameters of the model with respect to AMSE, which means the estimated parameters are close to the real ones used in the generated data. It can be seen that the AMSE is decreased as sample sizes are increased. However, the ACSR is decreased when component parameters are close to each other. The graphs of the PDF of the model with respect to true parameters and estimators are shown in Figure 1, which shows the curves of probability distribution functions given the estimated parameters compared to the one of the true parameters values that has been used in the simulation. The results of applying the EM algorithm to the three-component model are shown in Table 1.Table 1
The results of AMSE and ACSR of applying the EM algorithm to simulated data from the two-component mixture model for three different scenarios and three different sample sizes with 1000 replications.
Nπ1=0.35,π2=0.65, θ1=1,θ2=4π1=0.35,π2=0.65, θ1=0.5,θ2=2π1=0.25,π2=0.75, θ1=2,θ2=3AMSEACSRAMSEACSRAMSEACSR900.0021120.82900.0074570.77240.0021090.67811200.0016400.84040.0065230.78340.0019230.69632000.001550.8750.0063120.79320.0017220.7261Figure 1
The graphs of the PDF of the model for true parameters and estimated ones by the EM algorithm. The columns of graphs represent the sample sizesn = 90, 120, and 200, respectively. The rows of graphs represent the parameter scenarios (π1, θ1, θ2) = (0.35, 1, 4), (0.35, 0.5, 2), and (0.25, 2, 3), respectively.The results of applying the EM algorithm to the three-component model are shown in Table2. It can be seen that the results of AMS of the large sample sizes are better than the small sample sizes. However, it is obvious that as the component parameters are close to each other, the ACSR is decreased, which means that the clustering of data becomes difficult. This is obvious in Figure 2, where the of the probability distribution functions of the third scenario does not show the distance between the tops.Table 2
The results of AMSE and ACSR of applying the EM algorithm to simulated data from the three-component mixture model for three different scenarios and three different sample sizes with 1000 replications.
Nπ1=0.3,π2=0.3,π3=0.4, θ1=1,θ2=4,θ3=7π1=0.3,π2=0.3,π3=0.4, θ1=0.5,θ2=2,θ3=4π1=0.2,π2=0.3,π3=0.5, θ1=1,θ2=2,θ3=3AMSEACSRAMSEACSRAMSEACSR900.0044270.695050.013160.631540.013880.562361200.004160.72360.012870.657090.013630.579032000.0040070.748380.012390.675960.013110.5921Figure 2
The graphs of the PDF of the model for true parameters and estimated ones by the EM algorithm. The columns of the graphs represent the sample sizesn = 90, 120, and 200, respectively. The rows of the graphs represent the parameter scenarios (π1, π2, θ1, θ2, θ3) = (0.3, 0.3, 1, 4, 7), (0.3, 0.3, 0.5, 2, 4), and (0.2, 0.3, 1, 2, 3), respectively.


## 8. Conclusion and Discussion

We derived the formulas for estimating the parameters of the finite mixture of Rayleigh distribution by the EM algorithm in a general form. The mathematical formulas of estimators are valid for thek-component mixture of Rayleigh distribution. However, the computations of the EM algorithm can be more accurate when sample sizes are increased. The underlying reason for this is the computation of the member variable as it is calculated with the E-step of the algorithm, given the observations. This in turn will result in great accuracy when it comes to estimating component parameters and weights. We apply the proposed framework to simulated data. The data were simulated according to two and three components as a special case from the general one. It can be seen from the results that the algorithm is performed in a good manner when sample sizes are large and parameters are not close enough to each other. We considered the AMSE and ACSR to be tools for measuring the accuracy of the estimation and clustering of data into their groups. The method was illustrated by graphs as well to show the estimated graph and compare it to one of the true parameters that are used in generating data for several scenarios.



---

*Source: 7596449-2022-12-30.xml*