| Library: | robtech |
| See also: | lts |
| Quantlet: | lms | |
| Description: | Computes the least median of squares estimate for the coefficients of a linear model. |
| Usage: | b = lms(x, y{, nrep}) | |
| Input: | ||
| x | n x p design matrix of explanatory variables. | |
| y | n x 1 vector, dependent variable. | |
| nrep | an optional scalar, identifying how many random subsamples will be analyzed for the LMS estimate; the higher, the more precise estimate is. By default, it is equal to 1000 times the number of explanatory variables. | |
| Output: | ||
| b | p x 1 vector of estimated coefficients. | |
library("robtech")
;
; simulate data
;
n = 100
randomize(1)
;
x = matrix(n)~uniform(n,2)
y = 3 - 2*x[,2] + x[,3] + normal(n)
;
; fit the data...
;
b = lms(x,y)
b
Contents of coefs - estimates of b = (3,-2,1)' coefficient vector [1,] 2.9133 [2,] -2.411 [3,] 0.79301
; procedure for graphical representation of the results estimated by lms 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 = lms(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="LMS regression with an outlier"
setgopt(d,1,1,"title",title)
; ; sets title
endp ; end of myquant
;
; load metrics library
;
library("robtech")
;
; 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.