NMFM301, ZS 2015/16

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

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


I. Štatistický software R (úvod a základy)

Š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)

Základné matematické operácie v R

## operace + - * / ^
(5 + 3^2)/7
## [1] 2
## vektorove funkce c(), rep(), seq()
1:5 + 3
## [1] 4 5 6 7 8
rep(pi, 4)
## [1] 3.141593 3.141593 3.141593 3.141593
(c(1,2,3,4) + c(4,3,2,1))/5
## [1] 1 1 1 1
## ukladani vysledku pomoci '<-' nebo '='
vysledek1 <- 1:5 + 3
vysledek2 <- pi^(1:5)
## indexovani prvku
vysledek1[1]
## [1] 4
vysledek2[2:3]
## [1]  9.869604 31.006277



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?


II. odhady charakteristik rozdělení

Generování náhodných veličin z rozdělení N(0,1)

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ú iný 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


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)

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)

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í?

Nestrannost a konsistence 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] 0.9448921 0.8434147 0.9701962 0.9296819 1.1232725 0.9614167 0.8934768
##  [8] 0.9215996 0.8295457 0.8539167 1.2935939 1.2315135 0.7628708 1.0861577
## [15] 1.0497368 0.9435954 1.2018081 0.8395649 1.1433947 1.0294736 1.0802827
## [22] 1.0303864 0.8708213 1.1421448 1.2494750 0.9998519 1.1734000 0.7910824
## [29] 1.2351750 0.6980485 0.8306335 0.9272069 1.1455639 0.7816000 0.8765481
## [36] 1.0411563 1.2810517 1.1428847 0.9619594 0.9801965 1.0259930 1.2374975
## [43] 1.1617587 0.7163770 1.2732380 0.8696839 1.0217924 0.8831469 0.9594164
## [50] 1.2756419
## Jaký je průměr z 1000 odhadů rozptylu?
mean(rozptyly)
## [1] 0.9991265

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 
## 0.5090051 0.9549169 0.9323312 0.9649256 0.9945406

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)