Zvláštní poděkování patří panu Profesoru doc. RNDr. Janu Pickovi, CSc. S nímž jsem konzultoval práci. Rád bych také poděkoval panu Mgr. Martinu Schindlerovi, Ph.D. Nesmím opomenout ani ty, kteří mě v době mého studia nejvíce podporovali: „ I Vám děkuji.”
ANOTACE Tato práce je určena pro učitele matematiky, ale i pro zájemce z řad široké veřejnosti. Cílem práce je prostřednictvím statistických analytických postupů získat poznatky o analyzovaných časových řadách a porovnat výsledky pro časové řady naměřené na jednotlivých meteorologických stanicích prostřednictvím programů Statgraphics Centurion XV. a Matlab R2009b. Text je doplněn názornými obrázky a tabulkami. Práce byla vytvořena programem TeXnicCenter.
Klíčová slova časové řady, distribuční funkce, nulová hypotéza, alternativní hypotéza, úhrny srážek, těžké chvosty, lehké chvosty, sféra přitažlivosti
ANNOTATION This thesis is focused mainly on mathematics teachers, but also on general public. Its target is to get knowledge about analyzed time series by methods of statistical analysis and to compare the results of time series captured at various climatological stations with Statgraphics Centurion XV. and Matlab R2009b applications. The text is complemented with illustrative pictures and tables. The thesis was created in TeXnicCenter.
Key words time series, distribution function, null hypothesis, alternative hypothesis, precipitation, heavy tails, light tails, domain of attraction
Obsah 1 Úvod
5
2 Typy stanic na území České republiky
7
2.1
Automatizovaná stanice s profesionální obsluhou (AMS) . . . . . . . . . .
7
2.2
Automatizované stanice s dobrovolnou obsluhou (AKS) . . . . . . . . . . .
8
2.3
Automatizované srážkoměrné stanice (ASS) . . . . . . . . . . . . . . . . .
8
2.4
Manuální klimatologická stanice (MKS) . . . . . . . . . . . . . . . . . . . .
8
2.5
Manuální srážkoměrné stanice(MSS) . . . . . . . . . . . . . . . . . . . . .
9
2.6
Totalizátory (TOTAL) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3 Grafické znázornění klimatologických časových řad
11
4 Statistická hypotéza
19
4.1
Chyby testu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5 Analýza rozptylu
21
5.1
Příklad 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2
Ověření normality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.3
Příklad 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
6 Hodnoty extrémních distribucí
37
7 Závěr
46
Seznam zkratek a použitých symbolů ANOVA analýza rozptylu atd.
a tak dále
cm
centimetr
d.f.
distribuční funkce
H0
nulová hypotéza
HA
alternativní hypotéza
max
maximálně
min
minimálně
mm
milimetr
Mn
maximální hodnota náhodného výběru
m n. m.
metr nad mořem
n
přirozené číslo nabývá hodnot 1,2,3,5,7
sup
suprémum
T
dobu, za kterou v daných letech překročí naměřená data určitou hranici
tzv.
takzvaný
x
aritmetický průměr
γ
extremální index
µ
stupeň volnosti
µ, EX
střední hodnota
σ2
rozptyl
1
Úvod
Cílem této bakalářské práce je porovnávání naměřených dat z jednotlivých meteorologických a klimatologických stanic v České republice za zvolené období. Porovnávají se zde jednodenní, dvoudenní, třídenní, pětidenní a sedmidenní srážkové úhrny za sledovaná léta, ve kterých byla provedena měření. Přihlíží se zde k nadmořské výšce. Ve většině případů jsou stanice seřazeny vzestupně od nejníže položeného místa k místu s nejvyšší nadmořskou výškou. Pro představu porovnání stanic naměřených dat slouží tabulka 1. Tabulka byla vytvořena pomocí dat stanic Churáňov, Lysá Hora, Praděd. Data jsou seřazena vzestupně podle nadmořské výšky, která je též v tabulce uvedena. Nejvíce naměřených srážek je ze stanice Lysá hora. Maximální hodnoty jednodenního úhrnu srážek zde dosahují hodnoty 233.80 mm. Průměrná hodnota převyšuje 80 mm v jednodenním srážkovém úhrnu. Důležité je, že stanice Lysá Hora má oproti stanici Praděd menší nadmořskou výšku. Rozdíl mezi těmito stanicemi je 168m. Což z hlediska nadmořské výšky u nás není zanedbávaný rozdíl. Města
m n. m.
max
min
průměr
Churáňov
1118,00
135,30
22,30
51,47
Lysá Hora
1322,00
233,80
30,50
82,67
Praděd
1490,00
139,40
31,00
59,72
Tabulka 1: Statistika jednodenních úhrnů srážek stanic Churáňov, Lysá Hora, Praděd.
Graf 1 popisuje úhrn srážek u stanic Praděd, Lysá Hora, Chůráňov v jednodenním srážkovém úhrnu. Podle grafu je zřejmé, že data, která byla v letech 1961-2007 naměřena ve stanici Lysá Hora, výrazně převyšují data ostatních stanic. Povšimněme si také let 19951998, kde podle grafu strmě roste jednodenní úhrn srážek pro data Lysá Hora. V těchto letech byly některé oblasti Moravy těžce zasaženy povodněmi. Provedením analýzy a zvolením vhodného statistického softwaru provádím výpočty. 5
2D Graph 250 Praded Lysa Hora Churanov
precipitation [mm]
200
150
100
50
0 1960
1970
1980
1990
2000
2010
years
Obrázek 1: Graf jednodenního úhrnu srážek stanic s nadmořskou výškou nad 1000m. (zdroj vlastní) Dle těchto grafů a propočtů lze optimálně předpovědět co nejpřesněji pravděpodobnost povodní v nejohroženějších lokalitách České republiky a doporučit vybudování protipovodňových opatření. Následující kapitola se zabývá typy stanic, v další kapitole se pomocí grafického znázornění klimatologických časových řad popisuje chování úhrnů srážek v některých stanicích v České republice. Nakonec se vysvětlí pojem statistická hypotéza, díky němuž se následně provedou statistické testy a odhady extremální hodnoty.
6
2
Typy stanic na území České republiky
Na území Čech, Moravy a Slezska se nachází šest typů stanic. Prostřednictvím nich se získávají srážková a jiná data, která slouží k analýze počasí na našem území. [1]
Obrázek 2: Příklad automatizované stanice (Zdroj: http://www.campbellsci.com).
2.1
Automatizovaná stanice s profesionální obsluhou (AMS)
Automatizovaná meteorologické stanice jsou obsluhovány pozorovateli, kteří jsou zaměstnanci ČHMÚ a jejich provoz a metodické řízení zajišťuje odbor profesionální staniční sítě (OPSS), odbor letecké meteorologie (OLM) ČHMÚ v Praze, popř. Jsou pod správou Armády ČR (AČR). Program klimatologických pozorování je stejný jako u klimatologických automatizovaných stanic, ale jejich hlavním úkolem jsou synoptická pozorování. Tato pozorování probíhají každou hodinu a jsou okamžitě předávána v kódované zprávě SYNOP do centra v Praze-Komořanech. Data z některých synoptických stanic jdou pak i do celosvětového systému výměny meteorologických dat. Kromě prvků, které jsou pozorovány na klimatologických stanicích, se v synoptickém pozorování sleduje i výška a druh oblačnosti, atmosférický tlak vzduchu, tlaková tendence, dohlednost, podobně se rozlišuje typ padajících srážek atd. 7
2.2
Automatizované stanice s dobrovolnou obsluhou (AKS)
Tyto stanice měří všechny meteorologické prvky jako stanice manuální. Měření se však provádí pomocí automatických přístrojů (čidel) s intervalem záznamu 10 minut, u srážek s intervalem záznamu 1 minuta. Stanice jsou vybaveny modemem, který přenáší naměřená data pomocí systému GPRS každých 5 minut na centrální počítač na pobočku. Pozorovatel je povinen sledovat a zapisovat atmosférické jevy v průběhu dne a v zimním období měřit sněhové charakteristiky klasickými přístroji.
2.3
Automatizované srážkoměrné stanice (ASS)
Stanice jsou vybaveny automatickým srážkoměrem MR3H-FC a jejich data jsou přenášena pomocí GPRS do centrálního počítače na pobočce.
Obrázek 3: přístroj MR3H-FC (Zdroj: http://www.meteoservis.cz)
Interval záznamu je 1 minuta, přenos dat probíhá každých 5 minut. Měření sněhových charakteristik a pozorování meteorologických jevů se provádí stejně jako u stanic manuálních.
2.4
Manuální klimatologická stanice (MKS)
Ve třech klimatologických termínech v 7, 14, 21 hodin SEČ se provádí měření základních meteorologických prvků teplota vzduchu, maximální a minimální teplota vzduchu, 8
přízemní minimální teplota vzduchu, relativní vlhkost vzduchu, směr a vlhkost větru. Pozorovatel zaznamenává oblačnost, stav počasí a stav půdy. V termínu 7 hodin měří úhrn srážek, výšku nového sněhu, celkovou výšku sněhové pokrývky a vodní hodnotu sněhové pokrývky. Sleduje výskyt atmosférických jevů v průběhu celého dne. Na některých stanicích se dále měří také sluneční svit, teploty půdy, popř. Výpar nebo atmosférický tlak vzduchu. Data jsou podobně jako u srážkoměrných stanic zapisovaná do klimatologických výkazů, popřípadě v souborech na disketě nebo e-mailem.
2.5
Manuální srážkoměrné stanice(MSS)
V klimatologickém termínu 7 hodin středoevropského času pozorovatel provádí měření denního úhrnu srážek a sněhových charakteristik, tj. Výška nového sněhu, celková výška sněhové pokrývky a vodní hodnota celkové sněhové pokrývky. V průběhu celého dne pozorovatel zaznamenává výskyt a průběh meteorologických jevů. Data jsou zapisována do měsíčních srážkoměrných výkazů a zasílána po skončení měsíce poštou nebo vkládána do pořizovacího programu na PC a zasílána elektronickou poštou, popřípadě na disketě.
Obrázek 4: Stanice MSS (zdroj http://www.tobaccovilleweather.com )
9
• Nahoře: Vyhřívací srážkový kolektor. • Uprostřed: Senzor propojovacího modulu s vysílačem. • Dole: Ventilátor odsávajícího radiačního štítu, který zachycuje venkovní teplotu a vlhkost vzduchu sondy.
• Vzadu: 4 ”(0,01”rozlišení) MSS.
2.6
Totalizátory (TOTAL)
Mezi doplňkové stanice můžeme zařadit rovněž zvláštní srážkoměrná zařízení, zvaná totalizátory. Tyto přístroje se umísťují do horských nebo špatně dostupných oblastí a tam, kde nelze zajistit pravidelné měření srážek. Měření se provádí zpravidla dvakrát ročně, provádí je pracovníci oddělení, a tyto „ půlroční” úhrny jsou pomocí referenčních srážkoměrných stanic přepočítávány na měsíční úhrny srážek.
Obrázek
5:
Totalizátor
(zdroj:
http://www.tisicovky.cz/cs/hory/moravskoslezskej-
beskydy/norici-hora-hlv268/foto/)
10
3
Grafické znázornění klimatologických časových řad Moravia 1961−2007
160 140 120
precipitation [mm]
200
100
150
2010
1990
50 0 1
80
2000
100
1980 2
60 40
1970 3
4
years 5
6
20
1960
towns
Obrázek 6: 3D Graf jednodenního úhrnu srážek v letech 1961-2007. Stanice 1-Brno, 2Olomouc, 3-Štěpánov, 4-Kunčice, 5-Tovačov, 6-Klášterní Hradisko (Zdroj: Vlastní)
Velmi vysoké srážky v letech 1961-2007 na Moravě zachycuje trojrozměrný graf. Stanice Kunčice (MKS) již v jednodenním období výrazně převyšuje úhrny srážek v ostatních stanicích. Je třeba se podrobněji zabývat tím, zda takový rozdíl nastal pouze v roce 1997. Byl tedy překvapujícím jevem, který zaskočil tuto oblast, nebo se již v předchozích letech objevovaly podobné rozdíly úhrnů srážek? Bylo by možné se na tyto situace předem připravit?
11
Moravia 1961−1997 240 220 200 250 precipitation [mm]
180 200
160
150
140
100
120
50
100 6
0 2000
80
5 1990
4 3
1980 years
40
2
1970 1960
1
60
towns
Obrázek 7: 3D graf pětidenního úhrnu srážek v letech 1961-1996. Stanice 1-Brno, 2Olomouc, 3-Štěpánov, 4-Kunčice, 5-Tovačov, 6-Klášterní Hradisko (Zdroj: Vlastní) Za zmínku stojí rok 1977, kdy v Kunčicích spadlo během pěti dnů 250,20 mm srážek Z pouhého grafu 7 však nelze vyčíst žádná závislost jednotlivých let, ba dokonce závislost úhrnu srážek mezi jednotlivými městy. K tomu, aby se některé závislosti prokázaly je zapotřebí použití statistických testů. Zajímavější může být porovnání dvou stanic s naprosto odlišnými nadmořskými výškami. Graf 8 charakterizuje jednodenní úhrny srážek pro stanici Děčín (MKS) a Filpovou Huť (ASS). Ve stanici s vyšší nadmořskou výškou, Filipově Huti napadlo v letech 19682007 více srážek než v nejníže položené stanici Děčín. V závorkách vedle názvu stanic jsou uvedeny typy stanic, ze kterých byla data pořízena.
12
2D Graph 120 Filipova Hut Decin (Hrensko)
precipitation [mm]
100
80
60
40
20
0 1965
1970
1975
1980
1985 1990 years
1995
2000
2005
2010
Obrázek 8: 2D graf jednodenního úhrnu srážek měst Filipova Huť a Děčín (Hřensko)(Zdroj: Vlastní)
Na dalším grafu 9 jsou zaznamenány podíly srážek v létě a v zimě v letech 1968-2007. n-denní úhrny Filipova Huť Hřensko 1
51.20
32.50
2
71.85
40.95
3
83.30
44.15
5
103.40
51.95
7
114.25
60.25
Tabulka 2: Střední hodnoty meteorologických stanic Filipova Huť, Děčín (Hřensko)
Průměrná hodnota jednodenního úhrnu srážek ve městě Děčín (Hřensko) je 32,50 mm, zatímco ve Filipově Huti se průměrná hodnota jednodenního úhrnu srážek vyšplhala až na 51,20 mm. Což je výrazný rozdíl. Zatímco rozdíl průměrných hodnot u jednodenního úhrnu srážek byl 18,70 mm, tak u dvoudenního úhrnu je již 30,90 mm. Pro třídenní 13
2D Graph 4.5 Filipova Hut Decin (Hrensko)
4
summer / winter
3.5 3 2.5 2 1.5 1 0.5 1965
1970
1975
1980
1985 1990 years
1995
2000
2005
2010
Obrázek 9: 2D graf podílu srážek v létě a v zimě v letech 1968-2007 pro města Filipova Huť a Děčín (Hřensko)(Zdroj: Vlastní) úhrny srážek je rozdíl středních hodnot 39,15 mm, pro pětidenní úhrny srážek 51.45 a pro sedmidenní úhrny srážek činí rozdíl středních hodnot 54 mm.
n-denní úhrny Filipova Huť Hřensko 1
350.69
151.06
2
608.80
241.31
3
755.30
270.14
5
943.59
354.23
7
1394
564.80
Tabulka 3: Rozptyly meteorologických stanic Filipova Huť, Děčín (Hřensko)
Zajímavé je také, jakých hodnot nabývaly u měst rozptyly (viz tabulka 3). Jelikož rozptyly charakterizují proměnlivost srážek u konkrétního města v letech 1968-2007. Pro jednodenní úhrny srážek nabývá rozdíl rozptylů obou měst hodnotu 199.63. A stejně jako v předchozím pozorování středních hodnot, s přibývajícím n-denním úhrnem srážek 14
narůstají i rozdíly rozptylů mezi oběma městy. Tabulky 4 a 5 srovnávají kvantily n-denních úhrnů srážek mezi městy Filipova Huť a Hřensko. n/p
80,00% 90,00% 91,00% 92,00% 93,00% 95,00% 99,00%
1
74.45
82.90
84.18
84.68
84.92
85.75
117.70
2
93.55
119.80
120.52
122.47
124.83
127.80
133.50
3
112.40
134.60
136.04
136.52
136.68
143.35
152.90
5
131.65
150.25
158.61
161.18
161.82
170.50
189.20
7
139.50
170.25
172.13
177.04
182.96
195.15
252.40
Tabulka 4: Kvantily města Filipova Huť
n/p
80,00% 90,00% 91,00% 92,00% 93,00% 95,00% 99,00%
1
43.35
58.20
58.28
58.39
58.51
59.00
67.10
2
55.05
64.25
66.69
69.43
72.27
75.60
96.40
3
61.55
71.35
72.27
74.72
77.68
80.1
108.40
5
77.00
86.25
88.05
88.86
89.34
92.55
118.30
7
83.10
95.25
98.17
99.89
101.21
106.40
160.00
Tabulka 5: Kvantily města Děčín (Hřensko)
Proměnná T , neboli T -letá úroveň charakterizuje dobu, za kterou v daných letech překročí naměřená data určitou hranici. 1/T je pravděpodobnost překročenní T -leté úrovně (viz tabulka 6). Například pravděpodobnost výskytu stoleté vody je pouhé jedno procento. Pravděpodobnost pětileté vody je již dvacet procent. V tabulce 7 se testují data 162 stanic. Data jsou z let 1961-2007 seřazena od nejnižší po nejvyšší umístěnou stanici. 15
p
0,95 0,90 0,85 0,80 0,70 0,50 0,40 0,30 0,20 0,10
0,05
0,02
0,01
T
1,05 1,11 1,18 1,25 1,43 2,00 2,50 3,33 5,00 10,00 20,00 50,00 100,00 Tabulka 6: pravděpodobnost překročení T -leté úrovně.
Nadmořská výška [m] Průměrná hodnota rozptyl minimum maximum 150-225
35,15
172,17
12,6
95,8
225-271
37,55
240,57
12,7
162,9
271-320
40,18
241,75
16,5
160,4
320-373
38,33
232,34
14,7
159,4
373-426
41,68
347,15
14,4
137,1
426-460
38,02
252,85
13,5
134
460-517
41,04
321,43
16,1
159
517-583
40,38
268,9
15,3
159
583-1118
42,99
440,92
15,2
233,8
Tabulka 7: Pořadové statistiky 162 stanic
Data byla pořízena ze 162 stanic v České republice za léta 1961 až 2007. Se zvyšující se nadmořskou výškou, se zvětšují také průměrné hodnoty a s nimi podle očekávání rostou i rozptyly. Tento fakt popisují v tabulce maxima. V nadmořské výšce 583-1118 m.n.m napadlo výrazněji více srážek než v ostatních lokalitách s nižší nadmořskou výškou. Nejvíce srážek dopadne na města s vyšší nadmořskou výškou. Skutečnost, že ve vyšších místech napadne o poznání více srážek než v oblastech s nižší nadmořskou výškou zachycuje graf na obrázku 10, který byl vytvořen pomocí programu MATLAB R2009b. Data měst na obrázku jsou seřazena vzestupně podle nadmořské výšky. Z grafu lze vyčíst, že u některých měst s vyšší nadmořskou výškou dopadlo v letech, kdy byly na území České republiky zaznamenány povodně, o poznání více srážek než 16
v oblastech s nižší nadmořskou výškou. Graf České republiky (viz obrázek 11) zobrazuje průměrný roční úhrn srážek v letech 1961-1990. Nejvíce množství srážek spadlo na severu Čech v Krkonoších, v místě, kde pramení řeka Labe. Také v povodí řeky Moravy, řeky Bečvy a na Šumavě si můžeme všimnout vyššího výskytu srážek. Nejvíce srážek spadlo v oblastech s vyšší nadmořskou výškou, nejčastěji v povodí velkých řek. Ve Středozemí již takové srážkové úhrny nezaznamenáváme ani v oblastech velkých řek.
Obrázek 10: 3D graf jednodenního úhrnu srážek měst. Česká republika v letech 1961 až 2007 (Zdroj: Vlastní)
Mapa na obrázku 12 popisuje srážkové rozpoložení v roce 2002. Tento rok byl jedním z nejhorších let z pohledu povodňové aktivity v historii Čech. U zrodu povodní stála tlaková níže, která způsobila masivní déšť v delším časovém rozmezí. Situace se později uklidnila. To však netrvalo dlouho a došlo k dalšímu přibývání srážek. Půda nedokázala absorbovat předchozí množství srážek a opět musela čelit dalšímu přívalu deště. To mělo za následek zvýšení hladin řek, nádrží na území Čech. Výsledkem toho bylo, že se nádrže, které se již v předchozím období dešťů musely vyrovnávat s velkým množstvím deště, nyní již nedokázaly čelit tak velkému množství srážek. Došlo k rozlití vody a zaplavení údolí, 17
míst, které byly obývané lidmi. Na snímku je vidět jak silné deště postihly jihozápadní Čechy. Tato oblast byla jednou z nejpostiženějších oblastí v Čechách při záplavách v roce 2002.
Obrázek 11: Průměrný roční úhrn srážek v letech 1961-1990 (zdroj: ČHMÚ)
Obrázek 12: Roční úhrn srážek v roce 2002(zdroj: ČHMÚ)
18
4
Statistická hypotéza
Statistickou hypotézou [2] myslíme určité předpoklady o rozdělení náhodných veličin. Jestliže se tyto předpoklady týkají hodnot parametrů rozdělení náhodné veličiny, pak se nazývají parametrické hypotézy. V opačném případě se jedná o hypotézy neparametrické. Jsou-li hypotézou specifikovány všechny parametry rozdělení sledované veličiny, tzn. rozdělení je určeno jednoznačně, pak je hypotéza jednoduchá. Pokud není některý parametr rozdělení specifikován jednoznačně, např. je vymezen Při testování statistických hypotéz vždy porovnáváme dvě hypotézy. Jedna hypotéza, tzv. nulová (testovaná), je hypotéza, kterou testujeme. Značíme ji obvykle H0 . Druhou hypotézou je tzv. alternativní hypotéza, kterou obvykle značíme H1 (někdy také HA ). Posuzování nulové hypotézy H0 je založeno na následující úvaze: • Předpokládáme, že platí hypotéza H0 . • Rozhodneme se, kterým náhodným pokusem (například založeném na náhodném výběru) hypotézu ověříme. Určíme, která náhodná veličina bude výsledkem pokusu.
• Určíme hladinu spolehlivosti α neboli pravděpodobnost (míru rizika) toho, že hypotézu H0 neoprávněně zamítneme, ačkoliv platí. Parametr α se přitom stanovuje
jako malé, obvykle 0,05 nebo 0,01.
• V oboru možných hodnot použité náhodné veličiny určíme takovou část, do níž
za platnosti H0 padne výsledek veličiny s pravděpodobností α. Tato část oboru
možných hodnot se nazve kritický obor a značíme V
• Pokud nyní hodnota náhodné veličiny padne do kritického oboru, hypotézu zamí-
táme, neboť nastal jev, který by za platnosti H0 měl jen velmi malou pravděpodobnost a jehož výskyt tudíž svědčí proti platnosti nulové hypotézy.
19
Výsledkem testu je rozhodnutí o nulové hypotéze. Přijetí hypotézy H0 znamená, že ji považujeme za možnou. Zamítnutí hypotézy H0 je ekvivalentní přijetí hypotézy HA . Testování hypotéz je tedy proces, při němž se na základě náhodného výběru rozhodne pro testovanou nebo alternativní hypotézu.
4.1
Chyby testu
Uvedený postup může také vést k chybnému zamítnutí testované hypotézy (tzv. chyba I. druhu) nebo k chybnému přijetí testované hypotézy (tzv. chyba II. druhu). Pravděpodobnost chyby I. druhu je označována jako hladina významnosti testu. α = P (T ∈ H0 )
(1)
Pravděpodobnost, že hodnota testovacího kritéria padne do oboru přijetí H0 , jestliže platí HA , tzn. pravděpodobnost chyby II. druhu, je β = P (T ∈ V | HA )
(2)
Doplněk k β se nazývá síla testu (hovoříme také o silofunkci) 1 − β = P (T ∈ W | H1 )
20
(3)
5
Analýza rozptylu
Analýzu rozptylu zkráceně ANOVU vymyslel Sir Ronald Fisher Aylmer. Byl to anglický statistik, evoluční biolog, eugenista a genetik, který žil v letech 1890–1962. Dnes je považován za zakladatele moderní statistické vědy. Analýza rozptylu pracuje se vzorci, které jsou uvedeny v mnoha knihách zabývajících se statistikou. µk = EX k
(4)
k = 1, 2, . . .
Rovnice (4) počítá obecný moment k-tého řádu. Existuje-li moment µ1 a je-li konečný, pak µk = E (E − EX)k
k = 0, 1 . . .
(5)
číslo k se nazývá centrální moment k-tého řádu. Pro k=2 se nazývá centrální moment rozptyl a označuje se symbolem σ 2 Můžeme jej také označit var X, variancí veličiny X. Nechť X je náhodná veličina s konečnou střední hodnotou EX. Potom rozptyl náhodné veličiny X definujeme jako DX = E (E − EX)2
(6)
Rozptyl je střední hodnota kvadrátů odchylek od střední hodnoty. Pokud hodnota na pravé straně existuje, rozptylem je myšleno, jak moc náhodná veličina X kolísá kolem střední hodnoty EX. Nechť σ > 0 Pak symbolem a3 =
µ , σ3
a4 =
µ σ4
rozumíme šikmost, respektive
špičatost. Chceme-li zjistit, zda číselná hodnota y závisí na číselné či slovní hodnotě x použije se tzv. jednofaktorovou analýzu variance - analysis of variance (ANOVA). Analýza rozptylu ve své parametrické podobě předpokládá normalitu rozdělení a tzv. homoskedasticitu (identické rozptyly). Tato metoda se dnes běžně používá k běžným statistickým účelům. Součet čtvercových odchylek n pozorovaných hodnot měřitelné proměnné y od jejich průměru rozkládá na součty čtvercových odchylek, které odrážejí jednotlivé zdroje měnlivosti hodnot proměnné y. Na základě tohoto poznatku se pak testuje, zda jsou tyto zdroje účinné. Při jednofaktorové analýze rozptylu se zkoumá, zda lze změny hodnot měřitelné proměnné y vysvětlovat faktorem x. Symbolem yij označujeme j-té pozorování z i-té skupiny. Kde i = 1, . . . , k a j = 1, . . . , n. yi je průměr i-té skupiny (i-té třídy), symbolem y se označuje celkový průměr (průměr průměrů yi ). Součty čtverců jednotlivých odchylek 21
všech n hodnot proměnné y od jejich průměrů (tj. Součet Sy ) se rozkládá na součty odchylek Sy,m (meziskupinová variabilita) a Sy,ν (vnitroskupinová variabilita). Předpokládá se, že meziskupinová variabilita je vysvětlena daným faktorem, zatímco vnitroskupinová nikoli. Součty čtverců jednotlivých odchylek mají určitý počet stupňů volnosti. Počet stupňů volnosti součtu čtverců m je určen tím, kolik je z těchto veličin nezávislých. Existuje-li mezi m veličinami c lineárních vztahů, má součet čtverců těchto m veličin m − c stupňů
volnosti. Počet stupňů volnosti pro Sy je ν = n − 1 součet čtverců odchylek od společného průměru dostaneme
k X n X i=1 j=1
(yij − y) = 0
(7)
Počet stupňů volnosti Sy,m má µ1 = k − 1 stupňů volnosti a pro součet čtverců odchylek
mezi skupinami dostaneme
k X i=1
(yi − y) ni = 0
(8)
Pro Sy,ν máme ν2 = n − k stupňů volnosti a konečně pro součet čtverců odchylek √ uvnitř skupin (reziduální) (yi − y) ni = 0 Jednofaktorová analýza rozptylu vychází z předpokladu, že každý z k nezávislých vý-
běrů hodnot proměnné y pochází z normálního rozdělení a že tato rozdělení se středními hodnotami µ1 , µ2 ,. . ., µk mají stejné rozptyly rovné konstantě. Normální rozdělení je úplně specifikováno stření hodnotou a rozptylem. Jestliže se předpokládá, že všech k normálních rozdělení proměnné y má stejný rozptyl, je možné účinnost faktoru x posuzovat pouze podle toho, zda na různých úrovních tohoto faktoru mají či nemají tato normální rozdělení stejné střední hodnoty. Hypotézu H0 o neúčinnosti faktoru lze tedy formulovat jako hypotézu, že rozdělení proměnné y mají na různých faktoru x stejné střední hodnoty, tzn. H0 : µ1 = µ2 = . . . = µk . Platí hypotéza H0 , lze dokázat, že střední hodnota výběrových charakteristik Sy,m /ν1 a Sy,ν /ν2 jsou rovny rozptylu σ 2 a že poměr
F =
Sy,m /ν1 Sy,ν /ν2
(9)
má rozdělení F o ν1 = k – 1 a ν2 = n – k stupních volností. Jestliže hypotéza H0 o neúčinnosti (o nezávislosti) neplatí, je střední hodnota charakteristiky Sy,m /ν1 větší než střední hodnota charakteristiky Sy,ν /ν2 . Na to, že H0 neplatí, lze tedy usuzovat 22
z vysokých hodnot poměru F . Tento poměr se volí za testové kritérium a kritický obor je při zvolené hladině významnosti vymezen nerovností (10)
F > F1−α (k − 1, n − k)
na jejíž pravé straně je 100(1 – α)% kvantil F –rozdělení o k – 1 a n – k stupních volnosti. Pokud hodnota testového kritéria padne do tohoto kritického oboru, přijímá se hypotéza o účinnosti faktoru, res. Hypotéza o závislosti proměnné y na proměnné x. [3] Variabilita
Součet čtverců Stupeň volnosti Průměrný součet Testové kritérium
Faktor x
Sy,m
ν1 = k − 1
Sy,m /ν1
Reziduální
Sy,ν
ν2 = n − k
Sy,ν /ν2
Celkový
Sy
ν =n−1
F =
Sy,m /ν1 Sy,ν /ν2
Tabulka 8: ANOVA
5.1
Příklad 1.
Cílem pokusu je srovnání několik středních hodnot vytvořených z úhrnů srážek v jednotlivých letech. V případě, že se střední hodnoty nerovnají, zamítá se hypotéza H0 . Libovolně jsem zvolil Moravská města: 1-České Budějovice, 2-Tábor, 3-Hranice u Nových Hradů, 4-Jistebnice, 5-Kamenice nad Lipou, 6-Kardašova Řečice, 7-Kovářov, 8-Lodhéřov, 9-Nadějkov-Větrov, 10-Nová Bystřice, 11-Ševětín-Švamberk, 12-Soběnov. Pro ověření hypotézy jsem zvolil náhodně data z jednodenního úhrnu srážek za léta 1961-2007. Tvar rozdělení, jeho střední hodnoty a variability zobrazuje na obrázku vrubový krabicový graf 13. Střední čáry v krabicích představují mediány, hranice krabic pak 1. a 3. kvartily. Oblast mezi 1. a 3. kvartilem se označuje jako interkvartilový interval. Rozdíl mezi 1. a 3. kvartilem se označuje jako mezikvartilové rozpětí. Extrémní hodnoty znázorňují koncové úsečky. Body, které se nacházejí ve větší vzdálenosti od mediánu jsou zobrazeny v podobě červených křížků a reprezentují extrémní hodnoty. Na dolních koncích žádné extrémy nejsou. Oproti tomu jich je mnoho na horních koncích. Tento rozdíl 23
značí, že distribuce není symetrická. U hodnoty 10, tedy u Ševětínu, dokonce extrém převyšuje 120 mm hranici.
120
100
80
60
40
20 1
2
3
4
5
6
7
8
9
10
11
12
Obrázek 13: Vrubový krabicový graf (Zdroj: Vlastní)
Výsledky tohoto softwaru odpovídají obrázku 14. Error zde označuje Reziduální (vnitroskupinovou) variabilitu. Slovo Total v programu označuje součet vnitroskupinové a meziskupinové variability. Columns je experimentální (meziskupinová) variabilita. SS jsou součty čtverců, df stupně volnosti, MS Průměrný součet čtverců. F je testovací statistika a poslední sloupeček označuje tabelovanou hodnotu (p-value). Jelikož testovací statistika je větší než tabelovaná hodnota (p-value), můžeme pro pětiprocentní hladinu významnosti a 1-denni úhrn srážek v daných letech 1961-2007 nulovou hypotézu zamítnout a tvrdit, že střední hodnoty jednotlivých měst nejsou stejné. To jinými slovy znamená, že úhrny srážek nejsou závislé na výběru jednotlivých měst. Stejně bychom zamítli nulovou hypotézu pro 2, 3, 5 a 7 denní úhrny srážek v daných letech. Nulová hypotéza byla zamítnuta. V takovém případě je dobré použít jiný test. Vhodným testem je Tukeyův test, který srovnává střední hodnoty naměřených dat a ukazuje na velké rozdíly mezi středními hodnotami. 24
Obrázek 14: Jednofaktorová analýza rozptylu (Zdroj: Vlastní)
25
Towns
Count
Mean
Homogeneous Groups
Tábor
47
63,9128
X
Kamenice nad Lípou
47
70,3149
XX
Kovářov
47
71,2894
XX
Kardašova Řečice
47
73,7426
XX
Jistebnice
47
74,8362
XX
Ševětín Švamberk
47
75,4298
XX
Lodhéřov
47
79,0745
XX
České Budějovice
47
79,7128
XX
Nadějkov Větrov
47
79,7277
XX
Nová Bystřice
47
84,6447
XX
Hranice u Nových Hradů
47
85,6447
X
Soběnov
47
91,3468
X
Tabulka 9: Method: 95,0 percent Tukey HSD
Contrast
Sig. Difference
České Budějovice – Hranice u Nových Hradů
+/- Limits
-5,93191
21,0719
České Budějovice – Jistebnice
4,8766
21,0719
České Budějovice – Kamenice nad Lípou
9,39787
21,0719
České Budějovice – Kardašova Řečice
5,97021
21,0719
České Budějovice – Kovářov
8,4234
21,0719
České Budějovice – Lodhéřov
0,638298
21,0719
-0,0148936
21,0719
České Budějovice – Nadějkov Větrov
26
České Budějovice – Nová Bystřice
-4,93191
21,0719
České Budějovice – Soběnov
-11,634
21,0719
České Budějovice – Ševětín Švamberk
4,28298
21,0719
15,8
21,0719
Hranice u Nových Hradů – Jistebnice
10,8085
21,0719
Hranice u Nových Hradů – Kamenice nad Lípou
15,3298
21,0719
Hranice u Nových Hradů – Kardašova Řečice
11,9021
21,0719
Hranice u Nových Hradů – Kovářov
14,3553
21,0719
Hranice u Nových Hradů – Lodhéřov
6,57021
21,0719
Hranice u Nových Hradů – Nadějkov Větrov
5,91702
21,0719
1,0
21,0719
Hranice u Nových Hradů – Soběnov
-5,70213
21,0719
Hranice u Nových Hradů – Ševětín Švamberk
10,2149
21,0719
21,7319
21,0719
Jistebnice – Kamenice nad Lípou
4,52128
21,0719
Jistebnice – Kardašova Řečice
1,09362
21,0719
Jistebnice – Kovářov
3,54681
21,0719
Jistebnice – Lodhéřov
-4,2383
21,0719
Jistebnice – Nadějkov Větrov
-4,89149
21,0719
Jistebnice – Nová Bystřice
-9,80851
21,0719
Jistebnice – Soběnov
-16,5106
21,0719
Jistebnice – Ševětín Švamberk
-0,593617
21,0719
10,9234
21,0719
České Budějovice – Tábor
Hranice u Nových Hradů – Nová Bystřice
Hranice u Nových Hradů – Tábor
*
Jistebnice – Tábor 27
Kamenice nad Lípou – Kardašova Řečice
-3,42766
21,0719
Kamenice nad Lípou – Kovářov
-0,974468
21,0719
Kamenice nad Lípou – Lodhéřov
-8,75957
21,0719
Kamenice nad Lípou – Nadějkov Větrov
-9,41277
21,0719
Kamenice nad Lípou – Nová Bystřice
-14,3298
21,0719
Kamenice nad Lípou – Soběnov
-21,0319
21,0719
Kamenice nad Lípou – Ševětín Švamberk
-5,11489
21,0719
Kamenice nad Lípou – Tábor
6,40213
21,0719
Kardašova Řečice – Kovářov
2,45319
21,0719
Kardašova Řečice – Lodhéřov
-5,33191
21,0719
Kardašova Řečice – Nadějkov Větrov
-5,98511
21,0719
Kardašova Řečice – Nová Bystřice
-10,9021
21,0719
Kardašova Řečice – Soběnov
-17,6043
21,0719
Kardašova Řečice – Ševětín Švamberk
-1,68723
21,0719
Kardašova Řečice – Tábor
9,82979
21,0719
Kovářov – Lodhéřov
-7,78511
21,0719
Kovářov – Nadějkov Větrov
-8,4383
21,0719
Kovářov – Nová Bystřice
-13,3553
21,0719
Kovářov – Soběnov
-20,0574
21,0719
Kovářov – Ševětín Švamberk
-4,14043
21,0719
7,3766
21,0719
Lodhéřov – Nadějkov Větrov
-0,653191
21,0719
Lodhéřov – Nová Bystřice
-5,57021
21,0719
Kovářov – Tábor
28
Lodhéřov – Soběnov
-12,2723
21,0719
Lodhéřov – Ševětín Švamberk
3,64468
21,0719
Lodhéřov – Tábor
15,1617
21,0719
Nadějkov Větrov – Nová Bystřice
-4,91702
21,0719
Nadějkov Větrov – Soběnov
-11,6191
21,0719
Nadějkov Větrov – Ševětín Švamberk
4,29787
21,0719
Nadějkov Větrov – Tábor
15,8149
21,0719
Nová Bystřice – Soběnov
-6,70213
21,0719
Nová Bystřice – Ševětín Švamberk
9,21489
21,0719
Nová Bystřice – Tábor
20,7319
21,0719
Soběnov – Ševětín Švamberk
15,917
21,0719
27,434
21,0719
11,517
21,0719
Soběnov – Tábor
*
Ševětín Švamberk – Tábor
Tabulka 10: Rozdíly středních hodnot
Tabulka 10 se vztahuje na více porovnávací procedury. Určuje, které střední hodnoty jsou výrazně odlišné od ostatních. Tabulka ukazuje odhadovaný rozdíl mezi každým párem střední hodnoty. Hvězdička byla umístěna vedle dvou párů, které vykazují statisticky významné rozdíly na 95,00% úrovni spolehlivosti. V tabulce 9 jsou 2 homogenní skupiny identifikovány pomocí sloupce X. V každém sloupci, úrovně obsahující formu X, skupinu středních hodnot, v němž nejsou žádné statisticky významné rozdíly. Tato metoda se v současné době používá k rozlišení mezi středními hodnotami v Tukeyově HSD (honestly significant difference) proceduře. Při použití metody je 5,0% riziko, že se jeden nebo více dvojic výrazně liší, i když jejich skutečný rozdíl je roven nule. 29
5.2
Ověření normality
K tomu aby se mohla provést analýza rozptylu je zapotřebí ověřit normalitu rozdělení a tzv. homoskedasticitu (identické rozptyly).
Obrázek 15: Testy generované pomocí programu Statgraphics Centurion XV. (Zdroj: Vlastní)
Na obrázku 15 jsou v pravé části graficky znázorněny data sedmidenních úhrnů srážek stanice České Budějovice. Modrá křivka charakterizuje normální rozdělení a histogram představuje data stanice České Budějovice. Pomocí tohoto grafu vidíme, jak moc se liší data stanice české Budějovice od normálního rozdělení. Pro hodnoty 50-100 Histogram přelézá křivku normálního rozdělení a podle tohoto grafu bychom mohli zamítnout nulovou hypotézu, data se kterými pracujeme nepochází z normálního rozdělení. Stejnou situaci znázorňuje qantilový graf, který je uveden níže. Porovnává hodnoty kvantilů s lineární křivkou, která znázorňuje tvar normálního rozdělení. V tomto případě qantily korelují s lineární křivkou pro hodnoty 0-150, proto bychom mohli usoudit, že data sedmidenních úhrnů srážek stanice české Budějovice pochází z normálního rozdělení. V levé části okénka programu Statgraphics Centurion XV. Se provádí testy. Jedním z nich je Shapiro-Wilkův 30
test. Shapiro-Wilkův test je založen na porovnávání dvou škálových statistik. Hodnota W je ”korelační koeficient”, který udává, jak těsně data korelují s křivkou normálního rozdělení. Hodnota p-value pak udává jaké chyby bychom se dopustili, kdybychom zamítli nulovou hypotézu H0 . Pokud je p-value větší než zvolená hladina alfa, přijímáme nulovou hypotézu H0 . Pro sedmidenní úhrny srážek v Českých Budějovicích v letech 1961 až 2007 dopadl ShapiroWilkův test zamítnutím hypotézy H0 . P-value provedených testů byla menší než hodnota 0,05. Můžeme tedy s 95% pravděpodobností odmítnout hypotézu H0 , že data Českých Budějovic pro sedmidenní úhrny srážek pořízených v letech 1961 až 2007 pochází z normálního rozdělení. Stejně tomu bylo u ostatních klasických testů ověřujících normální rozdělení. Test špičatosti, test šikmosti, test chí-kvadrát dopadly stejně jako ShapiroWilkův test. Kolmogorův-Smirnovův Test tvrdil, že nejmenší p-value mezi provedenými testy je 0,0599798, tedy větší skoro o jednu desetinu než hodnota 0,05. Z tohoto hlediska nebylo možné s 95% pravděpodobností zamítnout hypotézu H0 . Tedy nemůžeme zamítnout hypotézu, že naměřená data Českých Budějovic pochází z normálního rozdělení. Další okénko v programu Statgraphics Centurion XV. znázorňuje koncové plochy, tak aby odpovídaly co možná nejlépe normálnímu rozdělení. Výstupní údaje ukazují, že pravděpodobnost získaná z hodnoty, která je menší nebo rovna 63,7702, je ze 32,3491% vhodná pro normální rozdělení. Ale ani z tohoto grafu zcela přesně nevyčteme, zda data pochází z normálního rozdělení. Další tabulka v programu Statgraphics Centurion XV. Ukazuje, která disribuční funkce nejlépe odpovídá datům. Data nejvíce odpovídají distribuční funkci loglogistic se třemi parametry. Dále pak následuje noncentral t distribuční funkce, která je dvouparametrová. Stejné testy byly provedeny pro sedmidenní úhrny srážek města Jistebnice v letech 1961 až 2007 (viz obrázek 16). Pomocí Shapiro-Wilkůvova testu byla zamítnuta hypotéza H0 . P-value nabývá hodnoty 1, 47344 · 10−8 , tedy je menší než hodnota 0,05. S 95% pravděpodobností odmítnoutáme hypotézu H0 , že data města Jistebnice pro sedmidenní úhrny srážek pořízených v letech 1961 až 2007 pochází z normálního rozdělení. Na rozdíl od dat Českých Budějovic dopadl Kolmogorův-Smirnovův test s 95% pravděpodobností zamítnutím nulové hypotézy. Nejvíce těmto datům jako v předchozím případě odpovídají 31
distribuční funkce loglogistic se třemi parametry a noncentral t distribuční funkce.
Obrázek 16: Testy generované pomocí programu Statgraphics Centurion XV. (Zdroj: Vlastní)
Proto, aby se prováděný test stal přesnějším bylo zapotřebí pracovat s větším množstvím dat. Místo dat českých Budějovic se použila data více měst. A vzhledem k tomu, že u programu Statgraphics Centurion XV. Není zřejmé, jak postupovali vývojáři tohoto programu při sestavování testů. Nikde v tomto programu není k dispozici zdrojový kód, pomocí něhož by bylo možné ověřit, že testy provedené v tomto programu pracují správně. Proto slouží program Statgraphics Centurion XV. pouze k ověření propočtů. Následně byl pomocí programu matlab vytvořen chi-kvadrát test dobré shody a výpočty prokázaly, že data těchto n-denních úhrnů srážek nepocházela z normálního rozdělení. Ukázalo se, že analýza rozptylu pro tato data není vhodná. Pro málo dat je vhodným testem Kolmogorův – Smirnův test. Jeho účelem je jako v předchozím případě analýzy rozptylu ověřit normalitu rozdělení. H0 : Výběr pochází z rozdělení N (µ, σ 2 ) proti alternativě HA : Výběr nepochází z rozdělení N (µ, σ 2 ). Testovým kritériem je náhodná veličina 32
Dn = sup |Fn (x) − F0 (x)| = max (D1∗ , D2∗ , . . . , Dn∗ )
(11)
F0 (x) je distribuční funkce testového rozdělení, Fn (x) je empirická (výběrová) distribuční
funkce
Di = sup |Fn (x) − F0 (x)|
(12)
x
Pozorovaná hodnota je maximální odchylkou mezi teoretickou a empirickou distribuční funkcí.
Di = max (D1∗ , D2∗ , . . . , Dn∗ )
(13)
kde Di = max i − F0 (x) ; 1−i − F0 (x) i = 1, . . . , n n n
(14)
Hypotézu zamítáme na hladině významnosti α, pokud pozorovaná hodnota testové statistiky překročí tabelovanou kritickou hodnotu Dn (α).[4]
5.3
Příklad 2.
Pomocí Kolmogorova-Smirnova testu určete, zda data stanic Filipova Huť, Děčín (Hřensko) patří do normálního rozdělení.
33
Obrázek 17: Část zdrojového kódu programu KLM-2012 (Zdroj: Vlastní)
Na obrázku 17 je část zdrojového kódu programu KLM-2012. Tento program testuje hypotézu, zda data, která jsou k dispozici, pochází z normálního rozdělení. Celý algoritmus byl vytvořen podle testovací statistiky Kolmogorov- Smirnov. Program si nejprve nalezne na disku počítače data, se kterými bude pracovat. V další části testu seřadí vzestupně n-denní úhrny srážek. Následně se spočítá směrodatné odchylka a průměrná hodnota. Dále se pomocí for cyklu vytvoří podprogram, který postupně spočítá teoretickou distribuční funkci. Pomocí hodnot teoretické distribuční funkce se v cyklu dále spočítají pomocné hodnoty pro výpočet maximálních odchylek teoretické a empirické distribuční funkce v bodech nespojitosti empirické distribuční funkce. Na závěr se mimo for cyklus určí pozorovaná hodnota testovacího kritéria. Která se dále bude porovnávat s tabelovanou kritickou hodnotou. Pokud pozorovaná hodnota testovací statistiky překročí tabelovanou kritickou hodnotu, zamítneme na hladině významnosti nulovou hypotézu, podle níž data pochází z normálního rozdělení. Test byl proveden ve dvou částech pro jednodenní úhrny srážek v letech 1968-2007. První část byla vytvořena pomocí vstupního příkazu, který program obsahuje a umí s tímto vstupem pracovat. Druhá část pak pomocí vlastního kódovacího jazyka, který 34
byl vytvořen v okénku Function M-File a část zdrojového kódu je uvedeno na obrázku 17. V obou případech dopadl test stejným výsledkem, zamítnutím hypotézy H0 a přijetím hypotézy HA . Červená křivka (viz obrázek 19) znázorňuje Standardní distribuční funkci a modrá křivka znázorňuje Empirickou distribuční funkci, která byla získána z dat jednodenního úhrnu srážek města Děčín (Hřensko). Testovací statistika k je maximální rozdíl mezi těmito křivkami. Na obrázku 20 byla empirická distribuční funkce získána z dat jednodenního úhrnu srážek města Filipova Huť. Ze dvou předešlých obrázku je patrné, že testovací statistika k bude výší u grafu charakterizující Kolmogorův-Smirnův test pro data Filipova huť. Jinými slovy, maximální vzdálenost mezi křivkami na tomto grafu je větší než vzdálenost na grafu, který byl vytvořen z dat města Děčín (Hřensko). Z toho dále můžeme usoudit, že opravdu ve městě Filipova Huť, které je u nás nejvýše položeným městem, napadne v daném časovém rozmezí více srážek než v nejníže položení městě Děčín (Hřensko).
35
Obrázek 18: Okno programu KLM-2012 (Zdroj: Vlastní)
Obrázek 19: Graf jednosměrného Kolmogorova-Smirnova testu pro město Děčín (Hřensko). (Zdroj: Vlastní) 36
Obrázek 20: Graf jednosměrného Kolmogorova-Smirnova testu pro město Filipova Huť. (Zdroj: Vlastní)
6
Hodnoty extrémních distribucí
V předchozí kapitole se pracovalo s daty X1 , . . . , Xn . s tzv. náhodným výběrem, který tvoří distribuční funkci (df ), označím ji velkým písmenem F . Malé f je pak derivací funkce F a nazývá se hustotou funkce F . Nechť Mn je maximální hodnota z náhodného výběru. Mn = max{X1, . . . , Xn} pro n ≥ 2. Pro minimum pak platí M in = min{X1 , . . . , Xn } =
−max{X1, . . . , Xn .}. Pravděpodobnost maxima P (Mn ≤ x) = P (X1 ≤ x, ..., Xn ≤ x) =
F n (x). Podle Fisher–Tippetovy věty lze vlastnosti těchto maxim, chvostů popsat pomocí jednoho parametru. Tomuto parametru se jinak říká extremální hodnota γ. Předpokládejme, že můžeme najít posloupnost reálných čísel an > 0 a bn tak, že posloupnost (Mn−bn ) /an konverguje v distribuci, t.j. F n (an x + bn ) → G (x) , n → ∞,
(15)
pro nějakou nedegenerovanou d.f. G (x) Jestliže podmínka platí, říkáme, že F je v ”domain of attraction of G” (sféře přitažlivosti), F ∈ M DA(G) (viz [5]) Jestliže F ∈ M DA(G) potom G je typu jedné z následujících tří distribučních funkcí.
37
Fréchet Φ1/γ (x) =
0, x ≤ 0
γ>0
exp −x−/γ , x > 0 Weibull Ψ1/γ (x) =
o n exp − (−x)1/γ , x ≤ 0
γ>0
1 x>0
Gumbel Λ (x) = exp (−e−x ) , x ∈ R Provedením parametrizace dostaneme [6]
1. Gumbel (EV0) g0 (x) = G0 (x) e−x , pro všechna x; 2. Fréchet (EV1), α > 0 : g1,α(x) = αG1,α (x) x−(1+α) , x ≥ 0, 3. Weibull (EV2), α < 0 : g2,α(x) = |α| G2,α (x) (−x)−(1+α) , x ≤ 0, Rozdělení s tzv. lehkými chvosty zahrnuje tyto dvě velké třídy rozdělení, tj. Gumbellovu třídu - velmi zjednodušeně řečeno hustota je exponenciální funkce (tj. patří sem např. normální a exponenciální rozděleni) a Weibullovu třídu, rozdělení na konečném nosiči (např. rovnoměrné rozdělení). Fréchetovu třídu označujeme jako rozdělení s těžkými chvosty. G0 (x) = exp e−x ∀x ∧ Gγ (x) = exp (1 + γx)−1/x , 1 + γx > 0, γ 6= 0.
(16)
Rovnici (16) označujeme jako zobecněné rozdělení extrémních hodnot GEV (Generalized extreme value distribution). Podle dobře známé formule (1 + γx)1/γ → exp (x) pro γ → 0
(17)
V takovém případě dostáváme rovnost Gγ (x) = G0 (x) ovšem jen pro nulový parametr γ. 38
• Gγ = G0 Gumbelova standardní distribuční funkce, • Gγ je Fréchetova distribuční funkce pro γ > 0, • Gγ je a Weibullova distribuční funkce pro γ < 0, K určení extremální hodnoty slouží mnoho odhadů. Nejznámějším a nejpoužívanějším z nich je Hillův odhad [7], který se používá v případě Fréchetovy třídy. Odhad byl vytvořen v programu Matlab R2009b, za předpokladu, že data pochází z rozdělení s těžkými chvosty (Fréchetovy třídy).
Obrázek 21: Část zdrojového kódu programu Hill-2012. (Zdroj: Vlastní)
Program Hill-2012 provádí Hillův odhad pro data meteorologických stanic. Část zdrojového kódu (viz obrázek 21) zobrazuje dva do sebe vnořené for cykly. Díky vnitřnímu a vnějšímu cyklu program provádí Hillův odhad pro více než jednu meteorologickou stanici. Graf na obrázku 22 popisuje závislost Hillova odhadu na počtu kroků k pro data stanic České Budějovice (AMS), Tábor (AKS), Hranice u Nových Hradů (MSS), Jistebnice (MSS), Kamenice nad Lipou (ASS), Kardašova Řečice (MSS), Kovářov (MSS), Lodhéřov (MSS), Nadějkov-Větrov (AKS), Nová Bystřice (ASS), Ševětín-Švamberk (MSS), Soběnov (MSS).
39
2D Hill Plot 0.045 CB HNH J KnL KR K L NV NB SS S T
0.04 0.035
estimation
0.03 0.025 0.02 0.015 0.01 0.005 0
0
10
20 30 number of steps
40
50
Obrázek 22: 2D Hill plot pro sedmidenní úhrny srážek a 46 kroků (Zdroj: Vlastní)
2D Hill Plot 0.12 CB HNH J KnL KR K L NV NB SS S T
0.1
estimation
0.08
0.06
0.04
0.02
0
1
2
3
4
5 6 number of steps
7
8
9
10
Obrázek 23: 2D Hill plot pro sedmidenní úhrny srážek a 10 kroků (Zdroj: Vlastní)
40
Pro Hillův odhad s 10 kroky (viz obrázek 23) je u některých měst možné pozorovat ustálení mezi 6 a 0 procentem. Těchto procent nabývá graf pro 0 až 10 kroků. Zřejmě se jedná o interval, na kterém by mohla ležet hledaná hodnota parametru γ.
2D Hill Plot 0.2 CB HNH J KnL KR K L NV NB SS S T
0.18 0.16
estimation
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0
1
1.5
2
2.5 3 3.5 number of steps
4
4.5
5
Obrázek 24: 2D Hill plot pro sedmidenní úhrny srážek a 5 kroků (Zdroj: Vlastní)
Pro Hillův odhad s 5 kroky (viz obrázek 24) se ustálení pohybuje v rozmezí 10 a 2 procent. Rozptyl se s ubývajícím počtem kroků zvyšuje. Proto by již nemělo smysl snižovat počet kroků. Naopak pro vyšší počet k by se nám zvyšovalo vychýlení. S malým k roste rozptyl a klesá vychýlení a naopak s velkým počtem kroků k roste vychýlení a klesá rozptyl. Je třeba zvolit optimální podíl dat, tak aby se například výrazně nezvětšil rozptyl na úkor vychýlení, pak se jen velmi těžce může odhadnou parametr. Zajímavější může být porovnán míst s různými nadmořskými výškami. Na obrázku 25 je znázorněn Hillův odhad jednodenních úhrnů srážek měst Filipova Huť (modře) a Hřensko (zeleně). Bylo provedeno deset kroků a z obrázku lze vyčíst, že hodnota γ parametru bude nejspíš mezi druhým až šestým krokem, tedy mezi 4,3% a 2,3%. Pro k rovno pěti, se γ parametr pohybuje mezi 2,5% a 1%. Ale pro příliš malé k se musí počítat s velkým rozptylem. 41
U následujících testů [8] se porovnávají dvě hypotézy. Pokud je první tzv. nulová hypotéza H0 splněna, pak data pochází z rozdělení s lehkými chvosty, v opačném případě přijímáme alternativní hypotézu HA , data pochází z rozdělení s těžkými chvosty.
2D Hill Plot 0.06 Filipova Hu Dìèín (Høensko) 0.05
estimation
0.04
0.03
0.02
0.01
0
1
2
3
4
5 6 number of steps
7
8
9
10
Obrázek 25: 2D Hill plot jednodenního úhrnu srážek stanic Filipova Huť a Děčín (Hřensko) (Zdroj: Vlastní)
Dalším statistickým postupem je W k test. Pokud je splněna rovnice (20) připouštíme nulovou hypotézu, že data na předem zvolené hladině významnosti pochází z rozdělení s lehkými chvosty.
2 k k X k − Xn−k+1:n 1X Xn−i+1:n Wk = 2 , X k = Pk k (k − 1) i=1 Xn−i+1:n − X k i=1
(18)
Statistika W k má asymptoticky normální rozdělení se střední hodnotou µk a rozptylem σk2
µk =
4 (k − 2) 1 1 , σk2 = . (k − 1) (k − 1)2 (k + 1) (k + 2) 42
(19)
2D Hill Plot 0.08 Filipova Hu Dìèín (Høensko)
0.07 0.06
estimation
0.05 0.04 0.03 0.02 0.01 0
1
1.5
2
2.5 3 3.5 number of steps
4
4.5
5
Obrázek 26: 2D Hill plot jednodenního úhrnu srážek měst Filipova Huť a Děčín (Hřensko) (Zdroj: Vlastní) Kritický obor pro oboustrannou alternativu je
|Wk∗ | > u1−α/2 ,
(20)
kde u1−α/2 je (1 − α/2) kvantil normálního rozdělení a
Wk∗ = (Wk − µk ) /σk ,
(21)
Jestli je splněna nerovnost, přijímáme alternativu, že dané rozdělení má lehké chvosty oproti alternativě rozdělení má těžké chvosty. Pro k = 20 nám W k statistika na 5% hladině významnosti konverguje ke Gumbelovu rozdělení. Rozdělení má tedy lehké chvosty. Jak však zvolit správné k, tak abychom se přiblížili správnému výsledku? Existuje spousta návrhů, jak zvolit k. k/n = 0.2 pro 50 ≤ n ≤ 500 a k/n = 0.1 pro 500 < n ≤ 5000. Někde se podle Ga-lambose doporučuje √ zvolit k = 2 n. Jinou alternativou, jak určit, zda má dané rozdělení těžké či lehké chvosty, je volba jiného statistického testu. Dalším testem je test od vývojářů Segers a Teugels.
43
2D Wk plot 1.4 Wk 1.2
1
Wk
0.8
0.6
0.4
0.2
0
0
5
10 number of steps
15
20
Obrázek 27: 2D Graf jednodenního úhrnu srážek pro stanici Praha-Klementinum v letech 1961-2005 (Zdroj: Vlastní)
Sm =
5 m
m X i=1
T (ξi )
!2
, T (x) := 1 −
6x (1 + x)2
(22)
Sm > χ21 (1 − α), kde χ21 (ε) označuje ε-kvantil χ2 rozdělení s jedním stupněm volnosti.
Pro m → ∞ testovací statistika konverguje k χ2 rozdělení
Na obrázku 28 je graf charakterizující Sm statistiku. Test byl proveden pro data města
Praha-Klementinum. Data zaznamenávají jednodenní úhrny srážek v letech 1961-2004. Modrá čára zde znázorňuje Sm statistiku. Oproti předchozímu testu je patrné, že má testovací statistika Sm velký rozptyl. Což by mohlo vést k domněnce, že tato data pochází z rozdělí s těžkými chvosty. Ovšem testovací statistika Sm byla menší než P-value o jednom stupni volnosti a pětiprocentní hladině významnosti. Zde se pracovalo s pouhými 15 skupinami, což je z pohledu testovací statistiky velmi malé číslo. Zapotřebí by bylo rozdělit data do více skupin, tak aby počet prvků v daných skupinách byl vyšší než 2. To ale s tak malým množstvím dat není možné.
44
2D plot Sm 0.25 Sm
0.2
Sm
0.15
0.1
0.05
0
1
1.5
2
2.5
3 m
3.5
4
4.5
5
Obrázek 28: 2D Graf jednodenního úhrnu srážek pro stanici Praha-Klementinum v letech 1961-2005 (Zdroj: Vlastní) Aby bylo testování ještě přesnější, přistupme k poslednímu testu T k, který vytvořili C. Neves, J. Picek a M. I. Fraga Alves. V případě, že se prokáže nulová hypotéza, test konverguje ke Gumbelovu rozdělení, v našem případě tedy k rozdělení s lehkými chvosty. Nulová hypotéza je zamítnuta na hladině významnosti alfa pokud platí (23)
∗ Tk,n < g(α/2) ,
kde ∗ Tk,n =
1 k
Pk
Xn:n − Xn−k:n
i=1 (Xn−i+1:n − Xn−k:n )
− log k.
(24)
Testovací statistika T k má oproti statistice W k klesající charakter. Ovšem stejně jako u předchozích testů přijímá nulovou hypotézu, že data města Praha - Klementinum pro jednodenní úhrny srážek v letech 1961-2007 pochází z rozdělení s lehkými chvosty, V tomto případě přesněji z Gumbelova rozdělení. Ve všech třech testech nám vyšlo, že data Města Praha - Klementinum pochází z rozdělení s lehkými chvosty. Pro tato data by nebylo vhodné provést Hillův odhad. 45
5
3
2D plot Tk
x 10
Tk 2.5
Tk
2
1.5
1
0.5
0
0
5
10
15
20 25 number of steps
30
35
40
45
Obrázek 29: 2D Graf jednodenního úhrnu srážek pro stanici Praha-Klementinum v letech 1961-2005 (Zdroj: Vlastní)
7
Závěr
Cílem práce bylo porovnat naměřená data na jednotlivých stanicích v České republice. Pomocí naměřených dat provést analýzu rozptylu, která neprokázala závislost mezi středními hodnotami vybraných měst. Posléze se neprokázala ani normalita, která je zapotřebí při provedení analýzy rozptylu. Z tohoto důvodu jsem provedl Tukeyův test, který porovnává střední hodnoty, které jsou výrazně odlišné od ostatních. středních hodnot. Následně provedl odhady, které mají určit k jakému rozdělení se naměřená data blíží. Při zjištění rozdělení by se dalo lépe charakterizovat chování úhrnů srážek v některých srážkami ohrožených oblastech. A nadále by se tento fakt mohl uplatnit v regresní analýze, bylo by o další velmi důležitý poznatek více. Pomocí programu Statgraphics Centurion se podařilo dokázat, že data měst České Budějovice a Jistebnice se blíží ke tří parametrické distribuční funkci loglogistic. Pozadu nezůstalo ani Gumbelovo rozdělení. Navíc se prokázalo, že data města Praha-Klementinum pocházejí z rozdělení s lehkými chvosty. To je důležitý fakt pro vznik dalších propočtů. Dále by bylo možné pokračovat v regresní analýze. Bylo by však zapotřebí použít i jiných faktorů ovlivňujících srážkové rozpoložení na našem území. Faktory by byly v podobě naměřených dat, například proudění větru, teplota vzduchu, přesná časová data a jiné. 46
Z vypracovaných grafů porovnávaných stanic, vyplývá přehled nejvíce ohrožených oblastí České republiky, jsou to stanice s vyšší nadmořskou výškou. Například Filipova Huť, Churáňov, Lysá Hora, Praděd a další horské oblasti především Šumava. Doporučoval bych v těchto oblastech věnovat větší pozornost k vybudování protipovodňových opatření. Například vybudovat atypické protipovodňové zátarasy. Omezit kácení lesů. Včas vypouštět nádrže přehrad a častěji čistit koryta řek.
47
Reference [1] Český hydrometeorologický ústav: Meteorologické stanice ČHMÚ. ČESKÝ HYDROMETEOROLOGICKÝ ÚSTAV. Český hydrometeorologický ústav: Meteorologické stanice ČHMÚ [online]. neznámé. neznámé, 32.01.2012 [cit. 2012-04-09]. Dostupné z: http://portal.chmi.cz/files/portal/docs/poboc/OS/stanice/ShowStations_CZ.html.
[2] Wikipedia: Cs.wikipedia.org, Testování statistických hypotéz [online]. neznámé, 201202-10 [cit. 2012-02-22]. Dostupné z: http://cs.wikipedia.org/wiki/Testova ´nı ´_statisticky ´ch_hypote ´z.
[3] HINDLS Richard, KAŇOKOVÁ Jara, NOVÁK Ilja, Metody statistické analýzy pro ekonomy. Praha 2004: Management Press. ISBN 80-86419-26-6. [4] MATEMATIKA PRO INŽENÝRY 21. STOLETÍ: Kolmogorovův-Smirnovův test řešený příklad. MATEMATIKA PRO INŽENÝRY 21. STOLETÍ [online]. VŠB - TU Ostrava, 1.9.2009 [cit. 2012-03-15]. Dostupné z: http://mi21.vsb.cz/flash-animace/kolmogorovuv-smirnovuv-test-reseny-priklad.
[5] EMBRECHTS Paul, KLÜPPELBERG Claudia, MIKOSCH Thomas. Modelling Extremal Events for Insurance and Finance. Berlín: Springer, 1997. ISBN 3-540-60931-8. [6] REISS Rolf-Dieter; THOMAS Michael. Statistical Analysis of Extreme Values : with Applications to insurance, Finance, Hydrology and Other Fields. Berlin : [s.n.], 2001. Modeling by Extreme Value Distributions, p. 15. ISBN 3-7643-6487-4. [7] HILL Bruce M. (1975). A simple general approach to inference about the tail of a distribution. Ann. Statist. 3, 1163 – 1174. [8] NEVES Cláudia., PICEK Jan., FRAGA Alves M. I. (2005). The contribution of the maximum to the sum of excesses for testing max-domains of attractions. J. Statist. Planning Infer.,
48
