## Hash oznacuje komentar :-) ## R lze volne stahnout z www.r-project.org ## dalsi odkazy na dokumentaci a manualy k R lze nalezt napr. na ## http://www.karlin.mff.cuni.cz/~kulich/vyuka/Rdoc/index.html ## ## R je - interaktivni vypocetni prostredi pro statisticke vypocty a grafiku ## - open-source implementace jazyka S ## ## case sensitive, zakladem je rozhrani prikazove radky - ale existuji ruzna GUI (jednoznacne nejhorsi je default ve Windows) ## ## zakladem R jsou vektory, seznamy a funkce (vse spolu souvisi, nekdy mirne neprehledne) ## ## data v R jsou typicky seznamy vektoru ruzneho typu, ktere navic mohou mit dalsi atributy (cely seznam i jednotlive vektory) ## ## (programovani casto pomoci funkci, ktere automaticky aplikuji jine funkce, ## ktere se mohou skladat z ruznych metod - podle typu objektu a jeho atributu) ## # krok 1: nainstalovani R z www.r-project.org ######## Napoveda # help(q) # napoveda ?q # jinak zavolana stejna napoveda ?help # napoveda k napovede :-) ######## Vektory, matice a zakladni operace # 1+5 # soucet se vytiskne a "zapomene" vysledek <- 1+5 # soucet se zapise do promenne vysledek, ale nevytiskne se vysledek # vypise hodnotu vysledek print(vysledek) # udela totez # # Nebo rychleji: prikazy lze sdruzovat print(vysledek <- 3+7) (vysledek <- 3+7) # udela totez # Zakladni typ objektu je vektor, vytvorime ho prikazem c() - c jako concatenate print(x <- c(2,4,5.4,2.3,3.113,4)) # bezne operace funguji po slozkach s prirozenou prioritou: # Prace s objekty: ls() # vypis seznamu rm("a","vysledek") # vymaz ls() ## typ objektu lze zjistit pomoci funkci class() a mode() class(a) mode(a) # atributy pomoci attributes() attributes(a) # strukturu objektu lze prozkoumat pomoci str() str(a) ## pro data se pouzivaji objekty typu data.frame # data() # vypis dostupnych souboru s defaultnimi daty data(women) # nacteni dat women # data ulozena jako typ data.frame women # tabulka - muze obsahovat vice typu dat soucasne # (numericke, znakove, logicke) str(women) women[,2] # k prvkum lze pristupovat jako pri praci s maticemi, # sikovnejsi je ale vyuzit jmena sloupcu weight # chyba - nezna objekt women$weight # "podobjekty" jsou pristupne pres znak "$", POZOR - ne "." jako nekde jinde :-) women[,"weight"] women[,names(women)=="weight"] attach(women) # zpristupni slozky 1. urovne pod jejich jmeny weight ## vetsina funkci ma metody pro ruzne typu vstupu summary(women) summary(wome$weight) methods(summary) ## jmeno funkce vypise jeji zdrojovy kod: summary ## ## R umoznuje kreslit hezke obrazky: hist(women$height) # histogram boxplot(women$height) ################################################################## ######### Grafy hustot nejpouzivanejsich rozdeleni: ################################################################## ### Exponencialni rozdeleni grid <- seq(from=0,to=15,by=.01) #posloupnost 0 az 15 s krokem 0.01 grid <- 0:1500/100 #jiny zpusob vyjadreni tehoz plot(grid,exp(-grid),type="l",lwd=2,col="red",main="Hustoty Exp(lambda) s ruznou hodnotou E(X)=1/lambda", xlab="x", ylab="Hustota f(x)") lines(grid,0.5*exp(-0.5*grid),lwd=2,col="green") lines(grid,10*exp(-10*grid),lwd=2, col="blue") lines(grid,.1*exp(-.1*grid),lwd=2,col="brown") legend(x=8, y=.8, legend=c("Exp(1)","Exp(0.5)","Exp(10)","Exp(0.1)"), lwd=2, col=c("red","green","blue","brown"), bty="n") # tisk do souboru pomoci presmerovani grafickeho vystupu (prikazy jpeg(), pdf(), ps() s ruznymi parametry, viz napr. help("pdf") nebo ?pdf) # pro kresleni grafu funkci se da vyuzit rovnou funkce curve() curve(dexp(x,rate=1),from=0,to=15,n=151,main="Hustoty Exp(lambda) s ruznou hodnotou E(X)=1/lambda",ylab="Hustota f(x)") curve(dexp(x,0.5),0,15,col="green",add=TRUE) curve(dexp(x,2),0,15,col="violet",add=TRUE) # k dispozici je cela rada barev barvy <- colors() barvy[grep("blue",barvy,ignore.case=TRUE)] # barvy obsahujici nazev 'blue' dev.off() #ukonci soucasne graficke zarizeni ######## Cviceni: ######## Nakreslete obdobne grafy pro distribucni a kvantilovou funkci F(x). Viz help(dexp). ######## Vypocet distribucni a kvantilove funkce. ######## Vypocet stredni hodnoty a rozptylu. # normalni rozdeleni curve(dnorm(x,mean=0,sd=1),-15,15,lty=1,main="Hustoty normalniho rozdeleni s ruznymi parametry", xlab="x", ylab="Hustota f(x)") curve(dnorm(x,0,sqrt(10)),lty=2,add=TRUE) curve(dnorm(x,0,10),lty=3,add=TRUE) curve(dnorm(x,5,1),lty=4,add=TRUE) legend(x=-15, y=.3, c("N(0,1)","N(0,10)","N(0,100)","N(5,1)"), lty=c(1,2,3,4), bty="o") ######## Cviceni: ######## pomoci funkce qnorm() spocitejte "95% predpovedni interval" pro hodnoty z normalniho rozdeleni ######## a pomoci funkce rnorm() vyzkousejte, ze skutecne asi 95% hodnot pada do prislusneho intervalu ######## Cviceni: ######## prozkoumejte "hustotu" diskretniho binomickeho rozdeleni - jak se meni s rostoucim parametrem "size"? n=2 curve(dbinom(x,size=n,prob=0.3),0,n) # pro zajimavost i distribucni funkce curve(pbinom(x,size=n,prob=0.3),0,n,ylim=c(0,1)) #################################################################### ######## generovani odhadu parametru theta pro X ~ R(0,theta) ##### ######## - a porovnavani jejich "statistickych" vlastnosti ##### #################################################################### theta <- 5 # skutecna hodnota parametru opak <- 100 # kolik "nahodnych" vyberu z R(0,theta) budu generovat rozsah <- 50 # pocet pozorovani v jednom vyberu ("n") #rozsah <- 500 mom <- vector() # Momentove odhady (zatim prazdny vektor) mle <- vector() # Max. verohodne odhady (zatim prazdny vektor) ## Ted 100x provedeme "nahodny" vyber z R(0,theta) a 100x odhadneme theta ## momentovou metodou a metodou max. verohodnosti: for (i in 1:opak) { u <- runif(rozsah,0,theta) mom[i] <- 2*mean(u) #### momentovy odhad mle[i] <- max(u) #### max.verohodny odhad } #### nakresli body mle a spoji je modrou carou #### 'expression' vyrobi z retezce vyraz, tj. spec. znaky, #### ten retezec se musi nejdrive zpastovat a az pak pouzit 'expression' #### (bez hacku a carek je to celkem snadne, jinak viz napr. embedFonts()) nadpis <- expression(paste("Odhad ", theta," pro X ~ R(0,",theta,")")) plot(mle, lty=2, type="o", ylim=c(theta-0.9,theta+1.17), col = "blue", pch="+", main = nadpis, xlab="cislo vyberu", ylab="odhad") #### nakresli body momentoveho odhadu a spoji je cervenou prerusovanou carou lines(1:opak,mom,lty=3, type="b", col="red", pch="x") ### skutecna hodnota ################ abline(theta,0, lwd=2) #### nestranny MLE ##### n.mle<-((rozsah + 1)/rozsah)*mle lines(1:opak,n.mle,lty=1, type="b", col="green", pch="*") #### nakresli legendu na souradnice x,y leg.txt <- c(expression(theta),"mom ", "mle","umle") leg.clrs <- c("black","red","blue","green") leg.pchrs <- c(".","x","+","*") leg.ltys <- c(1,3,2,1) legend(x=10, y=6.15, legend = leg.txt, col = leg.clrs, pch = leg.pchrs, lty = leg.ltys, merge = TRUE, trace = FALSE, horiz = FALSE, ncol = 2, bty="n") ### charakteristiky odhadu # odhad stredni hodnoty (skutecna hodnota je theta=5) mean(mom) mean(n.mle) mean(mle) # skutecna stredni hodnota odhadu je theta*rozsah/(rozsah+1) # odhad rozptylu var(mom) theta^2/(3*rozsah) # skutecna hodnota var(n.mle) theta^2/(rozsah*(rozsah+2)) # skutecna hodnota var(mle) theta^2*rozsah/((rozsah+2)*(rozsah+1)^2) # skutecna hodnota ## ## Odhady jsou pekne, ale v praxi je potreba znat i jejich presnost - proto ## se pocitaji konfidencni intervaly a testuji hypotezy... ##