Improving the weighted least squares estimation of parameters in errors-in-variables models

Available online at www.sciencedirect.com

Journal of the Franklin Institute 356 (2019) 8785–8802 www.elsevier.com/locate/jfranklin

Improving the weighted least squares estimation of parameters in errors-in-variables models Peiliang Xu Disaster Prevention Research Institute, Kyoto University, Uji, Kyoto 611-0011, Japan Received 18 October 2018; received in revised form 17 March 2019; accepted 20 June 2019 Available online 12 July 2019

Abstract Although the weighted least squares (LS) method is straightforward to deal with errors-in-variables (EIV) models, it results in the biased estimates of parameters and the variance of unit weight. The total least squares (TLS) method is statistically rigorous and optimal but is computationally much more demanding. This paper aims at constructing an improved weighted LS estimate of parameters from the perspective of multiplicative error models, which is expected to mainly maintain the advantages of computational simplicity of the weighted LS method and the optimality of the weighted TLS method to some extent but almost remove the bias of the weighted LS estimate and free the computational burden of the TLS method. The statistical aspects of the weighted LS estimate developed from the perspective of multiplicative error models have been analyzed and an almost unbiased estimate of the variance of unit weight proposed. Finally, an N-calibrated weighted LS estimate has been constructed from the perspective of multiplicative error models. © 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

1. Introduction Given a measurement system, if both the responses y and the independent variables A are contaminated with random errors, then such a measurement system is well known as an errors-in-variables (EIV) model, which is usually written as follows: y = (A − EA )β + ,

E-mail address: [email protected]. https://doi.org/10.1016/j.jfranklin.2019.06.016 0016-0032/© 2019 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

(1a)

8786

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

(see e.g., [4,11,24,27,30]), where y is an (n × 1) vector of measurements with the random errors , A is a measurement matrix of size (n × t), and EA is the random matrix of A, and β is an unknown vector of parameters. The stochastic models of both  and EA are assumed as follows:  E () = 0, D() = Qy σ 2 , (1b) E (EA ) = 0, D(A ) = Qa σa2 , where E(·) and D(·) stand for the expectation and variance-covariance operators, respectively, A stands for the vector by sequentially placing the columns of the random error matrix EA underneath each other or symbolically A = vec(EA ) [16]. Here Qy and Qa are the cofactor matrices of y and A, respectively. Qy is generally assumed to be positive definite and often alternatively replaced with the weighting matrix W, namely, Qy = W−1 . For technical simplicity but without loss of generality, y and A are assumed to be stochastically independent in this paper, namely, E (, A ) = 0. To estimate the unknown parameters β, one often makes an additional assumption that both σ 2 and σa2 are given or σ 2 = σa2 in case that σ 2 and σa2 are unknown. In this paper, we follow the second part of the assumption, namely, σ 2 = σa2 but unknown. σ 2 is also called the unknown variance of unit weight. There can be five different classes of methods to estimate the unknown parameters β in the EIV model (1): least squares (LS), orthogonal regression, minimizing the Frobenius norm of the residuals to both A and y, maximum likelihood and total least squares (TLS). Except for the LS method (see e.g., [6,15,29,30]), all the other four classes of methods apply the corrections to both the measurements y and A from the point of view of approximation theory by minimizing the Frobenius norm of the residuals to both A and y [13,27] or applying orthogonal regression [1,2,5,19,20]. The corrections to both y and A can also be adjusted from the statistical points of view of maximum likelihood and the weighted TLS method [3,7–9,12,17,21–23,25,31,33]. For some more details, the reader is referred to [30]. The weighted LS method is easy to compute and estimate the unknown parameters β in the EIV model (1), if the matrix A is treated as if it were not random, namely, βˆ wls = (AT WA )−1 AT Wy,

(2)

where βˆ wls is the weighted LS estimate of β. Nevertheless, the estimate βˆ wls of β in association with the EIV model (1) can result in a number of unwanted consequences. The weighted LS estimates of β itself and the variance of unit weight can be significantly biased (see e.g., [6,15,29,30]). The effect of errors-in-variables can also significantly affect the estimation of variance components of y [29,30]. Although the weighted TLS method is most rigorous statistically to estimate the parameters β in the EIV model (1), we have to numerically solve the following nonlinear equalityconstrained minimization problem: T min: S(va , r, β) = vaT Q−1 a va + r Wr,

(3a)

subject to the nonlinear equality conditions y + r = (A + VA )β,

(3b)

if Qa is regularly invertible, where r and VA are the correction vector and matrix of y and A, respectively, and va = vec(VA ). The weighted TLS problem (3) can also be reformulated as a mathematically equivalent unconstrained nonlinear adjustment model, as can be seen in [33] (see also [30]). The simultaneous estimation of both β and the true but

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8787

unknown values of A in the minimization problem (3) has substantially increased the computational complexity of TLS estimation (see e.g., [29,30,35]). For example, in this big data era, let us assume n = 10, 000, 000 and t = 1, 000. If all the elements of A are stochastically independent variables and if we apply the weighted TLS method to the EIV model (1), we will have to simultaneously estimate a total number of (1.0 × 1010 + 1, 000) unknown parameters from the measurements y and A. This computational burden can still be rather heavy and beyond our computational capability, even though the modern fastest supercomputers in the world are now able to handle 2.0 × 1017 flops per second (https: // www.en.wikipedia.org/ wiki/ Summit _(supercomputer)). From this point of view, a simple LS estimate still has its computational advantage over a statistically strict TLS solution in a certain sense. The major purpose of this paper is to construct an improved weighted LS estimate of β in the EIV model (1), which aims at making a best possible compromise between the weighted LS and TLS methods by maintaining the advantages of both methods but minimizing their statistical and computational disadvantages. In other words, the improved weighted LS estimate of parameters in the EIV model (1) should maintain the computational simplicity of the weighted LS method but avoid its poor statistical performances. The statistical aspects of the improved estimate to be constructed should also be TLS-comparable but without having to simultaneously estimate the true values of A. Xu [29,30] directly started with the weighted LS solution of β, analyzed the source of bias and then constructed an N-calibrated almost unbiased estimator for β. Accordingly, he was also able to successfully construct some almost unbiased estimators for the variance of unit weight. In this paper, we will first reformulate the EIV model (1) into a multiplicative error model and then estimate the parameters β from the perspective of multiplicative error models. More precisely, the paper is organized as follows. Section 2 will reformulate the EIV model (1) into a multiplicative error model. We will then apply the weighted LS method to estimate the parameters β. Section 3 will focus on the statistical aspects of the weighted LS estimate developed from the perspective of multiplicative error models. We will analyze the bias of the weighted LS solution, derive the accuracy formula up to the first order approximation and then follow [32] (see also [29]) to construct a bias-corrected weighted LS estimate of β. An almost unbiased estimate of the variance of unit weight will also be given here. Section 4 will use examples to illustrate the estimation methods developed in this paper. 2. Weighted LS estimation of parameters in the perspective of multiplicative error models To reformulate the EIV model (1) as a multiplicative error model, we construct the new pseudo-measurements by subtracting Aβ from y, which is denoted by l and given as follows: l = y − Aβ.

(4)

With the stochastic model (1b) in mind, we have the expectation of l: E (l) = E (y) − E (A )β = 0.

(5a)

Applying the variance-covariance propagation law to Eq. (4), we obtain the variancecovariance matrix of l as follows: D(l) = Qy σ 2 + (βT  In )Qa (β  In )σ 2

8788

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

= Ql σ 2 .

(5b)

where Ql = Qy + ( β T  In )Qa ( β  In ),

(5c)

Ql is the cofactor matrix of l,  stands for Kronecker product of two matrices (see e.g., [16]). Since Qy is assumed to be positive definite, Ql is positive definite as well and regularly invertible. If the column vectors of A are further assumed to be stochastically independent, with each column being of a different cofactor matrix Qia , then the stochastic model (5b) becomes   t  D(l) = σ 2 Qy + Qia βi2 . (6) i=1

The stochastic model (6) of the pseudo-measurements l is of a typical multiplicative error nature (see e.g., [25,26,34]). The random errors of l are obviously proportional to the model parameters β. The larger the parameters β, the noisier the pseudo-measurements l. We should also note that although Ql of l depends on β, the accuracy of l is not proportional to the magnitudes of the measurement signals themselves, namely, E(A)β. We are now in a position to estimate the parameters β in the EIV model (1) from the perspective of the multiplicative error models (4), (5a) and (5b). Although quasi-likelihood has been well known as a standard statistical method to handle multiplicative error models (see e.g., [14,18,28]), a major disadvantage of the method is that its statistical foundation is formally established on the assumption that the corresponding likelihood-like function is the solution to the differential equation as derived by maximizing the conventional likelihood, which even cannot necessarily be normalized into a conventional distribution function, as explained in details in [34]. Actually, from the computational point of view, the quasi-likelihood solution is simply to solve a set of nonlinear estimation equations. If the nonlinear system of equations is of multiple solutions, we do not know either which solution is the correct quasilikelihood solution due to the lack of a clearly defined objective function. Thus, in this paper, we will follow the LS-based approach of [32] (see also [26,34]) to estimate the parameters β in the EIV model (1) from the perspective of multiplicative error models. 2.1. The ordinary LS method The simplest method to estimate the unknown parameters β in the multiplicative error models (4), (5a) and (5b) is to ignore the complicated stochastic model (5b) but directly apply the ordinary LS principle to the pseudo-measurements (4). Since the pseudo-measurements l of Eq. (4) are actually random errors, the ordinary LS principle is equivalent to minimizing the following objective function: min : F1 (β) = lT l = (y − Aβ)T (y − Aβ).

(7)

Thus, we can readily obtain the ordinary LS solution to the minimization problem (7) as follows: βˆ ls = (AT A )−1 AT y,

(8)

which, though derived from the perspective of multiplicative error models, is clearly the same as that by directly applying the ordinary LS principle to the original EIV model (1).

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8789

Accordingly, the ordinary LS estimates of the pseudo-measurements l and the variance of unit weight are given as follows: ˆl = y − Aβˆ ls ,

(9)

and 2 = (y − Aβˆ )T (y − Aβˆ )/(n − t ), σ ls ls ls

(10)

respectively. One may also alternatively estimate the variance of unit weight with the analysis method of [25]. The advantage of the ordinary LS estimate (8) is easy to compute. However, since all or some of the elements of A are measured but not exactly given, βˆ ls can be significantly affected by the random errors of A. Hodges and Moore [15] studied the sensitivity of the LS estimate βˆ ls to the elements of A. They further went ahead to derive the bias and approximate accuracy of βˆ ls after the assumption that the elements of A are all stochastically independent. Davies and Hutton [6] further studied the effect of the random errors of A on the estimate of variance of unit weight. The asymptotic properties of the ordinary LS estimates of the parameters and the variance of unit weight were also investigated. Recently, Xu [29,30] systematically investigated the effects of EA on the weighted LS estimates of β, σ 2 and the MINQUE variance component estimation. Bias-corrected weighted LS estimates were also constructed for both β and σ 2 . 2.2. The weighted LS method in the perspective of multiplicative error models Since the elements of l are clearly not of the same accuracy, the ordinary LS principle cannot be statistically optimal, as can be seen from the stochastic model (5b). A more natural way to estimate β in the multiplicative error models (4), (5a) and (5b) would have to take the variance-covariance matrix (5b) into account. In this case, the weighted LS principle is to minimize the following objective function: min : F2 (β) = lT Q−1 l l

 −1 = (y − Aβ)T Qy + (βT  In )Qa (β  In ) (y − Aβ).

(11)

Differentiating F2 (β) of the objective function (11) with respect to β and letting the derivatives equal to zero, we have ˆ −1 (Aβˆ − y) + Gl (Aβˆ − y) = 0, 2AT Q l

(12)

ˆ the ˆ l is the estimate of Ql at the point of β, where βˆ is the weighted LS estimate of β, Q matrix Gl is given as follows: ⎡ ⎤ (Aβˆ − y)T Pˆ 1 ⎢ ˆ ⎥ ⎢(Aβ − y)T Pˆ 2 ⎥ ⎢ ⎥, Gl = ⎢ (13a) .. ⎥ ⎣ ⎦ . (Aβˆ − y)T Pˆ t and

  ∂Q−1  l Pˆ i =  ∂βi 

β=βˆ

(i = 1, 2, . . . , t ).

(13b)

8790

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

The derivative matrices Pˆ i of (13b) can be computed by following the definition of differentiation of an inverse matrix, i.e. Theorem 3 of ([16], p.151). As a result, we can readily obtain ˆ More precisely, the matrices Pˆ i the derivatives of Q−1 with respect to β i at the point of β. l in (13b) are given as follows:   ˆ −1 ˆ −1 (eiT  In )Qa (βˆ  In ) + (βˆ T  In )Qa (ei  In ) Q Pˆ i = −Q l l ⎧ ⎫ t ⎨ ⎬ ˆ −1 , ˆ −1 ˆ j (Qi j + Q ji ) Q = −Q β (13c) l ⎩ ⎭ l j=1

where Qij stands for the cofactor matrix between the ith and the jth column vectors of A, Qi j = QTji and Qii = Qia , and ei is the ith natural basis vector of dimension t. With the formulae (13a) to (13c) in hands, we can now rewrite the nonlinear system (12) of normal equations as follows: ˆ −1 A )βˆ − AT Q ˆ −1 y − G(Aβˆ − y) = 0, (AT Q l l

(14a)

where ⎡

⎤ t  −1 T ˆ −1 ˆ ˆ ˆ β j (Q1 j + Q j1 )Ql /2⎥ ⎢(Aβ − y) Ql j=1 ⎢ ⎥ ⎢ ⎥ t  ⎢(Aβˆ − y)T Q −1 −1 ˆ ˆ /2⎥ βˆ j (Q2 j + Q j2 )Q ⎢ ⎥ l l ⎥. j=1 G=⎢ ⎢ ⎥ .. ⎢ ⎥ ⎢ ⎥ . ⎢ ⎥ t  ⎣ ˆ −1 −1 ˆ ˆ /2 ⎦ ˆ j (Qt j + Q jt )Q (Aβ − y)T Q β l l

(14b)

j=1

Because the elements of the cofactor matrix Ql are the functions of the parameters β, the normal equations (14a) are essentially nonlinear. As a result, we will have to use numerical ˆ Since the nonlinear normal equations are methods to compute the weighted LS estimate β. derived from the nonlinear weighted LS principle, the Gauss-Newton method has been commonly used to solve such nonlinear equations ([10], chapter 10), in particular, if residuals are small. As a result, we can iteratively solve for the solution βˆ to Eq. (14a) as follows: ˆ −1 A )−1 {AT Q ˆ −1 y + Gk (Aβˆ k − y)}, k = 0, 1, . . . βˆ k+1 = (AT Q lk lk

(15)

ˆ lk and Gk stand for Ql and G but computed at the point of (compare [10], p.90), where Q βˆ k , respectively, and the script k stands for the kth iteration. Nevertheless, if the residuals are ˆ −1 A ) is ill-conditioned, one may have to use the damped Gauss-Newton too large or if (AT Q lk method or global optimization techniques (see e.g., [33]). The variance of unit weight can now be estimated as follows:  2 = (y − Aβ) ˆ TQ ˆ /(n − t ). ˆ −1 (y − Aβ) σwls l

(16)

In particular, if the columns of A are stochastically independent mutually and are of the stochastic model (6), the formulae (13c) and (14b) become ˆ −1 Qia Q ˆ −1 , Pˆ i = −2βˆi Q l l

(17)

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

and

⎤ ˆ −1 Q1a Q ˆ −1 βˆ1 (Aβˆ − y)T Q l l ⎢ˆ ˆ −1 Q2a Q ˆ −1 ⎥ ⎢β2 (Aβˆ − y)T Q l l ⎥ ⎥, G=⎢ .. ⎢ ⎥ ⎣ ⎦ . ˆ −1 Qta Q ˆ −1 βˆt (Aβˆ − y)T Q

8791

⎡

l

(18)

l

respectively. 3. Statistical and computational analysis The nonlinear normal equations (14a) obviously reflect two mixed features: one as a result of multiplicative error models in the form of the third term on the left hand side of Eq. (14a) and the other attributed to the effect of errors in variables, namely, the random matrix A, on the weighted LS estimate, as is immediately clear by comparing (14a) with formula (16) of [34] and the results in Section 4 of [29]. Thus, the weighted LS estimate βˆ of β is expected to be significantly affected by these two different types of error sources. In this section, we will follow [34] and [29] to study the statistical aspects of the weighted LS ˆ as computed iteratively with Eq. (15). More specifically, we will work out its estimate β, bias, its accuracy of first order approximation and then with these results, further construct a bias-corrected weighted LS estimator for β. 3.1. Bias analysis Let us first denote l = l, namely, l = l =  − EA β,

(19)

since the expectation of l is a zero vector. We further denote the left hand side of Eq. (14a) by ˆ = (AT Q ˆ −1 A )βˆ − AT Q ˆ −1 y − G(Aβˆ − y). L(y, A, β) l l

(20)

For convenience, we also denote the true values of y and A by y and A, respectively, which are actually equal to the expectations of y and A, namely, y = E (y) and A = E (A ), if their elements are directly observable. Obviously, if the true values of y, A and βˆ are all used in the expression (20), it is trivial to show that L(y, A, β) = 0,

(21)

since the true values of direct observables are the same as their expectations. Following [29], we can now write the weighted LS estimate βˆ in terms of the first and second differentials as follows: 1 ˆ βˆ = β + d βˆ + d 2 β, (22) 2 up to the second order approximation, where d βˆ and d 2 βˆ are the first and second differentials ˆ To represent the weighted LS estimate βˆ in terms of the measurement errors EA and , of β. we mean that we expand βˆ at the point β and truncate the expansion up to the second order approximation, which is equivalent to computing the first and second differentials of d βˆ and

8792

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

d 2 βˆ at the true values (y, A, β). As a result, Eq. (22) can be rewritten in terms of the random errors of measurements: βˆ = β + β + b,

(23)

where β and b are unknown error vectors, representing the linear and second order approximations of βˆ in terms of EA and , respectively. Obviously, to connect d 2 βˆ with b, we have to replace d 2 βˆ with 2b. The basic technique to derive the biases of the weighted LS estimate βˆ is to work out the expressions for both β and b. All the terms higher than the second order will be truncated. To ˆ in terms of dy, dA and d βˆ by computing start with, we represent the functionals L(y, A, β) ˆ the differentials of L(y, A, β) up to the second order approximation, namely, ˆ = dL(y, A, β) ˆ + 1 d 2 L(y, A, β) ˆ . L(y, A, β) 2

(24)

ˆ dA and dy with We then respectively replace y, A and βˆ with their true values, and d β, their error representations, namely, β , EA ,  and/or l , and finally obtain the second order ˆ in terms of random errors. representation of the differentials of L(y, A, β) ˆ of (20), we have Applying the differential operator to L(y, A, β) ˆ = dAT Q ˆ −1 (Aβˆ − y) + AT d Q ˆ −1 (Aβˆ − y) dL(y, A, β) l l T ˆ −1 ˆ ˆ + A Q (d Aβ + Ad β − d y) l

− dG(Aβˆ − y) − G(d Aβˆ + Ad βˆ − d y).

(25)

ˆ −1 can be According to the differential of an inverse matrix ([16], Theorem 3, p.151), d Q l written as follows: ˆ −1 = −Q ˆ lQ ˆ −1 ˆ −1 d Q dQ l l l   ˆ −1 ˆ −1 (d βˆ T  In )Qa (βˆ  In ) + (βˆ T  In )Qa (d βˆ  In ) Q = −Q l l ˆ −1 = −Q l

t  t 

ˆ −1 . (d βˆi βˆ j + βˆi d βˆ j )Qi j Q l

(26)

i=1 j=1

ˆ : ˆ −1 of (26) into (25) yields the first differential of L(y, A, β) Inserting d Q l ˆ = dAT Q ˆ −1 (Aβˆ − y) + AT Q ˆ −1 (d Aβˆ + Ad βˆ − d y) dL(y, A, β) l l ˆ −1 − AT Q l

t  t 

ˆ −1 (Aβˆ − y) (d βˆi βˆ j + βˆi d βˆ j )Qi j Q l

i=1 j=1

− dG(Aβˆ − y) − G(d Aβˆ + Ad βˆ − d y).

(27)

On the basis of the differential (27), we can continue to compute the second differential of ˆ by applying the differential operator to the first differential dL(y, A, β) ˆ , which L(y, A, β)

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8793

can be symbolically written, after neglecting or omitting all the terms higher than the second order, as follows: ˆ = d (d L(y, A, β)) ˆ d 2 L(y, A, β) T ˆ −1 ˆ −1 (dAβˆ + Ad βˆ − dy) = d A d Q (Aβˆ − y) + dAT Q l

l

ˆ −1 (d Aβˆ + Ad βˆ − d y) + AT d Q ˆ −1 (d Aβˆ + Ad βˆ − d y) + dA Q l l T ˆ −1 2ˆ ˆ ˆ + A Q (d Ad β + d Ad β + Ad β) T

l

ˆ −1 (d βˆ T  In )Qa (βˆ  In )Q ˆ −1 }(Aβˆ − y) − d {AT Q l l ˆ −1 (d βˆ T  In )Qa (βˆ  In )Q ˆ −1 (d Aβˆ + Ad βˆ − d y) − AT Q l l T ˆ −1 ˆ T ˆ ˆ −1 }(Aβˆ − y) − d {A Ql (β  In )Qa (d β  In )Q l T ˆ −1 ˆ T −1 ˆ ˆ − A Ql (β  In )Qa (d β  In )Ql (dAβˆ + Ad βˆ 2 ˆ ˆ ˆ

− dy)

− d G(Aβ − y) − dG(dAβ + Ad β − dy) ˆ . − d G(d Aβˆ + Ad βˆ − d y) − G(d Ad βˆ + d Ad βˆ + Ad 2 β)

(28)

We can now expand the left hand side of the nonlinear normal equations (14a) at the ˆ dA and dy are replaced with point (y, A, β) up to the second order approximation. d β, their corresponding error representations, namely, β , EA ,  and/or l . Bearing in mind that Aβ − y = 0 and G = 0 at (y, A, β) (compare (14b)), we have ˆ = L(y, A, β) + dL(y, A, β) + 1 d 2 L(y, A, β) L(y, A, β) 2 T T −1 = A Q−1 { E β + A  − } + E A β A Ql {EA β + Aβ − } l  T  1 T (β  In )Qa (β  In ) + (βT  In )Qa (β  In ) Q−1 − A Q−1 l l (EA β+ Aβ −) 2 T + A Q−1 l (EA β + Ab ) 1 T −1 T − A Ql (β  In )Qa (β  In )Q−1 l (EA β + Aβ − ) 2 1 T T −1 − A Q−1 l (β  In )Qa (β  In )Ql (EA β + Aβ − ) 2 − G (EA β + Aβ − ) T

T −1 = A Q−1 l {EA β + Aβ − } + EA Ql {EA β + Aβ − } T

+ A Q−1 l (EA β + Ab ) T

T −1 − A Q−1 l (β  In )Qa (β  In )Ql (EA β + Aβ − ) T

T −1 − A Q−1 l (β  In )Qa (β  In )Ql (EA β + Aβ − )

− G (EA β + Aβ − ),

(29a)

ˆ of (27) and d 2 L(y, A, β) ˆ of (28) where dL(y, A, β) and d 2 L(y, A, β) stand for dL(y, A, β) but are computed at the point (y, A, β), respectively. G is given as follows:

8794

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

⎡

⎤ t  −1 T −1 (E β + A  − ) Q β (Q + Q ) Q / 2 β j 1j j1 l l ⎢ A ⎥ j=1 ⎢ ⎥ ⎢ ⎥ t ⎢(E β + A − )T Q−1  β (Q + Q )Q−1 /2⎥ β j 2 j j2 ⎢ A ⎥ l l ⎥. j=1 G = ⎢ ⎢ ⎥ .. ⎢ ⎥ ⎢ ⎥ . ⎢ ⎥ t  ⎣ ⎦ −1 (EA β + Aβ − )T Q−1 β (Q + Q ) Q / 2 j t j jt l l

(29b)

j=1

ˆ of (20) must satisfy the normal equations (14a), namely, L(y, A, β) ˆ = Because L(y, A, β) 0, all the linear terms of β , EA ,  and/or l in (29a) should satisfy the zero constraint, namely, T

A Q−1 l {EA β + Aβ − } = 0, or equivalently,  T −1 T β = A Q−1 A A Q−1 l l ( − EA β)  T −1 T = A Q−1 A Q−1 l A l l −1

T

= Nl A Q−1 l l ,

(30)

where T

Nl = A Q−1 l A. With β of (30) in hands, we can write   −1 T l EA β + Aβ −  = − I − ANl A Q−1 l = −(I − Hl )l ,

(31)

where −1

T

Hl = ANl A Q−1 l . In the similar manner, all the second order terms of β , EA ,  and/or l in (29a) must satisfy the zero constraint of the normal equations (14a) as well. Bearing β of (30) in mind, we obtain the following equations with the second order terms of random errors, namely, T

T

−1 −1 0 = −EAT Q−1 l (I − Hl )l + A Ql EA β + A Ql Ab T

+ A Q−1 l

t  t 

(iβ β j + βi  jβ )Qi j Q−1 l (I − Hl )l − Gl (I − Hl )l ,

(32a)

i=1 j=1

from which we can now represent the second order term b of the weighted LS estimate βˆ in terms of the random errors as follows: T

−1

−1 Nl b = EAT Q−1 l (I − Hl )l − A Ql EA β + Nl Gl (I − Hl )l T

− A Q−1 l

t  t  i=1 j=1

(iβ β j + βi  jβ )Qi j Q−1 l (I − Hl )l ,

(32b)

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

where  iβ is the ith element of β and the error matrix Gl is given below: ⎡ ⎤ t  −1 T T −1 β j (Q1 j + Q j1 )Ql /2⎥ ⎢l (I − Hl ) Ql j=1 ⎢ ⎥ ⎢ ⎥ t ⎢T (I − H )T Q−1  β (Q + Q )Q−1 /2⎥ l j 2 j j2 ⎢ l ⎥ l l ⎥. j=1 Gl = ⎢ ⎢ ⎥ .. ⎢ ⎥ ⎢ ⎥ . ⎢ ⎥ t  ⎣ T ⎦ −1 l (I − Hl )T Q−1 β (Q + Q ) Q / 2 j tj jt l l

8795

(32c)

j=1

ˆ which can be obtained by We can now derive the biases of the weighted LS estimate β, ˆ , we applying the expectation operator to b of Eq. (32b). Denoting the bias of βˆ by bias(β) have ˆ = Nl E (b ) Nl bias(β) T

−1 = E {EAT Q−1 l (I − Hl )l } − A Ql E (EA β ) + E {Gl (I − Hl )l } T

− A Q−1 l

t  t    E (iβ β j + βi  jβ )Qi j Q−1 l (I − Hl )l i=1 j=1

= m1 + m2 + m3 + m4 ,

(33a)

where m1 = E {EAT Q−1 l (I − Hl )l }, T

(33b)

m2 = −A Q−1 l E (EA β )

(33c)

m3 = E {Gl (I − Hl )l },

(33d)

T

m4 = −A Q−1 l

t  t    E (iβ β j + βi  jβ )Qi j Q−1 l (I − Hl )l .

(33e)

i=1 j=1

Since the second order statistics of EA ,  and/or l are all given, all the four terms of m1 , m2 , m3 and m4 are computable, which are given in the appendix. Thus, the biases of the weighted LS estimate βˆ can be symbolically written as follows: ˆ = N−1 (m1 + m2 + m3 + m4 ). bias(β) l

(34)

ˆ of the weighted LS estimate βˆ are It is clear from formula (34) that the biases bias(β) attributed to four terms m1 , m2 , m3 and m4 . These four terms (33b) to (33e) clearly come from two different sources of errors. When compared with the results of biases in [29], the first two terms m1 and m2 have a source in the random errors of the design matrix A. The third term m3 is solely due to the multiplicative error nature of the random errors of l, as can be seen from [34] and [26]. An interesting point is that the multiplicative nature (5c) of the transformed observational model (4) is simply caused by the errors in variables. As a result, this mixed effect has directly been reflected in the last term m4 as a combination of

8796

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

ˆ −1 A ) of the differential of Ql of (5c) and the random errors EA in the normal matrix (AT Q l Eq. (14a). On the other hand, it is also clear from the appendix that the terms m1 and m4 are proportional to the noise level σ 2 , while the other two terms m2 and m3 are not only proportional to the noise level σ 2 but also the magnitudes of the parameters β. 2 of the estimate (16) up to the With the error representation (31) in hands, we expand  σwls second order approximation of random errors, apply the expectation operator to the expansion and obtain ˆ TQ ˆ /(n − t )} ˆ −1 (y − Aβ) E ( σ 2 ) = E {(y − Aβ) l

wls

= E {( − EA β − Aβ )T Q−1 l ( − EA β − Aβ )}/ (n − t ) = E {Tl (I − Hl )T Q−1 l (I − Hl )l }/ (n − t ) 2 =σ ,

(35)

indicating that the weighted LS estimate of the variance of unit weight from the perspective of multiplicative error models is unbiased up to the second order approximation. This result of unbiasedness is surprisingly different from the corresponding biased estimates of the variance of unit weight with the ordinary or naive weighted LS estimates of parameters, as reported in [6,15] and [29]. 3.2. First order accuracy and a bias-corrected weighted LS estimator To derive the first order accuracy of the weighted LS estimate βˆ for the parameters in the EIV model (1) from the perspective of multiplicative error models, we can directly apply the error propagation law to the linear representation (30) of the random errors of βˆ in terms of l , which is given as follows: −1

ˆ = N σ 2. D(β) l

(36)

By comparing formula (36) with formula (11) of [29], we can see that these two formulae of the first order accuracy are different. Actually, unlike the weighted LS estimate βˆ as derived from the normal equations (14a) in this paper, the naive weighted LS estimate of β [29] or the ordinary LS estimate (see e.g., [6,15]) does not fully take both the errors in y and A into account to construct a proper weighting matrix and then to estimate β in the EIV model (1). Theoretically, the weighted LS estimate βˆ from the perspective of multiplicative error models in this paper should be more appropriate than the weighted and/or ordinary LS estimates reported in [6,15] and [29]. Since the third bias term m3 of βˆ is completely attributed to the nature of multiplicative error models, following the strategy of constructing a bias-corrected weighted LS estimator in multiplicative error models proposed by Xu and Shimada [32] and Xu et al. [34], we can first eliminate this characteristic bias by directly deleting the third term on the left hand side of the normal equations (14a) and obtain the remaining part of the normal equations: ˆ −1 A )βˆ pb − AT Q ˆ −1 y = 0, (AT Q l l

(37)

where βˆ pb stands for the partially biased weighted LS estimate of β after removing the effect of multiplicative error nature. As a special stochastic model (6), if Qy = In and Qia = qi In with qi being a given scalar, then it is trivial to show that βˆ pb of Eq. (37) essentially becomes the ordinary LS estimate. From this point of view, βˆ pb may not be a good estimator of β.

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8797

On the other hand, since the matrix A is random, it will create extra bias of the weighted LS estimate, as is clear from [29]. Although one may go ahead to compute the bias due to the errors EA of A, as given in the formulae (36) and (37) of [29], we prefer to follow ˆ −1 A ) of Eq. (37) proposed by Xu [29] to the strategy of calibrating the normal matrix (AT Q l construct an N-calibrated bias-corrected weighted LS estimate of β, which is denoted by βˆ nc and given as follows: ˆ −1 A − WN )−1 AT Q ˆ −1 y, βˆ nc = (AT Q (38a) l

l

where the ijth element of the matrix WN is given by ˆ −1 Q ji )σ 2 . wi j = tr (Q N

l

(38b)

Here the matrix Qji has been defined in (13c) (also compare formula (41) of [29]). In particular, if the columns of the coefficient matrix A are supposed to be stochastically independent mutually, all the off-diagonal elements of WN are equal to zero, namely, wiNj = 0 for all i = j, and the diagonal elements of WN become ˆ −1 Qia ). wiiN = σ 2 tr (Q l

3.3. A simplified analysis of computational advantage over TLS Both the weighted LS estimator (11) or its iterative version (15) in the perspective of multiplicative error models and the weighted TLS solution to the EIV model (1) are essentially two different variants of the Gauss-Newton method applied to unconstrained nonlinear LS problems, as may be readily inferred from [10]. Thus, following Dennis and Schnabel [10], we can conclude that the iterative solution (15) is linearly convergent, as in the case of TLS (see [33]). The weighted LS solution (15) is significantly different from its TLS counterpart in the sense that our method does not need to adjust and update the measurements (or elements) of A. As a result, our method should be expected to be computationally more advantageous than TLS. More precisely, since both methods need to estimate the model parameters β, and since they are basically involved, more or less, numerically with the same matrix computation, the weighted TLS method is computationally more expensive due to the extra estimation to update the elements of A. To estimate the computational flops for this extra step of estimation, we assume a full variance-covariance matrix of A. Based on the weighted TLS solution of A in [33] and [37] (see also [22], but with a variance-covariance of special structure), under the implicit assumption that both the related inverse matrix and (Aβˆ tls − y) are readily available, we will need 2n2 flops of arithmetic operations to compute λ in formula (18) of [33] and 2n2 t2 flops of arithmetic operations to update the elements of A (compare formula (16a) in [33]), where βˆ tls stands for the weighted TLS solution of β. More importantly, since statistical quality measures have been essential for any statistical estimation, we will also have to compute, at the very least, the first order accuracy of the weighted LS solution (15) and the weighted TLS solution. Because the weighted LS solution (15) is only concerned about using the stochastic information on A to construct a suboptimal estimator of β without updating A, the cost of computing its first order variance-covariance matrix requires only O(t3 ) arithmetic operations. However, as is clear from [33], we need O((n + t )3 ) arithmetic operations to complete the first order accuracy of the weighted (nonlinear) TLS estimation for the EIV model (1). When n and t are sufficiently large, this computational complexity can become a very heavy burden. For the example of n = 10, 000, 000

8798

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

Table 1 ˆ TrueV — the true values of α and β; βˆ ls — the LS estimate; Parameter estimates and the standard deviations of β. βˆ — the weighted LS estimate from the perspective of multiplicative error models; βˆ nc — the N-calibrated biascorrected weighted LS estimate; βˆ tls — the TLS estimate; STD βˆ — the standard deviations of βˆ up to the first order approximation. Parameters

TrueV

βˆ ls

βˆ

βˆ nc

βˆ tls

STD βˆ

α β

2.20 1.00

1.012956 0.482063

2.111783 0.778811

2.172361 0.793592

2.111782 0.778811

0.238290 0.134371

and t = 1, 000 mentioned in the introduction, our weighted LS solution (15) and the weighted TLS method require O(109 ) and O(1021 ) arithmetic operations to compute the first order accuracy, respectively. 4. Examples This section will use examples to illustrate an application and the performances of the improved weighted LS estimate in the nonlinear normal equations (14a) and its first order accuracy, the N-calibrated bias-corrected weighted LS estimate (38a), and the weighted LS estimate of the variance of unit weight in (16). We will use part of the real climate data in Mann and Emanuel [36] to simulate the examples. One of the regression models used by Mann and Emanuel [36] is written as follows: T (t ) = αG(t ) + βS(t ) + (t ).

(39)

The main purpose of Mann and Emanuel [36] is to find an approximate relationship among the annual averaged tropical North Atlantic sea surface temperature (SST) T(t), the global mean SST G(t) and the annual tropospheric aerosol forcing/cooling S(t). The total number of data collected from 1870 to 1999 is 130. There are two unknown model parameters α and β in the model (39). (t) is the measurement error of T(t), which is assumed to be of mean zero and of the same variance. Although Mann and Emanuel [36] treated G(t) and S(t) as if they were free of errors, we will treat the model (39) as an EIV model, since both G(t) and S(t) are also measurements. We will follow [29] to simulate our examples and set the true values of α and β to 2.2 and 1.0, respectively. Since the standard deviation for the T(t) series was found to be equal to 0.19 by applying an ordinary LS fitting of the original data to the regression model (39), we will assume this value of standard deviation for all the three data series T(t), G(t) and S(t) in the following simulations. We further assume that the data series T(t), G(t) and S(t) are stochastically independent mutually and identically distributed. We can now use the true values of α and β and the standard deviation of 0.19 to generate the measurement series for T(t), G(t) and S(t). As a result, the model (39) becomes a standard EIV model for each set of simulated data series of T(t), G(t) and S(t). With the simulated example, we have computed the ordinary LS estimate βˆ ls of (8), the weighted LS estimate βˆ of Eq. (14a) and the N-calibrated bias-corrected weighted LS estimate βˆ nc of (38a). The results of parameter estimation are listed in Table 1, together with the true values and the weighted total least squares (TLS) estimates of the two parameters, which are denoted by βˆ tls (see e.g., [33]). The theoretical precisions (or standard deviations) of the first order approximation of βˆ are also listed in Table 1. Since the weights of the

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8799

measurements T(t) in the model (39) are assumed to be identical, the LS estimates βˆ ls of (8) in the third column of Table 1 are exactly the same as those of βˆ wls in (2). They are obviously biased from comparing the LS estimation results in the third column with the true values of parameters in the second column, which, however, can be theoretically expected from [29,30]. Although βˆ and βˆ nc are not statistically as rigorous as the TLS estimate βˆ tls , they perform almost as well as βˆ tls in this example, as can be numerically seen from Table 1. Actually, the differences among these three estimators are small and well within the stanˆ Nevertheless, since βˆ and βˆ nc require no estimation of the true values dard deviations of β. of the elements of A, they are computationally more advantageous than the TLS estimate βˆ tls . We have also computed the estimates of the variance of unit weight by using the three parameter estimators βˆ ls , βˆ and βˆ nc , respectively. The corresponding estimates of σ are equal to 0.354372, 0.189026 and 0.189082, respectively. When compared with the true value of 0.19, we can conclude that the LS estimate βˆ ls results in a significant bias, as can be theoretically expected by [15] (see also [6,29,30]). Both βˆ and βˆ nc have produced an almost unbiased estimate of σ 2 , as has been theoretically expected in this paper and [25]. 5. Conclusions Errors-in-variables models have been extensively investigated and applied to science and engineering problems, in particular, for the recent 40 years or so, thanks to the publication of Golub and van Loan [13]. Linear regression analysis or the weighted LS method can be most straightforwardly applied to deal with errors-in-variables (EIV) models, if the random errors in variables are treated as if they were deterministic or not random (see e.g., [6,15,29]). Although the method is computationally simple, it is well known to result in biased estimates of parameters and the variance of unit weight (see e.g., Hodges and Moore [15]; Davies and Hutton [6]). The effects of random errors in variables on variance component estimation have also been found to be significant by Xu [29,30]. The weighted TLS method is statistically rigorous and optimal. However, this method simultaneously estimates the parameters and the corrections to all the stochastically independent elements or functionals of A. In this digital and highly automatic era when data can be readily collected in an unprecedented large amount, the weighted TLS method can still be computationally too demanding. For example, let us assume n = 109 and t = 106 for the EIV model (1), which would roughly correspond to the problem of estimating an earth gravity model of 1000 degrees and orders with data of orbital positions, ranges and range-rates collected from a pair of tracking satellites at a sampling rate of 100 Hz for 15 days. Even with the fastest supercomputer of today, it is not reasonable to simultaneously estimate a total number of (1015 + 106 ) unknowns. This paper has strived to find an improved compromise between the weighted LS and TLS estimates of parameters in the EIV model (1), which would be expected to maintain the advantages of both the methods but keep away their disadvantages as far as possible. We have reformulated and constructed the weighted LS criterion from the perspective of multiplicative error models. The proposed weighted LS method is shown through the numerical example to maintain the advantages of computational simplicity of the weighted LS method and the optimality of the weighted TLS method but reduce the biasedness of the weighted LS method and completely avoid the computational burden of the weighted TLS method. We have carried out the analysis of statistical aspects of the weighted LS estimate from the perspective of

8800

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

multiplicative error models. We have provided an iterative formula to compute the weighted LS estimate and derived its accuracy of the first order approximation. Finally, we have proposed an almost unbiased estimate of the variance of unit weight and constructed an N-calibrated weighted LS estimate of parameters from the perspective of multiplicative error models. Acknowledgements The author thanks two anonymous reviewers for the constructive comments, which have helped improve the paper on the computational efficiency. This work is partially supported by the National Natural Science Foundation of China (41874012 and 41674013) through Prof. Y Shi. Appendix A We will now derive the four terms m1 , m2 , m3 and m4 required for the computation of biases of the weighted LS estimate βˆ in formula (34). We first rewrite the ith element of m1 of expression (33b) as follows: T −1 m1i = E {Eic Ql (I − Hl )l } −1 T = tr {Ql (I − Hl )E (l Eic )},

(40)

where m1i is the ith element of m1 , Eic is the ith column error vector of EA , and tr(·) stands for the trace of a square matrix. Inserting l of (19) into (40) yields T m1i = tr {Q−1 l (I − Hl )E (l Eic )} T = tr {Q−1 (I − Hl )E [( − EA β)Eic ]} ⎧l ⎡ ⎤⎫ t ⎨ ⎬  T⎦ ⎣ = −tr Q−1 (I − H ) E E β E l jc j ic ⎩ l ⎭ j=1

= −σ 2 tr {Q−1 l (I − Hl )

t 

β j Q ji },

(41)

j=1

since  is assumed to be stochastically independent of EA . To compute m2 of expression (33c), we only need to compute E(EA β ), which is denoted by m2b and given by, m2b = E (EA β )   −1 T = E EA Nl A Q−1 l l   −1 T = E EA Nl A Q−1 (  − E β) A l   −1 T −1 = −E EA Nl A Ql EA β = −Cm2 β, where

  −1 T Cm2 = E EA Nl A Q−1 l EA ,

(42)

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802

8801

whose ijth element is given by  −1 T   −1 T  ij 2 −1 cr cm2 = tr Nl A Q−1 l E (E jc Eir ) = σ tr Nl A Ql Q ji . Here Qcrji is the cofactor matrix between the jth column vector and the ith row vector of EA . Thus, m2 of (33c) becomes T

T

−1 m2 = −A Q−1 l E (EA β ) = A Ql Cm2 β.

(43)

Bearing (32c) in mind, we can compute the ith element of m3 , which is denoted by m3i and given by ⎧ ⎫ t ⎨ ⎬  −1 m3i = E Tl (I − Hl )T Q−1 β (Q + Q ) Q / 2(I − H )  j i j ji l l l l ⎩ ⎭ j=1 ⎧ ⎫ t ⎨ ⎬  −1 T = tr (I − Hl )T Q−1 β (Q + Q ) Q (I − H ) E (   ) /2 j i j ji l l l l l ⎩ ⎭ j=1 ⎧ ⎫ t ⎨ ⎬ = σ 2 tr β j (Qi j + Q ji )Q−1 (I − H ) /2. (44) l l ⎩ ⎭ j=1

Now we will compute the last term, namely, m4 of (33e), which, after slight re-arrangement, can be written as follows: T

t  t 

T

t  t 

T

t  t 

m4 = −A Q−1 l

  Qi j Q−1 l (I − Hl )E l (iβ β j + βi  jβ )

i=1 j=1

= −A Q−1 l

  Qi j Q−1 l (I − Hl )E l (iβ β j + βi  jβ )

i=1 j=1

= −A Q−1 l

 T Qi j Q−1 l (I − Hl )E l l vβ

i=1 j=1 T

= −σ 2 A Q−1 l

t  t 

Qi j (I − Hl )T vβ ,

(45)

i=1 j=1

where the vector vβ is equal to zero, except for the hth and kth elements, which take the values of β k and β h , respectively, with k = max (i, j) and h = min (i, j). References [1] R.J. Adcock, Note on the method of least squares, Analyst 4 (1877) 183–184. [2] R.J. Adcock, A problem in least squares, Analyst 5 (1878) 53–54. [3] A.R. Amiri-Simkooei, F. Zangeneh-Nejad, J. Asgari, On the covariance matrix of weighted total least-squares estimates, J. Surv. Eng. 142 (2016). Art 04015014. [4] R.J. Carroll, D. Ruppert, L.A. Stefanski, C.M. Crainiceanu, Measurement Error in Nonlinear Models – a Modern Perspective, 2nd edn, Chapman and Hall, London, 2006. [5] J.L. Coolidge, Two geometrical applications of the method of least squares, Am. Math. Mon. 20 (1913) 187–190. [6] R.B. Davies, B. Hutton, The effect of errors in the independent variables in linear regression, Biometrika 62 (1975) 383–391.

8802 [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

[27] [28] [29] [30] [31] [32] [33] [34]

[35] [36] [37]

P. Xu / Journal of the Franklin Institute 356 (2019) 8785–8802 W.E. Deming, The application of least squares, Phil. Mag. 11 (1931) 146–158. W.E. Deming, On the application of least squares — II, Phil. Mag. 17 (1934) 804–829. W.E. Deming, Statistical Adjustment of Data, Dover Publications, New York, 1964. J.E. Dennis Jr., R.B. Schnabel, Numerical methods for unconstrained optimization and nonlinear equations, in: SIAM Classics in Applied Mathematics, SIAM, Philadelphia, 1996. W.A. Fuller, Measurement Error Models, Wiley Interscience, New York, 1987. G.A. Gerhold, Least-squares adjustment of weighted data to a general linear equation, Am. J. Phys. 37 (1969) 156–161. G.H. Golub, C.F. van Loan, An analysis of the total least squares problem, SIAM J. Numer. Anal. 17 (1980) 883–893. C.C. Heyde, Quasi-Likelihood and Its Applications, Springer, New York, 1997. S.D. Hodges, P.G. Moore, Data uncertainties and least squares regression, Appl. Stat. 21 (1972) 185–195. J.R. Magnus, H. Neudecker, Matrix Differential Calculus with Applications in Statistics and Econometrics, Wiley, New York, 1988. I. Markovsky, S. van Huffel, On weighted structured total least squares, 2006, I. Lirkov, S. Margenov, J. Wa´sniewski, Proceedings of the LSSC 2005, LNCS, vol. 3743, 695–702. P. McCullagh, J. Nelder, Generalized Linear Models, (2nd edition), Chapman and Hall, London, 1989. Y. Nievergelt, Total least squares: state-of-the-art regression in numerical analysis, SIAM Rev. 36 (1994) 258–264. K. Pearson, On lines and planes of closest fit to systems of points in space, Phil Mag. 2 (1901) 559–572. B. Schaffrin, A note on constrained total least-squares estimation, Linear Algebra Appl. 417 (2006) 245–258. B. Schaffrin, A. Wieser, On weighted total least-squares adjustment for linear regression, J. Geod. 82 (2008) 415–421. B. Schaffrin, Y.A. Felus, An algorithmic approach to the total least-squares problem with linear and quadratic constraints, Stud. Geophys. Geod. 53 (2009) 1–16. G. Seber, C. Wild, Nonlinear Regression, John Wiley & Sons, New York, 1989. Y. Shi, P.L. Xu, Comparing the estimates of the variance of unit weight in multiplicative error models, Acta Geod. Geophys. 50 (2015) 353–363, doi:10.1007/s40328- 014- 0096- y. Y. Shi, P.L. Xu, J. Peng, C. Shi, J.N. Liu, Adjustment of measurements with multiplicative errors: error analysis, estimates of the variance of unit weight, and effect on volume estimation from Lidar-type digital elevation models, Sensors 14 (2014) 1249–1266, doi:10.3390/s140101249. S. van Huffel, J. Vandewalle, The Total Least Squares Problem: Computational Aspects and Analysis, SIAM, Philadelphia, 1991. R. Wedderburn, Quasi-likelihood functions, generalized linear models, and the Gauss–Newton method, Biometrika 61 (1974) 439–447. P.L. Xu, The effect of errors-in-variables on variance component estimation, J. Geod. 90 (2016) 681–701. P.L. Xu, Parameter estimation, variance components and statistical analysis in errors-in-variables models, in: W. Freeden (Ed.), Handbuch der Geodäsie, Springer, Berlin, 2018, doi:10.1007/978- 3- 662- 46900- 2_99- 1. P.L. Xu, Measurement-based perturbation theory and differential equation parameter estimation with applications to satellite gravimetry, Commun. Nonlinear Sci. Numer. Simul. 59 (2018) 515–543. P.L. Xu, S. Shimada, Least squares parameter estimation in multiplicative noise models, Commun. Stat. B29 (2000) 83–96. P.L. Xu, J.N. Liu, C. Shi, Total least squares adjustment in partial errors-in-variables models: algorithm and statistical analysis, J. Geod. 86 (2012) 661–675. P.L. Xu, J.N. Liu, Variance components in errors-in-variables models: estimability, stability and bias analysis. invited talk, in: Proceedings of the VIII Hotine-Marussi Symposium on Mathematical Geodesy, 2013. Rome, June 17–21, 2013. P.L. Xu, J.N. Liu, Variance components in errors-in-variables models: estimability, stability and bias analysis, J. Geod. 88 (2014) 719–734. M.E. Mann, K.A. Emanuel, Atlantic hurricane trends linked to climate change, EOS 87 (2006) 233–241. Y. Shi, P.L. Xu, J.N. Liu, C. Shi, Alternative formulae for parameter estimation in partial errors-in-variables models, J. Geod. 89 (2015) 13–16, doi:10.1007/s00190- 014- 0756- 2.