STP097 STATISTIKA – CVIČENÍ 12.12.2007 EMPIRICKÁ DISTRIBUČNÍ FUNKCE, JEDNOVÝBĚROVÉ TESTY
Postupujte podle zadání. Vše potřebné k dnešnímu cvičení natáhnete z webu do R příkazy: adr="http://artax.karlin.mff.cuni.cz/~kraud8am/stp097/stp097_cvic_2007-12-12.RData" load(url(adr))
Vyzkoušejte, že se vše povedlo: příkaz ls() musí mezi vypsanými objekty ukázat obrazek.edf (předdefinovaná funkce) a doly, rychlost (datové vektory).
1. Empirická distribuční funkce Empirickou distribuční funkci počítá funkce ecdf(x). Její argument x je vektor představující náhodný výběr. Jejím výsledkem je objekt, který můžeme pod nějakým jménem uschovat a dále zpracovávat. Spočtěte a nakreslete empirickou distribuční funkci náhodného výběru z normovaného normálního rozdělení o rozsahu 25 následujícím způsobem: (1) Vygenerujte a uschovejte náhodný výběr z normovaného normálního rozdělení (x = rnorm(25)). (2) Spočtěte empirickou distribuční funkci tohoto výběru (F = ecdf(x)). (3) Objekt F se chová jako funkce, tj. můžeme spočítat jeho hodnoty v libovolném bodě nebo bodech. Zkuste spočítat F(-0.5). (4) Udělejte obrázek: stačí napsat plot(F). Nyní zopakujte tento postup pro beta rozdělení s parametry α = β = 0.5 a rozsah výběru n = 35 [rbeta(n,alpha,beta)]. Podívejme se, jak se při vzrůstajícím počtu pozorování přibližuje empirická distribuční funkce skutečné distribuční funkci. Nakreslíme si empirickou distribuční funkci pro čtyři výběry z N(0, 1) o rozsahu 10, 50, 500, 2000. Každým obrázkem proložíme skutečnou distribuční funkci a dáme si je na jeden list. Taktéž spočítáme maximální absolutní rozdíl mezi skutečnou a empirickou distribuční funkcí. Jádrem výpočtu je připravená funkce obrazek.edf, kterou jste si natáhli z webové adresy. Vypište si, jak vypadá [print(obrazek.edf)]. Jejím jediným argumentem je rozsah výběru n. Funkce vygeneruje data z N(0, 1) a vyrobí obrázek empirické d.f a skutečné d.f., přitom vrací maximální rozdíl mezi empirickou a skutečnou distribuční funkcí. 1
Obrázky zakreslíme na jeden graf tak, že před voláním funkce obrazek.edf napíšeme příkaz par(mfrow=c(2,2)) (kreslící oblast se rozdělí na 2 × 2 obrázky). Pak čtyřikrát zavoláme obrazek.edf s argumentem 10, 50, 500 a 2000 a uvedeme grafiku do původního stavu [par(mfrow=c(1,1))].
2. Jednovýběrové testy: hladina, síla, p-hodnota Uvažujme náhodný výběr X1 , . . . , Xn z rozdělení N(0.5, 2) o rozsahu n = 60. Vygenerujte jeden takový výběr příkazy n = 60 x = rnorm(n,mean=0.5,sd=sqrt(2))
Nyní provedeme jednovýběrový Kolmogorovův-Smirnovův test hypotézy H0 : Xi ∼ N(0.5, 2) proti alternativě, že Xi mají libovolné jiné rozdělení. V R se takový KS test na výběru x provede příkazem ks.test(x,y="pnorm",mean=0.5,sd=sqrt(2))
(pnorm znamená distribuční funkci normálního rozdělení). Prozkoumejte výstup z této funkce: kde je testová statistika, kde je p-hodnota? Rozhodněte, zdali došlo k zamítnutí nulové hypotézy. Nakreslíme si obrázek empirické distribuční funkce spolu s distribuční funkcí za hypotézy: od = min(x)*0.9 do = max(x)*1.1 plot(ecdf(x),xlim=c(od,do)) body = seq(od,do,length=500) lines(body,pnorm(body,mean=0.5,sd=sqrt(2)),col="blue")
Dá se z obrázku okem odhadnout hodnota testové statistiky KS testu? Nyní na ten samý výběr x proveďte postupně test hypotézy H0 : Xi ∼ N(µ, 2) pro µ = 0.6, 0.65, 0.7, 0.75, . . . . Jak se mění výsledek testu? Pokračujte v oddalování střední hodnoty hypotetického rozdělení od střední hodnoty skutečného rozdělení dat po stejných krůčcích, dokud nedojde k zamítnutí hypotézy. Pak si nakreslete obrázek empirické distribuční funkce dat a distribuční funkce za platnosti hypotézy. Zopakujte to samé zadání s t-testem: nejdříve proveďte t-test hypotézy H0 : E Xi = 0.5 proti alternativě H1 : E Xi 6= 0.5 na původním datovém souboru x příkazem t.test(x,mu=0.5) a prozkoumejte výstup z funkce t.test. Pak zkoumejte výsledky t-testu hypotéz H0 : E Xi = µ pro µ = 0.6, 0.65, 0.7, 0.75, . . . , dokud nedojde k zamítnutí H0 . Došlo k němu dříve nebo později než u KS testu? 2
3. Jednovýběrové testy – simulace Nyní budeme simulovat hladinu a sílu jednovýběrových testů. Vyrobíme si jednoduchou funkci, která provede test na data x a vrátí pouze p-hodnotu. vem.ph = function(x,test,...) { test(x,...)$p.value }
Vyzkoušejte si ji na původní data s KS testem a t-testem: vem.ph(x,ks.test,y="pnorm",mean=0.5,sd=sqrt(2)) vem.ph(x,t.test,mu=0.5)
Teď vygenerujeme M = 1000 výběrů o rozsahu n = 60 z rozdělení N(0.5, 2) a uspořádáme do matice n × M : n = 60 M = 1000 x.vyb = matrix(rnorm(n*M,0.5,sqrt(2)),nrow=n,ncol=M)
Na každý výběr (každý sloupec) provedeme KS test hypotézy H0 : Xi ∼ N(0.5, 2) a získáme jeho p-hodnotu: ks.ph = apply(x.vyb,2,vem.ph,ks.test,y="pnorm",mean=0.5,sd=sqrt(2))
Aplikovali jsme funkci vem.ph s testem ks.test na sloupce matice x.vyb a získali vektor p-hodnot pro těchto 1000 výběrů. Nakreslete si jejich histogram. Výsledné p-hodnoty jsou náhodné veličiny, jaké je v tomto případě (tj. když platí hypotéza) jejich rozdělení? Váš úsudek si teoreticky zdůvodněte a ověřte provedením KS testu na výběr 1000 p-hodnot (tzn. na vektor ks.ph). Spočtěte, jaký podíl p-hodnot je menších než 0.05, tj. mean(ks.ph<0.05). Co odhaduje toto číslo? Zopakujte tuto úlohu za následujících podmínek: • Generujte výběry z N(µ, 2) pro µ = 0.7 a µ = 0.9, testujte stále hypotézu H0 : Xi ∼ N(0.5, 2) KS testem. Jak se mění rozdělení p-hodnot? Jak se mění počet p-hodnot menších než 0.05 a co to znamená? • Generujte výběry z N(0.5, 2) a provádějte t-test hypotézy H0 : E Xi = 0.5. Interpretujte výsledky. • Generujte výběry z N(µ, 2) pro µ = 0.7 a µ = 0.9, testujte hypotézu H0 : E Xi = 0.5 t-testem. Interpretujte výsledky. 3
4. KS test: data o důlních nehodách Proměnná doly, kterou jste si na začátku načetli, obsahuje okamžiky významných důlních neštěstí ve Velké Británii mezi lety 1875 a 1951. Zajímá nás, zda během sledovaného období docházelo k důlním nehodám rovnoměrně v průběhu času. Proměnná obsahuje přepočítaná data. Původní data1 udávala intervaly (ve dnech) mezi jednotlivými nehodami. V proměnné doly jsou skutečné okamžiky událostí (tj. kumulativní součty intervalů mezi událostmi) vydělené celkovou dobou pozorování, která byla 26263 dní. Pro představu si data vypište. Jsou-li události náhodně rozloženy v čase, měla by tato data představovat (uspořádaný) náhodný výběr z rovnoměrného rozdělení. Vykreslete si histogram a empirickou distribuční funkci. Z histogramu se zdá, že události jsou rozloženy nerovnoměrně, že dříve byly nehody častější. Je to jen optický dojem, nebo je to opravdu významné? Otestujte rovnoměrnost použitím testu Kolmogorov–Smirnov [ks.test(doly,"punif")] a rozmyslete si význam výsledků.
5. Jednovýběrový t-test, znaménkový test: měření rychlosti světla Kolem roku 1880 provedl A. A. Michelson pokusy za účelem stanovení rychlosti světla. Dnes známe “skutečnou hodnotu” rychlosti světla, takže můžeme posoudit úspěšnost tehdejších měření. Skutečná rychlost světla ve vakuu je 299 792.5 km/s. V prostředí, kde Michelson měřil, je ”správná” rychlost 299 710.5 km/s. K disposici máme výsledky 23 pokusů v proměnné rychlost. Od všech hodnot je odečteno 299 000 km/s. Popisné statistiky a grafy. Pomocí příkazu summary(rychlost) získáme základní informace o datech. (Další popisné statistiky lze získat pomocí funkcí quantile, median, mean, var, sd a podobně.) Nakreslete si histogram a empirickou distribuční funkci. Představu o rozdělení dat dává rovněž krabicový diagram (boxplot), který získáme příkazem boxplot(rychlost). Prostřední vodorovná čára je ve výši mediánu, horní a dolní hranice prostřední ”krabice” ve výši kvartilů (přibližně). Krajní vodorovné čáry ukazují největší/nejmenší pozorování ležící do vzdálenosti 1.5 krát výška krabice od nejbližší hranice krabice (tzn. od nejbližšího kvartilu, přibližně). Odlehlá pozorování jsou znázorněna samostatnými body. t-test. Testujme hypotézu, že Michelson měřil správně, proti alternativě, že měřil špatně. Jinými slovy, nulová hypotéza je, že naměřená data odpovídají skutečné střední hodnotě 710.5 km/s, alternativa je, že nikoli. Použijeme tedy jednovýběrový t-test, který porovnává průměr pozorování s hypotetickou střední hodnotou: t.test(rychlost,mu=710.5) 1Tabulka 1 v článku Maguire, Pearson & Wynn (1952). 4
Rozmyslete si, co znamenají jednotlivé vypsané údaje. Jak souvisí testování hypotézy s intervalovým odhadem, který je ve výpisu rovněž uveden? Znaménkový test. Připomeňme si, o co jde. Testujme hypotézu, že medián rozdělení, z něhož pocházejí data, je roven 710.5, proti alternativě, že tomu tak není. Pokud je 710.5 skutečně medianem, měl by počet hodnot nad 710.5 být blízký 23/2 (23 je počet pozorování length(rychlost)). Neobvykle mnoho (blízko 23) nebo neobvykle málo (blízko 0) hodnot nad 710.5 bude svědčit proti tomu, že median je 710.5. V našem případě je počet hodnot nad 710.5 (tedy testová statistika) roven 17: stat = sum(rychlost>710.5) stat Testová statistika má za platnosti hypotézy binomické rozdělení s parametry 23 a 0.5. Znaménkový test provedeme pomocí funkce binom.test: binom.test(stat,23,.5) Prostudujte si výpis a rozmyslete si závěr.
5
