Š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
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
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)
## 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
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.
c()
, seq()
, rep()
, repliate()
a pod.;
cbind()
, rbind()
, matrix()
, a pod.;
?c()
, ?seq()
, ?rep()
;
%*%
, 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?
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.
rbinom()
, rpois()
, rgeom()
, a ďalšie;
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)
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ý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)
## 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)
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í?
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 \]
## 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?
## 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)
apply()
, tapply()
, aggregate()
a by()
;
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()
.
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;
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)
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)