###########################################################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 rozdelenirnorm(100) # generuje nahodny vyber o rozsahu 100 z N(0,1)hustota normalniho rozdeleni N(2,1) v bode 2dnorm(2, mean = 2, sd = 1) Hodnota distribucni funkce rozdeleni N(2,1) v bode 2pnorm(2, mean = 2, sd = 1) kvantil rozdeleni N(2,1) v bode 2qnorm(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 rozdeleniModelujme si rozdeleni IQ v populaci studentu MFF. O IQ se predpoklada, ze ma normalni rozdelenise stredni hodnotou 115 a rozptylem 64. Pocetstudentu MFF je 2000. # Normalni rozdelenimff <- rnorm(n = 2000, mean = 115, sd = 8);Pokud jsme to nepopletli, tak prumer by mel byt kolem 115mean(mff)a rozptyl kolem 64var(mff)Histogram by mel pripominat normalni rozdelenihist(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 nasledovnesit hodnot x, pro ktera chceme pocitat hustotougrid <- 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 reseni1-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 opakovaniopak <- 1000;velikost vyberu n <- 50prumery <- 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 rozptylvar(...)Muzeme si take zkusit nakreslity histogramyhist(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 vyberfarby <- colors()farby;napr "modre" farby (nazvy z anglictiny):farby[grep("blue",farby,ignore.case=TRUE)] -hranate zatvorky - indexovanie pola# BonusJaka je pst, ze kdyz se potka se potka student z MFFse 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 jakopbinom(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 mincix <- rbinom(10000, size = 1, prob = 0.5);Spocteme si kumulativni pocty uspechu...x <- cumsum(x);...a podily uspechupodil.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 pouzijemedev.copy2eps(file="konvergence.eps")pokud preferujeme jiny format, napr. "pdf" pak muzeme pouzitdev.print(pdf, file = "konvergence.pdf")Dalsi prikazy pro graficka zarizenidev.list() # zoznam vsetkych graf. zariadenidev.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-rozdeleniecurve(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-rozdeleniecurve(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 oknadev.set(2) # nastavi aktualni okno na device cislo 2###########################################################REALNE DATA - ZAKLADNE CHARAKTERISTIKY #################################################################pruzkum IQ a znamekPokud data nemate od minula, tak si stahnete (do sveho pracovniho adresare) soubordata.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 sloupeckynames(DATA) <- c("index","skola","pohlavi","IQ","znamka1","znamka2","znamka3")Prejmenuju pohlavi 0,1 na Z, Mlevels(DATA$pohlavi) <- c("Z", "M");summary(DATA);Zajistime si pristup k promennym attach(DATA);Spocteme si korelacni koeficient mezi IQ a znamka1cor(cbind(IQ, znamka1))BonusZjistete jak se lisi korelacni koeficient mezi IQ a znamkou1 u divek a kluku zvlast# Korelacni koeficient mnohonasobne korelacecormat <- 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 korelaceryz <- cor(znamka1,znamka3);ryx <- cor(znamka1,IQ);rzx <- cor(znamka3,IQ);Koeficient parcialni korelace potom jeryz.x <- (ryz - ryx * rzx)/sqrt((1-ryx^2)*(1-rzx^2));porovnejte tuto hodnost s kor. koeficientem IQ a znamka1
