; 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))