Modified MM Algorithm and Bayesian Expectation Maximization-Based Robust Localization Under NLOS Contaminated Environments

Robust localization methods that employ distance measurements to predict the position of an emitter are proposed in this paper. The occurrence of outliers due to the non-line-of sight (NLOS) propagation of signals can drastically degrade the localization performance in crowded urban areas and indoor situations. Hence, robust positioning methods are considered to mitigate the effects of outliers. Specifically, localization methods based on robust statistics are considered. Modified multi-stage ML-type method (MM) based weighted least squares (WLS), maximum a posteriori (MAP) expectation maximization (EM) WLS and variational Bayes (VB) EM WLS algorithms are developed under various outlier-contaminated environments. Simulation results show that the position estimation accuracy of the proposed modified MM WLS method, which uses the novel weight, is higher than that of the other methods under most outlier-contaminated conditions. Furthermore, the MAP-EM WLS and VB-EM WLS methods are the most accurate among algorithms that do not require statistical testing. Additionally, the mean square error (MSE) and asymptotic unbiasedness of the proposed algorithms are analyzed.


I. INTRODUCTION
Emitter positioning involves predicting the coordinates of a point emitter using the observations from each receiver. It is of considerable interest in diverse areas of study, such as telecommunication, signal processing and remote sensing. Emitter positioning algorithms are classified into direct and indirect localization methods. The direct localization method estimates the location information directly from a signal model based on the grid search method. In contrast, the indirect localization method estimates the source position using information such as the time difference of arrival (TDOA), time of arrival (TOA), received signal strength indicator (RSSI) and angle of arrival (AOA). The direct localization method has better accuracy in the low signal-to-noise ratio (SNR) condition, but incurs a higher computational The associate editor coordinating the review of this manuscript and approving it for publication was Chao Tan . burden compared with that of the indirect algorithm. The indirect positioning algorithm is examined in this study. The localization methods based on TOA, TDOA and AOA require expensive and energy-consuming hardware. The TDOA and TOA-based positioning methods also require synchronization between sensors or between sensors and the emitter. However, the TOA-based positioning method shows the best location performance among the aforementioned algorithms. The localization systems using RSSI can be implemented at a reasonable price, but their localization accuracy using the RSSI is lower than that of the techniques that use the TOA, TDOA and AOA. Moreover, this technique suffers from problems such as model-mismatch, multipath effects and background interferences, which render the distance observation inaccurate. Therefore, the TOA is utilized as the measurement in this paper. Previous studies have extensively investigated the localization problems in the line-of-sight (LOS) context [1]- [3]. However, certain problems remain to be solved. VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ A critical theme among position prediction problems is the localization in mixed LOS/non-line-of-sight (NLOS) environments. For example, the LOS path between a source and sensors may be blocked in urban or indoor situations. Robust statistics are required because the conventional localization methods based on a sample mean are quite vulnerable to outliers in mixed LOS/NLOS environments. The outlier is an observation that is distant from most gathered observations. In general, studies on position prediction for the mixed LOS/NLOS problems are classified into three main areas: 1) robust statistics, 2) LOS and LOS/NLOS mixture sensor identification and 3) mathematical programming. Some examples of localization methods based on robust statistics are the M and least median squares (LMedS) algorithms [4]- [9]. The M estimator exploits a robust cost function to minimize the effects of outliers [4], [5] and the LMedS estimator yields the smallest median of squared residuals for the candidate estimates [6], [7]. The identification and discard method identifies the outlier-contaminated sensors and discards them prior to localization [10], [11]. Quadratic programming, linear programming and interior point optimization are mathematical programming methods [12]- [14], wherein, some constraints associated with NLOS measurements are employed. Some papers deal with the TOA localization under asynchronous receivers formulation in wireless networks using the expectation maximization (EM) and factor graph. In this study, we focus on robust-statistics-based localization.
The MM estimator has been a popular method because it attenuates the adverse influences of NLOS observation [4], [5]. The existing MM estimation-based localization algorithm uses a transformed range observation to obtain the position estimate and the weight to reflect the different accuracy of each sensor [15]. The transformed observation and transformed distance estimate based on the MM algorithm are assumed statistically independent in the existing MM-WLS algorithm, but this is not the case. Therefore, novel weights are derived using the fact that the transformed range measurements and transformed distance estimate based on MM estimate are correlated. From now on, the modified MM weighted least squares (WLS) estimation is abbreviated as the MM WLS-M algorithm. Moreover, the MM WLS-M method requires statistical testing to discriminate the inliers and outliers. M denotes the modified. This necessitates the determination of an optimal threshold for discriminating inliers and outliers. However, this is a non-trivial task under rapidly varying environments or channels, as the threshold needs to be determined considering the changing conditions. Thus, maximum a posteriori (MAP) EM and variational Bayesian (VB) EM WLS localization without statistical testing are proposed to avoid this limitation. Statistical testing is not required as the parameters, such as mean and variance, are estimated from the outlier-contaminated observations. In the remainder of this paper, the MAP EM and VB EM combined with the WLS method are denoted as the MAP-EM WLS and VB-EM WLS algorithms, respectively. We utilize the EM algorithm because it has been widely used in the estimation of parameters of mixture model. The MAP EM method is the Bayesian version of EM algorithm in which the posterior is maximized and can outperform the standard EM algorithm in the low SNR or small sample size condition. The VB learning method has the advantages compared to the MAP algorithm. The posterior to be maximized is difficult to be known in certain case. In this situation, the posterior can be estimated using the mean-field theorem of the VB EM algorithm. Also, the VB learning is less prone to the overfitting than MAP estimation [16]. The essential of the VB EM algorithm is that the posterior can be represented as the multiplication of the variational distributions. Therefore, the performance of the VB EM algorithm may be degraded if the variational distribution is far from the true posterior. Although the Monte-Carlo-based algorithms such as Gibbs sampling and Metropolis-Hastings algorithms can be used for estimating parameters, the proposed MAP-EM and VB-EM algorithms are computationally efficient than the Monte-Carlo methods. Also, the convergence of the algorithm can be confirmed more easily than the Monte-Carlo methods. The novelty and main contributions of this study are summarized as follows: • We develop a robust MM WLS-M localization method in which the weight calculation is different from that of existing MM WLS localization algorithm. Namely, the statistical dependence between transformed range measurements and the MM estimate are considered in the proposed method for derivation of the weight. In general, the weight is determined as the variance of response variables. As the weight is more accurate, the accuracy of the estimation algorithm is better.
• The existing statistical testing method requires an optimal threshold for discriminating the inliers and outliers. However, this is not an easy task under rapidly changing environments or channels, as the threshold should be predetermined according to the changing conditions. To address these drawbacks, we propose a localization technique that does not require hypothesis testing and an empirical threshold.
• To the best of the authors' knowledge, the MM estimation-based localization using the novel weight proposed in this paper has not been investigated yet. Moreover, the combination of Bayesian EM and WLS methods has not been reported yet. The remainder of this paper is organized as follows. Section II describes the mixed LOS/NLOS position determination problem covered in this paper. Section III describes the existing estimation algorithms. Section IV describes the proposed robust MM WLS-M, MAP-EM WLS and VB-EM WLS positioning algorithms. Section V analyses the mean square error (MSE) and asymptotic unbiasedness of the proposed localization algorithms. Section VI demonstrates the superiority of the proposed techniques over the conventional algorithms via simulation results. Finally, Section VII presents our conclusions and future work.

II. PROBLEM FORMULATION
The localization method based on distance measurements aims to determine the location of an emitter accurately, where the error risk function, e.g., the MSE, should be minimized. In the problem of mixed LOS/NLOS emitter localization, the observation equation is expressed as where n i,j is a Gaussian mixture random variable defined as . . , M, j = 1, 2, . . . , P with M and P denoting the number of receivers and observations in each receiver, respectively. The inliers follow N (0, σ 2 i,1 ) and the outliers follow N (µ i,2 , σ 2 i,2 ) (N (µ i,l , σ 2 i,l ) denotes a Gaussian probability density function (PDF) with mean µ i,l and variance σ 2 i,l (l=1, 2)). The inliers have a probability of 1−ν i,2 and the outliers have a probability of ν i,2 . The Gaussian mixture noise distribution with two-mode has been utilized extensively for modeling mixed LOS/NLOS situations [17]- [21]. More specifically, the authors of [17], [19] obtained the experimental result that the mixed LOS/NLOS noise distribution is consistent with two-mode mixed Gaussian distribution. It is assumed that the statistics such as the mean and variance of the LOS and NLOS noise distributions cannot be obtained. Here, ν i,2 ∈ [0, 1] represents a contamination ratio, which is commonly lower than 0.1 [8]. Some papers also consider the exponential distribution for modeling the NLOS bias. However, the applicability of that model in the proposed algorithm is out of scope because it requires the thorough investigation. Moreover, [x y] T denotes the target location to be estimated and [x i y i ] T represents the known coordinates of the ith receiver. Additionally, r i,j is the range observation from the emitter to the ith receiver at the jth time instance and d i is the true distance from the source to the ith receiver. Squaring (1) and rearranging the expressions yield the following equation: where We can obtain the following equation by expressing (2) in a matrix form where . .
The objective is to determine the parameter x by fusing the transformed distance measurements optimally, i.e., [b 1 , · · · , b P ]. The first-step solution is obtained using the aforementioned equations. Furthermore, the second-step procedure is employed to enhance the accuracy of the first-step solution. In this paper, a vector is represented by a lowercase boldface letter and a matrix by an uppercase boldface letter.
The operator [·] T denotes a transpose.

III. REVIEW OF THE CONVENTIONAL APPROACHES UNDER GAUSSIAN MIXTURE MODELS
The MM estimator has been employed in the presence of outliers. The rationale of this algorithm is to mitigate the adverse effects of outliers by utilizing a weight function, that is inversely proportional to the squared residual. The mixture EM method has been used to estimate the parameters (mixing coefficient, mean and variance) of the mixture model. Moreover, these parameters have been used for robust positioning in the presence of outliers [17], [21]. The lower bound of the marginal likelihood function is maximized when the marginal likelihood function is difficult to be directly maximized. The mixture EM method can be extended to the MAP-EM and VB-EM algorithms using a prior distribution. The principle of the MAP-EM and VB-EM algorithms is to estimate the responsibility of the observation via conditional expectation.
A. MM ESTIMATOR [8], [9] The MM estimation involves the calculation of the scale parameter based on the S estimation, followed by location estimation using the M algorithm. The details of MM estimator can be found in [8], [9], [15].

B. MAP-EM ESTIMATOR [22]
The EM method is a recursive method to find the local maximum of a complete likelihood function, where the model depends on unobserved latent variables [23], [24]. Also, EM method can be employed for finding the parameters maximizing the log-posterior function and it is termed as MAP-EM technique in this study. The MAP-EM method is divided into an expectation (E) step and a maximization (M) step. The expectation of the log-posterior evaluated using the current estimate is determined in the E-step and the parameters that maximize the log-posterior function are determined in the M-step. The objective of the MAP-EM algorithm is to find the MAP estimate of the posterior function using the following two-step. E-step: where ln(·) denotes the natural logarithm. In the M-step, the expected log-posterior at the current estimate is maximized as follows: M-step: C. VB-EM ESTIMATOR [16], [22], [25], [26] The VB-EM algorithm estimates the random parameter θ and the latent variable z of a model, given the observation VOLUME 9, 2021 r. The Bayes' rule can be applied for inferring the posterior as follows: The evidence p(r) is determined to be The main limitation of the Bayesian estimation is that the above integral is intractable, except for simple models. The variational algorithm is designed to detour infeasible integrals, approximating true posterior distribution p(θ , z|r) using the variational posterior distribution q(θ , z).
The Kullback-Leibler (KL) divergence is used for measuring the discrepancy between the true and variational posterior distributions, and is defined as where E q {·} denotes the expectation over distribution q. The marginal distribution ln p(r) can be rewritten from (10) as follows: becomes the lower bound of the marginal distribution because D KL (q||p) ≥ 0 and B(q) is maximized for estimating θ , which is equivalent to minimizing the KL divergence. The variational posterior q(θ , z) can be factorized as q(θ)q(z) in the variational Bayesian inference. The variational Bayesian EM method alternates between E-step and M-step, i.e., the optimal q(z) is determined by assuming that q(θ ) is fixed and the optimal q(θ ) is determined with q(z) fixed. The solutions are represented as follows:

IV. PROPOSED MM WLS-M, MAP-EM, VB-EM WLS LOCALIZATION METHODS
In this section, the MM WLS-M method is described in detail. Moreover, the MAP-EM and VB-EM algorithms are combined with the WLS method.

A. MM WLS-M ALGORITHM
A detailed explanation of the modified MM estimation-based WLS localization algorithm is presented in this subsection.
T is the transformed distance set and merged employing the MM estimation as follows: where is the estimate of true transformed distance that combines b i,1:P andσ i is the S-estimate for the i th sensor. The initial point ofb MM i is selected as the sample median of b i,1:P [8], [9]. Mostly, it can not be known whether the measurement is an inlier or outlier; thus, the inliers and outliers should be discerned utilizing the statistical testing given in the below of (14). The estimate of true transformed distancê b MM i is then found recursively since both sides of (14) contain unknownb MM i . The variance of the statisticb MM i is then calculated based on algebraic methods as follows: where which are determined as inliers in the i th sensor based on the statistical testing appeared in the below of (14), and Q i is the total number of measurements decided to be inliers. In The output of this algorithm is given by: The y i,1:P are averaged using the sample mean to increase the accuracy, i.e.,b i = P j=1 y i,j P . The variance of the statisticb i was determined using algebraic methods in existing method as follows [15]: where S i = q i s (y i,q i − me i2 ) 2 , R i is the number of samples determined as outliers in the i th sensor, P = (Q i + R i ) is the total number of samples in the i th sensor, me i2 = (15). The statistical independence between b i,q i andb MM i is assumed in existing MM WLS algorithm. However, this is not true becauseb MM i includes b i,q i . Therefore, the variance calculation should be modified using the property that b i,q i andb MM i are correlated as follows: where q i s and r i s are the LOS sensor index set and LOS/NLOS sensor index set, respectively. In (18), the third line is derived using the property that the weights that belong to the outlier set are zeros. Then, the transformed measurement equation is determined as follows: where the superscript f is the acronym of the first-step estimate and ε is the additive noise ofb f . The first-step WLS estimate based on the M estimation is determined as follows: where Furthermore, the estimation performance of the first-step estimate is enhanced utilizing the second-step approach [2], [3]. Because the two-step localization method is well known, its explanation is omitted. For more details of the two-step localization method, it can be referred to [2], [3].

B. MAP-EM WLS ALGORITHM
An EM estimate of the parameters of mixture distribution can be obtained by maximizing the complete data log-likelihood, i.e., the lower bound of the marginal likelihood function. The objective is the maximization of the marginal likelihood function with respect to the parameters of the mixture model, which is obtained as follows: where However, the optimization is difficult because the summation exists in the inner part of the logarithm. Instead of solving (21) directly, we maximize the complete log-likelihood function in which the latent variable, z j,l (j = 1, · · · , P, l = 1, 2), is appended to the incomplete range observation. z j,l indicates the mixture component responsible for the corresponding distance observation. Moreover, the EM method can be used for optimizing the complete log-posterior and the conditional expectation for the complete log-posterior of the i th sensor is expressed as where z i is a vector enclosing latent variables z j,l and superscript (k) denotes the k th iteration-step. We assume that the prior of the mixing coefficient, ν i , is modeled as the Dirichlet distribution with the parameter α i , where α i consists of α i,l = α i,0 + N i,l (i = 1, · · · , M , l = 1, 2), N i,l = P j=1 γ j,l , γ j,l = p(z j,l |r i,j , θ (k) i ) and α i,0 is the hyper-parameter. The prior mean (µ p ) is modeled as a Gaussian distribution with zero mean and the precision prior (λ p ) is modeled as a gamma distribution with the parameter (a 0 , b 0 ). Then, (22) can be further derived as follows: (2α i,l ) and (·) denotes the gamma function. Rearranging (23) with respect to ν i yields (24): where const represents the irrelevant terms of ν i = {ν i,1 , ν i,2 }. Differentiating (24) with respect to ν i,l and setting it to zero yield (25): where β represents the Lagrange multiplier and the constraint 2 l=1 ν i,l = 1 is used. Then, ν i,l is obtained by solving (25) VOLUME 9, 2021 as follows: can be obtained similarly as follows: In the following, the superscript (k) is omitted for brevity of notation. Subsequently, the WLS method is derived using the mean and covariance of r i = [r i,1 · · · r i,P ] T (i = 1, · · · , M ). The mean of r i is represented as follows: Also, the variance of r i is expressed as follows: Subsequently, the mean and variance of the transformed observation are obtained as follows: where It should be noticed that the mixture EM algorithm has not been combined with the WLS method in existing literatures. The MAP-EM WLS method does not require statistical testing because the mean and variance of the observations are calculated with the outliers included. This property is also applied to the VB-EM WLS method introduced in the following section. The transformed measurement equation for the MAP-EM WLS method is expressed as follows: where ε = [ε 1 , · · · , ε M ] T , b m = [µ b 1 , · · · , µ b M ] T , the superscript m is the abbreviation of the mean and ε is the noise component of b m . Then, the first-step WLS estimate is determined as follows: Furthermore, the accuracy of the first-step estimate can be improved using the second-step approach and it is identical to that of the MM WLS-M algorithm.

C. VB-EM WLS ALGORITHM
The VB-EM-based method has been used for approximating the posterior probability using unobserved latent variables given some observed data. The approximated distribution is restricted to belong to an exponential family to make the approximated distribution similar to the true posterior. The joint distribution of r i , θ i = [ν i , µ i , λ i ] T and z i is expressed as follows: The variational distribution is factorized with respect to the latent variables and the parameters so that First, the variational distribution of z i is derived using the mean field theory [25], [26] as follows: Absorbing any terms independent of z i into the constant yields: ρ j,l , β i,0 and ν i,0 are hyper-parameters and const is the term independent of z i . Then, q(z i ) can be expressed as follows: where u j,l = ρ j,l 2 k=1 ρ j,k . Second, the variational distribution of ν i is derived as follows: where const is the irrelevant term of ν i .
Finally, the variational posterior distribution q(µ i , λ i ) are obtained similarly using the above derivation as follows: . The hyper-parameters of the gamma distribution were set to zero without loss of generality in the derivation of q(µ i , λ i ). The variational posterior means of ν i,l , µ i,l , λ i,l and d i are obtained as follows: The range is estimated by maximizing the complete log-likelihood (22) as follows: Then, the mean and variance of the transformed observations are obtained using (30)− (33). The remaining procedures are identical to those of the MAP-EM WLS algorithm.

V. MSE AND ASYMPTOTIC UNBIASEDNESS ANALYSIS
This section analyzes the MSE and asymptotic unbiasedness of the proposed methods.

A. MSE PERFORMANCE ANALYSIS
The MSE of the estimator is the sum of the squared bias and error variance. The root MSE (RMSE) can be obtained as the square root of the MSE. The estimation error x e can be expressed as cov The MSEs of the MAP-EM and VB-EM WLS methods can be obtained similarly using the above procedure.

B. ASYMPTOTICAL UNBIASEDNESS ANALYSIS
The unbiasedness of the MM WLS-M estimator can be proven as follows. The expected value of b MM (i) can be expressed as follows: where n MM i denotes the estimation error ofr MM i . The second-order error term, (n MM i ) 2 , can be neglected under sufficiently low noise conditions. Furthermore, E[n MM i ] is zero because the MM estimator is an asymptotic unbiased estimator [8], [9]. Therefore, the bias of b MM (i) is approximately zero and the MM WLS-M estimator ( x e ) is asymptotic unbiased estimator. Moreover, lim  [27]. Furthermore, the unbiasedness of the VB-EM method can be verified. The Gaussian mixture parameters can be estimated using the EM algorithm and they are unbiased estimates when each Gaussian component is separated sufficiently [28]. Moreover, the second-order error terms can be neglected VOLUME 9, 2021 when the estimation is performed well. Then, E[µ b i ] can be expressed as follows: The following property is used in the derivation of (50): ν i,1 + ν i,2 = 1 and µ i,1 = 0. Usually, ν i,2 is close to zero. Therefore, the estimation bias is smaller when ν i,2 is closer to zero. However, the bias term cannot be neglected as ν i,2 is closer to one.

VI. NUMERICAL RESULTS
The localization performances of the proposed robust source localization algorithms are compared with those of MM WLS [15], robust EKF [29], robust extrapolated single propagation unscented Kalman filter (ESPUKF) [15], [30], [31] and semi-definite programming (SDP) [32] algorithms. The SDP method showed accurate performance under the heavily contaminated environment, thus it was included for comparing the localization performance. Also, the Kalman filter has been widely utilized for localization, the robust version of Kalman filter was compared with the proposed algorithm. We used Matlab R2019b to obtain the simulation results.

A. SIMULATION ENVIRONMENTS
The simulation settings are summarized as follows. The RMSE is defined as the estimation performance measure as follows:

B. RMSEs AS A FUNCTION OF STANDARD DEVIATION OF LOS AND NLOS NOISE
The positioning accuracy as a function the standard deviation of the outlier noise is illustrated in Fig. 2. The contamination ratio (ν i,2 ) was 20%, the sample length was 20 and the standard deviation of the inlier noise (σ i,1 ) was 0.3 m in Fig. 2(a). The bias of the outlier observation noise (µ i,2 ) was set to 20 m, such that the outlier-contaminated distance observation was consistently larger than zero. The standard deviation of the outlier noise (σ i,2 ) was set as a remarkably larger value than that of the inlier noise such that the observation noise is always positive. The RMSE of the proposed MM WLS-M method was lower than that of the other existing methods in Fig. 2(a). The MM WLS-M method outperformed the MM WLS algorithm because the correlation between the transformed range observations and the MM estimate was considered in the derivation of the weight. The RMSE of the VB-EM WLS method was slightly inferior to that of MM WLS method. However, the VB-EM algorithm did not require statistical testing to distinguish the inliers and outliers unlike the MM WLS method. Moreover, the Cramér-Rao lower bound (CRLB) was not displayed since its algebraic derivation for a nonlinear function and non-Gaussian distribution is non-trivial problem [17], [33]. As shown in Fig. 2(b), the contamination ratio was 40%, whereas the remaining environments were identical to those in Fig. 2(a). The RMSE of the proposed MM WLS-M method was the lowest among the localization methods. Further, the RMSEs of the VB-EM and MAP-EM WLS methods were lower that that of the MM WLS algorithm. Moreover, the proposed VB-EM and MAP-EM WLS methods outperformed the SDP and Kalmanfilter-based methods. In particular, the RMSE of the VB-EM WLS method was nearly identical to that of the MM WLS-M method in Fig. 2(c). The square root of error variances of the VB-EM WLS and MM WLS-M algorithms were similar with their RMSEs in Figs. 2(a)−(c). Fig. 3 assumes the identical situation with that appeared in Fig. 2, except that sensors 4−7 were outlier-contaminated sensors. The RMSE of the proposed MM WLS-M algorithm was lower than that of the existing and EM-based WLS algorithms. The RMSEs of the EM-based WLS algorithms were degraded. In Fig. 3(b), the contamination ratio  was 40%, whereas the remaining environments were same as those of Fig. 3(a). The RMSE of the proposed MM WLS-M algorithm was lower than that of the existing methods. As shown in Fig. 3(c), the contamination ratio was 80%, whereas the remaining environments were identical to those in Fig. 3(a). This situation modeled the severely outlier-contaminated environments. The performance of the SDP method was similar to that of the EM-based WLS methods under the severely outlier-contaminated conditions. Fig. 4(a) illustrates the RMSEs as a function of the standard deviations of inliers. Sensor 7 was postulated to be a LOS/NLOS mixture sensor, whereas the other sensors belong to LOS sensors. The contamination ratio was 40% and the accuracy of the proposed EM-based WLS and MM WLS-M methods were slightly higher than those of the other methods in Fig. 4(a). The RMSEs of the EM-based WLS methods were nearly identical to that of the MM WLS-M method. As shown in Fig. 4(b) and (c), the MM WLS-M algorithm outperformed the existing methods and the proposed EM-based WLS algorithms. The performances of all the robust methods became inferior as the standard deviation of the LOS noise was higher. In particular, the accuracy of the EM-based WLS methods was very sensitive to the number of LOS sensors as shown in Fig. 4(c). Hence, the EM-based WLS method is not recommended unless the number of LOS sensors is more than three.

C. RMSEs AS A FUNCTION OF SEVERAL PARAMETERS
The RMSEs with respect to bias are shown in Fig. 5. The estimation accuracies of the proposed methods were nearly unchanged with the variation in the bias and the proposed MM WLS-M method was superior to the other algorithms. The RMSEs of the proposed WLS-based algorithms were not influenced by the bias because the localization was dependent on the LOS sensors. Furthermore, comparing the localization methods that do not require statistical testing, the accuracy of the proposed MAP-EM and VB-EM WLS algorithms was higher than those of the other localization algorithms. Fig. 6 compares the accuracy as a function of the sample size. As the sample size increased, the RMSEs of all the methods decreased. The RMSE of the MM WLS-M method was lower than that of the other localization methods for all sample sizes. Moreover, comparing the localization methods that do not require statistical testing, the proposed EM-based WLS algorithms outperformed the existing SDP and Kalman-filter-based methods. Fig. 7 displays the RMSEs as a function of the number of sensors. The number of LOS/NLOS mixture sensors was set as two and the number of LOS sensors was increased one by one. As the number of sensors increased, the RMSEs of all the methods decreased. The RMSE of the MM WLS-M method was lower than that of the MM WLS method for the entire  number of sensors. However, the RMSEs of the MAP-EM and VB-EM WLS methods were severely degraded when the number of LOS sensors was less than three. This observation coincides with the results shown in Fig. 4(c), i.e., the EM-based WLS algorithms are not effective when the number of LOS sensors is less than three.
The RMSEs with respect to the contamination ratio (ν i,2 ) are shown in Fig. 8. The localization performances of all the methods were degraded as the contamination ratio increased. Specifically, the RMSEs of the localization algorithms were evidently degraded when the contamination ratio was higher than 0.5. The RMSE of the MM WLS-M method was lower than those of the other methods when the contamination ratio was less than 0.5, but it was similar to that of the MM WLS method when the contamination ratio exceeded 0.5. Fig. 9 provides the ECDFs with respect to estimation errors. The contamination ratio (ν i,2 ) was 90% and sensor 7 was an outlier-contaminated sensor, whereas the other sensors were LOS sensors. The other conditions were same as those of Fig. 2(a). This situation represents a mildly contaminated environment. The RMSEs of the MM WLS-M, MM WLS and EM-based WLS methods were nearly identical as shown in Fig. 9(a). Note that the VB-EM and MAP-EM WLS methods do not require statistical testing, where the optimal threshold must be determined. This advantage can be important under rapidly varying channel conditions or environments. The conditions in Fig. 9(b) were postulated to be identical to those of Fig. 9(a), except that sensors 4−7 were outlier-contaminated sensors. These simulation settings represent moderately outlier-contaminated situations. The RMSEs of the proposed methods were lower than those of the other existing techniques. All the sensors besides sensor 1 were outlier-contaminated sensors in Fig. 9(c). This situation represents a severely outlier-contaminated condition. In Fig. 9(c), the RMSEs of the EM-based WLS methods were lower than those of the other methods, although they were degraded significantly compared with those of the former cases. The RMSEs of all the robust methods were increased as the number of outlier-contaminated sensors was incremented. The EM-based WLS methods were significantly superior to the other methods under seriously outlier-contaminated situations. The simulation results illustrated that the MM WLS-M method was superior to the other algorithms in most environments and the EM-based WLS algorithms showed effective performances when the number of LOS sensors was more than three.

VII. CONCLUSION
Robust localization algorithms were developed using the MM WLS-M, MAP-EM and VB-EM WLS techniques. The proposed MM WLS-M method utilized a weighting matrix estimated using the property that transformed measurements are correlated with the MM estimate. The RMSE of the MM WLS-M method was lower than those of the other localization methods under all the conditions, excluding severely outlier-contaminated environments. Furthermore, the RMSEs of the MAP-EM and VB-EM WLS algorithms were the lowest among the localization algorithms that do not require statistical testing. However, the accuracies of the MAP-EM WLS and VB-EM WLS algorithms were severely degraded when the number of LOS sensors was less than three. Moreover, the MSE and asymptotic unbiasedness of the proposed methods were analyzed. The proposed localization methods are iterative methods, i.e., they may diverge under severely low-SNR or heavily outlier-contaminated noise conditions because the initial position estimate cannot be appropriately determined. Moreover, the computational complexity of the iterative method was higher than that of the closed-form technique. Therefore, algebraic robust localization algorithms shall be developed in the future work.