Cvičenie 7 (praktické cvičenie 1)

(Základy práce so statistickým softwarom R)

Zdrojový soubor pro Knit: Rmd soubor

I. Štatistický software R (stručný úvod)

Štatistický softwar R (R Development Core Team, 2016) je výpočetný násttroj určeným k spracovaniu, príprave a následnej štatistickej analýze dat. Jedná sa o GNU projekt, podobný programu S a primárne je navrhnutý a určený pre štatisticku analýzu dat.

Program R (dostupný pod GNU GPL licenciou) je k dispozícii k stiahnutiu (free of charge) na adrese

https://www.r-project.org

K dispozícii sú distribúcie s priamou podporou pre OS Windows, Linux aj Macintosh.

Základnú inštaláciu programu R je možne jednoducho rozšíriť pomocou dodatočných knižníc (balíčkov), ktoré sú k dispozícii na rôznych online repozitároch (zoznam hlavných repozitárov je na adrese https://cran.r-project.org/mirrors.html). Jednotlivé R knižnice sú tvorené samotnými užívateľmi softwaru R a ich správne fungovanie nie je garantované - je preto namieste určitá opatrnosť a hlavne aktívne premýšľanie pri ich používaní.

Pre užívateľov programu R sú k dispozícii aj rôzne grafické rozhrania, ktoré je možne dodatočne nainštalovať a umožňujú (v určitých smeroch) jednoduchšiu a prehľadnejšiu prácu. Najznámejší a pravdepodobne aj jeden z najlepších R interfacov je RStudio.

Užitočné materiály pre prácu so štatistickým softwarom R

Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (PDF súbor)
Komárek, A.: Základy práce s R. (PDF súbor)
Kulich, M.: Velmi stručný úvod do R. (PDF súbor)
De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808)

V nasledujúcom uvedieme niekoľko jednoduchých príkladov a úloh, ktoré maju primárne slúžiť ako názorné ukážky toho, ako funguje program R, ako sa v ňom implementuju jednotlivé príkazy a akym spôsobom načítať vlastny datový súbor (nejedná sa ale o kompletný súhrn jednotlivých možnosti, ktoré sú v programe R k dispozícii).

V prípade potreby/záujimu je možné využiť podrobnejšieho “sprievodcu” na cestu základmi Rka: Hrátky s R (autor: doc. Arnošt Komárek);
(k Hrátkam s R sú potrebné datove súbory: auta2004.dat, auta2004.csv a auta2004.xls)

1. Základné matematické operácie v R

2 + (4 * 5)^2 ## za '#' nasleduje komentar, ktory R ignoruje
1:5 * 2
rep(4,3)

Výsledky je možné vždy uložiť do nejakého vhodne zvoleného objektu, na ktorý sa následne môžeme odvolávať.

v1 <- 2 + (4 * 5)^2 ## za '#' nasleduje komentar, ktory R ignoruje
v2 <- 1:5 * 2
v3 <- seq(1,5, length = 3)
v4 <- rep(4,3)

Jednotlivé výsledky je možné kombinovať a opäť uložiť ako nový výsledok (resp. prepísať pôvodný).

v1 + v2
v5 <- v3 - 2
v5 <- v5 * 2

Užitočné

Program R je objektový nástroj a zvláda aj matematické operacie, ktoré štandardne nie sú definované. Napr. súčet dvoch vektorov s rôznou dĺžkou, operácie s maticami, ktoré maju rôzne typy/dimenzie. Je podstatné správne pochopiť, ako R v takýchto prípadoch funguje. Zabráni sa možným problémom v budúcnosti a v mnohých prípadoch to zase na druhu stranu može zjednodušiť výpočet a ušetriť čas.

V programe R existuje niekoľko prikazov, ktoré vytvoria vektor: napr. c(), seq(), rep(), repliate() a pod.;
Podobne existuje niekoľko príkazov, ktoré vytvoria maticu: napr. cbind(), rbind(), matrix(), a pod.;
Ku každému príkazu v R je k dispozícii návod, ako príkaz funguje - help sa zavolá zadaním samotného príkazu, ktorému predchádza otazník: napr. ?c(), ?seq(), ?rep();
Pre maticové násobenie je užitočné poznať operátor %*%, ktorý výnasobi dve matice (ak majú správne rozmery). Vhodnosť rozmerov sa jednoducho skontroluje pomocou príkazu dim() (užitočne je poznať tiež príkazy nrow() a ncol()).

Štandardne binárne operátory +, -, *, /, ^ fungujú ‘po zložkách’ - toto ‘fungovanie’ je ale R-špecifické - porovnajte nasledujúce výstupy:

1:4 + 1 ### sucet po zlozkach
1:4 * 1:2 ### nasobenie po zlozkach v paroch (dlzka jedneho vektora je celociselny nasobok dlzky druheho)
1:4 / 1:4 ### nasobenie v paroch (rovnake dlzky oboch vektorov)
1:5 + 1:2 ## Error - preco?

2. Práca s vektormi a maticami

V prípade, že máme v objektoch uložené nejake štruktúry (vektor, maticu, pole, list, atď.), je možné sa v Rku odkazovať aj na jednotlivé elementy v daných objektoch.

v3[1] ### zobrazi prvu zlozku vektora v3
v3[-c(2,3)] ### rovnaky vystup, vektor v3 bez druhej a tretej zlozky

A to isté platí aj pre matice (a analogicky aj pre pole).

m1 <- rbind(v3, v3 * 2, v3 * 3)
m1[1:2,1:2] ## elementy z prveho a druheho riadku a stlpca matice m1

Príkazy je možné ľubovoľne kombinovať a vytvárať nové, zložitejšie objekty.

m2 <- rbind(cbind(m1[1:2,1:2], c(0,0)), seq(10,11, length = 3))
l1 <- list(m1, m2) ## list, ktory obsahuje dve rozne matice, m1 a m2
l1[[2]][1,] # vypise na obrazovku prvy riadok matice m2

Užitočné

Pre maticové násobenie je užitočné poznať operátor %*%, ktorý výnasobi dve matice (ak majú správne rozmery). Vhodnosť rozmerov sa jednoducho skontroluje pomocou príkazu dim() (užitočne je poznať tiež príkazy nrow() a ncol()).

Štandardne binárne operátory +, -, *, /, ^ fungujú ‘po zložkách’ - toto ‘fungovanie’ je ale R-špecifické - porovnajte nasledujúce výstupy:

1:4 + 1 ### sucet po zlozkach
1:4 * 1:2 ### nasobenie po zlozkach v paroch (dlzka jedneho vektora je celociselny nasobok dlzky druheho)
1:4 / 1:4 ### nasobenie v paroch (rovnake dlzky oboch vektorov)
1:5 + 1:2 ## Error - preco?

3. Náhodné data - generátory náhodných čísel

Program R je primárne štatistický nástroj - umožňuje preto aj prácu s náhodnymi hodnotami. Implementované je celá rada rôznych generátorov náhodných čísel. Pre rôzne pravdepodobnostné rozdelenia sú k dispozícii príslušne generátory, ktoré simulujú náhodné hodnoty z daného pravdepodobnostného rozdelenia.

Diskrétne rozdelenia: rbinom(), rpois(), rgeom(), a ďalšie;
Spojité rozdelenia: runif(), rexp(), rnorm(), a ďalšie;

Pre správne použitie geneátorov je potrebné preštudovať prislušný návod - napr. ?rnorm zobrazí help pre generátor náhodných čísel s normálnym rozdelením. Je vhodné si tiež všimnúť, že existujú podobné príkazy, označené ako pnorm(), dnorm() a rnorm(), ktoré postupne značia distribučnú funkciu daného rozdelenia, hustotu a kvantilovú funkciu.

x <- rnorm(n=100,mean=0,sd=1)
y <- rexp(n = 1000)

Užitočné

Je nutné dôsledne si preštudovať implementáciu jednotlivých generátorov a funkcií, ktoré súvisia s pravdepodobnostnými rozdeleniami - je totiž pomerne časté, že parametre, ktoré špecifikujú nejaké rozdelenie, majú různy význam. Napríklad:
```
rnorm(10, 0, 2) ## simuluje 10 pozorovani z N(0, 4) - sigma^2 = 4
rexp(10, 4) ## simuluje 10 pozorovani z Exp(lambda = 4), kde EX = 1/4
```
z důvodu umožnění replikovatelnosti celého experimentu je nutné pri používaní náhodných generátoru používat príkaz set.seed(), kde uvnitř závorky je nutne specifikovat nějakou počáteční numerickou hodnotu. Porovnejte mezi sebou vústup z predchozích dvou příkazu a výstup z nasledujících dvou přikazu:
```
set.seed(1234)

rnorm(10, 0, 2) ## ted ma kazdy stejne nasimulovane hodnoty
rexp(10, 4) ## 
```

4. Reálne data - načítanie vlastných datových súborov

V prvom rade, v programe R je k dispozícii niekoľko skutočných datových súborov. Ich zoznam je možne zobrazit zavolaním príkazu data(). Konkrétny datovy súbor následne zobrazíme pomocou názvu prislušného súboru - napr. Orange. Vo väčšíne sa jedná o data.frame - resp. datovu maticu. Na rozdiel od klasickej matice umožňuje data.frame rôzne datové typy v jednotlivých stĺpcoch.

stlpec1 <- seq(1:5)
stlpec2 <- c(rep("male", 3), rep("female", 2))
stlpec3 <- c("a", "b", "c", "d", "e")

### nie je mozne vytvorit maticu
cbind(stlpec1, stlpec2, stlpec3)
### ale je mozne vytvorit data.frame
data.frame(stlpec1, stlpec2, stlpec3)

K načítaniu vlastných datových súborov slúži niekoľko príkazov. V programe R je technicky možné načítať vpodstate akýkoľvek datovy súbor (txt, csv, xls, ale aj jpg, png, mp3 a mnoho ďalších). Najčastejšie používané sú read.table() a read.csv(), ktoré vyžaduju súbor uložený vo formáte txt resp. csv. Dodatočné parametre je možné využiť pri volani príkazu - tieto parametre upresňujú format súboru a spôsob načítania do R. Viac podrobnosti v príslušnom helpe ?read.table prípadne ?read.csv.

Užitečné

Program R umožňuje priamo načítať aj online zdroje. Potrebné je iba špecifikovať link a použiť vhodny príkaz. Napríklad:
```
test.data <- read.table("http://www.karlin.mff.cuni.cz/~maciak/NMSA407/NMSA407-1617-HW1.txt", header=T)
```

5. Príprava a programovanie vlastných R skriptov

Štatistický program R slúži aj ako výborný programovací jazyk: pomocou tohto programu je možné naprogramovať takmer akúkoľvek úlohu. Je pri tom užitočné poznať rôzne nástroje a konkrétne funkcie, ktoré k tomuto účelu môžu poslúžiť (napr. implementácia cyklov for a while, alebo overovanie podmienok pomocou if).

for cyklus:
```
for (i in 1:10){
  print(1:i)
}
```

while cyklus:

i <- 1
while (i < 10){
  print(1:i)
  i <- i + 1
}

overovanie podmienky pomocou if

N <- rpoiss(1)
if (N %% 2 == 0){
  print("Sude cislo")
} else {
  print("Liche cislo")
}

Pri používaní cyklov je dôležité dôsledne skontrolovať skript - môže sa totiž veľmi jednoducho stáť, že program R neúmyselne zacyklíme. V takom prípade je dobré poznať klávesovú skratku CTR + C (v konzole programu R), prípadne ESC v interface RStudio.

Užitočné

Jednotlivé príkazy je možné v programe R vzájomne kombinovať a vytvárať z nich komplexné Rkové skripty. K dispozícii je samozrejme mnoho ďalších nástrojov určených k pokročilému programovaniu (viď napr. help a rôzne tutoriály).

II. odhady charakteristik rozdělení

Primárně je program R určený pro statistickou analýzu dat. Každá analýza libovolného datového souboru by měla pozostávat z dvou časti: popisná část a analyticka čast. V popisnej časti se pomoci popisnych charakteristik (průměr, výberový rozptyl, vyběrový median, a pod.) a obrázků/grafů snážíme získat dostatečne kvalifikovaný pohled na data - jakou maji struktúru, co s čím vzájemne souvisí a pod. Na základne této popisné části pak volíme vhodný pravděpodobnostný/statistický model, kterým data formálne analyzujem (napr. dvouvýběrový t-test, alebo analýza rozptulu, nebo regrese, a pod.).

1. Výpočet základních charakteristik náhodného výběru

## Výběrový průměr
mean(x)
## výběrový medián
median(x)
## výběrový rozptyl
var(x)

## minimum
min(x)
## maximum
max(x)

## summary
summary(x)

## pořádkové statistiky
sort(x)
## 22. pořádková statistika
sort(x)[22]
## vektor pořadí
rank(x)

2. Tabulace náhodného výběru, empirické relativní četnosti

## najdi, v kterých intervalech pozorování jsou
tbl=cut(x,breaks=c(-Inf,-3:3,Inf))

## kolik pozorování je v kterém intervalu
table(tbl)
## empirické relativní četnosti
table(tbl)/length(x)

3. Grafické zobrazování náhodného výběru

Boxplot (krabicový diagram) zobrazuje krabici ohraničenou horním a dolním výběrovým kvartilem, výběrový medián tlustou vodorovnou čárou uvnitř krabice, vousy, které jdou do jedenapůlnásobku délky boxu (mezikvartilového rozmezí) nahoru i dolů. Pozorování, která se nevejdou do vousů, jsou nakreslena zvlášť (odlehlá pozorování/outliers).

boxplot(x)

Stejný boxplot, ale s väčším dôrazom na estetickú formu.

boxplot(x, col = "lightblue", xlab = "Popis pro os x", ylab = "Popis pro os y", main = "Boxplot pro X")

Histogram ukazuje absolutní nebo relativní četnosti pozorování v určitých intervalech (lze měnit jejich počet i umístění). Histogram s relativními četnosti odhaduje hustotu pomocí po částech konstantní funkce. Histogram ale není konsistentním odhadem hustoty.

## histogram s absolutními četnostmi
hist(x, col = "lightblue")

Stanovíme si rozsah osy x, namaluje e reativní četnosti místo absolutních a proložíme teoretickou hustotu rozdělení N(0,1).

## Stanovíme rozsah osy x na (-3.5,3.5), 
## malujeme relativní četnosti místo absolutních
hist(x,freq=FALSE,xlim=c(-3.5,3.5), col = "lightblue")

## Spočítej a domaluj hustotu N(0,1) na intervalu (-3.5,3.5)
xpts=seq(-3.5,3.5,length=500)
lines(xpts,dnorm(xpts,mean=0,sd=1), col = "red", lwd = 2)

Nyní zvětšíme rozsah výběru na n=5000.

x=rnorm(n=5000,mean=0,sd=1)

hist(x,freq=FALSE,xlim=c(-3.5,3.5), ylim = c(0, 0.6), col = "lightblue")
xpts=seq(-3.5,3.5,length=500)
lines(xpts,dnorm(xpts,mean=0,sd=1), col = "red")

Podobá se tvar histogramu hustotě rozdělení, z kterého data pocházejí?

Ešte raz zväčšíme rozsah vyberu na n=50000 a zjemníme historam.

x=rnorm(n=50000,mean=0,sd=1)

hist(x,freq=FALSE,xlim=c(-3.5,3.5), ylim = c(0, 0.6), breaks = 40, col = "lightblue")
xpts=seq(-3.5,3.5,length=500)
lines(xpts,dnorm(xpts,mean=0,sd=1), col = "red")

Empirická distribuční funkce je nestranným a konsistentním odhadem skutečné distribuční funkce (viz větník a poznámky z přednášek).

\[ \widehat{F}_n(u)=\frac{1}{n}\sum_{i=1}^n\mathbb{I}_{(-\infty,u\rangle}(X_i) \]

Ukažme si empirickou distribuční funkci pro výběr z N(0,1) o rozsahu n=30.

## Rozsah výběru 30
x=rnorm(n=30,mean=0,sd=1)

## Výpočet a malování empirické distribuční funkce
plot(ecdf(x))

## Stanovíme rozsah osy x na (-3.5,3.5), zmenšíme puntíky, přepíšeme titulek
plot(ecdf(x),xlim=c(-3.5,3.5),cex.points=0.5, 
     main="Empiricka distribucni funkce N(0,1), n=30")

## Spočítej a namaluj skutečnou d.f. N(0,1) na intervalu (-3.5,3.5)
xpts=seq(-3.5,3.5,length=500)
lines(xpts,pnorm(xpts,mean=0,sd=1), col = "red")

Nyní zvětšíme rozsah výběru na n=200.

## Rozsah výběru 200
x=rnorm(n=200,mean=0,sd=1)

## Empirická distribuční funkce
plot(ecdf(x),xlim=c(-3.5,3.5),cex.points=0.1, 
     main="Empiricka distribucni funkce N(0,1), n=200")

## Spočítej a namaluj d.f. N(0,1) na intervalu (-3.5,3.5)
xpts=seq(-3.5,3.5,length=500)
lines(xpts,pnorm(xpts,mean=0,sd=1), col = "red")

Podobá se empirická distribuční funkce distribuční funkci, z které data pocházejí?

4. Nestrannost a konsistence bodových odhadů

Ilustrujme si pojmy nestrannosti a konsistence. Budeme pracovat s výběrovým rozptylem na simulovaných datech. Máme

\[ S^2_n=\frac{1}{n-1}\sum_{i=1}^n (X_i-\overline{X}_{n})^2 \]

Nestrannost výběrového rozptylu

## Spočítej 1000x výběrový rozptyl výběru o rozsahu 100
nsim=1000
nobs=100
xmat=rnorm(n=nsim*nobs,mean=0,sd=1)
xmat=matrix(xmat,nrow=nobs,ncol=nsim)

dim(xmat)

## [1]  100 1000

## V i-tém řádku je 1000 opakování pozorování X_i.
## V j-tém sloupci je náhodný výběr o rozsahu 100.
## Z každého výběru spočítáme výběrový rozptyl, bude jich 1000.
rozptyly=apply(xmat,2,var) 

## Podívejte se na prvních 50 odhadů rozptylu.
## Podobají se skutečnému rozptylu?
rozptyly[1:50]

##  [1] 1.0509725 1.0816750 0.9677750 0.9558585 0.8989324 1.0163461 0.9378098
##  [8] 1.0201024 0.9408067 0.7454111 1.1168161 0.8532366 1.0300669 1.0941115
## [15] 0.9098931 0.9240618 0.8768127 1.2364632 0.8937028 0.8620224 0.8911758
## [22] 0.9946921 0.7888259 1.1455473 1.2220950 1.0283665 1.0628913 1.0256976
## [29] 1.0790981 1.0966200 0.8489592 1.1657755 1.1143616 0.9035604 1.0481938
## [36] 1.0126741 0.9682610 0.9762675 0.9424606 0.8422652 0.8623573 0.9168330
## [43] 1.1372099 0.8371186 0.8534834 0.8095092 1.1058609 1.0576522 0.9421184
## [50] 1.0370796

## Jaký je průměr z 1000 odhadů rozptylu?
mean(rozptyly)

## [1] 0.9969071

Odpovídá tento výsledek tomu, co bychom očekávali od nestranného odhadu?

Konsistence výběrového rozptylu

## Spočítej výběrový rozptyl z výběru o rozsahu 10, 100, 200, 1000
rozsahy=c(10,100,1000,10000, 100000)

## Sem dáme výsledky
vybrozp=rep(0,length(rozsahy))
names(vybrozp)=paste("n",rozsahy,sep="=")

x=rnorm(max(rozsahy),mean=0,sd=1)

for(n in rozsahy)
{
  vybrozp[paste("n",n,sep="=")]=var(x[1:n])
}

vybrozp

##      n=10     n=100    n=1000   n=10000   n=1e+05 
## 1.7384124 1.1655157 1.0873766 1.0077570 0.9991455

Tohle lze zobrazit graficky napr. pomocou nasledujúceho obrázku.

plot(0,0, xlim = c(1,length(rozsahy)), ylim = c(0.5,1.5), pch = "", xlab = "Logaritmus rozsahu výběru", ylab = "Vyběrový rozptyl")
lines(vybrozp ~ seq(1,length(rozsahy)), col = "red")
abline(1,0, col = "black", lty = 2)

Užitočné

Pomocou helpu naštudujte, k čomu slúžia funcie apply(), tapply(), aggregate() a by();
Použijte tieto funkcie v súvislosti s datovým súborom ‘test.data’ a pomocou nich spočítajte niektoré empirické charakteristiky;
V pripade, že pracujeme v programe R len s jednym datovým súborom, je možné prácu (resp. implementáciu príkazov) mierne zjednodušiť - pomocou príkazu attach(test.data) program R zprístupni jednotlivé premenné, ktoré je potom možne volať priamo pomocou ich názvu:
```
attach(test.data)
mean(waitingTime[driverGender == "male"])
mean(waitingTime[driverGender == "female"])
```
V prípade, že takto zprístupnene premenné už neplánujeme využívať, je dobre využiť komplementárny príkaz detach().
Je dôležité správne pochopiť fungovanie Rka a odkazovanie na premenné po použití príkazu attach(). Program R totiž zprístupni premenne v stave, v akom sa v momente volania príkazu attach() aktuálne nachádzaju v data.frame. Akýkoľvek následný zásah do data.framu a zmena hodnôt v ňom, nebude reflektovaná ak sa na premennu odkážeme pomocou jej názvu a nie priamo cez odkaz na data.frame;

Doplňující úloha 1.

Zopakujte vše pro rozdělení Exp(2) namísto N(0,1)

## Generování veličin z Exp(2): 
rexp(n,2)
## Hustota Exp(2): 
dexp(x,2)
## Distribuční funkce Exp(2): 
pexp(x,2)

Doplňující úloha 2.

Zopakujte vše pro Cauchyovo rozdělení C(0,1) namísto N(0,1)

## Generování veličin z C(0,1): 
rcauchy(n,0,1)
## Hustota C(0,1): 
dcauchy(x,0,1)
## Distribuční funkce C(0,1): 
pcauchy(x,0,1)

NMFM301, ZS 2018/19

Cvičenie 7 (praktické cvičenie 1)

(Základy práce so statistickým softwarom R)

I. Štatistický software R (stručný úvod)

1. Základné matematické operácie v R

Užitočné

2. Práca s vektormi a maticami

Užitočné

3. Náhodné data - generátory náhodných čísel

Užitočné

4. Reálne data - načítanie vlastných datových súborov

Užitečné

5. Príprava a programovanie vlastných R skriptov

Užitočné

II. odhady charakteristik rozdělení

1. Výpočet základních charakteristik náhodného výběru

2. Tabulace náhodného výběru, empirické relativní četnosti

3. Grafické zobrazování náhodného výběru

4. Nestrannost a konsistence bodových odhadů

Nestrannost výběrového rozptylu

Konsistence výběrového rozptylu

Užitočné

Doplňující úloha 1.

Doplňující úloha 2.