NMFM301, ZS 2013/14
Praktické cvičení 1
Odhady 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.01273## výběrový medián
median(x)## [1] -0.01943## výběrový rozptyl
var(x)## [1] 0.8556## minimum
min(x)## [1] -2.309## maximum
max(x)## [1] 2.167## pořádkové statistiky
sort(x)## [1] -2.30901 -1.84637 -1.80453 -1.64941 -1.46733 -1.44817 -1.43497
## [8] -1.23285 -1.12030 -1.11952 -1.07592 -1.04727 -1.00805 -1.00170
## [15] -0.94656 -0.94204 -0.87523 -0.86999 -0.82348 -0.81663 -0.80082
## [22] -0.79181 -0.76535 -0.75929 -0.75742 -0.74604 -0.70577 -0.62315
## [29] -0.61297 -0.58055 -0.57370 -0.56953 -0.53707 -0.52240 -0.50836
## [36] -0.50544 -0.43008 -0.38644 -0.32310 -0.32303 -0.31763 -0.30527
## [43] -0.28852 -0.21339 -0.17925 -0.15174 -0.12196 -0.07830 -0.05101
## [50] -0.02386 -0.01501 -0.01413 0.01907 0.09127 0.09418 0.12700
## [57] 0.12894 0.17464 0.19465 0.22875 0.24226 0.28496 0.30535
## [64] 0.31625 0.36083 0.39791 0.40417 0.45466 0.45850 0.50693
## [71] 0.55052 0.57403 0.60910 0.62278 0.62332 0.65510 0.66161
## [78] 0.67543 0.70061 0.70400 0.72218 0.77710 0.82950 0.88918
## [85] 0.93755 1.01318 1.03567 1.05203 1.06347 1.10742 1.23589
## [92] 1.27641 1.39410 1.66098 1.66898 1.68328 1.79483 1.79938
## [99] 1.87362 2.16673## 22. pořádková statistika
sort(x)[22]## [1] -0.7918## vektor pořadí
rank(x)## [1] 73 69 100 86 26 62 77 78 71 85 16 47 93 35 22 12 70
## [18] 55 83 82 46 67 59 50 64 57 48 14 41 33 66 80 92 25
## [35] 96 36 9 29 65 45 49 30 88 24 1 99 28 32 8 31 27
## [52] 72 81 58 95 2 44 38 53 34 3 6 13 94 84 75 42 54
## [69] 97 56 68 76 63 23 52 74 15 89 51 20 91 79 60 61 39
## [86] 19 18 11 37 4 5 10 7 87 17 40 21 90 98 43Tabulace 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 1 13 38 33 14 1
## (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.01 0.13 0.38 0.33 0.14 0.01
## (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)
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)
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))
## 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))
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))
xpts=seq(-3.5,3.5,length=500)
lines(xpts,dnorm(xpts,mean=0,sd=1))
Podobá se tvar histogramu hustotě rozdělení, z kterého data pocházejí?
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))
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))
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] 1.1176 1.0134 0.8345 0.9791 0.9902 1.0245 1.2892 0.9802 1.0766 0.9212
## [11] 0.9574 0.7845 1.0316 0.8180 0.8505 0.9456 0.9156 1.0057 0.9158 0.8355
## [21] 1.0991 1.1130 0.9327 0.9032 0.9488 0.8600 1.0785 0.9008 0.8918 1.0713
## [31] 0.9457 1.2102 0.9083 0.8164 1.0609 1.0235 1.0327 1.0901 0.8878 1.1528
## [41] 0.8474 1.0625 1.0164 0.8794 0.8377 0.9451 0.9451 0.9040 0.9216 0.7081## Jaký je průměr z 1000 odhadů rozptylu?
mean(rozptyly)## [1] 1.003Odpoví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
## 1.453 1.093 1.093 1.020Doplň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)