1

Wyrównanie odporne na

bł

ę

dy grube. Wpływ danych

na modele prostej regresji

Wykonawca:

Grzegorz Kruczek

Grupa

ć

wiczeniowa 1

Rok akademicki 2015/2016

Numeryczne algorytmy in

ż

ynierskie

2

Spis tre

ś

ci

Dane wej

ś

ciowe...........................................................................................................3

Tok post

ę

powania........................................................................................................4

3

Dane wej

ś

ciowe

x

y

x

y

1

0.8159985

0.845318602

21 20.29533811 20.280916953

2

1.11419088

1.116836595

22 21.61858242 21.620956198

3

2.7380354

2.733633973

23 22.62647535 22.679587383

4

3.19497711

3.232829583

24 23.66728777 23.642666225

5 28.04798287

4.031829307

25 24.18528434 24.155247356

6

5.50007597

5.536825322

26 25.01542013 25.007648092

7

6.20189719

6.180988191

27 26.34964609 26.369400241

8

7.10734958

7.109668801

28 27.05127128 27.064629659

9

8.36876885

8.327941511

29 28.09764659 28.032069515

10

9.83088126

9.773894612

30 29.56409045

29.45576344

11 10.78951177 10.721589685

31 30.44990471 30.491950953

12 11.20967706 11.229551574

32 31.46518095 31.494271634

13 12.01225721 12.018654188

33 32.31663455 32.371260719

14 13.25754969 39.311144271

34 33.46146807 33.459222876

15 14.59213905 14.611905216

35 34.85470867 13.898829534

16 15.41742361 15.389581572

36 35.01600899 35.038755136

17 16.45590928 16.430955864

37 36.65967839 36.720278796

18 17.97281482 17.929151577

38 37.32772208 37.289466864

19 18.16006091 18.157476287

39 30.78828089 38.781149378

20

0.27834333 19.279055692

40 39.36120673 39.434106016

4

Tok post

ę

powania:

1: Wybrany program: R v3.2.0

2. Wczytanie współrz

ę

dnych x i y punktów z oddzielnych plików Mxn i Myn:

x=scan(file="Mxn.txt")
Read 40 items
y=scan(file="Myn.txt")
Read 40 items

Program potwierdził wczytanie zało

ż

onej liczny elementów.

3. Wyplotowanie danych:

plot(y~x)

4. Utworzenie modelu regresji liniowej:

ramka danych: ds=data.frame(x=x,y=y)
regresja liniowa: m.lm=lm(y~x,data=ds)
podsumowanie modelu: summary(m.lm)

Call:
lm(formula = y ~ x, data = ds)

5

Residuals:
Min 1Q Median 3Q Max
-22.7221 -2.1387 -0.0017 1.8158 24.6576

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.80728 2.25089 1.691 0.0989 .
x 0.81812 0.09815 8.335 4.16e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.15 on 38 degrees of freedom
Multiple R-squared: 0.6464, Adjusted R-squared: 0.6371
F-statistic: 69.48 on 1 and 38 DF, p-value: 4.158e-10

Model ma du

ż

e odchylenie standardowe.

Wyplotowanie danych z wpasowan

ą

funkcj

ą

sklejan

ą

:

layout(1)
plot(x~y)
lines(lowess(y,x))

6

definicja wydruku: op=par(mfrow=c(2,2),pty="s")
wyplotowanie wykresów dotycz

ą

cych modelu: plot(m.lm)

Te oraz poprzednie wykresy wyra

ź

nie wskazuj

ą

punkty odstaj

ą

ce od modelu i

zaburzaj

ą

ce go.

7

wczytanie nowej biblioteki: library(MASS)
obliczenie odległo

ś

ci Cooka: dl=cooks.distance(m.lm)

obliczenie warto

ś

ci zestandaryzowanych reszt: r=stdres(m.lm)

wytypowanie punktów odległych o wi

ę

cej ni

ż

jedn

ą

jednostk

ę

odległo

ś

ci od

wykresu funkcji obrazuj

ą

cego model:

a=cbind(ds,dl,r)
a[dl>1/40,]

x

y

dl

r

5 28,04798 4,031829 0,205732 -3,23971

14 13,25755 39,31114 0,210838 3,507321
20 0,278343 19,27906 0,270538 2,243748
35 34,85471 13,89883 0,257879 -2,66853
39 30,78828 38,78115 0,049177 1,402413

Znalezienie w zbiorze 10 punktów o najwi

ę

kszych odległo

ś

ciach Cooka

w kolejno

ś

ci od najwi

ę

kszej do najmniejszej:

rabs=abs(r)
a=cbind(ds,dl,r,rabs)
asorted=a[order(-rabs),]
asorted[1:10,]

x

y

dl

r

rabs

14 13,25755 39,31114 0,210838 3,507321 3,507321

5 28,04798 4,031829 0,205732 -3,23971

3,23971

35 34,85471 13,89883 0,257879 -2,66853 2,668527
20 0,278343 19,27906 0,270538 2,243748 2,243748
39 30,78828 38,78115 0,049177 1,402413 1,402413

1 0,815999 0,845319 0,014593 -0,53308 0,533076
2 1,114191 1,116837 0,013979 -0,52841 0,528411

40 39,36121 39,43411 0,013625 0,504005 0,504005

3 2,738035 2,733634

0,01016 -0,48321 0,483214

4 3,194977

3,23283 0,009008

-0,4642 0,464202

wczytanie nowej biblioteki: library(car)
wyplotowanie wykresu przedstawiaj

ą

cego odległo

ś

ci Cooka dla ka

ż

dego

punktu: influencePlot(lm(y~x),main="influencePlot",sub="Promie

ń

Kółka jest

proporcjonalny do odległo

ś

ci Cooka")

StudRes Hat CookD
14 4.208437 0.03314286 0.4591710
20 2.377084 0.09704551 0.5201324

8

Automatyczne wyszukanie punktów zaburzaj

ą

cych model:

influence.measures(lm(y~x))
Influence measures of
lm(formula = y ~ x) :

dfb,1_

dfb,x

dffit

cov,r

cook,d

hat

inf

1 -0,169179

0,144728 -0,16921

1,146

1,46E-02

0,0931

2 -0,165549

0,141037

-0,1656

1,143

1,40E-02

0,091

3 -0,140781

0,11701

-0,1411

1,133

1,02E-02

0,0801

4 -0,132403

0,109203 -0,13282

1,13

9,01E-03

0,0772

5

0,069456 -0,432087 -0,74399

0,574

2,06E-01

0,0377

*

6 -0,102106

0,080501 -0,10328

1,117

5,45E-03

0,0637

7 -0,095814

0,074323 -0,09731

1,113

4,84E-03

0,06

8 -0,084745

0,064231 -0,08664

1,109

3,84E-03

0,0555

9 -0,072901

0,053235 -0,07547

1,104

2,92E-03

0,0498

10 -0,059443

0,041188 -0,06283

1,098

2,02E-03

0,0438

9

11 -0,051424

0,03418 -0,05535

1,095

1,57E-03

0,0404

12

-0,04568

0,029745 -0,04964

1,094

1,26E-03

0,039

13 -0,039898

0,024865 -0,04427

1,092

1,01E-03

0,0365

14

0,673845 -0,386214

0,77917

0,499

2,11E-01

0,0331

*

15 -0,022849

0,011605 -0,02803

1,086

4,03E-04

0,0302

16 -0,019314

0,008878 -0,02481

1,085

3,16E-04

0,0287

17 -0,014228

0,005517 -0,01962

1,083

1,98E-04

0,0271

18 -0,008367

0,002103 -0,01321

1,082

8,95E-05

0,0257

19 -0,007126

0,001647 -0,01147

1,082

6,76E-05

0,0255

20

0,779276 -0,671452

0,77929

0,879

2,71E-01

0,097

*

21 -0,001365 -0,000118 -0,00292

1,082

4,38E-06

0,025

22

0,001047

0,000442

0,00288

1,082

4,26E-06

0,0256

23

0,002366

0,001964

0,00833

1,083

3,56E-05

0,0265

24

0,002273

0,003534

0,01118

1,085

6,42E-05

0,0278

25

0,002211

0,004768

0,01348

1,085

9,34E-05

0,0286

26

0,001868

0,007441

0,01813

1,087

1,69E-04

0,0301

27

0,000291

0,012838

0,02606

1,09

3,49E-04

0,033

28 -0,001015

0,015966

0,03006

1,091

4,64E-04

0,0348

29 -0,003325

0,020152

0,03456

1,095

6,13E-04

0,0379

30 -0,007618

0,028172

0,04365

1,099

9,77E-04

0,0429

31 -0,011977

0,03743

0,05522

1,102

1,56E-03

0,0463

32 -0,016608

0,04519

0,06359

1,106

2,07E-03

0,0505

33 -0,021365

0,053224

0,07241

1,11

2,69E-03

0,0544

34 -0,027577

0,062627

0,08199

1,116

3,44E-03

0,06

35

0,299273 -0,623883 -0,78613

0,747

2,58E-01

0,0675

*

36

-0,03868

0,079953

0,10035

1,124

5,15E-03

0,0685

37 -0,052771

0,101289

0,12274

1,134

7,70E-03

0,0784

38 -0,056954

0,106631

0,12766

1,139

8,33E-03

0,0827

39 -0,073774

0,219041

0,31779

0,996

4,92E-02

0,0476

40 -0,080038

0,140778

0,16344

1,152

1,36E-02

0,0969

Punkty zaburzaj

ą

ce model zostały oznaczone symbolem *

5. Regresja liniowa z "r

ę

cznym" usuni

ę

ciem danych:

Usuni

ę

cie punktów o najwi

ę

kszych odległo

ś

ciach Cooka:

m.lm2=lm(y~x,subset=-c(5,14,20,35,39))
Podsumowanie modelu: summary(m.lm2)

Call:
lm(formula = y ~ x, subset = -c(5, 14, 20, 35, 39))

Residuals:
Min 1Q Median 3Q Max
-0.111728 -0.024608 0.005648 0.025206 0.064585

Coefficients:

10

Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0114228 0.0137374 -0.832 0.412
x 1.0005014 0.0006073 1647.588 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04058 on 33 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 2.715e+06 on 1 and 33 DF, p-value: < 2.2e-16

Odchylenie standardowe modelu po usuni

ę

ciu punktów o najwi

ę

kszych

odległo

ś

ciach Cooka zmniejszyło si

ę

o dwa rz

ę

dy wielko

ś

ci.

wyplotowanie wykresu izolinii odległo

ś

ci Cooka: plot(m.lm2)

obliczenie odległo

ś

ci Cooka: dl.2=cooks.distance(m.lm2)

obliczenie warto

ś

ci zlinearyzowanych reszt: r2=stdres(m.lm2)

s=seq(1,40,by=1)
s=s[-c(5,14,20,35,39)]

11

a2=cbind(s,dl.2,r2)
wytypowanie punktów odległych o wi

ę

cej ni

ż

jedn

ą

jednostk

ę

odległo

ś

ci od

wykresu funkcji obrazuj

ą

cego model:

a2[dl.2>1/35,]

s

dl,2

r2

1

1 0,0667249 1,052129

4

4 0,0738404 1,230724

6

6 0,0532742 1,162406

10 10 0,0428335 -1,276564
11 11 0,0587549

-1,56191

29 29 0,069328 -1,720609
30 30 0,2137143

-2,82598

33 33 0,0558723

1,27012

37 37 0,0997115 1,388527
38 38 0,0767688 -1,182474
40 40 0,1880271 1,692757

Ż

aden z dziesi

ę

ciu najbardziej odległych punktów nie spełnił powy

ż

szego kryterium.

obliczenie warto

ś

ci bezwzgl

ę

dnych reszt zestandaryzowanych: rabs2=abs(r2)

a2=cbind(s, dl.2, r2, rabs2)
posortowanie malej

ą

co: asorted2=a2[order(-rabs2), ]

asorted2[1:10,]

s

dl,2

r2

rabs2

30 0,21371426

-2,82598

2,82598

29 0,06932798 -1,720609

1,720609

40 0,18802708

1,692757

1,692757

11 0,05875494

-1,56191

1,56191

37 0,09971154

1,388527

1,388527

23 0,02799205

1,331256

1,331256

10 0,04283348 -1,276564

1,276564

33 0,05587234

1,27012

1,27012

4 0,07384035

1,230724

1,230724

38 0,07676884 -1,182474

1,182474

6. Regresja z modelem Hubera:

library(MASS)
m.hub=rlm(y~x)
summary(m.hub)

Call: rlm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-2.402e+01 -3.000e-02 8.793e-04 3.024e-02 2.605e+01

12

Coefficients:
Value Std. Error t value
(Intercept) 0.0008 0.0167 0.0508
x 1.0000 0.0007 1373.7983

wyplotowanie wykresu regresji z modelem Hubera: plot(m.hub)

Wykres wyra

ź

nie przedstawia odstajace punkty.

13

wyplotowanie wykresu dla wag Hubera: plot(m.hub$w, ylab="Wagi Hubera")

Wykres wyra

ź

nie przedstawia punkty odstaj

ą

ce od modelu.

a=cbind(ds, dl, r, m.hub$w)
asorted=a[order(m.hub$w),]
asorted3=a[order(m.hub$w),]
asorted3[1:10,]

x

y

dl

r

m,hub$w

14 13,2575497 39,3111443 0,210837967 3,5073213 0,002330115

5 28,0479829

4,0318293 0,205732127 -3,2397097 0,002527485

35 34,8547087 13,8988295 0,257878892

-2,668527 0,002896538

20 0,2783433 19,2790557 0,270537736 2,2437482 0,003195056
39 30,7882809 38,7811494 0,049176639 1,4024133 0,007596463
30 29,5640904 29,4557634 0,000977148 0,2089273 0,551844067
40 39,3612067

39,434106 0,013624927 0,5040054 0,854401701

11 10,7895118 10,7215897 0,001570374

-0,273112

0,87640733

29 28,0976466 28,0320695

0,00061289 0,1764515 0,902976017

1 0,8159985

0,8453186 0,014592603 -0,5330759

1

Najmniejsze wagi Hubera otrzymały punkty najbardziej odstaj

ą

ce od modelu.

14

7. Regresja z modelem Bisquare

m.bi=rlm(y~x,method='MM')
m.bi

Call:
rlm(formula = y ~ x, method = "MM")
Converged in 3 iterations

Coefficients:
(Intercept) x
-0.01149484 1.00057661

Degrees of freedom: 40 total; 38 residual
Scale estimate: 0.0507

wyplotowanie wykresu dla wag Bisquare:
plot(m.bi$w, ylab="Wagi Bis",type="p",lwd=10)

Wykres ten podobnie jak poprzedni pozwala na łatw

ą

, cho

ć

niej jednoznaczn

ą

,

identyfikacje punktów odstaj

ą

cych od modelu.