Ústav matematiky a statistiky Přírodovědecká fakulta Masarykova univerzita Statistická inference I Zadáni domácího úkolu - rok 2016 Stanislav Katina, Veronika Bendová katina@math.muni.cz, 375612@math.muni.cz 12. prosince 2016 Katina, S., Bendová, V., 2016: Statistická inference I 1 Instrukce k domácímu úkolu: Odevzdává se jeden pdf soubor nazvaný prijmeni-jmeno-text-statinf-l-2016.pdf (obsahuje řešení příkladů, obrázky, Cií-kód napsaný v TjíjXu), jeden zdrojový soubor naprogramovaných funkcí prijmeni-jmeno-source-statinf-l-2016.R a jeden soubor ^-kódu konkrétních zadání z DÚ prijmeni-jmeno-priklady-statinf-l-2016. R, který používá tento zdrojový kód. Dejte si záležet na přehlednosti programovaného kódu, na doplnění komentářů a vhodné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ž číselných nebo grafických. I to bude součástí celkového hodnocení úkolu. Na psaní 'Sít-kódu doporučuji Tj5]Xovský 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}. DÚ je nutné odevzdat 7 dní před termínem zkoušky, na který se přihlásíte. (12. prosince 2016) Katina, S., Bendová, V., 2016: Statistická inference I 2 Příklad 1. Jana s Bárou a Vojtíškem dostali hořko-mléčný adventní kalendář, ve kterém je polovina čokolád hořkých a polovina čokolád mléčných, přičemž příchutě čokolád jsou v kalendáři rozmístěny náhodně. O čokolády se děti rozhodly podělit rovným dílem, ale protože je Vojtíšek nejmenší, dovolily mu sestry, aby svůj díl čokolád snědl jako první. 1. Vypočítejte, jaká je pravděpodobnost, že Vojtíšek, který vůbec nemá rád hořkou čokoládu, bude mít ve svém dílu: a) všechny čokolády mléčné, b) maximálně jednu čokoládu hořkou, c) více než polovinu přidělených čokolád mléčných. Získané pravděpodobnosti interpretujte, vždy uveďte odpověď. 2. Nakreslete graf pravděpodobnostní a distribuční funkce rozdělení, které popisuje rozložení počtu mléčných čokolád ve Vojtíškově přídělu. Grafy řádně okomentujte. Příklad 2. Načtěte datový soubor 01-one-sample-mean-skull-mf.txt obsahující údaje o délce a šířce lebky ze starověké egyptské populace. 1. Vytvořte tabulku základních charakteristik (průměr, rozptyl, směrodatná odchylka, medián, l.kvartil, 3.kvar-til, IQR, koeficient šikmosti a špičatosti) pro délku a šířku lebky starověké egyptské populace. 2. Nakreslete histogramy (v rel. škále) pro délku a šířku lebky starověké egyptské populace. 3. Vypočítejte Pearsonův korelační koeficient pro vztah délky a šířky lebky egyptské populace a závislost obou veličin demonstrujte pomocí tečkového grafu. 4. Vypočítejte dvourozměrný jádrový odhad hustoty délky a šířky lebky starověké egyptské populace, zakreslete jej pomocí funkce image() a superponujte jej a) konturovými křivkami dvourozměrného jádrového odhadu, b) teoretickými konturami dvourozměrného normálního rozdělení, kde střední hodnoty /ii a /i2, rozptyly a\ a a\ a korelační koeficient p nahraďte jejich MLE odhady získanými z dat. Poznámka: Pohlídejte si, aby kontury přesně seděly s barevnými přechody image grafu. 5. Dvourozměrný jádrový odhad hustoty vykreslete také pomocí funkce persp(). Hustotu rozsekejte na 12 intervalů, kde hodnoty v těchto intervalech budou odpovídat barvám terrain.colors(12). 6. Bylo v tomto případě vhodné použít Pearsonův korelační koeficient k určení závislosti mezi délkou a šiřkou lebky? Kterou informaci z datového souboru jsme při počítání korelačního koeficientu zanedbali a jak se to projevilo v grafech dvourozměrného jádrového odhadu? 7. Zohledněte informaci, kterou jsme na začátku příkladu zanedbali, a vygenerujte nový tečkový graf, oddělující barevně zanedbanou vlastnost, a superponujte jej konturami dvourozměrného jádrového odhadu. Určete nové korelační koeficienty zohledňující zanedbanou informaci z datového souboru. Všechny výsledky a grafy okomentujte, své závěry a postupy zdůvodněte. Nezapomeňte na přesnou interpretaci Pearsonových korelačních koeficientů. (12. prosince 2016) Katina, S., Bendová, V., 2016: Statistická inference I 3 Příklad 3. 1. Nakreslete škálovaný logaritmus profilové funkce věrohodnosti normálního rozdělení pro jj,. Na ose x bude jj, a na ose y hiCp(p\x.) = lp(p\x.) — max(Zp(/i|x)). Porovnejte lnCp(p\x.) s kvadratickou aproximací vypočítanou pomocí Taylorova rozvoje ln£p(/i|x) = ln(^|^j) w — ^X(JÍ)(fi — Jí)'2. 2. Nechť skóre funkce S(p) = ln Lp(fi\x). Vezmeme-li derivaci kvadratické aproximace uvedené výše, dostaneme S(fi) ~ — X(Jl)(fi — Jí) nebo —X~1^2(Jľ)S(fí) ~ X1/2(/i)(/i — Jí). Potom zobrazením pravé strany na ose x a levé strany na ose y dostaneme asymptoticky lineární funkci s jednotkovým sklonem. Je postačující mít rozsah osy x rovný (—2,2), protože funkce je asymptoticky (lokálně) lineární na tomto intervalu. Rozumně škálujte osu y. Zobrazte pro (a) n = 10, (b) n = 100 a (c) n = 1000. Použijte (1) X ~ N(0,1) a (2) X ~ (1 — p)N(0,1) +pN(0, 2), kde p = 0.05. Okomentujte rozdíly mezi (a), (b) a (c), stejně jako rozdíly mezi (1) a (2). Příklad 4. maximálně věrohodný odhad /i a a2 Vygenerujte pseudonáhodná čísla z X ~ iV(4,1), n = 1000. 1. Napište logaritmus odhadnuté a profilové funkce věrohodnosti pro /i a a2 a porovnejte maximálně věrohodné odhady parametrů jj, a a2 získané maximalizací těchto funkcí. Nakreslete grafy Ze(/i|x), Zp(/i|x), Ze(cr2|x) a Zp(cr2|x), kde zvýrazníte polohu maxim těchto funkcí. 2. MLE odhady parametrů jj, a a2 z bodu 1 najděte pomocí (a) funkce optimize(), (b) Newton-Raphsonovy metody (naprogramujte), (c) metody sečen (naprogramujte). Získané odhady zaokrouhlete na šest desetinných míst, uspořádejte do přehledné tabulky a vzájemně porovnejte. optimizeQ Newton-Raphson metoda.secen Jip 3% Porovnejte také numerické metody, s jejichž pomocí byly parametry odhadnuty. Která z použitých numerických metod je pro odhadování parametrů efektivnější a proč? Porovnání získaných odhadů a numerických metod slovně popište, rozeberte a své závěry zdůvodněte. 3. Napište logaritmus funkce věrohodnosti pro 9 = (/i, a2)T a prověřte, zda je maximálně věrohodný odhad 0 dostatečně blízko k jeho skutečné hodnotě. Nakreslete graf l(0\x), kde na ru-ové ose bude parametr jj, a na y-ové ose parametr a2, použitím funkce image() a superponujte ho konturovým grafem použitím funkce contour(). Zvýrazněte polohu maxima. (12. prosince 2016) Katina, S., Bendová, V., 2016: Statistická inference I 4 Příklad 5. Simulační studie Na základě simulační studie ověřte, že pokud a) X ~ N(n, cr2), kde fi = 0, a2 = 1, b) X ~ [(1 -p)N(fj,,a%) +pN(fi,a'i)}, kde ^ = 0, of = \,o\ = 4, p = 0.1, potom 1. testovací statistika x — fj, /— tw = -Vn s má Studentovo rozdělení o n — 1 stupních volnosti TV ~ tn-i, 2. testovací statistika má Fisherovo rozdělení o 1 a n — 1 stupních volnosti T^r ~ FiiTJ_i. Použijte rozsahy náhodných výběrů n = 15 a n = 500. Pro každou simulaci X vypočítejte t\y,m (resp. í^m), m = 1,2,..., M, kde M = 1000. Superponujte histogram vygenerovaných testovacích statistik v relativní škále s teoretickou křivkou hustoty Studentova (resp. Fisherova rozdělení). Vygenerované histogramy pro n = 15 a n = 500 jakožto i pro rozdělení (a) a (b) vzájemně srovnejte a srovnání slovně okomentujte. Příklad 6. Poissonovo rozdělení: část 1 — odvození Nechť náhodná veličina X pochází z Poissonova rozdělení, X ~ -Po(A) a realizace X = x. Pro náhodou veličinu X odvoďte 1. pravděpodobnostní funkci f(x,\), 2. věrohodnostní funkci (+ uveďte, jak vypadá jádro) L(\\x), 3. logaritmus věrohodnostní funkce (+ uveďte, jak vypadá jádro) l(\\x), 4. skóre funkci S(X), 5. maximálně věrohodný odhad parametru A, ozn. A, 6. pozorovanou Fisherovu míru informace X(A), 7. rozptyl odhadu parametru A, ozn. Var [A]. (12. prosince 2016) Katina, S., Bendová, V., 2016: Statistická inference I 5 Příklad 7. Poissonovo rozdělení: část 2 — maximálně věrohodné odhady Mějme početnosti úrazů mezi dělníky v továrně, kde početnosti dělníků mn při daném počtu úrazů n jsou uvedeny v následující tabulce (Greenwood a Yule (1920)). n II 0 I 1 I 2 I 3 I 4 I >5 mn || 447 | 132 | 42 | 21 | 3 | 2 1. Vypočítejte maximálně věrohodný odhad parametru A (a) pomocí maximalizace věrohodnostní funkce; výsledek zobrazte do grafu s křivkou věrohodnostní funkce Poissonova rozdělení, (b) pomocí maximalizace logaritmu věrohodnostní funkce; výsledek zobrazte do grafu s křivkou logaritmu věrohodnostní funkce Poissonova rozdělení. 2. Naprogramujte metodu sečen a s její pomocí vypočítejte maximálně věrohodný odhad parametru A. Výsledek zobrazte do grafu s (a) křivkou věrohodnostní funkce Poissonova rozdělení, (b) křivkou logaritmu věrohodnostní funkce Poissonova rozdělení. 3. Vypočítejte odhad rozptylu odhadu parametru A získaného pomocí (a) maximalizace věrohodnostní funkce, (b) maximalizace logaritmu věrohodnostní funkce, (c) metody sečen. 4. Všechny tři odhady parametru A a k nim příslušné odhady rozptylů zaokrouhlete na šest desetinných míst a uspořádejte do přehledné tabulky maximalizace L(A|x) maximalizace Z(A|x) metoda sečen A Var [A] Získané odhady porovnejte s explicitně vyjádřenými odhady (zaokrouhlenými též na šest desetinných míst). Porovnání slovně opište, stejně jako vygenerované grafy. Příklad 8. Poissonovo rozdělení: část 3 — kvadratická aproximace logaritmu funkce věrohodnosti 1. Pro data z příkladu 7 nakreslete škálovaný logaritmus funkce věrohodnosti Poissonova rozdělení. Na a;-ové ose bude A a na y-ové ose ln£(A|x) = Z(A|x) — max(Z(A|x)). Porovnejte ln£(A|x) s kvadratickou aproximací vypočítanou pomocí Taylorova rozvoje ln£(A|x) = ln ^^|x^ w —^I(A)(A — A)2. 2. Nechť skóre funkce S(X) = ln L(A|x). Vezmeme-li derivaci kvadratické aproximace uvedené výše, dostaneme 5(A) = -X(A)(A - A) anebo -X"1/2(A)5(A) «XX/2(A)(A-A). (1) Potom zobrazením pravé strany na a;-ové ose a levé strany na y-ové ose dostaneme asymptoticky lineární funkci s jednotkovým sklonem. Nakreslete graf, kde na a;-ové ose bude vynesena pravá strana rovnosti 1 a na y-ové ose levá strana rovnosti 1. Křivku superponujte lineární křivkou x = y. Rozumně škálujte oj-vou a y-vou osu. Vygenerované grafy řádně okomentujte. (12. prosince 2016)