4.1 Regresní úloha a regresní funkce

Lekce 4




Metoda nejmenších čtverců

Metoda nejmenších čtverců je další z „výkladních skříní“ statistiky. My se seznámíme pouze s její nejzákladnější verzí, kdy regresní funkce, měřící průběh závislosti, je funkcí jedné proměnné lineární v parametrech, případně pro ni existuje linearizující transformace. Budeme pomocí ní řešit regresní úlohu, což je úloha, v níž vysvětlující proměnná je řízená veličina a náhodnou veličinou je pouze vysvětlovaná proměnná. Součástí metody nejmenších čtverců je rovněž rozklad součtu čtverců odchylek pozorovaných hodnot vysvětlované proměnné od jejich průměru na složku vysvětlenou regresí a nevysvětlenou — reziduální — složku. Rozklad součtu čtverců odchylek nás přivede k charakteristice intenzity závislosti — indexu korelace. Čtvercem indexu korelace je index determinace, který vyjadřuje (zpravidla v %), jakou část rozptylu pozorovaných hodnot vysvětlované proměnné se podařilo objasnit regresí.

adjustace; funkce lineární v parametrech; Gaussovy normální rovnice; index determinace; index korelace; kritérium nejmenších čtverců; linearizující transformace; matice regresorů; metoda nejmenších čtverců; regresní funkce; reziduální součet čtverců; řízená proměnná; vysvětlovaná proměnná; vysvětlující proměnná

4.1 Regresní úloha a regresní funkce Datový soubor je tvořen n uspořádanými dvojicemi hodnot [xi , yi ] . Zatímco hodnoty yi jsou realizace (pozorované hodnoty) náhodné veličiny Y, hodnoty xi jsou hodnoty řízené proměnné, jejichž velikost je určena experimentátorem. Náhodná veličina Y se nazývá vysvětlovaná proměnná, hodnoty xi představují hodnoty vysvětlující proměnné. Cílem úlohy je určit průběh závislosti vysvětlované proměnné na měnících se hodnotách vysvětlující proměnné pomocí spojité regresní funkce a změřit intenzitu této závislosti. Pokud jde o regresní funkce, omezíme se na funkce jedné proměnné, které jsou navíc lineární v parametrech. Takovouto funkci můžeme zapsat jako součet součinů neznámých parametrů b j a známých funkcí — regresorů — vysvětlující proměnné. Jeden z parametrů, označovaný b0 , je absolutní člen regresní funkce. Regresní funkci jedné proměnné lineární v parametrech zapíšeme jako m

y ′ = ∑ b j f j ( x ) , kde [b0 b1 ... bm ] je vektor parametrů a funkce f j (x ) , pro j = 0,1,2,..., m , které T

j =0

neobsahují žádné další neznámé parametry, jsou regresory. Regresor f 0 ( x ) = 1( = x 0 ), další regresory jsou zpravidla elementární funkce hodnot x, např. x 1 = x, x −1 =

1 2 , x , log x, x ,.... . Mezi typické x

regresní funkce patří např. polynom 1. stupně (lineární funkce, přímka), 2. stupně (kvadratická funkce, parabola), 3. stupně (kubická parabola) a mnohé jiné. Existuje ovšem mnoho funkcí, které bychom potenciálně chtěli použít při řešení regresní úlohy, které však nejsou lineární v parametrech, např. funkce b0 b1x , b0 x b1 ,

1 , b0 + b1x . b0 + b1 x

Řadu těchto funkcí lze podrobit linearizující transformaci (např. všechny uvedené s výjimkou poslední). Pro některé však linearizující transformace, jak jsme právě viděli, neexistuje.

22

Příklad 4.1

1 . b0 + b1 x Logaritmická transformace exponenciální funkce: log y ′ = c0 + c1 x , kde c0 = log b0 , c1 = log b1 . 1 Lomenou funkci transformujeme jako = b0 + b1 x . y′

Budeme linearizovat regresní funkce y ′ = b0 b1x (exponenciální funkce) a y ′ =

Linearizovanou funkci obecně označíme ϕ ( y ′) =

m

∑b j =0

j

f j ( x ) , kde ϕ je linearizující trans-

formační funkce.




y ′ = b0 x b1 , y ′ = b0 + b1 x . (4–1)

Po vzoru příkladu 4.1 linearizujeme funkce

Příklad 4.2 Neexistence linearizující transformace. Např. pro funkci y ′ = b0 + b1x neexistuje transformace, která by umožnila vyjádřit ji ve výše uvedex

ném tvaru. Podobně je tomu např. u funkce

y ′ = b0b1

.

Případy funkcí, pro které neexistuje linearizující transformace, se nebudeme zabývat.

4.2 Metoda nejmenších čtverců Sestavíme vektor pozorovaných hodnot závisle proměnné y T = [ y1 y 2 ... y n ] a matici regresorů

1 1 F= ...  1

f 1 ( x1 ) f1 ( x2 ) ... f1 ( xn )

... ... ... ...

f m ( x1 )  f m ( x2 )  (první sloupec souvisí s tím, že funkce má absolutní člen). ...   f m ( xn )

Vektor parametrů regresní funkce určíme jako b = (F T F) −1 F T y . Příklad 4.3 Nalezení vektoru parametrů pro kvadratickou regresní funkci y ′ = b0 + b1 x + b2 x 2 .

 1 x1  1 x2 Matice regresorů F =  ... ...   1 xn

x12   n  x22   T −1 , matice ( F F) =  ∑ xi ...  ∑ xi2   xn2 

 ∑ yi    tor F T y =  ∑ yi xi  . ∑ yi xi2   

23

∑x ∑x ∑x ∑x ∑x ∑x i 2 i 3 i

2 i 3 i 4 i

−1

   , a konečně vek 

∑y ∑y x ∑y x

b0  Vektor b =  b1  je tedy řešením soustavy tzv. Gaussových normálních rovnic   b2  − nb0 − b1 ∑ xi − b2 ∑ xi2 = 0

i i

i

− b0 ∑ xi − b1 ∑ xi2 − b2 ∑ xi3 = 0

i

2 i

− b0 ∑ xi2 − b1 ∑ xi3 − b2 ∑ xi4 = 0 n

n

(všechny součty jsou

∑

), jejíž j–tá rovnice je dána jako

i =1




∑y i =1

i

m

n

j =0

i =1

f j ( xi ) − ∑ b j ∑ x i f j ( x i ) = 0 .

Po vzoru příkladu 4.3 ukážeme matici regresorů a sestavíme soustavu normálních rovnic pro regresní funkci

y ′ = b0 +

b1 . (4–2) x

Lze ukázat, že tímto způsobem získané parametry b0 , b2 ,..., bm po dosazení do regresní funkce lineární v parametrech y ′ =

m

∑b j =0

j

f j ( x ) vedou ke splnění dvou podmínek

n n m    ∑ ( yi − yi′ ) = ∑  yi − ∑ f j ( xi ) = 0 (z toho plyne y = y ′ ), i =1 i =1  j =0  2

m    ∑ ( yi − yi′ ) = ∑  yi − ∑ f j ( xi ) → min . i =1 i =1  j =0  n

n

2

Druhý z výrazů představuje kritérium nejmenších čtverců, tzv. reziduální součet čtverců. To, že výraz nabývá své minimální hodnoty právě pro parametry b0 , b2 ,..., bm , vyplývá z toho, že j-tá normální rovnice

n

m

n

i =1

j =0

i =1

∑ yi f j ( xi ) − ∑ b j ∑ xi f j ( xi ) = 0 vznikne jako parciální derivace kritéria podle parametru

b j , položená rovná nule. Odtud název metoda nejmenších čtverců. Zatímco existuje nekonečně mnoho parametrů b0 , b2 ,..., bm , pro které je splněna první z obou podmínek, existuje pro každý typ funkce (lineární, kvadratická, lomená apod.) jen jediná sada parametrů, která splňuje druhou podmínku. Určení konkrétního typu funkce není součástí úlohy a je v rukou uživatele. Vektor parametrů regresní funkce b0 , b2 ,..., bm , které byly stanoveny výše popsaným postupem, jednoznačným způsobem minimalizuje součet čtverců odchylek pozorovaných a naměřených hodnot závisle proměnné pro předem zvolený typ regresní funkce. Příklad 4.4 Při zkoušení turbíny byly experimentátorem voleny otáčky a měřena frekvence chvění2. Pro každou zvolenou hodnotu měrných otáček (vysvětlující proměnná) bylo provedeno měření jedné až tří hodnot měrné frekvence chvění (vysvětlovaná proměnná). Naměřené hodnoty jsou znázorněny na obr. 4.1. Z obrázku vyplývá, že vhodnou regresní funkcí by mohla být kvadratická funkce. Stanovíme tedy její rovnici. 2

V úloze jsou použity tzv. měrné veličiny. Každou rozměrnou fyzikální veličinu můžeme převést na

bezrozměrnou a na intervalu

0;1 normovanou měrnou veličinu pomocí vztahu

24

xi − x min . x max − x min .

Obr. 4.1 Naměřené hodnoty měrných otáček a měrné frekvence chvění a uvažované regresní funkce

K určení regresní paraboly (přímka je evidentně nevhodná) stanovíme z naměřených hodnot tyto veličiny

n = 16; ∑ xi = 8,35; ∑ xi2 = 5,0525; ∑ xi3 = 3,365875; ∑ xi4 = 2,37453125;

∑y

i

= 6,515; ∑ xi yi = 3,27005; ∑ xi2 yi = 3,365875

Tyto použijeme k sestavení soustavy normálních rovnic z příkladu 4.3

6,515

− 16b0

3,27005 − 8,35b0

− 8,35b1

− 5,0525b2

=0

− 5,0525b1

− 3,365875b2

=0,

3,365875 − 5,0525b0 − 3,365875b1 − 2,37453125b2 = 0 jejímž řešením obdržíme b0 = 0,9745; b1 = −2,3130; b2 = 2,0261 ; z čehož rovnice kvadratické funkce je y ′ = 0,9745 − 2,3130 x + 2,0261x 2 (viz obr. 4.2) Obr. 4.2 Naměřené hodnoty a vypočtená regresní parabola




Určete při jakých měrných otáčkách bude měrné chvění minimální a stanovte jeho hodnotu.

25

4.3 Rozklad součtu čtverců a index korelace Reziduální součet čtverců lze vyjádřit rovněž jako rozdíl součtu čtverců pozorovaných hodnot yi kolem aritmetického průměru a součtu čtverců vyrovnaných (vypočtených) hodnot kolem aritmetického průměru (při tom y = y ′ ):

n

n

n

i =1

i =1

i =1

∑ ( yi − yi′) 2 =∑ ( yi − y ) 2 − ∑ ( yi′ − y ) 2 .

Rovnice rozkladu součtu čtverců pozorovaných hodnot vysvětlované proměnné rozkládá variabilitu pozorovaných hodnot vysvětlované proměnné na složku vysvětlenou regresí a složku reziduální. Na tomto principu je založeno měření intenzity závislosti pomocí indexu korelace Obr. 4.2 Rozklad součtu čtverců vysvětlované proměnné

( yi − y)

y

n

I yx =

∑ ( y′ − y)

2

∑( y

2

i =1 n

2

i =1

i

i

− y)

, 0 ≤ I yx ≤ 1 .

Je-li regresí vysvětlena veškerá variabilita pozorovaných hodnot vysvětlované proměnné, je I yx = 1 a všechny pozorované hodnoty leží přesně

y

na regresní čáře. Není-li regresí vysvětlena žádná variabilita, je I yx = 0 (např. mění-li se hodnoty xi , zůstává yi = konst. ) — viz obr. 4.3.

x

a) čtverce odchylek pozorovaných hodnot kolem jejich aritmetického průměru

y

y′

y′ ( yi − yi′ ) 2

y

( y i′ − y ) 2

y

x

x b) čtverce odchylek vyrovnaných hodnot kolem jejich aritmetického průměru

c) čtverce odchylek pozorovaných hodnot kolem hodnot vyrovnaných

n

n

n

i =1

i =1

i =1

∑ ( y i − y ) 2 = ∑ ( y i′ − y ) 2 + ∑ ( y i − y i′ ) 2

26

Obr. 4.3 Extrémní případy indexu korelace

y

y

I yx = 1

y ′ = konst.

I yx = 0

y i = y i′

x




x

n

K nákresům na obr. 4.3 přiřaďte

∑(y i =1

i

n

n

i =1

i =1

− y ) 2 = ∑ ( y i − y i′ ) 2 , ∑ ( y i − y i′ ) 2 = 0 .

Příklad 4.5 Výpočet indexu korelace Rozklad rovnice součtu čtverců odchylek pro vysvětlovanou proměnnou vedoucí k indexu korelace v příkladu 4.4 (jen v náznaku): 0,0907244 = 0,081770 + 0,0089544 a pak

I yx =

0,081770 = 0,9494 . Podle očekávání je závislost poměrně těsná. 0,0907244

Testování významnosti indexu korelace Hypotézu o nulové hodnotě indexu korelace proti alternativní hypotéze ověříme pomocí testového kri-

I yx2 téria F =

m 1 − I yx2

, kde m je počet parametrů regresní funkce (tedy s výjimkou b0 ). Testové krité-

n − m −1 rium má Fisherovo–Snedecorovo rozdělení s m a (n – m – 1). Jde o jednostranný test, takže vypočtená hodnota se porovnává s tabelovanou hodnotou pro F1−α [m; n − m − 1] . Příklad 4.6 Test významnosti indexu korelace Ověříme významnost indexu korelace z příkladu 4.5

0,9013 2 F= = 59,36 . Hypotézu tedy na všech myslitelných hladinách významnosti zamítáme, 1 − 0,9013 16 − 3 neboť pravděpodobnost, že index korelace nabyl své hodnoty náhodou je přibližně 3 ⋅ 10 −8 . Při posuzování statistické významnosti indexu korelace je třeba vzít v úvahu, že jeho nízká, nevýznamná, hodnota nemusí být nutně způsobena neexistencí závislosti vysvětlované proměnné, ale příčina může být v nevhodné volbě typu funkce (závislost existuje, jen zvolená funkce její průběh měří špatně). Pokud bychom např. v příkladu 4.4 použili k vyrovnání přímku, získali bychom index korelace I yx = 0,52 , což je zavádějící hodnota ( F = 4,39 a pravděpodobnost je 0,0349, takže hypotézu

27

o nulové hodnotě indexu korelace nelze na hladině α = 0,01 zamítnout). Poslední výpočty byly provedeny s použitím specializovaného programového vybavení a čtenář si je nemůže zkontrolovat. Adjustace indexu korelace Dalším problémem je nesrovnatelnost hodnot indexu korelace u funkcí s různým počtem parametrů. Vzhledem k tomu, že regresní funkce s větším počtem parametrů má (při volbě vhodného typu funkce) lepší předpoklady vystihnout průběh funkce, má přidání dalšího parametru za následek zvýšení hodnoty indexu korelace. Proto se provádí tzv. adjustace indexu korelace, která snižuje jeho hodnotu v zá2 vislosti na počtu parametrů funkce (s výjimkou parametru b0 ): I adj = 1 − (1 − I yx )

n −1 . Ze n − (m + 1)

vzorce je zřejmé, že adjustace má význam zejména při malém počtu měření.

Příklad 4.7 Adjustace indexu korelace Provedeme adjustaci vypočteného indexu korelace z příkladu 4.5

I adj = 1 − (1 − 0,9013)

15 = 0,9413 . Hodnota adjustovaného indexu se tedy nepatrně snížila. 13

Index determinace n

2 2 Výše zmíněná hodnota I yx , často vyjadřovaná 100 I yx =

∑ ( y ′ − y)

2

i

i =1 n

∑(y i =1

i

− y) 2

100 , kde 0 ≤ 100 I yx2 ≤ 100 ,

udává, jaká část variability (v %) byla vysvětlena regresí a nazývá se index determinace.

4.4

Problémy linearizující transformace

Použití linearizující transformace u funkcí, které nejsou lineární v parametrech (řada z nich se běžně používá), není bez problémů. Je-li funkce lineární v parametrech ve tvaru ϕ ( y ′) =

m

∑b j =0

j

f j ( x ), kde ϕ je linearizující trans2

n

formační funkce, je reziduální součet čtverců odchylek roven

∑ [ϕ ( y ) − ϕ ( y ′ )] i =1

n

∑ (log y i =1

n

i

− log y i′ ) 2 → min ., ∑ ( i =1

i

i

→ min . (tj. např.

1 1 − ) 2 → min . apod.). Je zřejmé, že všechna tato kritéria y i y i′

jsou nesrovnatelná s původním kritériem z odst. 4.2. Proto nejsou srovnatelné ani indexy korelace, vypočtené při použití různých linearizujících funkcí. n

Tento problém se týká rovněž

∑ [ϕ ( y ) − ϕ ( y ′ )] = 0 , které je splněno pouze pro transformoi =1

i

i

vané hodnoty vysvětlované proměnné. Pokud vypočtenou funkci podrobíme inverzní transformaci, kladné a záporné odchylky pozorovaných hodnot kolem regresní funkce se již nekompenzují.

28

Příklad 4.8 Datový soubor, v němž vysvětlovaná proměnná je exponenciální funkcí vysvětlující proměnné, tj. y ′ = b0 b1x , byl po logaritmické transformaci vyrovnán funkcí log y ′ = 0,2312 + 0,3144 x . V tomto n

případě je

∑ (log y i =1

i

− log y i′ ) = 0 a tedy log y = log y ′ .

Po odlogaritmování je exponenciální funkce y ′ = 1,7029 ⋅ 2,0625 x , přičemž ovšem zatímco y ′ = 41,86 a

n

∑(y i =1

i

y = 41,17 ,

− y i′ ) = −4,15 . Rovnice rozkladu rozptylu pro exponenciální funkci ne-

platí, protože 10926,83 ≠ 11947,28 + 37,28 (index korelace by byl dokonce větší než jedna), zatímco pro logaritmovanou funkci je I log y . x = 0,996 . Vzhledem k úskalím linearizujících transformací bylo vyvinuto několik algoritmů pro přímou aplikaci metody nejmenších čtverců na funkce, které nejsou lineární v parametrech. Jedním z těchto algoritmů lze např. stanovit y ′ = 1,9393 ⋅ 2,0047 x , I yx = 0,999 . Tyto algoritmy nelineární regrese ovšem přesahují rámec tohoto textu a nebudeme se jimi zabývat.




K příkladu 4.8 máte k dispozici data

6  1 2 3 4 5 3 9 14 33 62 126 . Na těchto datech prověříme  

všechny výpočty z příkladu (samozřejmě s výjimkou nelineární regrese). (4–3)

V technické praxi se využívá rovněž řada funkcí, které mají často složité rovnice s mnoha parametry, které nelze žádným způsobem linearizovat. Rovněž v tomto případě se využívají (v příkladu 4.8 zmíněné) procedury nelineární regrese, které jsou však výpočetně mnohem náročnější.

Σ

1. Další úlohou o měření závislosti je regresní úloha, kdy pouze vysvětlovaná proměnná je náhodnou veličinou, zatímco vysvětlující proměnná je řízena experimentátorem. Charakteristikou intenzity závislosti je v tomto případě index korelace a průběh závislosti charakterizuje spojitá regresní funkce. 2. Stěžejní metodou řešení regresních úloh je metoda nejmenších čtverců, pomocí níž určujeme sadu parametrů regresní funkce, která minimalizuje hodnotu reziduálního součtu čtverců. 3. Metodu nejmenších čtverců jsme ukázali jednak pro funkce lineární v parametrech, jednak pro funkce, pro něž existuje vhodná linearizující transformace. 4. Metoda nejmenších čtverců umožňuje rozložit součet čtverců pozorovaných hodnot vysvětlované proměnné na složku objasněnou regresí a reziduum. Na tomto principu je založen index korelace. 5. Závěrem je třeba zdůraznit, že jsme se nezabývali ani měřením závislostí více než dvou proměnných, ani případem, kdy v regresní úloze figuruje regresní funkce, nelineární v parametrech, pro níž neexistuje linearizující transformace.

29



(4–1)

log y ′ = log b0 + b1 log x, y ′ 2 = b0 + b1 x 1 =0 x yi 1 1 − b0 ∑ − b1 ∑ 2 = 0 xi xi xi

∑y (4–2)

∑

i

b1 ∑

− nb0 −

(4–3) Výsledky jsou součástí příkladu 4.8. Stačí je pouze zkontrolovat.




1.

Pokuste se pokud možno co nejpřesněji specifikovat a porovnat regresní a korelační úlohu.

2.

Sestavte soustavy normálních rovnic pro regresní funkce (a) (b)

3. (a)

(b)

Upřesněte okolnosti, za kterých bude n

n

i =1 n

i =1 n

∑ ( yi − y) 2 > ∑ ( yi′ − y ) 2 , ∑ ( y ′ − y) i =1 n

(c)

5. 6.

2

i

∑ ( y′ − y) i =1

4.

2 y ′ = b0 + b1 log x, (c) y ′ = (b0 + b1 x ) .

y ′ = b0 + b1 x ,

i

= ∑ ( y i − y i′ ) 2 , i =1

2

= 0 . Tam, kde je to možné, určete i index korelace.

Patří index korelace mezi symetrické nebo asymetrické charakteristiky závislosti? Kolik procent variability chvění se podařilo objasnit závislostí na otáčkách turbíny v příkladu 4.4? Jak se nazývá příslušná charakteristika? Porovnejte výsledek adjustace s indexem korelace vypočteným na základě vyrovnání údajů v úloze 4.4 polynomem 3. stupně

y ′ = 0,9223 − 1,9531x + 1,2783x 2 + 0,4777 x 3 pro I yx = 0,9497 .

100 − 100 I yx2 ?

7.

Jak nazveme a vysvětlíme rozdíl

8.

Při zkoušení traktoru byla experimentátorem stanovena rychlost ( xi v km.hod-1) jako vysvětlující proměnná a měřena síla na hnací nápravě ( y i v kN) jako vysvětlovaná proměnná. Zjištěná data jsou v tabulce. Ke každé hodnotě vysvětlující proměnné jsou naměřeny shodně tři hodnoty vysvětlované proměnné (a tedy n = 18 ).

xi

2

4

6

8 10 12

65 49 44 39 36 35

y i 63 51 43 37 39 36 64 48 42 38 37 34 Ze zkušeností je známo, že vhodnou funkcí je lomená funkce y ′ = b0 +

b1 . Řešte x

úlohu po vzoru příkladů 4.4 až 4.7. Soustavu normálních rovnic pro tuto funkci jste určovali ve (4–2).

30