NMFM301, ZS 2016/17

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

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



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)

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