NMFM301, ZS 2015/16
Cvičenie 7 (praktické cvičenie 1)
(Základy práce so statistickým softwarom R)
Základné matematické operácie, jednoduché funkcie
## operace + - * / ^
(5 + 3^2)/7## [1] 2## vektorove funkce c(), rep(), seq()
1:5 + 3## [1] 4 5 6 7 8rep(1, 4) * pi## [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] 4vysledek2[2:3]## [1] 9.869604 31.006277Odhady charakteristik rozdělení
Generování náhodných veličin z rozdělení N(0,1)
x <- rnorm(n=100,mean=0,sd=1)Výpočet základních charakteristik náhodného výběru
## Výběrový průměr
mean(x)## [1] 0.1158581## výběrový medián
median(x)## [1] 0.137608## výběrový rozptyl
var(x)## [1] 1.384921## minimum
min(x)## [1] -2.849593## maximum
max(x)## [1] 2.743918## summary
summary(x)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.8500 -0.7520 0.1376 0.1159 0.9896 2.7440## pořádkové statistiky
sort(x)## [1] -2.84959273 -2.35229475 -2.21621282 -2.18552273 -1.96318340
## [6] -1.88703172 -1.79768193 -1.79157715 -1.45022967 -1.24676190
## [11] -1.21843672 -1.21666663 -1.17379962 -1.09755476 -1.08526652
## [16] -1.07067880 -1.06879250 -0.93210326 -0.90170645 -0.89634344
## [21] -0.82845758 -0.81843549 -0.79943706 -0.76513594 -0.76377083
## [26] -0.74807848 -0.71028095 -0.60561881 -0.60075770 -0.51087803
## [31] -0.49899954 -0.48165786 -0.46892471 -0.43673714 -0.42347463
## [36] -0.29080839 -0.28907381 -0.27967735 -0.27803688 -0.27147523
## [41] -0.24521107 -0.11443179 -0.11034552 -0.10407092 -0.06031791
## [46] -0.05633260 -0.01246909 0.07578795 0.10791224 0.12569792
## [51] 0.14951806 0.19326820 0.22158921 0.22801341 0.23275280
## [56] 0.32242095 0.33213604 0.33428333 0.33773370 0.47778354
## [61] 0.48042316 0.49091362 0.55484748 0.58286549 0.61614864
## [66] 0.63192199 0.68564368 0.68928482 0.70959616 0.73889414
## [71] 0.75695647 0.77940383 0.85031993 0.92944788 0.97691131
## [76] 1.02774267 1.07328137 1.09417837 1.12568787 1.14646590
## [81] 1.16831402 1.19711176 1.23335872 1.27330815 1.28508253
## [86] 1.29258717 1.30992332 1.34361925 1.36573209 1.40048960
## [91] 1.54909434 1.70790283 1.88787226 2.03389788 2.05385980
## [96] 2.05570736 2.39458354 2.50201576 2.68192948 2.74391827## 22. pořádková statistika
sort(x)[22]## [1] -0.8184355## vektor pořadí
rank(x)## [1] 89 42 13 26 1 79 43 95 11 99 88 48 76 23 57 87 31
## [18] 18 22 19 77 37 49 36 92 8 81 50 83 51 73 6 58 17
## [35] 28 98 82 60 32 56 64 21 61 33 15 46 62 100 84 85 2
## [52] 68 47 30 41 91 27 34 7 55 69 25 5 12 63 39 97 3
## [69] 75 29 65 16 67 44 80 71 66 93 70 78 90 54 74 94 86
## [86] 96 9 20 14 4 59 52 45 24 53 10 40 72 38 35Tabulace 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)## tbl
## (-Inf,-3] (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3]
## 0 4 13 30 28 18 7
## (3, Inf]
## 0## empirické relativní četnosti
table(tbl)/length(x)## tbl
## (-Inf,-3] (-3,-2] (-2,-1] (-1,0] (0,1] (1,2] (2,3]
## 0.00 0.04 0.13 0.30 0.28 0.18 0.07
## (3, Inf]
## 0.00Grafické 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 o hodně krajší…
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 me reativní četnosti místo absolutních a proložíme 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.9936585 0.9855783 1.1037764 1.1130100 0.8277200 1.0589077 0.9729988
## [8] 0.7700253 1.0478645 1.0910309 1.0523408 1.1080988 1.0991173 0.8797521
## [15] 1.0092145 1.1230577 0.9787103 0.9605295 0.9691735 1.2000659 0.7766102
## [22] 1.1149297 1.0316581 0.9083414 0.9415115 0.9477570 1.2572628 0.8654387
## [29] 0.9917328 0.7659439 0.7832118 0.7775869 1.1478719 1.0811024 0.9288772
## [36] 0.8621170 0.9092884 1.0651611 0.9199779 1.1866927 0.7135960 0.8467047
## [43] 0.9622396 1.0327821 1.0956775 0.7997090 0.9336769 0.9229877 0.7034714
## [50] 0.9098362## Jaký je průměr z 1000 odhadů rozptylu?
mean(rozptyly)## [1] 0.9996575Odpoví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,200,1000)
## 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=200 n=1000
## 0.6150442 1.1492391 1.0812933 0.9807934Doplň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)