MAI061 cvičení 2

Z ωικι.matfyz.cz
Verze z 21. 3. 2007, 16:30, kterou vytvořil 195.113.20.80 (diskuse)

(rozdíl) ← Starší verze | zobrazit aktuální verzi (rozdíl) | Novější verze → (rozdíl)
Přejít na: navigace, hledání
###########################################################
#####         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
Citováno z „http://wiki.matfyz.cz/index.php?title=MAI061_cvičení_2&oldid=18969“