Modern Stochastics: Theory and Applications

Consider a problem of estimation of two lines by perturbed observations of points that lie on the lines. Let the true points $(\xi _{i},\eta _{i})$ lie on the union of two different lines $\eta =k_{1}\xi +h_{1}$ and $\eta =k_{2}\xi +h_{2}$, that is, Let these points be observed with perturbations $(\delta _{i},\varepsilon _{i})$, $i=1,\dots ,n$, that is, the observed points are $(x_{i},y_{i})$, $i=1,\dots ,n$, with The perturbations are assumed to be independent and identically normally distributed, where I is the $2\times 2$ identity matrix.

The parameters $k_{1}$, $h_{1}$, $k_{2}$, $h_{2}$, and ${\sigma }^{2}$ are to be estimated.

We consider both functional and structural models. In functional model, the true points are assumed to be nonrandom. In structural model, the true points are assumed to be independent and identically distributed (i.i.d.). The errors $(\delta _{i},\varepsilon _{i})$ are i.i.d. and independent of the true points.

In the structural model, $(\xi _{i},\eta _{i},\delta _{i},\varepsilon _{i})$ are i.i.d. random vectors, and thus, the observed points $(x_{i},y_{i})$ are i.i.d. In the functional model, the observed points are independent, Gaussian, with different means but with common covariance matrix.

Let the true points $(\xi _{i},\eta _{i})$ lie on the second-order algebraic curve Hereafter, a second-order algebraic curve is called a “conic section” or a “conic.”

The points are observed with Gaussian perturbations, and the perturbed points are denoted as $(x_{i},y_{i})$. We have the same equations as (2)–(4) in the two-line fitting model.

The vector of coefficients in (5) is denoted by $\boldsymbol{\beta }={(A,2B,C,2D,2E,F)}^{\top }$. The nonzero vector $\boldsymbol{\beta }$ and the error variance ${\sigma }^{2}$ are the parameters of interest.

Similarly to the two-line fitting model, the functional and the structural models are distinguished.

A couple of lines is a degenerate case of a conic section. Therefore, the conic section fitting model is an extension of the two-line fitting model.

We consider the adjusted least squares (ALS) estimator for unknown ${\sigma }^{2}$. The estimator is constructed in [7]. Introduce the $6\times 6$ symmetric matrix Asterisks are typed instead of some entries above the diagonal of a symmetric matrix. The entries of the matrix $\psi (x,y;v)$ are generalized Hermite polynomials in x and y. The matrix $\psi (x,y;v)$ is constructed such that $\operatorname{\mathsf{E}}\psi (x_{i},y_{i};{\sigma }^{2})=\psi (\xi _{i},\eta _{i};0)$ in the functional model and $\psi (\xi _{i},\eta _{i};0)\boldsymbol{\beta }=0$ for the true points and true parameters.

Denote

The estimator ${\hat{\sigma }}^{2}$ of the error variance ${\sigma }^{2}$ is obtained from the equation

Equation (6) always has a unique nonnegative solution. If $n\ge 6$, then the solution to (6) is positive almost surely.

The matrix $\boldsymbol{\varPsi }_{n}({\hat{\sigma }}^{2})$ is singular. Define the estimator $\hat{\boldsymbol{\beta }}$ of the vector $\boldsymbol{\beta }$ as a nonzero solution to the equation

The strong consistency of the ALS2 estimator is proved in [7] and [11] under somewhat different conditions. The asymptotic normality is proved in [12] for the functional model and in [13] for the structural model. Two consistent estimators of the asymptotic covariance matrix are constructed in [13].

Denote Under the conditions of Proposition 1 stated further, the latter limits exist almost surely. See [12] for explicit expressions of the matrices ${\psi ^{\prime }_{v}}(x,y;v)$, ${\boldsymbol{\varPsi }^{\prime\prime }_{1}}$, $\overline{\boldsymbol{\varPsi }}_{\infty }$, and ${\overline{\boldsymbol{\varPsi }}^{\prime }_{\infty }}$. Note that ${\boldsymbol{\varPsi }^{\prime\prime }_{1}}$ is a constant matrix.

“Eventually” in the previous statement means that almost surely there exists $n_{0}$ such that ${\hat{\boldsymbol{\beta }}}^{\top }{\boldsymbol{\varPsi }^{\prime }_{n}}({\hat{\sigma }}^{2})\hat{\boldsymbol{\beta }}<0$ for all $n\ge n_{0}$. In other words, almost surely, ${\hat{\boldsymbol{\beta }}}^{\top }{\boldsymbol{\varPsi }^{\prime }_{n}}({\hat{\sigma }}^{2})\hat{\boldsymbol{\beta }}\ge 0$ holds only for finitely many n.

Denote the normalized version of the true parameter

Normalize the estimator of $\boldsymbol{\beta }$ in such a way that ${\widetilde{\boldsymbol{\beta }}}^{\top }{\boldsymbol{\varPsi }^{\prime }_{n}}({\hat{\sigma }}^{2})\widetilde{\boldsymbol{\beta }}=-n$ and ${\boldsymbol{\beta }}^{\top }\widetilde{\boldsymbol{\beta }}\ge 0$. Therefore, denote

The methods of fitting an algebraic curve (or surface) to observed points can be classified as follows.

Algebraic distance methods, where the residuals in the equations for the algebraic curve are minimized. For example, the minimum point of the sum of squared residuals ${\sum _{i=1}^{n}}{(A{x_{i}^{2}}+2Bx_{i}y_{i}+C{y_{i}^{2}}+2Dx_{i}+2Ey_{i}+F)}^{2}$ (with some normalizing constraint in order to avoid $A=B=\dots =F=0$) in the conic fitting problem and ${\sum _{i=1}^{n}}{(k_{1}x_{i}+h_{1}-y_{i})}^{2}{(k_{2}x_{i}+h_{2}-y_{i})}^{2}$ in the two-line fitting problem is called the ordinary least squares (OLS) estimator.

The criterion function for the OLS estimator is simple enough and can be adjusted so that the resulting estimator is consistent (under some conditions). Such an estimator is called the adjusted least squares (ALS) estimator. The OLS and ALS estimators are method-of-moments estimators, meaning that the criterion functions for the estimators are polynomials whose coefficients are sample moments of coordinates of the observed points. Hence, the OLS and ALS estimators can be computed efficiently.

In order to obtain parameters of two lines, the observed points are fitted with a conic section, and then the parameters of the conic section are used to obtain the parameters of two lines. There are some papers where this idea is used.

The problem of estimating the fundamental matrix for two-camera view is considered in [6]. The fundamental matrix is a singular matrix whose left and right null-vectors are the coordinates of each camera in the coordinate system of the other camera. Initially, the ALS estimator of the fundamental matrix is evaluated. Then it is projected so that the estimated fundamental matrix is singular.

In [14] the problem of segmentation of a finite-dimensional vector space onto linear subspaces is considered, and the generalized principal component analysis method is introduced. The sample is fitted with an algebraic cone (a set of points that satisfy a homogeneous algebraic equation) by the OLS method. Then subspaces are extracted from the algebraic cone with use of a small learning sample. An application of segmentation of a vector space onto hyperplanes for searching planes on binocular image is given in [16].

In [15] an ellipsoid fitting problem with a constraint such that a center of the ellipsoid lies on a given line is considered. The algebraic distance with embedded constraint is minimized. The analytical (behavioral) properties of the optimization problem are studied. We consider a conic section fitting problem but with different constraint—the conic is degenerated to a couple of straight lines.

Geometric distance methods, where distances between the estimated curve and each point are minimized. The sum of squares of those distances is minimized, and the orthogonal regression (OR) estimator is obtained.

A numerical algorithm for evaluation of the orthogonal regression estimator is presented in monograph [1].

The orthogonal regression is consistent in the single straight line fitting problem [3, Section 1.3.2(a)]. In nonlinear models, the estimator may be inconsistent. There is a one-step correction procedure in explicit and implicit models [5, 9] with application in the ellipsoid fitting model [9]. However, in the two-line fitting model, the correction from [9] is unstable.

Probabilistic methods. They are used to obtain the maximum likelihood (ML) estimator and Bayes estimators.

Let $\{A_{n},\hspace{0.2778em}n=1,\hspace{0.1667em}2,\hspace{0.1667em}\dots \}$ be a sequence of random events. The random event $A_{n}$ is said to hold eventually if almost surely there exists $n_{0}$ such that $A_{n}$ occurs for all $n\ge n_{0}$. In other words, the random event $A_{n}$ holds eventually if and only if it does not occur only for finitely many n almost surely.

The estimator $\hat{\beta }$ is called asymptotically normal if $\sqrt{n}(\hat{\beta }-\beta _{\mathrm{true}})\to N(0,\varSigma )$ in distribution, were the asymptotic covariance matrix Σ may be singular, and n is the sample size. This definition differs from the conventional one adopted in asymptotic theory because here only $\sqrt{n}$-asymptotic normality is considered.

Let $\boldsymbol{\zeta }\sim N(\boldsymbol{\mu },\varSigma )$ be a bivariate random vector. Then $E=\{\boldsymbol{z}:{(\boldsymbol{z}-\boldsymbol{\mu })}^{\top }{\varSigma }^{-1}(\boldsymbol{z}-\boldsymbol{\mu })\le 1\}$ is called the 40% ellipsoid of the normal distribution because $\operatorname{\mathsf{P}}(\boldsymbol{\zeta }\in E)\approx 0.3935$. This is the ellipsoid where the probability density function is at least 0.3679 of its maximum.

In Section 2, we construct five estimators for parameters of the two line fitting model. In Section 3, we propose two definitions of the equivariance of an estimator and state that all of the five estimators are equivariant. The estimators are compared numerically in Section 4. The proofs are given in Appendix A.

The two-line fitting model is a restriction of the conic section fitting model. A couple of lines defined by the equation $(k_{1}\xi -\eta +h_{1})(k_{2}\xi -\eta +h_{2})=0$ is a degenerate conic section with coefficients with a constraint $C\ne 0$.

The conic section ALS2 estimator provides estimation of the error variance ${\sigma }^{2}$ and the coefficients $A,B,\dots ,F$.

Denote by $\nu (i)\in \{1,2\}$ the indicator of a line which the true point $(\xi _{i},\eta _{i})$ belongs to. Equation (1) can be rewritten as The indicator $\nu (i)$ is nonrandom in the functional model, and it is a random variable in the structural model.

There are two cases where the structural model is not identifiable. If the common distribution of the true points is concentrated on a straight line and on a single point (presumably not on the line), that is, then there are many ways to fit the true points with two lines. If the common distribution of the true points is concentrated in four points, that is, then there are three ways to fit the true points with two lines (unless three of the four points lie on a straight line, which is a particular case of (13)).

In order to estimate the parameters $k_{1}$, $h_{1}$, $k_{2}$, and $h_{2}$, we can solve Eqs. (12). With ignoring the last equation $F=Ch_{1}h_{2}$, the solution is

Substituting the elements of the ALS2 estimator $\hat{\boldsymbol{\beta }}={(\hat{A},2\hat{B},\hat{C},2\hat{D},2\hat{E},\hat{F})}^{\top }$ into the right-hand side of (15)–(17), we obtain an “ignore-$\widehat{F}$” estimator:

If the conic section estimated by the ALS2 estimator is a hyperbola, then the “ignore-$\widehat{F}$” estimate of the two lines comprises the asymptotes of the hyperbola.

Choose the sign ± in (18) such that $\hat{k}_{1}<\hat{k}_{2}$.

We need the notation for the function that expresses the line parameters $k_{1}$, $h_{1}$, $k_{2}$, $h_{2}$ in elements of $\boldsymbol{\beta }$ and is defined by (15)–(17). With this notation, we can write

Now, we state the asymptotic normality of the “ignore-$\widehat{F}$” estimator.

The matrix $K\varSigma _{\tilde{\beta }}{K}^{\top }$ is nonsingular.

Equation (11) represents a couple of intersecting straight lines if and only if Denote where the function $\Delta (\boldsymbol{\beta })$ and its derivative ${\Delta ^{\prime }}(\boldsymbol{\beta })$ are evaluated at the point $\boldsymbol{\beta }={(A,2B,C,2D,2E,F)}^{\top }$.

Perform one-step update of the estimator $\widetilde{\boldsymbol{\beta }}$ to make it closer to the surface $\Delta (\boldsymbol{\beta })=0$: Then use expressions (15)–(17) to estimate $k_{1},h_{1},k_{2},h_{2}$:

The sum of squared distances between each observed point and the closer of two lines is equal to The orthogonal regression estimator is a Borel-measurable function of observations such that

In the functional model, the orthogonal regression estimator is the maximum likelihood estimator. However, because the dimension of parameter space grows as the sample size is increasing, the orthogonal regression estimator may be inconsistent.

The estimator is constructed in the structural model, so it should be called the structural maximum likelihood estimator.

If a Gaussian distribution of a random point $(\xi ,\eta )$ is concentrated on a straight line $\eta =k\xi +h$, then it is a singular normal distribution: where $\mu _{\xi }$ and ${\sigma _{\xi }^{2}}$ are the expectation and variance of the random variable ξ. Note that the covariance matrix ${\sigma _{\xi }^{2}}\big(\begin{array}{c@{\hskip3.33pt}c}1& k\\{} k& {k}^{2}\end{array}\big)$ is singular and positive semidefinite.

If the distribution of a random point $(\xi _{i},\eta _{i})$ is concentrated on two straight lines $\eta =k_{1}\xi +h_{1}$ and $\eta =k_{2}\xi +h_{2}$ and the distribution on each line is Gaussian, then, due to (25), the conditional distributions are for $j=1,\hspace{0.1667em}2$. The matrices $\varSigma _{0j}$ are positive semidefinite and singular, that is, $\lambda _{\min }(\varSigma _{01})=\lambda _{\min }(\varSigma _{02})=0$, and the points $\boldsymbol{\mu }_{j}$ are the centers of Gaussian distribution of the points on each line.

The distribution of $(\xi _{i},\eta _{i})$ is a mixture of two singular normal distributions where $p=\operatorname{\mathsf{P}}(\nu (i)=1)=\operatorname{\mathsf{P}}(\nu (1)=1)$ is the probability that the point $(\xi _{i},\eta _{i})$ lies of the first line.

The distribution of the observed points is also a mixture of two Gaussian distributions with $\varSigma _{j}=\varSigma _{0j}+{\sigma }^{2}\mathbf{I}$, where ${\sigma }^{2}$ is the error variance; see (4). Note that $\lambda _{\min }(\varSigma _{1})=\lambda _{\min }(\varSigma _{2})={\sigma }^{2}$.

The likelihood function for the sample of points with a mixture of two normal distributions is where is the density of a bivariate normal distribution.

One method of evaluating the maximum likelihood estimator is as follows:

The denominator $\hat{\sigma }_{jxx}-{\hat{\sigma }}^{2}$ may be equal to 0 with some positive probability. Occurrence of this event means that the estimated figure is a straight line and a single point outside the line rather than two straight lines.

In order to make the statement of consistency easier, assume that $k_{1}<k_{2}$ and choose the estimator such that $\hat{k}_{1}\le \hat{k}_{2}$.

The regular best asymptotically normal (RBAN) estimators were developed by Chiang [4]. Our RBAN moment estimator differs from the original RBAN so that not only the observed points $(x_{i},y_{i})$, but also monomials ${x_{i}^{p}}{y_{i}^{q}}$, $p+q\le 4$, are averaged.

Introduce the 14-dimensional vectors whose elements are the monomials of coordinates of observed points:

Evaluate the average and sample covariance matrix of the vectors $\boldsymbol{m}_{i}$:

Denote where ξ, δ, and ε are independent random variables such that Basically, the function $f_{1}$ is defined for all $\mu _{p}$, $p=1,\dots ,4$, that comprise possible 4-tuples of moments of a random variable, that is, satisfy see [10]. However, since the elements of the vector-function $f_{1}$ are polynomials of its arguments, it can be extended to ${\mathbb{R}}^{7}$.

Denote

In the structural model,

Consider the equation It is a system of 14 equations in 14 variables. If (33) has a solution, then the moment estimator can be defined as one of the solutions. However, (33) may have no solution.

In the rest of Section 2.4, ${\mu _{\bullet }^{(\bullet )}}={({\mu _{q}^{(j)}})_{j=1,}^{2}}{_{q=1}^{4}}$ is a $2\times 4$ matrix.

The estimator is defined as a point where ${(f_{2}(\dots )-\overline{\boldsymbol{m}})}^{\top }{\varSigma _{m}^{-1}}(f_{2}(\dots )-\overline{\boldsymbol{m}})$ attains its minimum: This minimization problem is similar to that in Theorem 6 in [4]. The minimum can be evaluated explicitly, and this allows us to reduce the dimension of minimization problem. The reduction of dimension of the optimization problem was used, for example, in [8].

The routines evaluating the RBAN-moment estimator and the estimator for its covariance matrix are developed without rigid theoretical basis; see Section 4.3.

The similarity transformation of ${\mathbb{R}}^{2}$ is where U is an orthogonal matrix, $K\ne 0$ is a scaling coefficient, and $\Delta \boldsymbol{z}\in {\mathbb{R}}^{2}$ is an intercept.

The transformation of a sample of points acting elementwise is also denoted $g(Z)$: if $Z=\{\boldsymbol{z}_{i},\hspace{0.2778em}i=1,\dots ,n\}$, then $g(Z)=\{g(\boldsymbol{z}_{i}),\hspace{0.2778em}i=1,\dots ,n\}$.

Hereafter, we use vector notation: the observed points are denoted $\boldsymbol{z}_{i}={(x_{i},y_{i})}^{\top }$, and the true points are denoted $\boldsymbol{\zeta }_{i}={(\xi _{i},\eta _{i})}^{\top }$.

The underlying statistical structure is $({\mathbb{R}}^{n\times 2},\mathcal{B}({\mathbb{R}}^{n\times 2}),P_{Z\mid \theta },\theta \in \varTheta )$, where $Z\in {\mathbb{R}}^{n\times 2}$ is the observed sample, $Z=\{\boldsymbol{z}_{i},\hspace{0.2778em}i=1,\dots ,n\}$, $\mathcal{B}({\mathbb{R}}^{n\times 2})$ is the Borel σ-field, and θ is a parameter that uniquely identifies the distribution of the observed points; $\theta =(\boldsymbol{\zeta }_{1},\dots ,\boldsymbol{\zeta }_{n};{\sigma }^{2})$ in the functional model, and $\theta =(P_{\zeta };{\sigma }^{2})$ in the structural model. Here $\boldsymbol{\zeta }_{1},\dots ,\boldsymbol{\zeta }_{n}$ are points located on two strait lines, and $P_{\zeta }$ is a probability measure concentrated on two straight lines.

The statistical structure is invariant with respect to transformation g if the change of the probability measure induced by the transformation of the sample can be obtained by some transformation $\tilde{g}$ of parameters, that is, if there exists a bijection $\tilde{g}:\varTheta \to \varTheta $ such that Here $P_{g(Z)\mid \theta }$ is the induced probability measure; it is sometimes denoted $P_{g(Z)\mid \theta }=P_{Z\mid \theta }{g}^{-1}$.

The statistical structure is similarity invariant if it is invariant with respect to all similarity transformations of the form (34).

In order to become similarity invariant, the underlying statistical structure needs some extension. We assume that the true points lie on two lines, which may be parallel to the y-axis. The following restrictions do not ruin the invariance:

With these restrictions, the statistical structure is invariant with in the functional model and in the structural model.

Let $\ell \ell (\theta )$ and $\operatorname{sigma2}(\theta )$ be functions that extract the parameters of interest. If, in the functional model, $\theta =(\boldsymbol{\zeta }_{1},\dots ,\boldsymbol{\zeta }_{n};{\sigma }^{2})$ and points $\boldsymbol{\zeta }_{1},\dots ,\boldsymbol{\zeta }_{n}$ lie on the lines $\ell _{1}$ and $\ell _{2}$ or if, in the structural model, $\theta =(P_{\zeta };{\sigma }^{2})$ and the probability measure is concentrated on the union of two lines $\ell _{1}\cup \ell _{2}$, then $\ell \ell (\theta )=\{\ell _{1},\ell _{2}\}$. If $\theta =(\dots ;{\sigma }^{2})$, then $\operatorname{sigma2}(\theta )={\sigma }^{2}$.

We treat $\{\ell _{1},\ell _{2}\}$ as an unordered couple, that is, $\{\ell _{1},\ell _{2}\}=\{\ell _{2},\ell _{1}\}$.

The transformation of the lines parameters and the transformation of ${\sigma }^{2}$ do not interfere each other, and these transformations are not interfered by a particular location or distribution of true points on the lines, that is, the parameter transformation $\tilde{g}$ is such that there exists transformations $\tilde{g}_{\ell \ell }$ and $\tilde{g}_{{\sigma }^{2}}$ such that These transformations are

The estimator is called equivariant with respect to the transformation g if, when the data are transformed, the estimator follows the inducing transformation of parameters. The estimator $\widehat{\ell \ell }(Z)$ for two lines and the estimator ${\hat{\sigma }}^{2}(Z)$ for error variance are equivariant with respect to similarity transformation g if The estimator is called similarity equivariant if it is equivariant with respect to any similarity transformation g.

In a fitting problem, an estimator for a “true figure” is called fitting equivariant with respect to transformation $g(\boldsymbol{z})$, $g:\mathbb{R}\to \mathbb{R}$ if, when the sample is transformed, the estimated “true figure” follows the same transformation g. An estimator is called similarity fitting equivariant if it is fitting equivariant with respect to any similarity transformation.

In the two-line fitting problem, denote by $\cup \{\ell _{1},\ell _{2}\}=\ell _{1}\cup \ell _{2}$ the union of a pair of two lines. An estimator $\widehat{\ell \ell }(Z)$ is similarity fitting equivariant if and only if for any similarity transformation $g(\boldsymbol{z})$,

The similarity fitting equivariant estimator depends on geometry of the plane and does not depend on the Cartesian coordinate system used.

Because of (35), in the two-line fitting model, the estimator for two lines $\widehat{\ell \ell }(Z)$ is similarity equivariant if and only if it is similarity fitting equivariant.

Some troubles, which may arise during estimation, are not addressed yet.

In order to define the equivariance of an unreliable estimator, we allow that the estimators fail simultaneously in both sides of (36), (37), or (38). Also, we allow that for fixed similarity transformation $g(\boldsymbol{z})$, equation (36), (37), or (38) does not hold with probability 0.

The equivariance of the ALS2 estimator in the conic section fitting problem is verified in [11, Section 5.5] (see Theorem 30 there for similarity fitting equivariance). That implies the equivariance of the “ignore-$\widehat{F}$” estimator.

In order to make the updated before ignore-$\widehat{F}$ step estimator equivariant, we use normalization of the ALS2 estimator (9) rather than $\| \hat{\boldsymbol{\beta }}\| =1$.

The orthogonal regression estimator and the parametric maximum likelihood estimator are maximum likelihood estimators, but in different models. Thus, they are equivariant.

The criterion function for the RBAN-moment estimator is similarity invariant. This means that the criterion function does not change when the data sample follows a similarity transformation and the parameters follow the inducing transformation. Thus, the RBAN-moment estimator is equivariant.

Consider a further restriction of the mixture-of-two-normal-distributions model from Section 2.3. Assume that covariance matrices of $\varSigma _{1}$ and $\varSigma _{2}$ have the same diagonal entries but additive inverse off-diagonal entries:

The statistical structure is invariant in scaling of the y-coordinate, $(x_{\mathrm{new}},y_{\mathrm{new}})=(x_{\mathrm{old}},\hspace{0.2222em}ry_{\mathrm{old}})$, $r>0$. This transformation maps the lines $y=-kx+h_{1}$ and $y=kx+h_{2}$ onto the lines $y=-rkx+rh_{1}$ and $y=rkx+rh_{2}$, respectively. The maximum likelihood estimator in this model is equivariant. However, this equivariance is somewhat strange. The transformation of parameters that induces the scaling of the y-coordinate of the observed points does not induce the same transformation of the true points nor the same mapping of the true lines. The estimated lines follow the transformation of parameters rather than the transformation of observed points. This is illustrated in Fig. 1.

Let $k_{\mathrm{old}}$ be the true value of the parameter k before the transformation. Then after the transformation, the value of the parameter is with $t=({r}^{2}{k_{\mathrm{old}}^{2}}-1)\hspace{0.1667em}{\sigma _{\xi \hspace{0.1667em}\mathrm{old}}^{2}}+({r}^{2}-1)\hspace{0.1667em}{\sigma _{\mathrm{old}}^{2}}$. If $0<{r}^{2}\ne 1$, $k_{\mathrm{old}}\ne 0$, and ${\sigma }^{2}>0$, then $k_{\mathrm{new}}\ne rk_{\mathrm{old}}$. Hence, the maximum likelihood estimator is not fitting equivariant with respect to scaling of the y-coordinate here.

A sample of the true points $(\xi _{i},\eta _{i})$, $i=1,\dots ,n$, is generated from a random distribution concentrated on (a subset of) two lines. Three distributions of the true points are used; see Fig. 2:

These three distributions of true points are concentrated on the same two lines which intersect one another at the point $(-0.08,0.31)$.

For the same sample of true points $\{(\xi _{i},\eta _{i}),\hspace{0.2778em}i=1,\dots ,n\}$, 100 samples of the measurement errors $\{(\delta _{i},\varepsilon _{i}),\hspace{0.2778em}i=1,\dots ,n\}$, ${(\delta _{i},\varepsilon _{i})}^{\top }\sim N(0,{\sigma }^{2})$, are simulated, and 100 samples of the observed points $(x_{i},y_{i})$ are obtained; see (2) and (3). For each sample of the observed points, the estimates of the parameters of the true lines were evaluated with the following five methods: two ALS2-based estimators (the ignore-$\widehat{F}$ estimator and the estimator with one-step update of the ALS2 estimator before the ignore-$\widehat{F}$ step), the orthogonal regression estimator, the parametric maximum likelihood estimator, and the RBAN moment estimator.

For each estimated couple of lines, the point of their intersection is found. The 100 estimates of intersection points are averaged, and their sample standard deviations are evaluated. For the ALS2-based estimators and the RBAN moment estimator, the standard errors of the estimators are also evaluated.

For computation of the orthogonal regression estimator, the k-means method is used. Initially, two lines were chosen randomly. Then classification and mean steps are alternated. On the classification step, the observed points are split into two clusters based on which line is closer to the point. (The first cluster contains all the observed points that are closer to the first line than to the second line, and the second cluster contains the other observed points.) On the means step, each cluster is fitted with a straight line by the orthogonal regression method (the two lines are updated). The algorithm is completed when the classification step does now change the clusters. The obtained parameters of the two lines deliver a local minimum to the criterion function $Q(k_{1},h_{1},k_{2},h_{2})$ (24). Trying to obtain the global minimum, the algorithm is restarted several times with different initial two lines.

For computation of the parametric maximum likelihood estimator, the expectation–maximization algorithm [2] is used. The equation for optimization problem of finding a minimum of the likelihood function $L(p,\boldsymbol{\mu }_{1},\varSigma _{1},\boldsymbol{\mu }_{2},\varSigma _{2})$ (28) such that $\lambda _{\min }(\varSigma _{1})=\lambda _{\min }(\varSigma _{2})$ is where $\textit{par}$ is a vector parameterization of $(\boldsymbol{\mu }_{1},\varSigma _{1},\boldsymbol{\mu }_{2},\varSigma _{2})$ such that $\lambda _{\min }(\varSigma _{1})=\lambda _{\min }(\varSigma _{2})$. Hence, the maximum likelihood estimator is a stationary point of the function with fixed

The EM algorithm is iterative. Once the $(m-1)$th approximation $({p}^{(m-1)},{\mu _{1}^{(m-1)}},{\varSigma _{1}^{(m-1)}},{\mu _{2}^{(m-1)}},{\varSigma _{2}^{(m-1)}})$ is obtained, the weights are evaluated: Then mth approximation $({p}^{(m)},{\mu _{1}^{(m)}},{\varSigma _{1}^{(m)}},{\mu _{2}^{(m)}},{\varSigma _{2}^{(m)}})$ is obtained by minimizing under the constraint $\lambda _{\min }(\varSigma _{1})=\lambda _{\min }(\varSigma _{2})$. The point where the minimum is attained can be explicitly expressed in ${w_{i}^{(m-1)}}$, $x_{i}$, and $y_{i}$, $i=1,\dots ,n$.

In case the criterion function has multiple minima, a consistent estimator—that is, the “ignore-$\widehat{F}$” estimator—is used as the initial point, and the criterion function $Q(\theta )$ is searched for a local minimum nearby. Here $\theta ={(k_{1},h_{1},k_{2},h_{2},{\sigma }^{2})}^{\top }$ is a vector meaning the parameters of interest.

The knowledge or misspecification of the parameter p does not affect the estimator for the parameters of interest $k_{1}$, …, ${\sigma }^{2}$. Thus, for estimation of the asymptotic covariance matrix, assume $p=0.5$ to be known. The estimator of the asymptotic covariance matrix of $(\theta ,\mathbf{M})$ is where The estimator of the asymptotic covariance matrix of θ is the principal submatrix of $\widehat{\varSigma }_{\theta ,\mathbf{M}}$.

Average of estimated centers over 100 simulations, standard deviations over 100 simulations, and medians of estimated standard errors are presented in Tables 1–3.

Using the estimator $\widehat{F}$ (by one-step update before ignore-$\widehat{F}$ step), we improve the precision of estimation. The precision of the RBAN-moment estimator approximates the precision of the updated before ignore-$\widehat{F}$ step estimator.

The parametric maximum likelihood estimator is the best when the normality condition, which was assumed during construction of the estimator, is satisfied. Otherwise, it is biased.

The orthogonal regression and the maximum likelihood estimators are good for small error variance (${\sigma }^{2}={0.02}^{2}$). For ${\sigma }^{2}=0.{1}^{2}$, the orthogonal regression estimator is broken down when the distribution of true points is a mixture of two normal distributions and is biased for the two other distributions of true points.

Mean-square deviance of the intersection of the estimated lines from the true intersection point is presented in Table 4.

For small errors, the RBAN-moment estimator is a bit less accurate than the updated before ignore-$\widehat{F}$ step estimator. For ${\sigma }^{2}=0.{1}^{2}$, the difference is negligible.

The parametric maximum likelihood estimator has the smallest deviation from the true value, except for the discrete distribution of true points and ${\sigma }^{2}=0.01$.

For small errors (${\sigma }^{2}={0.02}^{2}$), the orthogonal regression estimator outperforms the consistent estimators and has the deviation approximately as small as the parametric maximum likelihood estimator.

Normalization of the estimator of $\boldsymbol{\beta }$ affects the ALS2-based estimator of two lines with one-step update before the ignore-$\widehat{F}$ step. With normalization $\| \hat{\boldsymbol{\beta }}\| =1$, the derived estimator of two lines is not equivariant, whereas with normalization ${\tilde{\boldsymbol{\beta }}}^{\top }{\boldsymbol{\varPsi }^{\prime }_{n}}({\hat{\sigma }}^{2})\tilde{\boldsymbol{\beta }}=-n$, the derived estimator is equivariant. Comparison of equivariant and nonequivariant versions of the estimator is displayed in Table 5.

There is a tendency that the equivariant version of the estimator is more accurate for small samples than the nonequivariant version. The two versions of the estimator are consistent and asymptotically equivalent. When the estimation is precise, the difference between the versions is negligible. When the estimation is imprecise, it is impossible to make inference which version is more accurate.

In [13] a conic section fitting model is considered, and two estimators ($\widehat{\varSigma }_{\mathrm{true}}$ and $\widehat{\varSigma }_{\mathrm{sample}}$) for the asymptotic covariance matrix of the ALS2 estimator are constructed.

The software developed here can be used to make numerical comparison of the estimates of the asymptotic covariance matrices. The data are generated as described in Section 4.1 with 1000 simulations for each set of true points. Thus, the true conic unnecessarily was chosen degenerate. For each simulation, the parameters of the conic section were estimated; its center is found, and two confidence ellipsoids for the center were constructed using two different estimators of the asymptotic covariance matrix.

The sample coverage probability and median (over 1000 ellipsoids) area of the confidence ellipsoids is presented in Table 6. The ellipsoids were constructed for confidence levels 0.8 and 0.95. The area of 95% confidence ellipsoids is displayed in Table 6, and the area of 80% confidence ellipsoid is $\log _{20}(5)=0.5372$ of the area of 95% confidence ellipsoids.

Note that standard errors for coverage probability are $1.3\% $ for 80% confidence ellipsoids and $0.7\% $ for 95% confidence ellipsoids. The simulations do not allow us to make an inference which estimator is better. Thus, $\widehat{\varSigma }_{\mathrm{sample}}$-based estimator updated before ignore-$\widehat{F}$ step is compared with other estimators in simulations in Section 4.4 because of simpler explicit expression for $\widehat{\varSigma }_{\mathrm{sample}}$.