# NMFM301 # Praktické cvičení 2 -- Instruktáž ************** ## Kvantilový diagram Mějme náhodný výběr $X_1,\ldots,X_n$ z rozdělení $F_X$. Kvantilový diagram je grafická metoda sloužící k neformálnímu porovnání rozdělení $F_X$ s hypotetickým rozdělením $F_0$. Na ose $y$ se kreslí pořádkové statistiky $X_{(i)}$. Na ose $x$ se vynášejí aproximace střední hodnoty pořádkových statistik rozdělení $F_0$, tj. $F_0^{-1}\bigl(\frac{i}{n+1}\bigr)$ *[proč zrovna tohle?]*. Jestliže data pocházejí z rozdělení $F_0$, vykreslené body leží přibližně na přímce $y=x$. Jestliže data pocházejí z lineární transformace rozdělení $F_0$, vykreslené body stále leží přibližně na přímce (ale jiné než $y=x$). Jestliže data pocházejí z rozdělení odlišného od $F_0$, vykreslené body tvoří rostoucí křivku. Z tvaru této křivky dokážeme poznat, v čem se rozdělení dat liší od $F_0$. V R se dá porovnávat rozdělení dat s rozdělením N(0,1) pomocí funkce `qqnorm`. Argumentem funkce je vektor obsahující náhodný výběr. Vygenerujme 50 pozorování z rozdělení N(10,4) a namalujme normální kvantilový diagram: ```{r} x <- rnorm(50,10,2) qqnorm(x) qqline(x) ``` Funkce `qqline` namaluje přímku procházející dolním a horním výběrovým kvartilem dat, která usnadnuje čtení grafu. Nyní vygenerujme 50 pozorování z rozdělení Exp(0.4) a namalujme opět normální kvantilový diagram: ```{r} x <- rexp(50,0.4) qqnorm(x) qqline(x) ``` Poznali byste, že tato data nepocházejí z normálního rozdělení? ### Úkol: Vytvořte normální kvantilový diagram pro data z rozdělení $\Gamma(a,p)$ a sledujte, co se děje, když se zvětšuje argument $p$. 1. Vygenerujte 50 pozorování z rozdělení $\Gamma(a,p)$, kde *a*=0.4 a *p*=0.5. Musíte dát pozor na správnou parametrizaci Gama rozdělení; použijte volání `rgamma(50,rate=0.4,shape=0.5)` (argument `rate` udává hodnotu parametru *a*, argument `shape` udává hodnotu parametru *p*, oba musíte pojmenovat). 2. Nakreslete normální kvantilový diagram pro tato data. 3. Zvětšujte hodnotu parametru *p* na 1, 2, 20, 50 a sledujte, co se děje s kvantilovým diagramem. Jaké to má vysvětlení? ## Práce s daty Nyní se naučíme pracovat s datovými soubory. Načtěte data z pracovního direktoráře nebo přímo z webu příkazem ```{r} load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/cv2.RData")) ?load ``` Nyní máte k dispozici datový soubor nazvaný `cv2`. Načteme ješte několik funkcí, s kterými dnes budeme pracovat (`plotCI`,`plotmeans`,`ci.asym`,`ci.t`,`sign.test`). ```{r} load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/functions.RData")) ``` Data mají 500 pozorování (řádků) a 6 veličin (sloupců). Jejich velikost zjistíme funkcí `dim`: ```{r} dim(cv2) ``` Pozorování tvoří náhodný výběr 500 obyvatel USA. V jednotlivých veličinách o nich máme následující údaje: Název | Význam ------|------- `id` | ID jedince `pohl` | Pohlaví `vzdelani` | Vzdělání `stav` | Rodinný stav `vaha` | Hmotnost [kg] `vyska` | Výška [cm] Data si prohlédněte poklepáním na příslušný řádek v okénku "Workspace". Strukturu dat vypíšeme pomocí funkce `str`: ```{r} str(cv2) ``` Základní popisné statistiky všech veličin získáme pomocí funkce `summary`: ```{r} summary(cv2) ``` S datovým souborem se zachází podobně jako s maticí. Můžeme vypsat jeden nebo více řádků pomocí hranaté závorky: ```{r} cv2[1:5,] ``` K jednotlivým veličinám máme přístup buď pomocí čísla sloupce (`cv2[,2]` je vektor hodnot pohlaví) nebo lépe pomocí názvu veličiny napsané za názvem dat a oddělené znakem $, např `cv2$pohl`. Výšky prvních pěti obyvatel vypíšeme takto: ```{r} cv2$vyska[1:5] ``` Hmotnosti lidí se základním vzděláním získáme pomocí podmínky v hranaté závorce (rovnítko musí být zdvojeno): ```{r} cv2$vaha[cv2$vzdelani=="Elementary"] ``` Zkusme se nyní podívat na normální kvantilový diagram výšky. ```{r} qqnorm(cv2$vyska) qqline(cv2$vyska) ``` Obecně se věří, že výška v populaci má normální rozdělení. Nasvědčuje tomu tento obrázek? Podívejme se ještě na výšky mužů a žen zvlášť. ```{r} qqnorm(cv2$vyska[cv2$pohl=="Male"]) qqline(cv2$vyska[cv2$pohl=="Male"]) ``` ```{r} qqnorm(cv2$vyska[cv2$pohl=="Female"]) qqline(cv2$vyska[cv2$pohl=="Female"]) ``` Jsou výšky pro každé pohlaví normální? ## Intervaly spolehlivosti pro střední hodnotu Funkce `ci.t` počítá 95% interval spolehlivosti pro střední hodnotu založený na t-rozdělení: \[ \left(\overline{X}_{n}-\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr), \ \overline{X}_{n}+\frac{S_n}{\sqrt{n}}t_{n-1}\Bigl(1-\frac{\alpha}{2}\Bigr)\right). \] (pro $\alpha=0.05$). Je to přesný interval pro normální rozdělení a přibližný pro jakékoli rozdělení s konečným druhým momentem. ```{r} smp <- rnorm(20,2,1) prum <- mean(smp) ci.t(smp) ``` Výsledkem funkce `ci.t` je vektor o třech složkách: dolní konec intervalu, průměr, horní konec intervalu. Vygenerujme nyní 100 náhodných výběrů o rozsahu 20 se střední hodnotou 2. Nejprve nastavíme požadované vstupy: ```{r} nobs <- 20 nvyb <- 100 str.h <- 2 ``` Nyní vygenerujeme 2000 veličin z rozdělení N(2,1) a uspořádáme je do matice 20 x 100. Ve sloupcích bude 100 výběrů o rozsahu 20 z rozdělení N(2,1). ```{r} vybery <- rnorm(nobs*nvyb,mean=str.h,sd=1) data.mat <- matrix(vybery,nrow=nobs,ncol=nvyb) ``` Na každý sloupec (výběr) spočítáme interval spolehlivosti pro střední hodnotu funkcí `ci.t`: ```{r} vs.ci <- apply(data.mat,2,ci.t) ``` Zjistíme, kolik intervalů pokrylo skutečnou střední hodnotu (2) a jaká byla průměrná šířka intervalu. Kolik intervalů by ji teoreticky mělo pokrýt? ```{r} (pokryti <- sum(vs.ci[1,]78),length(cv2$vaha)) ``` Otestujte nyní přesným znaménkovým testem hypotézy, že medián hmotnosti mužů (žen) je 78 kg. Liší se nějak záSadně výsledky přesného testu od testu asymptotického? Proč asi?