Masarykova univerzita v Brně Přírodovědecká fakulta sbírka príkladu k predmetu aplikovaná statistika i Marie Budíková, Veronika Bendová Brno, 2016 1 Bodové a intervalové rozdělení četností Přehled použitých funkcí: read.delim, source, head, names, factor, data.frame, sum, cumsum, row.names, va-riacnLrada, barplot, abline, plot, lines, points, axis, table, cbind, rbind, prop.table, dim, round, range, min, max, hist, paste, text. Bodové rozdělení četností Bodové rozdělení četností procvičíme pomocí datového souboru znamky.txt, který obsahuje údaje o známkách z matematiky, angličtiny a pohlaví 20 studentů 1. ročníku. Příklad 1.1. Načtěte soubor znamky.txt. Znakům X, Y, Z vytvořte návěští (X - známka z matematiky, Y - známka z angličtiny, Z - pohlaví studenta). Popište, co znamenají jednotlivé varianty (u znaků X a Y: 1 - výborně, 2 -chvalitebně, 3 - dobře, 4 - neprospěl, u znaku Z: 0 - žena, 1 - muž). setwd("C:/Disk D/ND-Skola/02-Vyuka/04-Aplikovana statistika 2016/Sbirka") source('AS-funkce.R') data <- read.table('znamky.txt', sep='\t', dec='.') head, (data) ## VI V2 V3 ## 1 2 2 0 ## 2 1 3 1 ## 3 4 3 1 ## 4 1 1 0 ## 5 1 2 1 ## 6 4 4 1 names(data) <- c('matematika', 'angličtina', 'pohlavi') head(data) ## matematika angličtina pohlavi ## 1 2 2 0 ##2 1 3 1 ##3 4 3 1 ##4 1 10 ##5 1 2 1 ##6 4 4 1 fl <- factor(data$matematika, levels=c(l,2,3,4), labels=c('výborne','chvalitebné','dobre','nedostatečné')) f2 <- factor(data$anglictina, levels=c(l,2,3,4), labels=c('výborne','chvalitebné','dobre','nedostatečné')) f3 <- factor(data$pohlavi, levels=c(0,1), labels=c('zena','muz')) data2 <- data.frame(fl, f2, f3) names(data2) = c('matematika','angličtina','pohlavi') head(data2) ## matematika ## 1 chvalitebné ## 2 výborne ## 3 nedostatečné ## 4 výborne ## 5 výborne ## 6 nedostatečné angličtina pohlavi chvalitebné zena dobre muz dobre muz výborne zena chvalitebné muz nedostatečné muz 1 Příklad 1.2. Vytvořte a) variační řadu známek z matematiky a známek z angličtiny; b) sloupkový diagram absolutních četností znaků X=Matematika a Y=Angličtina; c) polygon absolutních četností znaků X=Matematika a Y=Angličtina, a) Variační řada známek z matematiky matematika <- data2$matematika ni <- sum(matematika=='výborne') n2 <- sum(matematika== n3 <- sum(matematika== n4 <- sum(matematika== chvalitebné') dobre') nedostatečné') nj <- c(nl,n2,n3,n4) n <- sum(nj) PJ <_ nj/n Nj <- cumsum(nj) Fj <- cumsum(pj) variacni.rada <- data.frame(nj=nj, Nj=Nj, pj=pj, Fj=Fj) row.names(variacni.rada) <- c('výborne', 'chvalitebné', 'dobre', 'nedostatečné') variacni.rada ## výborne ## chvalitebné ## dobre ## nedostatečné nJ NJ PJ FJ 7 7 0.35 0.35 3 10 0.15 0.50 2 12 0.10 0.60 8 20 0.40 1.00 Variační řadu můžeme také získat použitím funkce variacni_rada(X, názvy), která je naprogramovaná ve skriptu AS-funkce. (VR.Mat <- variacni_rada(X=matematika, nazvy=c('výborne', 'chvalitebné', 'dobre', 'nedostatečné'))) ## výborne ## chvalitebné ## dobre ## nedostatečné nJ NJ PJ FJ 7 7 0.35 0.35 3 10 0.15 0.50 2 12 0.10 0.60 8 20 0.40 1.00 Variační řada známek z angličtiny angličtina <- data2$anglictina (VR.Ang <- variacni_rada(X=anglictina, nazvy=c('výborne', 'chvalitebné', 'dobre', 'nedostatečné'))) ## výborne ## chvalitebné ## dobre ## nedostatečné nj Nj pj Fj 4 4 0.20 0.20 4 8 0.20 0.40 7 15 0.35 0.75 5 20 0.25 1.00 2 b) Sloupkový diagram absolutních četností znaků X=Matematika a Y=Angličtina názvy.známek <- c('výborne', 'chvalitebné', 'dobre', 'nedostatečné') # Matematika barplot(VR.Mat$nj, col='white', border='white', axes=T, xlab='Známka', ylab='Pocet pozorováni', names=nazvy.známek, main='Sloupkový diagram pro predmet matematika') abline(h=0:9, col='grey80', lty=2) barplot(VR.Mat$nj, col='blue', axes=F, density=20, border='darkblue', add=T, names=F) #Anglictina barplot(VR.Ang$nj, col='white', border='white', axes=T, xlab='Známka', ylab='Pocet pozorováni', names=nazvy.známek, main='Sloupkový diagram pro předmět angličtina') abline(h=0:9, col='grey80', lty=2) barplot(VR.Ang$nj, col='blue', axes=F, density=20, border='darkblue', add=T, names=F) Sloupkový diagram pro predmet matematika Sloupkový diagram pro predmet angličtina o výborne chvalitebné dobre nedostatečné výborne chvalitebné dobre nedostatečné Známka Známka c) Polygon četností # Matematika plot(l:4, VR.Mat$nj, type='n', xlim=c(0.5,4.5), ylim=c(l,9), xlab='Známka',ylab='Absolútni četnost', main='Polygon četnosti pro predmet matematika', axes=F) abline(h=0:9, col='grey80', lty=2) abline(v=0:9, col='grey80', lty=2) lines(l:4, VR.Mat$nj, col='darkblue', lwd=2 ) points(l:4, VR.Mat$nj,col='darkblue',pch=20, cex=1.2) axis(l, at=0:5, lab=c(11,názvy.známek,'')) axis(2, at=0:10) # Angličtina 3 plot(1:4, VR.Ang$nj, type='n', xlim=c(0.5,4.5), ylim=c(3.5,7.5), xlab='Známka',ylab='Absolutní četnost', main='Polygon četnosti pro predmet angličtina', axes=F) abline(h=seq(3.5, 7.5, by=0.5), col='grey80', lty=2) abline(v=0:9, col='grey80', lty=2) lines(l:4, VR.Ang$nj, col='darkblue', lwd=2 ) points(l:4, VR.Ang$nj,col='darkblue',pch=20, cex=1.2) axis(l, at=0:5, lab=c(11,názvy.známek,'')) axis(2, at=0:10) Polygon četnosti pro predmet matematika Polygon četnosti pro predmet angličtina výborne chvalitebné dobre nedostatečné výborne chvalitebné dobre nedostatečné Známka Známka Příklad 1.3. Vytvořte variační řady známek z matematiky a angličtiny pouze a) pro ženy, b) pro muže. a) Variační řada známek z matematiky pro ženy pohlavi<-data2$pohlavi variacni_rada(X=matematika[pohlavi=='zena'], nazvy=nazvy.známek) ## nj Nj pj Fj ## výborne 5 5 0.50.5 ## chvalitebné 2 7 0.2 0.7 ## dobre 1 8 0.1 0.8 ## nedostatečné 2 10 0.2 1.0 Variační řada známek z angličtiny pro ženy variacni_rada(X=anglictina[pohlavi=='zena'], nazvy=nazvy.známek) ## nj Nj pj Fj ## výborne 4 4 0.4 0.4 ## chvalitebné 2 6 0.2 0.6 ## dobre 1 7 0.1 0.7 ## nedostatečné 3 10 0.3 1.0 4 b) Variační řada známek z matematiky pro muže variacni_rada(X=matematika[pohlavi=='muz'], nazvy=nazvy.známek) ## nj Nj pj Fj ## výborne 2 2 0.20.2 ## chvalitebné 1 3 0.1 0.3 ## dobre 1 4 0.1 0.4 ## nedostatečné 6 10 0.6 1.0 Variační řada známek z angličtiny pro muže variacni_rada(X=anglictina[pohlavi=='muz'], nazvy=nazvy.známek) ## nj Nj pj Fj ## výborne 0 0 0.00.0 ## chvalitebné 2 2 0.2 0.2 ## dobre 6 8 0.60.8 ## nedostatečné 2 10 0.2 1.0 Příklad 1.4. Nadále budeme pracovat s celým datovým souborem. Vytvoříme kontingenční tabulku simultánních absolutních četností znaků X a Y. K.Tab <- table(matematika, angličtina) K.Tab2 <- cbind(K.Tab, suma=apply(K.Tab, 1, sum)) (K.Tab3 <- rbind(K.Tab2, suma=apply(K.Tab2, 2, sum))) ## výborne chvalitebné dobre nedostatečné suma ## výborne 4 12 0 7 ## chvalitebné 0 2 1 0 3 ## dobre 0 0 1 12 ## nedostatečné 0 13 4 8 ## suma 4 4 7 5 20 Vidíme, že ve výběrovém souboru byli 4 studenti, kteří měli z obou předmětů "výborně", jeden student, který měl z matematiky "výborně"a z angličtiny "chvalitebně"atd. až 4 studenti, kteří z obou předmětů neprospěli. Příklad 1.5. Vytvořte kontingenční tabulku řádkově a sloupcově podmíněných relativních četností znaků X=Ma-tematika a Y=Angličtina. Tab <- table(matematika, angličtina) # Radkové podmíněné relativní četnosti round(prop.table(Tab, margin=l), digits=3) ## angličtina ## matematika výborne chvalitebné dobre nedostatečné ## výborne 0.571 0.143 0.286 0.000 ## chvalitebné 0.000 0.667 0.333 0.000 ## dobre 0.000 0.000 0.500 0.500 ## nedostatečné 0.000 0.125 0.375 0.500 Interpretace např. 2. sloupce ve 4. řádku: V souboru bylo 8 studentů, kteří neprospěli z matematiky. Mezi nimi byl jeden, který měl chvalitebně z angličtiny, což představuje 1/8 = 12.5%. 5 # Sloupcové podmíněné relativní četnosti round(prop.table(Tab, margin=2), digits=3) ## angličtina ## matematika výborne chvalitebné dobre nedostatečné ## výborne 1.000 0.250 0.286 0.000 ## chvalitebné 0.000 0.500 0.143 0.000 ## dobre 0.000 0.000 0.143 0.200 ## nedostatečné 0.000 0.250 0.429 0.800 Interpretace např. 4. řádku ve 2. sloupci: V souboru byli 4 studenti, kteří měli chvalitebně z angličtiny. Mezi nimi byl jeden, který neprospěl z matematiky, což představuje 1/4 = 25%. Intervalové rozdělení četností Práci s intervalovým rozdělením četností si ukážeme na datovém souboru lebky.txt. Popis datového souboru: Máme k dispozici údaje o rozměrech lebek staroegyptské populace. Jedná se o 216 mužů a 109 žen. Znak X ... nej větší délka mozkovny v mm (tj. přímá vzdálenost kraniometrických bodů glabella a opisthocranion) Znak Y ... nej větší šířka mozkovny v mm (tj. přímá vzdálenost kraniometrických bodů euryon dx a euryon sin) Znak Z ... pohlaví osoby (1-muž, 0-žena) Příklad 1.6. Načtěte soubor lebky.txt. Podle Sturgersova pravidla najděte optimální počet třídicích intervalů pro znaky X a Y a vhodně stanovte meze třídicích intervalů, a to zvlášť pro muže a zvlášť pro ženy. data <- read. delim('lebky.txt', sep='\ť, dec= ' . ' , header=F) names(data) <- c('délka', 'sirka', 'pohlavi') head(data) ## délka ## 1 188 ## 2 172 ## 3 176 ## 4 184 ## 5 183 ## 6 177 sirka pohlavi 145 muz 139 muz 138 muz 128 muz 139 muz 143 muz # Muzi data.M <- data[data$pohlavi=='muz',] n.M <- dim(data.M)[1] (Sturges.M <- round(l+3.3*loglO(n.M), digits=0)) ## [1] 9 Protože mužů je 216, podle Sturgersova pravidla je optimální počet třídicích intervalů 9. Musíme zjistit minimum a maximum, abychom vhodně stanovili meze třídicích intervalů: # Znak X = Délka lebky délka.M <- data.M$delka range(délka.M) ## [1] 164 199 max(delka.M) - min(delka.M) ## [1] 35 6 round((max(delka.M) - min(delka.M))/Sturges.M, digits=0) ## [1] 4 Pro znak X = Délka lebky je minimum 164 a maximum 199, rozsah těchto hodnot je 35 a ideální délka jednoho třídicího intervalu vyšla jako 199~164 « 4. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 163, horní mez posledního třídicího intervalu 199. Celkem třídicí intervaly pro znak X budou: (163,167), (167,171), ..., (195,199). # Znak Y = Sirka lebky sirka.M <- data.M$sirka range(sirka.M) ## [1] 124 149 max(sirka.M) - min(sirka.M) ## [1] 25 round((max(sirka.M) - min(sirka.M))/Sturges.M, digits=0) ## [1] 3 Pro znak Y = šířka lebky je minimum 124 a maximum 149, rozsah těchto hodnot je 25 a ideální délka jednoho třídicího intervalu vyšla jako 149~124 « 3. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 123, horní mez posledního třídicího intervalu 150. Celkem třídicí intervaly pro znak X budou: (123,126), (126,129), ..., (147,150). # Zeny data.F <- data[data$pohlavi=='zena',] n.F <- dim(data.F)[1] (Sturges.F <- round(l+3.3*loglO(n.F), digits=0)) ## [1] 8 Protože žen je 109, podle Sturgersova pravidla je optimální počet třídicích intervalů 8. Postup je analogický jako u mužů. # Znak X = Délka lebky délka.F <- data.F$delka range(délka.F) ## [1] 157 188 max(délka.F) - min(delka.F) ## [1] 31 round((max(délka.F) - min(delka.F))/Sturges.F, digits=0) ## [1] 4 Pro znak X = Délka lebky je minimum 157 a maximum 188, rozsah těchto hodnot je 31 a ideální délka jednoho třídicího intervalu vyšla jako 188~157 « 4. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 156, horní mez posledního třídicího intervalu 188. Celkem třídicí intervaly pro znak X budou: (156,160), (160,164), ..., (184,188). # Znak Y = Sirka lebky sirka.F <- data.F$sirka range(sirka.F) 7 ## [1] 118 146 max(sirka.F) - min(sirka.F) ## [1] 28 round((max(sirka.F) - min(sirka.F))/Sturges.F, digits=0) ## [1] 4 Pro znak Y = šířka lebky je minimum 118 a maximum 146, rozsah těchto hodnot je 28 a ideální délka jednoho třídicího intervalu vyšla jako 146~118 « 4. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 116, horní mez posledního třídicího intervalu 148. Celkem třídicí intervaly pro znak X budou: (116,120), (120,124), ..., (144,148). Příklad 1.7. Vytvořte histogram pro X a pro Y (s uvedenými absolutními a relativními četnostmi jednotlivých třídicích intervalů), a to zvlášť pro muže a zvlášť pro ženy. # Muzi # X=Delka lebky hist(delka.M, breaks=seq(163, 199, by=4), ylim=c(0,52), main='Histogram délky lebky u múzu1, xlab='Délka lebky', ylab='Početnosti', col='white', border='white', density=20, axes=F) abline(h=seq( 0, 60, by=10), col='grey80', lty=2) hist(delka.M, breaks=seq(163, 199, by=4), col='blue', border='darkblue', density=20, add=T) axis(l, at=seq(163, 199, by=4)) axis(2, at=seq( 0, 50, by=10)) abs.c <- hist(delka.M, breaks=seq(163, 199, by=4), plot=F)$counts stred <- hist(delka.M, breaks=seq(163, 199, by=4), plot=F)$mids rel.c <- round(abs.c/sum(abs.c)*100, 0) četnosti <- paste(abs.c, '; ', rel.c, sep='') text(stred, abs.c+2, četnosti, cex=0.8) #------------------------------------------------------------------------------------- # Y=Sirka lebky hranice <- seq(123, 150, by=3) hist(sirka.M, breaks=hranice, ylim=c(0,52), main='Histogram sirky lebky u múzu1, xlab='sirka lebky', ylab='Početnosti', col='white', border='white', density=20, axes=F) abline(h=seq(0, 60, by=10), col='grey80', lty=2) hist(sirka.M, breaks=hranice, col='blue', border='darkblue', density=20, add=T) axis(l, at=hranice) axis(2, at=seq( 0, 50, by=10)) abs.c <- hist(sirka.M, breaks=hranice, plot=F)$counts stred <- hist(sirka.M, breaks=hranice, plot=F)$mids rel.c <- round(abs.c/sum(abs.c)*100, 0) četnosti <- paste(abs.c, '; ', rel.c, sep='') text(stred, abs.c+2, četnosti, cex=0.8) 8 Histogram delky lebky u muzu Histogram sirky lebky u muzu 163 171 —I— 179 187 195 Délka lebky Pro ženy je postup analogický jako pro muže. Histogram delky lebky u zen i-1-1-1-1-1-1-1-1 156 160 164 168 172 176 180 184 188 46; 21%j_í£F-1 123 —I— 129 T "T 135 sirka lebky 141 —I— 147 Histogram sirky lebky u zen 33; 30% I-1-1-1-1-1-1-1-1 116 120 124 128 132 136 140 144 148 Délka lebky sirka lebky 9 99 Příklady k samostatnému řešení Příklad 1.8. Bodové rozdělení četností V severozápadním Skotsku byla provedena studie, která zkoumala výskyt krevních skupin. V oblasti Eskdale bylo náhodně vybráno 100 osob, v oblasti Annandale 125 osob a v oblasti Nithsdale 253 osob. Výsledky jsou uvedeny v tabulce: oblast Krevní skupina n j A B 0 AB Eskdale 33 6 56 5 100 Annandale 54 14 52 5 125 Nithsdale 98 35 115 ■5 253 a.k 185 55 223 15 478 Jako znak X označíme oblast (má 3 varianty: Eskdale, Annandale, Nithsdale) a jako znak Y označíme krevní skupinu (má 4 varianty: A, B, AB a 0). Data jsou uložena v souboru krevni_skupiny.txt. a) Vytvořte variační řadu znaku Y, a to pro všechny tři oblasti dohromady a pak pro každou zvlášť. b) Nakreslete sloupkový diagram a polygon absolutních četností znaku Y. c) Nakreslete výsečový diagram pro znak X. d) Vytvořte kontingenční tabulku sloupcově a poté řádkově podmíněných relativních četností znaků X, Y. Řešení: a) Variační řady: • Všechny tři oblasti dohromady ## nj Nj pj Fj ## A 185 185 0.3870 0.3870 ## B 55 240 0.1151 0.5021 ## 0 223 463 0.4665 0.9686 ## AB 15 478 0.0314 1.0000 • Eskdale ## nj Nj PJ Fj ## A 33 33 0 33 0.33 ## B 6 39 0 06 0.39 ## 0 56 95 0 56 0.95 ## AB 5 100 0 05 1.00 Annandale ## nj Nj PJ Fj ## A 54 54 0 432 0.432 ## B 14 68 0 112 0.544 ## Q 52 120 0 416 0.960 ## AB 5 125 0 040 1.000 • Nithsdale ## nj Nj pj Fj ## A 98 98 0.3874 0.3874 ## B 35 133 0.1383 0.5257 ## Q 115 248 0.4545 0.9802 ## AB 5 253 0.0198 1.0000 10 b) Sloupkový graf a polygon četností Sloupkový diagram - krevní skupiny Polygon četnosti - krevní skupiny A B O AB A B O AB Krevni skupina Krevni skupina c) Výsečový graf Vysecovy graf - zastoupeni oblasti Nithsdale; 52.93% d) Řádkově a sloupcově podmíněné četnosti • Tabulka řádkově podmíněných četností ## Krev.Skupina ## Oblast A B 0 AB ## Annandale 0.432 0.112 0.416 0.040 ## Eskdale 0.330 0.060 0.560 0.050 ## Nithsdale 0.387 0.138 0.455 0.020 11 4694 • Tabulka sloupcově podmíněných četností ## Krev.Skupina ## Oblast A B 0 AB ## Annandale 0.292 0.255 0.233 0.333 ## Eskdale 0.178 0.109 0.251 0.333 ## Nithsdale 0.530 0.636 0.516 0.333 Příklad 1.9. Intervalové rozdělení četností U 50 studentů a studentek byla zjišťována jejich hmotnost (znak X, v kg), výška (znak Y, v cm) a pohlaví (znak Z, 0.. .žena, 1.. .muž). Data jsou uložena v souboru vys_vah.txt. a) Podle Sturgesova pravidla najděte optimální počet třídicích intervalů pro znaky X a Y a vhodně stanovte meze třídicích intervalů. b) Vytvořte histogram pro X a pro Y (s uvedenými absolutními a relativními četnostmi jednotlivých třídicích intervalů). Řešení a) Protože osob je 50, podle Sturgesova pravidla je optimální počet třídicích intervalů 7. Musíme zjistit minimum a maximum, abychom vhodně stanovili meze třídicích intervalů. Pro znak X je minimum 51 a maximum 90. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 50, horní mez posledního třídicího intervalu 92. Délka třídicích intervalů je tedy 92~50 = 6. Celkem třídicí intervaly pro znak X budou: (50; 56), (56; 62), (86; 92). Pro znak Y je minimum 160 a maximum 192. Jeví se vhodné dolní mez prvního třídicího intervalu zvolit 159, horní mez posledního třídicího intervalu 194. Délka třídicích intervalů je tedy 194~159 = 5. Celkem třídicí intervaly pro znak Y budou: (159; 164), (164; 169), (189; 194). b) Histogramy pro váhu a výšku Histogram - vaha studentu a studentek Histogram - vyska studentu a studentek 1-1-1-1-1-1-1-1 1-1-1-1-1-1-1-1 50 56 62 68 74 80 86 92 159 164 169 174 179 184 189 194 vaha vyska Poznámka: Tytéž úkoly lze řešit zvlášť pro muže a zvlášť pro ženy. 12 2 Výpočet číselných charakteristik jednorozměrného a dvourozměrného datového souboru Přehled použitých funkcí: data.frame, apply, library, round, cramersV, read.delim, source, head, names, factor, quantile, boxplot, cor, dotplot, abline, length, mean, var, sqrt, skewness, kurtosis, cbind. Příklad 2.1. U 100 náhodně vybraných domácností byl zjiťován způsob zásobování bramborami (znak X, varianty 1 = vlastní sklep, 2 = jinde, 3 = nákup) a bydliště (znak Y, varianty 1 = velké město, 2 = malé město, 3 = vesnice). velké město malé město vesnice vlastní sklep 13 15 14 jinde 11 7 2 nákup 19 9 10 a) Pro oba znaky určíme modus. b) Vypočteme Cramérův koeficient znaků X, Y. a) Stanovení modu (data <- data.frame(velke.mesto=c(13,11,19), male.mesto=c(15,7,9), vesnice=c(14,2,10), row.names=c('sklep','j inde','nakup'))) ## velké.město male.město vesnice ## sklep 13 15 14 ## jinde 11 7 2 ## nakup 19 9 10 apply(data,1,sum) ## sklep j inde nakup ## 42 20 38 apply(data,2,sum) ## velké.město male.město vesnice ## 43 31 26 Znak X má modus 1, tj. nejvíce domácností skladuje brambory ve vlastním sklepě a znak Y má také modus 1, tj. nejvíce domácností bydlí ve velkém městě. b) Výpočet Cramérova koeficientu Hodnotu Cramérova koeficientu vypočítáme pomocí funkce cramersV, která je součástí knihovny Isr. Nejprve tedy musíme nainstalovat tuto knihovnu (Packages —> Install —> Isr —> Install) a následněji načíst (library(lsr)). Teprve potom můžeme funcki cramersV() použít na naši datovou tabulku a Cramérův koeficient dopočítat. library(Isr) round(cramersV(data), digits=3) ## [1] 0.179 13 Cramérův koeficient nabývá hodnoty 0.179, tedy mezi způsobem zásobování bramborami a bydlištěm domácnosti existuje jen slabá závislost - viz následující tabulka: Cramérův koeficient interpretace 0 - 0.1 zanedbatelná závislost 0.1 - 0.3 slabá závislost 0.3 - 0.7 střední závislost 0.7 - 1 silná závislost Příklad 2.2. Otevřeme datový soubor znamky.txt. a) Pro známky z matematiky a angličtiny vypočteme medián, dolní a horní kvartil, kvartilovou odchylku a vytvoříme krabicový diagram. b) Vypočteme Spearmanův korelační koeficient známek z matematiky a angličtiny pro všechny studenty, pak zvlášť pro muže a zvlášť pro ženy. Získané výsledky budeme interpretovat. a) data <- read. delim('známky .txť , sep='\ť, dec= ' . ' ,header=F) source('AS-funkce.R') head(data) ## VI V2 V3 ## 1 2 2 0 ## 2 1 3 1 ## 3 4 3 1 ## 4 1 1 0 ## 5 1 2 1 ## 6 4 4 1 names(data) <- c('matematika', 'angličtina', 'pohlavi') f3 <- factor(data$pohlavi, levels=c(0,1), labels=c('zena','muz')) data[,3] <- f3 head(data) ## matematika angličtina pohlavi ##1 2 2 zena ##2 1 3 muz ## 3 4 3 muz ##4 1 1 zena ##5 1 2 muz ## 6 4 4 muz matematika <- data$matematika angličtina <- data$anglictina pohlavi <- data$pohlavi q.M <- quantile(matematika, probs=c(0.5,0.25,0.75), type=2) #type=5 q.A <- quantile(angličtina, probs=c(0.5,0.25,0.75), type=2) iqr.M <- q.M[3]-q.M[2] iqr.A <- q.A[3]-q.A[2] (tabulka<-data.frame(median=c(q.M[l],q.A[l]), kvl=c(q.M[2],q.A[2]), kv3=c(q.M[3],q.A[3]), IQR=c(iqr.M, iqr.A), row.names=c('matematika','angličtina'))) 14 ## medián kvl kv3 IQR ## matematika 2.5 1 4.0 3.0 ## angličtina 3.0 2 3.5 1.5 boxplot(matematika, angličtina, main='Krabicový graf dvou proměnných', names=c('matematika','angličtina'), ylab='známka', ylim=c(0,5), border='darkgreen', col='darkolivegreenl') Krabicový graf dvou proměnných m - -1- matematika -1- angličtina b) cor(matematika, angličtina, method='spearman') ## [1] 0.6884422 cor(matematika[pohlavi==lzena'], angličtina[pohlavi=='zena'], method='spearman') ## [1] 0.8603138 cor(matematika[pohlavi=='muz'], angličtina[pohlavi=='muz'], method='spearman') ## [1] 0.3735437 Vidíme, že nejsilnější přímá pořadová závislost mezi známkami z matematiky a angličtiny je u žen, r s = 0.86. U mužů je tato závislost mnohem slabší, r s = 0.37. U žen tedy dochází k tomu, že se sdružují podobné známky z obou předmětů, zatímco u mužů se projevuje spíše tendence k různým známkám. Je to zřetelně vidět na dvourozměrných tečkových diagramech. dotplot(matematika[pohlavi==lzena'], angličtina[pohlavi=='zena'], main='Tečkový graf známek - Zeny', xlab='matematika', ylab='angličtina', col='darkgreen', bg='darkolivegreenl', xlim=c(l,4), ylim=c(l,4)) abline(v=seq(l,4,by=0.5), col='grey80', lty=2) abline(h=seq(l,4,by=0.5), col='grey80', lty=2) 15 dotplot(matematika[pohlavi=='muz'], angličtina[pohlavi=='muz'], main='Tečkový graf známek - Muzi1, xlab='matematika', ylab='angličtina', col='darkgreen', bg='darkolivegreenl', xlim=c(l,4), ylim=c(l,4)) abline(v=seq(l,4,by=0.5), col='grey80', lty=2) abline(h=seq(l,4,by=0.5), col='grey80', lty=2) Tečkový graf známek - Zeny - a> i- i- T- i- i- 1.0 1.5 2.0 2.5 3.0 3.5 4.0 matematika Tečkový graf známek - Muzi [ 1.0 1.5 2.0 2.5 3.0 3.5 4.0 matematika Význam hodnot Spearmanova (i Pearsonova) koeficientu korelace je popsán v tabulce: Abs.hod. korel.koef. Interpretace hodnoty 0 pořadová (lineární) nezávislost (0;0.1) velmi nízký stupeň závislosti [0.1; 0.3) nízký stupeň závislosti [0.30; 0.50) mírný stupeň závislosti [0.50; 0.70) význačný stupeň závislosti [0.70; 0.90) vysoký stupeň závislosti [0.90; 1) velmi vysoký stupeň závislosti 1 úplná pořadová (lineární) závislost Podle výše uvedené tabulky existuje mezi známkami z matematiky a známkami z angličtiny význačný stupeň přímé pořadové závislosti (rs = 0.69), dále v případě žen existuje mezi známkami z matematiky a z angličtiny vysoký stupeň přímé pořadové závislosti (rs = 0.86), zatímco u mužů existuje mezi známkami z matematiky a z angličtiny pouze mírný stupeň přímé pořadové závislosti (rs = 0.37). Příklad 2.3. Otevřeme datový soubor lebky.txt. a) Pro největší délku a největší šířku mozkovny mužů vypočteme aritmetický průměr, rozptyl, směrodatnou odchylku, koeficient variace, šikmost a špičatost. b) Vypočítejte Pearsonův koeficient korelace největší délky a největší šířky mozkovny mužů. Dále vypočtěte kova-rianci těchto dvou znaků a nakreslete dvourozměrný tečkový diagram. a) library(el071) data <- read. delim( ' lebky .txt' , sep='\ť, dec= ' . ' , header=F) names(data) <- c('délka', 'sirka', 'pohlavi') head(data) 16 ## délka ## 1 188 ## 2 172 ## 3 176 ## 4 184 ## 5 183 ## 6 177 sirka pohlavi 145 muz 139 muz 138 muz 128 muz 139 muz 143 muz délka.M <- data$delka[data$pohlavi=='muz'] n <- length(délka.M) prumer.D <- mean(délka.M) rozptyl.D <- l/n*sum((délka.M-prumer.D)~2) sm.odch.D <- sqrt(rozptyl.D) koef.var.D <- sm.odch.D/mean(delka.M)*100 sikmost.D <- skewness(délka.M, type=2) spicatost.D <- kurtosis(delka.M, type=2) (tab.D <- round(data.frame(n=n, prumer=prumer.D, rozptyl=rozptyl.D, sm.odch=sm.odch.D, koef.var=koef.var.D, sikmost=sikmost.D, spicatost=spicatost.D), digits=4)) ## n prumer rozptyl sm.odch koef.var sikmost spicatost ## 1 216 182.0324 40.5777 6.3701 3.4994 -0.0551 -0.4511 Analogický postup zvolíme pro výpočty základních charakteristik pro šířku mozkovny mužů. Výsledné charakteristiky pro obě proměnné sloučíme do jedné tabulky. ## n prumer rozptyl sm.odch koef.var sikmost spicatost ## délka 216 182.0324 40.5777 6.3701 3.4994 -0.0551 -0.4511 ## sirka 216 137.1852 23.1694 4.8135 3.5087 0.0853 -0.2485 b) Výpočet Pearsonova korelačního koeficientu cor(delka.M, sirka.M, method='pearson') ## [1] 0.168157 Vidíme, že mezi délkou mozkovny a šířkou mozkovky u mužů existuje nízký stupeň přímé lineární závislosti. Výpočet kovariance kovariance <- sum((delka.M-prumer.D)*(sirka.M-prumer.S))/n round(kovariance, 4) ## [1] 5.156 Tečkový diagram dotplot(delka.M, sirka.M, main='Tečkový graf délky a sirky lebky múzu1, xlab='délka lebky', ylab='sirka lebky', col='darkgreen', bg='darkolivegreenl') abline(v=seq(160,200,by=5), col='grey80', lty=2) abline(h=seq(120,145,by=5), col='grey80', lty=2) 17 Tečkový graf délky a sirky lebky muzu oo o o o o.o o ---•--! oo o! o! o ■ O O' oo !o o ! ooqo o o o ooo o oo o O O O OO O I---O--0O-G OO OO O OOOOOOOOOO O OOO OO OO 0000 o oo oo oOo oo! O OO OOOOOO O O OiOO ■-e©--6e-©o-!---©coco ■--©-: OO o o o o , o o oo ooooooo o OOO O ' o ooo o í o o! O i o o oo o o -- 0 je střední počet těchto událostí. Píšeme X ~ Po(A). Pravděpodobnostní funkce Poissonova rozdělení má tvar ÍXX —-e"A pro x=0,l,...; , , x\ (13) 0 jinak. Distribuční funkce Poissonova rozdělení má tvar ^) = E7Te"A- (14) t=o l- Příklad 4.2. Nakreslete graf pravděpodobnostní funkce a distribuční funkce náhodné veličiny X ~ Po(5). #graf pravděpodobnostní funkce x <- 0:16 px <- dpois(x, lambda=5) plot(x, px, type='h', main='X~Po(5)', ylab='pravděpodobnostní funkce', xlab='x') points(x, px, col='reď, pch=19, cex=0.8) #graf distribuční funkce Fx <- ppois(x, lambda=5) n <- length(Fx) plot(x, Fx, type='n', main='X~Po(5)', ylab='distribuční funkce', xlab='x', xlim=c(-l,n), ylim=c(0,1)) segments(x, Fx, x+1, Fx) arrows(0, 0, -1, 0, length=0.1) 26 arrows(n-1, 1, n, 1, length=0.1) points(x, Fx, col='red', pch=19, cex=0.8) points(x, c(0, Fx[-n]), col='red', bg='white', pch=21, cex=0.8) X~Po(5) X~Po(5) 0 5 10 15 0 5 10 15 x x Príklad 4.3. Při provozu balicího automatu vznikají během směny náhodné poruchy, které se řídí rozdělením Po(2). Jaká je pravděpodobnost, že během směny dojde k alespoň jedné poruše? X ... počet poruch během směny; X ~ Po(2); 2° P(X > 1) = 1-P(X < 1) = 1-P(X < 0) = 1 - P(X = 0) = 1 - — e"2 = 0.8647. l-dpois(0, lambda=2) ## [1] 0.8646647 l-ppois(0, lambda=2) ## [1] 0.8646647 27 4.3 Exponenciální rozdělení Exp(A) Náhodná veličina X udává dobu čekání na příchod nějaké události, která se může dostavit každým okamžikem se stejnou šancí bez ohledu na dosud pročekanou dobu. (Jde o tzv. čekání bez paměti.) Přitom j vyjadřuje střední dobu čekání. Náhodná veličina X ~ Exp(A) má hustotu \e pro x > 0; 0 jinak. Příklad 4.4. Nakreslete graf hustoty a distribuční funkce náhodné veličiny X ~ Exp(2). (15) #graf hustoty x <- seq(from=0, to=2.5, length=512) fx <- dexp(x, rate=2) plot(x, fx, type='ľ, main='X~Exp(2) ' , xlab='x', col='reď) #graf distribuční funkce Fx <- pexp(x, rate=2) plot(x, Fx, main='X~Exp(2)', xlab='x', ylab='distribuční funkce1, type='ľ, col='reď) X~Exp(2) X~Exp(2) Příklad 4.5. Doba do ukončení opravy v opravně obuvi je náhodná veličina, která se řídí exponenciálním rozdělením se střední dobou opravy 3 dny. Jaká je pravděpodobnost, že oprava bude ukončena do dvou dnů? X ~ Exp(l/3); 1 P(X<2) = / -e"*da; = [-e"*]Q = 1 - e"3 = 0.4866. /o 3 pexp(2, rate=l/3) ## [1] 0.4865829 28 Příklad 4.6. Doba (v hodinách), která uplyne mezi dvěma naléhavými příjmy v jisté nemocnici, se řídí exponenciálním rozdělením se střední dobou čekání 2h. Jaká je pravděpodobnost, že uplyne více než 5 h bez naléhavého příjmu? X ~ Exp(l/2); ľ°° 1 2 5 P(X > 5) = P(X > 5) = / -e-*dx = [-e"f ] = e" = 0.0821. l-pexp(5, rate=l/2) ## [1] 0.082085 4.4 Normální rozdělení N(fi, a2) Náhodná veličina X ~ (/i, cr2) má hustotu 1 (*-n)2 /(x) = —=e-^~. (16) CrV27T Pro /i = 0 a cr2 = 1 se jedná o standardizované normální rozdělení, píšeme [7 ~ N(0,1). Hustota pravděpodobnosti má v tomto případě tvar f(x) = J=e-^. (17) y žír Příklad 4.7. Nakreslete graf hustoty a distribuční funkce náhodné veličiny X ~ N(—l, 2). #graf hustoty x <- seq(from=-5, to=3, length=512) fx <- dnorm(x, mean=-l, sd=sqrt(2)) plot(x, fx, type='ľ, main=bquote(paste('X " N( ',mu,'=-l, 1,sigma~2,1=2 )')), xlim=c(-4.5,2.5), col=1red1) #graf distribuční funkce Fx <- pnorm(x, mean=-l, sd=sqrt(2)) plot(x, Fx, type='ľ, main=bquote(paste('X ~ N( ',mu,'=-1, ',sigma~2,'=2 )')), xlim=c(-4.5,2.5), ylab= 'distribuční funkce', col='reď) 29 Příklad 4.8. Výsledky u přijímacích zkoušek na jistou VŠ jsou normálně rozděleny s parametry jj, = 550 bodů, a = 100 bodů. S jakou pravděpodobností bude mít náhodně vybraný uchazeč alespoň 600 bodů? X ... výsledek náhodně vybraného uchazeče; X ~ 7V(550,1002). l-pnorm(600, mean=550, sd=100) ## [1] 0.3085375 l-pnorm(0.5, mean=0, sd=l) ## [1] 0.3085375 Příklad 4.9. Životnost baterie v hodinách je náhodná veličina, která má normální rozdělení se střední hodnotou 300 hodin a směrodatnou odchylkou 35 hodin. Jaká je pravděpodobnost, že náhodně vybraná baterie bude mít životnost a) alespoň 320 hodin? b) nejvýše 310 hodin? a) ## [1] 0.2838546 b) ## [1] 0.6124515 Příklad 4.10. Na výrobní lince jsou automaticky baleny balíčky rýže o deklarované hmotnosti 1000 g. Působením náhodných vlivů hmotnost balíčků kolísá. Lze ji považovat za náhodnou veličinu, která se řídí normálním rozdělením se střední hodnotou 996 g a směrodatnou odchylkou 18 g. Jaká je pravděpodobnost, že náhodně vybraný balíček rýže neprojde výstupní kontrolou, jestliže je povolená tolerance ±30 g od deklarované hmotnosti 1000 g? ## [1] 0.1037604 Příklad 4.11. Nakreslete graf hustoty dvourozměrného standardizovaného normálního rozdělení. source('AS-funkce.R') x <- seq(from=-3, to=3, length=40) y <- seq(from=-3, to=3, length=40) nx <- length(x) ny <- length(y) z <- matrix(NA, nrow=nx, ncol=ny) for(i in l:nx){ for(j in l:ny){ 600 - 550 ) = 1 - P(Č7 < 0.5) = 1 - $(0.5) = 1-0.69146 = 0.3085. P(X £ (970; 1030)) = 1 - P(970 < X < 1030) = 1 - P(970 < X < 1030) 30 z[i,j] <- norm2(x[i], y[j], mul=0, mu2=0, sigmal=l, sigma2=l) } } color <- terrain.colors(12) stredy <- (z[-l, -1] + z[-l, -ncol(z)] + z[-nrow(z), -1] + z[-nrow(z), -ncol(z)])/4 stredy.col <- cut(stredy, 12) persp(x, y, z, col = color [stredy.col], phi = 30, theta = -45) 31 5 Výpočet číselných charakteristik náhodných veličin, aplikace Moi-vreovy — Laplaceovy věty 5.1 Kvantily vybraných spojitých rozdělení a-kvantil náhodné veličiny X značíme xa. Normální rozdělení N(p,, a2) Náhodná veličina X ~ N(fi, a2) má hustotu f(x) = —E=e ■ 0\j žít Pro /i = 0, a2 = 1 se jedná o standardizované normální rozdělení, píšeme U ~ ÍV(0,1). Hustota standardizovaného normálního rozdělení má v tomto případě tvar 1 _«íi y Zít a-kvantil standardizovaného normálního rozdělení značíme ua. Standardizované normální rozdělení je symetrické okolo nuly, proto pro kvantily tohoto rozdělení platí vztah Poznámka: Vyjádření hustot následujících tří rozdělení je příliš složité, lze ho najít např. v příloze A skript Marie Budíková, Pavel Osecký, Štěpán Mikoláš: Teorie pravděpodobnosti a matematická statistika. Sbírka příkladů. MU Brno 2007. X2 rozdělení s n stupni volnosti x'2(n) Nechť Xi,..., Xn jsou stochasticky nezávislé náhodné veličiny, Xi ~ ÍV(0,1), i = 1,..., n. Pak náhodná veličina x = X2 +... + x2n má x'2 rozdělení s n stupni volnosti X~X\n). a-kvantil x'2 rozdělení s n stupni volnosti značíme Xa(.n)-Studentovo rozdělení s n stupni volnosti t(n) Nechť X\, X'2 jsou stochasticky nezávislé náhodné veličiny, X\ ~ ÍV(0,1), X2 ~ x'2(n)- Pak náhodná veličina X- Xl má Studentovo rozdělení s n stupni volnosti X ~ t(n). a-kvantil Studentova rozdělení s n stupni volnosti značíme ta(n). Studentovo rozdělení je symetrické okolo nuly, proto pro kvantily tohoto rozdělení platí vztah ta(n) = -íi_a(n). 32 Fisherovo-Snedecorovo rozdělení s n\ a n2 stupni volnosti F(ni,n2) Nechť X\,..., Xn jsou stochasticky nezávislé náhodné veličiny, Xi ~ X2(ni)i i= lj 2. Pak náhodná veličina má Fishorovo rozdělení s ni a 112 stupni volnosti X ~ F(ni,n2). a-kvantil Fisherova rozdělení sni a «2 stupni volnosti značíme F„(ni,n2) Pro kvantily Fisherova rozdělení platí následující vztah Fa(n1,n2) Fi-a(n1,n2y Příklad 5.1. Najděte medián a horní a dolní kvartil náhodné veličiny U ~ N(0,1). qnorm(0.5) ## [1] 0 qnorm(0.25) ## [1] -0.6744898 qnorm(0.75) ## [1] 0.6744898 Příklad 5.2. Najděte dolní kvartil náhodné veličiny X ~ iV(3,5). qnorm(0.25, 3, sqrt(5)) ## [1] 1.491795 Příklad 5.3. Určete kvantil Xo.025(25)-qchisq(0.025, 25) ## [1] 13.11972 Příklad 5.4. Určete kvantily ío.gg(30) a í0.05(14). qt(0.99, 30) ## [1] 2.457262 qt(0.05, 14) ## [1] -1.76131 Příklad 5.5. Určete kvantily F0.975(5,20) a F0.05(2,10). 33 qf(0.975, 5,20) ## [1] 3.289056 qf(0.05, 2,10) ## [1] 0.0515573 5.2 Výpočet střední hodnoty a rozptylu diskrétních náhodných veličin Příklad 5.6. Postupně se zkouší spolehlivost čtyř přístrojů. Další se zkouší jen tehdy, když předchozí je spolehlivý. Každý z přístrojů vydrží zkoušku s pravděpodobností 0.8. Náhodná veličina X udává počet zkoušených přístrojů. Vypočtěte střední hodnotu a rozptyl náhodné veličiny X. X nabývá hodnot 1, 2, 3, 4 a její pravděpodobnostní funkce je vr(l) = 0.2, tt(2) = 0.8*0.2 = 0.16, tt(3) = 0.82*0.2 = 0.128, tt(4) = 0.83 * 0.2 + 0.84 = 0.512, tt(0) = 0 jinak. E(X) = 1 * 0.2 + 2 * 0.16 + 3 * 0.128 + 4 * 0.512 = 2.952 D(X) = l2 * 0.2 + 22 * 0.16 + 32 * 0.128 + 42 * 0.512-2.95222 = 1.4697 x <- 1:4 pi <- c(0.2, 0.8*0.2, 0.8~2*0.2, 0.8~3*0.2+0.8~4) (EX <- sum(x*pi)) ## [1] 2.952 (DX <- sum(x~2*pi)-EX~2) ## [1] 1.469696 Příklad k samostatnému řešení Příklad 5.7. Náhodná veličina X udává počet ok při hodu kostkou. Vypočtěte její střední hodnotu a rozptyl, x <- 1:6 pi <- rep(l/6, 6) (EX <- sum(x*pi)) ## [1] 3.5 (DX <- sum(x~2*pi)-EX~2) ## [1] 2.916667 5.3 Výpočet koeficientu korelace diskrétních náhodných veličin Příklad 5.8. Náhodná veličina X udává příjem manžela (v tisících dolarů) a náhodná veličina Y příjem manželky (v tisících dolarů). Je známa simultánní pravděpodobnostní funkce ir(x,y) diskrétního náhodného vektoru (X,Y): tt(10, 10) = 0.2, tt(10, 20) = 0.04, vr(10, 30) = 0.01, vr(10, 40) = 0, vr(20,10) = 0.1, vr(20, 20) = 0.36, vr(20, 30) = 0.09, vr(20,40) = 0, 7r(30,10) = 0, vr(30,20) = 0.05, vr(30,30) = 0.1, vr(30,40) = 0, vr(40,10) = 0, vr(40,20) = 0, 34 7r(40, 30) = 0, 7r(40, 40) = 0.05, ir(x, y) = 0 jinak. Vypočtěte koeficient korelace příjmů manžela a manželky. Náhodná veličina X i náhodná veličina Y nabývají hodnot 10, 20, 30, 40. Stanovíme hodnoty marginálních pravděpodobnostních funkcí: 7ri(10) = 0.25, 7ri(20) = 0.55, 7ri(30) = 0.15, 7ri(40) = 0.05, tti(x) = 0 jinak, vr2(10) = 0.3, vr2(20) = 0.45, vr2(30) = 0.2, 7r2(10) = 0.05, tt2(í/) = 0 jinak. Všechny hodnoty si zapíšeme do tabulky simultánních a marginálních pravděpodobností. Tabulka pravděpodobnostních funkcí n(X, Y) X - příjem manžela Y - příjem manželky 10 20 30 40 E 10 0.2 0.04 0.01 0 0.25 20 0.1 0.36 0.09 0 0.55 30 0 0.05 0.1 0 0.15 40 0 0 0 0.05 0.05 E 0.3 0.45 0.2 0.05 1 Spočteme E(X) = 20, E(Y) = 20, D(X) = 60, D(Y) = 70. Dosazením do vzorce pro výpočet kovariance zjistíme, 49 že C(X,Y) = 49, tedy koeficient korelace R(X,Y) = ^_ = 0.76. x <- c(10, 20, 30, 40) y <- c(10, 20, 30, 40) n <- length(x) pi <- data.frame(Deset= c(0.2, 0.1, 0, 0), Dvacet= c(0.04, 0.36, 0.05, 0), Tricet= c(0.01, 0.09, 0.1, 0), Ctyricet=c(0, 0, 0, 0.05), row.names=c('Deset', 'Dvacet', 'Třicet', 'Čtyřicet')) pix <- apply(pi, 1, sum) piy <- apply(pi, 2, sum) (EX <- sum(x*pix)) ## [1] 20 (EY <- sum(y*piy)) ## [1] 20 (DX <- sum(x~2*pix)-EX~2) ## [1] 60 (DY <- sum(y~2*piy)-EY~2) ## [1] 70 (CXY <- sum(c((x-EX)*(y-EY)[1], (x-EX)*(y-EY)[2], (x-EX)*(y-EY)[3], (x-EX)*(y-EY) [4]) * c(as.matrix(pi)))) ## [1] 49 (RXY <- CXY/(sqrt(DX)*sqrt(DY))) ## [1] 0.7560864 35 Příklady k samostatnému řešení Příklad 5.9. Objektem zájmu rozsáhlé studie bylo sledování pohřebního rituálu dnes již vymřelého, ale v minulosti velmi dlouho přetrvávajícího a rozsáhlého jihoamerického kmene. Součástí pohřebního rituálu tohoto kmene bylo odsekávání článků prstů na rukou a nohou zemřelého a jejich následné obětování bohům jako dar, aby zemřelého přijali mezi sebe. Zemřelému tak byl na ruce odetnut buď jeden nebo dva prsty a na noze tři nebo čtyři prsty. Dále bylo zjištěno, že domorodci odtínali jeden prst na ruce a tři prsty na noze zemřelého s pravděpodobností 0.1, dva prsty na ruce a tři prsty na noze s pravděpodobností 0.3, jeden prst na ruce a čtyři prsty na noze s pravděpodobností 0.35 a dva prsty na ruce a čtyři prsty na noze s pravděpodobností 0.25. Určete korelaci znaků X - počet odetnutých prstů na rukou &Y- počet odetnutých prstů na nohou. ## EX EY DX DY CXY RXY ## Charakteristiky 1.6 3.55 0.24 0.2475 -0.08 -0.3282 Z tabulky výsledků vidíme, že střední hodnota počtu prstů odetnutých na rukou je 1.6, zatímco střední hodnota počtu prstů odetnutých na nohou je 3.55. Rozptyl počtu prstů odetnutých na rukou je 0.24 a rozptyl počtu prstů odetnutých na nohou je 0.2475. Kovarince mezi znaky X & Y nabývá hodnoty -0.08. Hodnota korelačního koeficientu vyšla -0.3282, což značí, že mezi počtem prstů odetnutých na rukou a počtem prstů odetnutých na nohou existuje mírný stupeň nepřímé lineární závislosti. Příklad 5.10. Diskrétní náhodný vektor (X, Y) má simultánní pravděpodobnostní funkci s hodnotami 7r(0, —1) = c, 7r(0, 0) = tt(0, 1) = tt(1, -1) = vr(2, -1) = 0, tt(1, 0) = tt(1, 1) = vr(2,1) = 2c, vr(2, 0) = 3c, tt(x, y) = 0 jinak. Určete konstantu c a vypočtěte R(X,Y). Tabulka pr. funkcí n (X, Y) X Y -1 0 1 E 0 c 0 0 c 1 0 2c 2c 4c 2 0 3c 2c 5c E c 5c 4c 10c=l ## EX EY DX DY CXY RXY ## Charakteristiky 1.4 0.3 0.44 0.41 0.18 0.4238 Střední hodnota náhodné veličiny X je 1.4, střední hodnota náhodné veličiny Y je 0.3. Rozptyl náhodné veličiny X nabývá hodnoty 0.44, rozptyl náhodné veličiny Y nabývá hodnoty 0.41. Kovariance mezi veličinami X a Y je 0.18 a korelační koeficient nabývá hodnoty 0.4238, což značí, že mezi znaky X a Y existuje mírný stupeň přímé lineární závislosti. Příklad 5.11. Zkoumali jsme potomky kosmanů. Náhodná veličina X udává počet manželských potomků, které samice porodila a náhodná veličina Y počet nemanželských potomků, které samice porodila. Je známa simultánní pravděpodobnostní funkce ir(x,y) diskrétního náhodného vektoru (X,Y): Tabulka simultánní pstní fce tt(X,Y) X - počet manž.p. Y - počet nemanž.p. 1 2 3 1 0.2 0.04 0.01 2 0.15 0.36 0.09 3 0.05 0.1 0.0 Vypočtěte koeficient korelace manželských a nemanželských potomků. 36 ## EX EY DX DY CXY RXY ## Charakteristiky 1.9 1.7 0.39 0.41 0.11 0.2751 Střední hodnota počtu manželských potomků kosmanů je 1.9, střední hodnota počtu nemanželských potomků je 1.6. Rozptyl manželských potomků je 0.39, rozptyl nemanželských potomků je 0.41. Kovariance mezi počtem manželských a nemanželských potomků je 0.11. Korelační koeficient nabývá hodnoty 0.2751, což znamená, že mezi počtem manželských a nemanželských potomků kosmanů existuje nízký stupeň přímé lineární závislosti. 5.4 Aplikace Moivreovy-Laplaceovy věty Příklad 5.12. Pravděpodobnost úspěchu při jednom pokusu je 0.3. S jakou pravděpodobností lze tvrdit, že počet úspěchů ve 100 pokusech bude v mezích od 20 do 40? Výpočet proveďte a) přesně; b) pomocí aproximace normálním rozdělením. # a) sum(dbinom(20:40, 100, 0.3)) ## [1] 0.9786144 pbinom(40, 100, 0.3)-pbinom(19, 100, 0.3) ## [1] 0.9786144 # b) pnorm(40, 100*0.3, sqrt(100*0.3*0.7))-pnorm(19, 100*0.3, sqrt(100*0.3*0.7)) ## [1] 0.9772632 Příklad 5.13. Pravděpodobnost, že zakoupený elektrospotřebič bude vyžadovat opravu během záruční doby, je rovna 0.2. Jaká je pravděpodobnost, že během záruční doby bude nutno ze 400 prodaných spotřebičů opravit více než 96? Výpočet proveďte a) přesně; b) pomocí aproximace normálním rozdělením. # a) l-pbinom(96, 400, 0.2) ## [1] 0.02138855 l-pnorm(96, 400*0.2, sqrt(400*0.2*0.8)) ## [1] 0.02275013 Výsledek: ad a) 0.0246, ad b) 0.0228 Příklad k samostatnému řešení Příklad 5.14. Pravděpodobnost, že určitý typ výrobku má výrobní vadu, je 0.05. Jaká je pravděpodobnost, že ze série 1000 výrobků bude mít výrobní vadu nejvýše 70? Výpočet proveďte 37 a) přesně; b) pomocí aproximace normálním rozdělením. # a) pbinom(70, 1000, 0.05) ## [1] 0.9976697 # b) pnorm(70, 1000*0.05, sqrt(1000*0.05*0.95)) ## [1] 0.9981455 6 Základní pojmy matematické statistiky 6.1 Bodové odhady parametrů Příklad 6.1. Ve 12-ti náhodně vybraných internetových obchodech byly zjištěny následující ceny deskriptoru artefaktů (v Kč): 102, 99,106,103, 96, 98,100,105,103, 98,104,107. Těchto 12 hodnot považujeme za realizace náhodného výběru X\,..., X\2 z rozdělení, které má střední hodnotu /i a rozptyl a2. a) Určete nestranné bodové odhady neznámé střední hodnoty jj, a neznámého rozptylu a2. b) Najděte výběrovou distribuční funkci Fi2(x) a nakreslete její graf. ad a) Vypočteme realizaci výběrového průměru m = -7j(102 + 99 + • • • + 107) = 101.75 Kč Vypočteme realizaci výběrového rozptylu: s2 = ij- [(102 - 101.75)2 + (99 - 101.75)2 + • • • + (107 - 101.75)2] = 12.39 Kč2 x <- c(96, 98, 98, 99, 100, 102, 103, 103, 104, 105, 106, 107) n <- length(x) (m <- mean(x)) ## [1] 101.75 (s2 <- var(x)) ## [1] 12.38636 ad b) Pro usnadnění výpočtu hodnot výběrové distribuční funkce F\2{x) uspořádáme ceny podle velikosti: 96, 98, 98, 99, 100, 102, 103, 103, 104, 105, 106, 107. Číselnou osu rozdělíme na 11 intervalů a v každém intervalu 38 stanovíme hodnotu výběrové distribuční funkce: x < 96 : F12(x) = 0 96 < x < 98 : F12(x) 12 3 0.083 98 < x < 99 : F12(x) = — = 0.25 99 < x < 100 : Fi2(a;) 100 < x < 102 : F12(x 102 < a; < 103 : F12(x 103 < x < 104 : F12(x 104 < a; < 105 : F12(x 105 < x < 106 : Fi2(a; 106 < x < 107 : F12(x x > 107 : Fi2(a;) = 1 4 12 " 5_ ~ 12 _6_ ~ 12 _8_ ~ 12 _9_ ~ 12 10 ~ 12 11 " 12 0.3 = 0.416 = 0.5 = 0.6 = 0.75 = 0.83 = 0.916 # Výberová distribuční funkce t <- unique(sort(x)) y <- sort(x) nt <- length(t) četnost <- NULL for(i in l:nt){ četnost[i] <- sum(y<=t [i])} Fx <- cetnost/n round(Fx, digits=4) ## [1] 0.0833 0.2500 0.3333 0.4167 0.5000 0.6667 0.7500 0.8333 0.9167 1.0000 # graf výběrové distribuční funkce F(x) x <- c(min(t)-l,t, max(t)+l) y <- c(0,Fx,l) plot(x, y, type='n', xlab='x', ylab='F(x)', main='Výberová distribucni funkce') abline(h=seq(0,1,by=0.1), col='grey85') abline(v=seq(95, 108,by=2), col='grey85') lines (x, y, type= ' s ' , col='reď, lwd=2) arrows(96,0,95,0, col='reď, lwd=2, length=0.1) arrows(107,1,108,1, col='reď, lwd=2, length=0.1) 39 Výberová distribucni funkce co o o o ô 96 98 100 102 104 106 108 Příklad 6.2. Přírůstky cen akcií v % na burze v New Yorku u 10 náhodne vybraných společností dosáhly těchto hodnot: 10,16, 5,10,12, 8,4, 6, 5,4. a) Odhadněte střední hodnotu, rozptyl a směrodatnou odchylku růstu cen akcií. b) Odhadněte pravděpodobnost růstu cen akcií aspoň o 8.5%. ad a) x <- c(10, 16, 5, 10, 12, 8, 4, 6, 5, 4) x <- sort(x) n <- lengtli(x) s2 <- var(x) s <- sd(x) Tab <- data.frame(m=m, s2=s2, s=s, row.names='akcie') round(Tab, digits=2) ## m s2 s ## akcie 101.75 15.78 3.97 ad b) P(X > 8.5) P(X > 8.5) = 1 nx>8.5 n nx=8.5) pst <- sum(x>=8.5)/lengtli(x) pst <- l-sum(x<8.5)/lengtri(x) round(pst,4) ## [1] 0.4 40 Průměrný růst cen akcií odhadujeme na 8 % se směrodatnou odchylkou 3.97%. Dále, u 40 % akcií vzrostla cena alespoň o 8.5 %. Příklad 6.3. Bylo zkoumáno 9 vzorků půdy s různým obsahem fosforu (veličina X). Hodnoty veličiny Y označují obsah fosforu v obilných klíčcích (po 38dnech), jež vyrostly na těchto vzorcích půdy. číslo vzorku 1 2 3 4 5 6 7 8 9 X 1 4 5 9 11 13 23 23 28 Y 64 71 54 81 76 93 77 95 109 Těchto 9 dvojic hodnot považujeme za realizace náhodného výběru {X\, Y\),..., (Xg, Yg) z dvourozměrného rozdělení s kovariancí a\2 a koeficientem korelace p. Najděte bodové odhady kovariance a\2 a koeficientu korelace p. x <- c(l, 4, 5, 9, 11, 13, 23, 23, 28) y <- c(64, 71, 54, 81, 76, 93, 77, 95, 109) cov(x,y) ## [1] 130 cor(x,y) ## [1] 0.8049892 Výběrová kovariance veličin X, Y se realizuje hodnotou 130. Výběrový koeficient korelace veličin X, Y nabyl hodnoty 0.805, tedy mezi veličinami X, Y existuje silná přímá lineární závislost. 6.2 Intervalové odhady parametrů Nechť Xi,..., Xn je náhodný výber z rozdělení L(9), 9 je sledovaný parametr a a € (0,1). Interval (D, H) se nazývá 100(1 — a)% oboustranný interval spolehlivosti parametru 9, pokud pro každé 9 € 0 platí P(D < 9 < H) = 1 - a. Interval (D, oo) se nazývá 100(1 — a)% levostranný interval spolehlivosti parametru 9, pokud pro každé 9 € 0 platí P(D <9) = l-a. Interval (—oo,H) se nazývá 100(1 — a)% pravostranný interval spolehlivosti parametru 9, pokud pro každé 9 € 0 platí P(0 < H) = 1 - a. Číslo a se nazývá riziko, číslo 1 — a se nazývá spolehlivost. Příklad 6.4. Při kontrolních zkouškách životnosti 16 žárovek byl stanoven odhad m = 3000 h střední hodnoty jejich životnosti. Z dřívějších zkoušek je známo, že životnost žárovky se řídí normálním rozdělením se směrodatnou odchylkou a = 20 h. Vypočtěte a) 99% empirický interval spolehlivosti pro střední hodnotu životnosti; b) 90% levostranný empirický interval spolehlivosti pro střední hodnotu životnosti; c) 95% pravostranný empirický interval spolehlivosti pro střední hodnotu životnosti. Výsledek zaokrouhlete na jedno desetinné místo a vyjádřete v hodinách a minutách. 41 d = m - -^=Ui_Q = 3000 - -==2.57583 = 2987.1 Vň VIE on h = m- -^tia = 3000 + -==2.57583 = 3012.9 m <- 3000 s <- 20 n <- 16 # a) alpha <- 0.01 (dh <- m-s/sqrt(n)*qnorm(l-alpha/2)) ## [1] 2987.121 (hh <- m-s/sqrt(n)*qnorm(alpha/2)) ## [1] 3012.879 2987 h a 6 min < fi < 3012 h a 54 min s pravděpodobnostní 0.99. ad b) d = m- -^=ui-« = 3000 - 42=1.28155 = 2993.6 y/ň vTě alpha <- 0.1 (dh <- m-s/sqrt(n)*qnorm(l-alpha)) ## [1] 2993.592 2993 h a 36 min < fi s pravděpodobnostní 0.9. h = m- -^=ua = 3000 + —==1.95996 = 3008.2 Vň y/m alpha <- 0.05 (hh <- m-s/sqrt(n)*qnorm(alpha)) ## [1] 3008.224 3008 h a 13 min > fi s pravděpodobnostní 0.95. Užitečný odkaz: na adrese http://www.prevody-jednotek.cz je program, s jehož pomocí lze převádět různé fyzikální jednotky, v našem případě hodiny na minuty. 42 6.3 Testování hypotézy Příklad 6.5. Víme, že výška hochů ve věku 9.5 až 10 let má normální rozdělení s neznámou střední hodnotou jj, a známým rozptylem a2 = 39.112 cm2. Dětský lékař náhodně vybral 15 hochů uvedeného věku, změřil je a vypočítal realizaci výběrového průměru m = 139.13 cm. Podle jeho názoru by výška hochů v tomto věku neměla přesáhnout 142 cm s pravděpodobností 0.95. Lze tvrzení lékaře akceptovat? Testujeme Hq : jj, > 142 proti H\ : jj, < 142 na hladině významnosti a = 0.05. a) Test provedeme pomocí kritického oboru. Pro úlohy o střední hodnotě normálního rozdělení při známém rozptylu používáme pivotovou statistiku M — u U = —^- ~ AT(0,1). Testovací statistika bude mít tedy tvar M -c a za platnosti nulové hypotézy Hq se tato statistika bude řídit standardizovaným normálním rozdělením T0 ~ N(0,1). Vypočítáme realizaci testového kritéria: 139.13 - 142 t0 =---— = -1.7773. ^39.112 sigma <- sqrt(39.112) n <- 15 m <- 139.13 c <- 142 alpha <- 0.05 # ad a) tO <- (m-c)/(sigma/sqrt(n)) qnorm(alpha) ## [1] -1.644854 Stanovíme kritický obor: W G (—oo,ua) = (—oo,Mo.os) = (—oo, — «0.95) = (—00,—1.6449). Protože —1.7773 G W, nulovou hypotézu Hq zamítáme na hladině významnosti a = 0.05. Tvrzení lékaře lze tedy akceptovat s rizikem omylu 5 %. b) Test provedeme pomocí intervalu spolehlivosti. Meze 100(1 — 0;)% empirického pravostranného intervalu spolehlivosti pro střední hodnotu jj, při známém rozptylu cr2 jsou (—00, h) = (—00, m--j=ua). J n V našem případě dostáváme: h = 139.13 - V39^12Mo 05 = 139.13 + V39A£121.645 = 141.79. 43 (hh <- m-(sigma/sqrt(n))*qnorm(alpha)) ## [1] 141.7861 Protože 142 ^ (—oo; 141.79), Hq zamítáme na hladine významnosti 0.05. c) Test provedeme pomocí p-hodnoty. p = P(T0 < í0) = $(-1.7773) = 0.0378 (pval <- pnorm(tO)) ## [1] 0.03775549 Jelikož 0.0378 < 0.05, nulovou hypotézu zamítáme na hladině významnosti 0.05. 7 Ověřování normality a parametrické úlohy o jednom náhodném výběru z normálního rozdělení a dvourozměrného rozdělení 7.1 Grafické ověřování normality Příklad 7.1. Při nanášení tenkých kovových vrstev stříbra na polymerní materiál se vyžaduje, aby tloušťka vrstvy byla 0.020 /im. Pomocí atomové absorpční spektroskopie se zjistily hodnoty, jež jsou uvedeny v tabulce a uloženy v souboru vrstva_stribra.txt. Posuďte Q-Q grafem, zda se výsledky měření řídí normálním rozdělením. tloušťka 0.0212 0.0186 0.0192 0.0207 0.0200 0.0200 0.0190 0.0188 0.0208 0.0194 0.0188 0.0193 0.0204 0.0185 0.0187 0.0195 0.0191 0.0195 0.0199 0.0205 0.0189 0.0188 0.0199 0.0202 0.0208 data <- read.delimCvrstva_stribra.txt') stribro <- data$tloustka_vrstvy qqnorm(stribro, col='darkblue', xlab='teoreticky kvantiľ, ylab='pozorovaný kvantil', main='Q-Q graf') qqline(stribro, col='reď) 44 Q-Q graf teoreticky kvantil Dle vzhledu Q-Q grafu lze soudit, že data vykazují jen lehké odchylky od normality. 7.2 Testy normality Příklad 7.2. U 48 studentek VŠE v Praze byla zjišťována výška a obor studia (1 - národní hospodářství, 2 — informatika). Hodnoty jsou uloženy v souboru vyska.txt. Pomocí Lilieforsovy modifikace K-S testu, pomocí S-W testu a pomocí A-D testu testujte na hladině významnosti a = 0.05 hypotézu, že data pochází z normálního rozdělení. Pomocí Q-Q grafu posuďte vizuálně předpoklad normality. library(nortest) data <- read.tableCvyska.txt', header=T, row.names=NULL) head(data) ## vyska obor ## 1 165 1 ## 2 170 1 ## 3 170 1 ## 4 179 1 ## 5 170 1 ## 6 168 1 vyska <- data$vyska shapiro.test(vyska) ## ## Shapiro-Wilk normality test ## ## data: vyska ## W = 0.966, p-value = 0.176 lillie.test(vyska) ## ## Lilliefors (Kolmogorov-Smirnov) normality test 45 ## ## data: vyska ## D = 0.15562, p-value = 0.005258 ad.test(vyska) ## ## Anderson-Darling normality test ## ## data: vyska ## A = 0.66099, p-value = 0.07933 Shrnutí výsledků st <- shapiro.test(vyska) It <- lillie.test(vyska) at <- ad.test(vyska) (tab <- data.frame(test.statistika=round(c(lt$statistic, st$statistic, at$statistic),6), p.value=round(c(lt$p.value, st$p.value, at$p.value),6), rozhodnuti=c('zamitame', 'nezamitame', 'nezamitame'), row .names=c( ' Lillie .test' , 'S-W.tesť, 'A-D.test'))) ## test.statistika p.value rozhodnuti ## Lillie.test 0.155621 0.005258 zamitame ## S-W.test 0.965996 0.176031 nezamitame ## A-D.test 0.660990 0.079330 nezamitame Q-Q graf qqnorm(vyska, col='darkblue', xlab='teoreticky kvantil', ylab='pozorovaný kvantil', main='Q-Q graf) qqline(vyska, col='reď) Q-Q graf LO 00 —l-1 T- O o o co - ■■- o o o i-1-1-1-r -2-10 1 2 teoreticky kvantil 46 Tečky se řadí podél ideální přímky, normalita je jen lehce porušena. Příklad k samostatnému řešení Příklad 7.3. Testy normality a grafické ověření normality proveďte a) pro výšku studentek oboru národní hospodářství; b) pro výšku studentek oboru informatika. ad a) ## test.statistika p.value rozhodnuti ## Lillie.test 0.167473 0.042926 zamitame ## S-W.test 0.970969 0.606793 nezamitame ## A-D.test 0.419238 0.305268 nezamitame Vidíme, že Lilieforsova varianta K-S testu zamítá hypotézu o normalitě na hladině významnosti a = 0.05 (p-hodnota je menší než 0.05), zatímco S-W test hypotézu o normalitě nezamítá (p-hodnota je větší než 0.05). A-D test poskytne hodnotu testové statistiky 0.4192, odpovídající p-hodnota je 0.3053, tedy A-D test nezamítá hypotézu o normalitě na hladině významnosti a = 0.05. ad b) ## test.statistika p.value rozhodnuti ## Lillie.test 0.172301 0.123974 nezamitame ## S-W.test 0.922747 0.111924 nezamitame ## A-D.test 0.566019 0.123709 nezamitame V případě b) ani jeden z testů hypotézu o normalitě nezamítá na hladině významnosti a = 0.05. 47 7.3 Parametrické úlohy o jednom náhodném výběru z normálního rozdělení Upozornění: Pokud to povaha úlohy vyžaduje, proveďte test normality dat. Příklad 7.4. Vlastnosti výběrového průměru z normálního rozdělení: Předpokládejme, že velký ročník na vysoké škole má výsledky ze statistiky normálně rozděleny kolem střední hodnoty 72 bodů se směrodatnou odchylkou 9 bodů. Najděte pravděpodobnost, že průměr výsledků náhodného výběru 10 studentů bude větší než 80 bodů. Návod: Xi,...,Xio je náhodný výběr z N(72, 81). Počítáme P(M > 80), přičemž výběrový průměr M má normální rozdělení se střední hodnotou E(M) = \i = 72 a rozptylem D(M) = ^ = f± = 8.1. Tedy P(M > 80) = 1 - P(M < 80) = l-$(80), kde $(80) je hodnota distribuční funkce rozdělení AT(72, 8.1) v bodě 80. l-pnorm(80, 72, sqrt(8.1)) ## [1] 0.002470053 Příklad 7.5. Intervaly spolehlivosti pro parametry /i, a2 normálního rozdělení: Z populace stejně starých selat téhož plemene bylo vylosováno šest selat a po dobu půl roku jim byla podávána táž výkrmná dieta. Byly zaznamenávány průměrné denní přírůstky hmotnosti v dkg. Z dřívějších pokusů je známo, že v populaci mívají takové přírůstky normální rozdělení, avšak střední hodnota i rozptyl se měnívají. Přírůstky v dkg: 62, 54, 55, 60, 53, 58. a) Najděte 95% empirický levostranný interval spolehlivosti pro neznámou střední hodnotu jj, při neznámé směrodatné odchylce a. b) Najděte 95% empirický interval spolehlivosti pro směrodatnou odchylku a. Shapirův - Wilkův test normality: x <- c(62, 54, 55, 60, 53, 58) shapiro.test(x)$p.val ## [1] 0.6194994 P-hodnota S-W testu je 0.6195 > 0.05, tedy nulovou hypotézu o normálním rozdělení náhodného výběru nezamítáme na hladině významnosti a = 0.05. Výpočet intervalů spolehlivosti: ad a) x <- c(62, 54, 55, 60, 53, 58) m <- mean(x) s <- sd(x) n <- length(x) alpha <- 0.05 dd <- m-s/sqrt(n)*qt(1-alpha, n-1) print(paste('dd =', round(dd, 4))) ## [1] "dd = 54.0568" IS= (54.0568; oo) /i > 54.06 dkg s pravděpodobností 0.95. 48 ad b) dh <- (n-l)*s~2/qchisq(l-alpha/2, n-1) hh <- (n-l)*s~2/qchisq(alpha/2, n-1) print(paste('dh =', round(sqrt(dh), 4))) ## [1] "dh = 2.2332" print(paste('hh =', round(sqrt(hh), 4))) ## [1] "hh = 8.7747" IS= (2.2332; 8.7747) 2.23 dkg < a < 8.77 dkg s pravděpodobností 0.95. Příklad 7.6. Testování hypotézy o střední hodnotě jj,: Systematická chyba měřícího přístroje se eliminuje nastavením přístroje a měřením etalonu, jehož správná hodnota je jj, = 10.00. Nezávislými měřeními za stejných podmínek byly získány hodnoty: 10.24, 10.12, 9.91, 10.19, 9.78, 10.14, 9.86, 10.17, 10.05, které považujeme za realizace náhodného výběru rozsahu 9 z rozdělení N(fi,a2). Je možné při riziku 0.05 vysvětlit odchylky od hodnoty 10.00 působením náhodných vlivů? Shapirův - Wilkův test normality: x <- c(10.24, 10.12, 9.91, 10.19, 9.78, 10.14, 9.86, 10.17, 10.05) shapiro.test(x)$p.val ## [1] 0.2873252 P-hodnota S-W testu je 0.2873 > 0.05, tedy nulovou hypotézu o normálním rozdělení náhodného výběru nezamítáme na hladině významnosti a = 0.05. Testování nulové hypotézy: Na hladině významnosti a = 0.05 testujeme hypotézu Hq K testování použijeme jednovýběrový í-test. a) Testování pomocí kritického oboru : /i = 10 proti oboustranné alternativě H\ : jj, ^ 10. x <- c(10.24, 10.12, 9.91, 10.19, 9.78, 10.14, 9.86, 10.17, 10.05) alpha <- 0.05 m <- mean(x) s <- sd(x) n <- length(x) c <- 10 tO <- (m-c)/(s/sqrt(n)) wl <- qt(l-alpha/2, n-1) w2 <- qt(alpha/2, n-1) print(paste('tO =', round(tO, 4))) ## [1] "tO = 0.9426" print(paste('wl =', round(wl, 4))) 49 ## [1] "wl = 2.306 print(paste('w2 =', round(w2, 4))) ## [1] "w2 = -2.306" Testovací statistika íg nabývá hodnoty 0.9426, kritický obor má tvar W = (-co ; -2.3060) U (2.3060 ; oo) Protože tg ^ W, Hq nezamítáme na hladině významnosti a = 0.05. b) Testování pomocí intervalu spolehlivosti dh <- m-s/sqrt(n)*qt(l-alpha/2, n-1) hh <- m-s/sqrt(n)*qt(alpha/2, n-1) print(paste('dh =', round(dh, digits=4))) ## [1] "dh = 9.9261" print(paste('hh =', round(hh, digits=4))) ## [1] "hh = 10.1761" 95% empirický interval spolehlivosti pro střední hodnotu jj, má tvar IS = (9.9261; 10.1761). Protože c = 10 G IS, nulovou hypotézu Hq nezamítáme na hladině významnosti a = 0.05. c) Testování pomocí p-hodnoty p.val <- 2*min(pt(t0, n-1),1-pt(tO, n-1)) print(paste('p-hodnota =', round(p.val, 4))) ## [1] "p-hodnota = 0.3735" Protože p-hodnota = 0.3735 > 0.05, nulovou hypotézu Hq nezamítáme na hladině významnosti a = 0.05. Poznámka: K otestování nulové hypotézy o střední hodnotě jj, můžeme použít funkci t.test(x) s argumentem mu=10 (hodnota c z nulové hypotézy) a argumentem alternative='two.sideď (oboustranná alternativa). x <- c(10.24, 10.12, 9.91, 10.19, 9.78, 10.14, 9.86, 10.17, 10.05) t.test(x, mu=10, alternative='two.sideď) ## ## One Sample t-test ## ## data: x ## t = 0.94261, df = 8, p-value = 0.3735 ## alternative hypothesis: true mean is not equal to 10 ## 95 percent confidence interval: ## 9.926073 10.176149 ## sample estimates: ## mean of x ## 10.05111 50 Příklad 7.7. Testování hypotézy o směrodatné odchylce a: U 25 náhodně vybraných dvoulitrových lahví s nealkoholickým nápojem byl zjištěn přesný objem nápoje. Výběrový průměr činil m = 1.991 a výběrová směrodatná odchylka s = 0.11. Předpokládejme, že objem nápoje v láhvi je náhodná veličina s normálním rozdělením. Na hladině významnosti a = 0.05 ověřte tvrzení výrobce, že směrodatná odchylka je 0.081. Na hladině významnosti a = 0.05 testujeme hypotézu Hq : a = 0.08 proti oboustranné alternativě H\ : a ^ 0.08 neboli Hq : a2 = 0.0064 proti oboustranné alternativě H\ : a2 ^ 0.0064. K testování nulové hypotézy použijeme test o rozptylu. a) Testování pomocí kritického oboru Vypočteme realizaci testového kritéria (n-l)s2 24 * 0.12 alpha <- 0.05 m <- 1.99 s <- 0.1 c <- 0.0064 n <- 25 (tO <- (n-l)*s"2/c) ## [1] 37.5 (wl <- qchisq(alpha/2, n-1)) ## [1] 12.40115 (w2 <- qchisq(l-alpha/2, n-1)) ## [1] 39.36408 Testovací statistika íg nabývá hodnoty 37.5, kritický obor má tvar W = (-00 ; 12.4012) U (39.3640 ; oo) Protože tg ^ W, nejsme oprávněni na hladině významnosti a = 0.05 zamítnout tvrzení výrobce. b) Testování pomocí intervalu spolehlivosti (dh <- (n-l)*s~2/qchisq(l-alpha/2, n-1)) ## [1] 0.006096929 (hh <- (n-l)*s~2/qchisq(alpha/2, n-1)) ## [1] 0.01935304 95% empirický interval spolehlivosti pro a má tvar IS= (0.0781; 0.1391). Protože c = 0.08 G IS, nejsme oprávněni na hladině významnosti a = 0.05 zamítnout tvrzení výrobce. 51 c) Testování pomocí p-hodnoty (p.val <- 2*min(pchisq(t0, n-1), l-pchisq(tO, n-1))) ## [1] 0.0779636 Protože p-hodnota = 0.078 > 0.05, nejsme oprávněni na hladině významnosti a = 0.05 zamítnout tvrzení výrobce. Příklad 7.8. Interval spolehlivosti pro rozdíl parametrů /ii — /i2 dvourozměrného rozdělení: Bylo vylosováno 6 vrhů selat a z nich vždy dva sourozenci. Jeden z nich vždy dostal náhodně dietu č. 1 a druhý dietu č. 2. Přírůstky v dkg jsou následující: (62,52), (54,56), (55,49), (60,50), (53,51), (58,50). Za předpokladu, že uvedené dvojice tvoří náhodný výběr z dvourozměrného rozdělení s vektorem středních hodnot (/íi,/í2) a jejich rozdíly se řídí normálním rozdělením, sestrojte 95% interval spolehlivosti pro rozdíl středních hodnot. dl <- c(62, 54, 55, 60, 53, 58) d2 <- c(52, 56, 49, 50, 51, 50) x <- dl-d2 # rozdil středních hodnot přírůstku u diety 1 a diety 2 Shapirův - Wilkův test normality: shapiro.test(x)$p.val ## [1] 0.3241142 P-hodnota S-W testu je 0.3241 > 0.05, tedy nulovou hypotézu o normálním rozdělení náhodného výběru nezamítáme na hladině významnosti a = 0.05. Výpočet intervalů spolehlivosti: ## [1] "dh = 0.6265" ## [1] "hh = 10.7069" 95% interval spolehlivosti pro rozdíl středních hodnot fii — /i2 má tvar: (0.6265; 10.7069). 0.63 dkg < jj, < 10.71 dkg s pravděpodobností 0.95. Poznámka: K nalezení hranic 95% empirického intervalu spolehlivosti pro rozdíl středních hodnot fii — /i2 dvourozměrného rozdělení můžeme použít funkci t.test(x.y) s argumentem paired=T (párová test) a argumentem alter-native='two.sideď (oboustranná alternativa). x <- c(62, 54, 55, 60, 53, 58) y <- c(52, 56, 49, 50, 51, 50) t.test(x, y, paired=T, alternative='two.sideď)$conf.int ## [1] 0.6264613 10.7068720 ## attr(,"conf.level") ## [1] 0.95 Příklad 7.9. Testování hypotézy o rozdílu parametrů /ii —/i2 dvourozměrného rozdělení: Bylo vybráno šest nových vozů téže značky a po určité době bylo zjištěno, o kolik mm se sjely jejich levé a pravé přední pneumatiky. Výsledky: (1.8,1.5), (1.0,1.1), (2.2,2.0), (0.9,1.1), (1.5,1.4), (1.6,1.4). Za předpokladu, že uvedené dvojice tvoří náhodný výběr z dvourozměrného rozdělení s vektorem středních hodnot (pi, /i2) a jejich rozdíly se řídí normálním 52 rozdělením, testujte na hladině významnosti a = 0.05 hypotézu, že obě pneumatiky se sjíždí stejně rychle. Označme \i = \i\ — Hi- Na hladině významnosti a = 0.05 testujeme hypotézu Hq : ji = 0 proti oboustranné alternativě H\ : fj, ^ 0. Testování provedeme párovým í-testem. a) Testování pomocí kritického oboru Shapirův - Wilkův test normality: shapiro.test(x)$p.val ## [1] 0.4522054 P-hodnota S-W testu je 0.4522 > 0.05, tedy nulovou hypotézu o normálním rozdělení náhodného výběru nezamítáme na hladině významnosti a = 0.05. Testování nulové hypotézy: m <- mean(x) n <- length(x) s <- sd(x) c <- 0 (tO <- (m-c)/(s/sqrt(n))) ## [1] 1.051758 (wl <- qt(l-alpha/2, n-1)) ## [1] 2.570582 (w2 <- qt(alpha/2, n-1)) ## [1] -2.570582 Testovací statistika ío nabývá hodnoty 1.0518, kritický obor má tvar W = (-oo ; -2.5706) U (2.5706 ; co) Protože ío 4- W, Hq nezamítáme na hladině významnosti a = 0.05. b) Testování pomocí intervalu spolehlivosti (dh <- m-s/sqrt(n)*qt(l-alpha/2, n-1)) ## [1] -0.1203401 (hh <- m-s/sqrt(n)*qt(alpha/2, n-1)) ## [1] 0.2870068 95% empirický interval spolehlivosti pro a má tvar IS= (-0.1203; 0.2870). Protože c = 0 G IS, Hq nezamítáme na hladině významnosti a = 0.05. 53 c) Testování pomocí p-hodnoty (p.val <- 2*min(pt(t0, n-1), l-pt(tO, n-1))) ## [1] 0.341062 Protože p-hodnota = 0.3411 > 0.05, Hq nezamítáme na hladině významnosti a = 0.05. Poznámka: K otestování nulové hypotézy o rozdílu parametrů /ii — /j.2 dvourozměrného rozdělení můžeme použít funkci t.test(x.y) s argumentem paired=T (párový test) a argumentem alternative='two.sideď (oboustranná alternativa). x <- c(1.8, 1.0, 2.2, 0.9, 1.5, 1.6) y <- c(1.5, 1.1, 2.0, 1.1, 1.4, 1.4) t.test(x, y, paired=T, alternative='two.sided') ## ## Paired t-test ## ## data: x and y ## t = 1.0518, df = 5, p-value = 0.3411 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -0.1203401 0.2870068 ## sample estimates: ## mean of the differences ## 0.08333333 54 8 Parametrické úlohy o dvou nezávislých náhodných výběrech z normálních rozdělení a jednom náhodném výběru z alternativního rozdělení Parametrické úlohy o dvou nezávislých náhodných výběrech z normálních rozdělení Příklad 8.1. Interval spolehlivosti pro parametrickou funkci yui — /j-2' Bylo vylosováno 11 stejně starých selat téhož plemene. Šesti z nich byla předepsána výkrmná dieta č. 1 a zbylým pěti výkrmná dieta č. 2. Průměrné denní přírůstky v dkg za dobu půl roku jsou následující: Zjištěné hodnoty považujeme za realizace dvou nezávislých náhodných výběrů pocházejících z rozdělení N(fii,a2) a N(p2, o"2)- Sestrojte 95% empirický interval spolehlivosti pro rozdíl středních hodnot yui — /i2- x <- c(62, 54, 55, 60, 53, 58) y <- c(52, 56, 49, 50, 51) ml <- mean(x) m2 <- mean(y) si <- sd(x) s2 <- sd(y) nl <- length(x) n2 <- length(y) alpha <- 0.05 sh <- sqrt(((nl-l)*sl~2 + (n2-l)*s2~2)/(nl+n2-2)) (dh <- ml-m2-sh*sqrt(l/nl+l/n2)*qt(l-alpha/2, nl+n2-2)) ## [1] 0.9919634 (hh <- ml-m2-sh*sqrt(l/nl+l/n2)*qt(alpha/2, nl+n2-2)) ## [1] 9.808037 S pravděpodobností alespoň 0.95 platí, že 0.99 dkg < fii — /J.2 < 9.81 dkg. Příklad 8.2. Testování hypotéz o parametrických funkcích /ii — /i2, o"i/o"2: 1. Pro datový soubor z příkladu 8.1 testujte na hladině významnosti a = 0.05 hypotézu, že a) rozptyly hmotnostních přírůstků selat při obou výkrmných dietách jsou shodné; b) obě výkrmné diety mají stejný vliv na hmotnostní přírůstky selat. 2. Výsledek testování podpořte krabicovým diagramem. Shapirův - Wilkův test normality Nejprve je potřeba otestovat normalitu obou náhodných výběrů. x <- c(62, 54, 55, 60, 53, 58) y <- c(52, 56, 49, 50, 51) shapiro.test(x)$p.value ## [1] 0.6194994 dieta č. 1 dieta č. 2 62 54 55 60 53 58 52 56 49 50 51 IS= (0.9920; 9.8080) 55 shapiro.test(y)$p.value ## [1] 0.4271986 P-hodnota S-W testu pro přírůstky selat krmených dietou č. 1 je 0.6195 > 0.05, p-hodnota S-W testu pro přírůstky selat krmených dietou č. 2 je 0.4272 > 0.05. V obou případech tedy nulovou hypotézu o normalitě dat nezamítáme na hladině významnosti a = 0.05. ad a) Testování hypotézy o shodě rozptylů. 2 Na hladině významnosti a = 0.05 testujeme nulovou hypotézu Hq : ^ = 1 oproti alternativní hypotéze CT2 cr2 Hi : —h ^ 1. K testování použijeme dvouvýběrový ŕ1-test. CT2 Testování nulové hypotézy: i. Testování pomocí kritického oboru ml <- mean(x) m2 <- mean(y) sl <- sd(x) s2 <- sd(y) ni <- length(x) n2 <- length(y) alpha <- 0.05 sh <- sqrt(((nl-l)*sl~2 + (n2-l)*s2~2)/(nl+n2-2)) (tO <- sl~2/s2~2) ## [1] 1.753425 (wl <- qf (alpha/2, nl-1, n2-l)) ## [1] 0.1353567 (w2 <- qf(l-alpha/2, nl-1, n2-l)) ## [1] 9.364471 Testovací statistika ío nabývá hodnoty 1.7534, kritický obor má tvar W = (-00 ; 0.1354) U (9.3645 ; oo) Protože ío ^ W, Hq o shodě rozptylů a'f a o\ nezamítáme na hladině významnosti a = 0.05. ii. Testování pomocí intervalu spolehlivosti (dh <- (sl~2/s2~2)/(qf(l-alpha/2, nl-1, n2-l))) ## [1] 0.1872423 (hh <- (sl~2/s2~2)/(qf(alpha/2, nl-1, n2-l))) ## [1] 12.9541 95% empirický interval spolehlivosti pro podíl o\ja\ má tvar IS = (0.1872; 12.9541). Protože c = 1 € IS, Hq o shodě rozptylů a'f a trf nezamítáme na hladině významnosti a = 0.05. 56 iii. Testování pomocí p-hodnoty (p.val <- 2*min(pf(tO, nl-1, n2-l), l-pf(tO, nl-1, n2-l))) ## [1] 0.6063451 Protože p-hodnota = 0.6063 > 0.05, Hq o shodě rozptylů a\ a a\ nezamítáme na hladině významnosti a = 0.05. Upozornění: V případě zamítnutí hypotézy o shodě rozptylů je zapotřebí použít test se samostatnými odhady rozptylu. ad b) Testování hypotézy o shodě středních hodnot Na hladině významnosti a = 0.05 testujeme nulovou hypotézu Hq : pi — p2 = 0 oproti alternativní hypotéze H\ : pi — p2 0. K testování použijeme dvouvýběrový test o rozdílu středních hodnot. Testování nulové hypotézy: i. Testování pomocí kritického oboru c <- 0 (tO <- ((ml-m2)-c)/(sri*sqrt(l/nl+l/n2))) ## [1] 2.771222 (wl <- qt (alpha/2, nl+n2-2)) ## [1] -2.262157 (w2 <- qt(l-alpha/2, nl+n2-2)) ## [1] 2.262157 Testovací statistika tg nabývá hodnoty 2.7712, kritický obor má tvar W = (-oo ; -2.2622) U (2.2622 ; oo) Protože tg G W, Hq o shodě středních hodnot pi a pg zamítáme na hladině významnosti a = 0.05. ii. Testování pomocí intervalu spolehlivosti V příkladu 8.1 jsme zjistili, že 95% oboustranný interval spolehlivosti pro rozdíl středních hodnot pi — pg má tvar IS = (0.9920; 9.8080). Protože c = 0 ^ IS, Hq o shodě středních hodnot pi a p2 zamítáme na hladině významnosti a = 0.05. iii. Testování pomocí p-hodnoty (p.val <- 2*min(pt(t0, nl+n2-2), l-pt(tO, nl+n2-2))) ## [1] 0.02171008 Protože p-hodnota = 0.0217 < 0.05, Hq o shodě středních hodnot pi a p2 zamítáme na hladině významnosti a = 0.05. Znamená to, že s rizikem omylu nejvýše 5 % se prokázalo, že obě výkrmné diety se liší účinností. Poznámka: K otestování nulové hypotézy o rozdílu středních hodnot pi — p2 dvou nezávislých náhodných výběrů z normálních rozdělení můžeme použít funkci t.test(x.y) s argumentem alternative='two.sideď (oboustranná alternativa) a argumentem var.equal=T (rozptyly obou náhodných výběrů si jsou rovné). 57 x <- c(62, 54, 55, 60, 53, 58) y <- c(52, 56, 49, 50, 51) t.test(x, y, alternative='two.sided', var.equal=T) ## ## Two Sample t-test ## ## data: x and y ## t = 2.7712, df = 9, p-value = 0.02171 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.9919634 9.8080366 ## sample estimates: ## mean of x mean of y ## 57.0 51.6 Upozornení: Pokud bychom na hladině významnosti a = 0.05 zamítli nulovou hypotézu o shodě rozptylů a\ a oj, mohli bychom k otestování nulové hypotézy o shodě středních hodnot \x\ a \xi použít opět funkci t.test s argumentem alternative='two.sideď (oboustranná alternativa) a argumentem var.equal=F. Tento argument modifikuje klasický í-test na í-test s Welschovou aproximací stupňů volnosti, která se používá v případě, že rozptyly obou náhodných výběrů nejsou shodné. Krabicový diagram Řešení příkladu doplníme ještě o krabicový diagram. boxplot(x, y, names=c('dieta cl', 'dieta c.2'), ylab='hmotnostni prirustky selat', main='Hmotnostni prirustky selat', col='blanchedalmond', border='coral4') mtext('Krabicový graf', line=0.4, cex=0.9) CD Hmotnostni prirustky selat Krabicový graf dieta c.1 dieta c.2 58 Příklad k samostatnému řešení Příklad 8.3. Načtěte datový soubor vyska.txt, který obsahuje údaje o výšce 48 studentek VŠE v Praze (proměnná vyska) a obor jejich studia (1 - národní hospodářství, 2 - informatika). a) Pomocí S-W testu ověřte na hladině významnosti a = 0.1 předpoklad o normalitě výšek v obou skupinách studentek. b) Na hladině významnosti a = 0.1 testujte hypotézu o shodě rozptylů výšek studentek v daných dvou oborech studia. c) Na hladině významnosti a = 0.1 testujte hypotézu o shodě středních hodnot výšek studentek v daných dvou oborech studia. d) Výpočet doplňte krabicovými diagramy. Shapirův - Wilkův test normality data <- read.tableCvyska.txt', header=T) x <- data$vyska[data$obor==l] y <- data$vyska[data$obor==2] shapiro.test(x)$p.value ## [1] 0.6067928 shapiro.test(y)$p.value ## [1] 0.1119235 P-hodnota S-W testu pro výšku studentek oboru národní hospodářství je 0.6068 > 0.1, p-hodnota S-W testu pro výšku studentek oboru informatika je 0.1119 > 0.1. V obou případech tedy nulovou hypotézu o normalitě dat nezamítáme na hladině významnosti a = 0.1. Testování hypotézy o shodě rozptylů i. Testování pomocí kritického oboru ## [1] "tO = 1.9873" ## [1] "wl = 0.5033" ## [1] "w2 = 2.0905" Protože ío ^ W, Ho o shodě rozptylů a\ a trf nezamítáme na hladině významnosti a = 0.1. ii. Testování pomocí intervalu spolehlivosti ## [1] "dh = 0.9506" ## [1] "hh = 3.9487" Protože c = 1 G IS, Ho o shodě rozptylů cr'( a a\ nezamítáme na hladině významnosti a = 0.1. iii. Testování pomocí p-hodnoty ## [1] "p.val = 0.1249" Protože p-hodnota = 0.1249 > 0.05, Hq o shodě rozptylů a'f a a\ nezamítáme na hladině významnosti a = 0.1. 59 Testování hypotézy o shodě středních hodnot i. Testování pomocí kritického oboru ## [1] "tO = 1.744" ## [1] "wl = -1.6787" ## [1] "w2 = 1.6787" Protože tg G W, Hq o shodě středních hodnot /ii a /j.2 zamítáme na hladině významnosti a = 0.1. ii. Testování pomocí intervalu spolehlivosti ## [1] "dh = 0.1095" ## [1] "hh = 5.7334" Protože c = 0 ^ IS, Hq o shodě středních hodnot /ii a /j.2 zamítáme na hladině významnosti a = 0.1. iii. Testování pomocí p-hodnoty ## [1] "p.val = 0.0878" Protože p-hodnota = 0.0878 < 0.1, Hq o shodě středních hodnot /ii a /j.2 zamítáme na hladině významnosti a = 0.1. Krabicový diagram boxplot(x, y, names=c('narodni hospodářství', 'informatika'), ylab='vyska studentek', xlab='obor', main='Vyska studentek VSE', col='blanchedalmond', border='coral4') mtext('Krabicový graf', line=0.4, cex=0.9) Vyska studentek VSE Krabicový graf narodni hospodářství obor informatika 60 8.1 Parametrické úlohy o jednom náhodném výběru z alternativního rozdělení Příklad 8.4. Asymptotický interval spolehlivosti pro parametr 9 alternativního rozdělení: Může politická strana, pro niž se v předvolebním průzkumu vyslovilo 60 z 1000 dotázaných osob, očekávat se spolehlivostí 0.95, že by v této době ve volbách překročila 5 % hranici pro vstup do parlamentu? Zavedeme náhodné veličiny X\,..., XiooOj přičemž Xi = 1, když se i-tá osoba vysloví pro danou politickou stranu, a Xi = 0 jinak; i = 1,..., 1000. Tyto náhodné veličiny tvoří náhodný výběr z rozdělení A(9). V tomto případě n = 1000, m = 60/1000 = 0.06, a = 0.05, = m0.95 = 1.645. Ověření podmínky n9(l — 9) > 9: parametr 9 neznáme, musíme ho tedy nahradit výběrovým průměrem. Pak 1000 * 0.06 * 0.94 = 56.4 > 9. 95% levostranný interval spolehlivosti pro parametr 9 má potom tvar S pravděpodobností přibližně 0.95 je tedy 9 > 0.0476. Protože tento interval zahrnuje i hodnoty nižší než 0.05, nelze vyloučit, že strana získá méně než 5% hlasů. x <- 60 n <- 1000 m <- x/n alpha <- 0.05 (dh <- m-sqrt(m*(1-m)/n)*qnorm(l-alpha)) ## [1] 0.04764716 Příklad k samostatnému řešení Příklad 8.5. Přírůstky cen akcií na burze (v %) u 10 náhodně vybraných společností dosáhly těchto hodnot: 10, 16, 5, 10, 12, 8, 4, 6, 5, 4. Sestrojte 95% asymptotický empirický interval spolehlivosti pro pravděpodobnost, že přírůstek ceny akcie překročí 8.5%. ## [1] "dh = 0.0964" ## [1] "hh = 0.7036" 0.096 < 9 < 0.704 s pravděpodobností aspoň 0.95. Znamená to, že pravděpodobnost, že přírůstek ceny akcie překročí 8.5%, je alespoň 9.6% a nanejvýš 70.4% (při spolehlivosti 95%.) Příklad 8.6. Testování hypotézy o parametru 9 alternativního rozdělení: Určitá cestovní kancelář organizuje zahraniční zájezdy podle individuálních přání zákazníků. Z několika minulých let ví, že 30 % všech takto organizovaných zájezdů má za cíl zemi X. Po zhoršení politických podmínek v této zemi se cestovní kancelář obává, že se zájem o tuto zemi mezi zákazníky sníží. Ze 150 náhodně vybraných zákazníků v tomto roce má 38 za cíl právě zemi X. Potvrzují nejnovější data pokles zájmu o tuto zemi? Volte hladinu významnosti a = 0.05. Máme náhodný výběr Xi,...,Xi5o z rozdělení .4(0.3). Testujeme Hq : 9 = 0.3 proti levostranné alternativě H\ : 9 < 0.3. V tomto případě je testovacím kritériem statistika V našem případě M - c 61 která se za platnosti nulové hypotézy asymptoticky řídí standardizovaným normálním rozdělením ÍV(0,1). Nejprve musíme ověřit splnění podmínky n9(l — 9) > 9: 150 * 0.3 * 0.7 = 31.5 > 9. a) Testování pomocí kritického oboru Vypočteme realizaci testovacího kritéria: x <- 38 n <- 150 c <- 0.3 m <- x/n alpha <- 0.05 (tO <- (m-c)/sqrt((c*(l-c)/n))) ## [1] -1.247219 (wl <- qnorm(alpha)) ## [1] -1.644854 Protože testovací kritérium nepatří do kritického oboru, Hq nezamítáme na asymptotické hladině významnosti a = 0.05. b) Testování pomocí intervalu spolehlivosti Proti levostranné alternativě H\ postavíme 95% pravostranný interval spolehlivosti. = -1.2472. Kritický obor má tvar: W = (-oo ; ua) U (-00 ; -1.645). (—oo; hh) kde realizace horní hranice = 0.3117 (hh <- m-sqrt(m*(1-m)/n)*qnorm(alpha)) ## [1] 0.3117439 Protože 0.3 G IS = (—oo; 0.3117), Hq nezamítáme na asymptotické hladině významnosti a = 0.05. 62 c) Testování pomocí p-hodnoty V části a) jsme spočítali hodnotu testovací statistiky ío = —1-2472. Protože máme levostrannou alternativní hypotézu Hi, vypočítáme p-hodnotu podle vzorce p-hodnota = P(T0 < t0) = P(T0 < t0) = pnorm(t0) = 0.1062 (p.val <- pnorm(tO)) ## [1] 0.1061586 Protože p-hodnota = 0.1062 > 0.05, H0 nezamítáme na asymptotické hladině významnosti a = 0.05. 63 9 Analýza rozptylu jednoduchého třídění Příklad 9.1. V jisté továrně se měřil čas, který potřeboval každý ze tří dělníků k uskutečnění téhož pracovního úkonu. Cas v minutách: 1.dělník: 3.6 3.8 3.7 3.5 2.dělník: 4.3 3.9 4.2 3.9 4.4 4.7 3.dělník: 4.2 4.5 4.0 4.1 4.5 4.4 Na hladině významnosti a = 0.05 testujte hypotézu, že výkony těchto tří dělníků jsou stejné. Zamítnete-li nulovou hypotézu, určete, výkony kterých dělníků se liší na dané hladině významnosti a = 0.05. Průzkumová analýza Úloha vede na analýzu rozptylu jednoduchého třídění. Načteme datový soubor cas_delniku.txt. Proměnná X obsahuje zjištěné časy, proměnná ID nabývá hodnoty 1 pro 1. dělníka, hodnoty 2 pro 2. dělníka a hodnoty 3 pro 3. dělníka. data <- read.delim('cas_delniku.txt', sep='', dec=',', header=T) názvy <- c('delnik ľ, 'delnik 2' , 'delnik 3') X <- data$X ID <- data$ID prum <- srn <- N <- NULL for(i in 1:3){ prum[i] <- mean(X[ID==i]) sm[i] <- sd(X[ID==i]) N [i] <- length(X[ID==i]) } prum[4] <-mean(X) srn[4] <- sd(X) N[4] <- length(X) tab <- data.frame(prumer=prum, N=N, sm.odch=sm, row.names=c(názvy, 'všichni')) round(tab,4) ## prumer N sm.odch ## delnik 1 3.6500 4 0.1291 ## delnik 2 4.2333 6 0.3077 ## delnik 3 4.2833 6 0.2137 ## všichni 4.1063 16 0.3530 Komentář: Na uskutečnění daného pracovního úkonu potřebuje nejkratší čas 1. dělník. Podává také nej vyrovnanější výkony - směrodatná odchylka proměnné X je u něj nejmenší. Naopak nejpomalejší je 3. dělník. 64 Krabicový graf boxplot(X[ID==1], X[ID==2], X[ID==3], names=nazvy, ylab='cas (v min)', xlab='delnik', main='Cas potřebný k provedeni pracovniho úkonu', col='blanchedalmond', border='coral4') mtext('Krabicový graf', line=0.4, cex=0.9) points(c(mean(X[ID==l]), mean(X[ID==2]), mean(X[ID==3])), pch=15, col='red') Cas potrebný k provedeni pracovniho úkonu Krabicový graf oj o 00 oó co OÓ delnik 1 delnik 3 Testování normality Na testování normality všech tří výběrů použijeme kvůli jejich malým rozsahům S-W test. ## [1] ## [1] ## [1] "S-W test, delnik 1 "S-W test, delnik 2 "S-W test, delnik 3 0.9719' 0.5819' 0.3313' Q-Q graf Cas pracovniho úkonu — delnik 1 ~i-1-r -0.5 0.0 0.5 teoreticky kvantil Q-Q graf Cas pracovniho úkonu — delnik 2 ~i-1-1-r -1.0 -0.5 0.0 0.5 1.0 teoreticky kvantil Q-Q graf Cas pracovniho úkonu — delnik 3 ~i-1-1-1-r -1.0 -0.5 0.0 0.5 1.0 teoreticky kvantil Komentář: Protože ve všech třech případech je p-hodnota S-W testu > 0.05, nulovou hypotézu o normalitě časů všech tří dělníků nezamítáme na hladině významnosti a = 0.05. Z Q-Q grafů dále vidíme, že tečky se ve všech třech případech jen málo odchylují od přímky, což podporuje náš závěr, že všechny tři výběry pochází z normálního rozdělení. 65 Test homogenity rozptylů Jelikož náhodné výběry pochází z normálního rozdělení, je vhodné na testování hypotézy o shodě rozptylů všech tří výběrů spoužít Levenův test. Na hladině významnosti a = 0.05 testujeme nulovou hypotézu Hq : o\ = o\ = a\ oproti alternativní hypotéze H\ : a2 ^ a2 pro alespoň jednu dvojici i, j. K otestování použijeme funkci levene.test z knihovny lawstat s argumentem location='mean', čímž zvolíme klasickou formu Levenova testu. Pokud bychom zadali funkci levene.test s argumentem location='median', získali bychom robustnější modifikaci Levenova testu, která ve svém výpočtu nahrazuje aritmetický průměr mediánem. Balíček lawstat disponuje také příkazem na výpočet Bartlettova testu, který je k dispozici ve funkci bartlett.test. library(lawstat) levene.test(X, ID, location='mean') ## ## classical Levene's test based on the absolute deviations from the mean ( none ## not applied because the location is not set to median ) ## ## data: X ## Test Statistic = 1.5142, p-value = 0.2564 #levene.test(X, ID, location=1 median 1) #bartlett.test(X, ID) Komentár: Testovací statistika Levenova testu nabývá hodnoty 1.5142, odpovídající p-hodnota = 0.256, tedy na hladině významnosti a = 0.05 nezamítáme hypotézu o shodě rozptylů a'f, cr| a cr|. Test o shodě středních hodnot: Na hladině významnosti a = 0.05 testujeme nyní nulovou hypotézu Ho ■ Mi = M2 = M3 oproti alternativní hypotéze H\ : /j-í fij pro alespoň jednu dvojici i,j. r <- 3 n <- length(X) Xi. <- Mi. <- ni <- NULL for(i in l:r){ Xi. [i] <- sum(X[ID==i]) ni[i] <- length(X[ID==i]) Mi. [i] <- sum(X[ID==i])/ni[i] } X. . <- sum(X) M. . <- mean(X) SA <- sum(ni*(Mi.-M. .)~2) fA <- r-1 ST <- sum((X-M..)~2) SE <- ST-SA fE <- n-r Fa <- (SA/fA)/(SE/fE) pl <- pf(Fa,r-1,n-r) 66 p2 <- 1-pl p.val <- min(pl, p2) (tab <- round(data.frame(SA=SA, fA=fA, SE=SE, fE=fE, ST=ST, fT=n, Fa=Fa, p.val=p.val), digits=5)) ## SA fA SE fE ST fT Fa p.val ## 1 1.11771 2 0.75167 13 1.86938 16 9.66533 0.00268 Skupinový součet čtverců Sa = 1.1177, počet stupňů volnosti f a = 2, reziduálni součet čtverců Se = 0.7517, počet stupňů volnosti f e = 13, testovací statistika p = Sa/fa Se/f e nabývá hodnoty 9.6653, počet stupňů volnosti čitatele = 2, jmenovatele = 13. p-hodnota = 0.00268, tedy na hladině významnosti a = 0.05 zamítáme nulovou hypotézu o shodě středních hodnot. Metoda mnohonásobného porovnávání Jelikož jsme na hladině významnosti 0.05 zamítli nulovou hypotézu o shodě středních hodnot, chceme nyní zjistit, které dvojice středních hodnot se od sebe významně liší. Stanovíme nulové a alternativní hypotézy pro dvojice středních hodnot H01 : Mi = M2 oproti Hn : >ui ^ fi2; H02 : Mi = M3 oproti H12 : Hi ^ ^3 5 Hoz ■ M2 = M 3 oproti H13 : fi2 M3- Protože v každé skupině máme různý počet pozorování, je vhodné použít na mnohonásobné porovnávání Scheffého metodu. Pravou a levou stranu Scheffého metody získáme použitím funkce Scheffe(X, group, names, alpha), která je k dispozici v RSkriptu AS-funkce.R. Argumenty funkce jsou X ... vektor hodnot, group ... vektor přiřazující každému pozorování skupinu, do níž náleží, names ... názvy jednotlivých skupin, alpha ... hladina významnosti. source('AS-funkce.R') Scheffe(X, ID, názvy, alpha=0.05) ## $R ## delnik 1 delnik 2 delnik 3 ## delnik 1 0.4690839 0.4282131 0.4282131 ## delnik 2 0.4282131 0.3830054 0.3830054 ## delnik 3 0.4282131 0.3830054 0.3830054 ## ## $L ## delnik 1 delnik 2 delnik 3 ## delnik 1 0.0000000 0.5833333 0.6333333 ## delnik 2 0.5833333 0.0000000 0.0500000 ## delnik 3 0.6333333 0.0500000 0.0000000 l*(Scheffe(X, ID, názvy, alpha=0.05)$R >= Scheffe(X, ID, názvy, alpha=0.05)$L) ## delnik 1 delnik 2 delnik 3 ## delnik 110 0 ## delnik 2 0 11 ## delnik 3 0 11 Komentář: Porovnáním pravé a levé strany Scheffého metody vidíme, že na hladině významnosti a = 0.05 zamítáme nulovou hypotézu o shodě středních hodnot yui a /i2 a středních hodnot yui a jj,3. Výsledek Scheffého metody tedy ukazuje, že na hladině významnosti a = 0.05 se liší výkony dělníků (1,2), (1,3), naopak výkony dělníků (2,3) se neliší. 67 Příklad 9.2. Na střední škole byl uskutečněn experminet zjišťující efektvitu jednotlivých pedagogických metod. Studenti byli rozděleni do pěti skupin a každá skupina byla vyučována pomocí jedné z pedagogických metod: tradiční způsob, programová výuka, audiotechnika, audiovizuální technika a vizuální technika. Z každé skupiny byl potom vybrán náhodný vzorek studentů a všichni byli podrobeni témuž písemnému testu. Výsledky testu jsou uvedeny v následující tabulce a v souboru pet_metod.txt: metoda počet bodů tradicni programová audio audiovizuálni vizuálni 76.2 48.3 85.1 63.7 91.6 87.2 85.2 74.3 76.5 80.3 67.4 67.9 72.1 60.4 67.3 60.1 55.4 72.3 40.0 75.8 81.6 90.3 78.0 67.8 57.6 50.5 70.2 88.8 67.1 77.7 73.9 Na hladině významnosti a = 0.05 testujte hypotézu, že znalosti všech studentů jsou stejné a nezávisí na použité pedagogické metodě. V případě zamítnutí hypotézy zjistěte, které výběry se liší na hladině významnosti 0.05. Průzkumová analýza Načteme datový soubor pet_metod.txt. Proměnná body obsahuje dosažené počty bodů a proměnná metoda označení příslušné pedagogické metody. Nejprve vypočítáme průměry, směrodatné odchylky a rozsahy všech tří výběrů: data <- read.delimCpet_metod.txt', sep='', dec=',', header=T) názvy <- c('tradicni', 'programová', 'audiotechnika', 'audiovizuálni', 'vizuálni') X <- data$body ID <- data$metoda r <- length(unique(ID)) prum <- srn <- N <- NULL for(i in l:r){ prum[i] <- mean(X[ID==i]) sm[i] <- sd(X[ID==i]) N [i] <- length (X [ID==i]) } prum[r+l] <-mean(X) sm[r+l] <- sd(X) N[r+1] <- length(X) tab <- data.frame(prumer=prum, N=N, sm.odch=sm, row.names=c(názvy, 'všechny')) round(tab,4) ## tradicni 75.3500 ## programová 73.0125 ## audiotechnika 59.0200 ## audiovizuálni 75.1833 ## vizuálni 71.3667 prumer N srn.odch 16.5390 7.8650 12.4594 11.3286 12.6920 ## všechny 71.3097 31 12.6953 Komentář: Nejlepších výsledků dosahují studenti vyučovaní tradiční metodou, podávají však nejméně vyrovnané výkony (počty bodů v této skupině mají největší směrodatnou odchylku). Naopak nejhoršího výsledku dosáhli studenti vyučovaní audio metodou. Nejvyrovnanější výkony pozorujeme u studentů vyučovaných programovou metodou. 68 Krabicový graf boxplot(X[ID==l], X[ID==2], X[ID==3], X[ID==4], X[ID==5], names=nazvy, ylab='body' xlab='metóda', main='Výukové metódy', col='khakil', border='orange4') mtext('Krabicový graf', line=0.4, cex=0.9) points(c(mean(X[ID==l]), mean(X[ID==2]), mean(X[ID==3]), mean(X[ID==4]), mean(X[ID Výukové metody Krabicový graf tradicni t-1-r programová audiotechnika audiovizuálni vizuálni metoda Testování normality Na testování normality všech tří výběrů použijeme z důvodu jejich malých rozsahů S-W test. print(paste('S-W test, metoda 1:', round(shapiro.test(X[ID==1])$p.value,4))) ## [1] "S-W test, metoda 1: 0.4177" print(paste('S-W test, metoda 2:', round(shapiro.test(X[ID==2])$p.value,4))) ## [1] "S-W test, metoda 2: 0.9966" print(paste('S-W test, metoda 3:', round(shapiro.test(X[ID==3])$p.value,4))) ## [1] "S-W test, metoda 3: 0.7663" print(paste('S-W test, metoda 4:', round(shapiro.test(X[ID==4])$p.value,4))) ## [1] "S-W test, metoda 4: 0.9577" print(paste('S-W test, metoda 5:', round(shapiro.test(X[ID==5])$p.value,4))) ## [1] "S-W test, metoda 5: 0.8814" for(i in l:r){ qqnorm(X[ID==i], col='darkblue', xlab='teoreticky kvantil', ylab='pozorovaný kvantil', main='Q-Q graf') qqline(X[ID==i] , col='reď) mtext(paste('Cas pracovniho úkonu — metoda', i), line=0.4, cex=0.8)} 69 Q-Q graf Q-Q graf Q-Q graf Cas pracovního úkonu — metoda 1 Cas pracovního úkonu — metoda 2 Cas pracovního úkonu — metoda 3 -1.0 -0.5 0.0 0.5 1.0 -1.5 -0.5 0.0 0.5 1.0 1.5 -1.0 -0.5 0.0 0.5 1.0 teoreticky kvantil teoreticky kvantil teoreticky kvantil Q-Q graf Q-Q graf Cas pracovního úkonu — metoda 4 Cas pracovního úkonu — metoda 5 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 teoreticky kvantil teoreticky kvantil Komentář: Protože ve všech pěti případech je p-hodnota S-W testu > 0.05, nulovou hypotézu o normalitě dat u každé metody nezamítáme na hladině významnosti a = 0.05. Vzhled Q-Q grafů potvrzuje náš závěr, že předpoklad normality je ve všech pěti případech oprávněný. Test homogenity rozptylů K testování hypotézy o shodě rozptylů opět použijeme klasický Levenův test. Na hladině významnosti a = 0.05 testujeme nulovou hypotézu Hq : a\ = a\ = cr| = a\ = a2 oproti alternativní hypotéze H\ : a2 ^ a2 pro alespoň jednu dvojici i, j. library(lawstat) levene.test(X, ID, location='mean') ## ## classical Levene's test based on the absolute deviations from the mean ( none ## not applied because the location is not set to median ) ## ## data: X ## Test Statistic = 0.81903, p-value = 0.5248 #levene.test(X, ID, location=1 median 1) #bartlett.test(X, ID) Komentár: Testovací statistika F se realizuje hodnotou 0.8190, odpovídající p-hodnota = 0.5248. Na hladině významnosti a = 0.05 tedy nezamítáme hypotézu o shodě rozptylů. Test o shodě středních hodnot: Na hladině významnosti a = 0.05 testujeme nyní nulovou hypotézu 70 Ho ■ Mi = M2 = M3 = IH = M5 oproti alternativní hypotéze Hi : /j-í nj pro alespoň jednu dvojici i,j. n <- length(X) Xi. <- Mi. <- ni <- NULL for(i in l:r){ Xi. [i] <- sum(X[ID==i]) ni [i] <- length(X[ID==i]) Mi. [i] <- sum(X[ID==i])/ni[i] } X. . <- sum(X) M. . <- mean(X) SA <- sum(ni*(Mi.-M. .)~2) fA <- r-1 ST <- sum((X-M.. )~2) SE <- ST-SA fE <- n-r Fa <- (SA/fA)/(SE/fE) pl <- pf(Fa,r-1,n-r) p2 <- 1-pl p.val <- min(pl, p2) (tab <- round(data.frame(SA=SA, fA=fA, SE=SE, fE=fE, ST=ST, fT=n, Fa=Fa, p.val=p.val), digits=5)) ## SA fA SE fE ST fT Fa p.val ## 1 966.3737 4 3868.773 26 4835.147 31 1.62362 0.19825 Komentář: Testovací statistika F se realizuje hodnotou 1.6236, počet stupňů volnosti čitatele = 4, jmenovatele = 26, odpovídající p-hodnota = 0.1983, na hladině významnosti a = 0.05 tedy nezamítáme hypotézu o shodě středních hodnot. Znamená to, že s rizikem omylu nejvýše 5 % se neprokázal rozdíl v účinnosti jednotlivých pedagogických metod. Příklad 9.3. Pan Novák může cestovat z místa bydliště do místa pracoviště třemi různými způsoby: tramvají, autobusem a metrem s následným přestupem na tramvaj. Máme k dispozici jeho naměřené časy cestování do práce v době ranní špičky (včetně čekání na příslušný spoj) v minutách: autobus: 32 39 42 37 34 38 tramvaj: 30 34 28 26 32 metro: 40 37 31 39 38 33 34 Pro všechny tři způsoby dopravy vypočtěte průměrné časy cestování. Na hladině významnosti a = 0.05 testujte hypotézu, že doba cestování do práce nezávisí na způsobu dopravy. V případě zamítnutí nulové hypotézy zjistěte, které způsoby dopravy do práce se od sebe liší na hladině významnosti a = 0.05. Průzkumová analýza Načteme datový soubor doby_cestovani.txt. Proměnná cas obsahuje zjištěné doby cestování a proměnná ID označení příslušného způsobu dopravy. Nejprve vypočteme průměry, směrodatné odchylky a rozsahy všech tří výběrů: 71 data <- read.delimCdoby_cestovani.txt', sep='', dec=',', header=T) názvy <- c('tramvaj', 'autobus', 'metro') X <- data$cas ID <- data$ID r <- length(unique(ID)) prum <- sm <- N <- NULL for(i in l:r){ prum[i] <- mean(X[ID==i]) sm[i] <- sd(X[ID==i]) N[i] <- length(X[ID==i]) } prum[r+l] <-mean(X) sm[r+l] <- sd(X) N[r+1] <- length(X) tab <- data.frame(prumer=prum, N=N, sm.odch=sm, row.names=c(nazvy, 'všechny')) round(tab,4) ## prumer N sm.odch ## tramvaj 37.0000 6 3.5777 ## autobus 30.0000 5 3.1623 ## metro 36.0000 7 3.3665 ## všechny 34.6667 18 4.3791 Komentář: Nejkratší průměrnou dobu do zaměstnání pan Novák cestuje, když použije autobus, naopak nejdéle cestuje tramvají. Variabilita dob jednotlivých způsobů cestování je vcelku vyrovnaná. 72 Krabicový graf boxplot(X[ID==1], X[ID==2], X[ID==3], names=nazvy, ylab='doba cestováni (v min)', xlab='způsob dopravy', main='Doprava do prace', col='ali ceblue',border='dodgerblue4') mtext('Krabicový graf', line=0.4, cex=0.9) points(c(mean(X[ID==l]), mean(X[ID==2]), mean(X[ID==3])), pch=15, col='blue3') lo .a o "O Doprava do prace Krabicový graf tramvaj autobus způsob dopravy metro Testování normality Na testování normality všech tří výběrů použijeme z důvodu jejich malých rozsahů S-W test. print(paste('S-W test, autobus:', round(shapiro.test(X[ID==1])$p.value,4))) ## [1] "S-W test, autobus: 0.9539" print(paste('S-W test, tramvaj:', round(shapiro.test(X[ID==2])$p.value,4))) ## [1] "S-W test, tramvaj: 0.9672" print(paste('S-W test, metro: ', round(shapiro.test(X[ID==3])$p.value,4))) ## [1] "S-W test, metro: 0.6294" for(i in l:r){ qqnorm(X[ID==i], col='darkblue', xlab='teoreticky kvantil', ylab='pozorovaný kvantil', main='Q-Q graf') qqline(X[ID==i], col='red') mtext(paste('Cas pracovniho úkonu — metoda', i), line=0.4, cex=0.8)} 73 Q-Q graf Q-Q graf Q-Q graf Cas pracovního úkonu — metoda 1 Cas pracovního úkonu — metoda 2 Cas pracovního úkonu — metoda 3 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 teoreticky kvantil teoreticky kvantil teoreticky kvantil Komentář: Protože ve všech třech případech je p-hodnota S-W testu > 0.05, nulovou hypotézu o normalitě dat u každé metody nezamítáme na hladině významnosti a = 0.05. Vzhled Q-Q grafů potvrzuje náš závěr, že předpoklad normality je ve všech třech uvedených případech oprávněný. Test homogenity rozptylů Nulovou hypotézu o shodě rozptylů Hq : a'f = = Scheffe(X, ID, názvy, alpha=0.05)$L) ## tramvaj autobus metro ## tramvaj 1 0 1 ## autobus 0 10 ## metro 1 0 1 #Scheffe(X, ID, alpha=0.05, názvy) Komentáře: Z tabulky vyplývá, že s rizikem omylu nejvýše 5 % se liší cestování tramvají a autobusem a dále cestování autobusem a metrem. 75 10 Neparametrické úlohy o mediánech Příklad 10.1. Párový znaménkový test a párový Wilcoxonův test Při zjišťování kvality jedné složky půdy se používají dvě metody označené A a B. Výsledky jsou uvedeny v následující tabulce: Vzorek 1 2 3 4 5 6 7 8 9 10 11 12 A B 0.275 0.28 0.312 0.312 0.284 0.288 0.3 0.298 0.365 0.361 0.298 0.307 0.312 0.319 0.315 0.315 0.242 0.242 0.321 0.323 0.335 0.341 0.307 0.315 Na hladině významnosti a = 0.05 testujte hypotézu, že metody A a B dávají stejné výsledky. K testování použijte jak párový znaménkový test, tak párový Wilcoxonův test. Pro lepší představu sestrojte krabicové diagramy pro obě metody. Testujeme Hq : zo.50 = 0 proti oboustranné alternativě H\ : zo.50 7^ Oj kde zo.50 Je medián rozdělení, z něhož pochází rozdílový náhodný výběr Z\ = X\—Y\,..., Z\i = X\i—Y\i. Vypočteme rozdíly mezi výsledky metod A a B: Xi-yi I -0.005 0 -0.004 0.002 0.004 -0.009 -0.007 0 0 -0.002 -0.006 -0.008 Párový znaménkový test Nenulových rozdílů je 9, testovací statistika = 2. Ve statistických tabulkách najdeme pro n = 9 a a = 0.05 kritické hodnoty k\ = 1, fc2 = 8. Protože kritický obor W = (0; 1) U 8 neobsahuje hodnotu 2, nemůžeme H0 zamítnout na hladině významnosti a = 0.05. Neprokázalo se tedy, že by metody A a B dávaly rozdílné výsledky. library(PASWR) xl <- c(0.275, 0.312,0.284,0.300,0.365,0.298,0.312,0.315,0.242,0.321,0.335,0.307) x2 <- c(0.280, 0.312,0.288,0.298,0.361,0.307,0.319,0.315,0.242,0.323,0.341,0.315) SIGN.test(xl-x2, alternative='two.sided') ## ## One-sample Sign-Test ## ## data: xl - x2 ## s = 2, p-value = 0.1797 ## alternative hypothesis: true median is not equal to 0 ## 95 percent confidence interval: ## -0.006893636 0.000000000 ## sample estimates: ## median of x ## -0.003 ## Conf.Level L.E.pt U.E.pt ## Lower Achieved CI 0.8540 -0.0060 0 ## Interpolated CI 0.9500 -0.0069 0 ## Upper Achieved CI 0.9614 -0.0070 0 Párový Wilcoxonův test Absolutní hodnoty nenulových rozdílů uspořádáme vzestupně podle velikosti: Usp. abs(a;i-yi) 0.002 0.002 0.004 0.004 0.005 0.006 0.007 0.008 0.009 Pořadí 123456789 Průměrné pořadí 1.5 1.5 3.5 3.5 5 6 7 8 9 76 S+ = 1.5 + 3.5 = 5, S~ = 1.5 + 3.5 + 5 + 6 + 7 + 8 + 9 = 40, n = 9, a = 0.05, tabelovaná kritická hodnota pro n = 9 a a = 0.05 je 5, testovací statistika je = min(5V+, S^) = min(5, 40) = 5. Protože 5 < 5, Hq zamítáme na hladině významnosti a = 0.05. wilcox.test(xl-x2, alternative='two.sideď, correct=F, exact=F) ## ## Wilcoxon signed rank test ## ## data: xl - x2 ## V = 5, p-value = 0.03781 ## alternativě hypothesis: true location is not equal to 0 V tomto případě je p-hodnota 0.03781, tedy nulová hypotéza se zamítá na asymptotické hladině významnosti a = 0.05. Nejsou však splněny předpoklady pro použití asymptotické varianty testu (příliš malý rozsah výběru), i když závěr je stejný jako při testování pomocí kritické hodnoty. Krabicový graf boxplot(xl, x2, names=c('A', 'B'),main='Kvalita složky pudy', ylab='výsledek testu',xlab='metoda',col='burlywoodl',border='chocolate4') Kvalita složky pudy CD "O CD CO CO O OJ o metoda Z krabicových diagramů je vidět, že obě metody se poněkud liší v úrovni, ale neliší se ve variabilitě. Příklad 10.2. Jednovýběrový znaménkový test a jednovýběrový Wilcoxonův test Vyráběné ocelové tyče mají kolísavou délku s předpokládanou hodnotou mediánu 10 m. Náhodný výběr 10-ti tyčí poskytl tyto výsledky: 9.83, 10.10, 9.72, 9.91, 10.04, 9.95, 9.82, 9.73, 9.81, 9.90. Na hladině významnosti 0.05 testujte hypotézu, že předpoklad o mediánu délky tyčí je oprávněný. K testování použijte jak jednovýběrový znaménkový test, tak jednovýběrový Wilcoxonův test. Pro lepší představu sestrojte krabicový diagram. Testujeme Hq : £0.50 = 10 proti oboustranné alternativě H\ : £0.50 7^ 10- Vypočteme rozdíly mezi naměřenými délkami a konstantou 10: 77 Xi-W -0.17 0.1 -0.28 -0.09 0.04 -0.05 -0.18 -0.27 -0.19 -0.1 Jednovýběrový znaménkový test Nenulových rozdílů je 10, testovací statistika = 2. Ve statistických tabulkách najdeme pro n = 10 a a = 0.05 kritické hodnoty k\ = 1, fc2 = 9. Protože kritický obor W = (0; 1) U (9; 10) neobsahuje hodnotu 2, Hq nezamítáme na hladině významnosti a = 0.05. x <- c(9.83, 10.10, 9.72, 9.91, 10.04, 9.95, 9.82, 9.73, 9.81, 9.90) SIGN.test(x, md=10, alternative='two.sided') ## ## One-sample Sign-Test ## ## data: x ## s = 2, p-value = 0.1094 ## alternative hypothesis: true median is not equal to 10 ## 95 percent confidence interval: ## 9.755956 10.010800 ## sample estimates: ## median of x ## 9.865 ## Conf.Level L.E.pt U.E.pt ## Lower Achieved CI 0.8906 9.810 9.9500 ## Interpolated CI 0.9500 9.756 10.0108 ## Upper Achieved CI 0.9785 9.730 10.0400 Jednovýběrový Wilcoxonův test Absolutní hodnoty nenulových rozdílů uspořádáme vzestupně podle velikosti: Usp. abs(xi - yi) 0.04 0.05 0.09 0.1 0.1 0.17 0.18 0.19 0.27 0.28 Pořadí 123456789 10 Průměrné pořadí 1 2 3 4.5 4.5 6 7 8 9 10 5+ = 1 + 4.5 = 5.5, S~ = 2 + 3 + 4.5 + 6 + 7 + 8 + 9+ 10 = 49.5, n = 10, a = 0.05, tabelovaná kritická hodnota pro n = 10 a a = 0.05 je 8, testovací statistika = min(5^r, S^) = min(5.5; 49.5) = 5.5. Protože 5.5 < 8, Hq zamítáme na hladině významnosti a = 0.05. Vidíme, že Wilcoxonův test dospěl k odlišnému závěru než znaménkový test. wilcox.test(x,mu=10, alternative = 'two.sided', exact=F, correct=F) ## ## Wilcoxon signed rank test ## ## data: x ## V = 5.5, p-value = 0.02484 ## alternative hypothesis: true location is not equal to 10 78 Krabicový diagram boxplot(x, main='Ocelové tyče', ylab='delka tyči (v m)', col='gray85') Ocelové tyče Q) "O O O 00 O) Příklad 10.3. Dvouvýběrový Wilcoxonův test Majitel obchodu chtěl zjistit, zda velikost nákupů (v dolarech) placených kreditními kartami Master/EuroCard a Visa jsou přibližně stejné. Náhodně vybral • 7 nákupů placených Master/EuroCard: 42, 77, 46, 73, 78, 33, 37; • 9 nákupů placených Visou: 39, 10, 119, 68, 76, 126, 53, 79, 102. Lze na hladině významnosti a = 0.05 tvrdit, že velikost nákupů placených těmito dvěma typy karet se shodují? Pro lepší představu sestrojte krabicové diagramy pro oba typy platebních karet. Na hladině významnosti a = 0.05 testujeme Hq : £0.50 — ž/0.50 = 0 proti oboustranné alternativě H\ : 2:0.50— yo.50 7^ 0. usp. hodnoty 10 33 37 39 42 46 53 68 73 76 77 78 79 102 119 126 pořadí xt 2 3 5 6 9 11 12 pořadí Vi 1 4 7 8 10 13 14 15 16 Ti = 2 + 3 + 5 + 6 + 9 + 11 + 12 = 48, T2 = 1 + 4 + 7 + 8 + 10 + 13 + 14 + 15 + 16 = 88, U-l = 7.9 + 7.8/2 - 48 = 43, (72 = 7.9 + 9.10/2 - 88 = 20, Kritická hodnota pro a = 0.05, min(7, 9) = 7, max(7, 9) = 9 je 12. Protože min(43, 20) = 20 > 12, nemůžeme na hladině významnosti a = 0.05 zamítnout hypotézu, že velikosti nákupů placených kreditními kartami Master/EuroCard a Visa se shodují. x <- c(42, 77, 46, 73, 78, 33, 37, 9) y <- c(39, 10, 119, 68, 76, 126, 53, 79, 102) wilcox.test(x, y, alternativě = 'two.sideď, exact=F, correct=F) ## ## Wilcoxon rank sum test ## 79 ## data: x and y ## W = 20, p-value = 0.1237 ## alternative hypothesis: true location shift is not equal to 0 Krabicový diagram boxplot(x,y,main='Platebni karty',xlab='platebni karta', ylab='cena za nakup',names=c('Master/Eurocard','Visa'), col='darkolivegreenl',border='darkolivegreen') Platebni karty o OJ Master/Eurocard platebni karta Visa Příklad 10.4. Kruskalův-Wallisův test Voda po holení jisté značky se prodává ve čtyřech různých lahvičkách stejného obsahu. Údaje o počtu prodaných lahviček za týden v různých obchodech jsou uvedeny v následující tabulce: l.typ: 50 35 43 30 62 52 43 57 33 70 64 58 53 65 39 2.typ: 31 37 59 67 44 49 54 62 34 42 40 3.typ: 27 19 32 20 18 23 4.typ: 35 39 37 38 28 33 Posuďte na 5 % hladině významnosti, zda typ lahvičky ovlivňuje úroveň prodeje. V případě zamítnutí nulové hypotézy zjistěte, prodeje kterých typů lahviček se od sebe významně liší. K testování použijte Kruskalův - Wallisův test; v případě zamítnutí nulové hypotézy použijte k zjištění významných rozdílů vhodnou metodu mnohonásobného porovnávání. Pro lepší představu sestrojte krabicové diagramy pro všechny typy lahviček. Všech 38 hodnot uspořádáme vzestupně podle velikosti a stanovíme součet pořadí hodnot patřících do 1. až 4. výběru. Ti = 379, T2 = 257, T3 = 24, T4 = 81 Q=-i2vS-3(n+l) = ^f^ + H^ + ^ + ^)>) -3*39 = 18.79 nin+l)jr[nj 38 * 39 V 15 11 6 6 7 80 W = (xLa(r - 1); oo) = <*£95(3); oo) = (7.815; oo) Protože Q € W, Hq zamítáme na asymptotické hladině významnosti a = 0.05. xl <- c(50, 35, 43, 30, 62, 52, 43, 57, 33, 70, 64, 58, 53, 65, 39) x2 <- c(31, 37, 59, 67, 44, 49, 54, 62, 34, 42, 40) x3 <- c(27, 19, 32, 20, 18, 23) x4 <- c(35, 39, 37, 38, 28, 33) x <- c(xl, x2, x3, x4) ni <- c(length(xl),length(x2),length(x3),length(x4)) group <- c(rep(l, ni[l]), rep(2, ni[2]), rep(3, ni[3]), rep(4, ni [4])) kruskal.test(x, group) ## ## Kruskal-Wallis rank sum test ## ## data: x and group ## Kruskal-Wallis chi-squared = 18.802, df = 3, p-value = 0.0003004 Krabicový graf boxplot(xl, x2, x3, x4, main='Voda po holeni',xlab='typ lahvičky', ylab='pocet prodaných kusu', names=c('Smart','Sport','Atractive','Mystic'),col='lightcyan',border='dodgerblue4') Voda po holeni o o _ co co J= o _ o en >. CO "O P o _ Q. ^í Smart Sport Atractive Mysti typ lahvičky Je vidět, že úroveň prodeje pro l.typ je nevyšší, zatímco pro 3. typ nejnižší. 81 Metoda mnohonásobného porovnávání Nyní provedeme mnohonásobné porovnávání, abychom zjistili, které dvojice typů lahviček se liší na hladině významnosti a = 0.05: library(PMCMR) posthoc.kruskal.nemenyi.test(x=x, g=group, method="Chisq") ## ## Pairwise comparisons using Tukey and Kramer (Nemenyi) test ## with Tukey-Dist approximation for independent samples ## ## data: x and group ## ## 1 2 3 ## 2 0.97310 - ## 3 0.00043 0.00334 - ## 4 0.12541 0.29841 0.44923 ## ## P value adjustment method: none Vidíme, že se liší typy (1, 3) a (2, 3). Příklady k samostatnému řešení Příklad 10.5. Ve skupině 12-ti studentů se sledovala srdeční frekvence při změně polohy z lehu do stoje. Získaly se tyto rozdíly počtu tepů srdce za 1 minutu: -2, 4, 8, 25, -5, 16, 3, 1, 12, 17, 20, 9. Za předpokladu, že tyto rozdíly mají symetrické rozdělení, testujte na hladině významnosti a = 0.05 hypotézu, že medián rozdílů obou tepových frekvencí je 15 proti oboustranné alternativě. Sestrojte krabicový diagram. ## [1] "Wilcoxonuv test: p-value = 0.0499" ## [1] "Znaménkový test: p-value = 0.3877" Krabicový graf srdecni frekvence Výsledek: Zaménkový test nulovou hypotézu nezamítá na hladině významnosti a = 0.05, avšak Wilcoxonuv test ano. 82 Příklad 10.6. Z produkce tří podniků vyrábějících televizory bylo vylosováno 10, 8 a 12 kusů. Byly získány následující výsledky zjišťování citlivosti těchto televizorů v mikrovoltech: 1.podnik: 420 560 600 490 550 570 340 480 510 460 2.podnik: 400 420 580 470 470 500 520 530 3.podnik: 450 700 630 590 420 590 610 540 740 690 540 670 Ověřte na hladině významnosti a = 0.05 hypotézu o shodě úrovně citlivosti televizorů v jednotlivých podnicích. Sestrojte krabicové diagramy pro všechny tři podniky. ## [1] "Kruskaluv-Wallisuv test: p-value= 0.0157" Krabicový graf Citilivost televizoru v různých podnicich > o o E > o o o o LO o o Metoda mnohonásobného porovnávání ## podnikl podnik2 ## podnik2 0.8728 NA ## podnik3 0.0672 0.0251 Výsledek: Na hladině významnosti a = 0.05 se liší televizory vyráběné ve 2. a 3. podniku. 83 11 Hodnocení kontingenčních tabulek Příklad 11.1. Testování hypotézy o nezávislosti, měření síly závislosti V roce 1950 zkoumali Yule a Kendall barvu očí a vlasů u 6800 mužů. Výsledky zkoumání jsou uvedeny v následující tabulce a v souboru vlasy_oci.txt. Barva očí Barva vlasů světlá kaštanová černá rezavá modrá šedá/zelená hnědá 1768 807 180 47 946 1387 746 53 115 438 288 16 Na asymptotické hladině významnosti a = 0.05 testujte hypotézu o nezávislosti barvy očí a barvy vlasů. Vypočtěte Cramérův koeficient. Simultánní četnosti znázorněte graficky. Ověření podmínek dobré aproximace library(scatterplot3d) data <- read.delim('vlasy_oci.csv',sep=';',dec='.',header=T) data <- data.frame(data[,2:5], row.names=data[,1]) názvy.v <- names(data) názvy.o <- row.names(data) # Overení podmnky dobre aproximace chisq.test(data)$expected ## svetla kaštanová cerna rezavá ## modra 1167.2593 1085.976 500.9024 47.86217 ## seda/zelena 1304.7310 1213.875 559.8952 53.49904 ## hneda 357.0097 332.149 153.2025 14.63879 Podmínky dobré aproximace jsou splněny. Všechny teoretické četnosti jsou větší než 5. Nyní budeme testovat hypotézu o nezávislosti proměnných oci a vlasy. chisq.test(data) ## ## Pearson's Chi-squared test ## ## data: data ## X-squared = 1088.1, df = 6, p-value < 2.2e-16 Ve výstupní tabulce najdeme mj. hodnotu testové statistiky K = 1088.149 s počtem stupňů volnosti (df= 6) a odpovídající p-hodnotou (p-value < 2.2e — 16). Protože p-hodnota je mnohem menší než 0.05, nulovou hypotézu o nezávislosti barvy očí a barvy vlasů zamítáme na asymptotické hladině významnosti a = 0.05. Pro zjištění míry závislosti v kontingenční tabulce použijeme Cramérův koeficient. library(lsr) cramersV(data) ## [1] 0.2830494 Hodnota Cramérova koeficientu je 0.283, což svědčí o slabé závislosti barvy očí a vlasů. 84 Grafické znázornění četností x=rep(c(l,2,3,4),c(3,3,3,3)) y=rep(l:3,4) z=c(as.matrix(data)) #cbind(x,y, z) scatterplot3d(x, y , z, type='h', pch=18, xlab='barva vlasu', ylab='barva oci', zlab='početnosti', main='Scatterplot', x.ticklabs = c('světla','','kaštanová','','cerna','','rezavá'), y.ticklabs = cCmodra', '', 'seda/zelena', '','hneda'),color=rep('brown4',12), label.tick.marks=T, axis=T, tick.marks=T, angle=40, lwd=2) barva vlasu Příklad k samostatnému řešení Příklad 11.2. Otevřete si soubor ped_hodnost.txt. Na hladině významnosti a = 0.05 testujte hypotézu o nezávislosti pedagogické hodnosti a pohlaví. Dále vypočtěte Cramérův koeficient vyjadřující intenzitu závislosti pedagogické hodnosti na pohlaví. Data v souboru mají následující tvar: pohlaví pedagogická hodnost odb. asistent docent profesor muž žena 32 15 8 34 8 3 ## [1] "Podminky dobre aproximace:" ## odb.asistent docent profesor ## muz 36.3 12.65 6.05 ## zena 29.7 10.35 4.95 Podmínky dobré aproximace jsou splněny, pouze jediná teoretická četnost klesne pod 5. 85 ## [1] "Chi-kvadratovy test:" ## ## Pearson's Chi-squared test ## ## data: data ## X-squared = 3.4988, df = 2, p-value = 0.1739 Testovací statistika K nabývá hodnoty 3.5, p-hodnota = 0.1739, tedy na asymptotické hladině významnosti a = 0.05 nezamítáme hypotézu o nezávislosti pedagogické hodnosti a pohlaví. ## [1] "Crameruv koeficient: V= 0.187" Hodnota Cramérova koeficientu je 0.187, což svědčí o slabé závislosti mezi pedagocickou hodností a pohlavím. Grafické znázornění četností Scatterplot odborný asistent docent profesor pedagogická hodnost Příklad 11.3. Fisherův faktoriálový test 100 náhodně vybraných mužů a žen bylo dotázáno, zda dávají přednost nealkoholickému nápoji A či B. Údaje jsou uvedeny ve čtyřpolní kontingenční tabulce. pref. nápoj pohlaví muž žena A B 20 30 30 20 Na hladině významnosti a = 0.05 testujte pomocí Fisherova faktoriálového testu hypotézu, že preferovaný typ nápoje nezáleží na pohlaví respondenta. V našem případě se jedná o oboustranný test (nevíme, zda muži více preferují nápoj A či nápoj B než ženy. data <- data.frame(muz=c(20,30), zena=c(30,20), row.names=c('A','B')) fisher.test(data, alternative='two.sideď) 86 ## ## Fisher's Exact Test for Count Data ## ## data: data ## p-value = 0.07134 ## alternative hypothesis: true odds ratio is not equal to 1 ## 95 percent confidence interval: ## 0.1846933 1.0640121 ## sample estimates: ## odds ratio ## 0.4481632 Ve výstupní zprávě funkce fisher.test() je uvedena p-hodnota = 0.07134. Protože p-hodnota je větší než 0.05, nezamítáme na hladině významnosti a = 0.05 hypotézu, že preferovaný typ nápoje nezáleží na pohlaví respondenta. Příklad 11.4. Podíl šancí Pro údaje z příkladu č.3 vypočtěte podíl šancí a sestrojte 95% asymptotický interval spolehlivosti pro logaritmus podílu šancí. Pomocí tohoto intervalu spolehlivosti testujte na asymptotické hladině významnosti a = 0.05 hypotézu, že preferovaný typ nápoje nezáleží na pohlaví respondenta. Nejprve zopakujme teorii: Ve čtyřpolních tabulkách používáme charakteristiku OR ad bč' která se nazývá podíl šancí (odds ratio). Můžeme si představit, že pokus se provádí za dvojích různých okolností a může skončit buď úspěchem nebo neúspěchem. výsledek pokusu okolnosti nj. I. II. úspěch a b a+b neúspěch c d c+d a+c b+d n Poměr počtu úspěchů k počtu neúspěchů (tzv. šance) za prvních okolností je ^, za druhých okolností je ^. Podíl šancí je OR = Považujeme ho za odhad skutečného podílu šancí op. Pomocí 100(1 — a)% asymptotického intervalu spolehlivosti pro logaritmus skutečného podílu šancí ln op lze na asymptotické hladině významnosti a testovat hypotézu o nezávislosti nominálních veličin X &Y. Upozornění: Musí být splněny podmínky dobré aproximace. Asymptotický 100(1 — a)% interval spolehlivosti pro přirozený logaritmus skutečného podílu šancí má tvar \nOR 1111 - + t H---1- jiii-Q/2 5 In O R abed ' 1111 - + T H---h -,Ua/2 a b c d ' Jestliže interval spolehlivosti nezahrne 0, pak hypotézu o nezávislosti zamítneme na asymptotické hladině významnosti a. Výpočet příkladu Ověříme splnění podmínek dobré aproximace a zjistíme, že všechny teoretické četnosti jsou rovny 25. 87 # Ověřeni podminek dobré aproximace a <- 20 b <- 30 c <- 30 d <- 20 alpha <- 0.05 data <- data.frame(muz=c(a,c), zena=c(b,d), row.names=c('A','B')) chisq.test(data)$expected ## muz zena ## A 25 25 ## B 25 25 Podíl šancí ad 20 * 20 4 ^"=30730 = 9=0-4 Dolní a horní mez 95% intervalu spolehlivosti pro ln op (-1,61108; -0,01078) Protože tento interval spolehlivosti neobsahuje 0, na asymptotické hladině významnosti a = 0.05 zamítáme hypotézu, že preferovaný typ nápoje nezáleží na pohlaví respondenta. # Podii šanci + IS (OR <- (a*d)/(b*c)) ## [1] 0.4444444 (dh <- log(0R)+sqrt(l/a+l/b+l/c+l/d)*qnorm(alpha/2)) ## [1] -1.611082 (hh <- log(0R)+sqrt(l/a+l/b+l/c+l/d)*qnorm(l-alpha/2)) ## [1] -0.01077827 Tento výsledek je v rozporu s výsledkem, ke kterému dospěl Fisherův přesný test. Je to způsobeno tím, že test pomocí asymptotického intervalu spolehlivosti je pouze přibližný. Ke stejnému závěru, jaký jsme dostali u testování pomocí podílu šancí, dospějeme, pokud použijeme Pearsonův chí-kvadrát test o nezávislosti. chisq.test(data, correct=F) ## ## Pearson's Chi-squared test ## ## data: data ## X-squared =4, df = 1, p-value = 0.0455 Ve funkci chisq.test() však můžeme zadat parametr correct=T, který provede korekci Pearsonova testu pro kontingenční tabulky typu 2x2. Výsledek takto provedeného testu je již v souladu s Fisherovým přesným testem. chisq.test(data, correct=T) ## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: data ## X-squared = 3.24, df = 1, p-value = 0.07186 88 Příklad 11.5. 36 mužů onemocnělo určitou chorobou. Někteří z nich se léčili, jiní ne. Někteří se uzdravili, jiní zemřeli. Údaje jsou uvedeny ve čtyřpolní kontingenční tabulce. přežití léčení ano ne ano ne 10 6 12 8 Vypočtěte a interpretujte podíl šancí. Pomocí intervalu spolehlivosti pro logaritmus podílu šancí testujte na asymptotické hladině významnosti a = 0.05 hypotézu, že přežití nezávisí na léčení, proti tvrzení, že léčení zvyšuje šance na přežití. ## muz zena ## A 9.777778 6.222222 ## B 12.222222 7.777778 ## [1] "0R= 1.1111" ## [1] "dolni hranice IS: -1.0283" Výsledek: OR = 1.1; nulovou hypotézu nezamítáme asymptotické hladině významnosti a = 0.05, protože le-vostranný 95% asymptotický interval spolehlivosti pro logaritmus podílu šancí je (-1.03; oo). Příklad k samostatnému řešení Příklad 11.6. V průzkumu o kuřáctví bylo dotázáno 92 osob. Z 64 mužů jich kouří 19 a z 28 žen jich kouří 6. a) Na hladině významnosti a = 0.05 testujte hypotézu, že kouření se vyskytuje stejně často u mužů a žen. Použijte Pearsonův chi-kvadrát test i Fisherův přesný test. b) Vypočtěte a interpretujte podíl šancí a stanovte meze 95% intervalu spolehlivosti pro podíl šancí. Výsledek ad a) Před provedením Pearsonova chí-kvadrát testu je zapotřebí ověřit splnění podmínek dobré aproximace. ## muz zena ## kurak 17.3913 7.608696 ## nekuřák 46.6087 20.391304 Jsou splněny, všechny čtyři teoretické četnosti jsou větší než 5. ## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: data ## X-squared = 0.31889, df = 1, p-value = 0.5723 Testovací statistika Pearsonova chí-kvadrát testu je K = 0.6714, p-hodnota je 0.5723 > 0.05, tedy nulovou hypotézu nezamítáme na asymptotické hladině významnosti a = 0.05. 89 ## ## Fisher's Exact Test for Count Data ## ## data: data ## p-value = 0.4576 ## alternative hypothesis: true odds ratio is not equal to 1 ## 95 percent confidence interval: ## 0.498056 5.398695 ## sample estimates: ## odds ratio ## 1.54109 Pro Fisheruv presný test vychází p-hodnota 0.4576, což je větší než hladina významnosti 0.05, nulovou hypotézu nezamítáme na hladině významnosti a = 0.05. ad b) ## [1] "0R= 1.5481" ## [1] "dolni hranice IS: 0.5418" ## [1] "horni hranice IS: 4.4239" Podíl šancí je 1.55, což znamená, že u mužů je šance na kouření 1.55x vyšší než u žen. 0.5418 < op < 4.4239 s pravděpodobností aspoň 0.95. 90 12 Jednoduchá korelační analýza Příklad 12.1. Testování nezávislosti ordinálních veličin 12 různých softwarových firem nabízí speciální programové vybavení pro vedení účetnictví. Jednotlivé programy byly posouzeny odbornou komisí složenou z počítačových odborníků a komisí složenou z profesionálních účetních. Úkolem bylo doporučit vhodný program na základě stanovení pořadí jednotlivých programů. Výsledky posouzení: Produkt firmy číslo 1 2 3 4 5 6 7 8 9 10 11 12 Pořadí dle odborníků Pořadí dle účetních 6 7 1 8 4 2.5 9 12 10 2.5 5 11 4 5 2 10 6 1 7 11 8 3 12 9 Vypočtěte Spearmanův koeficient pořadové korelace a na hladině významnosti a = 0.05 testujte hypotézu, že hodnocení obou komisí jsou nezávislá. Data jsou uložena v souboru ucetnictvi.txt. X <- c(6, 7, 1, 8, 4, 2.5, 9, 12, 10, 2.5, 5, 11) Y <- c(4, 5, 2, 10, 6, 1, 7, 11, 8, 3, 12, 9) cor.test(X, Y, method='spearman', exact=F) ## ## Spearman's rank correlation rho ## ## data: X and Y ## S = 81.642, p-value = 0.009024 ## alternativě hypothesis: true rho is not equal to 0 ## sample estimates: ## rho ## 0.714537 Spearmanův koeficient pořadové korelace nabývá hodnoty r s = 0.7145, tedy mezi hodnocením obou komisí existuje vysoký stupeň přímé pořadové závislosti. Testovací statistika se realizuje hodnotou 81.642, odpovídající p-hodnota je 0.009024, tedy na asymptotické hladině významnosti a = 0.05 zamítáme hypotézu o pořadové nezávislosti hodnocení dvou komisí ve prospěch oboustranné alternativy. Upozornění: Pokud rozsah výběru nepřesáhne 20, měli bychom testování provést pomocí tabelované kritické hodnoty. V našem případě pro n = 12 a a = 0.05 je kritická hodnota 0.5804. Vidíme, že nulovou hypotézu zamítáme na hladině významnosti a = 0.05, protože 0.7145 > 0.5804. Příklad 12.2. Testování nezávislosti intervalových a poměrových veličin Zjišťovalo se, kolik mg kyseliny mléčné je ve 100 ml krve matek prvorodiček (veličina X) a u jejich novorozenců (veličina Y) těsně po porodu. Byly získány tyto výsledky: Číslo matky i 2 3 4 5 6 40 64 34 15 57 45 Ví 33 46 23 12 56 40 Nakreslete dvourozměrný tečkový diagram, vypočtěte výběrový korelační koeficient, sestrojte 95% interval spolehlivosti pro korelační koeficient a na hladině významnosti a = 0.05 testujte hypotézu o nezávislosti výsledků obou měření. Data jsou uložena v souboru kyselina_mlecna.txt. 91 Dvourozměrný tečkový diagram data <- read.delimCkyselina_mlecna.txt', header=F) names(data) <- c('matka', 'dite') X <- data$matka Y <- data$dite plot(X, Y, xlab='množství kyseliny mlecne v krvi matky (mg/100 ml)', ylab='mnozstvi kyseliny mlecne v krvi novorozence (mg/100 ml)', main='Mnozstvi kyseliny mlecne v krvi', xlim=c(0,70), ylim=c(0,70), pch=21, col='darkreď , bg='reď) mtext('Graficke ověřeni normality', line=0.4) Mnozstvi kyseliny mlecne v krvi Grafické ověřeni normality 10 20 T" "T 50 T" 30 40 50 60 70 mnozstvi kyseliny mlecne v krvi matky (mg/100 ml) Testování hypotézy o nezávislosti: cor.test(X, Y, type='pearson') ## ## Pearson's product-moment correlation ## ## data: X and Y ## t = 5.2653, df = 4, p-value = 0.006232 ## alternativě hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.5108072 0.9930174 ## sample estimates: ## cor ## 0.9348324 Ve výstupní zprávě funkce cor.test() je mj. uvedena hodnota výběrového korelačního koeficientu r 12 = 0.9348, tzn. že mezi X a Y existuje velmi vysoký stupeň přímé lineární závislosti. Hodnota testovací statistiky t = 5.2653 a p-hodnota pro test hypotézy o nezávislosti vychází 0.006232, Hq tedy zamítáme na hladině významnosti a = 0.05. S rizikem omylu nejvýše 5% jsme tedy prokázali, že mezi oběma koncentracemi existuje závislost. 95% interval spolehlivosti pro p mající meze 0.5108 a 0.9930 nepokrývá hodnotu 0, a tudíž hypotézu o nezávislosti veličin X, Y zamítáme na hladině významnosti a = 0.05. 92 Příklad 12.3. Porovnání dvou korelačních koeficientů V psychologickém výzkumu bylo vyšetřeno 426 hochů a 430 dívek. Ve skupině hochů činil výběrový koeficient korelace mezi verbální a performační složkou IQ 0.6033, ve skupině dívek činil 0.5833. Za předpokladu dvourozměrné normality dat testujte na hladině významnosti a = 0.05 hypotézu, že korelační koeficienty se neliší. Rl <- 0.6033 R2 <- 0.5833 ni <- 426 n2 <- 430 ksi <- 0 Zl <- l/2*log((l+Rl)/(l-Rl)) Z2 <- l/2*log((l+R2)/(l-R2)) Zw <- (Zl-Z2-ksi)/sqrt(l/(nl-3)+l/(n2-3)) (p.val <- 2*min(pnorm(Zw), l-pnorm(Zw))) ## [1] 0.6527169 Výsledek: Příslušná p-hodnotaje 0.6528, tedy nezamítáme nulovou hypotézu o shodě dvou koeficientů korelace na asymptotické hladině významnosti a = 0.05. Příklady k samostatnému řešení Příklad 12.4. Načtěte datový soubor IQ.txt. Za předpokladu dvourozměrné normality dat (orientačně ověřte pomocí dvourozměrného tečkového diagramu) testujte na hladině významnosti a = 0.1 hypotézu, že korelační koeficienty mezi verbální a performační složkou IQ jsou stejné u dětí z města a venkova. ## [1] 0.0780111 Bydliště = mesto Grafické ověřeni normality verbálni složka IQ (v bodech) E Bydliště = venkov Grafické ověřeni normality 80 90 verbálni složka IQ (v bodech) Výsledek p-hodnota= 0.07801, tedy s rizikem omylu nejvýše 10% jsme prokázali, že korelační koeficienty se liší. 93 Příklad 12.5. V náhodném výběru 10 dvoučlenných domácností byl zjišťován měsíční příjem (veličina X, v tisících Kč) a vydání za potraviny (veličina Y, v tisících Kč). 15 21 34 35 39 42 58 64 75 90 Ví 3 4.5 6.5 6 7 8 9 8 9.5 10.5 Vypočtěte a interpretujte výběrový koeficient korelace. Na hladině významnosti a = 0.05 testujte hypotézu o nezávislosti veličin X, Y. Sestrojte 95% asymptotický interval spolehlivosti pro p. Data jsou uložena v souboru prijem_vydani.txt. Grafické ověření normality Závislost přijmu a vydajú domácnosti Grafické ověřeni normality prijem domácnosti (v tisicich Kc) Výsledek ri2 = 0.9405, mezi měsíčními příjmy a výdaji tedy existuje velmi vysoký stupeň přímé lineární závislosti, p-hodnota= 5.095 e — 05, tedy Hq zamítáme na hladině významnosti a = 0.05. S pravděpodobností alespoň 0.95 platí: 0.7623 < p < 0.9862. Příklad 12.6. Bylo sledováno 10 žáků. Na základě psychologického vyšetření byli tito žáci seřazeni podle nervové lability (čím byl žák labilnější, tím dostal vyšší pořadí Ri). Kromě toho sledování žáci dostali pořadí Qi na základě svých výsledků v matematice (nejlepší žák v matematice dostal pořadí 1). Výsledky jsou uvedeny v tabulce: Pořadí Ri 1 2 3 4 5 6 7 8 9 10 Pořadí Qi 9 3 8 5 4 2 10 1 7 6 Vypočtěte vhodný korelační koeficient a jeho hodnotu řádně interpretujte. Na hladině významnosti a = 0.05 testujte hypotézu, že nervová labilita a výsledky v matematice jsou nezávislé. Data jsou uložena v souboru nervova_labilita.txt 94 Závislost mezi labilitou zaka a výsledky v matematice Grafické ověřeni normality 0 2 4 6 8 10 poradi zaka podle nervové lability Výsledek: Spearmanův koeficient pořadové korelace r s = —0.127, tedy mezi nervovou labilitou žáka a jeho výsledky v matematice existuje nízký stupeň nepřímé pořadové závislosti. p-hodnota= 0.7329, a tedy Hq nezamítáme na hladině významnosti a = 0.05. 95