Ústav matematiky a statistiky Přírodovědecká fakulta Masarykova univerzita Statistická inference II Zadáni domácího úkolu - rok 2017 1. část Stanislav Katina, Veronika Bendová katinaSmath.muni.cz, xbendovavSmath.muni.cz 25. května 2017 Katina, S., Bendová, V., 2017: Statistická inference II 1 Instrukce k domácímu úkolu: Odevzdává se jeden pdf soubor nazvaný prijmeni-jmeno-text-statinf-ll-2017.pdf (obsahuje řešení příkladů, obrázky, ^-kód napsaný v TjíjXu), jeden zdrojový soubor naprogramovaných funkcí prijmeni-jmeno-source-statinf-11-2017. R a jeden soubor ^ť-kódu konkrétních zadání z D Ú prijmeni-jmeno-priklady-statinf-11-2017. R, který používá tento zdrojový kód. Dejte si záležet na přehlednosti programovaného kódu, na doplnění komentárů avhodného užití zavedených pravidel, které máte k dispozici v prezentaci Standards of programming in R: R style guide. Také věnujte svou pozornost a čas dostatečným popisům vašich úvah a zvolených postupů a interpretacím výsledků, ať už slovních nebo grafických. I to bude součástí celkového hodnocení úkolu. Na psaní ^ť-kódu doporučuji TgXovský balíček listings a vytvoření prostředí v hlavičce dokumentu pomocí následujícího kódu: \lstset{1anguage=R, % nastavenie jazyka R basicstyle=\footnotesize\ttfamily, % typ pisma R-kodu comment st y le =\ttf amily \ color {f arbal }■, % farba komentára k funkciám numbersty1 e = \color{farba2}\f ootnotesize , % farba a velkost cislovania numbers=left, % cislovanie vlavo stepnumber=1, % cislovanie po krokoch jedna frame=leftline, % vytvorenie lavej hraničnej čiary breaklines=true} % zalomenie riadkov V textu potom kód vkládáme do prostředí \begin{lstlisting} a \end{lstlisting}. Kompletní řešení domácího úkolu je nutné nahrát do odevzdávárny v IS nejpozději 7 dní před termínem zkoušky, na který se přihlásíte. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 2 Příklad 1. Vylepšená věrohodnost pomocí g(0): 1. Nakreslete logaritmus relativní funkce věrohodnosti parametru p binomického rozdělení Bin(7V, p), kde N = 10 a n = 8, superponovaný jeho kvadratickou aproximací. 2. Nakreslete logaritmus relativní funkce věrohodnosti g(p) = logit(p) = ln j^- (při stejném zadání N a n jako v (1)), superponovaný jeho kvadratickou aproximací. 3. Nakreslete graf porovnávající vzájemně logaritmus relativní funkce věrohodnosti s její kvadratickou aproximací získanou na základě parametru p (ad 1) a aproximací získanou za základě parametrické funkce g(p) (ad 2). 4. Vypočítejte Waldův a věrohodnostní 100 x (1 — a)% empirický DIS pro p. 5. Vypočítejte Waldův a věrohodnostní 100 x (1 — a)% empirický DIS pro g(p) z bodu (2) a transformujte jej zpět do originální škály. 6. Vzájemně porovnejte Waldovy empirické DIS pro p a pro g(p) po zpětné transaformaci do originální škály a věrohodnostní empirické DIS pro p a pro g(p) po zpětné transaformaci do originální škály. Který z intervalů vykazuje lepší vlastnosti a proč? 7. Naprogramujte dvě numerické metody: metodu bisekce (funkce bisekce()) a metodu sečen (funkce metoda.secen()) ke zpřesnění hranic věrohodnostních intervalů spolehlivosti. Požadovaná forma výstupu příkladu: • dvě samostatně použitelné funkce bisekce() a metoda.secen() s naimplementovanými iteračními metodami; • trojice grafů: (i) graf s parametrem p na ose x a log. rel. věroh. funkcí + její kvadratickou aproximací na ose y; (ii) graf s parám, funkcí g(p) na ose x a příslušnou log. rel. věroh. funkcí + její kvádr, aproximací na ose y; (iii) graf s parametrem p na ose x a log. rel. věrohodnostní funkcí + její kvádr, aproximací pomocí parametru p a pomocí parám, funkce g(p) na ose y (na základě tohoto grafu porovnejte kvalitu obou kvádr, aproximací); • tabulka hranic intervalů spolehlivosti: parametr/p.funkce Waldův IS - dh Waldův IS - hh Věroh. IS - dh Věroh. IS - hh P g{p) g(p) zpětně transf. do škály p • tabulka přesnějších hranic věrohodnostních intervalů spolehlivosti získaných pomocí vlastnoručně naprogramovaných funkcí bisekceQ a metoda.secenQ: parametr/p.funkce pův. - dh pův. - hh m.bisekce - dh m.bisekce - hh m.sečen - dh m.sečen - hh P g{p) g(p) z.t. do šk.p Poznámka: Hodnoty ve výsledných tabulkách i textu zaokrouhlete na šest desetinných míst. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 3 Příklad 2. Test o směrodatné odchylce a: Z archivních materiálů máme k dispozici původní kraniometrické údaje o délce lebky mužů a žen ze starověké egyptské populace (soubor one-sample-mean-skull-mf.csv). Současně máme k dispozici průměrné hodnoty délky lebky a hodnoty směrodatných odchylek pro muže a ženy novověké egyptské populace (délka lebky mužů xm = 177.568 mm se směrodatnou odchylkou sm = 7.526 mm; délka lebky žen x f = 171.962 mm se směrodatnou odchylkou Sf = 7.052 mm; rozsah datového souboru nm = 88, rif = 52). Načtěte datový soubor one-sample-mean-skull-mf.csv, kde proměnná skuli.L označuje délku lebky (v mm) starověké egyptské populace a proměnná sex označuje pohlaví měřeného jedince. Zaměřte se na délku lebky žen, o které předpokládáme že má normální rozdělení N(fi, a2). 1. Otestujte nulovou hypotézu, že směrodatná odchylka délky lebky žen u starověké egyptské populace je rovna směrodatné odchylce délky lebky žen u novověké egyptské populace. 2. Vypočítejte 100 x (1 — a)% empirický DIS, tj. (au ; an) pro směrodatnou odchylku délky lebky žen starověké egyptské populace, kde koeficient spolehlivosti 1 — a = 0.95. Jak v části (1), tak i v části (2) použijte k otestování nulové hypotézy a ke stanovení příslušných DIS: (a) Waldovu testovací statistiku Uw', (b) skóre testovací statistiku Us; (c) věrohodnostní testovací statistiku Ulr- Požadovaná forma výstupu příkladu: • H0; • íři; • tabulka výsledků: Statistika a test. stat. &h p-hodnota Uw Us ULR • komentář k výsledkům uvedeným v tabulce + zdůvodněné rozhodnutí o tom, kterou testovací statistiku byste v praxi pro konečnou analýzu využili. Poznámka: Hodnoty ve výsledné tabulce i textu zaokrouhlete na čtyři desetinná místa. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 4 Příklad 3. Pokračování příkladu 2: Vraťme se nyní k datovému souboru one-sample-mean-skull-mf.csv, kde proměnná skuli.L označuje délku lebky (v mm) starověké egyptské populace a proměnná sex označuje pohlaví měřeného jedince. 1. Na hladině významnosti a = 0.05 otestujte nulovou hypotézu, že směrodatná odchylka délky lebky žen u starověké egyptské populace je větší než směrodatná odchylka délky lebky žen u novověké egyptské populace. Testování proveďte pomocí: (a) kritického oboru; (b) intervalu spolehlivosti; (c) p-hodnoty. Požadovaná forma výstupu příkladu: • H0; • íři; • tabulka výsledků: Název statistiky a statistika krit.obor au cth p-hodnota • komentář k výsledkům uvedeným v tabulce + zdůvodněné rozhodnutí o nulové hypotéze. Poznámka: Hodnoty ve výsledné tabulce i textu zaokrouhlete na čtyři desetinná místa. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 5 Příklad 4. Simultánní oblasti spolehlivosti + elipsa spolehlivosti pro střední hodnotu a směrodatnou odchylku (dokončení ze cvičení): 1. Nakreslete simultánní množinu spolehlivosti pro 6 = (/i, a)T použitím asymptotického intervalu spolehlivosti pro jj, a exaktního intervalu spolehlivosti pro a. 2. Nakreslete simultánní množinu spolehlivosti pro 6 = (/i, a)T použitím asymptotických intervalů spolehlivosti pro jj, a pro a. 3. Do obrázku dokreslete 100(1 — a) % elipsu spolehlivosti pro 6 = (/í, tjg (pravostranná); c) Hqs : a2 > ctq oproti H13 : a2 < crg (levostranná); kde cr2, = 2. Vypočítejte minimální rozsah náhodného výběru pro testy hypotéz (a)-(c) při a = 0.05 a 1 — /3 = 0.8, pro a2 € {1.01,1.03,1.05,... 2.99, 3.01} (ad (a)), a2 e {2.01, 2.03, 2.05,... 2.99, 3.01} (ad (b)), a2 e {1.01,1.03,1.05,... 1.97,1.99} (ad (c)). Závislost minimálního rozsahu náhodného výběru na hodnotě a2 zakreslete do grafu pomocí křivky (na osu x vyneste parametr a2, na osu y minimální rozsah náhodného výběru.). V grafech barevně odlište minimální rozsah náhodného výběru pro: (i) a2 = 1.4; (ii) cr2 = 1.85; (iii) cr2 = 2.2; (iv) cr2 = 2.8; je-li to možné. Poznámka: Čím více se budeme s hodnotou a2 blížit k hodnotě crg, tím větší minimální rozsah souboru budeme potřebovat. V souladu s touto informací a s vhodným vykreslením grafu rozumně stanovte maximální rozsah souboru (např. TV = 3000, příp. TV = 5000). K našim potřebám nám stačí vědět, že pro a2 blízká hodnotě crg potřebujeme více než 3000 pozorování. Požadovaná forma výstupu příkladu: • funkce min.rozsah(), která pro stanovenou spolehlivost 1 — a a sílu j3* vypočítá pro libovolnou alternativu minimální rozsah náhodného výběru; • tři grafy závislosti TV na cr2; • tabulka výsledků: Alt. hypotéza cr2 = 1.4 cr2 = 1.85 cr2 = 2.2 cr2 = 2.8 oboustranná pravostranná levostranná • komentář ke grafům a k výsledkům uvedeným v tabulce. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 7 Příklad 6. Rozdělení testovací statistiky, síla a silofunkce pro test o rozdílu středních hodnot /ii — /j.2 1. Nechť náhodný výběr X pochází z normálního rozdělení, X ~ N(jj,i,a'f), kde /ii = 4 a a'f = 2.52, a nechť náhodný výběr Y pochází z normálního rozdělení, Y ~ N(/j,2, cr'2), kde \xi = 2 a a\ = 2.52. Rozsahy náhodných výběrů ri\ = 122 = 20. Pomocí simulační studie v ® porovnejte rozdělení testovací statistiky pro test nulové hypotézy Hq: /j-i — /j-2 = Mo oproti alternativní hypotéza -H11 : /ii — /i2 7^ Meh kde jj,q = 0. a. Nasimulujte M = 10 000 pseudonáhodných výběrů laľ takových, že X ~ 7V(/i!,2.52) aľ~ 7V(/i2, 2.52). Pro každé m = 1,..., 10 000 vypočítejte realizaci testovací statistiky íj^ pro nulovou hypotézu Hq: \i2 = 0. Vykreslete histogram testovacích statistik a superponujte jej jednak křivkou hustoty necentrálního í-rozdělení a jednak křivkou hustoty centrálního í-rozdělení. b. Vypočítejte empirickou sílu testu za platnosti alternativní hypotézy Hn : fii — /j.2 = 2. c. Vytvořte animaci zobrazující odchýlení rozdělení testovací statistiky od centrálního í-rozdělení. Zvolte /ii = 4 a/i26 {1,1.5,..., 6, 6.5, 7}. Použijte i. klasický dvouvýběrový í-test; ii. Welchův dvouvýběrový í-test. 2. Nechť nyní X pochází ze směsi dvou normálních rozdělení, t.j. X ~ [pJV(4, 2.52) + (1 — p)N(A, 4.52)], kde p = 0.9 a nechť Y ~ iV(2, 2.52). Proveďte simulační studii popsanou v bodě (1.) pro tento náhodný výběr. Požadovaná forma výstupu příkladu: • funkce cn.rozdeleni(), která pro zadané hodnoty yui, [12, o"ii, cr22, o"2i> a22i ni> n2, M, a a p vykreslí dvojici grafů histogramů (jeden pro klasický í-test a druhý pro Welchův í-test); • animace (ad 1.) zachycující posun necentrálního rozdělní k centrálnímu rozdělení /ii — /j.2 pro pevně zvolené /ii = 4 a měnící se /j.2 G {1,1.5,..., 6.5, 7}. V animaci budou vedle sebe dva grafy histogramů (jeden pro klasický í-test a druhý pro Welchův í-test) superponované křivkou centrálního a necentrálního í-rozdělení; pod popiskem osy x bude uvedena hodnota empirické síly příslušného testu (í-test resp. Welchův test); • animace (ad 2.) zachycující posun necentrálního rozdělní k centrálnímu rozdělení /ii — /j.2 pro pevně zvolené /ii = 4 a měnící se /j.2 G {1,1.5,..., 6.5, 7}. V animaci budou vedle sebe dva grafy histogramů (jeden pro klasický í-test a druhý pro Welchův í-test) superponované křivkou centrálního a necentrálního í-rozdělení; pod popiskem osy x bude uvedena hodnota empirické síly příslušného testu (í-test resp. Welchův test); • komentář ke grafům + popis okometrického srovnání dvojic histogramů; • odpovědi na následující otázky — Co mají společného tvary histogramů klasického a Welchova í-testu a proč? — Jak se liší síla klasického a Welchova í-testu a proč? — Jak se mění síla testů se snižující se vzdáleností fii — /i2? Jak si tento jev vysvětlujete? (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 8 Příklad 7. Pravděpodobnost pokrytí Waldova 95% empirického DIS pro rozdíl středních hodnot \x\—\xi' 1. Nechť A. X3 ~ N(iXj,a2), kde j = 1, 2, >ui = 20, li2 = 25 a a2 = 100; B. X, ~ N(iXj,a^), kde j = 1, 2, >ui = 20, ii2 = 25, cr2 = 100 a cr2 = 150; C. Xj ~ pN(fXj,a2) + (1 - p)N(fXj,al), kde j = 1, 2, ^1 = 20, ^2 = 25, cr2 = 100 a^ = 400; D. X, ~ pN(1x3,(7?) + (1 -p)N(fxj,a?a), kde j = 1,2, ^1 = 20, ii2 = 25, cr2 = 100, a\ = 150, afa = 400 a ^ = 450. Pomocí simulační studie (M = 5 000) vypočítejte pravděpodobnost pokrytí 95% oboustranného Waldova empirického intervalu spolehlivosti pro ii\ — ii2 jako podíl 'Y^?=\ ^-(\^w,m\ < tdf(l — &/2))/M, kde tw,m, m = 1, • • •, 5 000 jsou a. testovací statistiky klasického dvouvýběrového í-testu; b. testovací statistiky Welchova dvouvýběrového í-testu. Rozsahy náhodných výběrů zvolte i. ni = n2 = 5; ii. ni = n2 = 50; iii. ni = n2 = 100 2. Celý postup uvedený v bodě 1. zopakujte také pro /xi = /x2 = 20. Ostatní hodnoty ponechte stejné jako hodnoty zadané v bodě 1. Požadovaná forma výstupu příkladu: • funkce pokryti() která pro zadané hodnoty /xi, /x2, a\, a2, a\a, a2a, M, ni, n2, a a p vypočítá empirickou pravděpodobnost pokrytí DIS pro rozdíl středních hodnot /xi — /x2; • dvě tabulky: — tabulka pravděpodobností pokrytí pro náhodné výběry A.-D., testovací statistiky a.-b. a rozsahy náhodných výběrů i.-iii. pro /xi = 20 a /x2 = 25 (ad 1.); — tabulka pravděpodobností pokrytí pro náhodné výběry A.-D., testovací statistiky a.-b. a rozsahy náhodných výběrů i.-iii. pro /xi = /x2 = 20 (ad 2.); n=5 n=50 n=100 A. klasicky í-test A. Welch í-test B. klasicky í-test B. Welch í-test C. klasicky í-test C. Welch í-test D. klasicky í-test D. Welch í-test • komentář k výsledkům uvedeným v tabulkách + porovnání výsledků získaných z obou tabulek; • odpověď na otázky: — S jakou pravděpodobností pokrývají Waldovy empirické DIS hodnotu 0? — Jak a proč se hodnota pravděpodobnosti pokrytí mění (ad 1.)? Poznámka: Hodnoty ve výsledné tabulce i textu zaokrouhlete na čtyři desetinná místa. (25. května 2017) Katina, S., Bendová, V., 2017: Statistická inference II 9 Příklad 8. Nechť početnosti úmrtí X jako následek kopnutí koněm v Pruských armádních jednotkách (Bortkiewicz, 1898) mají Poissonovo rozdělení s parametrem A, tj. X ~ Poiss(X). Pravděpodobnost, že někdo bude smrtelně zraněný v daném dni, je extrémně malá. Mějme 10 vojenských jednotek za 20-letou periodu (rozsah M = 10 x 20 = 200), kde, při početnostech úmrtí n = 1, 2, 3, 4, 5+ v dané jednotce a v daném roce, zaznamenáváme také početnosti vojenských jednotek mn při daném n, kde M = ~^mn = 200 (viz tabulka). n || 0 | 1 | 2 | 3 | 4 | 5+ mn || 109 | 65 | 22 | 3 | 1 | 0 1. Vypočítejte a. Waldův 95% empirický DIS pro A; b. skóre 95% empirický DIS pro A; c. věrohodnostní 95% empirický DIS pro A. 2. Na hladině významnosti a = 0.05 otestujte nulovou hypotézu Hq : A = 0.6 oproti H\ : A ^ 0.6. Testování proveďte Waldovým, skóre i věrohodnostním přístupem, a to pomocí i. kritického oboru; ii. intervalu spolehlivosti; iii. p-hodnoty. Požadovaná forma výstupu příkladu: • H0 : A = A0, kde A0 = 0.61; • Hi : A ^ A0, kde A0 = 0.61; • tabulka výsledků: A, test.stat. W \d p-hodnota Zw Us ULR • komentář k výsledkům uvedeným v tabulce + zdůvodněné rozhodnutí o nulové hypotéze; • Zamyslete se, jakým způsobem by se dalo zjistit, který typ testu (Waldův skóre, věrohodnostní) je nejvhodnější k otestování hypotézy o parametru A Poissonova rozdělení. Z jakých hledisek můžeme testy vzájemně porovnávat? Vyjmenujte alespoň dvě hlediska, která vás napadají. Poznámka: Hodnoty ve výsledné tabulce i textu zaokrouhlete na čtyři desetinná místa. (25. května 2017)