MAI061 cvičení 2
Z ωικι.matfyz.cz
###########################################################
##### 2. CVICENIE - 14.03.2007 ######
###########################################################
# Podrobnosti k cviceniu na mojej web stranke, na adrese:
# http://www.mmatthew.matfyz.cz/school.php
# Podrobnosti k prednaske na stranke Dr. Hlavku, na adrese:
# http://www.karlin.mff.cuni.cz/~hlavka
###########################################################
## ROZDELENIA, NAHODNY VYBER, ZAKLADNE CHARAKTERISTIKY ##
###########################################################
# # Normalni rozdeleni
rnorm(100) # generuje nahodny vyber o rozsahu 100 z N(0,1)
# hustota normalniho rozdeleni N(2,1) v bode 2
dnorm(2, mean = 2, sd = 1)
# Hodnota distribucni funkce rozdeleni N(2,1) v bode 2
pnorm(2, mean = 2, sd = 1)
# kvantil rozdeleni N(2,1) v bode 2
qnorm(2, mean = 2, sd = 1)
# # dalsi rozdeleni
# (rdpq)unif ... rovnomerne
# (rdpq)exp ... exponencialni
# (rdpq)t .... t-rozdeleni
# (rdpq)cauchy ... cauchyho rozdeleni
# (rdpq)f .... F-rozdeleni
# (rdpq)binom ... binomicke rozdeleni
# (rdpq)pois ... poissonovo rozdeleni
# Modelujme si rozdeleni IQ v populaci studentu MFF.
# O IQ se predpoklada, ze ma normalni rozdeleni
# se stredni hodnotou 115 a rozptylem 64. Pocet
# studentu MFF je 2000.
# # Normalni rozdeleni
mff <- rnorm(n = 2000, mean = 115, sd = 8);
# Pokud jsme to nepopletli, tak prumer by mel byt kolem 115
mean(mff)
# a rozptyl kolem 64
var(mff)
# Histogram by mel pripominat normalni rozdeleni
hist(mff, freq = FALSE);
# Prolozime teoretickou hustotou..
curve(dnorm(x,mean = 115, sd = 8), col="blue", add=TRUE);
# Pokud bychom neznali prikaz curve, mohli bychom
# postupovat nasledovne
# sit hodnot x, pro ktera chceme pocitat hustotou
grid <- seq(80,150,by= 0.1);
lines(x=grid, y = dnorm(grid,mean = 115, sd = 8), col="blue");
# # Otazka: Jake procento studentu MFF ma IQ vyssi nez 130?
# Simulacni reseni:
# 1.Generujeme si dostatecne velky vyber
N <- 1000000;
x <- rnorm(N, mean = 115, sd = 8);
# Podil studentu s IQ vyssim nez 130 pak spocteme jako
sum(x > 130)/N;
# 2.Teoreticke reseni
1-pnorm(130, mean = 115, sd = 8);
# Kazdy rok nahodne vybereme 100 studentu za ucelem
# odhadnuti prumerneho IQ (tj. chceme odhadnout
# teoreticke prumerne IQ, ktere je 115 a ktere nyni nezname).
# Je lepe pouzivat prumer ci median vyberu?
# Jak kvalitni jsou tyto odhady? Tj. jak dobre odhaduji 115?
# pocet opakovani
opak <- 1000;
# velikost vyberu
n <- 50
prumery <- rep(0, opak);
mediany <- rep(0, opak);
for(i in 1:opak){
# generujeme si nahodny vyber
x <- rnorm(50, mean = 115, sd = 8);
prumery[i] <- mean(x);
mediany[i] <- median(x);
}
# Ktery odhad se Vam zda vyhodnejsi pouzivat?
# Odhaduji tyto v prumeru odhady to, co maji odhadovat?
mean(prumery);
mean(mediany);
# A ktery z odhadu ma mensi rozptyl
var(...)
# Muzeme si take zkusit nakreslity histogramy
hist(prumery, freq = FALSE, breaks = 20)
hist(mediany, freq = FALSE, breaks = 20, add = TRUE, border = "blue")
# Prehlednejsi je ale asi si nakreslit "vyhlazene histogramy",
# tj. "jadrove odhady"
plot(density(prumery), main = "");
lines(density(mediany), col = "blue");
# # Porovnejme jeste studentyu z MFF a z CVUT.
# O rozdeleni IQ studentu z CVUT prepokladejme,
# ze je normalni se stredni hodnotu 110 a rozptylem 10.
curve(dnorm(x,mean=115,sd=8),from=80,140, lty=1, main="Hustoty rozdeleni IQ", xlab="IQ", ylab="Hustota f(x)")
curve(dnorm(x,mean=110,sd=10),lty=2,add=TRUE, col = "blue")
# a udelame legendu, viz help(legend)
legend(x=85, y=0.04, c("mff","CVUT"), lty=c(1,2), bty="o", col = c("black", "blue"))
# Komu by se nelibily pouzite barvy, tak ma nasledujici vyber
farby <- colors()
farby;
### napr "modre" farby (nazvy z anglictiny):
farby[grep("blue",farby,ignore.case=TRUE)]
### -hranate zatvorky - indexovanie pola
# # Bonus
# Jaka je pst, ze kdyz se potka se potka student z MFF
# se studentem z CVUT, tak matfyzak bude mit vyssi IQ?
# # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # Alternativni rozdeleni
# # # # # # # # # # # # # # # # # # # # # # # # # # #
# Q: Hazime (nezfalsovanou) minci. Jaka je pst, ze pocet
# rubu v 1000 pokusech bude mezi 450 a 550?
# A: Jedna se o binomicke rozdeleni, kdy pocet pokusu
# je n = 1000 a pst uspechu je p = 0.5. Tudiz hledanou
# pst spocteme analyticky jako
pbinom(550, size = 1000, prob = 0.5) - pbinom(449, size = 1000, prob = 0.5)
# Ze zakona velkych cisel vime, ze podil rubu by mel
# konvergovat k jedne polovine.
# Generujme si 10 000 hodu minci
x <- rbinom(10000, size = 1, prob = 0.5);
# Spocteme si kumulativni pocty uspechu...
x <- cumsum(x);
# ...a podily uspechu
podil.uspechu <- x/seq(1,10000,1);
plot(podil.uspechu[-(1:50)], type = "l")
abline(h = 0.5, lty = "dotted")
# Pokud bychom si chteli dany obrazek ulozit,
# tak bud se proklikame (coz pod Linuxem nejde),
# anebo pouzijeme
dev.copy2eps(file="konvergence.eps")
# pokud preferujeme jiny format, napr. "pdf" pak
# muzeme pouzit
dev.print(pdf, file = "konvergence.pdf")
# Dalsi prikazy pro graficka zarizeni
dev.list() # zoznam vsetkych graf. zariadeni
dev.off() # ukonci aktualne graficke zariadenie
graphics.off() # ukonci vsetky graficke zariadenia
# # Bonus: 200 studentu jde na zkousku z analyzy.
# Kazdy z nich uspeje s pravdepodobnosti 0.23.
# a) Jaka je pst, ze zkousku udela mene nez 50 studentu?
# b) Predpokladejme, ze pst udelani zkousku zavisi na
# IQ studenta pomoci vzorce p = (IQ - 70)/200.
# IQ studentu ma normalni rozdeleni N(115,64).
# Jaka je pst, ze zkousku udela mene nez 50 studentu?
###########################################################
##### PRIKLAD (BONUS) ######
###########################################################
### Studentovo t-rozdelenie
curve(dnorm(x,0,1),-5,5,lwd=2,col="red",main="Hustoty N(0,1), Studentova a Cauchyho rozdeleni", xlab="x", ylab="Hustota f(x)")
curve(dt(x,15),lwd=2,col="blue",add=TRUE)
curve(dt(x,5),lwd=2,col="green",add=TRUE)
curve(dcauchy(x),lwd=2,col="brown",add=TRUE) # odpovida t(1)
legend(x=-4, y=.3, legend=c("N(0,1)","t(15)","t(5)","Cauchy"), lwd=2, col=c("red","blue","green","brown"), bty="n")
### Chi^2-rozdelenie
curve(dchisq(x,df=2),from=0.0001,to=15,lwd=2,col="blue",main=expression(Hustoty~~chi[(df)]^{2}~~rozdeleni), xlab="x",ylab="Hustota f(x)")
curve(dchisq(x,df=3),lwd=2,col="green",add=TRUE)
curve(dchisq(x,df=5),lwd=2,col="brown",add=TRUE)
curve(dchisq(x,df=10),lwd=2,col="black",add=TRUE)
curve(dchisq(x,df=1),lwd=2,col="red",add=TRUE)
legend(x=8, y=.4, legend=c("df = 1","df = 2","df = 3","df = 5","df = 10"), lwd=2, col=c("red","blue","green","brown","black"), bty="n")
### F-rozdeleni - skusime v novom okne (pozri "help("X11")")
X11() # nove graficke okno...
curve(df(x,20,25),0.0001,10,n=500,lwd=2,col="blue",main=paste("Hustoty ",expression(F(m,n))," rozdeleni"), xlab="x", ylab="Hustota f(x)")
curve(df(x,2,10),0.0001,10,n=500,lwd=2,col="green",add=TRUE)
curve(df(x,10,2),0.0001,10,n=500,lwd=2,col="brown",add=TRUE)
curve(df(x,15,20),0.0001,10,n=500,lwd=2,col="black",add=TRUE)
curve(df(x,1,1),0.0001,10,n=500,lwd=2,col="red",add=TRUE)
legend(x=5, y=.8, legend=c("m = 1, n = 1","m = 2, n=10","m = 10, n = 2","m = 15, n = 20","m = 20, n = 25"), lwd=2, col=c("red","green","brown","black","blue"), bty="n")
dev.cur() # vypise cislo aktualniho okna
dev.set(2) # nastavi aktualni okno na device cislo 2
###########################################################
##### REALNE DATA - ZAKLADNE CHARAKTERISTIKY ######
###########################################################
# pruzkum IQ a znamek
# Pokud data nemate od minula, tak si stahnete
# (do sveho pracovniho adresare) soubor
# data.txt a a nactete tento soubor do R-ka pomoci
DATA <- read.table("data.txt", colClasses=c("factor", "factor", "factor", "numeric", "numeric", "numeric", "numeric"), header=F);
# pojmenuju sloupecky
names(DATA) <- c("index","skola","pohlavi","IQ","znamka1","znamka2","znamka3")
# Prejmenuju pohlavi 0,1 na Z, M
levels(DATA$pohlavi) <- c("Z", "M");
summary(DATA);
# Zajistime si pristup k promennym
attach(DATA);
# Spocteme si korelacni koeficient mezi IQ a znamka1
cor(cbind(IQ, znamka1))
# Bonus
# Zjistete jak se lisi korelacni koeficient mezi IQ a znamkou1
# u divek a kluku zvlast
# # Korelacni koeficient mnohonasobne korelace
cormat <- cor(cbind(IQ,znamka1, znamka3));
#
Rxx <- cormat[2:3,2:3];
Ryx <- cormat[1,2:3];
rho.xy <- sqrt(t(Ryx)%*%solve(Rxx)%*%Ryx);
# Porovnejte rho.xy s prostym korelacnim koeficientem
# cor(IQ, znamka1). Co z toho vyplyva
# # Koeficient parcialni korelace
ryz <- cor(znamka1,znamka3);
ryx <- cor(znamka1,IQ);
rzx <- cor(znamka3,IQ);
# Koeficient parcialni korelace potom je
ryz.x <- (ryz - ryx * rzx)/sqrt((1-ryx^2)*(1-rzx^2));
# porovnejte tuto hodnost s kor. koeficientem IQ a znamka1
