A maximum likelihood estimator has been applied to find regression parameters of a straight line in case of different error models. Assuming Gaussian-type noise for the measurement errors, explicit results for the parameters can be given employing Mathematica. In the case of the ordinary least squares (), total least squares (), and least geometric mean deviation () approaches, as well as the error model of combining ordinary least squares ( and ) in the Pareto sense, simple formulas are given to compute the parameters via a reduced Gröbner basis. Numerical examples illustrate the methods, and the results are checked via direct global minimization of the residuals.

Introduction

Generally, to carry out a regression procedure one needs to have a model , an error definition , and the probability density function of the error . Considering the set as measurement points, the maximum likelihood approach aims at finding the parameter vector that maximizes the likelihood of the joint error distribution. Assuming that the measurement errors are independent, we should maximize (see eg. [1])

(1)

Instead of maximizing this objective, we minimize

(2)

Consider the Gaussian-type error distribution as ; then our estimator is

(3)

In our case the model is a line,

(4)

It can be seen that (in the case of Gaussian-type measurement noise) only the type of the error model determines the parameter values, since we should always minimize the least squares of the errors. There are different error models, which can be applied to fitting a line in a least-squares sense. The error model frequently employed, assuming an error-free independent variable , is the ordinary least squares model ()

(5)

Similarly, one may also consider an error-free dependent variable . Then the error model () is

(6)

These approaches are called the algebraic approach.

Another error model considers the geometrical distance between the data point and the line to be fitted. This type of fitting is also known as orthogonal regression, since the distances of the sample points from the line are evaluated by computing the orthogonal projection of the measurements on the line itself. The error in this case [2] is

(7)

This geometrical approach or total least squares () approach can also be considered as an optimization problem with constraints; namely, one should minimize the errors in both variables [3]:

(8)

under the conditions

(9)

In addition, one can also combine and to construct an error model. The first possibility is to consider the geometric mean of these two types of errors,

(10)

These error models are illustrated in Figure 1.

Figure 1. The different error models in the case of fitting a straight line.

This model is also called the least geometric mean deviation approach or model (see [4)]. As a second possibility, one may consider and as competing functions of the parameters and find their Pareto-front representing a set of optimal solutions for the parameters . Since this multi-objective problem is convex, the objective can be expressed as a linear combination of these error functions, namely

(11)

where is a parameter, , and the set of optimal solutions of the parameters belonging to the Pareto-front is . You can choose the value of depending on your trade-off preference between and [5].

Application of Symbolic Computation

SuperLog Function

Symbolic computation can be used to avoid direct minimization and to get an explicit formula for the estimated parameters. We apply the Mathematica function SuperLog developed in [6], which uses pattern matching that enhances Mathematica‘s ability to simplify expressions involving the natural logarithm of a product of algebraic terms.

Let us activate this function.

Then this is the ML estimator for Gaussian-type noise.

Ordinary Least Squares ()

Now let us consider the problem.

Here are the necessary conditions for the optimum.

Let us introduce the following constants:

(12)

(13)

(14)

(15)

(16)

In those terms, here are the necessary conditions for the optimum.

Then this is the optimal solution of the parameters.

Total Least Squares ()

Although the equation system for the parameters of is linear, for other error models we get a multivariable algebraic system. Now consider the problem. Here is the maximum likelihood function.

Therefore here is the equation system to be solved.

Since , the conditions are as follows.

A Gröbner basis solves this system, eliminating .

Since the second equation is linear, it is reasonable to compute first, then .

Least Geometric Mean Deviation ()

The error model also leads to a second-order polynomial equation system. Now here is the ML estimator.

Consequently, here is the system to be solved for the parameters.

Assume .

Again a Gröbner basis gives a second-order system.

When is known, the other parameter can be computed.

Pareto Approach

In the case of the Pareto approach, the system is already fourth order.

Here is the system.

Here is the system in compact form.

Here is the Gröbner basis for the first parameter.

Assume that .

After solving this polynomial for , the other parameter can be solved from the second equation, which is linear in .

Numerical Example

Data Samples

Consider some data on rainfall (in mm) and the resulting groundwater level changes (in cm) from a landslide along the Ohio River Valley near Cincinnati, Ohio [7].

There are 14 measurements.

This displays the measured data.

Figure 2. The measured data: rainfall versus water level change in dimensional form.

Computation of the Constants for the Equation Systems

The constants , , , , and in equations (12) to (16) are needed.

This separates the data.

This transforms the data into dimensionless form.

Figure 3. The measured data: rainfall versus water level change in dimensionless form.

Now the constants can be computed.

Computation of the Parameters of the Fitted Line

Model

Here are the estimated parameters employing the explicit solutions.

This checks the result.

Figure 4 shows the estimated line with the sample points.

Figure 4. The sample points with the line estimated with .

Model

Here are the first and second parameters.

Here is a check of this result on the basis of the definition. Equation (8) gives the objective function.

The constraints are .

The unknown variables are not only the parameters, but the adjustments as well.

This uses a built-in global optimization method. (This takes a long time to compute.)

The estimation gives a result quite different from the model; see Figure 5.

Figure 5. The lines estimated with the (red) and (green) models.

Since the constraints are linear, the optimization can be written in unconstrained form, reducing the original number of variables to .

Model

Now here is the first parameter.

This uses the result.

Here is a numerical check of the objective.

Figure 6 shows this result together with the and models.

Figure 6. The lines estimated with the (red), (green), and (blue) models.

Pareto Approach

The first parameter is a fourth-order polynomial.

The best trade-off between and is to let .

This is the real positive solution.

Using this value gives the second parameter.

We compute the solution using direct global minimization. Here is the objective.

This gives the result.

Figure 7 shows this solution with the results of the other models.

Figure 7. The lines estimated with the (red), (green), and (blue) models, and the Pareto approach with (magenta).

Conclusion

The numerical computations show that the formulas developed by an ML estimator via symbolic computation to determine the parameters of a straight line to be fitted provide correct results and require considerably less computation time than the direct methods based on global minimization of the residuals. Our examples also illustrate that the , , and Pareto approaches give more realistic solutions than the traditional , since Figure 7 shows there are at least two outliers in the sample set.

References

[1]	W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, 2nd ed., Cambridge: Cambridge University Press, 1992.
[2]	M. Zuliani. “RANSAC for Dummies.” (Jan 10, 2014) vision.ece.ucsb.edu/~zuliani/Research/RANSAC/docs/RANSAC4Dummies.pdf.
[3]	B. Schaffrin, “A Note on Constrained Total Least-Squares Estimation,” Linear Algebra and Its Applications, 417(1), 2006 pp. 245-258. doi:10.1016/j.laa.2006.03.044.
[4]	C. Tofallis, “Model Fitting for Multiple Variables by Minimising the Geometric Mean Deviation,” in Total Least Squares and Errors-in-Variables Modeling: Analysis, Algorithms and Applications (S. Van Huffel and P. Lemmerling, eds.), Dordrecht: Kluwer, 2002.
[5]	B. Paláncz and J. L. Awange, “Application of Pareto Optimality to Linear Models with Errors-in-All-Variables,” Journal of Geodesy, 86(7), 2012 pp. 531-545. doi:10.1007/s00190-011-0536-1.
[6]	C. Rose and M. D. Smith, “Symbolic Maximum Likelihood Estimation with Mathematica,” The Statistician, 49(2), 2000 pp. 229-240. www.jstor.org/stable/2680972.
[7]	W. C. Haneberg, Computational Geosciences with Mathematica, Berlin: Springer, 2004.
B. Paláncz, “Fitting Data with Different Error Models,” The Mathematica Journal, 2014. dx.doi.org/doi:10.3888/tmj.16-4.

About the Author

Béla Paláncz received his D.Sc. degree in 1993 from the Hungarian Academy of Sciences and has wide-ranging experience in teaching and research (RWTH Aachen, Imperial College London, DLR Köln, and Wolfram Research). His main research fields are mathematical modeling and symbolic-numeric computation.

Béla Paláncz
Department of Photogrammetry and Geoinformatics,
Budapest University of Technology and Economics
1521 Budapest, Hungary
palancz@epito.bme.hu