| Library: | metrics |
| See also: |
| Quantlet: | lts | |
| Description: | Computes the least trimmed squares estimate for the coefficients of a linear model. |
| Usage: | b = lts(x, y{, h, all, mult}) | |
| Input: | ||
| x | n x p design matrix of explanatory variables. | |
| y | n x 1 vector, dependent variable. | |
| h | optional, trimming constant; default value is [n/2]+[(p+1)/2], where n is the number of observations and p is the number of parameters. h should belong to {[(n+1)/2],...,n} and must be bigger than p. | |
| all | optional, logical flag for the exact computation (nonzero = TRUE), default value is 0 (FALSE). If ci = 0, an approximation of the estimator is computed. If ci != 0, the estimate is computed precisely - this can be rather time demanding and applicable only for small sample sizes. The number of iterations corresponds in this case to n over h. | |
| mult | optional, affects the maximal number of iterations, after which the algorithm is stoped (if convergence was not achieved earlier), default value equals to 1. The maximal number of iterations is currently '(600 * h) * times', so this variable allow to reduce or extend this number. | |
| Output: | ||
| b | p x 1 vector of estimated coefficients. | |
library("metrics")
;
; simulate data
;
randomize(555)
x = uniform(50,3)
y = x[,1] + 2*x[,2] - x[,3] + normal(50)
;
; fit the data...
;
b = lts(x,y)
b
Contents of coefs - estimates of b = (1,2,-1)' coefficient vector [1,] 0.95542 [2,] 2.0811 [3,] -0.25978
; procedure for graphical representation of the results estimated by rqfit quantlet
; parameter obs1 is an outlier added to randomly generated data points
proc() = myquant(obs1)
;
; initialization
;
n = 10 ; number of observations
randomize(17654321) ; sets random seed
beta = #(1, 2) ; defines intercept and slope
x = matrix(n)~sort(uniform(n))
; ; creates design matrix
;
; new x-observation is added
;
x = x|(1~obs1[1])
m = x*beta ; defines regression line
eps = 0.05*normal(n) ; creates obs error
;
; new y-observation is added
;
y = m[1:n] + eps ; noisy line
y = y|obs1[2]
;
; create graphical display and draw data points
;
d = createdisplay(1,1)
dat = x[,2]~y
outl = obs1[1]~obs1[2]
setmaskp(outl,1,12,15) ; outlier is blue big star
tdat = x[,2]~m
setmaskl(tdat,(1:rows(tdat))', 1, 1, 1)
; ; thin blue line
setmaskp(tdat, 0, 0, 0) ; reduces point size to min
;
; estimation of the model using rqfit
;
beta1 = lts(x,y)
;
; draw estimated regression line
;
yhat = x*beta1
hdat = x[,2]~yhat
setmaskp(hdat, 0, 0, 0)
setmaskl(hdat,(1:rows(hdat))', 4, 1, 3)
; ; thick red line
show(d, 1, 1, dat[1:n], outl, tdat, hdat)
title="LTS regression with an outlier"
setgopt(d,1,1,"title",title)
; ; sets title
endp ; end of myquant
;
; load metrics library
;
library("metrics")
;
; call estimation function with outlier #(0.9,4.5)
;
myquant(#(0.9,4.5))
As a result, you should see a graph, in which observations are denoted by black circles and an outlier (observation #(0.9,4.5)) is represented by the blue big star in the right upper corner of the graph. The blue line depicts the true regression line (beta = #(1,2)), while the thick red line shows the estimated regression line.