Total Msplit estimation

Msplit estimation is a method that enables the estimation of mutually competing versions of parameters in functional observation models. In the presented study, the classical functional models found in it are replaced by errors-in-variables (EIV) models. Similar to the weighted total least-squares (WTLS) method, the random components of these models were assigned covariance matrix models. Thus, the proposed method, named Total Msplit (TMsplit) estimation, corresponds to the basic rules of WTLS. TMsplit estimation objective function is constructed using the components of squared Msplit and WTLS estimation objective functions. The TMsplit estimation algorithm is based on the Gauss–Newton method that is applied using a linear approximation of EIV models. The basic properties of the method are presented using examples of the estimation of regression line parameters and the estimation of parameters in a two-dimensional affine transformation.


Introduction
proposed a method for estimating parameters in split functional models of geodetic observations. Such a split occurs when two functional models that differ from each other in terms of mutually competing versions of the same parameter can correspond to a single observation. For example, in a network deformation analysis carried out based on an aggregate set of observations obtained during two measurement epochs, any given observation from this set corresponds to one of the following models: a functional model of observations from the first measurement epoch or a functional model of observations obtained during the second measurement epoch. Another example concerns the sets containing outliers. In that case, any given observation from this set can be a "good" observation or a wrong observation with a functional model appropriate for it.
The proposed method called M split estimation assumes that if a particular observation already occurs, it brings two mutually competing pieces of f -information (Jones and Jones 2000) determined in relation to two versions of the same parameter (Wiśniewski 2009(Wiśniewski , 2010. M split estimators of B Zbigniew Wiśniewski zbyszekw@uwm.edu.pl 1 Department of Geodesy, Faculty of Geoengineering, University of Warmia and Mazury, 10-719 Olsztyn, Poland these versions are the quantities that minimise the aggregate information being the product of competing information. Similar assumptions are also adopted in the maximum likelihood method (ML-method) (e.g. Rao 1973;Wiśniewski 2017). However, this method does not allow the existence of several versions of the parameters in a functional model relating to the same observation. From this perspective, M split estimation can be regarded as a particular kind of development of the ML method. In the absence of competing versions of the parameter, M split estimators become ML-estimators. Since the study conducted by Huber (1964), a generalisation of the ML-method has been very popular, known as M-estimation, in which f -information is replaced by certain arbitrary functions. A similar substitution can also be observed in M split estimation. M split estimation based on L1 norm condition was developed to take advantage of this possibility (Wyszkowska and Duchnowski 2019).
The general theory of M split estimation was developed without detailed assumptions about probabilistic observation models, which enables the creation of M split estimation varieties corresponding to specific models of this nature. The most commonly accepted probabilistic model of geodetic observations is the normal distribution. The family of normal distributions corresponds to the basic variant of M split estimation called "squared M split estimation". This variant of M split estimation can be regarded as a particular type of expansion of the least-squares (LS) method (Wiśniewski 2009). In the absence of competing parameter versions, squared M split estimators become LS-estimators. Whenever M split estimation is mentioned further on, this will indicate this particular variant.
M split estimation was applied inter alia in the analysis of geodetic network deformation Wiśniewski 2012, 2014;Zienkiewicz 2014;Zienkiewicz et al. 2017;Wiśniewski and Zienkiewicz 2016). In these problems, M split estimation is particularly effective in identifying stable potential reference points (PRPs) (Nowel 2019). Janowski and Rapiński (2013) applied M split estimation in 3D modelling, primarily for the detection of surface structures (e.g. roof planes) of engineering structures. The modelling was carried out based on laser scanning data. A similar problem is also analysed in a study by Janicka et al. (2020), where M split estimation was proposed as a means of detecting and determining the displacements of adjacent planes. Laser scanning data also provided the basis for determining the terrain profiles using M split estimation (Błaszczak-Bąk et al. 2015;Wyszkowska et al. 2021).
M split estimation can provide an alternative to Mestimation which is robust to gross errors. The possibility of such applications was indicated in the following studies (Wiśniewski 2009(Wiśniewski , 2010Yang et al. 2010;Ge et al. 2013;Janicka and Rapinski 2013;Amiri-Simkooei et al. 2017). Wiśniewski and Zienkiewicz (2021a, b) demonstrated that with properly established, competitive functional models, the robustness of M split estimators to gross errors is their inherent feature. The robustness of these estimators in a wider context (e.g. to poorly chosen models) was analysed in detail by Wiśniewski (2019, 2020).
In both the traditional models and split functional models constructed on their basis (in M split estimation), it is assumed that only observations are affected by random errors. Currently, an errors-in-variables (EIV) model, in which design matrix elements are also affected by random errors, is applied in many geodetic problems. For example, this model was applied in geodetic datum transformation (Teunissen 1988;Davis 1999;Acar et al. 2006;Akyilmaz 2007;Schaffrin and Felus 2008;Mahboub 2012;Fang 2015;Aydin et al. 2018;Mercan et al. 2018)) as well as in remote sensing (Felus and Schaffrin 2005), in a function approximation (Wang and Zhao 2019), in linear regression (Schaffrin and Wieser 2008;Amiri-Simkooei and Jazaeri 2012;Zeng et al. 2018;Lv and Sui 2020) and in the least-squares collocation (Schaffrin 2020;Wiśniewski and Kamiński 2020). The effect of the random design matrix on the weighted LS estimate is presented in . This paper also proposed a bias-corrected weighted LS estimate for the EIV model. A developed EIV stochastic model and an estimation of the model's components (using, inter alia, the MINQUE method) are presented in .
The estimation of parameters in functional models extended to the EIV form is most commonly carried out using the total least-squares (TLS) method. The optimisation problem of this method as well as its solution based on the singular value decomposition (SVD) was presented by Golub and Loan (1980). TLS using the SVD procedure was developed and adapted to geodetic purposes as well (e.g. Felus 2004;Akyilmaz 2007;Schaffrin and Felus 2008). Another way to solve the TLS optimisation problem, based on a nonlinear Lagrange function, is proposed in Schaffrin et al. (2006).
In the practical applications of the TLS method, besides having effective algorithms at one's disposal, the possibility for taking into account random weights of EIV model components is also very important. The basic solutions in this regard were presented by van Huffel and Vandewalle (1991), who established the generalised total least-squares (GTLS) method. On the other hand, Schaffrin and Wieser (2008) proposed an expansion of the TLS method, in which weights were derived from the adopted covariance matrix models (stochastic models). In the method proposed in the cited study, called "the weighted total least-squares (WTLS) method", stochastic models can apply to both observation vectors and the vectors created from random errors affecting the design matrix.
WTLS is still developed and analysed. For example, Fang (2013) analysed necessary and sufficient conditions for WTLS optimality. Amiri-Simkooei (2017, 2018 presented the theory behind the constrained weighted total least-squares (CWTLS) method. The WTLS optimisation problem was also formulated and solved using a secondorder approximation function (Wang and Zhao 2019). Due to their nonlinear nature, the WTLS estimators are biased. Biascorrected versions of these estimators are presented in studies by (Xu et al. 2012;Tong et al. 2015). Moreover, an important problem in WTLS theory and practice is the assessment of the accuracy of the determined estimators. What might be helpful in this regard are the strategies for determining the covariance matrix of WTLS estimates (Amiri-Simkooei et al. 2016) and methods for estimating the variance components in EIV models .
In the optimisation problem of the WTLS method based on the Lagrange approach, the objective function is minimised with the conditions defined by the nonlinear EIV model. An iterative algorithm to solve this problem was proposed in (Schaffrin and Wieser 2008). Shen et al. (2011), based on Newton-Gauss algorithm of nonlinear LS adjustment (Pope 1974), proposed another iterative method for solving WTLS problems, which is easier in practical applications. In this model, a nonlinear EIV model is replaced with a linear approximation, which significantly facilitates the organisation of a corresponding computational algorithm.
The origin of TLS or WTLS estimation is the LS-method which is neutral for all observations. The WTLS estimators' lack of robustness to gross errors is, therefore, an inherent feature, which may restrict the scope of the practical application of these estimators. Therefore, a robust estimation of EIV model parameters is of interest to many authors. For example, Wang et al. (2016) proposed a robust total least-squares (RTLS) method in which the robustness of WTLS estimators was obtained by means of the application of weight functions adopted in robust M-estimation. Another proposal, based on the least trimmed squares (LTS) method, was presented by Lv and Sui (2020). In that method, the authors used the inherent robustness of estimators minimising the sum of the squared orthogonal errors. They called the LTS version adjusted to EIV models "total least trimmed squares" (TLTS).
The current study will apply the EIV model in a basic M split estimation variant. When referring to the WTLS method, models of covariance matrices of random components of this model will also be taken into account. According to the basic principles of M split estimation, it will be assumed that the parameters in EIV model will have two mutually competing versions (which consequently leads to the split of this model). The objective function of the proposed method, called the "Total M split (TM split ) estimation", will be created through the application of the Lagrange approach (Schaffrin and Wieser 2008) using the approach adopted in (Shen et al. 2011), i.e. the split EIV models will be replaced with their linear approximations.
The paper is organised as follows. As the proposed method is an expansion of M split estimation, which takes into account the basic assumptions used in the WTLS method, it appears necessary to review the theoretical foundations of both these methods. These foundations, set in the context relevant to this study, are provided in Sect. 2. The theory behind TM split estimation and its algorithm are provided in Sect. 3. In Sect. 4, examples of the method application will be provided. TM split estimation will be applied to estimate parameters in competing bias models (Sect. 4.1). The obtained results will be compared with classical M split estimators calculated in Wiśniewski (2010). In Sect. 4.2, the data provided in Neri et al. (1989) and also used, inter alia, in studies by Schaffrin and Wieser (2008), Shen et al. (2011) andMahboub (2012), will be used to determine TM split estimators of linear regression parameters. It will be assumed that the basic set of "good" observations is disturbed by "strange" observations for which a corresponding regression line also exists. Moreover, in this Chapter, the behaviour of TM split estimators will be checked in the event that the set contains one observation affected by gross error with different values. The determined TM split estimators will be compared with the WTLS estimators published in the cited studies. TM split estimators' robustness to gross errors is additionally analysed using an example of a two-dimensional affine transformation (Sect. 4.3). The data for this example are derived from (Lv and Sui 2020). The RTLS and TLTS estimators showed in the cited paper will be compared with TM split estimators. The paper concludes with a summary.

M split estimation
Let y = AX + v be a functional model of y = [y 1 , . . . , y n ] T observation vector, where A is the n × m coefficient matrix (rank(A) = m), X is the m-vector of unknown parameters to be estimated, and v is n-vector of random observation errors. In M split estimation, this model is split into two models: where X α and X β are mutually competing versions of the same vector of parameters X. The vectors v α ,v β are respective versions of the vector v, which result from the observation errors and the errors of the functional models. M split estimators of parameters X α and X β are quantities that minimise the following general objective function (Wiśniewski 2009).
where ρ α and ρ β are arbitrary functions. In the context of cross-weighting that is natural in M split estimation, function (2) can also be expressed in the following form: where w α (v iβ ) = ρ β (v iβ ) and w β (v iα ) = ρ α (v iα ) are now regarded as a special type of weight function. The specific character of weighting is that the contribution of function ρ α (v iα ) to the optimisation problem is enhanced (or weakened) by the weight function whose argument is quantity v iβ competing in relation to v iα (and vice versa). The weight functions are not like those in M-estimation, which are modified to make the estimator robust. Mutual "cross weighting" functions w α (v iβ ) and w β (v iα ) are applied to determine mutually competitive estimates related to the same observation set (Wiśniewski 2009). In the case of M split estimation, one supposes that the observation set might be a mixture of realisations of two different random variables that differ from each other in the parameters of the functional models. One of those variables might be regarded as a "strange" one and its realisations as outliers in a particular case. Then results of M split estimation are estimates of the parameters of the "good" variable (like in robust M-estimation) but also estimates of the parameters of the "strange" variable. This study will use the basic variant of M split estimation, in which ρ(v α ) = v 2 iα q −1 i and ρ(v β ) = v 2 iβ q −1 i . The quantities q i are diagonal elements of the Q y cofactor matrix occurring in the C y = σ 2 0 Q y covariance matrix model (σ 2 0 -unknown variance component). The adopted functions can be associated (although this is not necessary) with normal distributions as probabilistic observation models. Taking these functions into account, based on Eqs.
(2) and (3), the following will be recorded (Wiśniewski 2009;Zienkiewicz 2018aZienkiewicz , 2018b The solution to the optimisation problem ϕ(X α , X β ) → min includes such quantitiesX α andX β (M split estimators) for which the following is true: whereṽ α = y − AX α andṽ β = y − AX β are residual vectors. The above equations are solved by means of iteration. The iterative procedure can be organised in such a manner that in the steps l = 1, . . . , s, the following quantities are determined Zienkiewicz 2021a, 2021b): (the iterative procedure using gradients and Hessians of the function ϕ(X α , X β ) is presented in Wiśniewski 2009Wiśniewski , 2010. The iterative process defined by Eq. (7) is convergent and ends for such l = s that X α(s) = X α(s−1) and X β(s) = X β(s−1) . Then,X α = X α(s) andX β = X β(s) . In M split estimation, the stochastic model C v = C y = σ 2 0 Q y , similar to the functional model, is split. The split results in covariance matrices C v α = σ 2 0α Q y and C v β = σ 2 0β Q y , which are two versions of the covariance matrix C v (Wiśniewski and Zienkiewicz 2021b). The invariant and unbiased estimators of variance coefficients σ 2 0α and σ 2 0β are the following quantities (Wiśniewski and Zienkiewicz, 2021a, b): and where

Weighted TLS method
The total least-squares method is applied where the classical model y = AX + v is replaced by the EIV model of the following form: where E is an n × m random matrix corresponding to the matrix A being observed. The random vector corresponding to observation vector y in the EIV model was denoted as υ.
If we assume that v = y − AX is the error vector in the classical functional model, then υ = v + EX. It should be considered that EX = (X T ⊗ I n )e, where e = vec(E) is a vector formed from successive columns of matrix E (⊗-the Kronecker product, I n -an identity matrix of dimensions n× n). Moreover, in order to simplify further notation, additional designations (10) can then be expressed in the following form: In the WTLS method, in addition to the stochastic model of the observation vector C y = σ 2 0 Q y , a stochastic model of e vector is also adopted. In the simplest case, it can be assumed that such a model is the expression C e = σ 2 0 Q e , where Q e is the known cofactor matrix. However, there are examples in which not all columns of matrix A are affected by random disturbances (e.g. in the linear regression analysis). In that case, matrix Q e can be subject to appropriate decomposition Q e = Q 0 ⊗ Q x , where Q x denotes a nonnegative definite diagonal matrix of size n × n (Schaffrin and Wieser 2008;Shen et al. 2011).
In the Lagrange approach applied in Schaffrin and Wieser (2008) and Lv and Sui (2020), the WTLS method is based on the following objective function: where λ denotes an n × 1 vector of "Lagrange multipliers". Matrix Q + e is the Moore-Penrose inverse of Q e . For example, Rao 1973;Felus 2004).The iteration procedures that solve the optimisation problem ϕ(υ, e, λ, X) → min with the nonlinear conditions y − AX − υ + X T ⊗ e = 0 are, in general, complicated. The application of a similar objective function in M split estimation will generate numerical procedures and computational algorithms of even greater complexity. From the perspective of computational process optimisation, however, the approach adopted in Shen et al. (2011), also applied in Wang et al. (2016), is particularly interesting for the purposes of the study. The iterative method proposed in the study is based on the Newton-Gauss algorithm of nonlinear LS adjustment, proposed by Pope (1974). In this method, the EIV model (11) is replaced, in the j-the iteration, by a linear approximation of the following form: where A j = A−Ẽ j , X = X j +δX and X j ⊗ = X j ⊗I n .Ẽ j are residual matrices built on the basis of the residual vectorẽ j determined in the j-th iteration. Vector δX is a small quantity to be solved in the iteration. After taking into account the (13), the objective function (12) can be rewritten in the following form (Shen et al. 2011).
The minimum of this function is obtained through satisfying the following Euler-Lagrange necessary conditions (Shen et al. 2011) The solution to the equations contained in Eq. (15) are the following quantities: and (υ-residual vector corresponding to the observation vector y). Shen et al. (2011) also apply the iterative process on the assumption that E j δX is a negligible quantity. Then, The estimation of variance coefficient σ 2 0 , common to stochastic models C y = σ 2 0 Q y and C e = σ 2 0 Q e , is also of interest in WTLS. Schaffrin and Wieser (2008) proposed a biased estimator in the following form: A correction of the estimator (19) by means of introducing a bias term δb to it was presented by Shen et al. (2011). More complex EIV stochastic models are also being under consideration. For example, Xu and Liu (2014) introduced a model containing variance components and proposed a way to estimate these components.
By applying the approach used in Shen et al. (2011) and Wang et al. (2016), see Eq. (13), the above models in the j-th iteration will be replaced with the following linear approximations: After determining the quantities E j , X j α and X j β , the models that are valid in the j-th iteration, contained in Eq. (20), take the following forms: Total M split estimators of parameters X α and X β are quan-titiesX α andX β which minimise the following objective function: This function [lake functions Eqs. (3) and (4)] is not designed to make the method robust and also does not "predict" outliers among observations of elements of matrix A. In view of the following conditions resulting from Eq. (21) the original objective function (23) will be supplemented to the following form: where λ α and λ β are Lagrange multiplier vectors corresponding to conditions (24). It is established that the Euler-Lagrange necessary conditions have the following forms in the optimisation problem ϕ(υ α , υ β , e, X α , X β , λ α , λ β ) → min: The set of simultaneous substitutions: υ α =υ α , υ β =υ β , e =ẽ, δX α = δX α , δX β = δX β , λ α =λ α , λ β =λ β , introduced to simplify the notation, was denoted as . Based on Eqs. (26a)-(26c), the following residual vectors are determined: By substituting the quantities obtained above to Eqs. (26f) and (26g), a system of normal equations relating to vectorŝ λ α andλ β is obtained. The system has the following form: After introducing block matrices (1 2 = [1, 1] T ) and after including the mutually competing quantities being determined in combined vectors, i.e. having introduced vectors system of Eqs. (28) can also be expressed as follows: The solution to Eq. (32) is the combined Lagrange multiplier vector of the following form: After substituting Eq. (33) to conditions (26d) and (26e), jointly recorded as the following normal equation is obtained: The solution to this equation is vector which represents an evaluation of the combined vector of δX increments in the ( j + 1) iteration. In order to determine the combined vector X = X j + δX that is valid in this iteration, it will be taken into account that A = A j + E j , and thus, (35) then takes the following form: which yields the following: Based on Eqs. (26a) and (26b), jointly recorded as the combined residual vector that is valid in the ( j + 1) iteration can be determined: On the other hand, based on Eq. (26c), expressed in the following form: the following residual vector is determined: Taking that vector, we can obtainẼ j+1 , and hence, the matrix A j+1 = A −Ẽ j+1 . Variance coefficient estimators in TM split estimation should be derived by referring to the split EIV models (e.g. by applying the theory presented in Wiśniewski and Zienkiewicz (2021b). However, this problem that requires additional, detailed theoretical and empirical analyses is beyond the scope of this paper. With minor random disturbances of matrix A, variance coefficient estimators appropriate for M split estimation can also be recommended here. Then, in Eq. (8), the vectorsṽ α andṽ β should be replaced by the vectorsυ α = y − (A −Ẽ)X α andυ β = y − (A −Ẽ)X β . These estimators of variance coefficients stay unbiased; however, they loose their invariance according to growing values of EX α andẼX β .

Algorithm
Total M split estimation algorithm contains the following basic elements: Step 0-a starting step, Step 1-iterative calculation of M split estimators for valid split EIV models (with internal iterations l = 0, . . . , s), Step 2-updating the EIV models' parameters, and return to Step 1 (until the adopted criterion for stopping the iterative process, j = 0, . . . , k), Step 3-adopting the final values of Total M split estimators. Each of these steps is described in more detail below.
Step 0: Similar to M split estimation, the iterative process can also be initiated here using the following classical leastsquares (LS) estimatorŝ where W = Q −1 y is the weight matrix. Therefore, the following are adopted: Step 1: Calculate M split estimatorsX j α andX j β : step 1 (1) : Based on the valid vector υ j β(l) , the following weight matrix is constructed and the following is calculated: step 1 (2) : Based on the vector υ j α(l+1) , the following weight matrix is constructed and the following is calculated: Once criterion (49) has been satisfied, the following M split estimators in the j-th iteration: and residual vectors are adopted.
Step 2: For the determined iterative M split estimatorsX j α , X j β , and residual vectorsυ j α ,υ j β , weight matrices that are valid in the j-th iteration are determined: Based on these matrices and takingX j ⊗α =X j α ⊗ I n ,X j ⊗β = X j β ⊗ I n , we can determine matrix j , Eq. (30); then, the combined vector of increments is calculated: This vector enables the calculation of Lagrange multiplier vector that is valid in the ( j + 1) iteration and then the calculation of combined residual vectors and Based on vectorẽ j+1 , it is necessary to build a matrix of disturbancesẼ j+1 and to calculate the vector of the following parameters which are valid in the ( j + 1) iteration: Step 3: Repeat Step 1 and Step 2 until Once criterion (58) has been satisfied, it is assumed that TM split estimator of the combined vector of X = The blocks of this vector are TM split estimatorsX α andX β , the "removal" of which from vectorX may be facilitated by relationshipŝ where Once the iterative process is complete, the final residual vectorsυ =υ k andẽ =ẽ k are also determined. Based on vectorυ = [υ T α ,υ T β ] T , two versions of the residual vector corresponding to the observation vector y, i.e.
where D υ α = I n , 0 n, n and D υ β = 0 n, n , I n , are obtained. On the other hand, vectorẽ provides the basis for the construction of residual matrixẼ corresponding to matrix A, hence alsoÃ = A −Ẽ.
The given above process of iterative determination of TM split estimates is also described by the flowchart presented in Fig. 1

Example 1: competitive models of systematic errors
In one of the examples provided in a study by Wiśniewski (2010), it was assumed that y i , i = 1, . . . , n were observations of a certain value of Y disturbed not only with random errors v i but also with systematic errors s i = s(t i ) = a + bt i (e.g. Wiśniewski 1985;Kubáčková and Kubáček 1991;Yang and Zhang 2005). The problem, however, is that two versions of this model can be used: whereas it is not known which of them concerns specific observation y i . For this reason, in M split estimation, the classical observation model is split into the following models: where X = Y + a. In these models, two mutually competing versions of the parameters occur, namely X α = Y + a α , X β = Y +a β . Observations were simulated with the assumption of theoretical values of the parameters X α , b α and X β , b β . Theoretical observations y i , i = 1, . . . , 10, were affected by Gaussian errors with the expected value of 0, and standard deviation of σ y . For theoretical values X α = 6.0, b α = 0.5, X β = 3.0, b β = 1.0, and standard deviation σ y = 0.14, the set of observations presented in Table 1 was obtained. The table also presents M split estimates, namelyX α ,b α ,X β ,b β , the mutual competitive residualsv iα ,v iβ and weights related to such residuals A graphical illustration of the set of observations, and a graphical interpretation of the obtained results (as compared to LS-estimators determined using model (61)), are shown in Fig. 2. For the sake of clarity, the figure shows the competitive residuals only of the observation y 10 . In this figure, the mutually competing results of M split estimation are conventionally denoted as M split(α) and M split(β) (when it is convenient, these notations will also be used further on in this paper).
The models contained in Eq. (62) will now be replaced with EIV models of the following form: where e t i is a random error affected to the variable t i .
(the first column of matrix A is not random). Vector e = vec(E), built from matrix E columns, has the following form: In view of the structure of this vector, Q e cofactor matrix, similarly as in Schaffrin and Wieser (2008) and Shen et al. (2011), will be expressed in the following form: where Q e t is the cofactor matrix of vector e t . Note that here Q e t is regular, but Q 0 is not; thus, Q e is not regular too (the similar situation is in the example given by Schaffrin and Wieser 2008). Random errors e t i are simulated as Gaussian quantities with the expected value of 0, and standard deviation of σ e . M split and TM split estimators of model (63) parameters will be determined for four values of standard deviation: σ e = 0 (variant I), σ e = 0.13 (variant II), σ e = 0.28(variant III), σ e = 0.37(variant IV). In each of these variants, observations y i , i = 1, . . . , 10, are as in the example cited above (Table 1). The data adopted for the calculations are listed in Table 2, while the M split and TM split estimators obtained for these data are presented in Table 3. Additionally, the residualsẽ t i are presented in Table 4. Table  5 shows the competitive residuals and the respective weights for variant III. In variant I, since matrix A is not disturbed by random disturbances, TM split estimators are equal to M split estimators. With an increase in the σ e standard deviation value, the differences between these estimators increase, with TM split estimators remaining close to the theoretical values (as do M split estimators for the variant without random errors of matrix A, σ e = 0). This is well illustrated by the norm values of the vector of differences between the vector of theoretical    Table  3.
In general, the iterative process involved in the determination of TM split estimators ended after 4 to 6 steps of "external" iteration. In each of these steps, 6 to 7 "internal" iterations were carried out resulting in M split estimators that are valid for this step. The course of the iterative process in Total M split estimation, based on the example of variant II, is presented in Fig. 3.
The given examples apply one observation set, respectively. The additional analyses will be based on Monte Carlo (MC) simulations. The main objective is to determine the empirical accuracy of Total M split estimation and the measures of its efficacy.
The accuracy of M split estimation can be determined by applying asymptotical covariance matrices proposed in Wiśniewski and Zienkiewicz (2021b). The diagonal elements allow us to compute the estimated standard deviation σX α , σb α , σX β , σb β of the respective M split estimates. To apply such an approach to Total M split estimation, one should develop the theory presented in the paper mentioned, which is beyond the scope of the present paper. Based on simulated observation sets, an empirical way could be an alternative for the analytical assessment in question. Total M split estimatesX k α , b k α ,X k β ,b k β are computed for each k = 1, . . . , N simulation,  I  II  III  IV  I  II  III  IV   1   which is a base for determining the following MC-estimates of the parameters X α , b α ,X β , b β : The Monte Carlo estimators of the parameter standard deviations can be computed in the following way (e.g. Koch 2013;Nowel 2016;Lv and Sui 2020).
M split estimation and its several developments (including Total M split estimation) focus on optimal fitting the competitive functional models in the observation set. The efficacy of M split estimation, like the efficacy of other methods, can be described by the differences between the parameter estimates obtained and the actual parameter values. Considering N simulations, the efficacy is usually measured by the root mean squared error (RMSE) (e.g. Kargoll et al. 2018;Lv and Sui 2020). Here, the efficacy of M split or Total M split estimates is determined by following RMSEs Additionally, the global root mean squared error, concerning the whole parameter vector, is determined (e.g. Wiśniewski 2014) where r is the number of estimated parameters (here r = 4). The simulated errors v k i , e k t, i , i = 1, . . . , 10 (for each k = 1, . . . , N ) affect the theoretical observations y i or the elements a i, 2 of the second column of the matrix A (the theoretical observations and the theoretical matrix A stay the same). The simulations are performed by applying the Gaussian random generators of the MatLab system, σ y randn(n, 1) or σ e randn(n, 1), respectively.
First, let us examine the efficacy of M split estimates, which apply the models of (62). Since matrix A is constant, only observation errors are simulated. The computations are determined in several variants of the observation standard deviation, i.e. σ y = 0.05, 0.1, 0.2, 0.3. Table 6 presents results obtained for N = 3000 and the theoretical parameter values X = [X α , b α , X β , b β ] T = [6.0, 0.5, 3.0, 1.0] T . Figure 4 shows M split estimates obtained in each simulation and MC estimates (for σ y = 0.1).
The accuracy and efficacy of the estimatesb α andb β are the most satisfying. The values ofσ MĈ b α ,σ MĈ b β and RMSEb α , RMSEb β are relatively small for all values of the observation standard deviations. The values obtained for the estimateŝ X α ,X β are higher; however, they are still acceptable.
Let us now use the models (63) to examine how random disturbances of matrix A might influence the accuracy and efficacy of the estimates. The observation errors are simulated assuming the constant standard deviation σ y = 0.1, whereas the errors e t, i in several variants, in which σ e = 0, 0.05, 0.1, 0.2, 0.3. Table 7 presents the results for N = 3000.
Both accuracy and efficacy of M split estimates decrease when the coefficient matrix is disturbed with random errors. That effect can be reduced by using Total M split estimation, for which the measures in questions are smaller. As in the previous case, TM split estimates of the parameters b α and b β have the most satisfying accuracy and efficacy. The efficacy of Total M split estimation is confirmed by the parameter estimates obtained in each simulation. The example estimates and the MC estimates are presented in Fig. 5 (for σ e = 0.2). Schaffrin and Wieser (2008) as well as Shen et al. (2011) and Mahboub (2012) applied WTLS for the estimation of the intercept ξ 1 and slope ξ 2 of the regression line

Example 2: linear regression
Let it now be assumed that the set of observations not only contains observations concerning model (71) but also observations for which the regression line differs in parameters ξ 1 and ξ 2 (in Total M split estimation, these will be parameters ξ 1β and ξ 2β ). These observations will hereinafter be referred to conventionally as outliers. In contrast to the classical approach, the deviation here is of a different nature, and is not necessarily related to the existence of gross errors. If the assignment of observation y i to its respective regression line is not known, then this observation may correspond to both the following model and to the model that is competing in relation to it, namely: For the estimation of parameters in models (72) and (73), Total M split estimation will be applied. The calculations will be carried out using the data provided in (Neri et al. 1989) and also used in (Schaffrin and Wieser 2008;Shen et al. 2011;Mahboub 2012). These data will be supplemented with two variants of bias. In the first of these variants, regression line (73) with theoretical parameters ξ 1β = 2.0, ξ 2β = 0.75 is adopted, while in the second variant, ξ 1β = 4.5, ξ 2β = −0.70 is adopted. For the outliers, weights W x i were determined based on the information on the weights of coordinates x i concerning the group of original observations (where necessary, also using interpolation). Moreover, the outliers are assigned equal weights W y i = 50, which corresponds to the standard deviation σ y = 0.14. The values of these weights are very important to the success of Total M split estimation. Too high accuracy (high weights) of the added observations may cause the original observations to be "ignored", while too low accuracy (low weights) may cause the added observations to be "ignored". In the presented example, satisfying results that were not much different from each other were obtained for 40 < W y i < 80. Original observations (i = 1, . . . , 10) (Neri et al. 1989), outliers (i = 11, 12, 13), and the  weights W x i and W y i , corresponding to these observations, are listed in Table 8. The results of model (71) parameter estimation using WTLS for original observations, transcribed from Shen et al. (2011), are provided in columns 2 and 3 of Table 9 [the estimateσ 0 is computed by applying Eq. (19)]. Based on original observations and using models (72) and (73), TM split estimators of the parameters occurring there were calculated (column 4, Table 9). A graphical interpretation of the original set of observations and the location of regression lines (determined based on the WTLS and TM split estimators) in this set are shown in Fig. 6. On the other hand, the WTLS and TM split estimators determined for the sets   extended to include outliers are provided in other columns of Table 9. These sets and the corresponding regression lines are shown in Fig. 7. The TM split estimators determined for the original set of observations may appear not wholly satisfactory (column 4, Table 9). Due to the lack of outliers, Total M split estimation, predicting the existence of two mutually competing regression lines, forces two regression lines to fit into the set of observations (see Fig. 6). It should be noted that the regression line established by WTLS estimators lies between these lines. The average values of TM split estimatorsξ 1M split = (ξ 1α +ξ β )/2 = 5.4069 andξ 2M split = (ξ 2α +ξ 2β )/2 = −0.4612, as compared to WTLS estimatorŝ ξ 1 = 5.4799 andξ 2 = −0.4805, can already be considered satisfactory. For both sets containing outliers, TM split estimatorsξ 1α andξ 2α are close to the corresponding WTLS estimators obtained for the original set of observations. On the other hand, estimatorsξ 1β andξ 2β are close to the true values of parameters ξ 1β = 2.0, ξ 2β = 0.75(Variant I) and ξ 1β = 4.5, ξ 2β = −0.70(Variant II). In such cases, WTLS estimators yielded no good answers. This is particularly true Fig. 6 A set of original observations (Shen et al. 2011) and the regression line position were determined based on WTLS and TM split estimators for Variant II for which the relevant comparisons are particularly unfavourable. The results obtained using WTLS are not surprising, as the lack of WTLS estimators' robustness to observations is their inherent feature.
The example presented above concerned a situation where outliers can be assigned a regression line that is appropriate for them. In practice, however, the outlying of observations may relate to single observations and result from, for example, the effect of gross errors. In order to check the response of TM split and WTLS estimators to such errors, one of the observations from the original set will be affected by gross error with a few value versions. For example, let it be assumed that such an observation is y 5 = 3.5 (x 5 = 3.3) with weights W x 5 = 200 and W y 5 = 20 (Table 8). This observation will be affected by gross error with values of g = 1, g = 2, g = 5, g = 10, respectively. The data adopted for the calculations are provided in Table 10. On the other hand, the TM split and WTLS estimators of parameters ξ 1 and ξ 2 are provided in Table 11. TM split estimatorsξ 1α andξ 2α for each accepted gross error value are satisfactory, especially in comparison with WTLS estimators being obtained, which are unacceptable even for small gross errors. The quantitiesξ 1β andξ 2β are competing in relation toξ 1α andξ 2α . These are estimators of the parameters of the y 5 = ξ 1β + ξ 2β x 5 regression line on which the observation affected by gross error should lie. By using the equationỹ 5 =ξ 1β +ξ 2β x 5 , it is possible to calculate the prediction of observation y 5 affected by gross error. The prediction of this observation for the adopted gross error values, as compared to its simulated values, is presented in Table 12. A graphical interpretation of the obtained results is provided in Fig. 8.

Example 3: Two-dimensional affine transformation
The estimators of parameters in the EIV model, which are robust to gross errors, include inter alia robust total leastsquares (RTLS) and total least trimmed squares (TLTS) estimators (Wang et al. 2016;Lv and Sui 2020). One of the examples provided in Lv and Sui (2020) concerns a two-dimensional affine transformation carried out based on observations affected by gross errors. In that case, the authors applied the following transformation model: where u s , υ s and u t , υ t are the coordinates of the common points in the start and target coordinate systems, while a 1 , b 1 , c 1 , a 2 ,b 2 , c 2 are the parameters being determined. Table 13 presents 15 observation points simulated in the start and target coordinate systems (Lv and Sui 2020). Some of these observations (highlighted in bold) are affected by gross errors. For these data, Lv and Sui (2020) applied TLTS estimation (using two different algorithms yielding the same results) as well as TLTS with RTLS as a starting point (hereafter denoted as TLTS/RTLS). The authors compared the obtained estimators with WTLS and RTLS estimators. These results may now be supplemented with TM split estimators (the seventh column of Table 14). TM split estimators will also be calculated for the data without gross errors. All the cited and calculated estimators are provided in Table 14.  Table 14 shows that TM split(α) estimates are generally the closest to TLTS estimates. However, in the case of parameters a 1 , b 1 , c 1 and b 2 , TM split(α) estimators are closest to TLTS/RTLS estimators and thus to the true parameter values. The interpretation of TM split(β) estimators is similar to that in the previous example, i.e. these are estimators of the parameters of the model for observations affected by gross errors. It is worth noting here that in the absence of gross errors, both versions of TM split estimators differ slightly from each other.

Summary
By using M split estimation, it is possible to determine the estimators of mutually competing parameters in classical functional models. However, there are cases in geodetic practice in which classical models need to be replaced with EIV models. The method proposed in this paper, called "Total M split estimation", is an development of M split estimation that accepts such models. The Total M split estimation objective function was determined by applying the Lagrange approach, as in the case of WTLS method. The mutually competing Gross error Prediction of an observation Simulated observation g = 1ỹ 5 =ξ 1β +ξ 2β x 5 = 6.42 − 0.60 · 3.3 = 4.44 y 5 = 4.5 g = 2ỹ 5 =ξ 1β +ξ 2β x 5 = 9.44 − 1.24 · 3.3 = 5.34 y 5 = 5.5 g = 5ỹ 5 =ξ 1β +ξ 2β x 5 = 11.85 − 1.04 · 3.3 = 8.41 y 5 = 8.5 g = 10ỹ 5 =ξ 1β +ξ 2β x 5 = 19.25 − 1.80 · 3.3 = 13.31 y 5 = 13.5 EIV models occurring in this function were replaced with their linear approximations. This enabled the construction of a relatively simple yet efficient algorithm for determining TM split estimators. The basis of this algorithm is the iterative updating of EIV models (external iterations) based on M split estimators (internal iterations) obtained in the previous iterative step. The proposed algorithm is efficient in all cases presented in the paper. It concerns both results and the flow of the iterative process. The problems with the convergence of the external iterations might occur because of the linear approximation of EIV models applied. It is especially evident when the errors disturbing the matrix A are too big.
The examples presented in the paper showed that the properties of TM split estimators are, in general, similar to the properties of M split estimators. If the elements of matrix A are not affected by random errors, then TM split and M split estimators are equal to each other. The possibility of determining estimators of competing parameters is of particular importance when the sets of observations are a mixture of the realisations of two random variables with mutually competing positional parameters. Such a situation occurs in Examples 1 and 2, where each observation group can be assigned a corresponding regression line. The problem, however, is that it is not known which of these lines is the best for a particular observation. TM split estimators, similarly like M split estimators for classical models, yield satisfactory results here. Due to their theoretical origin (neutral LSmethod), WTLS estimators are not robust to gross errors. Where the sets are realisations of only a single random variable, Total M split estimation further offers two mutually competing solutions. These are forced solutions, yet so close to each other that even in such a situation, it is possible to evaluate the functional model parameters (e.g. after the calculation of average values of respective estimators). Total M split estimation can also be applied for the estimation of EIV model parameters in the case where the outlying of observations results from their being affected by gross errors. From the perspective of the M split and TM split estimation theory, such a case is not significantly different from that discussed earlier. In the second part of Example 2, it was shown that the determined TM split estimators enabled the determination of not only the regression line appropriate for "good" observations (TM split(α) solution), but also of the regression line on which the observation affected by gross error lies (TM split(β) solution).
In EIV models, it is assumed that matrix A is observed as well. Therefore, its elements can also be affected by gross errors arising from various reasons. Such a situation is the case in Example 3, in which certain coordinates are affected by gross errors, both in the start and target systems in twodimensional affine transformation. In this example, TM split estimators are close to robust TLTS estimators and, for certain transformation parameters, they are also close to the results of TLTS estimation with RTLS estimation as a starting step.
The estimates' accuracy is an important issue (especially in comparing the estimation methods). The accuracy of M split estimation can be determined by applying asymptotical covariance matrices. To apply such an approach to Total M split estimation, one should develop the theory, which is beyond the scope of the present paper. Section 4.1 presents the assessments of the M split and Total M split estimates accuracy (empirical standard deviations) obtained from the Monte Carlo simulations. Generally, the accuracy of M split estimates decreases with the growing standard deviation of errors disturbing the matrix A. Total M split estimates have smaller standard deviations than respective M split estimates in such a context. Similar relations concern the measures of efficacy,  namely values of RMSEs. Thus, when the competitive functional models are supplemented with EIV models, then M split estimation should be replaced by Total M split estimation. It is especially advisable when disturbances of matrix A have large standard deviations (here, the application of Total M split estimation is justified for σ e = σ y ).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.