# Improved Expectation Maximization Algorithm for Gaussian Mixed Model Using the Kernel Method

Authors: Mohd Izhan Mohd Yusoff; Ibrahim Mohamed; Mohd Rizam Abu Bakar
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Fraud activities have contributed to heavy losses suffered by telecommunication companies. In this paper, we attempt to use Gaussian mixed model, which is a probabilistic model normally used in speech recognition to identify fraud calls in the telecommunication industry. We look at several issues encountered when calculating the maximum likelihood estimates of the Gaussian mixed model using an Expectation Maximization algorithm. Firstly, we look at a mechanism for the determination of the initial number of Gaussian components and the choice of the initial values of the algorithm using the kernel method. We show via simulation that the technique improves the performance of the algorithm. Secondly, we developed a procedure for determining the order of the Gaussian mixed model using the log-likelihood function and the Akaike information criteria. Finally, for illustration, we apply the improved algorithm to real telecommunication data. The modified method will pave the way to introduce a comprehensive method for detecting fraud calls in future work.


---

## Body

## 1. Introduction

Every year telecommunication companies register loses amounting to millions of dollars due to fraud activities. Examples of such activities given by [1, 2] include the use of the customer’s line in Premium Rate Service without their knowledge or autodialers with no intention to pay for the outgoing calls; PABX for international calls; an unregistered user with an assigned number accessing the network (such activity is called stolen line unknown); and international roaming manipulation. Vendors, seeing the above as an opportunity not to be missed, compete to provide data mining applications that can detect the said activities effectively using various methods such as OLAP, deviation based outlier detection, and hidden Markov model. The focus of this paper is Gaussian mixed model, henceforth, GMM.A GMM is best known for providing a robust speaker representation for the difficult task of speaker identification on short-time speechspectra [3]. Its function is extended to detect fraud activities on the number (as well as length) of domestic and international calls made on a daily basis during office, evening, and night hours. Reference [4] presented three approaches to fraud detection in communication networks: neural networks with supervised learning, probability density estimation methods, and Bayesian networks. Information describing a subscriber’s behavior kept in toll tickets was used. For example, supervised learning used summary statistics over the whole observed time period (especially the number of times fraud activities were recorded in the data). The two latter approaches used a subscriber’s daily behavior. To improve the fraud detection system, they recommended the combination of the three presented methods together with the incorporation of rule based systems.The maximum likelihood estimation for a GMM is generally difficult to obtain directly, but it is made easier with the availability of the Expectation Maximization (EM) algorithm which was first introduced by [5]. Since then, there has been a significant increase in its use especially in finding the maximum likelihood for probabilistic models. For example, [6, 7] developed an online system for detecting fraud calls using a hierarchical switching generative model. The model is trained by using the EM algorithm on an incomplete data set and is further improved by using a gradient-based discriminative method. In this paper, we propose an improved EM algorithm for GMM in detecting fraud calls in telecommunication.This paper is organized as follows. Section2 gives an introduction to GMM. We then describe the EM algorithm for a GMM, the kernel method, and eventually the proposed modified EM algorithm for GMM in Section 3. In Section 4, we study the performance of the modified algorithm in estimating the parameters and the effect of overlapping areas of Gaussian components in GMM. In the next section, we propose graphical plots to identify the “best” number of components in a GMM. For illustration, an application of the improvement on a real data set is presented in Section 6.


## 2. Gaussian Mixed Model

Letx∈Rd and K be the number of components where each component has its own prior probability ai and probability density function with mean μi and covariance Σi,i=1,…,K. A Gaussian mixed model is then given by
(1)∑i=1Kaiϕ(x∣μi,Σi)=∑i=1Kai1(2π)d|Σi|exp(-(x-μi)tΣi-1(x-μi)2),
where ∑i=1Kai=1. We next define the likelihood function and the log-likelihood function by L(X∣θ)=∏j=1nf(xj∣θ) and l(X∣θ)=∑j=1nlog(∑i=1Kaiϕ(xj∣μi,Σi)) where X=(x1t,…,xnt)t, respectively. The maximum likelihood estimation (m.l.e) method aims at finding θ^ that maximizes l(X∣θ); see [8]. The expression log(∑i=1Kaiϕ(xj∣μi,Σi)) in l(X∣θ) is difficult to compute. We use the Expectation Maximization (EM) algorithm to overcome this problem.


## 3. Expectation Maximization (EM) Algorithm

### 3.1. EM for Gaussian Mixed Model

In a general setup of the EM algorithm given in [5], the authors considered an unobservable variable X in sample space 𝒳, which is indirectly observed through observed variable Y in sample space 𝒴. Assuming that f(x∣θ) is the sampling density depending on the parameter θ∈Ω, the corresponding family of sampling densities for Y, say g(y∣θ), can be derived from
(2)g(y∣θ)=∫χ(y)f(x∣θ)dx,
where χ(y) is a subset of 𝒳 under the mapping x→y(x) from 𝒳 to 𝒴. The main objective of the EM algorithm is to find the value of θ that maximizes (2). Consider the expected value of logf(x∣θ′) given y and θ, denoted by Q(θ′∣θ), where
(3)Q(θ′∣θ)=E(logf(x∣θ′)∣y,θ)
with the expectation assumed to exist for all pairs (θ′,θ) and f(x∣θ)>0 for θ∈Ω. According to [5], the EM iteration consists of two steps, namely, the E-step and the M-step. At the pth iteration with the estimate of θ denoted by θ(p), the E-step will give the value of Q(θ∣θ(p)) and the M-step will find a new estimate of θ, say θ(p+1), that maximizes Q(θ∣θ(p)). The steps are repeated until convergence is achieved.For the case of a GMM, we defineQ(θ′∣θ)=E[log∏i=1nayi′ϕ(xi∣μyi′,Σyi′)∣X,θ], where yi∈{1,2,…,K} and yi=k, if the ith sample is generated by the kth mixture component. It is simplified, by applying, amongst others, the Bayes formula f(θ∣x)∝f(x∣θ)P(θ) where f(θ∣x) is the posterior probability, f(x∣θ) is the likelihood function, and P(θ) is the prior probability to the following equations (see [9, 10]):
(4)Q(θ′∣θ)=∑i=1n∑k=1Kpi,klogak′+∑i=1n∑k=1Kpi,klogϕ(xi∣μk′,Σk′),
where
(5)pi,k=akϕ(xi∣μk,Σk)∑lalϕ(xi∣μl,Σl),(6)ϕ(xi∣μk,Σk)=1(2π)d|Σk|exp(-(xi-μk)tΣk-1(xi-μk)2).
Hence, the EM iteration for a GMM is defined by the following.E-step: use (5).M-step: use the formulas
(7)aj=1n∑ipij,μj=∑ipijxi∑ipij,Σj=∑ipij(xi-μj)(xi-μj)t∑ipij
which are derived from the Lagrange multipliers, ∂Q/∂μj=0 and ∂Q/∂Σj-1=0 (for details, see the appendices).The above steps (i.e.,E-step and M-step) are repeated until convergence is achieved.


### 3.2. The Kernel Method

The kernel method can be used to find the probability density estimate for univariate data; see, for example, [11]. Let α<min(xi)-3h, β>max(xi)+3h, M=2r; let h be the bandwidth for some integer r, δ=(β-α)/M; and let tk=α+kδ be the kth grid point where k=0,1,…,M-1. The density estimate at grid point tk is represented by the following equation:
(8)f^(tk)≈∑l=-M/2M/2(1M∑k=0M-1ξkexp(2πklMi))×exp((-2πklMi)-12h2(2πlβ-α)2),
where i2=-1. For x∈[tk,tk+1], the density estimate f^(x) is defined by f^(x)=(1/nh)∑i=1nK((x-xi)/h) where K(t)=(1/2π)exp(-(1/2)t2). To compute f^(x) at a grid of points, a method which makes use of the Fourier transform is employed. Let f~(s) be the Fourier transform of the kernel density estimate f^(x). It can be shown that f~(s)=(2π)1/2K~(hs)u(s)=exp(-(1/2)h2s2)u(s), where K~(s) is the Fourier transform of the Gaussian kernel and u(s)=(2π)-1/2n-1∑j=1nexp(isxj) is the Fourier transform of the data. Thus, f^(x)=(2π)-1/2∫e-isx(2π)1/2K~(hs)u(s)ds is the convolution of the data with the kernel.We will use the following algorithm by [11] to discretize the data to very fine grids and to find f^(x) by convolving the data with the kernel.Step A. Discretize the data to find the weight sequence {ξk} with M=28. If x∈[tk,tk+1], it is split into a weight (1/nδ2)(tk+1-x) at tk and a weight (1/nδ2)(x-tk) at tk+1; these weights are accumulated over all the data points xi to give a sequence of (ξk) weights summing up to 1/δ.Step B. Find the sequence {Yl} defined by Yl=M-1∑k=0M-1ξkexp((2πkl/M)i), where -(M/2)≤l≤M/2. It can be shown that when α=0, Yl≈(2π)1/2(β-α)-1u(sl), where sl=2πl/(β-α).Step C. Find the sequence {ζl*}, where ζl*=exp(-(1/2)h2sl2)Yl. Here, h=0.9An-1/5, where A=min(sd,IQR/1.34), sd is the standard deviation, and IQR is the interquartile range. The IQR is chosen here by [11], who claimed that the bandwidth is useful for a wide range of densities.Step D. Let ζk be the inverse discrete Fourier transform of ζl*, that is, ζk=∑l=-M/2M/2ζl*exp(-(2πkl/M)i).It can be shown that whenα=0, f^(tk)≈ζk. We then identify xi where its density estimate, denoted by f^(xi), is greater than those of its nearest neighbors xi+1 and xi-1. In other words, f^(xi)>f^(xi+1) and f^(xi)>f^(xi-1); refer to Figure 1 where the vertical line that touches tk and f^(tk) shows the location of the peak. Note that we may obtain more than one maximum points, which means that the data may consist of more than one Gaussian distribution. These results form a very important component of the improved EM algorithm for GMM to be described next.Figure 1
Plot off^(tk) against tk.


### 3.3. Improved EM Algorithm for GMM

A number of authors highlighted the importance of identifying the right number, sayk, of components in a GMM and subsequently choosing good initial values for the model parameters μi and σi2, i=1,2,…k, in the EM algorithm. Reference [12] noted the difficulty of using log-likelihood-ratio statistics to test the number of components and subsequently suggested using a nonparametric bootstrapping approach. Similarly, [13] pointed out the same concerns and introduced an algorithm called the stepwise-split-and-merge EM algorithm to solve the said problem. In addition, [14] investigated the possibility of using the minimization of the Kullback-Leiber distance between fitted mixture model and the true density as a method for estimating k where the said distance was estimated using cross validation. Reference [15] viewed the mixture distribution as a contaminated Gaussian density and proposed a recursive algorithm called the Gaussian Mixture Density Decomposition algorithm for identifying each Gaussian component in the mixture. Other works on this topic can also be found, for example, in [16, 17].In this paper, we propose an improved EM algorithm for GMM which can perform both tasks: identifying the initial number of components and providing automatic initial values for the EM algorithm. The full improved EM algorithm for a GMM is now presented.Step 1.
The kernel method as described in Section3.2 is used to determine the number, say K0, of components and also the corresponding means μi of each component, where i=1,2,…,K0. The initial estimates of the standard deviations σii are set to unity while the prior weights ai are set to be 1/K0.Step 2.
The EM algorithm for a GMM as described in Section3.1 is executed to give the final estimates of parameters μi, σii, and ai, i=1,2,…,K0. The log-likelihood function and Akaike information criteria (AIC) are calculated using the said parameters.Step 3.
Step2 is repeated for other possible number K of components with μi=0, σii=1 for the other K-K0  components, and ai=1/K.Step 4.
The log-likelihood function and AIC values forK=1,2,…,10 are plotted. The final number of components Kf is chosen when adding extra components in the model does not significantly increase or decrease the values of the log-likelihood function and the AIC, respectively.


## 3.1. EM for Gaussian Mixed Model

In a general setup of the EM algorithm given in [5], the authors considered an unobservable variable X in sample space 𝒳, which is indirectly observed through observed variable Y in sample space 𝒴. Assuming that f(x∣θ) is the sampling density depending on the parameter θ∈Ω, the corresponding family of sampling densities for Y, say g(y∣θ), can be derived from
(2)g(y∣θ)=∫χ(y)f(x∣θ)dx,
where χ(y) is a subset of 𝒳 under the mapping x→y(x) from 𝒳 to 𝒴. The main objective of the EM algorithm is to find the value of θ that maximizes (2). Consider the expected value of logf(x∣θ′) given y and θ, denoted by Q(θ′∣θ), where
(3)Q(θ′∣θ)=E(logf(x∣θ′)∣y,θ)
with the expectation assumed to exist for all pairs (θ′,θ) and f(x∣θ)>0 for θ∈Ω. According to [5], the EM iteration consists of two steps, namely, the E-step and the M-step. At the pth iteration with the estimate of θ denoted by θ(p), the E-step will give the value of Q(θ∣θ(p)) and the M-step will find a new estimate of θ, say θ(p+1), that maximizes Q(θ∣θ(p)). The steps are repeated until convergence is achieved.For the case of a GMM, we defineQ(θ′∣θ)=E[log∏i=1nayi′ϕ(xi∣μyi′,Σyi′)∣X,θ], where yi∈{1,2,…,K} and yi=k, if the ith sample is generated by the kth mixture component. It is simplified, by applying, amongst others, the Bayes formula f(θ∣x)∝f(x∣θ)P(θ) where f(θ∣x) is the posterior probability, f(x∣θ) is the likelihood function, and P(θ) is the prior probability to the following equations (see [9, 10]):
(4)Q(θ′∣θ)=∑i=1n∑k=1Kpi,klogak′+∑i=1n∑k=1Kpi,klogϕ(xi∣μk′,Σk′),
where
(5)pi,k=akϕ(xi∣μk,Σk)∑lalϕ(xi∣μl,Σl),(6)ϕ(xi∣μk,Σk)=1(2π)d|Σk|exp(-(xi-μk)tΣk-1(xi-μk)2).
Hence, the EM iteration for a GMM is defined by the following.E-step: use (5).M-step: use the formulas
(7)aj=1n∑ipij,μj=∑ipijxi∑ipij,Σj=∑ipij(xi-μj)(xi-μj)t∑ipij
which are derived from the Lagrange multipliers, ∂Q/∂μj=0 and ∂Q/∂Σj-1=0 (for details, see the appendices).The above steps (i.e.,E-step and M-step) are repeated until convergence is achieved.


## 3.2. The Kernel Method

The kernel method can be used to find the probability density estimate for univariate data; see, for example, [11]. Let α<min(xi)-3h, β>max(xi)+3h, M=2r; let h be the bandwidth for some integer r, δ=(β-α)/M; and let tk=α+kδ be the kth grid point where k=0,1,…,M-1. The density estimate at grid point tk is represented by the following equation:
(8)f^(tk)≈∑l=-M/2M/2(1M∑k=0M-1ξkexp(2πklMi))×exp((-2πklMi)-12h2(2πlβ-α)2),
where i2=-1. For x∈[tk,tk+1], the density estimate f^(x) is defined by f^(x)=(1/nh)∑i=1nK((x-xi)/h) where K(t)=(1/2π)exp(-(1/2)t2). To compute f^(x) at a grid of points, a method which makes use of the Fourier transform is employed. Let f~(s) be the Fourier transform of the kernel density estimate f^(x). It can be shown that f~(s)=(2π)1/2K~(hs)u(s)=exp(-(1/2)h2s2)u(s), where K~(s) is the Fourier transform of the Gaussian kernel and u(s)=(2π)-1/2n-1∑j=1nexp(isxj) is the Fourier transform of the data. Thus, f^(x)=(2π)-1/2∫e-isx(2π)1/2K~(hs)u(s)ds is the convolution of the data with the kernel.We will use the following algorithm by [11] to discretize the data to very fine grids and to find f^(x) by convolving the data with the kernel.Step A. Discretize the data to find the weight sequence {ξk} with M=28. If x∈[tk,tk+1], it is split into a weight (1/nδ2)(tk+1-x) at tk and a weight (1/nδ2)(x-tk) at tk+1; these weights are accumulated over all the data points xi to give a sequence of (ξk) weights summing up to 1/δ.Step B. Find the sequence {Yl} defined by Yl=M-1∑k=0M-1ξkexp((2πkl/M)i), where -(M/2)≤l≤M/2. It can be shown that when α=0, Yl≈(2π)1/2(β-α)-1u(sl), where sl=2πl/(β-α).Step C. Find the sequence {ζl*}, where ζl*=exp(-(1/2)h2sl2)Yl. Here, h=0.9An-1/5, where A=min(sd,IQR/1.34), sd is the standard deviation, and IQR is the interquartile range. The IQR is chosen here by [11], who claimed that the bandwidth is useful for a wide range of densities.Step D. Let ζk be the inverse discrete Fourier transform of ζl*, that is, ζk=∑l=-M/2M/2ζl*exp(-(2πkl/M)i).It can be shown that whenα=0, f^(tk)≈ζk. We then identify xi where its density estimate, denoted by f^(xi), is greater than those of its nearest neighbors xi+1 and xi-1. In other words, f^(xi)>f^(xi+1) and f^(xi)>f^(xi-1); refer to Figure 1 where the vertical line that touches tk and f^(tk) shows the location of the peak. Note that we may obtain more than one maximum points, which means that the data may consist of more than one Gaussian distribution. These results form a very important component of the improved EM algorithm for GMM to be described next.Figure 1
Plot off^(tk) against tk.


## 3.3. Improved EM Algorithm for GMM

A number of authors highlighted the importance of identifying the right number, sayk, of components in a GMM and subsequently choosing good initial values for the model parameters μi and σi2, i=1,2,…k, in the EM algorithm. Reference [12] noted the difficulty of using log-likelihood-ratio statistics to test the number of components and subsequently suggested using a nonparametric bootstrapping approach. Similarly, [13] pointed out the same concerns and introduced an algorithm called the stepwise-split-and-merge EM algorithm to solve the said problem. In addition, [14] investigated the possibility of using the minimization of the Kullback-Leiber distance between fitted mixture model and the true density as a method for estimating k where the said distance was estimated using cross validation. Reference [15] viewed the mixture distribution as a contaminated Gaussian density and proposed a recursive algorithm called the Gaussian Mixture Density Decomposition algorithm for identifying each Gaussian component in the mixture. Other works on this topic can also be found, for example, in [16, 17].In this paper, we propose an improved EM algorithm for GMM which can perform both tasks: identifying the initial number of components and providing automatic initial values for the EM algorithm. The full improved EM algorithm for a GMM is now presented.Step 1.
The kernel method as described in Section3.2 is used to determine the number, say K0, of components and also the corresponding means μi of each component, where i=1,2,…,K0. The initial estimates of the standard deviations σii are set to unity while the prior weights ai are set to be 1/K0.Step 2.
The EM algorithm for a GMM as described in Section3.1 is executed to give the final estimates of parameters μi, σii, and ai, i=1,2,…,K0. The log-likelihood function and Akaike information criteria (AIC) are calculated using the said parameters.Step 3.
Step2 is repeated for other possible number K of components with μi=0, σii=1 for the other K-K0  components, and ai=1/K.Step 4.
The log-likelihood function and AIC values forK=1,2,…,10 are plotted. The final number of components Kf is chosen when adding extra components in the model does not significantly increase or decrease the values of the log-likelihood function and the AIC, respectively.


## 4. Simulation

We use simulation to investigate the performance of the proposed modified algorithm.


### 4.1. Simulation Scheme

Simulation data were generated using the Box and Muller Transformation [18] as defined by (9) below:
(9)zj=μ+(-2σ2loguj)1/2cos2πuj+1,zj+1=μ+(-2σ2loguj)1/2sin2πuj+1,
where uj,uj+1~U(0,1). For the case of two components, we start by generating a random number u1~U(0,1). If 0<u1<a1, we generate two random numbers u2~U(0,1) and u3~U(0,1) and calculate z2+z3 using the first and second equations of (9) with μ*=μ1/2 and σ*=σ1/2. Otherwise, we use μ*=μ2/2, and σ*=σ2/2. The process continues until the required sample size is obtained. The scheme is easily extended to any number of components. For further details, refer to [19].


### 4.2. Study of Performance Based on the Log-Likelihood Function

We first look at the performance of the standard method, called Method 1, followed by that of the modified method, called Method 2. For Method 1, in place of Step1 of the modified method, we assign values zero and unity, respectively, to the means and variances of all components. We compare the performance by looking at the log-likelihood function via simulation study.Following [20], we consider two cases with two and one case with three components with the true values of the parameters given in Table 1. For each case, we generate 100 samples of size 1000 where the chosen sample size reflects the large size of data sets found in the telecommunication industry, the focus of our interest. We then apply Method 1 and Method 2 on the simulated data. For each case, for better quality viewing, we plot only 50 values of the log-likelihood function for both methods on the same plot, as given in Figure 2. It can be seen that, for Samples 1 and 3, the proposed Method 2 clearly outperforms the standard Method 1 with the values of the log-likelihood function corresponding to Method 2 being always larger than those of Method 1. However, we see that some values overlap for Sample 2, though the proposed Method 2 still generally performs better. In this case, the prior probabilities ai are distinctly different from the chosen values of ai in Sample 1 while other true values remain the same which leads to different percentages of overlapping of the Gaussian components in the GMM. Hence, we will investigate the performance of the improved EM algorithm in estimating the parameters of the GMM by taking into account the effect of different percentages of overlapping between the components observed in the data.Table 1
List of true values ofa’s, and
μ’s, σ’s.
Sample name and size (in bracket)
Prior probability
Mean
Variance
Sample 1Two components
a
1
=
0.4 a2=0.6
μ
1
=
0.0 μ2=2.0
σ
1
2
=
1.0 σ22=0.25
Sample 2Two components
a
1
=
0.85 a2=0.15
μ
1
=
0.0 μ2=2.0
σ
1
2
=
1.0 σ22=0.25
Sample 3Three components
a
1
=
0.33 a2=0.33, a3=0.34
μ
1
=
0.0 μ2=-1.0, μ3=4.0
σ
1
2
=
1.0 σ22=0.25, σ32=4.0Plots of values of log-likelihood function.
(a)
Sample 1
(b)
Sample 2
(c)
Sample 3


### 4.3. The Effects of Different Overlapping Percentages on Performance

The main objective here is to investigate the performance of the modified EM algorithm for different overlapping percentages of the components in the GMM. For simplicity, we restrict our attention to two components so thatθ=(θ1,…,θ6)=(a1,a2,μ1,μ2,σ11,σ22) are to be estimated. Data is simulated using the simulation scheme described in Section 4.1.After performing Steps1 and 2, we find Di=θi-θ^i, where θi is the true value of the ith parameter and θ^i is the EM estimate of the parameter, i=1,2,…,n. The sample mean and standard deviation of Di are computed using the formulas D-=(1/n)∑i=1nDi and SD=1/(n-1)∑i=1n(Di-D-)2. The estimates are considered good if D- is close to zero, indicating small biases observed in the simulation results, and SD is also close to zero, indicating that the parameter estimates are concentrated around their respective true values.We determine the area of overlapping between the two components for each model by using the misclassification concept given in [21], the details of which are provided in Appendix B. The formula to estimate the overlapping areas depends on the mean and standard deviation of the components. The effects of prior probabilities should not affect the estimates greatly as their sum equals unity.We consider three cases for different combinations of parameterθ which give different percentages of overlapping of the GMM components. The results are tabulated in Tables 2–4. Table 2 deals with case 1, where the true values of μ1=0, μ2=3.0, and σ11=σ22=0.316 are fixed but the true values of a1 and a2 are varied. In all cases, the percentage of overlapping is 0% as the separation of the means is rather large with small values of dispersion. We can see that the values of the mean are close to zero with the small standard errors less than unity for all parameters considered. On the other hand, Table 3 gives the results for case 2 where μ1=0, μ2=1.0, σ22=0.707, and σ11=0.447 are fixed but a1 and a2 are varied to give 25% of overlapping. The bias is still considered small but generally larger than that for case 1. In addition, the values are also more dispersed here. Finally Table 4 shows the results of case 3 where μ1=0, μ2=0.25, σ11=0.577, and σ22=1.414 are fixed with 45% of overlapping. As expected, the results deteriorate when the percentage of overlapping increases. We conclude that the modified EM algorithm for GMM performs well when the percentages of overlapping are small, but its performance is affected when the percentages increase. More comprehensive simulation results can be obtained from the authors upon request.Table 2
Simulation results for the caseμ1=0, μ2=3.0, and σ12=σ22=0.316.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
-
0.001
0.010
0.001
0.010
0.007
0.033
-
0.003
0.013
0.000
0.014
0.001
0.004
0.2
0.8
0.002
0.013
-
0.002
0.013
0.002
0.022
-
0.003
0.011
0.000
0.011
0.002
0.007
0.3
0.7
-
0.002
0.014
0.002
0.014
0.004
0.017
-
0.002
0.010
0.002
0.009
0.001
0.005
0.4
0.6
0.003
0.020
-
0.003
0.020
-
0.004
0.019
-
0.006
0.012
0.000
0.009
0.001
0.005Table 3
Simulation results for the caseμ1=0, μ2=1.0, σ22=0.707, and σ12=0.447.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
-
0.112
0.148
0.112
0.148
-
0.130
0.203
-
0.067
0.139
-
0.030
0.103
0.034
0.068
0.2
0.8
-
0.014
0.073
0.014
0.073
-
0.006
0.087
0.006
0.091
0.022
0.056
0.005
0.057
0.3
0.7
0.020
0.087
-
0.020
0.087
0.031
0.069
0.023
0.104
0.044
0.048
-
0.006
0.075
0.4
0.6
0.067
0.075
-
0.067
0.075
-
0.112
0.444
0.145
0.232
-
0.002
0.085
-
0.021
0.099Table 4
Simulation results for the caseμ1=0, μ2=0.25, σ12=0.577, and σ22=1.414.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
0.089
0.007
-
0.089
0.007
-
0.817
3.770
0.045
0.048
-
0.341
0.823
0.305
0.095
0.2
0.8
0.157
0.108
-
0.157
0.108
0.291
3.521
0.075
0.073
-
0.169
0.374
0.413
0.256
0.3
0.7
0.237
0.100
-
0.237
0.100
-
0.205
2.924
0.105
0.121
-
0.278
0.389
0.578
0.278
0.4
0.6
0.245
0.187
-
0.245
0.187
0.602
2.520
0.092
0.108
-
0.008
0.258
0.508
0.435


## 4.1. Simulation Scheme

Simulation data were generated using the Box and Muller Transformation [18] as defined by (9) below:
(9)zj=μ+(-2σ2loguj)1/2cos2πuj+1,zj+1=μ+(-2σ2loguj)1/2sin2πuj+1,
where uj,uj+1~U(0,1). For the case of two components, we start by generating a random number u1~U(0,1). If 0<u1<a1, we generate two random numbers u2~U(0,1) and u3~U(0,1) and calculate z2+z3 using the first and second equations of (9) with μ*=μ1/2 and σ*=σ1/2. Otherwise, we use μ*=μ2/2, and σ*=σ2/2. The process continues until the required sample size is obtained. The scheme is easily extended to any number of components. For further details, refer to [19].


## 4.2. Study of Performance Based on the Log-Likelihood Function

We first look at the performance of the standard method, called Method 1, followed by that of the modified method, called Method 2. For Method 1, in place of Step1 of the modified method, we assign values zero and unity, respectively, to the means and variances of all components. We compare the performance by looking at the log-likelihood function via simulation study.Following [20], we consider two cases with two and one case with three components with the true values of the parameters given in Table 1. For each case, we generate 100 samples of size 1000 where the chosen sample size reflects the large size of data sets found in the telecommunication industry, the focus of our interest. We then apply Method 1 and Method 2 on the simulated data. For each case, for better quality viewing, we plot only 50 values of the log-likelihood function for both methods on the same plot, as given in Figure 2. It can be seen that, for Samples 1 and 3, the proposed Method 2 clearly outperforms the standard Method 1 with the values of the log-likelihood function corresponding to Method 2 being always larger than those of Method 1. However, we see that some values overlap for Sample 2, though the proposed Method 2 still generally performs better. In this case, the prior probabilities ai are distinctly different from the chosen values of ai in Sample 1 while other true values remain the same which leads to different percentages of overlapping of the Gaussian components in the GMM. Hence, we will investigate the performance of the improved EM algorithm in estimating the parameters of the GMM by taking into account the effect of different percentages of overlapping between the components observed in the data.Table 1
List of true values ofa’s, and
μ’s, σ’s.
Sample name and size (in bracket)
Prior probability
Mean
Variance
Sample 1Two components
a
1
=
0.4 a2=0.6
μ
1
=
0.0 μ2=2.0
σ
1
2
=
1.0 σ22=0.25
Sample 2Two components
a
1
=
0.85 a2=0.15
μ
1
=
0.0 μ2=2.0
σ
1
2
=
1.0 σ22=0.25
Sample 3Three components
a
1
=
0.33 a2=0.33, a3=0.34
μ
1
=
0.0 μ2=-1.0, μ3=4.0
σ
1
2
=
1.0 σ22=0.25, σ32=4.0Plots of values of log-likelihood function.
(a)
Sample 1
(b)
Sample 2
(c)
Sample 3


## 4.3. The Effects of Different Overlapping Percentages on Performance

The main objective here is to investigate the performance of the modified EM algorithm for different overlapping percentages of the components in the GMM. For simplicity, we restrict our attention to two components so thatθ=(θ1,…,θ6)=(a1,a2,μ1,μ2,σ11,σ22) are to be estimated. Data is simulated using the simulation scheme described in Section 4.1.After performing Steps1 and 2, we find Di=θi-θ^i, where θi is the true value of the ith parameter and θ^i is the EM estimate of the parameter, i=1,2,…,n. The sample mean and standard deviation of Di are computed using the formulas D-=(1/n)∑i=1nDi and SD=1/(n-1)∑i=1n(Di-D-)2. The estimates are considered good if D- is close to zero, indicating small biases observed in the simulation results, and SD is also close to zero, indicating that the parameter estimates are concentrated around their respective true values.We determine the area of overlapping between the two components for each model by using the misclassification concept given in [21], the details of which are provided in Appendix B. The formula to estimate the overlapping areas depends on the mean and standard deviation of the components. The effects of prior probabilities should not affect the estimates greatly as their sum equals unity.We consider three cases for different combinations of parameterθ which give different percentages of overlapping of the GMM components. The results are tabulated in Tables 2–4. Table 2 deals with case 1, where the true values of μ1=0, μ2=3.0, and σ11=σ22=0.316 are fixed but the true values of a1 and a2 are varied. In all cases, the percentage of overlapping is 0% as the separation of the means is rather large with small values of dispersion. We can see that the values of the mean are close to zero with the small standard errors less than unity for all parameters considered. On the other hand, Table 3 gives the results for case 2 where μ1=0, μ2=1.0, σ22=0.707, and σ11=0.447 are fixed but a1 and a2 are varied to give 25% of overlapping. The bias is still considered small but generally larger than that for case 1. In addition, the values are also more dispersed here. Finally Table 4 shows the results of case 3 where μ1=0, μ2=0.25, σ11=0.577, and σ22=1.414 are fixed with 45% of overlapping. As expected, the results deteriorate when the percentage of overlapping increases. We conclude that the modified EM algorithm for GMM performs well when the percentages of overlapping are small, but its performance is affected when the percentages increase. More comprehensive simulation results can be obtained from the authors upon request.Table 2
Simulation results for the caseμ1=0, μ2=3.0, and σ12=σ22=0.316.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
-
0.001
0.010
0.001
0.010
0.007
0.033
-
0.003
0.013
0.000
0.014
0.001
0.004
0.2
0.8
0.002
0.013
-
0.002
0.013
0.002
0.022
-
0.003
0.011
0.000
0.011
0.002
0.007
0.3
0.7
-
0.002
0.014
0.002
0.014
0.004
0.017
-
0.002
0.010
0.002
0.009
0.001
0.005
0.4
0.6
0.003
0.020
-
0.003
0.020
-
0.004
0.019
-
0.006
0.012
0.000
0.009
0.001
0.005Table 3
Simulation results for the caseμ1=0, μ2=1.0, σ22=0.707, and σ12=0.447.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
-
0.112
0.148
0.112
0.148
-
0.130
0.203
-
0.067
0.139
-
0.030
0.103
0.034
0.068
0.2
0.8
-
0.014
0.073
0.014
0.073
-
0.006
0.087
0.006
0.091
0.022
0.056
0.005
0.057
0.3
0.7
0.020
0.087
-
0.020
0.087
0.031
0.069
0.023
0.104
0.044
0.048
-
0.006
0.075
0.4
0.6
0.067
0.075
-
0.067
0.075
-
0.112
0.444
0.145
0.232
-
0.002
0.085
-
0.021
0.099Table 4
Simulation results for the caseμ1=0, μ2=0.25, σ12=0.577, and σ22=1.414.
Prior prob.
Bias,Di
a
1
a
2
μ
1
μ
2
σ
1
2
σ
2
2
a
1
a
2
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
D
-
S
D
0.1
0.9
0.089
0.007
-
0.089
0.007
-
0.817
3.770
0.045
0.048
-
0.341
0.823
0.305
0.095
0.2
0.8
0.157
0.108
-
0.157
0.108
0.291
3.521
0.075
0.073
-
0.169
0.374
0.413
0.256
0.3
0.7
0.237
0.100
-
0.237
0.100
-
0.205
2.924
0.105
0.121
-
0.278
0.389
0.578
0.278
0.4
0.6
0.245
0.187
-
0.245
0.187
0.602
2.520
0.092
0.108
-
0.008
0.258
0.508
0.435


## 5. Determination of the Final Number of Components in the GMM

In the last two steps of the modified algorithm, we intend to confirm that the choice of the initial numberK0 of components in the GMM using kernel method is final. This can be done by considering extra components in the model. For that, as stated in Section 3.3, we repeat Step 2 for other possible number K of components, by setting μi=0, σii=1 for the other K-K0 components, and ai=1/K. The final number of components Kf is determined when adding extra component neither increases the log-likelihood nor decreases the AIC values significantly. The changes can easily be seen on a line plot of the values.


## 6. Real Example: Phone Call Data

The call detail record, which was supplied by Telekom Malaysia Berhad (henceforth, TM), consists of calls made by customers that fell victim to fraud activities. Table5 shows the format of the call detail record for each TM’s customer. We performed several steps on the original data in order to have the data in a desired format, that is, group the real data according to service number, find the country that matches with the country code as well as dialed digits, and sort the real data according to seize time. The column entitled “Seize time” gives the time when the call was made; the fourth column details the duration of the calls in the following format: hour (hh), minute (mm), and second (ss); and the fifth column is the result of converting the information in the fourth column into day format.Table 5
An extract from Customer A’ call detail record.
Service number
Dialed digits
Seize time
Duration (hhmmss)
Converted duration
Xxx
yyy
8:41:37
000339
2.53
E
-
03
Xxx
yyy
9:27:03
000035
4.05
E
-
04
Xxx
yyy
9:43:46
000048
5.56
E
-
04
Xxx
yyy
9:50:21
000031
3.59
E
-
04
Xxx
yyy
10:54:30
000138
1.13
E
-
03We consider real data consisting of the converted duration of each call made by Customer A (referring to the fifth column of Table5), whose identity is not revealed to ensure confidentiality, on March 31, 2011. Seventeen (17) calls were made and the data are displayed in Figure 3. Step 1 of the improved EM algorithm for GMM identifies two initial components. The plots of the log-likelihood function and AIC in Figures 4(a) and 4(b) are the results from performing Steps 2, 3, and 4 of the improved EM algorithm for GMM, which reveal that the EM algorithm fails to achieve convergence when the number of components equals to five or above. It can also be seen that a GMM with 2 components is identified as the “best” model, since the inclusion of more components not only fails to increase the value of the log-likelihood but also fails to decrease the values of the AIC. The final EM estimates for the two-component GMM are a^1=0.64, a^2=0.36, μ^1=-0.66, μ^2=1.17, σ^11=0.07, and σ^22=0.35, and they represent the behavior of calls made by Customer A on March 31, 2011. In a future paper, we will show how the above information produced from the improved EM algorithm for GMM can be used in the process of detecting fraud activities in the telecommunication industry.Figure 3
Duration (in day format) is displayed in the histogram.Plots of log-likelihood and AIC values.
(a)
Log-likelihood
(b)
AIC


## 7. Conclusion

In this paper, we proposed a modified EM algorithm which can numerically identify the number of components of a GMM and estimate the parameters of the model using the kernel method. We showed via simulation that the performance of the algorithm is generally good but, as expected, is affected by increasing percentages of overlapping of the Gaussian components. We then used the line plots of the log-likelihood and AIC values to identify the final number of GMM components. They could clearly be determined via the concave-like shape of the AIC plot which indicates that the AIC decreases to a minimum value and then increases as the number of components increases. Finally, the modified EM algorithm for GMM was tested on real telecommunication data. The results serve as testimony to the effectiveness of the improved EM algorithm for GMM and should be useful when considering the problem of fraud calls faced by the telecommunication companies.



---

*Source: 757240-2013-08-22.xml*