NMFM301, ZS 2021/22

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

(Základy práce so štatistickým programom R)

Zdrojový soubor pro Knit: Rmd soubor



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

Štatistický softwar R (R Development Core Team, 2019) je výpočetný nástroj určený k spracovaniu, príprave a následnej štatistickej analýze dat. Jedná sa o GNU projekt, podobný programu S a primárne je navrhnutý pre praktickú aplikáciu rôznych štatistických metód.

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 najčastejšie používaných R interfacov je RStudio.

Užitočné doplňujúce materiály pre prácu so štatistickým programom 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úcich riadkoch je uvedených niekoľko jednoduchých príkladov a úloh, ktoré majú primárne slúžiť ako názorné ukážky fungovania programu R (implementácia jednotlivých príkazov a práca s vlasnými datovými súbormi). Nejedná sa ale o kompletný súhrn možnosti, ktoré sú v programe R k dispozícii.

V prípade potreby 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

Najjednoduchšie využitie programu R ako jednoduchá kalkulačka:

2 + (4 * 5)^2 ## za '#' nasleduje komentar, ktory R ignoruje
sqrt(256)
1:5 * 2
1 + 2 * (2 * (2 - 1) + 2)

Výsledky je možné 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
print(v1)

Uložené objekty (výsledky, alebo aj komplexnejšie datové štruktúry ako vektor, matica, pole, list) je možné vhodné kompbinovať a ďalej s nimi pracovať.

v1 + v2
v3 <- log(v1 + v2)
v4 <- v3 - 2
v5 <- sin(pi * v2)



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, ktoré majú rôznu dĺžku, operácie s maticami, ktoré majú rôzne typy/dimenzie. Je podstatné správne pochopiť, ako program 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.

  • Štandardné operácie v programe R získame pomocou štandardných znamienok +, -, *, \, a symbolu ^ pro mocniny;
  • Komplexnejšie matematické výrazy môžeme formulovať pomocou zátvoriek. V Rku slúžia k tomuto účelu pouze okrúhle zátvorky (). Hranaté zátvorky [] a zátvorky {} maju v programe R vlastnú, špecifickú úlohu;
  • Ku každému príkazu v R je k dispozícii návod, ako príkaz správne funguje - help sa zavolá zadaním samotného príkazu, ktorému predchádza otazník: napr. ?c zobrazí help pre funkciu c(), analogicky ?seq zobrazí help k funkcii seq(). Help obsahuje popis samotnej funkcie, vysvetlenie implementácie a tiež niekoľko ilustratívnych príkladov a odkazov na dalšie príbuzne funkcie, alebo funkcie podobného charakteru.
  • Podrobný návod je k dispozícii aj na internete a v prípade problému stačí použiť google a riešenie sa určite ľahko nájde;


2. Práca s vektormi a maticami

V prípade, že máme v objektoch uložené nejaké hodnoty (vektor, maticu, pole, list, atď.), je možné sa v Rku odkazovať aj na jednotlivé elementy v daných objektoch, používať ich v dalších výpočtoch, alebo s nimi manipulovať. K tomuto účelu práve slúžia hranaté zátvorky [].

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 

Objekty rôzneho typu môžu byť pohromade uložené napr. v objekte typu list()

li <- list(v3, m1) ## list, ktory obsahuje vektor v3 a maticu m1

Na jednotlivé elementy v liste sa odkazuje následovne:

li[[1]] ## zobrazi prvy element listu, vektor v3
li[[2]] ## zobrazi druhy element listu, maticu m1

Vždy je vhodné voliť kompaktnú a čo najjednoduchšiu formu zápisu. Pre porovnanie, všetky nasledujúce príkazy vytvoria vo výsledku jednotkovú maticu typu 4x4, ale najkompaktnejši a najstručnejší je pouze posledný zápis. Vyžaduje ale znalosť príkazu diag() a jeho spávnu implementáciu.

m1 <- rbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))

m2 <- cbind(c(1,0,0,0), c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))

m3 <- matrix(rep(0,16), nrow = 4)
diag(m3) <- 1

m4 <- NULL
for (i in 1:4){
  m4 <- rbind(m4, c(rep(0, i - 1), 1, rep(0, 4 - i)))
}

m5 <- diag(1, 4, 4)



Užitočné


  • V programe R existuje niekoľko prikazov, ktoré vytvoria vektor: napr. c(), seq(), rep(), replicate() a pod.;
  • Podobne existuje niekoľko príkazov, ktoré vytvoria maticu: napr. cbind(), rbind(), matrix(), diag(), upper.tri(), lower.tri(), data.frame() a pod.;
  • Užitočné príkazy pre vytváranie a pracovanie s vektormi a maticami sú as.numeric(), as.matrix(), array().
  • Pre prácu s objektom typu list() sú užitočné príkazy unlist(), as.list(), pairlist() a pod.
  • Pre maticové násobenie je nutné použiť 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(), ncol(), alebo dimnames()).
  • Inverznú maticu lze získať príkazom solve(), vlastne čísla/vektory matice pomocou príkazu eigen(), rozklady pomocou príkazov eiv(), svd(), qr(), alebo chol().
  • Š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é sú rôzne generátory (pseudo) 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.

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:

    x <- rnorm(10, 0, 2) ## simuluje 10 pozorovani z N(0, 4) - sigma^2 = 4
    y <- rexp(10, 4) ## simuluje 10 pozorovani z Exp(lambda = 4), kde EX = 1/4
  • V prípade použítia nejakého generátora náhodných čísel je vhodné predom nastaviť tzv. seed. Umožní to následnu replikáciu generovania s garanciou rovnakých výsledkov.

    set.seed(1234)
  • Výsledok použitia generátora náhodných čisel môže byť napríklad vektor, alebo matica (každý stĺpec ma stejný počet hodnôt), alebo list (umožňuje rôzny počet hodnôt v jednotlivých elemntoch).
  • set.seed(1234)
    g1 <- rnorm(20, 160, 20) ## napr. vyska 20 nahodnych lidi s prumernou vyskou 160 cm a rozptylom 400
    
    g2 <- cbind(g1, rexp(20, 1), rbinom(20, size = 1, prob = 1/2)) ## matica
    
    g3 <- list(g1, rexp(10, 1), rbinom(50, size = 1, prob = 1/2)) ## list


4. Reálne data a vlastné datové súbory

V programe R je k dispozícii niekoľko vzorových ilustračných datových súborov. Ich zoznam je možne zobraziť zavolaním príkazu data(). Konkrétny datovy súbor je následne zobrazený zavolaním mena prislušného súboru z príkazoveho riadku v R - napr. Orange, alebo mtcars. Vo väčšíne prípadov sa jedná o objekt typu data.frame - resp. špeciálny typ matice (pozri help ?data.frame). Na rozdiel od klasickej matice umožňuje type data.frame()použiť rôzne datové typy (integer, boolean, character, etc.) v jednotlivých stĺpcoch - premenných. Je tak možné vytvoriť datový súbor, ktorý bude obsahovať premenné rôznych typov - spojité veličiny, disrétne hodnoty, ale aj faktory, alebo štandardný popis vo forme textu.

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)

S ľubovolným datovým súborom sa následne pracuje ako so štandardnou maticou - môžeme sa odkazovať ne jednotlivé elementy, celé stĺpce (väčšinou jednotlivé premenné), alebo konkrétne řádky (väčšinou jednotlivé pozorovania). Nasledujúce príkazy dávajú stejné výstupy, spôsob implementácie príkazu je ale zakaždým iný. ´

head(Orange) ### zaujima nas pouze promena circumference (treti sloupec)

attach(Orange)
circumference

Orange$circumference

Orange[,3]

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.

Analogicky je možné stejne súbory pomocou programu R uložiť a zapísať konkrétny subor na disk.

Užitoč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)
  • Pomocou príkazov dim(), str(), head(), alebo summary() preštudujte štruktúru datového súboru a pomocou vhodných indexov, názvu premenných, alebo pomocou príkazov subset() a select() sa pokúste definovať rôzne podsúbory (v zmysle pozorovaní–riadkov, aj v zmysle premenných–stĺpcov). Napríklad:

    ### pouze vzdalenosti do 50 km
    data50_1 <- test.data[test.data$travelDistance == "< 50", ] 
    data50_2 <- test.data[test.data[,2] == "< 50", ] 
    data50_3 <- subset(test.data, travelDistance == "< 50")
    
    ### jedna sa o tri ekvivalentne datove subory
    all.equal(data50_1, data50_2)
    ## [1] TRUE
    all.equal(data50_1, data50_3)
    ## [1] TRUE
    ### nebo pouze data, kde stopovali vzdy dvaja a trip byl v CR
    dataCR2_1 <- test.data[test.data$hitchhikers == "two guys" & test.data$country == "CZ", ]
    dataCR2_2 <- test.data[test.data[,1] == "two guys" & test.data[,4] == "CZ", ]
    dataCR2_3 <- subset(test.data, hitchhikers == "two guys" & country == "CZ")
    
    ### opat sa jedna o tri ekvivalentne datove podsubory
    all.equal(dataCR2_1, dataCR2_2)
    ## [1] TRUE
    all.equal(dataCR2_1, dataCR2_3)
    ## [1] TRUE
    ### vyberieme data, kde krajina je CZ nebo "other", stopuju pouze dva a zaujima nas  waiting time a vzdalenost
    dataCZother_1 <- test.data[(test.data$country == "CZ" | test.data$country == "other") & test.data$hitchhikers == "two guys", c(2,5)]
    dataCZother_2 <- test.data[(test.data$country == "CZ" | test.data$country == "other") & test.data$hitchhikers == "two guys", -c(1,3,4,6,7)]
    dataCZother_3 <- test.data[(test.data$country == "CZ" | test.data$country == "other") & test.data$hitchhikers == "two guys", names(test.data) %in% c("travelDistance", "waitingTime")]
    dataCZother_4 <- test.data[(test.data[,4] == "CZ" | test.data[,4] == "other") & test.data[,1] == "two guys", names(test.data) %in% c("travelDistance", "waitingTime")]
    dataCZother_5 <- subset(test.data, (country == "CZ" | country == "other") & hitchhikers == "two guys", select = c("travelDistance", "waitingTime"))
    
    ### opat sa jedna o ekvivalentne datove subory
    all.equal(dataCZother_1, dataCZother_2)
    ## [1] TRUE
    all.equal(dataCZother_1, dataCZother_3)
    ## [1] TRUE
    all.equal(dataCZother_1, dataCZother_4)
    ## [1] TRUE
    all.equal(dataCZother_1, dataCZother_5)
    ## [1] TRUE
  • Použijte google a help v programe R a pokúste sa načítať do prostredia R napr. obrázok (jpg súbor). Podívajte sa, čo ste do Rka načítali (resp. aký objekt je v Rku po načítaní uložený)


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

Štatistický program R je možné plnohodnotne využiť aj ako štandardný 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ť (najdôležitejšie sú asi 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 <- rpois(1, lambda = 10)
    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 různych charakteristik rozdělení

Primárne 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ájemně souvisí a pod. Na základne popisné části pak volíme vhodný pravděpodobnostní/statistický model, pomocou kterého data formálne analyzujeme (napr. dvouvýběrový t-test, analýza rozptylu, nebo model lineární egrese, a pod.).

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

Budeme uvažovať náhodný výběr z normálního rozdělení, který jsme vygenerovali v časti 3 a který máme v programe R uložený jako objekt x.

## nagenerovane hodnoty nehodneho vyberu
print(x)
## 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 (prezentovaný grafický výstup by nikdy neměl byt vytvořen pouze defaultním volaním nějaké funkce). Dobrý obrázek musí obsahovat všechny podstatné informace: popisky na na oboch osích (ve stejném jazyku jako je text celého dokumentu) a případně informaci o použitých jednotkách.

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ů

Pre nagenerovaný náhodný výber (uložený v objekte x) si spočítame jednu konkrétnu výberovú charakteristiku (napr. výberový rozptyl) a pomocou simulácii a obrázkov si ilustrujme pojmy nestrannosť a konsistence. Pre výberový rozptyl platí nasledujúci vzťah:

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

a) 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.8675656 0.9122906 1.0901116 0.9604818 1.1001118 1.1261891 0.9087724
##  [8] 0.8608041 1.1313290 0.9360728 0.9651728 0.9528824 0.9924122 1.1569610
## [15] 0.8494274 0.9810858 0.8821455 0.9446420 0.9454784 0.9861092 0.8545362
## [22] 0.8735921 0.9794742 0.8471861 1.0901251 0.9875349 0.9221278 1.1470627
## [29] 1.0365498 0.9000493 1.0818617 1.0283511 0.8801763 0.7772892 1.2684647
## [36] 1.0117876 0.9150579 0.8866460 0.5919318 0.8918105 0.9026967 0.8810698
## [43] 0.9929394 0.9541029 0.9023404 0.9886412 1.4267272 1.0943251 1.1008315
## [50] 1.0211727
## Jaký je průměr z 1000 odhadů rozptylu?
mean(rozptyly)
## [1] 0.9936741

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


b) 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.3015907 1.3448177 0.9984664 1.0257253 1.0054709

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úca samostatná úloha

Zopakujte Časť II pre náhodný výber generovaný z exponenciálneho rozdelenia Exp(2) namiesto pôvodného štandardného normálneho rozdelenia N(0,1) a potom z Cauchyho rozdelenia C(0,1). Zamerajte sa na nestrannosť a konzistenci výberového rozptylu (alebo výberovej strednej hodnoty, čo si vyžaduje trochu výraznejšie prispôsobenie vyššie uvedeného Rkového kódu)-

## Generování veličin z Exp(2): 
rexp(n,2)
## Hustota Exp(2): 
dexp(x,2)
## Distribuční funkce Exp(2): 
pexp(x,2)
## 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)