The LeastTrimmedSquares function computes a robust linear estimator from a collection of input data. It is a quantity first described by Rousseeuw [1], and is computed using the FAST-LTS algorithm by Rousseeuw and Van Driessen [2]. The article [1] that introduced least trimmed squares regression also introduced least median of squares regression, but least trimmed squares regression was later considered the better choice by many, including the original author - see e.g. [2].

Least trimmed squares regression yields a robust estimator. This means that it will continue to perform well if some points are replaced by outliers. Least-squares regression, the type of regression most commonly used and implemented by LinearFit and NonlinearFit, is very susceptible to outliers.

Conceptually, one performs least trimmed squares regression by selecting an arbitrary subset of h of the input data points, where h is an input parameter (which can be set using the include option). One then performs regular least-squares regression on these h data points, and records the sum of squared residuals. This sum of squared residuals is now minimized over all subsets of h data points.

The FAST-LTS algorithm implemented in this command is an approximation algorithm: it will in all reasonable cases return the optimal regression model, but this is not guaranteed. An exact algorithm that does guarantee the optimal result can be selected with the subsets = all option, but this algorithm will take a number of steps that grows very quickly with the size of the problem.

Least trimmed squares regression with parameter h on n input points has a breakdown point of h / n. This means that when fewer than n - h of the data points in a sample are changed, the result of least trimmed squares regression can only change a limited amount. See the discussion of the include option below for limits on the values h can assume.

In the first and third calling sequences, the parameter X is a Matrix containing the values of the independent variables; each row corresponds to one data point and each column corresponds to one independent variable. Alternatively, you can specify X as a list of lists. If there is only a single independent variable, you can also specify X as a column Vector or a list.

In those same calling sequences, Y is a Vector containing the values of the dependent variable, in the same order as the data points occur in X. Alternatively, you can specify Y as a list, or a Matrix with a single column.

In the second and fourth calling sequences, XY is a Matrix containing both the dependent and independent data, with the dependent data being in the right-most column.

The returned value is an expression representing the estimator evaluated at the value v. By supplying a variable name here, you obtain the general expression for the estimator in terms of variables v[1], v[2], etc., corresponding to the independent variables. In particular, if there is a single independent variable, the result will be in terms of v[1]. If instead you supply a list (or Vector) of values, corresponding to the values of the independent variables, you obtain the value of the estimator at those values.

include = h

This option determines how many data points are included (and consequently, how many are ignored). The value h should either be an integer greater than 1, in which case it is directly used as the number of data points to include, or a number that is less than or equal to one, in which case it is multiplied by n (the number of data points) and rounded to obtain the number of data points to include.

In all cases, the resulting number of data points to include (which will be denoted by h in the rest of this discussion) needs to be at least equal to (n + p + 1)/2, where n is the number of data points and p is the number of independent variables, including (if present) the intercept. The default value for h is this minimum number (rounded up to an integer).

intercept = tf

This option determines if a column for the intercept (a constant term in the resulting regression model) should be added to X. The default value is true; you can supply the option intercept = false to turn this off. The intercept is internally handled by adding a column filled with the value 1 to the right of X.

Using this option is necessary if there is already an independent variable assuming only one constant value, or if the all-one vector is otherwise already in the column space of X: the matrix used internally needs to have full column rank, so it will not work if the all-one vector is joined to a Matrix that already contains that vector in its column space.

subsets = s

This option can be used by experts to control the details of how the approximation process used for the FAST-LTS algorithm works. The discussion of this option assumes knowledge of the Rousseeuw and Van Driessen paper [2] that describes the algorithm.

If p = 1 (after including the intercept), this command uses a fast exact algorithm and ignores the subsets option.

By submitting the subsets = all option, the user can select a variant that tries all p-subsets of the data points. It runs two C-steps on all regression lines generated this way, selects the 10 best results, and runs C-steps on these results until convergence. This should in all cases return the optimal result, but the number of regression lines to process grows very quickly as a function of n and p. This algorithm is selected by default if binomial(n, p) < 1500.

Otherwise, the algorithm runs in a number of rounds in which candidate regression lines are generated and subjected to C-steps; the details of these depend by default on n but can be fine-tuned by submitting an option of the form subsets = s, where s is a list of lists. Each sublist corresponds to one round of the algorithm, and contains a number of options specified as equations, listed below.

subsets = c and size = s

These options specify that the data points should be subdivided into c disjoint subsets of size s each, all treated independently. Any points not included are ignored for the time being. The other options then apply to each subset; for example, the generate option specifies how many candidates are generated for each subset. The default values are c = 1 and s equal to the total number of data points. These options can only be specified for the first round.

generate = g

This specifies how to obtain the candidates at the start of the round:

If g is a positive integer, Maple will randomly create that many candidates. This can only be specified for the first round.

If g is the value merge, then the data points selected in all subsets in previous rounds are merged, and the candidates found are merged as well.

If g is the value all, then all original data points are used, including those potentially ignored due to subsets and size options in the first round. The candidates found are merged as well.

steps = t

This specifies the number of C-steps done on each candidate in this round. The default is 2.

select = k

This specifies that the top k candidates go on to the next round. This is the only required option: there is no default, because this option needs to be included for every round (if the subsets option is used at all).

The underlying computation is done in floating-point; therefore, all data points must have type realcons and all returned solutions are floating-point, even if the problem is specified with exact values.  For more information about numeric computation in the Statistics package, see the Statistics/Computation help page.

Another method for robust regression is implemented in the RepeatedMedianEstimator procedure.

$\mathrm{with}\left(\mathrm{Statistics}\right)\&colon;$

$x≔\mathrm{Sample}\left(\mathrm{Uniform}\left(0\&comma;10\right)\&comma;1000\right)$

$\mathrm{noise}≔\mathrm{Sample}\left(\mathrm{Cauchy}\left(0\&comma;1\right)\&comma;1000\right)$

$y≔\mathrm{`~`}\left[\mathrm{`+`}\right]\left(5x+\mathrm{noise}\&comma;\mathrm{`\; \$`}\&comma;10\right)$

$\mathrm{plots}\left[\mathrm{setoptions}\right]\left(\mathrm{size}=\left[0.7\&comma;''golden''\right]\right)\&colon;$

$\mathrm{pp}≔\mathrm{PointPlot}\left(y\&comma;\mathrm{xcoords}=x\right)\&colon;$$\mathrm{pp}$

$\mathrm{ls\_regression\_result}≔\mathrm{Fit}\left(aX+b\&comma;x\&comma;y\&comma;X\right)$

$\mathrm{ls\_regression\_result}≔3.44970682383807X+10.6816568681413$

$\mathrm{ls\_deviation\_from\_model}≔{\left(\mathrm{coeff}\left(\mathrm{ls\_regression\_result}\&comma;X\&comma;1\right)-5\right)}^2+{\left(\mathrm{coeff}\left(\mathrm{ls\_regression\_result}\&comma;X\&comma;0\right)-10\right)}^2$

$\mathrm{ls\_deviation\_from\_model}≔2.86806501793850$

$\mathrm{lts\_regression\_result}≔\mathrm{LeastTrimmedSquares}\left(x\&comma;y\&comma;\left[X\right]\right)$

$\mathrm{lts\_regression\_result}≔5.03537530551575X+9.82475419561271$

$\mathrm{lts\_deviation\_from\_model}≔{\left(\mathrm{coeff}\left(\mathrm{lts\_regression\_result}\&comma;X\&comma;1\right)-5\right)}^2+{\left(\mathrm{coeff}\left(\mathrm{lts\_regression\_result}\&comma;X\&comma;0\right)-10\right)}^2$

$\mathrm{lts\_deviation\_from\_model}≔0.0319625041956812$

$\mathrm{lts\_900\_regression\_result}≔\mathrm{LeastTrimmedSquares}\left(x\&comma;y\&comma;\left[X\right]\&comma;\mathrm{include}=900\right)$

$\mathrm{lts\_900\_regression\_result}≔5.00862730339998X+10.0156318668695$

$\mathrm{lts\_900\_deviation\_from\_model}≔{\left(\mathrm{coeff}\left(\mathrm{lts\_900\_regression\_result}\&comma;X\&comma;1\right)-5\right)}^2+{\left(\mathrm{coeff}\left(\mathrm{lts\_900\_regression\_result}\&comma;X\&comma;0\right)-10\right)}^2$

$\mathrm{lts\_900\_deviation\_from\_model}≔0.000318785625780940$

$\mathrm{rme\_regression\_result}≔\mathrm{RepeatedMedianEstimator}\left(x\&comma;y\&comma;X\right)$

$\mathrm{rme\_regression\_result}≔10.0306661686300+5.00564125476873X$

$\mathrm{rme\_deviation\_from\_model}≔{\left(\mathrm{coeff}\left(\mathrm{rme\_regression\_result}\&comma;X\&comma;1\right)-5\right)}^2+{\left(\mathrm{coeff}\left(\mathrm{rme\_regression\_result}\&comma;X\&comma;0\right)-10\right)}^2$

$\mathrm{rme\_deviation\_from\_model}≔0.000972237653807886$

$\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{pp}\&comma;\mathrm{plot}\left(\left[\mathrm{ls\_regression\_result}\&comma;\mathrm{lts\_regression\_result}\&comma;\mathrm{lts\_900\_regression\_result}\&comma;\mathrm{rme\_regression\_result}\right]\&comma;X=0..10\&comma;\mathrm{legend}=\left[''Least\; squares''\&comma;''Least\; trimmed\; squares''\&comma;''Least\; trimmed\; squares\; (900\; points)''\&comma;''Repeated\; median\; estimator''\right]\right)\&comma;\mathrm{view}=\left[0..10\&comma;-10..110\right]\right)$

$\mathrm{xy}≔\mathrm{Sample}\left(\mathrm{Normal}\left(0\&comma;1\right)\&comma;\left[1000\&comma;2\right]\right)$

$\mathrm{noise}≔\mathrm{Sample}\left(\mathrm{Normal}\left(0\&comma;1\right)\&comma;1000\right)$

$\mathrm{xy}\left[1..10\right]≔\mathrm{Sample}\left(\mathrm{Normal}\left(0\&comma;10\right)\&comma;\left[10\&comma;2\right]\right)$

${\mathrm{xy}}_{1..10}≔\left[\begin{array}{cc}\mathrm{-7.25835457416460}& 4.46546251109489\\ \mathrm{-12.8116871270602}& 6.95701976222452\\ 5.99013830391918& \mathrm{-0.0711220669337460}\\ 5.75154152981187& \mathrm{-24.7628908901018}\\ 9.09548316783929& 3.96336774375604\\ 8.73959176796014& 10.8028281369976\\ 10.7436650987406& \mathrm{-17.8756884460622}\\ 4.35976236509024& \mathrm{-11.7928910523180}\\ \mathrm{-5.77171156203078}& 4.33504635234648\\ 2.59252474066801& \mathrm{-10.1924397735107}\end{array}\right]$

$\mathrm{noise}\left[1..10\right]≔\mathrm{Sample}\left(\mathrm{Normal}\left(0\&comma;100\right)\&comma;10\right)$

${\mathrm{noise}}_{1..10}≔\left[\begin{array}{cccccccccc}100.800809666339& 8.96726548728334& \mathrm{-46.8119150928935}& \mathrm{-259.969027745269}& 19.0856519868389& \mathrm{-153.183461105396}& \mathrm{-25.7014908594142}& 94.1751893325151& 75.9056085428043& \mathrm{-120.552449126529}\end{array}\right]$

$z≔\mathrm{`~`}\left[\mathrm{`+`}\right]\left(\mathrm{xy}·⟨5\&comma;-3⟩\&comma;\mathrm{`\; \$`}\&comma;2\right)+{\mathrm{noise}}^{\mathrm{\%T}}$

$\mathrm{pp}≔\mathrm{plots}\left[\mathrm{pointplot3d}\right]\left(⟨\mathrm{xy}|z⟩\&comma;\mathrm{color}=\mathrm{black}\right)\&colon;$$\mathrm{pp}$

$\mathrm{ls\_regression\_result}≔\mathrm{Fit}\left(aX+bY+c\&comma;\mathrm{xy}\&comma;z\&comma;\left[X\&comma;Y\right]\right)$

$\mathrm{ls\_regression\_result}≔2.79396059575048X-0.846671122546459Y+1.81658520523718$

$\mathrm{lts\_regression\_result}≔\mathrm{LeastTrimmedSquares}\left(\mathrm{xy}\&comma;z\&comma;\left[X\&comma;Y\right]\right)$

$\mathrm{lts\_regression\_result}≔4.89975550872497X-3.01061932803024Y+2.16309086233403$

$\mathrm{plots}\left[\mathrm{display}\right]\left(\mathrm{pp}\&comma;\mathrm{plot3d}\left(\left[\mathrm{ls\_regression\_result}\&comma;\mathrm{lts\_regression\_result}\right]\&comma;X=-2..2\&comma;Y=-2..2\&comma;\mathrm{color}=\left[\mathrm{red}\&comma;\mathrm{blue}\right]\right)\&comma;\mathrm{view}=\left[-2..2\&comma;-2..2\&comma;-5..20\right]\&comma;\mathrm{orientation}=\left[25\&comma;80\&comma;-60\right]\right)$

$x≔\mathrm{Sample}\left(\mathrm{Uniform}\left(0\&comma;1\right)\&comma;\left[100000\&comma;10\right]\right)$

$y≔x·⟨1\&comma;2\&comma;3\&comma;4\&comma;5\&comma;6\&comma;7\&comma;8\&comma;9\&comma;10⟩$

$y\left[1..1000\right]≔\mathrm{`~`}\left[\mathrm{`+`}\right]\left(y\left[1..1000\right]\&comma;\mathrm{`\; \$`}\&comma;10\right)$

$\mathrm{CodeTools}\left[\mathrm{Usage}\right]\left(\mathrm{LeastTrimmedSquares}\left(x\&comma;y\&comma;X\right)\right)$

memory used=298.15MiB, alloc change=45.50MiB, cpu time=4.43s, real time=1.80s, gc time=1.56s

$1.00000000000002X_1+1.99999999999999X_2+2.99999999999999X_3+4.00000000000000X_4+5.00000000000000X_5+5.99999999999999X_6+7.00000000000000X_7+8.00000000000002X_8+8.99999999999999X_9+10.0000000000000X_{10}+8.20789049648978\times {10}^{\mathrm{-16}}$

[1] Rousseeuw, Peter J. Least Median of Squares Regression. Journal of the American Statistical Association 79, 1984, pp.871-880.

[2] Rousseeuw, Peter J., and Van Driessen, Katrien. Computing LTS Regression for Large Data Sets. Data Mining and Knowledge Discovery 12, 2006, pp.29-45.

The Statistics[LeastTrimmedSquares] command was introduced in Maple 2019.

For more information on Maple 2019 changes, see Updates in Maple 2019.