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

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.

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. Section 2 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

Let 𝐱∈𝑅^𝑑 and 𝐾 be the number of components where each component has its own prior probability 𝑎_𝑖 and probability density function with mean 𝝁_𝑖 and covariance Σ_𝑖,𝑖=1,…,𝐾. A Gaussian mixed model is then given by ^𝐾_𝑖=1𝑎_𝑖𝜙𝐱∣𝝁_𝑖,Σ_𝑖=^𝐾_𝑖=1𝑎_𝑖1(2𝜋)^𝑑Σ_𝑖−exp𝐱−𝝁_𝑖^𝑡Σ_𝑖⁻¹𝐱−𝝁_𝑖2, (1) where ∑^𝐾_𝑖=1𝑎_𝑖=1. We next define the likelihood function and the log-likelihood function by 𝐿(𝐗∣𝜃)=∏^𝑛_𝑗=1𝑓(𝐱_𝑗∣𝜃) and 𝑙(𝐗∣𝜃)=∑^𝑛_𝑗=1log(∑^𝐾_𝑖=1𝑎_𝑖𝜙(𝐱_𝑗∣𝝁_𝑖,Σ_𝑖)) where 𝐗=(𝐱^𝑡₁,…,𝐱^𝑡_𝑛)^𝑡, respectively. The maximum likelihood estimation (m.l.e) method aims at finding 𝜃 that maximizes 𝑙(𝐗∣𝜃); see [8]. The expression log(∑^𝐾_𝑖=1𝑎_𝑖𝜙(𝐱_𝑗∣𝝁_𝑖,Σ_𝑖)) in 𝑙(𝐗∣𝜃) 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 𝑋 in sample space 𝒳, which is indirectly observed through observed variable 𝑌 in sample space 𝒴. Assuming that 𝑓(𝑥∣𝜃) is the sampling density depending on the parameter 𝜃∈Ω, the corresponding family of sampling densities for 𝑌, say 𝑔(𝑦∣𝜃), can be derived from 𝑔𝑦∣𝜃=_𝜒𝑦𝑓(𝑥∣𝜃)𝑑𝑥, (2) where 𝜒(𝑦) is a subset of 𝒳 under the mapping 𝑥→𝑦(𝑥) from 𝒳 to 𝒴. The main objective of the EM algorithm is to find the value of 𝜃 that maximizes (2). Consider the expected value of log𝑓(𝑥∣𝜃^′) given 𝑦 and 𝜃, denoted by 𝑄(𝜃^′∣𝜃), where 𝑄𝜃^′∣𝜃=𝐸log𝑓𝑥∣𝜃^′∣𝑦,𝜃 (3) with the expectation assumed to exist for all pairs (𝜃^′,𝜃) and 𝑓(𝑥∣𝜃)>0 for 𝜃∈Ω. According to [5], the EM iteration consists of two steps, namely, the 𝐸-step and the 𝑀-step. At the 𝑝th iteration with the estimate of 𝜃 denoted by 𝜃^(𝑝), the 𝐸-step will give the value of 𝑄(𝜃∣𝜃^(𝑝)) and the 𝑀-step will find a new estimate of 𝜃, say 𝜃^(𝑝+1), that maximizes 𝑄(𝜃∣𝜃^(𝑝)). The steps are repeated until convergence is achieved.

For the case of a GMM, we define 𝑄(𝜃^′∣𝜃)=𝐸[log∏^𝑛_𝑖=1𝑎^′_𝑦𝑖𝜙(𝐱_𝑖∣𝝁^′_𝑦𝑖,Σ^′_𝑦𝑖)∣𝐗,𝜃], where 𝑦_𝑖∈{1,2,…,𝐾} and 𝑦_𝑖=𝑘, if the 𝑖th sample is generated by the 𝑘th mixture component. It is simplified, by applying, amongst others, the Bayes formula 𝑓(𝜃∣𝑥)∝𝑓(𝑥∣𝜃)𝑃(𝜃) where 𝑓(𝜃∣𝑥) is the posterior probability, 𝑓(𝑥∣𝜃) is the likelihood function, and 𝑃(𝜃) is the prior probability to the following equations (see [9, 10]): 𝑄𝜃^′∣𝜃=^𝑛𝐾_𝑖=1𝑘=1𝑝_𝑖,𝑘log𝑎^′_𝑘+^𝑛𝐾_𝑖=1𝑘=1𝑝_𝑖,𝑘log𝜙𝐱_𝑖∣𝝁^′_𝑘,Σ^′_𝑘, (4) where 𝑝_𝑖,𝑘=𝑎_𝑘𝜙𝐱_𝑖∣𝝁_𝑘,Σ_𝑘∑_𝑙𝑎_𝑙𝜙𝐱_𝑖∣𝝁_𝑙,Σ_𝑙𝜙𝐱_𝑖,∣𝝁_𝑘,Σ_𝑘=1(2𝜋)^𝑑Σ_𝑘−𝐱exp^𝑖−𝝁_𝑘^𝑡Σ_𝑘⁻¹𝐱_𝑖−𝝁_𝑘2. (5)(6) Hence, the EM iteration for a GMM is defined by the following. 𝐸-step: use (5).  𝑀-step: use the formulas 𝑎_𝑗=1𝑛_𝑖𝑝_𝑖𝑗,𝝁_𝑗=∑_𝑖𝑝_𝑖𝑗𝐱_𝑖∑_𝑖𝑝_𝑖𝑗,Σ_𝑗=∑_𝑖𝑝_𝑖𝑗𝐱_𝑖−𝝁_𝑗𝐱_𝑖−𝝁_𝑗^𝑡∑_𝑖𝑝_𝑖𝑗 (7) which are derived from the Lagrange multipliers, 𝜕𝑄/𝜕𝝁_𝑗=0 and 𝜕𝑄/𝜕Σ_𝑗⁻¹=0 (for details, see the appendices).

The above steps (i.e., 𝐸-step and 𝑀-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(𝑥_𝑖)−3ℎ, 𝛽>max(𝑥_𝑖)+3ℎ, 𝑀=2^𝑟; let ℎ be the bandwidth for some integer 𝑟, 𝛿=(𝛽−𝛼)/𝑀; and let 𝑡_𝑘=𝛼+𝑘𝛿 be the 𝑘th grid point where 𝑘=0,1,…,𝑀−1. The density estimate at grid point 𝑡_𝑘 is represented by the following equation: 𝑓𝑡_𝑘≈^𝑀/2_{𝑙=−𝑀/2}1𝑀^𝑀−1_𝑘=0𝜉_𝑘exp2𝜋𝑘𝑙𝑀𝑖×exp−2𝜋𝑘𝑙𝑀𝑖−12ℎ²2𝜋𝑙𝛽−𝛼², (8) where 𝑖²=−1. For 𝑥∈[𝑡_𝑘,𝑡_𝑘+1], the density estimate 𝑓(𝑥) is defined by 𝑓(𝑥)=(1/𝑛ℎ)∑^𝑛_𝑖=1𝐾((𝑥−𝑥_𝑖)/ℎ) where 𝐾(𝑡)=(1/√2𝜋)exp(−(1/2)𝑡²). To compute 𝑓(𝑥) at a grid of points, a method which makes use of the Fourier transform is employed. Let 𝑓(𝑠) be the Fourier transform of the kernel density estimate 𝑓(𝑥). It can be shown that 𝑓(𝑠)=(2𝜋)^1/2𝐾(ℎ𝑠)𝑢(𝑠)=exp(−(1/2)ℎ²𝑠²)𝑢(𝑠), where 𝐾(𝑠) is the Fourier transform of the Gaussian kernel and 𝑢(𝑠)=(2𝜋)^−1/2𝑛⁻¹∑^𝑛_𝑗=1exp(𝑖𝑠𝑥_𝑗) is the Fourier transform of the data. Thus, 𝑓(𝑥)=(2𝜋)^−1/2∫𝑒^{−𝑖𝑠𝑥}(2𝜋)^1/2𝐾(ℎ𝑠)𝑢(𝑠)𝑑𝑠 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 𝑓(𝑥) by convolving the data with the kernel.

Step A. Discretize the data to find the weight sequence {𝜉_𝑘} with 𝑀=2⁸. If 𝑥∈[𝑡_𝑘,𝑡_𝑘+1], it is split into a weight (1/𝑛𝛿²)(𝑡_𝑘+1−𝑥) at 𝑡_𝑘 and a weight (1/𝑛𝛿²)(𝑥−𝑡_𝑘) at 𝑡_𝑘+1; these weights are accumulated over all the data points 𝑥_𝑖 to give a sequence of (𝜉_𝑘) weights summing up to 1/𝛿.

Step B. Find the sequence {𝑌_𝑙} defined by 𝑌_𝑙=𝑀⁻¹∑^𝑀−1_𝑘=0𝜉_𝑘exp((2𝜋𝑘𝑙/𝑀)𝑖), where −(𝑀/2)≤𝑙≤𝑀/2. It can be shown that when 𝛼=0, 𝑌_𝑙≈(2𝜋)^1/2(𝛽−𝛼)⁻¹𝑢(𝑠_𝑙), where 𝑠_𝑙=2𝜋𝑙/(𝛽−𝛼).

Step C. Find the sequence {𝜁^∗_𝑙}, where 𝜁^∗_𝑙=exp(−(1/2)ℎ²𝑠²_𝑙)𝑌_𝑙. Here, ℎ=0.9𝐴𝑛^−1/5, where 𝐴=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 𝜁_𝑘 be the inverse discrete Fourier transform of 𝜁^∗_𝑙, that is, 𝜁_𝑘=∑^𝑀/2_{𝑙=−𝑀/2}𝜁^∗_𝑙exp(−(2𝜋𝑘𝑙/𝑀)𝑖).

It can be shown that when 𝛼=0, 𝑓(𝑡_𝑘)≈𝜁_𝑘. We then identify 𝑥_𝑖 where its density estimate, denoted by 𝑓(𝑥_𝑖), is greater than those of its nearest neighbors 𝑥_𝑖+1 and 𝑥_𝑖−1. In other words, 𝑓(𝑥_𝑖)>𝑓(𝑥_𝑖+1) and 𝑓(𝑥_𝑖)>𝑓(𝑥_𝑖−1); refer to Figure 1 where the vertical line that touches 𝑡_𝑘 and 𝑓(𝑡_𝑘) 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 of 𝑓(𝑡𝑘) against 𝑡𝑘.

3.3. Improved EM Algorithm for GMM

A number of authors highlighted the importance of identifying the right number, say 𝑘, of components in a GMM and subsequently choosing good initial values for the model parameters 𝜇_𝑖 and 𝜎²_𝑖, 𝑖=1,2,…𝑘, 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 𝑘 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 Section 3.2 is used to determine the number, say 𝐾₀, of components and also the corresponding means 𝜇_𝑖 of each component, where 𝑖=1,2,…,𝐾₀. The initial estimates of the standard deviations 𝜎_𝑖𝑖 are set to unity while the prior weights 𝑎_𝑖 are set to be 1/𝐾₀.

Step 2. The EM algorithm for a GMM as described in Section 3.1 is executed to give the final estimates of parameters 𝜇_𝑖, 𝜎_𝑖𝑖, and 𝑎_𝑖, 𝑖=1,2,…,𝐾₀. The log-likelihood function and Akaike information criteria (AIC) are calculated using the said parameters.

Step 3. Step 2 is repeated for other possible number 𝐾 of components with 𝜇_𝑖=0, 𝜎_𝑖𝑖=1 for the other 𝐾−𝐾₀  components, and 𝑎_𝑖=1/𝐾.

Step 4. The log-likelihood function and AIC values for 𝐾=1,2,…,10 are plotted. The final number of components 𝐾_𝑓 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: 𝑧_𝑗=𝜇+−2𝜎²log𝑢_𝑗^1/2cos2𝜋𝑢_𝑗+1,𝑧_𝑗+1=𝜇+−2𝜎²log𝑢_𝑗^1/2sin2𝜋𝑢_𝑗+1, (9) where 𝑢_𝑗,𝑢_𝑗+1∼𝑈(0,1). For the case of two components, we start by generating a random number 𝑢₁∼𝑈(0,1). If 0<𝑢₁<𝑎₁, we generate two random numbers 𝑢₂∼𝑈(0,1) and 𝑢₃∼𝑈(0,1) and calculate 𝑧₂+𝑧₃ using the first and second equations of (9) with 𝜇^∗=𝜇₁/2 and 𝜎^∗=𝜎₁/√2. Otherwise, we use 𝜇^∗=𝜇₂/2, and 𝜎^∗=𝜎₂/√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 Step 1 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 𝑎_𝑖 are distinctly different from the chosen values of 𝑎_𝑖 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 of a’s, and 𝜇’s, 𝜎’s.

(a) Sample 1

(b) Sample 2

(c) Sample 3

Figure 2 

Plots of values of log-likelihood function.

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 𝜽=(𝜃₁,…,𝜃₆)=(𝑎₁,𝑎₂,𝜇₁,𝜇₂,𝜎₁₁,𝜎₂₂) are to be estimated. Data is simulated using the simulation scheme described in Section 4.1.

After performing Steps 1 and 2, we find 𝐷_𝑖=𝜃_𝑖−𝜃_𝑖, where 𝜃_𝑖 is the true value of the 𝑖th parameter and 𝜃_𝑖 is the EM estimate of the parameter, 𝑖=1,2,…,𝑛. The sample mean and standard deviation of 𝐷_𝑖 are computed using the formulas 𝐷=(1/𝑛)∑^𝑛_𝑖=1𝐷_𝑖 and 𝑆_𝐷=1/(𝑛−1)∑^𝑛_𝑖=1(𝐷_𝑖−𝐷)². The estimates are considered good if 𝐷 is close to zero, indicating small biases observed in the simulation results, and 𝑆_𝐷 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 𝜇₁=0, 𝜇₂=3.0, and √𝜎₁₁=√𝜎₂₂=0.316 are fixed but the true values of 𝑎₁ and 𝑎₂ 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 𝜇₁=0, 𝜇₂=1.0, √𝜎₂₂=0.707, and √𝜎₁₁=0.447 are fixed but 𝑎₁ and 𝑎₂ 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 𝜇₁=0, 𝜇₂=0.25, √𝜎₁₁=0.577, and √𝜎₂₂=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 𝜇₁=0, 𝜇₂=3.0, and 𝜎²₁=𝜎²₂=0.316.

Table 3 

Simulation results for the case 𝜇₁=0, 𝜇₂=1.0, 𝜎²₂=0.707, and 𝜎²₁=0.447.

Table 4 

Simulation results for the case 𝜇₁=0, 𝜇₂=0.25, 𝜎²₁=0.577, and 𝜎²₂=1.414.

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 number 𝐾₀ 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 𝐾 of components, by setting 𝜇_𝑖=0, 𝜎_𝑖𝑖=1 for the other 𝐾−𝐾₀ components, and 𝑎_𝑖=1/𝐾. The final number of components 𝐾_𝑓 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. Table 5 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.

We consider real data consisting of the converted duration of each call made by Customer A (referring to the fifth column of Table 5), 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 𝑎₁=0.64, 𝑎₂=0.36, 𝜇₁=−0.66, 𝜇₂=1.17, 𝜎₁₁=0.07, and 𝜎₂₂=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.

(a) Log-likelihood

(b) AIC

Figure 4 

Plots of log-likelihood and AIC values.

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.

Appendices

A. Derivation of the First, Second, and Third Equations of (7)

(A.1) Using Lagrange multipliers defined by max/min 𝐹(𝑥,𝑦,𝑧) subject to Φ(𝑥,𝑦,𝑧)=0, 𝐺(𝑥,𝑦,𝑧)=𝐹(𝑥,𝑦,𝑧)+𝜆Φ(𝑥,𝑦,𝑧), 𝜕𝐺/𝜕𝑥=0,𝜕𝐺/𝜕𝑦=0,𝜕𝐺/𝜕𝑧=0 [22] on max∑_𝑖∑_𝑗𝑝_𝑖𝑗log(𝑎_𝑗) subject to ∑_𝑗𝑎_𝑗=1 (or (∑_𝑗𝑎_𝑗−1)=0), we get the first equation of (7).(A.2) From (𝜕/𝜕𝝁_𝑗)((1/2)∑_𝑖∑_𝑗𝑝_𝑖𝑗(𝐱_𝑖−𝝁_𝑗)^𝑡Σ_𝑗⁻¹(𝐱_𝑖−𝝁_𝑗))=0, we get the second equation of (7) by using the following matrix properties: 𝜕𝐱^𝑡𝐀𝐲/𝜕𝐱=𝐀𝐲,𝜕𝐚^𝑡𝐱/𝜕𝐱=𝐚.(A.3) The first and second expressions of 𝜕𝜕Σ_𝑗⁻¹12_𝑖𝑗𝑝_𝑖𝑗𝐱_𝑖−𝝁_𝑗^𝑡−1_𝑗𝐱_𝑖−𝝁_𝑗+𝜕𝜕Σ_𝑗⁻¹12_𝑖𝑗𝑝_𝑖𝑗log⁻¹_𝑗=0 (A.1) use the following matrix properties: 𝜕tr(𝐱𝐲)/𝜕𝐱=𝐲+𝐲^𝑡−Diag(𝐲) and ∑𝐱^𝑡_𝑖𝐀𝐱_𝑖=tr(𝐀∑𝐱_𝑖𝐱^𝑡_𝑖) to get the third equation of (7) [8].

B. The Value of Intersections

For the case when 𝜇₁≠𝜇₂ and 𝜎₁≠𝜎₂, 𝑓₁(𝑥)=(1/𝜎₁√2𝜋)𝑒^{−(1/2)((𝑥−𝜇1)/𝜎1)2} and 𝑓₂(𝑦)=1/𝜎₂√2𝜋𝑒^{−(1/2)((𝑦−𝜇2)/𝜎2)2} are obtained from 𝑥₁₁=(−𝑏+𝑏²−4𝑎𝑐)/2𝑎, and 𝑥₁₂=(−𝑏−𝑏²−4𝑎𝑐)/2𝑎 where 𝑎=(𝜎²₂−𝜎²₁), 𝑏=2(𝜎²₁𝜇₂−𝜎²₂𝜇₁) and 𝑐=(𝜎²₂𝜇²₁−𝜎²₁𝜇²₂)−2𝜎²₁𝜎²₂log(𝜎₂/𝜎₁).

Firstly, using the above formula as well as 𝑃((𝑥−𝜇)/𝜎)=∫^{(𝑥−𝜇)/𝜎}_−∞(1/√2𝜋)𝑒^{−(1/2)𝑡2}𝑑𝑡, we find the area between 𝑥₁₁ and 𝑥₁₂ (and convert it into percentage) for each component; refer to Figure 5(a). Secondly, we find the minimum between the areas of the two components. This value represents the percentage of overlapping between two components (which is an approximation). For the case when 𝜇₁≠𝜇₂, 𝜎₁=𝜎₂, 𝑓₁(𝑥)=(1/𝜎₁√2𝜋)𝑒^{−(1/2)((𝑥−𝜇1)/𝜎1)2}, and 𝑓₂(𝑦)=(1/𝜎₂√2𝜋)𝑒^{−(1/2)((𝑦−𝜇2)/𝜎2)2}, let 𝑑=(𝜇₁+2𝜎₁)−(𝜇₂−2𝜎₂). The value of the intersection, say 𝑥₁, is obtained from the following formula (which is an approximation): 𝑥₁=𝜇0,𝑑<0,₁+2𝜎₁𝜇,𝑑=0,₁+2𝜎₁−𝑑2,𝑑>0. (B.1) Taking similar steps, the area for the component on the left hand side of Figure 5(b) is obtained from 1−𝑃((𝑥₁−𝜇₁)/𝜎₁)=1−∫^{(𝑥1−𝜇1)/𝜎1}_−∞(1/√2𝜋)𝑒^{−(1/2)𝑡2}𝑑𝑡 and that of the component on the right hand side of Figure 5(b) from 𝑃((𝑥₁−𝜇₂)/𝜎₂)=∫^{(𝑥1−𝜇2)/𝜎2}_−∞(1/√2𝜋)𝑒^{−(1/2)𝑡2}𝑑𝑡. We convert them into percentages before adding them up to represent the percentage of overlapping between two components (which is an approximation).

(a)

(b)

Figure 5 

∫𝑓1(𝑥)𝑑𝑥 and ∫𝑓2(𝑥)𝑑𝑥 are used to find the shaded areas.

Conflict of Interests

Mohd Izhan Mohd Yusoff, the main author of the paper worked in Telekom Research & Development Sdn Bhd. (a subsidiary of Telekom Malaysia Berhad) for almost 20 years, the same organization that is sponsoring his Ph.D. study at the Institute of Mathematical Sciences, University of Malaya. The title of his study was proposed by two supervisors, who are the coauthors of the paper, namely, Ibrahim Mohamed and Mohd Rizam Abu Bakar. Due to the affiliation of the first author to Telekom Malaysia Berhad, the authors are able to obtain the relevant data sets to illustrate the development of the theory considered in this paper. In other words, Telekom Malaysia Berhad, as mentioned in the paper, is merely supplying the customer’s call detail record and does not provide any financial contribution towards the project. The data analysis and the paper preparation have been carried out independently. In conclusion, there is no potential conflict of interests in the study.

Acknowledgments

The authors would like to express their thanks to Telekom Malaysia Berhad for their cooperation in this study. This research is financially supported by the Fundamental Research Grant Scheme (no. FP012-2013A).