MASARYKOVA UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA ÚSTAV MATEMATIKY A STATISTIKY DIPLOMOVÁ PRÁCE VYHLAZOVÁNÍ A REGRESE

MASARYKOVA UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA ÚSTAV MATEMATIKY A STATISTIKY

DIPLOMOVÁ PRÁCE

VYHLAZOVÁNÍ A REGRESE

BRNO 2009

MARCELA HENZLOVÁ

Anotace Diplomová práce „Vyhlazování a regreseÿ se zabývá rozborem dvou přístupů k neparametrické regresi a neparametrickému vyhlazování, konkrétně jádrových odhadů a vyhlazovacích splajnů. Obě tyto metody jsou v práci zasazeny do kontextu regresní analýzy a čtenář je seznámen s jejich základními vlastnostmi a parametry, které ovlivňují fungování těchto metod. Dále jsou odvozeny asymptoticky optimální hodnoty těchto parametrů a také algoritmické přístupy k hledání těchto optimálních hodnot. Praktická část práce se věnuje implementaci zkoumaných metod v jazyce MATLAB a srovnáním jejich výsledků s užitím simulovaných dat.

Annotation Diploma thesis „Smoothing and regressionÿ is focused on two approaches to nonparametric regression, particularly kernel estimates and smoothing splines. Both of these methods are introduced and set up within the context of regression analysis. The reader is also introduced to their basic properties and the parameters, which affect their performance and output. Asymptotically optimal values of these parameters are derived and algorithms for finding these values are also presented. Practicall part of this thesis focuses on implementation of these methods in MATLAB language and comparing their performance using simulated data.

Poděkování Chtěla bych tímto poděkovat vedoucí své diplomové práce prof. RNDr. Ivance Horové, CSc. z Ústavu matematiky a statistiky PřF MU v Brně za cenné rady a připomínky při psaní práce, za čas strávený při jejím čtení a za trpělivost.

Prohlášení Prohlašuji, že jsem zadanou diplomovou práci vypracovala somostatně pod vedením prof. RNDr. Ivanky Horové, CSc. a veškerou použitou literaturu uvedla v seznamu literatury.

V Brně dne

................. Marcela Henzlová

Obsah Úvod

6

1 Regresní model

7

2 Jádrové odhady 2.1 Úvodní pojmy . . . . . . . . . . . . . . . . . . 2.2 Hledání optimální šířky vyhlazovacího okna h 2.3 Asymptotické vlastnosti jádrových odhadů . . 2.4 Hledání optimálního jádra K . . . . . . . . . 2.5 Hledání optimálního řádu jádra k . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

9 9 16 21 31 39

3 Vyhlazovací splajny 3.1 Úvodní pojmy . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Hledání optimálního parametru vyhlazení λ a parametru k 3.3 Hledání báze prostoru N S 2k (x1 , . . . , xn ) . . . . . . . . . . 3.4 Asymptotické vlastnosti vyhlazovacích splajnů . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

45 45 56 60 65

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

4 Simulační studie

70

Závěr

78

5

Úvod Regresní analýza tvoří nedílnou součást matematické statistiky, která nachází v praxi četná využití. Obory, pro něž je nutností zpracovávat množství dat, jako je například analýza finančních trhů, predikce časových řad, průzkumové studie, nebo různé technické obory, by si bez metod regresní analýzy jen těžko věděly rady. Regresní analýza představuje nástroj ke zkoumání dat, na které nahlížíme jako na realizace náhodné veličiny. Umožňuje nám v datech hledat trendy, různé periodicity, zkoumat vlastnosti a chování procesů, které tato data vygenerovaly, a také jejich chování předvídat. Základním přístupem regresní analýzy je předpoklad, že naměřené nebo jinak získané hodnoty náhodné veličiny lze rozdělit na deterministickou část a aditivní nekorelovanou chybu střední hodnoty nula. Tento předpoklad lze skutečně v mnoha případech v praxi uplatnit. Právě deterministickou složku jsme schopni metodami regresní analýzy odhadovat, zkoumat a modelovat. Z pohledu základního dělení rozlišujeme regresní analýzu parametrickou a neparametrickou. Parametrická regrese zkoumá případy, kdy je charakter regresní funkce dopředu znám a její tvar ovlivňuje předem daná skupina parametrů. Příkladem budiž polynomiální regrese, jejíž cílem je odhadnout koeficienty polynomu, který regresní funkci modeluje. Není-li však charakter regresní funkce dopředu znám, používáme metody neparametrické regrese, jíž je věnována tato práce. Konkrétně se budeme zabývat jádrovými odhady a vyhlazovacími splajny. V první kapitole si zopakujeme několik základních pojmů z oblasti regresní analýzy. Druhá kapitola se bude věnovat jádrovým odhadům, konkrétně jejich definicí a představením základních vlastností. Dále se budeme věnovat rozborem parametrů, které kvalitu jádrových metod ovlivňují a hledání jejich optimálních hodnot. Třetí kapitola se zabývá podobným způsobem vyhlazovacími splajny. Ve čtvrté kapitole budou srovnány výsledky obou metod na simulovaných datech za využití asymptotických vlastností obou metod. Tato kapitola se opírá o implementaci obou metod v jazyce MATLAB, která je k dispozici na přiloženém CD. 6

Kapitola 1 Regresní model Mějme regresní model s pevným plánem: yi = m(xi ) + εi,

i = 1, . . . , n,

(1.0.1)

kde y¯ = (y1 , . . . , yn )T je vektor závislých proměnných, xi body plánu takové, že platí 0 ≤ x1 < x2 < · · · < xn ≤ 1, a ε¯ = (ε1 , . . . , εn )T vektor chyb, přičemž předpokládáme, že εi , i = 1, . . . , n jsou nezávislé, stejně rozdělené a že platí E(εi) = 0, var(εi ) = σ 2 > 0, i = 1, . . . , n. Proces hledání vhodné aproximace m b neznámé funkce m bývá označován jako vyhlazování. Přístupy k problematice aproximace funkce m na základě uvedeného regresního modelu lze rozdělit do dvou kategorií. 1. Parametrické přístupy: Hledání aproximace neznámé funkce m se provádí v třídě funkcí zvoleného tvaru F (x; p1 , . . . , pm ), závislých na parametrech p1 , . . . , pm z nějaké množiny M. K výběru vhodné funkce ze zvolené třídy se často používá metoda nejmenších čtverců. To znamená, že se hledají takové hodnoty parametrů p1 , . . . , pm ∈ M, pro něž je výraz n X

(yi − F (xi ; p1 , . . . , pm ))2

i=1

minimální. Tuto metodu používáme, známe-li tvar aproximované funkce. 2. Neparametrické přístupy: Nesprávnou volbou třídy funkcí F (x; p1 , . . . , pm ) můžeme při parametrickém postupu obdržet velmi špatné výsledky. V takových případech je lepší volit některý z neparametrických přístupů, kdy se využívá pouze obecných vlastností funkce m, například její hladkosti. Mezi nejběžnější postupy neparametrické regrese patří jádrové odhady či splajnové vyhlazování, kterým se budeme podrobněji věnovat v následujících kapitolách. 7

8

KAPITOLA 1. REGRESNÍ MODEL

Vhodným nástrojem pro posouzení kvality odhadu je střední kvadratická chyba. 1.0.1 Definice. Nechť m(x) b je odhad funkce m(x). Střední kvadratická chyba odhadu m(x) b je definována vztahem Věta 1.0.2. Platí

2 MSE(m(x)) b = E(m(x) − m(x)) b .

MSE(m(x)) b = Em b 2 (x) − E 2 m(x) b + (E m(x) b − m(x))2 . | {z } | {z } rozptyl m(x) b

Důkaz.

2 (vychýlení m(x)) b

2 MSE(m(x)) b = E(m(x) − m(x)) b = E(m2 (x) − 2m(x)m(x) b +m b 2 (x)) = = m2 (x) − 2m(x)E m(x) b + Em b 2 (x) + E 2 m(x) b − E 2 m(x) b = 2 2 2 2 = Em b (x) − E m(x) b + m (x) − 2m(x)E m(x) b + E m(x) b = 2 2 2 = Em b (x) − E m(x) b + (E m(x) b − m(x)) .

Kapitola 2 Jádrové odhady 2.1

Úvodní pojmy

V této kapitole se budeme zabývat jádrovými odhady. Definujme tedy nyní jádro. 2.1.1 Definice. Nechť reálná funkce K definovaná na R splňuje tyto vlastnosti: 1. funkce K splňuje Lipschitzovu podmínku na intervalu [−1, 1], tj. pro každé x, y ∈ [−1, 1] platí kde L je konstanta,

|K(x) − K(y)| ≤ L|x − y|,

L > 0,

2. pro nosič funkce K platí supp(K) = [−1, 1], 3. pro daná celá nezáporná čísla ν, k mající stejnou paritu a 0 ≤ ν < k platí   pro 0 ≤ j < k, j 6= ν Z1 0 j ν x K(x) dx = (−1) ν! pro j = ν .   −1 βk 6= 0 pro j = k

Pak funkci K nazýváme jádrem řádu k a třídu všech takových jader značíme Kνk . Poznámka. Třetí vlastnost z definice 2.1.1 nazýváme momentové podmínky. Poznámka. Označíme-li jádro symbolem Kh , rozumíme tím, že 1 x , kde h ∈ R, h > 0, Kh (x) = K h h

přičemž jádru K přísluší nosič [−1, 1], zatímco jádru Kh přísluší nosič [−h, h]. 9

10

KAPITOLA 2. JÁDROVÉ ODHADY

Příklady jader (viz. Obr. 2.1): 1 2

1. Obdélníkové jádro: K(x) =

pro x ∈ [−1, 1]

3 2. Epanečnikovo jádro: K(x) = (1 − x2 ) 4 15 (1 − x2 )2 16

3. Kvartické jádro: K(x) =

pro x ∈ [−1, 1]

Epanecnikovo jadro

kvarticke jadro

2

2

2

1.5

1.5

1.5

1

1

0.5

0.5

y

1

0.5

y

y

obdelnikove jadro

pro x ∈ [−1, 1]

0

0

0

−0.5

−0.5

−0.5

−1 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(a) obdélníkové jádro

−1 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−1 −2

(b) Epanečnikovo jádro

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(c) kvartické jádro

Obr. 2.1: Druhy jader Pro odhad regresní funkce používáme jader třídy K0k , pro odhad ν-té derivace regresní funkce jsou vhodné jádra třídy Kνk . Nyní definujme některé pojmy a dokažme pomocná tvrzení. 2.1.2 Definice. Nechť β¯ = (β0 , . . . , βp )T a a ¯ = (a0 , . . . , ap )T jsou vektory takové, ¯ je skalární funkce vektoru β. ¯ Derivací funkce g(β) ¯ že platí a ¯, β¯ ∈ Rp+1. Nechť g(β) podle vektoru β¯ rozumíme vektor ¯ ∂g(β) = ∂ β¯

¯ ¯ T ∂g(β) ∂g(β) ,..., . ∂β0 ∂βp

Lemma 2.1.3. Platí tyto vztahy: 1.

∂¯ aT β¯ ∂ β¯T a ¯ = =a ¯, ¯ ¯ ∂β ∂β

∂ β¯T Aβ¯ ¯ kde A je symetrická matice typu (p + 1) × (p + 1). 2. = 2Aβ, ∂ β¯

11


¯ =a Důkaz. 1. Nechť g(β) ¯T β¯ = a0 β0 + a1 β1 + · · · + ap βp = β¯T a¯, pak ¯ ∂g(β) = (a0 , a1 , . . . , ap )T = a ¯. ¯ ∂β ¯ = β¯T Aβ¯ = Pp aii β 2 + 2 Pp 2. Nechť g(β) i i,j=0 aij βi βj , pak pro i = 0, . . . , p platí i=0 i6=j

p p X X ¯ ∂g(β) =2 aii βi + 2 aij βj , ∂βi i=0 i,j=0 i6=j

tedy

¯ ∂g(β) ¯ = 2Aβ. ∂ β¯

Vraťme se k modelu (1.0.1). Neznámou funkci m(x) lze aproximovat pomocí polynomu β0 + β1 (xi − x) + · · · + βp (xi − x)p stupně p na lokálním intervalu šířky h, tedy p X m(x) = βj (xi − x)j , pro x ∈ [xi + h, xi − h]. j=0

Odhad m(x; b p, h) = βb0 této funkce nalezneme pomocí vážené metody nejmenších ¯ kde čtverců, tedy hledáme argminβ¯ S(β), ¯ = S(β)

n X

[yi − β0 − β1 (xi − x) − · · · − βp (xi − x)p ]2 Kh (xi − x),

i=1

přičemž roli vah plní jádro Kh . Předpokládejme, že K je nezáporné jádro. Označme:  . . . (x1 − x)p . . . (x2 − x)p   , .. ..  . . p . . . (xn − x)  0  0  . ..  . . . . Kh (xn − x)

   1 (x1 − x) β0 1 (x2 − x) β1     β¯ =  ..  , X =  .. .. . . . 1 (xn − x) βp  Kh (x1 − x) 0 ...  0 Kh (x2 − x) . . .  W = .. .. ..  . . .

  y1  y2    y¯ =  ..  , . yn

0

0


12

Můžeme tedy zapsat váženou metodu nejmenších čtverců maticově takto: ¯ = (¯ ¯ T W (¯ ¯ = y¯T W y¯ − y¯T W X β¯ − (X β) ¯ T W y¯ + (X β) ¯ T W X β. ¯ S(β) y − X β) y − X β) ¯ přes β, ¯ položíme její derivaci podle proAbychom minimalizovali funkci S(β) měnné β¯ rovnu 0 a s využitím lemmatu 2.1.3 získáme 0=

¯ ¯ T W y¯ ∂(X β) ¯ T W X β¯ ∂S(β) ∂ y¯T W y¯ ∂ y¯T W X β¯ ∂(X β) = − + = − ∂ β¯ ∂ β¯ ∂ β¯ ∂ β¯ ∂ β¯ | {z } =0

¯ = −2X W y¯ + 2X T W X β. T

b¯ parametru β, ¯ tj. Nyní z rovnice X T W X β¯ = X T W y¯ získáme odhad β b¯ = (X T W X)−1 X T W y¯. β

Věta 2.1.4. Nechť X T W X je regulární matice a nechť e¯1 = (1, 0, . . . , 0)T . Pak m(x; b p, h) = βb0 = e¯T1 (X T W X)−1 X T W y¯. Důkaz. Víme, že platí:

Tedy

  βb0 βb  1 b¯ =  β  ..  = (X T W X)−1 X T W y¯. . βbp

    βb0 βb0 βb  βb   1 T T −1 T T  1 e¯1 (X W X) X W y¯ = e¯1  .  = (1, 0, . . . , 0)  .  = βb0 . .  ..  . βbp βbp

Důsledek 2.1.5. Je-li X T W X pozitivně definitní matice, pak je minimum jediné. Důkaz. Předpokládejme, že existují vzájemně různé β¯• , β¯◦ , jenž obě minimalizují ¯ tj. S(β), X T W X β¯• = X T W y¯ a X T W X β¯◦ = X T W y¯. Odečteme-li druhou rovnici od první, dostaneme, že (X T W X)(β¯• − β¯◦ ) = 0, a jelikož je podle předpokladu X T W X pozitivně definitní matice, musí nutně platit β¯• ≡ β¯◦ . Nyní se seznámíme s nejznámějšími typy jádrových odhadů.

13


1. Nadaraya-Watsonovy odhady: Jsou speciálním případem, kdy p = 0, tedy m b N W (x; h) := m(x; b 0, h).    Kh (x1 − x) 0 ... 0 1     0 K (x − x) . . . 0 h 2   1 T X W X = (1, 1, . . . , 1)    ..  = .. .. .. ..   . . . . . 0 0 . . . Kh (xn − x) 1 n X = Kh (xi − x). i=1

1 . Dále platí i=1 Kh (xi − x)    y1 Kh (x1 − x) 0 ... 0     0 K (x − x) . . . 0 h 2   y2   T X W y¯ = (1, 1, . . . , 1)    ..  = .. .. .. ..  .   . . . . 0 0 . . . Kh (xn − x) yn n X = Kh (xi − x)yi . Inverzní matice pak má tvar (X T W X)−1 = Pn

i=1

Výsledný odhad je tedy Pn Kh (xi − x)yi m b N W (x; h) = Pi=1 . n i=1 Kh (xi − x)

Na Obr. 2.2 je pro lepší představu znázorněna konstrukce Nadaraya-Watsonova odhadu. 2.5

2

1.5

1.5

1.5

1

1

1 y

2

y

2.5

2

y

2.5

0.5

0.5

0

0

0

−0.5

−0.5

−0.5

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(a) jádra vynásobená hodnotami bodů pozorování

−1

0.5

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

(b) součet jader

0.8

0.9

1

−1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(c) jádrový odhad

Obr. 2.2: Konstrukce Nadaraya-Watsonova odhadu

0.9

1

14


2. Lokální lineární odhady: Jsou speciálním případem, kdy p = 1, tedy m b LL (x; h) := m(x; b 1, h). XT W X =

   Kh (x1 − x) . . . 0 1 x1 − x 1 ... 1   ..  .. .. ..  = .. =   . . . . .  x1 − x . . . xn − x 0 . . . Kh (xn − x) 1 xn − x   1 x1 − x Kh (x1 − x) ... Kh (xn − x)  .. ..  = = . .  (x1 − x)Kh (x1 − x) . . . (xn − x)Kh (xn − x) 1 xn − x Pn Pn K (x − x) (x − x)K (x − x) h i i h i i=1 i=1 Pn = Pn . 2 i=1 (xi − x)Kh (xi − x) i=1 (xi − x) Kh (xi − x)

Definujme „pomocnéÿ funkce: n

1X (xi − x)r Kh (xi − x), sbr (x; h) = n i=1

Pak platí

r = 0, 1, 2, . . .

sb0 (x; h) sb1 (x; h) X WX = n , s1 (x; h) sb2 (x; h) b 1 sb2 (x; h) −b s1 (x; h) T −1 , (X W X) = s1 (x; h) sb0 (x; h) n(b s2 (x; h)b s0 (x; h) − sb21 (x; h)) −b T

kde

sb2 (x; h)b s0 (x; h) − sb21 (x; h) =

1 X (xi − xj )2 Kh (xi − x)Kh (xj − x). 2 n i,j i6=j

Dále platí X T W y¯ =

1 ... 1 x1 − x . . . xn − x

 Kh (x1 − x) . . .  .. ..  . . 0

0 .. .

  y1   ..   .  =

yn . . . Kh (xn − x)   y1 Kh (x1 − x) ... Kh (xn − x)  ..  = .= (x1 − x)Kh (x1 − x) . . . (xn − x)Kh (xn − x) yn Pn y K (x − x) = Pni=1 i h i . i=1 yi (xi − x)Kh (xi − x)

15


Výsledný odhad je n

1 X (b s2 (x; h) − sb1 (x; h)(xi − x)) Kh (xi − x)yi m b LL (x; h) = . n i=1 sb2 (x; h)b s0 (x; h) − sb1 (x; h)2 3. Gasser-Müllerovy odhady: m b GM (x; h) =

kde s0 = 0,

si =

n X

(xi + xi+1 ) 2

i=1

yi

Zsi

Kh (t − x) dt,

si−1

pro i = 1, . . . , n − 1,

sn = 1.

Pro odhad derivace regresní funkce je vhodné použít Gasser-Müllerova odhadu ve tvaru: Zsi n 1 X (ν) yi m b GM (x; h) = ν Kh (t − x) dt, kde K ∈ Kνk . h i=1 si−1

16


2.2

Hledání optimální šířky vyhlazovacího okna h

Volba šířky vyhlazovacího okna má výrazný vliv na kvalitu odhadu. Zvolíme-li h příliš malé, odhad bude méně vychýlen, avšak na úkor velké variability. Říkáme, že výsledný odhad je podhlazený. Při příliš velkém h, bude rozptyl odhadu malý, což však bude mít za následek nárust jeho vychýlení od skutečné hodnoty. Říkáme, že výsledný odhad je přehlazený. Pro volbu optimální šířky vyhlazovacího okna se obvykle využívá metody křížového ověřování, kterou popíšeme v následující podkapitole. Odhady regresních funkcí lze zapsat ve tvaru m(x; b h) =

n X

wi(x; h)yi ,

i=1

kde wi (x; h) je tzv. váhová funkce.

Příslušná váhová funkce je u Nadaraya-Watsonových odhadů Pn Kh (xi − x)yi m b N W (x; h) = Pi=1 n i=1 Kh (xi − x) tvaru

Kh (xi − x) wiN W (x; h) = Pn , i=1 Kh (xi − x)

u lokalních lineárních odhadů n

tvaru

m b LL (x; h) =

1 X (b s2 (x; h) − sb1 (x; h)(xi − x)) Kh (xi − x)yi n i=1 s2 (x; h)b b s0 (x; h) − sb1 (x; h)2

wiLL(x; h) =

(b s2 (x; h) − sb1 (x; h)(xi − x)) Kh (xi − x) n(b s2 (x; h)b s0 (x; h) − sb1 (x; h)2 )

a u Gasser-Müllerových odhadů

tvaru

m b GM (x; h) =

n X i=1

wiGM (x; h) =

yi

Zsi

Kh (t − x) dt,

si−1

Zsi

si−1

Kh (t − x) dt,

17


kde s0 = 0,

si =

(xi + xi+1 ) 2

pro i = 1, . . . , n − 1,

sn = 1.

Odhad funkce m v bodě xj , j = 1, . . . , n je roven m(x b j ; h) =

tedy

n X

wi (xj ; h)yi,

j = 1, . . . , n,

i=1

m(x b 1 ; h) = w1 (x1 ; h)y1 + w2 (x1 ; h)y2 + · · · + wn (x1 ; h)yn m(x b 2 ; h) = w1 (x2 ; h)y1 + w2 (x2 ; h)y2 + · · · + wn (x2 ; h)yn .. . m(x b n ; h) = w1 (xn ; h)y1 + w2 (xn ; h)y2 + · · · + wn (xn ; h)yn .

To lze maticově zapsat takto:    w1 (x1 ; h) w2 (x1 ; h) m(x b 1 ; h)    m(x  b 2 ; h)   w1 (x2 ; h) w2 (x2 ; h) m(¯ b x; h) =  = .. .. ..    . . . w1 (xn ; h) w2 (xn ; h) m(x b n ; h)

  y1 . . . wn (x1 ; h)   . . . wn (x2 ; h)   y2     ..  = Sh y¯, .. ..   . . . yn . . . wn (xn ; h)

kde Sh = {sijh }ni,j=1 je vyhlazovací (klobouková) matice (viz. Obr. 2.3), pro kterou sijh = wj (xi ; h).

vyhlazovaci matice jadra, h = 0.2

vyhlazovaci matice jadra, h = 0.6

0.35

0.14

0.3

0.12 0.1 hodnoty matice

hodnoty matice

0.25 0.2 0.15 0.1 0.05

0

0 0 10 10

0.06 0.04 0.02

0

0 0

5 5

0.08

5 10

5 10

15

15

15

15 20

sloupce matice Sh

20

radky matice S

20

h

sloupce matice Sh

(a) h = 0, 2

(b) h = 0, 6

Obr. 2.3: Vyhlazovací matice

20

radky matice S

h

18


Definujme nyní odhad funkce m(x) v bodě xi bez použití bodu xi : m b −i (xi ; h) =

n X j=1 j6=i

sijh yj , 1 − siih

kde

n X j=1 j6=i

sijh = 1. 1 − siih

2.2.1 Definice. Funkce křížového ověřování je definovaná vztahem n

CV (h) =

1X (yi − m b −i (xi ; h))2 . n i=1

Věta 2.2.2. Nechť m(x; b h) je odhad m(x) a nechť platí, že m b −i (xi ; h) =

n X j=1 j6=i

sijh yj + siih m b −i (xi ; h).

Předpokládejme, že siih 6= 1 pro i = 1, . . . , n. Pak funkci CV (h) můžeme vyjádřit ve tvaru 2 n 1 X yi − m(x b i ; h) CV (h) = . (2.2.1) n i=1 1 − siih

Důkaz. Vyjdeme nejprve z odhadu funkce m v bodě xi s vynecháním tohoto bodu, přičemž požadujeme, aby konstanta byla zachována. Toto vyjádříme následujícím způsobem: m b −i (xi ; h) =

n X j=1 j6=i

sijh yj , 1 − siih

přičemž platí n X j=1 j6=i

n sijh 1 X 1 = sijh = (1 − siih ) = 1. 1 − siih 1 − siih j=1 1 − siih j6=i

Malou úpravou dostaneme z rovnosti (2.2.2): (1 − siih )m b −i (xi ; h) = m b −i (xi ; h) =

n X

sijh yj ,

j=1 j6=i

sijh yj + siih m b −i (xi ; h).

j=1 j6=i n X

(2.2.2)

19


Nyní

yi − m b −i (xi ; h) = yi − =

n X j=1 j6=i

1 1 − siih





n X   sijh 1  = yj = y − s y − s y i ii i ij j h h  1 − siih 1 − siih  j=1

yi −

n X

sijh yj

j=1

Odtud plyne, že

n

1X CV (h) = n i=1

!

j6=i

=

1 (yi − m(x b i ; h)) . 1 − siih

yi − m(x b i ; h) 1 − siih

2

.

Chceme-li odhadnout optimální šířku vyhlazovacího okna, hledáme argument mi1 1 nima funkce CV (h) na množině Hn = [ak n− 2k+1 , bk n− 2k+1 ], kde 0 < ak < bk < ∞, tj. b hCV,0,k = argmin CV (h), K ∈ K0,k . (2.2.3) h∈Hn

Poznámka. Vztah (2.2.1) můžeme zobecnit na tvar n

1X GCV (h) = n i=1

yi − m(x b i ; h) 1 − tr(Sh )/n

2

,

P kde tr(Sh ) je stopa matice Sh = {sijh }ni,j=1, tj. platí tr(Sh ) = ni=1 siih . Metoda se pak v tomto případě nazývá zobecněná metoda křížového ověřování. Odhad šířky vyhlazovacího okna v případě odhadu 1. derivace funkce m dostaneme ze vztahu n−1 2 1 X (1) (1) (1) CV (1) (h) = ∆i − m b −{i,i+1} (xi ; h) , n − 1 i=1 kde

(0)

xi = xi , (0)

∆i = yi , (1)

(1)

1 (0) (0) xi+1 − xi , 2 (0) (0) ∆ − ∆i = i+1 , (0) (0) xi+1 − xi

xi = (1)

∆i (1)

(1)

pro i = 1, . . . , n − 1 a m b −{i,i+1} (xi ; h) je Gasser-Müllerův odhad v bodě xi konstruovaný na datech (x1 , y1), . . . , (xi−1 , yi−1), (xi+2 , yi+2), . . . , (xn , yn ).

20


V tomto případě symbolem b hCV (1) ,1,k označme číslo, které minimalizuje CV (1) (h) (1) (1) na vhodné množině Hn (zpravidla Hn = [ n1 , 2]), tj. b hCV (1) ,1,k = argmin CV (1) (h),

K ∈ K1,k .

(2.2.4)

(1)

h∈Hn

Číslo b hCV (1) ,1,k je tedy odhad optimálního vyhlazovacího parametru pro odhad 1. derivace funkce m.

21


2.3

Asymptotické vlastnosti jádrových odhadů

Pro zjednodušení se v této podkapitole zabývejme jádry 2. řádu. Připomeňme, že jádrem 2. řádu nazveme reálnou funkci K definovanou na R splňující tyto podmínky: 1. funkce K splňuje Lipschitzovu podmínku na intervalu [−1, 1], tj. pro každé x, y ∈ [−1, 1] platí |K(x) − K(y)| ≤ L|x − y|, kde L je konstanta a L > 0, 2. pro nosič funkce K platí supp(K) = [−1, 1], 3.

R1

K(x) dx = 1,

−1

R1

xK(x) dx = 0 a

−1

R1

x2 K(x) dx = β2 6= 0.

−1

Zaveďme nyní označení, které budeme nadále používat. Předpokládejme, že an a bn jsou posloupnosti reálných čísel. Pak 1. an = O(bn ), pro n → ∞ tehdy a jen tehdy, jestliže lim supn→∞ | abnn | < ∞, 2. an = o(bn ), pro n → ∞ tehdy a jen tehdy, jestliže lim supn→∞ | abnn | = 0. Následující lemma a důsledek je užitečný pro důkaz věty o tvaru vychýlení a rozptylu. Lemma 2.3.1. Jestliže 0 < h < x < 1 − h a xi = ni , i = 1, . . . , n, pak pro všechna r ∈ N0 platí: 1. r

sbr (x; h) = h

kde

Z1

xr K(x) dx + O(n−1 ),

−1

n

1X sr (x; h) = b (xi − x)r Kh (xi − x), n i=1

2. 1 h

Z1

−1

3. r−1

h

Z1

−1

n

1X 2 K (x) dx + O(n ) = K (xi − x), n i=1 h −1

2

n

1X x K (x) dx + O(n ) = (xi − x)r Kh2 (xi − x). n i=1 r

2

−1

(2.3.1)

22


Důkaz. Viz. [2]. Důsledek 2.3.2. s0 (x; h) = b

Z1

K(x) dx + O(n−1 ) = 1 + O(n−1 ),

−1

s1 (x; h) = h b

Z1

xK(x) dx + O(n−1 ) = O(n−1 ),

−1 2

s2 (x; h) = h b

3

s3 (x; h) = h b

Z1

x2 K(x) dx + O(n−1 ) = h2 β2 + O(n−1 ),

Z1

x3 K(x) dx + O(n−1 ) = O(n−1 ).

−1

−1

.. .

Důkaz. Plyne z podmínky 1. předchozího lemmatu a z definice jádra 2. řádu. Věta 2.3.3 (O tvaru vychýlení a rozptylu.). Nechť jsou splněny následující předpoklady: 1. Šířka vyhlazovacího okna h = hn splňuje limn→∞ h = 0 a limn→∞ nh = ∞. 2. Pro bod odhadu uvnitř intervalu [0, 1] platí, že h < x < 1 − h pro všechna n ≥ n0 , kde n0 je pevné. 3. Předpokládejme, že pro body plánu platí: xi = ni , i = 1, . . . , n. 4. Funkce m splňuje m ∈ C 2 [0, 1], tj. m má spojité derivace až do řádu 2 včetně. 5. Funkce K je jádro řádu 2. Pak pro odhad m b LL (x; h) := m(x, b 1, h) v bodě x platí a

1 vychýlení m b LL (x; h) = h2 β2 m′′ (x) + o(h2 ) + O(n−1 ) 2

σ 2 V (K) rozptyl m b LL (x; h) = + o((nh)−1 ), nh

kde

V (K) =

Z1

−1

K 2 (x) dx.

23


Důkaz. a) vychýlení odhadu: Víme, že Dokážeme, že

vychýlení m b LL (x; h) = E m b LL (x; h) − m(x).

1 Em b LL (x; h) − m(x) = h2 β2 m′′ (x)(1 + o(1)) + O(n−1 ). 2

Označme m(¯ x) = (m(x1 ), . . . , m(xn ))T . Podle věty 2.1.4 platí m b LL (x; h) = eT1 (X T W X)−1 X T W y¯.

Jelikož E y¯ = m(¯ x), pak

Em b LL (x; h) = eT1 (X T W X)−1 X T W E y¯ = eT1 (X T W X)−1X T W m(¯ x).

Taylorův rozvoj funkce m v bodě xi je ve tvaru:

1 m(xi ) = m(x) + (xi − x)m′ (x) + (xi − x)2 m′′ (x) + . . . . 2

Pomocí Taylorova rozvoje můžeme tedy maticově vyjádřit m(¯ x) takto:       m(x1 ) 1 x1 − x (x1 − x)2     ..  m(x) + 1 m′′ (x)  .. m(¯ x) =  ...  =  ...  . .  m′ (x) . 2 2 m(xn ) 1 xn − x (xn − x) | {z } X

Tedy

 (x1 − x)2 1 m(x)    .. Em b LL (x; h) = (1, 0)(X T W X)−1 X T W X + m′′ (x)   = . m′ (x) 2 2 (xn − x)   (x1 − x)2 1 m(x)   .. = (1, 0) + (1, 0)(X T W X)−1 X T W m′′ (x)  = . ′ m (x) 2 2 (xn − x)   (x1 − x)2 1   .. = m(x) + (1, 0)(X T W X)−1 X T W m′′ (x)  . . 2 2 (xn − x) 

Dále, jelikož



 (x1 − x)2 sb2 (x; h)   . T . X W , =n . sb3 (x; h) 2 (xn − x)



24


pak Em b LL (x; h) − m(x) =



 (x1 − x)2 1   .. = m′′ (x)(1, 0)(X T W X)−1X T W  = . 2 (xn − x)2 1 1 ′′ s2 (x; h) −b b s1 (x; h) s2 (x; h) b = m (x)(1, 0) = s1 (x; h) s0 (x; h) b s3 (x; h) b 2 sb2 (x; h)b s0 (x; h) − b s21 (x; h) −b m′′ (x) sb22 (x; h) − sb3 (x; h)b s1 (x; h) = . 2 sb2 (x; h)b s0 (x; h) − sb21 (x; h) Podle důsledku 2.3.2 platí 2 sb22 (x; h) − b s3 (x; h)b s1 (x; h) = h2 β2 + O(n−1 ) − O(n−1 )O(n−1 ), 2 sb2 (x; h)b s0 (x; h) − sb21 (x; h) = 1 + O(n−1 ) h2 β2 + O(n−1 ) − O(n−1 ) . Zanedbáme-li členy O(n−1 ), dostaneme Em b LL (x; h) − m(x) =

2

m′′ (x) (h2 β2 + O(n−1 )) m′′ (x) h2 β2 + O(n−1 ) = = 2 (1 + O(n−1 )) (h2 β2 + O(n−1 )) 2 1 + O(n−1 ) m′′ (x) h2 β2 (1 + O(n−1 )) − h2 β2 O(n−1 ) + O(n−1 ) = = 2 1 + O(n−1 ) m′′ (x) 2 O(n−1 ) O(n−1 ) = h β2 − h2 β2 + = 2 1 + O(n−1 ) 1 + O(n−1 ) 1 = h2 β2 m′′ (x) + o(h2 ) + O(n−1 ). 2 b) rozptyl odhadu: =

Víme, že rozptyl m b LL (x; h) = E(m b LL (x; h))2 −E 2 m b LL (x; h) = E(m b LL (x; h) − E m b LL (x; h))2 . Dokážeme, že

σ 2 V (K) E(m b LL (x; h) − E m b LL (x; h)) = +o((nh)−1 ), nh 2

kde V (K) =

Z1

K 2 (x) dx.

−1

Platí

m b LL (x; h) − E m b LL (x; h) = (1, 0)(X T W X)−1 X T W (¯ y − m(¯ x)) = = (1, 0)(X T W X)−1 X T W ε¯ = n X = wiLL (x; h)εi . i=1

25


Dále rozptyl m b LL (x; h) = E

=E

n X n X

n X

=

i=1

wi (x; h)wj (x; h)εi εj

i=1 j=1

n X

wi (x; h)εi

wi2 (x; h)Eε2i = σ 2

i=1

n X

!

!2

=

=

wi2 (x; h) =

i=1

2 n σ X (b s2 (x; h) − b s1 (x; h)(xi − x)) Kh (xi − x) = 2 = n i=1 sb2 (x; h)b s0 (x; h) − b s21 (x; h) 2

n X 2 σ2 sb2 (x; h)Kh2 (xi − x) − = 2 2 2 n (b s2 (x; h)b s0 (x; h) − sb1 (x; h)) i=1

− 2b s1 (x; h)b s2 (x; h)(xi − x)Kh2 (xi − x) + sb21 (x; h)(xi − x)2 Kh2 (xi − x) = " n X σ2 2 = 2 sb (x; h) Kh2 (xi − x) − n (b s2 (x; h)b s0 (x; h) − sb21 (x; h))2 2 i=1 − 2b s1 (x; h)b s2 (x; h)

n X i=1

(xi − x)Kh2 (xi − x) + sb21 (x; h)

n X

#

(xi − x)2 Kh2 (xi − x) .

i=1

S využitím vztahu (2.3.1), důsledku 2.3.2 a zanedbáme-li členy O(n−1 ), dostaneme rozptyl m b LL (x; h) =

" σ2 1 2 −1 = sb (x; h) V (K) + O(n ) − n(b s2 (x; h)b s0 (x; h) − sb21 (x; h))2 2 h  1  Z − 2b s1 (x; h)b s2 (x; h)  xK 2 (x) dx + O(n−1 ) + 

−1

+ sb21 (x; h) h

Z1

−1



x2 K 2 (x) dx + O(n−1 ) =

2 (h β2 + O(n−1 )) h1 V (K) + O(n−1 ) σ σ 2 h1 V (K) + O(n−1 ) = = = n [(1 + O(n−1 )) (h2 β2 + O(n−1 )) − (O(n−1 ))]2 n 1 + O(n−1 ) σ 2 h1 V (K) (1 + O(n−1 )) − h1 V (K)O(n−1 ) + O(n−1 ) = = n 1 + O(n−1 ) σ2 σ2 O(n−1 ) σ 2 O(n−1 ) σ 2 V (K) = V (K) − V (K) + = + o((nh)−1 ). −1 −1 nh nh 1 + O(n ) n 1 + O(n ) nh 2

2

26


Důsledek 2.3.4. Nechť jsou splněny předpoklady věty 2.3.3. Pro střední kvadratickou chybu odhadu m b LL (x; h) v bodě x platí 2 σ 2 V (K) 1 2 −1 ′′ 2 −1 MSE(m b LL (x; h)) = + o((nh) ) + h β2 m (x) + o(h ) + O(n ) . nh 2 Poznámka. Hlavní člen MSE se značí MSE, tedy 2 σ 2 V (K) 1 2 σ 2 V (K) 1 4 2 ′′ ′′ MSE(m b LL (x; h)) = + h β2 m (x) = + h β2 [m (x)]2 . nh 2 nh 4

Poznámka. Lze ukázat, že odhady m b N W (x; h), m b LL (x; h) a m b GM (x; h) jsou asymptoticky ekvivalentní, tj. MSE(m b N W (x; h)) = MSE(m b LL (x; h)) = MSE(m b GM (x; h)).

Věta 2.3.5. Nechť jsou splněny předpoklady věty 2.3.3. Pak asymptoticky optimální šířka vyhlazovacího okna je dána vztahem 2 1 σ V (K) 5 − 15 . hopt,2 = n β22 [m′′ (x)]2

Důkaz. Pro pevné x spočítáme derivaci střední kvadratické chyby odhadu podle h a položíme ji rovnu 0, tedy ∂MSE(m(x; b h)) σ 2 V (K) =− + h3 β22 [m′′ (x)]2 = 0. ∂h nh2 Vyjádříme-li nyní z poslední rovnosti h, dostaneme 2 1 σ V (K) 5 − 15 h=n . β22 [m′′ (x)]2 Důsledek 2.3.6.

4 1 5 4 MSE(m(x; b hopt,2 )) = n− 5 σ 2 V (K) 5 β22 [m′′ (x)]2 5 , 4 4

tj. asymptotická rychlost konvergence MSE → 0 je n− 5 . Důkaz. MSE(m(x; b hopt,2 )) =

2 15 !4 1 σ 2 V (K) 1 σ V (K) − = β22 [m′′ (x)]2 = 2 15 ! + 4 n 5 β 2 [m′′ (x)]2 2 σ V (K) 1 n n− 5 2 ′′ 2 β2 [m (x)] 4 1 1 4 4 1 4 = n− 5 σ 2 V (K) 5 β22 [m′′ (x)]2 5 + n− 5 σ 2 V (K) 5 β22 [m′′ (x)]2 5 = 4 1 45 2 ′′ 5 −4 2 = n 5 σ V (K) β2 [m (x)]2 5 . 4

27


Důsledek 2.3.7. Pro hopt,2 platí 1 (vychýlení m(x; b hopt,2 ))2 = . rozptyl m(x; b hopt,2 ) 4

Důkaz. Plyne z důkazu předchozí věty, neboť

MSE(m(x; b hopt,2 )) = 4 1 1 4 4 1 4 = n− 5 σ 2 V (K) 5 β22 [m′′ (x)]2 5 + n− 5 σ 2 V (K) 5 β22 [m′′ (x)]2 5 . | {z } |4 {z } 2 (vychýlení m(x;h b opt,2 ))

rozptyl m(x;h b opt,2 )

Střední kvadratická chyba určuje chybu odhadu v konkrétním bodě, je tedy pouze lokálního charakteru. Pokud chceme zjistit, jak se chová „celýÿ odhad, je vhodné použít integrální střední kvadratickou chybu. 2.3.8 Definice. Integrální střední kvadratická chyba odhadu m(x; b h) je definována vztahem Z1 MISE(m(h)) b = MSE(m(x; b h)) dx. 0

Poznámka. Hlavní člen MISE značíme MISE, obdobně jako u MSE. Věta 2.3.9. σ 2 V (K) 1 4 2 MISE(m(h)) b = + h β2 nh 4

Z1

[m′′ (x)]2 dx.

0

Důkaz. Plyne přímo z definice MISE.

Dosud jsme se zabývali jádry třídy K02 . Nyní se zabývejme obecněji jádry třídy K0k . Věta 2.3.10. Nechť K ∈ K0k , m ∈ C k [0, 1] a nechť h → 0, nh → ∞ pro n → ∞. Pak 2 σ 2 V (K) k 1 k (k) MSE(m(x; b h)) = + (−1) h βk m (x) , nh k! Z1 σ 2 V (K) βk2 2k MISE(m(h)) b = + 2h [m(k) (x)]2 dx. nh k! 0

28


Důkaz. Viz. [2]. Věta 2.3.11. Asymptoticky optimální hodnota vyhlazovacího parametru je dána vztahem 1   2k+1 hopt,k

  = 

 σ 2 V (K)   1  βk2 R (k) 2 2kn k!2 [m (x)] dx

.

0

Důkaz. Spočítáme derivaci integrální střední kvadratické chyby odhadu podle h a položíme ji rovnu 0, tedy σ 2 V (K) βk2 2k−1 ∂MISE(m(h)) b =− + 2k h ∂h nh2 k!2

Z1

[m(k) (x)]2 dx = 0.

0

Vyjádříme-li nyní z poslední rovnosti h, dostaneme 

Důsledek 2.3.12.

  h= 

1  2k+1

 σ 2 V (K)   1 R  βk2 (k) 2 2kn k!2 [m (x)] dx 0

2k

MISE(m(h b opt,k )) = n− 2k+1

kde

Dk =

.

1 + 2k (2k)

Z1 0

2k 2k+1

2k

(σ 2 V (K)) 2k+1 βk2 Dk

m(k) (x) k!

2

1 2k+1

,

dx, 2k

tj. asymptotická rychlost konvergence MISE → 0 je n− 2k+1 . Důkaz. σ 2 V (K) MISE(m(h b opt,k )) = + βk2 1 2k+1 2 σ V (K) n 2knβk2 Dk 1

σ 2 V (K) 2knβk2 Dk

2k 2k+1

Dk =

2k 2k 1 1 (2kn) 2k+1 2 1 σ 2 V (K) 2k+1 βk2 Dk 2k+1 = = σ V (K) 2k+1 βk2 Dk 2k+1 + 2k n (2kn) 2k+1 1 2k+1 2k 2k 1 + 2k 2 2k+1 β 2 D = n− 2k+1 . 2k (σ V (K)) k k (2k) 2k+1

29


Důsledek 2.3.13. Pro hopt,k platí (vychýlení m(x; b hopt,k ))2 1 = . rozptyl m(x; b hopt,k ) 2k

Důkaz. Plyne z důkazu předchozí věty, neboť: MISE(m(h b opt,k )) = 1

2k 1 2k+1 2k 1 2k+1 (2kn) 2k+1 2 1 2 2 = σ V (K) σ V (K) 2k+1 βk2 Dk 2k+1 + β D k 2k k n (2kn) 2k+1 | {z } | {z } 2 (vychýlení m(x;h b opt,k ))

rozptyl m(x;h b opt,k )

a

1

1 2k

(2kn) 2k+1

(2kn) 2k+1 1 = . : n 2k

Pro odhad derivace regresní funkce platí následující věty a důsledky. Věta 2.3.14. Nechť K ∈ Kνk , m ∈ C k [0, 1] a nechť h → 0, nh2ν+1 → ∞ pro n → ∞, 0 ≤ ν < k. Pak 2 σ 2 V (K) (ν) k 1 (k−ν) (k) + (−1) h βk m (x) , MSE(m b GM (x; h)) = nh2ν+1 k! Z1 2 2 σ V (K) β (ν) MISE(m b GM (h)) = + k2 h2(k−ν) [m(k) (x)]2 dx. nh2ν+1 k! 0

Důkaz. Viz. [2].

Věta 2.3.15. Asymptoticky optimální hodnota vyhlazovacího parametru je dána vztahem 1  2k+1    hopt,ν,k =  

 (2ν + 1)σ 2 V (K)   1  βk2 R (k) 2 2(k − ν)n k!2 [m (x)] dx

.

0

Důkaz. Spočítáme derivaci integrální střední kvadratické chyby odhadu podle h a položíme ji rovnu 0, tedy (ν)

∂MISE(m b GM (h)) (2ν + 1)σ 2 V (K) βk2 2(k−ν)−1 =− +2(k−ν) h ∂h nh2(ν+1) k!2

Z1 0

[m(k) (x)]2 dx = 0.

30


Vyjádříme-li nyní z poslední rovnosti h, dostaneme 

Důsledek 2.3.16.

  h= 

1  2k+1

 (2ν + 1)σ 2 V (K)   1  βk2 R (k) 2 2(k − ν)n k!2 [m (x)] dx

.

0

(ν)

kde

MISE(m b GM (hopt,ν,k )) = 2ν+1 2ν+1 2(k−ν) 2(k−ν) 2ν + 1 2(k − ν) 2k+1 2 − 2k+1 1+ (σ V (K)) 2k+1 βk2 Dk 2k+1 , =n 2(k − ν) 2ν + 1 Dk =

Z1

m(k) (x) k!

0

2

dx,

tj. asymptotická rychlost konvergence MISE → 0 je n−

2(k−ν) 2k+1

.

Důkaz. (ν)

=

MISE(m b GM (hopt,ν,k )) = σ 2 V (K)

βk2

(2ν + 1)σ 2 V (K) 2(k − ν)nβk2 Dk

2(k−ν) 2k+1

Dk = 2ν+1 + (2ν + 1)σ 2 V (K) 2k+1 n 2(k − ν)nβk2 Dk " 2ν+1 2(ν−k) # 2(k−ν) 2ν+1 2(k−ν) 2(k − ν) 2k+1 2(k − ν) 2k+1 = + n− 2k+1 σ 2 V (K) 2k+1 βk2 Dk 2k+1 = 2ν + 1 2ν + 1 2ν+1 2ν+1 2(k−ν) 2(k−ν) 2ν + 1 2(k − ν) 2k+1 2 − 2k+1 =n 1+ (σ V (K)) 2k+1 βk2 Dk 2k+1 . 2(k − ν) 2ν + 1

Důsledek 2.3.17. Pro hopt,ν,k platí

(vychýlení m(x; b hopt,ν,k ))2 2ν + 1 = . rozptyl m(x; b hopt,ν,k ) 2(k − ν)

Důkaz. Plyne z důkazu předchozí věty, neboť:

2(k − ν) 2ν + 1

2(ν−k) 2ν+1 2k+1 2(k − ν) 2k+1 2ν + 1 : = . 2ν + 1 2(k − ν)

31


2.4

Hledání optimálního jádra K

2.4.1 Definice. Řekneme, že reálná funkce f definovaná na konečném nebo nekonečném intervalu [a, b] má r znaménkových změn na [a, b], jestliže existuje r + 1 subintervalů [xi−1 , xi ] ⊂ [a, b], i = 1, . . . , r + 1, x0 = a, xr+1 = b takových, že platí: 1. f (x)f (y) ≥ 0 pro každé x, y ∈ [xi−1 , xi ], i = 1, . . . , r+1, přičemž f (x)f (y) > 0 nastává pro x, y ∈ Di ⊂ [xi−1 , xi ], kde Di má nenulovou Lebesgueovou míru, 2. f (x)f (y) ≤ 0 pro každé x ∈ [xi−1 , xi ] a y ∈ [xi , xi+1 ], i = 1, . . . , r. Označme ch(f ) počet znaménkových změn na [a, b]. Věta 2.4.2. Pro K ∈ Kνk , 0 ≤ ν < k platí ch(K) ≥ k − 2. Důkaz. Viz. [2]. 2.4.3 Definice. Legendreovy polynomy Pn jsou ortogonální polynomy na intervalu [−1, 1] s váhou w(x) = 1, přičemž platí 1.

Z1

Pn (x)Pm (x) dx =

−1

( 0

m 6= n , m=n

2 2n+1

2. Pn+1 (x) =

2n + 1 n xPn (x) − Pn−1 (x), n+1 n+1

kde P0 (x) = 1,

P1 (x) = x.

2.4.4 Definice. Nechť K ∈ Kν,k . Funkce, které jsou řešením minimalizačního problému Z1 min V (K), kde V (K) = K 2 (x) dx, K∈Kν,k

−1

se nazývají jádra s minimálním rozptylem. Věta 2.4.5. Jádra s minimálním rozptylem jsou jednoznačně určené polynomy stupně k − 2 omezené na intervalu [−1, 1]. Tyto polynomy jsou sudé funkce pro k sudé a liché funkce pro k liché. Mají právě k − 2 různých reálných kořenů na intervalu (−1, 1) a jejich explicitní formule je dána vztahem k−2

kde

(−1)ν ν! X e K(x) = (2r + 1)prν Pr (x), 2 r=ν Pr (x) =

r X i=0

je Legendreův polynom stupně r.

pri xi

32


Důkaz. Viz. [2]. 2.4.6 Definice. Nechť K ∈ Kν,k ∩ Nk−2 , kde Nk−2 je třída funkcí integrovatelných s druhou mocninou, které mají právě k −2 znaménkových změn na [−1, 1]. Funkce, které jsou řešením minimalizačního problému min

K∈Kν,k ∩Nk−2

T (K),

kde 2

T (K) = (V (K)k−ν |βk |2ν+1 ) 2k+1

2  2ν+1  2k+1 k−ν 1 Z Z1   xk K(x) dx =  K 2 (x) dx , 

−1

se nazývají optimální jádra.

−1

Věta 2.4.7. Optimální jádra jsou polynomy stupně k. Tyto polynomy jsou sudé funkce pro k sudé a liché funkce pro k liché. Mají právě k − 2 různých reálných kořenů na intervalu (−1, 1) a body −1, 1 jsou rovněž kořeny. Explicitní vyjádření je ve tvaru k (−1)ν ν! X (2r + 1)prν (Pr (x) − Pk (x)) , Kopt (x) = 2 r=ν kde Pr (x) je Legendreův polynom stupně r a Pk (x) Legendreův polynom stupně k. Důkaz. Viz. [2]. Poznámka. Optimální jádra jsou jádra, která minimalizují integrální střední kvadratickou chybu odhadu funkce m (viz. podkapitola 2.5). Definujme nyní ekvivalentní jádro Kδ k jádru K takto: x 1 Kδ (x) = ν+1 K , kde δ ∈ R, δ δ

δ > 0,

přičemž jádru K přísluší nosič [−1, 1], zatímco jádru Kδ přísluší nosič [−δ, δ]. Věta 2.4.8. Funkcionál T (K) je invariantní vzhledem k transformaci K(x) → Kδ (x) =

1 δ ν+1

x K( ). δ

33


Důkaz. 2  2ν+1  2k+1 k−ν δ Z Zδ   tk Kδ (t) dt T (Kδ ) =  Kδ2 (t) dt = 

−δ

 Zδ  = −δ

−δ

2 2ν+1  2k+1 k−ν δ Z 1 t  tk 1 K( t ) dt K 2 ( ) dt =  2ν+2 ν+1 δ δ δ δ

−δ

t = x 1 δ dt = dx δ

 Z1  = −1

2 2ν+1  2k+1 k−ν 1 Z k k x δ 1  2  = K (x)δ dx K(x)δ dx  2ν+2 ν+1 δ δ

−1

2  2ν+1  2k+1 k−ν 1 2 Z (k−ν)(2ν+1) 2k+1 Z1 δ   xk K(x) dx = K 2 (x) dx =  δ (2ν+1)(k−ν)

−1

−1

= T (K).

Pro důkaz věty 2.4.10 je užitečné následující lemma: Lemma 2.4.9. Nechť Pr je Legendreův polynom stupně r, 0 ≤ ν ≤ k − 2 a ν, k mají stejnou paritu. Pak 1. νpkν

=

k−1 X

(2r + 1)prν−1 ,

(2.4.1)

r=ν−1

2. Pr (x) =

1 ′ (P ′ (x) − Pr−1 (x)). 2r + 1 r+1

(2.4.2)

Důkaz. Viz. [2]. e Věta 2.4.10. Nechť K(x) ∈ Kν+1,k+1 je jádro s minimálním rozptylem a Kopt (x) ∈ Kν,k ∩ Nk−2 je optimální jádro. Pak d e Kopt (x) = K(x), dx

x ∈ (−1, 1).

34


Důkaz. Jádro s minimálním rozptylem je tvaru k−1 (−1)ν+1 (ν + 1)! X e K(x) = (2r + 1)prν+1 Pr (x), 2 r=ν+1

(2.4.3)

zatímco optimální jádro má tvar

k

Kopt (x) =

(−1)ν ν! X (2r + 1)prν (Pr (x) − Pk (x)) . 2 r=ν

Derivujeme-li předešlou rovnici podle x, pak k

(−1)ν ν! X d Kopt (x) = (2r + 1)prν (Pr′ (x) − Pk′ (x)) . dx 2 r=ν

(2.4.4)

Dále využijeme-li vztah (2.4.2), dostaneme Pr′ (x) − Pk′ (x) =

′ ′ ′ ′ = (Pr′ (x) − Pr+2 (x)) + (Pr+2 (x) − Pr+4 (x)) + · · · + (Pk−2 (x) − Pk′ (x)) = k−r

=−

2 X

(2(r + 2i − 1) + 1)Pr+2i−1 (x)

i=1

a dosazením do (2.4.4) máme d (−1) Kopt (x) = dx 2

ν+1

ν!

k X r=ν

 k−r  2 X (2r + 1)prν  (2(r + 2i − 1) + 1)Pr+2i−1 (x) . i=1

(2.4.5)

Porovnáme-li koeficienty u Pν+2j+1 (x) ve výrazech (2.4.3) a (2.4.5), tj.

a

(−1)ν+1 (ν + 1)! (2(ν + 2j + 1) + 1)pν+2j+1 ν+1 2 2j+ν

X (−1)ν+1 ν! (2(ν + 2j + 1) + 1) (2r + 1)prν , 2 r=ν

zjistíme, že jsou si rovny, neboť podle (2.4.1)

2j+ν

(ν +

1)pν+2j+1 ν+1

=

X

(2r + 1)prν .

r=ν

Na Obr. 2.4 a Obr. 2.6 je vyobrazeno jádro s minimálním rozptylem pro různé hodnoty k a ν, zatímco na Obr. 2.8 a Obr. 2.10 je vyobrazeno optimální jádro pro různé hodnoty k a ν, přičemž Epanečnikovo jádro je speciálním případ optimálního jádra pro k = 2 a ν = 0.

35


jadro s minimalnim rozptylem k = 2, ν = 0 1

0.8

0.6

0.4

0.2

y

0

−0.2

−0.4

−0.6

−0.8

−1 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(a) k = 2, ν = 0 jadro s minimalnim rozptylem k = 4, ν = 2

jadro s minimalnim rozptylem k = 4, ν = 0

10

2

8 1.5

6 1

4 0.5

2

0

y

y

0

−2 −0.5

−4 −1

−6 −1.5

−2 −2

−8

−1.5

−1

−0.5

0 x

0.5

1

1.5

−10 −2

2

(b) k = 4, ν = 0

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(c) k = 4, ν = 2



3


40

800

30

600

20

400

10

200

2

0

y

0

y

y

1

0

−10

−200

−20

−400

−1

−2 −30

−3 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

−40 −2

2

−600

−1.5

(e) k = 6, ν = 0

−1

−0.5

0 x

0.5

1

1.5

2

−800 −2

−1.5

(f) k = 6, ν = 2

−1

−0.5

0 x

0.5

1

(g) k = 6, ν = 4

Obr. 2.4: Jádro s minimálním rozptylem pro různé parametry k a ν sudé

k\ν

0

2

1 2 2

4 6

−3 (5x − 3) 8 15 (63x4 − 70x2 + 15) 128

2

4

15 (3x2 − 1) 4

− 105 (45x4 − 42x2 + 5) 32

945 (35x4 − 30x2 + 3) 16

Tab. 2.5: Jádro s minimálním rozptylem pro různé parametry k a ν sudé

1.5

2

36


jadro s minimalnim rozptylem k = 3, ν =1 2

1.5

1

0.5

y

0

−0.5

−1

−1.5

−2 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(a) k = 3, ν = 1 jadro s minimalnim rozptylem k = 5, ν =1


6

100

80

4 60

40

2

0

y

y

20

0

−20

−2 −40

−60

−4 −80

−6 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

−100 −2

2

−1.5

(b) k = 5, ν = 1

−1

−0.5

0 x

0.5

1

1.5

2

(c) k = 5, ν = 3

jadro s minimalnim rozptylem k = 7, ν =1


10


300

5000

8

4000

200

6

3000

4

2000

100

1000

0

y

0

y

y

2

−2

0

−1000

−100

−4

−2000

−6

−3000

−200

−8

−10 −2

−4000

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−300 −2

(e) k = 7, ν = 1

−1.5

−1

−0.5

0 x

0.5

1

1.5

(f) k = 7, ν = 3

2

−5000 −2

−1.5

−1

−0.5

0 x

0.5

1

(g) k = 7, ν = 5

Obr. 2.6: Jádro s minimálním rozptylem pro různé parametry k a ν liché

k\ν

1

3

−3 x 2

5

15 (14x3 − 10x) 16

7

− 105 (99x5 − 126x3 + 35x) 128

3

5

− 105 (5x3 − 3x) 4 945 (77x5 − 90x3 + 21x) 32

− 10395 (63x5 − 70x3 + 15x) 16

Tab. 2.7: Jádro s minimálním rozptylem pro různé parametry k a ν liché

1.5

2

37


optimalni jadro k = 2, ν = 0 2

1.5

1

0.5

y

0

−0.5

−1

−1.5

−2 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(a) k = 2, ν = 0 optimalni jadro k = 4, ν = 0

optimalni jadro k = 4, ν = 2

4

10

3

8

6

2 4

1

0

y

y

2

0

−2

−1

−4

−2 −6

−3 −8

−4 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−10 −2

(b) k = 4, ν = 0

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(c) k = 4, ν = 2



4

40

3

30

2

20

1

10


200

−1

0

y

0

y

y

400

0

−10 −200

−2

−20

−3

−30

−400

−4 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−40 −2

(e) k = 6, ν = 0

−1.5

−1

−0.5

0 x

0.5

1

1.5

(f) k = 6, ν = 2

2

−600 −2

−1.5

−1

−0.5

0 x

0.5

1

(g) k = 6, ν = 4

Obr. 2.8: Optimální jádro pro různé parametry k a ν sudé

k\ν

0

2

−3 (x2 − 1) 4

4

15 (7x4 − 10x2 + 3) 32

6

− 105 (33x6 − 63x4 + 35x2 − 5) 256

2

4

− 105 (5x4 − 6x2 + 1) 16 315 (77x6 − 135x4 + 63x2 − 5) 64

− 10395 (21x6 − 35x4 + 15x2 − 1) 32

Tab. 2.9: Optimální jádro pro různé parametry k a ν sudé

1.5

2

38



3

2

y

1

0

−1

−2

−3

−4 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(a) k = 3, ν = 1 optimalni jadro k = 5, ν = 1


10

60

8 40 6

4 20

0

0

y

y

2

−2 −20 −4

−6 −40 −8

−10 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−60 −2

−1.5

(b) k = 5, ν = 1

−1

−0.5

0 x

0.5

1

1.5

2

(c) k = 5, ν = 3



15


300

4000

3000

10

200

5

100

0

0

2000

y

y

y

1000

0

−1000

−5

−100

−10

−200

−2000

−3000

−15 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

−300 −2

(e) k = 7, ν = 1

−1.5

−1

−0.5

0 x

0.5

1

1.5

2

(f) k = 7, ν = 3

−4000 −2

−1.5

−1

−0.5

0 x

0.5

1

1.5

(g) k = 7, ν = 5

Obr. 2.10: Optimální jádro pro různé parametry k a ν liché

k\ν

1

3

15 (x3 − x) 4 − 105 (9x5 − 14x3 + 5x) 32 7 5 3

5 7

315 (143x 256

− 297x + 189x − 35x)

3

5

945 (7x5 − 10x3 + 3x) 16

− 10395 (39x7 − 77x5 + 45x3 − 7x) 64

135135 (33x7 − 63x5 + 35x3 − 5x) 32

Tab. 2.11: Optimální jádro pro různé parametry k a ν liché

2

39


2.5

Hledání optimálního řádu jádra k

Na kvalitu výsledného odhadu má vliv nejen šířka vyhlazovacího okna, ale i řád jádra. V následující podkapitole je popsán algoritmus pro nalezení optimálního řádu jádra. Nejprve ale dokážeme několik tvrzení s užitím Gasser-Müllerova odhadu. Věta 2.5.1. „Množstvíÿ vyhlazení s jádrem K a parametrem h je stejné jako „množstvíÿ vyhlazení s jádrem Kδ a parametrem h∗ = hδ , tj. (ν)

(ν)

m b GM (x; K, h) = m b GM (x; Kδ , h∗ ).

Důkaz. (ν)

m b GM (x; K, h) = =

=

n 1 X

hν+1

Yi

i=1

1 h∗ν+1

n X i=1

Zsi

si−1 Zsi

Yi

K

t−x h

1 δ ν+1

si−1

K

Zsi n X 1 t−x dt = ∗ν+1 ν+1 Yi dt = K h δ h∗ δ i=1

t−x δ h∗

dt =

1 h∗ν+1

n X i=1

Yi

si−1 Zsi

si−1

Kδ

t−x h∗

dt =

(ν)

=m b GM (x; Kδ , h∗ ). (ν)

Věta 2.5.2. Tvar integrální střední kvadratické chyby odhadu m b GM s jádrem Kδ a s parametrem h∗ je (ν)

MISE(m b GM (Kδ , h∗ )) =

2

σ nh∗2ν+1

Zδ

−δ

kde

Dk =

 δ 2 Z Kδ2 (x) dx + h∗2(k−ν) Dk  xk Kδ (x) dx ,

Z1

−δ

m(k) (x) k!

0

2

dx.

Důkaz. Plyne přímo z věty 2.3.14. Věta 2.5.3. Nechť K ∈ Kνk a nechť δ 2k+1 = oběma částem chyby je stejný, tj. Zδ

−δ

V (K) . βk2

Pak „příspěvekÿ jádra K k

 δ 2 Z Kδ2 (x) dx =  xk Kδ (x) dx . −δ

40


Důkaz. Zδ

Kδ2 (x) dx =

−δ

Zδ

1 δ 2ν+2

K2

−δ

=

Z1

1 δ 2ν+2

x δ

x =t dx = 1 δ = dx = dt δ

K 2 (t)δ dt = δ −(2ν+1)

−1

Z1

−1

|

K 2 (t) dt, {z

V (K)

}

 δ 2  δ 2 Z Z x =t x 1  xk Kδ (x) dx =  xk K dx = 1 δ = dx = dt δ ν+1 δ δ −δ −δ 2  1 2  1 Z k k Z t δ K(t)δ dx = δ 2(k−ν)  tk K(t) dt . = δ ν+1 −1 −1 | {z } βk2

V případě, že by platila rovnost δ −(2ν+1) V (K) = δ 2(k−ν) βk2 , pak by také platila rovnost δ 2k+1 = V β(K) 2 , což je původní předpoklad. k

Poznámka. Číslo δ =

V (K) βk2

1 2k+1

se nazývá kanonický faktor a značí se γ.

Věta 2.5.4. Nechť K ∈ Kνk , m ∈ C k [0, 1] a nechť h → 0, nh2ν+1 → ∞ pro n → ∞. Pak σ2 (ν) ∗ ∗2(k−ν) +h Dk , MISE(m b GM (Kγ , h )) = T (K) nh∗2ν+1 kde

Dk =

Z1

m(k) (x) k!

0

2

dx

a γ je kanonický faktor. Důkaz. Užitím věty 2.5.3 dostaneme (ν)

MISE(m b GM (Kγ , h∗ )) = =

2

σ nh∗2ν+1

Zγ

−γ

Zγ

−γ

|

 γ 2 Z Kγ2 (x) dx + h∗2(k−ν) Dk  xk Kγ (x) dx = −γ

Kγ2 (x) dx {z

V (Kγ )

}

σ2 ∗2(k−ν) +h Dk . nh∗2ν+1

41


Nyní upravíme V (Kγ ): V (Kγ ) =

Zγ

Kγ2 (x) dx =

−γ

=

−γ

1

γ

Zγ

V (K) = 2ν+1

x = t Z1 1 x K 2 ( ) dx = 1 γ K 2 (t) dt = = dx = dt γ 2ν+2 γ γ 2ν+1 1

γ

βk2 V (K)

2ν+1 2k+1

−1

2ν+1

V (K) = (βk2 ) 2k+1 V (K)

2(k−ν) 2k+1

= T (K).

Důsledek 2.5.5. Pro optimální hodnotu h∗opt,ν,k platí: 1 (2ν + 1)σ 2 2k+1 ∗ . hopt,ν,k = 2(k − ν)nDk

Důkaz. Spočítáme derivaci integrální střední kvadratické chyby odhadu podle h a položíme ji rovnu 0, tedy (ν) b GM (Kγ , h∗ )) (2ν + 1)σ 2 ∂MISE(m ∗2(k−ν)−1 = T (K) − + 2(k − ν)h Dk = 0. ∂h∗ nh∗2ν+2 Nyní vyjádříme-li z poslední rovnosti h, dostaneme 1 (2ν + 1)σ 2 2k+1 ∗ h = . 2(k − ν)nDk

1. Nechť ν, k jsou sudá čísla, 0 ≤ ν ≤ k − 2, pak

Důsledek 2.5.6.

(h∗opt,ν,k )2k+1 =

(2ν + 1)k ∗ (hopt,0,k )2k+1 . k−ν

(2.5.1)

2. Nechť ν, k jsou lichá čísla, 0 ≤ ν ≤ k − 2, pak (h∗opt,ν,k )2k+1 = Důkaz.

(2ν + 1)(k − 1) ∗ (hopt,1,k )2k+1 . 3(k − ν)

1. (h∗opt,0,k )2k+1 =

tedy h∗opt,ν,k h∗opt,0,k

!2k+1

=

(2.5.2)

σ2 , 2nkDk

(2ν + 1)σ 2 2nkDk (2ν + 1)k = , 2 2(k − ν)nDk σ k−ν

2. (h∗opt,1,k )2k+1

3σ 2 = , 2n(k − 1)Dk

tedy h∗opt,ν,k h∗opt,1,k

!2k+1

=

(2ν + 1)σ 2 2n(k − 1)Dk (2ν + 1)(k − 1) = . 2 2(k − ν)nDk 3σ 3(k − ν)

42


Důsledek 2.5.7. Platí (ν)

Důkaz.

MISE(m b GM (Kγ , h∗opt,ν,k )) = T (K) (2ν + 1)σ 2 2(k − ν)nDk

(h∗opt,ν,k )2k+1 =

2k + 1 σ2 . 2n(k − ν) (h∗opt,ν,k )2ν+1

Dk =

⇒

(2ν + 1)σ 2 2(k − ν)n(h∗opt,ν,k )2k+1

,

odtud σ

(ν)

MISE(m b GM (Kγ , h∗opt,ν,k )) = T (K)

= T (K)

σ

2

n(h∗opt,ν,k )2ν+1 2

T (K)σ = n(h∗opt,ν,k )2ν+1

n(h∗opt,ν,k )2ν+1

+ (h∗opt,ν,k )2(k−ν)

2ν + 1 1+ 2(k − ν)

2

+ (h∗opt,ν,k )2(k−ν) Dk

(2ν + 1)σ

2

2(k − ν)n(h∗opt,ν,k )2k+1

= T (K)

2k + 1 σ2 . 2n(k − ν) (h∗opt,ν,k )2ν+1

Definujme veličinu

a označme

2k + 1 1 b L(k) = T (Kopt,ν,k ) 2ν+1 ∗ 2n(k − ν) h opt,ν,k Iν (k0 ) =

Pro ν = 0 dostaneme

k0 − ν ν + 2j; j = 1, . . . , 2

L(k) = T (Kopt,ν,k )

2k + 1 1 2nk h∗opt,0,k

a I(k0 ) =

!

k0 2j; j = 1, . . . , 2

.

.

=

!

=

43


Algoritmus pro nalezení optimální hodnoty k: 1. pro sudé ν: Nechť k ∈ Iν (k0 ) a ν, k jsou sudá. (a) Pro každé k ∈ I(k0 ) najdeme optimální jádro Kopt,ν,k . (b) Pro každé k ∈ I(k0 ) a každé optimální jádro Kopt,ν,k najdeme metodou křížového ověřování (3.2.2) optimální hodnotu b hCV,0,k . b hCV,0,k Dále b h∗opt,0,k = a ze vztahu (2.5.1) vypočteme b h∗opt,ν,k . γ0,k (c) Jako kritérium pro volbu optimálního řádu k jádra budeme minimalib zovat funkci L(k) vzhledem ke k, tj. b b k = argmin L(k). k∈Iν (k0 )

2. pro liché ν: Nechť k ∈ Iν (k0 ) a ν, k jsou lichá.

(a) Pro každé k ∈ I(k0 ) najdeme optimální jádro Kopt,ν,k . (b) Pro každé k ∈ I(k0 ) a každé optimální jádro Kopt,ν,k najdeme metodou křížového ověřování (2.2.4) optimální hodnotu b hCV (1) ,1,k . b hCV (1) ,1,k Dále b h∗opt,1,k = a ze vztahu (2.5.2) vypočteme b h∗opt,ν,k . γ1,k (c) Jako kritérium pro volbu optimálního řádu k jádra budeme minimalib zovat funkci L(k) vzhledem ke k, tj. b b k = argmin L(k). k∈Iν (k0 )

Odhad m b funkce m pak sestrojíme pomocí b k, Kopt,ν,bk a b hopt,ν,bk .

Vyhlazení optimálním jádrem pro různý jádrový odhad a různý řád k při stejné hodnotě h můžeme porovnat z Obr. 2.12.

44


vyhlazeni optimalnim jadrem k=4, ν=0, h=0.0509

vyhlazeni optimalnim jadrem k=6, ν=0, h=0.0509 30

20

20

20

10

10

10

0

0

0

y

30

y

y


−10

−10

−10

−20

−20

−20

−30

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

0.1

0.2

(a) k = 2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

−30

1

0

0.1

0.2

(b) k = 4

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

(c) k = 6

1. Nadaraya-Watsonovy odhady vyhlazeni optimalnim jadrem k=4, ν=0, h=0.0509


20

20

20

10

10

10

0

0

0

y

30

y

y


−10

−10

−10

−20

−20

−20

−30

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

0

0.1

0.2

(d) k = 2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

−30

1

0

0.1

0.2

(e) k = 4

0.3

0.4

0.5 x

0.6

0.7

(f) k = 6

2. Lokální-lineární odhady vyhlazeni optimalnim jadrem k=4, ν=0, h=0.0509


20

20

20

10

10

10

0

0

0

y

30

y

y


−10

−10

−10

−20

−20

−20

−30

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

(g) k = 2

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(h) k = 4

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

(i) k = 6

3. Gasser-Müllerovy odhady Obr. 2.12: Vyhlazení optimálním jádrem pro různý jádrový odhad, různý řád k a stejnou šířku okna h

Kapitola 3 Vyhlazovací splajny 3.1

Úvodní pojmy

V této kapitole budeme předpokládat, že odhadovaná funkce m je dostatečně hladká, přesněji, že funkce m je prvkem Sobolevova prostoru. 3.1.1 Definice. Nechť L2 [a, b] je množina všech funkcí f , jejichž druhá mocnina je integrovatelná na intervalu [a, b]. Pak Sobolevovým prostorem nazveme množinu W2k [a, b] = {f ; f, f (1) , . . . , f (k−1) absolutně spojité, f (k) ∈ L2 [a, b]}. 3.1.2 Definice. Nechť k ∈ N a λ > 0. Funkci, která minimalizuje funkcionál Z n X 2 (yi − f (xi )) +λ (f (k) (x))2 dx b

−1

Φn,λ,k (f ) = n |

i=1

{z

}

ΦSU M (f )

|a

{z

ΦINTk (f )

}

přes všechna f ∈ W2k [a, b], nazýváme vyhlazovacím splajnem hodnot y1 , . . . , yn . Parametr λ se nazývá parametr vyhlazení. Funkcionál Φn,λ,k (f ) je součtem dvou členů, přičemž první z nich udává kvadrát vzdálenosti mezi vektorem dat y1 , . . . , yn a vektorem hodnot funkce f (x) v bodech x1 , . . . , xn a druhý vyjadřuje určitou míru porušení podmínky hladkosti. Pomocí parametru λ lze oběma členům přisoudit různou váhu. Při větších hodnotách λ bude preferována hladkost, při menších hodnotách bude výsledný splajn více kopírovat data. 3.1.3 Definice. Nechť θ0 , . . . , θr−1 , δ1 , . . . , δl ∈ R. Splajnem řádu r s uzly v ξ1 , . . . , ξl nazveme funkci tvaru s(x) =

r−1 X i=0

i

θi x +

l X i=1

45

δi (x − ξi )r−1 + ,

(3.1.1)

46

KAPITOLA 3. VYHLAZOVACÍ SPLAJNY

kde

( x − ξi (x − ξi)+ = 0

pro x − ξi ≥ 0 pro x − ξi < 0

se nazývá useknutá funkce. Věta 3.1.4. Funkce s je splajn řádu r s uzly v ξ1 , . . . , ξl , tj. s lze vyjádřit ve tvaru (3.1.1), právě tehdy, když platí následující podmínky: 1. funkce s je po částech polynomiální stupně r − 1 na podintervalu [ξi , ξi+1), 2. funkce s má r − 2 spojitých derivací na celém definičním oboru, 3. funkce s má derivaci r − 1, která je skoková funkce se skoky v ξ1 , . . . , ξl . Důkaz. Dokažme nejprve, že za platnosti všech třech podmínek, má funkce s(x) tvar (3.1.1). Platí-li podmínka 1., pak existují polynomy p0 (x), p1 (x), . . . , pl (x) definované na celém definičním oboru stupně r − 1 takové, že s(x) lze vyjádřit na podintervalech (−∞, ξ1 ), [ξ1 , ξ2), . . . , [ξl , ∞) takto:   p0 (x) pro x ∈ (−∞, ξ1 ),    p0 (x) + p1 (x) pro x ∈ [ξ1 , ξ2 ), s(x) = . ..     p (x) + p (x) + · · · + p (x) pro x ∈ [ξ , ∞). 0

1

l

l

Z podmínky 2. plyne, že

lim− s(i) (x) = lim+ s(i) (x)

x→ξ1

pro všechna

i = 0, . . . , r − 2,

x→ξ1

tedy (i)

(i)

(i)

lim− p0 (x) = lim+ p0 (x) + lim+ p1 (x)

x→ξ1

x→ξ1

pro všechna i = 0, . . . , r − 2.

x→ξ1

(i)

To znamená, že limx→ξ1+ p1 (x) = 0 pro všechna i = 0, . . . , r − 2, neboť funkce (i)

(i)

(i)

p0 (x) je spojitá v ξ1 . Analogicky limx→ξ2+ p2 (x) = 0, . . . , limx→ξ+ pl (x) = 0 pro l všechna i = 0, . . . , r − 2. Jelikož polynom p1 (x) má r − 2 derivací zprava rovných nule v bodě ξ1 , pak p1 (x) má r − 1-násobný kořen v ξ1 . Analogicky p2 (x) má r − 1násobný kořen v ξ2 , p3 (x) má r − 1-násobný kořen v ξ3 , atd. Dále, nechť platí podmínka 3., pak lim− s(r−1) (x) 6= lim+ s(r−1) (x), x→ξ1

tedy

(r−1)

lim p0

x→ξ1−

x→ξ1

(r−1)

(x) 6= lim+ p0 x→ξ1

(r−1)

(x) + lim+ p1 x→ξ1

(x).

47


(r−1)

Pak ale limx→ξ1+ p1

(r−1)

(x) 6= 0, neboť funkce p0

(r−1)

(x) je spojitá v ξ1 . Analogicky

(r−1)

(r−1)

limx→ξ2+ p2 (x) 6= 0, . . . , limx→ξ+ pl (x) 6= 0. Nyní, jelikož limx→ξ1+ p1 (x) 6= l 0 a p1 (x) má r − 1-násobný kořen v ξ1 , pak funkci p1 (x) můžeme napsat ve tvaru p1 (x) = δ1 (x − ξ1 )r−1 , kde δ1 je nějaká konstanta. Analogicky můžeme napsat polynomy p2 (x), . . . , pl (x). Definujme nyní funkce f1 , . . . , fl takto: ( 0 pro x ∈ (−∞, ξ1 ) f1 (x) = , p1 (x) pro x ∈ [ξ1 , ∞) ( 0 pro x ∈ (−∞, ξ2 ) f2 (x) = , p2 (x) pro x ∈ [ξ2 , ∞) .. . ( 0 pro x ∈ (−∞, ξl ) fl (x) = . pl (x) pro x ∈ [ξl , ∞) Pak platí, že s(x) = p0 (x) + f1 (x) + · · · + fl (x)

pro všechna x ∈ R.

Na Obr. 3.1 je ilustrační příklad, kde se k polynumu p0 (x) přičítají dvě useknuté funkce. 300 p (x) 0 f1(x) f (x) 2

s(x)

200

y

100

0

−100

−200

−300

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Obr. 3.1: Znázornění funkce s(x) jako součet funkcí p0 (x), f1 (x) a f2 (x)

48


Funkce fi , i = 1, . . . , l lze přepsat do tvaru fi (x) = δi (x − ξi )r−1 + , i = 1, . . . , l, a jelikož je p0 (x) polynom stupně r − 1, můžeme ho vyjadřit takto: p0 (x) =

r−1 X

θi xi .

i=0

Funkci s(x) tedy můžeme vyjádřit ve tvaru r−1 X

s(x) =

i

θi x +

i=0

l X

δi (x − ξi )r−1 + .

i=1

Nyní dokažme opačný směr ekvivalence, tj. je-li funkce s(x) tvaru (3.1.1), platí podmínky 1. až 3.. Jelikož fi (x) = δi (x − ξi)r−1 + , i = 1, . . . , l, zřejmě (r−2)

lim− fi

x→ξi

x→ξi

(r−2)

lim+ fi

x→ξi (r−2)

tedy fi

(x) = lim− 0 = 0, (x) = lim+ δi (r − 1)!(x − ξi ) = 0, x→ξi

(x) je spojitá v bodě ξi , i = 1, . . . , l a musí nutně platit i 2.. Dále (r−1)

lim− fi

x→ξi

x→ξi

(r−1)

lim fi

x→ξi+ (r−1)

a zároveň je fi podmínky 3..

(x) = lim− 0 = 0, (x) = lim+ δi (r − 1)! = δi (r − 1)!, x→ξi

(x) = δi (r − 1)! pro všechna x ∈ [ξi , ∞), což dokazuje platnost

Na Obr. 3.2 je ilustrační příklad, kde se k polynumu p0 (x) přičítají dvě useknuté funkce. 5

10

x 10

9

p(3)(x) 0 f(3)(x) 1 (3) f2 (x)

8

s(3)(x)

7 6

y

5 4 3 2 1 0 −1

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

Obr. 3.2: Znázornění derivace r − 1 funkce s(x) jako součet derivací r − 1 funkcí p0 (x), f1 (x) a f2 (x)

49


Podmínka 1. plyne přímo z toho, že pro každé x ∈ [ξj , ξj+1) platí r−1 X

s(x) =

i

θi x +

i=0

j X

δi (x − ξi )r−1 ,

i=1

tedy, že s(x) je polynom stupně r − 1 na tomto intervalu. Jinými slovy, splajn je po částech polynomiální funkce, jejíž jednotlivé polynomiální segmenty jsou na sebe navázany v uzlových bodech ξ1 , . . . , ξl takovým způsobem, že zaručují určitá kritéria hladkosti. Pokud bychom požadovali, aby i derivace řádu r − 1 byla spojitá, pak by tento splajn byl jediným polynomem stupně r − 1. Pro lepší pochopení funkce δi (x − ξi)r−1 je vhodný Obr. 3.3, který znázorňuje + konstrukci splajnu 4. řádu s uzlem v ξi, který je po částech polynomiální stupně 3.

3

−30(x−0.5)+ + 20(x−0.25)(x−0.5)(x−0.75) 2

1

1

1

0

0

0

−1

−1

−1

y

2

y

y

3

−30 * (x−0.5)+

20(x−0.25)(x−0.5)(x−0.75) 2

−2

−2

−2

−3

−3

−3

−4

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

3 (a) y = 20(x3 − 32 x2 + 11 16 x− 32 )

−4

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

(b) y = −30(x − 12 )3+

0.9

1

−4

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(c) součet obou funkcí

Obr. 3.3: Význam funkce δi (x − ξi )3+ při konstrukci splajnu 4. řádu s uzlem v ξi Označme symbolem S r (ξ1 , . . . , ξl ) množinu všech funkcí tvaru (3.1.1). Poznámka. Není obtížné dokázat, že S r (ξ1 , . . . , ξl ) je vektorový prostor. Jelikož r−1 funkce 1, x, . . . , xr−1 , (x − ξ1 )r−1 jsou lineárně nezávislé, má tento + , . . . , (x − ξl )+ prostor dimenzi r + l. Nadále se budeme zabývat přirozenými splajny, přičemž za uzly ξ1 , . . . , ξl budeme brát body plánu x1 , . . . , xn z regresního modelu (1.0.1), tedy l = n. 3.1.5 Definice. Nechť platí podmínky z věty 3.1.4 pro r = 2k s uzly v x1 , . . . , xn a navíc i čtvrtá podmínka, že funkce s je polynonom stupně k − 1 mimo interval [x1 , xn ], pak splajn (3.1.1) nazveme přirozeným splajnem řádu 2k s uzly v x1 , . . . , xn . Poznámka. Název přirozeného splajnu plyne z toho, že platí-li čtvrtá podmínka, jsou také splněny přirozené hraniční podmínky s(i) (a) = 0 = s(i) (b),

i = k, . . . , 2k − 1.

50


Označme symbolem N S 2k (x1 , . . . , xn ) množinu všech přirozených splajnů řádu 2k s uzly v x1 , . . . , xn . Poznámka. Prostor N S 2k (x1 , . . . , xn ) má dimenzi n a je podprostorem vektorového prostoru S 2k (x1 , . . . , xn ). Navíc ze čtvrté podmínky plyne, že je-li s přirozený splajn, musí platit θk = · · · = θ2k−1 = 0 (3.1.2) v (3.1.1), jelikož s musí být polynom stupně k − 1 pro x < x1 . Na Obr. 3.4 vidíme přirozený splajn pro k = 2, který nazýváme kubický splajn, neboť je po částech polynomiální stupně 3 na intervalu [ξ1 , ξl ]. Jelikož se jedná o přirozený splajn, tak na intervalech [a, ξ1 ] a [ξl , b] je polynomiální stupně 1, tedy lineární. Zejména si všimněme, že platí přirozené hraniční podmínky. prirozeny splajn, k = 2

prirozeny splajn − 1. derivace, k = 2

20

80

15

60

40

10

20

5

y

y

0

0

−20

−5 −40

−10 −60

−15

−80

−20 −0.2

0

0.2

0.4

0.6

0.8

1

−100 −0.2

1.2

0

0.2

0.4

0.6

x

0.8

1

1.2

1

1.2

x

(a) kubický splajn

(b) 1. derivace kubického splajnu 4

5



x 10

1500

4 1000

3

2

500

1

y

y

0

−500

0

−1

−2 −1000

−3 −1500

−4

−2000 −0.2

0

0.2

0.4

0.6

0.8

x

(c) 2. derivace kubického splajnu

1

1.2

−5 −0.2

0

0.2

0.4

0.6

0.8

x

(d) 3. derivace kubického splajnu

Obr. 3.4: Přirozený splajn pro k = 2 (kubický splajn)

51


Dokažme nyní následující pomocné lemma. Lemma 3.1.6 (Lyche a Schumaker). Nechť funkce g1 , . . . , gn tvoří bázi prostoru N S 2k (x1 , . . . , xn ). Pak existují koeficienty θ0,j , . . . , θk−1,j , δ1,j , . . . , δn,j takové, že platí k−1 n X X i gj (x) = θi,j x + δi,j (x − xi )2k−1 . (3.1.3) + Jestliže s(x) = Zb

f

i=0

Pn

j=1 βj gj (x)

(k)

i=1

af∈

W2k [a, b],

(k)

pak platí

k

(x)s (x) dx = (−1) (2k − 1)!

n X

f (xi )

i=1

a

n X

(3.1.4)

βj δi,j .

j=1

Důkaz. Rovnost (3.1.3) je pouhé přeformulování toho, že N S 2k (x1 , . . . , xn ) ⊂ S 2k (x1 , . . . , xn ) a že platí (3.1.2). Zbývá tedy dokázat, že platí (3.1.4). Nejprve si všimněme, že derivace s(k+j), j = 0, . . . , k − 1 jsou nulové mimo interval [x1 , xn ], což plyne ze čtvrté podmínky. Dále položme x0 = a, xn+1 = b a upravme následující integrál: Zb

f

(k)

(k)

k−1

(x)s (x) dx = (−1)

Zb

′

f (x)s(2k−1) (x) dx =

a

a

= (−1)k−1

n X

x n−1 Zi+1 X ′ f (x)s(2k−1) (x) dx = (−1)k−1 f (x)s(2k−1) (x) dx =

x Zi+1

i=0 x i

′

i=1 x i

" !# n−1 n i X X X = (−1)k−1 (f (xi+1 ) − f (xi )) βj (2k − 1)! δr,j = i=1

= (−1)k−1 (2k − 1)!

n X j=1

k

= (−1) (2k − 1)!

n X

βj

j=1

r=1

j=1

"

βj

n X

f (xi )

i=2

n X

i−1 X

δr,j −

r=1

n−1 X

f (xi )

i=1

i X r=1

δr,j

!#

=

f (xi )δi,j .

i=1

V posledním kroku jsme využili toho, že pro x > xn platí (2k−1)

gj

(x) = (2k − 1)!

n X

δi,j = 0,

j=1

což plyne ze čtvrté podmínky.

Nyní se zabývejme konstrukcí odhadu m(x; b λ) funkce m, který minimalizuje funkcionál Φn,λ,k (f ) definovaný v 3.1.2. Přitom budeme předpokládat, že parametr λ je pevně zvolen a platí 0 ≤ x1 < · · · < xn ≤ 1, tedy a = 0 a b = 1.

52


Nechť ω je nějaká konstanta a f libovolná funkce taková, že f ∈ W2k [0, 1]. Zaveďme následující označení: 1 Ψ (m(x; b λ), f (x); ω) = [ΦSU M (m(x; b λ) + ωf (x)) + λΦIN Tk (m(x; b λ) + ωf (x))] = 2   Z1 n X 1 = n−1 (yi − m(x b i ; λ) − ωf (xi ))2 + λ (m b (k) (x; λ) + ωf (k) (x))2 dx , 2 i=1 0

(3.1.5)

kde ΦSU M a ΦIN Tk jsou definovány v 3.1.2. Jelikož m b je minimalizátor, musí platit Ψ (m(x; b λ), f (x); 0) ≤ Ψ (m(x; b λ), f (x); ω)

pro všechna ω 6= 0.

Je zřejmé, že (3.1.5) je diferencovatelná funkce podle ω a platí ∂Ψ (m(x; b λ), f (x); ω) = ∂ω Z1 n X −1 = −n f (xi )(yi − m(x b i ; λ) − ωf (xi)) + λ f (k) (x)(m b (k) (x; λ) + ωf (k) (x)) dx. ′

Ψ (m(x; b λ), f (x); ω) = i=1

0

(3.1.6)

Nutnou podmínkou, aby funkce m(x; b λ) minimalizovala funkcionál Φn,λ,k (f ) definovaný v 3.1.2, je ∂Ψ (m(x; b λ), f (x); ω) ′ = 0. (3.1.7) Ψ (m(x; b λ), f (x); 0) = ∂ω ω=0 Nyní můžeme dokázat následující větu.

Věta 3.1.7. Nechť g1 , . . . , gn je báze prostoru N S 2k (x1 , . . . , xn ) a předpokládejme, že n ≥ k. Pro pevné 0 < λ < ∞ existuje jediná funkce m(x; b λ) minimalizující funk2k k cionál Φn,λ,k b λ) ∈ N S (x1 , . . . , xn ), a proto lze napsat P(f ) na W2 [0, 1]. Navíc m(x; m(x; b λ) = nj=1 βjλ gj . Koeficienty β¯λ = (β1λ , . . . , βnλ )T jsou řešením rovnice (G + nλH)β¯λ = y¯,

kde G = {gj (xi )}ni,j=1 , a

n H = (−1)k (2k − 1)!δi,j i,j=1 ,

kde δi,j jsou koeficienty pro bázové funkce g1 , . . . , gn definované v (3.1.3).

(3.1.8)

53


Důkaz. Nejprve si všimneme, že pro libovolné funkce f1 , f2 ∈ W2k [0, 1], obdobně jako v (3.1.6), platí ′

Ψ (f1 , f2 ; 0) = −n−1

n X

f2 (xi )(yi − f1 (xi )) + λ

i=1

Z1

(k)

(k)

f1 (x)f2 (x) dx.

0

Tedy, nutná podmínka, aby funkce f1 minimalizovala funkcionál Φn,λ,k (f ) definovaný v 3.1.2, je n

1 X f2 (xi )(yi − f1 (xi )) = nλ i=1

Z1

(k)

(k)

f1 (x)f2 (x) dx

pro všechny

f2 ∈ W2k [0, 1].

0

P Kdyby f1 byl přirozený splajn, pak bychom ho mohli napsat ve tvaru f1 = nj=1 βj gj pro nějaké koeficienty β1 , . . . , βn . Užitím lemmatu 3.1.6 můžeme nutnou podmínku přepsat ve tvaru n n n n X X X 1 X k f2 (xi )(yi − βj gj (xi )) = (−1) (2k − 1)! f2 (xi ) βj δi,j . nλ i=1 j=1 i=1 j=1

Jelikož toto musí platit pro všechny funkce f2 , je tato podmínka ekvivalentní tomu, že βj splňuje n X gj (xi ) + nλ(−1)k (2k − 1)!δi,j βj = yi , j=1

což je (3.1.8). Jestliže dokážeme, že existuje jediný vektor β¯λ , který je řešením (3.1.8), pak budeme moci ukázat, že existuje přirozený splajn, který splňuje nutnou podmínku (3.1.7) pro minimalizaci funkcionálu Φn,λ,k (f ). Abychom ověřili, že vektor β¯λ je jednoznačně definován rovnicí (3.1.8), využijeme poznatku, že systém A¯ x = ¯b má jediné řešení právě tehdy, když x¯ = 0 je jediné řešení systému A¯ x = 0. Jelik kož m(x; b λ) ∈ W2 [0, 1], proto také platí nutná podmínka (3.1.7) pro minimalizaci funkcionálu Φn,λ,k (f ), v tomto případě ′

Ψ (m(x; b λ), m(x; b λ); 0) = n−1

n X i=1

2

(m(x b i ; λ)) + λ

Z1 0

m(x; b λ)(k) (x)

2

dx = 0.

Proto (m(x b i ; λ))2 = 0, i = 1, . . . , n a m b (k) (x; λ) = 0 skoro všude. To znamená, že m(x; b λ) je polynom stupně k − 1, který je roven nule v n ≥ k bodech a tedy m(x; b λ) ≡ 0. Nicméně, jelikož g1 , . . . , gn je báze prostoruPN S 2k (x1 , . . . , xn ), jsou g1 , . . . , gn lineárně nezávislé, což implikuje, že m(x; b λ) = nj=1 βjλ gj (x) = 0 právě tehdy, když β1λ = · · · = βnλ = 0. Z toho plyne, že (3.1.8) má jediné řešení.

54


P Nyní zbývá dokázat, že m(x; b λ) = nj=1 βjλ gj minimalizuje Φn,λ,k (f ). Užitím toho, ′ že Ψ (m(x; b λ), f (x); 0) = 0 pro f ∈ W2k [0, 1] dostaneme n−1

n X

2

(yi − f (xi )) + λ

i=1

−1

=n

n X i=1

−1

+n

n X

n X i=1

f (k) (x)

0

2

(yi − m(x b i ; λ)) + λ

i=1

≥ n−1

Z1

Z1 0

(m(x b i ; λ) − f (xi )) + λ Z1 0

dx =

m b (k) (x; λ)

2

(yi − m(x b i ; λ))2 + λ

2

Z1

2

′

dx + 2Ψ (m(x; b λ), m(x; b λ) − f (x); 0)+

m b (k) (x; λ) − f (k) (x)

0

m b (k) (x; λ)

2

dx,

2

dx ≥

(3.1.9)

′

neboť m(x; b λ) − f (x) ∈ W2k [0, 1] a Ψ (m(x; b λ), m(x; b λ) − f (x); 0) = 0. Všechny ostatní sčítance musí být nezáporné, proto pro žadnou jinou funkci P z W2k [0, 1] n nemůže být hodnota funkcionálu Φn,λ,k (f ) menší. Tedy m(x; b λ) = j=1 βjλ gj opravdu minimalizuje Φn,λ,k (f ). Abychom se přesvědčili, že tento minimalizátor je jediný, předpokládejme, že funkce m(x; b λ) a f (x) dávají stejnou hodnotu funkcionálu. Pak z (3.1.9) plyne, že Z1 0

a

m b (k) (x; λ) − f (k) (x)

n X i=1

2

dx = 0

(m(x b i ; λ) − f (xi ))2 = 0.

První z těchto vztahů má za důsledek, že k-tá derivace m(x; b λ) − f (x) je nulová skoro všude. To znamená, že m(x; b λ) − f (x) musí být polynom stupně k − 1. Druhá rovnost zaručuje, že m(x b i ; λ) − f (xi ) = 0, i = 1, . . . , n a jelikož n ≥ k platí m(x; b λ) ≡ f (x). Předchozí věta je pro naše účely poměrně důležitá, neboť zaručuje existenci a jednoznačnost odhadu m(x; b λ) vyhlazovacím splajnem a také nám dává návod na jeho konstrukci. Konkrétně, odhad m b vyhlazovacím splajnem je přirozený splajn řádu 2k, tj. n X m(x; b λ) = βjλ gj (x) j=1

s βjλ získaných z (3.1.8).


55

Podívejme se blíže na vyhlazovací parametr λ ve větě 3.1.7. Mohou nastat dva extrémní případy, kdy λ = 0 a λ = ∞. Je-li λ = 0, pak hledaný přirozený splajn splňuje podmínku n X yi = βj0 gj (xi ), i = 1, . . . , n. j=1

Takovému splajnu, říkáme přirozený interpolační splajn pro y1 , . . . , yn . Z důkazu věty 3.1.7 lze odvodit následující důsledek.

Důsledek 3.1.8. Jestliže n ≥ k, přirozený interpolační splajn pro y1 , . . . , yn je funkce f , která jednoznačně minimalizuje ΦIN Tk (f ) na množině W2k [0, 1] za podmínky f (xi ) = yi , i = 1, . . . , n. Situace, kdy λ = ∞ převede problém hledání odhadu m(x; b λ) na běžnou polynomiální regresi řádu k, neboť v tomto případě platí ΦIN Tk (f ) = 0 a minimalizujeme tedy pouze ΦSU M (f ), kde ΦSU M a ΦIN Tk jsou definovány v 3.1.2.

56


3.2

Hledání optimálního parametru vyhlazení λ a parametru k

Vraťme se opět k větě 3.1.7. Rovnice (3.1.8) lze přepsat do vhodnějšího tvaru. Vynásobíme-li ji zleva maticí GT dostaneme (GT G + nλGT H)β¯λ = GT y¯, a tedy GT H =

(

(−1)k (2k − 1)!

n X

gr (xi )δi,j

r=1

.

i,j=1

To lze ovšem podle lemmatu 3.1.6 zapsat jako  1 n Z  (k) (k) gi (x)gj (x) dx   0

)n

.

i,j=1

Důsledek 3.2.1. Pro ) na W2k [0, 1] Pn n ≥ k je minimalizátor¯ funkcionálu Φn,λ,k (f T funkce m(x; b λ) = jsou řešením j=1 βjλ gj , kde koeficinty βλ = (β1λ , . . . , βnλ ) rovnice (GT G + nλΩ)β¯λ = GT y¯, (3.2.1) kde Ω=

 1 Z 

0

n  (k) (k) gi (x)gj (x) dx 

.

i,j=1

Tvar rovnice (3.2.1) můžeme dále upravit následovně: (GT G + nλΩ)β¯λ = GT y¯ β¯λ = (GT G + nλΩ)−1 GT y¯ Gβ¯ = G(GT G + nλΩ)−1 GT y¯ |{z}λ m(x;λ) b

Odhad m(¯ b x; λ) = (m(x b 1 ; λ), . . . , m(x b n ; λ))T můžeme tedy napsat ve tvaru kde

m(¯ b x; λ) = Sλ y¯,

Sλ = G(GT G + nλΩ)−1 GT . Matici Sλ se nazývá, stejně jako u jádrových odhadech, vyhlazovací (klobouková) matice (viz. Obr. 3.5). Klobouková matice Sλ = {sijλ }ni,j=1 je symetrická a pozitivně definitní.

57


vyhlazovaci matice splajnu, λ = 0.0001

vyhlazovaci matice splajnu, λ = 0.01

0.25

0.5

0.2

0.4

0.15 hodnoty matice

hodnoty matice

0.6

0.3 0.2 0.1 0

0

−0.1 0

0.1 0.05 0 −0.05

5 10

5 10

0

−0.1 0

5 10

5 10

15

15

15

15

20

20

radky matice S

20

λ

sloupce matice Sλ

20

radky matice S

λ

sloupce matice S

λ

(a) λ = 0, 0001

(b) λ = 0, 01

Obr. 3.5: Vyhlazovací matice Pozitivní definitnost matice Sλ = G(GT G + nλΩ)−1 GT není na první pohled zřejmá, tudíž si tuto vlastnost nyní ověříme. O pozitivně definitních maticích platí následující vlastnosti: 1. nechť R je regulární matice a S je pozitivně definitní matice, pak RT SR je také pozitivně definitní matice, 2. součet pozitivně definitních matic je také pozitivně definitní matice, 3. je-li S pozitivně definitní matice, pak inverzní matice S −1 je také pozitivně definitní matice. Stačí tedy dokázat, že Ω je pozitivně definitní. Definujeme-li na prostoru W2k [0, 1] skalární součin funkcí f1 a f2 následovně: hf1 , f2 i =

Z1

f1 (x)f2 (x) dx,

0

pak Z1

(k)

(k)

gi (x)gj (x) dx

0

(k) gi (x)

(k)

je skalární součin funkcí a gj (x), kde gi , i = 1, . . . , n jsou bázové funkce k prostoru W2 [0, 1]. Jelikož funkce gi , i = 1, . . . , n jsou lineárně nezávislé, jsou i (k) funkce gi , i = 1, . . . , n lineárně nezávislé. Dále  1 n Z  (k) (k) Ω= gi (x)gj (x) dx ,   0

i,j=1

58


(k)

tedy Ω je Grammova matice pro systém funkcí gi , i = 1, . . . , n, která je vždy pozitivně definitní, a proto Ω je pozitivně definitní matice. K volbě optimálního parametru vyhlazení, můžeme opět využít metody křížového ověřování, kterou jsme definovali již v podkapitole 2.2. Chceme-li odhadnout vyhlazovací parametr λ, hledáme argument minimuma funkce křížového ověřování CV (λ) na množině všech nezáporných reálných čísel R+ 0. bCV = argmin CV (λ), λ

(3.2.2)

λ∈R+ 0

kde

přičemž

n

1X (Yi − m b −i (xi ; λ))2 , CV (λ) = n i=1 m b −i (xi ; λ) =

n X j=1 j6=i

sijλ yj , 1 − siiλ

je odhad funkce m(x) v bodě xi bez použití bodu xi . Věta 3.2.2. Nechť m(x; b λ) je odhad m(x) a nechť platí, že m b −i (xi ; λ) =

n X j=1 j6=i

sijλ Yj + siiλ m b −i (xi ; λ).

Předpokládejme, že siiλ 6= 1 pro i = 1, . . . , n. Pak funkci CV (λ) můžeme vyjádřit ve tvaru 2 n b i ; λ) 1 X Yi − m(x . (3.2.3) CV (λ) = n i=1 1 − siiλ Důkaz. Viz. důkaz věty 2.2.2.

Poznámka. Vztah (3.2.3) můžeme opět zobecnit na tvar n

1X GCV (λ) = n i=1

Yi − m(x b i ; λ) 1 − tr(Sλ )/n

2

kde tr(Sλ ) je stopa matice Sλ = (sijλ ), tj. platí tr(Sλ ) = zobecněnou metodu křížového ověřování.

, Pn

i=1

siiλ , a získat tak

59


Nás ovšem nezajímá pouze parametr λ, ale také hodnota parametru k, která ovlivňuje řád splajnu. Tuto dvojici parametrů odhadneme jako argument minima funkce GCV (λ, k) na množině R+ 0 × {1, . . . , n}, tedy bGCV , b (λ kGCV ) =

Zavedeme-li množinu

argmin

GCV (λ, k).

λ∈R+ 0 , k∈{1,...,n}

bGCV ; k = 1, . . . , n}, Λb = {λ k

kde

bGCV = argmin GCV (λ, k), λ k λ∈R+ 0

pak

b bGCV , k). k = argmin GCV (λ k bGCV ∈Λb λ k

Poznámka. Volbou k = 1 získáme lineární splajny, volbou k = 2 kubické splajny. Ty se v praxi používají převážně. Vyhlazení splajnem pro různý řád 2k a stejný vyhlazovací parametr λ můžeme porovnat z Obr. 3.6 . vyhlazeni splajnem k = 2, λ = 0.0001

vyhlazeni splajnem k = 3, λ = 0.0001 30

20

20

20

10

10

10

0

0

0

y

30

y

y

vyhlazeni splajnem k = 1, λ = 0.0001 30

−10

−10

−10

−20

−20

−20

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

(a) k = 1

0.7

0.8

0.9

1

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

(b) k = 2

0.7

0.8

0.9

1

−30

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

(c) k = 3

Obr. 3.6: Vyhlazení splajnem pro různý řád 2k a stejný vyhlazovací parametr λ

1

60


3.3

Hledání báze prostoru N S 2k (x1, . . . , xn)

Ve větě 3.1.7 jsme předpokládali, že máme vhodnou bázi prostoru N S 2k (x1 , . . . , xn ). Existuje spousta možností, jak můžeme zkonstrulovat tyto báze. Mezi nejznámější z nich patří např. K-W báze (Kimeldorf-Wahba) či D-R báze (Demmler-Reinsch). Mějme funkci −2

q(s, t) = ((k − 1)!)

Z1

k−1 (s − u)k−1 + (t − u)+ du,

0 ≤ s, t ≤ 1,

0

a matici typu n × n

Q = {q(xi , xj )}ni,j=1 ,

která je pozitivně definitní. Dále definujme matici

typu n × k a matici

n,k X = xj−1 i i,j=1 M(λ) = Q + nλIn

typu n × n, kde In je jednotková matice. Pak K-W báze prostoru N S 2k (x1 , . . . , xn ) vypadá takto: ! A , (3.3.1) (g1 (x), . . . , gn (x)) = (1, x, . . . , xk−1 , q(x, x1 ), . . . , q(x, xn )) B kde A = (X T M(λ)−1 X)−1 X T M(λ)−1 a B = M(λ)−1 − M(λ)−1 X(X T M(λ)−1 X)−1 X T M(λ)−1 .

(3.3.2)

Není na první pohled zřejmé, že takto definované funkce tvoří bázi. Musíme tedy tuto vlastnost ověřit. Všimněme si, že po opakované integraci per partes a s využitím binomické věty můžeme psát (q(x, x1 ), . . . , q(x, xn )) = u¯(x)T X T +

(−1)k ((x − x1 )2k−1 , . . . , (x − xn )2k−1 ), + + (2k − 1)!

kde u¯(x) je vektor typu k × 1, ve kterém vystupuje x, ale ne xi . Vynásobíme-li rovnost (3.3.2) zleva maticí X T , dostaneme X T B = X T M(λ)−1 − X T M(λ)−1 X(X T M(λ)−1 X)−1 X T M(λ)−1 = 0. | {z } =Ik

61


Jelikož X T B = 0, platí (q(x, x1 ), . . . , q(x, xn ))B = v¯(x)T X T B + =

(−1)k ((x − x1 )2k−1 , . . . , (x − xn )2k−1 )B + + (2k − 1)!

(−1)k ((x − x1 )2k−1 , . . . , (x − xn )2k−1 )B, + + (2k − 1)!

můžeme napsat 2k−1 (g1 (x), . . . , gn (x)) = (1, x, . . . , xk−1 , (x − x1 )+ , . . . , (x − xn )2k−1 ) +

A (−1)k B (2k−1)!

!

.

(3.3.3)

Funkce 1, x, . . . , xk−1 , (x − x1 )2k−1 , . . . , (x − xn )2k−1 jsou součástí báze pro splajny + + řádu 2k, a jsou tedy lineárně nezávislé. ! Funkce g1 , . . . , gn proto budou lineárně A nezávislé právě tehdy, když matice bude mít hodnost n. Abychom tuto B skutečnost ověřili, všimneme si, že A má hodnost k a B má hodnost n − k a že BM(λ)AT = 0, neboť BM(λ)AT = B M(λ)M(λ)−1 X(X T M(λ)−1 X)−1 = {z } | =In −1

= M(λ) X(X T M(λ)−1 X)−1 −

− M(λ)−1 X(X T M(λ)−1 X)−1 X T M(λ)−1 X(X T M(λ)−1 X)−1 = {z } | =Ik

= 0.

Víme tedy, že g1 (x), . . . , gn (x) jsou lineárně nezávislé. Pokud tedy g1 (x), . . . , gn (x) náleží prostoru N S 2k (x1 , . . . , xn ), pak musí nutně tvořit bázi. Abychom ověřili, že g1 (x), . . . , gn (x) jsou přirozené splajny, stačí ukázat, (k+j) (k+j) že jsou splněny hraniční podmínky gi (0) = 0 = gi (1), j = 0, . . . , k − 1. Napíšeme-li (k+j)

(g1

kde

e (x), . . . , gn(k+j)(x)) = (2k − 1) . . . (k − j)((x − x1 )k−j−1 , . . . , (x − xn )k−j−1 )B, + + (k+j)

e= B

(−1)k B, (2k − 1)!

(k+j)

vidíme, že podmínka gi (0) = 0 je splněna. Pro ověření podmínky gi (1) = 0 k−j−1 použijeme binomickou větu, tedy vyjádříme (1 − xi )+ jako k−j−1 k−j−1 k−j−1 k−j−1 k − j − 1 (1−xi )+ = 1− xi +· · ·+(−1) xk−j−1 . 0 1 k−j−1 i

62


Zavedemeli-li vektor v¯ takto: v¯j = (v0 , v1 , . . . , vk−1 )T =  T k − j − 1 k − j − 1    k−j−1 k − j − 1 = ,− , . . . , (−1) , 0, 0, . . . , 0 , 0 1 k−j−1   | {z } | {z } j

k−j

pak můžeme psát

((1 − x1 )k−j−1 , . . . , (1 − xn )k−j−1 ) = v¯jT X T , + +

j = 0, . . . , k − 1,

a využijeme skutečnosti, že X T B = 0. Nyní jsme ukázali, že funkce g1 (x), . . . , gn (x) definované v (3.3.1) jsou bází prostoru N S 2k (x1 , . . . , xn ). Koeficienty β¯λ definované v (3.1.7) mají při použití této báze velice jednoduchý tvar. Z (3.3.1) a (3.3.3) dostaneme tvar matice G a H:   g1 (x1 ) . . . gn (x1 )  . ..  .. .. G= . .   = g1 (xn ) . . . gn (xn )   ! 1 x1 . . . xk−1 q(x , x ) . . . q(x , x ) 1 1 1 n 1 . . .  A . . . .  . . .. .. .. .. .. = . .  B = 1 xn . . . xk−1 q(xn , x1 ) . . . q(xn , xn ) n A = XA + QB = X|Q B

a

H = (−1)k (2k − 1)!

(−1)k B = B. (2k − 1)!

Proto G + nλH = XA + QB + nλB = XA + M(λ)B = = XA + M(λ)M(λ)−1 − M(λ)M(λ)−1 X (X T M(λ)−1 X)−1 X T M(λ)−1 = | {z } | {z } | {z } =In

=In

= In ,

což znamená, že β¯λ = y¯, neboť (G + nλH) β¯λ = y¯. | {z } =In

Koeficienty u této báze jsou tedy pozorované hodnoty y1 , . . . , yn .

A

63

1.5

1.5

1

1

1

0.5

0.5

0.5

0

−0.5

y

1.5

y

y


0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

−0.5

1

0

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(b) bázové funkcé vynásobené hodnotami bodů pozorování

(a) body s bázovými funkcemi

−0.5

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(c) vyhlazovací splajn

Obr. 3.7: Konstrukce přirozeného splajnu pomocí K-W báze Na Obr. 3.7 je pro lepší představu znázorněna konstrukce vyhlazovacího splajnu pomocí K-W báze. Nyní se seznámíme ještě s tvarem D-R báze, kterou později využijeme při formulaci asymptotického tvaru integrální střední kvadratické chyby odhadu funkce m. 3.3.1 Definice. Nechť f je libovolná funkce a x1 , . . . , xn jsou různá. Poměrnou diferencí k-tého řádu funkce f rozumíme [xi , . . . , xi+k ]f =

[xi+1 , . . . , xi+k ]f − [xi , . . . , xi+k−1 ]f , xi+k − xi

přičemž [xi ]f = f (xi ). 3.3.2 Definice. Nechť f je libovolná funkce a x1 , . . . , xn jsou různá. Normalizovanou poměrnou diferencí k-tého řádu funkce f rozumíme (xi+k − xi )[xi , . . . , xi+k ]f = [xi+1 , . . . , xi+k ]f − [xi , . . . , xi+k−1 ]f . [i]

[i]

Označme α0,k , . . . , αk,k koeficienty u f (xi ) ve vyjádření k-té normalizované diference funkce f . Definujme nyní matici U o rozměrech n × (n − k) tak, že i-tý sloupec matice má tvar [i] [i] (0, . . . , 0, α0,k , . . . , αk,k , 0, . . . , 0), kde i = 1, . . . , n − k, | {z } | {z } i−1

n−k−i

a dále matice V a W takto:

V = (X T Q−1 X)−1 X T Q−1 , W = U T QU.


64

Pro vlastní čísla ζ1 , . . . , ζn matice UW −1 U T platí, že ζ1 , . . . , ζk jsou nulová a ζk+1 ≤ · · · ≤ ζn (viz. [1]). Nechť z¯1 , . . . , z¯n jsou vlastní vektory příslušející vlastním číslům ζ1 , . . . , ζn . Označíme-li Z = [¯ z1 , . . . , z¯n ], pak D-R báze má tvar ! V Z (g1 (x), . . . , gn (x)) = (1, x, . . . , xk−1 , q(x, x1 ), . . . , q(x, xn )) UW −1 U T Z a odhad m(x; b λ) má v této bázi tvar m(x; b λ) =

n X i=1

k n X X z¯iT y¯ z¯iT y¯ T gi (x) = z¯i y¯gi (x) + gi (x). 1 + nλζi 1 + nλζi i=1 i=k+1

Jednotlivé kroky pro odvození tvaru D-R báze jsou podrobněji popsány v [1].

65


3.4

Asymptotické vlastnosti vyhlazovacích splajnů

Věta 3.4.1. Nechť jsou splněny následující předpoklady: 1. Vyhlazovací parametr λ = λn splňuje následující podmínky limn→∞ λ = 0, limn→∞ nλ = ∞. 2. Předpokládejme, že pro body plánu platí: xi = ni , i = 1, . . . , n. 3. Funkce m splňuje m ∈ C k [0, 1], tj. m má spojité derivace až do řádu k včetně. Pak pro hlavní člen MISE(m(λ)) b platí

MISE(m(λ)) b ≤ MISE max (m(λ)), b

kde

kde

MISE max (m(λ)) b =λ V∗ =

Z∞

Z1

(m(k) (x))2 dx +

0

σ2 V ∗ 1

nπλ 2k

,

1 dx. (1 + x2k )2

−∞

Důkaz. Integrální střední kvadratickou chybu lze vyjádřit ve tvaru MISE(m(λ)) b = =

Z1 0

Z1 0

E(m(x) − m(x; b λ))2 dx = 2

(m(x) − E m(x; b λ)) dx + | {z } 2 (vychýlení m(x;λ)) b

kde E m(x; b λ) je minimalizátor −1

n

n X

2

(m(x) − f (xi )) + λ

i=1

Z1 0

Z1

(E m b 2 (x; λ) − E 2 m(x; b λ)) dx, | {z } rozptyl m(x;λ) b

(f (k) (x))2 dx

0

pro f (x) ∈ W2k [0, 1]. Pro n → ∞ platí E m(x; b λ) = fλ (x), kde fλ (x) je minimalizátor Z1 Z1 2 (m(x) − f (x)) dx + λ (f (k) (x))2 dx 0

0

66


pro f (x) ∈ W2k [0, 1]. Pro k ≥ 2 platí (viz. [1]) Z1 0

2

(vychýlení m(x; b λ)) dx =

Z1

(m(x) − fλ (x))2 dx(1 + o(1)),

0

a označme tedy hlavní část Z1 0

Jelikož platí Z1 0

2

(vychýlení m(x; b λ)) dx =

2

b λ)) dx ≤ (vychýlení m(x;

Z1

Z1

(m(x) − fλ (x))2 dx.

0

2

(m(x) − f (x)) dx + λ

0

Z1

(f (k) (x))2 dx,

0

pro libovolnou funkci f (x) ∈ W2k [0, 1], a položíme-li f (x) = m(x), dostaneme Z1 0

2

(vychýlení m(x; b λ)) dx ≤ λ

Z1

(m(k) (x))2 dx.

0

Nyní se podívejme na rozptyl m(x; b λ). Připomeňme, že vyjádříme-li odhad m(x; b λ) v D-R bázi, obdržíme n X z¯iT y¯ m(x; b λ) = gi (x). 1 + nλζi i=1 Jelikož systém bázových funkcí g1 (x), . . . , gn (x) je ortonormální, pak pro kovarianci z¯jT y¯ a z¯jT′ y¯ platí C(¯ zjT y¯, z¯jT′ y¯)

=σ

2

n X i=1

a tedy

( σ2 gj (xi )gj ′ (xi ) = 0

′

j=j ′ , j 6= j

n X gi2 (x) rozptyl m(x; b λ) = σ . 1 + nλζi i=1 2

Nyní Z1 0

rozptyl m(x; b λ) dx == σ 2

n X i=1

R1

gi2 (x) dx

0

1 + nλζi

.

67


Pro dostatečně velké n platí n

1 . 1X 2 . = gj (xi ) = n n i=1

Z1

gj2 (x) dx

0

(viz. [1]), a proto můžeme napsat Z1 0

Nyní

n 1 σ2 X rozptyl m(x; b λ) dx = + O(n−1 ). n i=1 1 + nλζi

n σ2 X 1 σ2 X 1 σ2 X 1 = + , n i=1 1 + nλζi n 1 + nλζi n 1 + nλζi 3 3 i≤n 4k

i>n 4k

přičemž podle [1] 1 1 X . = O(λ−2 n−3 ) = o(n−1), n 1 + nλζi 3 i>n 4k

1 X 1 . V∗ = 1 , n 1 + nλζi nπλ 2k 3 i≤n 4k

kde ∗

V =

Z∞

1 dx. (1 + x2k )2

−∞

Hlavní část

R1 0

rozptyl m(x; b λ) dx je Z1 0

a tedy

rozptyl m(x; b λ) dx =

MISE max (m(x; b λ)) = λ

Z1

V∗ 1

nπλ 2k

,

(m(k) (x))2 dx +

0

σ2V ∗ 1

nπλ 2k

.

Z věty 3.4.1 plyne, že integrální střední kvadratická chyba odhadu m b je zhora ohraničena hodnotou MISE max (m(λ)) b =λ

Z1 0

(m(k) (x))2 dx +

σ2 V ∗ 1

nπλ 2k

.

68


Jelikož neznáme přesnou hodnotu MISE(m(λ)), b nemůžeme ani přesně spočítat hodnotu λopt , která MISE(m(λ)) b minimalizuje. Místo λopt budeme používat hode b notu λopt , která minimalizuje MISE max (m(λ)).

Věta 3.4.2. Nechť jsou splněny předpoklady věty 3.4.1. Pak přibližná hodnota asymptoticky optimálního vyhlazovacího parametru je dána vztahem 2k 2 ∗ 2k+1 σ V e λopt = , 2knπDk∗

kde

∗

V =

Z∞

1 dx (1 + x2k )2

Dk∗

a

=

−∞

Z1

[m(k) (x)]2 dx.

0

eopt minimalizuje MISE max (m(λ)), Důkaz. Jelikož λ b musí platit b ∂MISE max (m(λ)) e = 0. ∂λ λ=λopt

Tedy

Nyní z rovnosti

σ2V ∗ ∂MISE max (m(λ)) b ∗ = Dk − 1+2k . ∂λ 2knπλ 2k 0 = Dk∗ −

získáme eopt = λ

Důsledek 3.4.3.

σ2V ∗

2knπλ

σ2 V ∗ 2knπDk∗

1+2k 2k

2k 2k+1

.

2k 1 2k 1 + 2k eopt )) = n− 2k+1 MISE max (m( b λ (Dk∗ ) 2k+1 σ 2 V ∗ 2k+1 , 2k (2kπ) 2k+1 2k

tj. asymptotická rychlost konvergence MISE max → 0 je n− 2k+1 . Důkaz. eopt )) = MISE max (m( b λ 2k 2 ∗ 2k+1 σ V σ2 V ∗ ∗ = D + = 1 k 2k+1 2knπDk∗ σ2 V ∗ nπ 2knπD∗ k

=

1

2

∗

2k 2k+1

1

(Dk∗ ) 2k+1

1

2k 1 (2kπn) 2k+1 2 ∗ 2k+1 + σ V (Dk∗ ) 2k+1 = πn

+ σ V 2k (2kπn) 2k+1 2k 1 2k 1 + 2k = n− 2k+1 (Dk∗ ) 2k+1 σ 2 V ∗ 2k+1 . 2k (2kπ) 2k+1

69


eopt platí Důsledek 3.4.4. Pro λ

eopt ))2 1 (vychýlení m(x; b λ = . eopt ) 2k rozptyl m(x; b λ

Důkaz. Plyne z důkazu předchozí věty, neboť:

=

eopt )) = MISE max (m( b λ 1

2

2k

(2kπn) 2k+1 |

σ V

∗

{z

2k 2k+1

1

1

eopt ))2 (vychýlení m(x; b λ

a

eopt ) rozptyl m(x; b λ

1

1 2k

(2kπn) 2k+1 kde V∗ =

Z∞

−∞

2k 1 (2kπn) 2k+1 2 ∗ 2k+1 + σ V (Dk∗ ) 2k+1 πn } | {z }

(Dk∗ ) 2k+1

1 dx (1 + x2k )2

1 (2kπn) 2k+1 = , : πn 2k

a

Dk∗ =

Z1 0

[m(k) (x)]2 dx.

Kapitola 4 Simulační studie Cílem této kapitoly bude ověřit některé vlastnosti jádrových odhadů a vyhlazovacích splajnů, zejména vliv šířky vyhlazovacího okna h, resp. vyhlazovacího parametru λ. Testovací data si připravíme následujícím způsobem: Zvolíme si dostatečně hladkou funkci m, u níž máme zaručeno, že náleží do prostoru W 2 [0, 1]. Vezmeme soubor dat o n vzorcích ležících na intervalu [0, 1] a k hodnotě každého bodu přičteme bílý šum s rozložením N(0, σ 2 ). Máme tedy dvojice (x1 , y1 ), . . . , (xn , yn ), pro něž navíc známe i regresní funkci m. Na Obr. 4.1 jsou simulovaná data a regresní funkci. regresni funkce m1 40

30

20

y

10

0

−10

−20

−30

0

0.2

0.4

0.6

0.8

1

x

Obr. 4.1: Simulovaná data s regresní funkcí m1 (x)

70

71

KAPITOLA 4. SIMULAČNÍ STUDIE

Pro následující studii byly zvoleny tyto parametry: m1 (x) = 6 cos(10x + 12 ) + 3 cos(20x + 1) + 32 , n = 200, σ 2 = 10. Jelikož nás v této části bude zajímat zejména vliv šířky vyhlazovacího okna h, resp. vyhlazovacího parametru λ na kvalitu odhadu, zvolíme si pevně jádrový odhad a jádro, včetně řádu, a také řád vyhlazovacího splajnu. V našem případě byl použit Nadarya-Watsonův odhad s Epanečnikovým jádrem (k = 2, ν = 0) a kubický splajn (k = 2). Nejprve se zaměřme na nalezení vyhlazovacích parametrů h a λ pomocí zobecněné metody křížového odhadu. Parametr h budeme hledat na intervalu [0, 1] a parametr λ na intervalu [10−16 , 102]. Na Obr. 4.2 vidíme průběhy funkce GCV.

funkce m1, prubeh GCV

funkce m1, prubeh GCV

135

200

130 180 125 160 GCV(λ)

GCV(h)

120 115

140

110 120 105 100 100 95

0

0.2

0.4

0.6

0.8

1

80 −16

−14

−12

h

(a) jádrové odhady

−10

−8

−6 log(λ)

−4

−2

0

2

(b) vyhlazovací splajny

Obr. 4.2: Průběh funkce GCV Z obrázku je zřejmé, že funkce GCV nabývá pro jádrový odhad svého minima v h = 0, 0483 a pro vyhlazovací splajn v λ = 3, 1 · 10−7 . Jelikož h = 0, 0483 a n = 200, není obtížné si spočítat, že jádro zasahuje přibližně do 19 datových bodů. Na přůběhu funkce GCV (λ) pro vyhlazovací splajn můžeme zřetelně rozeznat následující tři oblasti: 1. interval [10−16 , 10−12], kde hodnota GCV (λ) se téměř nemění a pro odpovídající λ je odhad m(x; b λ) interpolačním splajnem,

72


2. interval [10−12 , 10−2 ], kde změna parametru λ probíhá ve prospěch hladkosti či záchování věrnosti dat, v této oblasti se také nachází minimum funkce GCV (λ), jemuž odpovídá λGCV , 3. interval [10−2, 102 ], kde hodnota GCV (λ) se téměř nemění a pro odpovídající λ vyhlazovací splajn přechází v regresní polynom. Pro Obr. 4.3 jsou použity zjištěné hodnoty hGCV a λGCV . Srovnáváme zde odhady e m(x; b hGCV ) a m(x; b λGCV ) s odhady m(x; b hopt ) a m(x; b λopt ), kde hodnoty hopt a e λopt jsme odvodili z asymptotických vlastností jádrových odhadů a vyhlazovacích splajnů. srovnani odhadu s regresni funkci m1

srovnani odhadu s regresni funkci m1

40

40 data m1(x) m1(x;h

30

data m1(x) m1(x;λ

30

)

opt

)

opt

m1(x;hGCV)

m1(x;λ

)

GCV

10

10 y

20

y

20

0

0

−10

−10

−20

−20

−30

0

0.2

0.4

0.6 x


0.8

1

−30

0

0.2

0.4

0.6

0.8

1

x


Obr. 4.3: Srovnání regresní funkce m1 (x) s jejími odhady s optimálním parametrem a s parametrem získaným GCV funkcí Všimněme si, že pro zvolené parametry se podařilo, jak u jádrového odhadu, tak u vyhlazovacího splajnu, najít zobecněnou metodou křížového ověřování hodnoty hGCV a λGCV velice blízké hodnotám hopt a e λopt . Jádrové i splajnové odhady regresní funkce jsou si v obou případech relativně podobné, přičemž v případě splajnu je výsledek hladší. Dále si všimněme podobnosti odhadů se samotnou regresní funkcí. Oba odhady velice věrně napodobují její průběh. Tedy, přestože data jsou značně zašuměná, můžeme přesto za daných podmínek při dostatečném počtu dat (n = 200) relativně dobře odhadnout průběh regresní funkce. Fakt, že zjištěné hodnoty hGCV a λGCV jsou v našem případě skutečně nejideálnější, lze posoudit na Obr. 4.4, kde odhad m(x; b hGCV ), resp. m(x; b λGCV ) porovnáváme s odhady m(x; b 5hGCV ) a m(x; b 5−1 hGCV ), resp. s odhady m(x; b 1000λGCV ) a −1 m(x; b 1000 λGCV ).

73


srovnani odhadu regresni funkce m1 pro ruzne hodnoty λ

srovnani odhadu regresni funkce m1 pro ruzne hodnoty h

40

40 data m1(x) m1(x;h

30

data m1(x) m1(x;λ

30

)

GCV

m1(x;h

20

)

GCV 3

m1(x;10 λ

m1(x;5hGCV )

GCV

/5)

20

GCV

−3

m1(x; 10

λ

) )

GCV

10

y

y

10

0

0

−10

−10

−20

−20

−30

0

0.2

0.4

0.6

0.8

1

−30

0

x


0.2

0.4

0.6

0.8

1

x


Obr. 4.4: Odhady funkce m1 (x) pro různé parametry Je zřejmé, že odhady s příliš vysokými hodnotami parametru, tj. m(x; b 5hopt ) a e m(x; b 1000λopt), jsou přehlazené, zatímco odhady s příliš nízkými hodnotami parametru, tj. m(x; b 5−1 hopt ) a m(x; b 1000−1 e λopt ), jsou podhlazené.

V druhé části simulační studie se budeme věnovat podrobněji hledání hodnot vyhlazovacích parametrů pomocí zobecněné metody křížového ověřování a porovnání výsledků s asymptoticky optimálními hodnotami těchto parametrů. Zároveň budeme zkoumat, jak se mění odhady parametrů vyhlazení při různé úrovni zašumění regresní funkce m2 . Pro tuto studii byly zvoleny následující parametry: m2 (x) = 5 cos(10x − π2 ) + 5 cos(18x), n = 100, σ 2 = 0, 8 a σ 2 = 0, 08.

Nejprve se podívejme na případ, kdy σ 2 = 0, 8. Data jsou zašuměná velice výrazně a pouhým okem v nich nelze rozeznat průběh regresní funkce m2 . Při vyhlazení eopt jsme schopni se regresní funkci přiblížit, ale ne do takové s hodnotami hopt a λ míry, jako v předchozím případě. To může být způsobeno skutečností, že zvolených 100 vzorků plánu ještě není dostatečné množství pro takto zašuměná data. Pokud bychom však chtěli zjistit hodnotu hGCV a λGCV , došli bychom k závěru, že v tomto případě zobecněná metoda křížového ověřování selhává, jak je vidět na Obr. 4.5, neboť nelze určit minimum funkce GCV .

74


funkce m2, ε ~ N(0,0.8), prubeh GCV

funkce m2, ε ~ N(0,0.8), prubeh GCV 1.15

2.2

1.1

2

1.05

1.8

1

GCV(λ)

GCV(h)

1.6 0.95

1.4 0.9

1.2 0.85

1

0.8

0.75

0

0.1

0.2

0.3

0.4

0.5 h

0.6

0.7

0.8

0.9

0.8 −16

1

−14

−12

−10

−8

−6

−4

−2

0

2

log(λ)



Obr. 4.5: Průběh funkce GCV Při pohledu na výsledek vyhlazení (viz. Obr. 4.6) je jasné, že odhady funkce m2 jsou příliš přehlazené a jakákoliv struktura dat se ztrácí. funkce m2(x), ε ~ N(0,0.8), srovnani odhadu s regresni funkci m2

funkce m2(x), ε ~ N(0,0.8), srovnani odhadu s regresni funkci m2

2.5

2.5 data m2(x) m2(x;h ) opt m2(x;h )

2

data m2(x) m2(x;λ ) opt m2(x;λ )

2

GCV

GCV

1

1

0.5

0.5 y

1.5

y

1.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5

−2

0

0.1

0.2

0.3

0.4

0.5 x

0.6


0.7

0.8

0.9

1

−2

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9


Obr. 4.6: Srovnání regresní funkce m2 (x) s jejími odhady s optimálním parametrem a s parametrem získaným GCV funkcí eopt )) a Srovnáním chyb MISE(m(h b opt )) a MISE(m(h b GCV )), resp. MISE(m( b λ MISE(m(λ b GCV )) vidíme, že se skutečně nejedná o kvalitní odhady, neboť je patrné, že jejich hodnoty se od sebe liší dvěma desetinými řády v případě jádrových odhadů (viz. Tab. 4.9) a třema desetinými řády v případě splajnového vyhlazování (viz. Tab. 4.10). Při dané úrovni šumu tedy nejsme schopni sestrojit věrný odhad regresní funkce.

1

75


Nyní se podívejme na případ, kdy σ 2 = 0, 08. Data nejsou zašuměná příliš a pouhým okem lze bez problémů rozpoznat jejich strukturu. Při vyhlazování s hodeopt oba odhady průběh funkce m2 relativně věrně napodobují. notami hopt a λ Pokusíme-li se opět najít hodnoty hGCV a λGCV pomocí zobecněné metody křížového ověřování, vidíme z Obr. 4.7, že v tomto případě je vše v pořádku, neboť obě eopt . GCV funkce mají minimum, které leží blízko vypočteným hodnotám hopt a λ funkce m2, ε ~ N(0,0.08), prubeh GCV

0.016

0.014

0.014

0.012

0.012

GCV(λ)

GCV(h)

funkce m2, ε ~ N(0,0.08), prubeh GCV

0.016

0.01

0.01

0.008

0.008

0.006

0.006

0.004

0

0.1

0.2

0.3

0.4

0.5 h

0.6

0.7

0.8

0.9

0.004 −16

1

−14

−12

−10

−8

−6

−4

−2

0

2

log(λ)



Obr. 4.7: Průběh funkce GCV Při pohledu na výsledek vyhlazení (viz. Obr. 4.8) je vidět, že odhady funkce m2 jsou velmi dobré a struktura dat je zachována. funkce m2(x), ε ~ N(0,0.08), srovnani odhadu s regresni funkci m2

funkce m2(x), ε ~ N(0,0.08), srovnani odhadu s regresni funkci m2

0.3

0.3 data m2(x) m2(x;h ) opt m2(x;h )

0.2

data m2(x) m2(x;λ ) opt m2(x;λ )

0.2

GCV

0

0 y

0.1

y

0.1

GCV

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3

−0.4

0

0.1

0.2

0.3

0.4

0.5 x

0.6


0.7

0.8

0.9

1

−0.4

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8


Obr. 4.8: Srovnání regresní funkce m2 (x) s jejími odhady s optimálním parametrem a s parametrem získaným GCV funkcí

0.9

1

76


eopt )) a Srovnáním chyb MISE(m(h b opt )) a MISE(m(h b GCV )), resp. MISE(m( b λ MISE(m(λ b GCV )) vidíme, že se jedná o kvalitní odhady, neboť je patrné, že rozdíly v jejich hodnotách jsou téměř zanedbatelné, jak u jádrových odhadů, tak v případě splajnového vyhlazování (viz. Tab. 4.9 a Tab. 4.10). Při dané úrovni šumu tedy jsme naopak schopni sestrojit velmi věrný odhad regresní funkce. Nyní shrňme získané výsledky do tabulek. σ2

hopt

0.8

0.176

0.7

0.1675 0.9

hGCV 0.9

MISE max (m(h b opt )) MISE(m(h b GCV )) 0.0136

2.787

0.0110

1.8287

0.6

0.2503 0.1575 0.0086

0.0152

0.5

0.1464 0.1065 0.0064

0.0074

0.4

0.1464 0.1602 0.0064

0.0065

0.3

0.1193 0.0963 0.0028

0.0030

0.2

0.1015 0.0999 0.0015

0.0015

0.1

0.0769 0.0893 5.1270·10−4

4.871·10−4

0.08

0.0703 0.0671 3.04·10−4

3.423·10−4

Tab. 4.9: Hodnoty pro jádrové vyhlazování s různým šumem

σ2

eopt λ

MISE max (m(λ b opt )) MISE(m(λ b GCV ))

λGCV

0.8

4.061·10−6

1.1·10−2

0.0061

5.8186

0.7

3.6495·10−6

1.1·10−2

0.0052

6.4005

0.6

3.2261·10−6

8.2309·10−5

0.5

2.7882·10

−6

0.0041

0.0489

−5

0.0032

0.0121

0.4

2.3324·10−6

9.7082·10−5

0.0024

0.0569

0.3

1.8529·10−6

8.412·10−5

0.0017

0.0492

0.2

1.3396·10−6

1.8926·10−6

0.0011

0.0014

0.1

7.694·10−7

8.871·10−7

5.3719·10−4

6.0254·10−4

0.08

6.436·10−7

1.121·10−6

4.34·10−4

6.355·10−4

1.9122·10

Tab. 4.10: Hodnoty pro splajnové vyhlazování s různým šumem


77

Z tabulek je patrné, že pro naše zvolené parametry (funkce m2 a n = 100) kovergují hodnoty parametrů nalezené zobecněnou metodou křížového ověřování snižováním úrovně šumu k asymptoticky optimálním hodnotám. Vidíme také, že při klesající úrovni šumu roste kvalita odhadu, jelikož hodnota integrální střední kvadratické chyby klesá k nule. Pro funkci m2 jsme zjistili, že úroveň šumu, pro níž začínají být odhady pomocí zobecněné metody křížového ověřování použitelné, je σ 2 = 0, 6. Dále si všiměme drobného kolísání hodnot v okolí σ 2 = 0, 4. Toto kolísání je s největší pravděpodobností zapříčiněno skutečností, že v těchto případech se u počítačového generátoru náhodných čísel s normálním rozložením projevila odchylka od normality pro n = 100, která by se pro větší n samozřejmě minimalizovala.

Závěr Tato práce představila dva základní přístupy k neparametrické regresi, konkrétně jádrové odhady a vyhlazovací splajny. Předvedli jsme si, jakým způsobem je možné tyto metody zavést a formulovat v kontextu regresní analýzy. Ke každé z metod jsme představili jejich nejpoužívanější varianty a seznámili jsme se s parametry, které ovlivňují výsledky těchto metod. Také jsme si odvodili postupy, pomocí nichž lze optimální hodnoty těchto parametrů přibližně najít. Praktická část této diplomové práce se zabývala implementací jádrových metod a vyhlazovacích splajnů v jazyce MATLAB. Pro testování těchto metod byla vygenerována simulovaná data, u nichž jsme ideální regresní funkci přímo znali, a byli jsme tedy schopni provést objektivní srovnání s jejím odhadem.

78

Literatura [1] Eubank, Randall L.: Spline smoothing and nonparametric regression. New York: Marcel-Dekker, 1988. [2] Horová, Ivanka: přednášky k předmětu Neparametrické vyhlazování, 2008. Masarykova univerzita. Nepublikováno. [3] Plch, Roman - Čechová, Lenka: Sázíme diplomovou práci z matematiky v TEXu. Brno: Masarykova univerzita, 2003. htpp://www.math.muni.cz/ ~plch/vyuka/b.ps. [4] Rybička, Jiří: LATEX pro začátečníky. Brno: Konvoj, 1999. [5] Wand, M. P. - Jones, M. C.: Kernel smoothing. London: Chapman & Hall, 1995.

79

MASARYKOVA UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA ÚSTAV MATEMATIKY A STATISTIKY DIPLOMOVÁ PRÁCE VYHLAZOVÁNÍ A REGRESE

Recommend Documents