########################################################## ######## PRIKLAD na procviceni zakladnich pojmu: ######## ######## generovani X ~ exp(lambda) ######## ########################################################## ## Hustota exponencialniho rozdeleni dexp(..) ## Distribucni funkce exp. rozdeleni pexp(..) ## Kvantilova funkce exp. rozdeleni qexp(..) ## Pseudonahodne veliciny s exp. r. rexp(..) lambda <- 1/3 # stredni hodnota je 3 ### s vyuzitim jeho vnitrni fce ex <- rexp(300,lambda) ### pres R(0,1) ### U <- runif(300) gex <- -(1/lambda)*log(U) ### Porovnani dvou vyberu - muzou pochazet ze stejneho rozdeleni? # popisne statistiky summary(ex) summary(gex) # obrazky plot(ex,gex) ### "nesmysl" ... par(mfrow=c(3,2)) plot(sort(ex),sort(gex)) ### ... musime to seradit abline(0,1, lty=2) # prikresli rovnou caru (parametry pocatek a smernice) ## porovnani obou vyberu pomoci krabicoveho diagramu boxplot(ex, gex, notch=TRUE, names = c("ex.orig", "U->ex"), horizontal=TRUE ) ## porovnani obou vyberu pomoci histogramu (tj. relativnich cetnosti) - odhad hustoty hist(ex, prob=TRUE) # histogram (na ose y relativni cetnosti) lines( density(ex), lty=1) # prikresli odhad do tehoz grafu, plnou carou hist(gex, prob=TRUE) # histogram (na ose y relativni cetnosti) lines( density(gex), lty=1) # prikresli odhad do tehoz grafu, plnou carou ## ## podobne muzeme odhadovat napr. i distribucni funkci (to je jeste trochu jednodussi) ## ## PRIKLAD: odvozeni odhadu distribucni funkce v bode "x", tj. F(x)=P(X<=x) pomoci "relativni cetnosti" ## eF.ex<-ecdf(ex) eF.gex<-ecdf(gex) plot(eF.ex) # emp.distr.fce plot(eF.gex) # emp.distr.fce # A jeste obe empiricke distribucni funkce do jednoho obrazku: par(mfrow=c(1,1)) ### vytvori opet 1-no polni obrazek plot(eF.ex,do.points=FALSE, verticals=TRUE, xlim=c(0,(max(ex,gex)+1)),main="Empiricke a teoreticka distr. fce Exp(1/3)") plot(eF.gex,do.points=FALSE,pch = "x", col.points="red",col.hor = "blue", col.vert= "blue",verticals=TRUE, add=TRUE) lines(0:25,pexp(0:25,lambda),type="l", col="red") legend(x=12, y=.6, legend=c("Interni funkce","Pres rovnomerne","Teoreticka F(x)"), lwd=1, col=c("black","blue","red"), bty="n") # ## jenom pro zajimavost: formalni statisticky text shody dvou rozdeleni... ks.test(ex,gex) #test Kolmogorov-Smirnov (testuje shodu dvou distr. funkci porovnanim jejich odhadu) dev.off() ####################################################################### ######### Intervaly spolehlivosti pro normalni rozdeleni - ilustrace ####################################################################### # n <- 50 # rozsah vyberu opak <- 100 # pocet generovanych vyberu mu <- 0 # stredni hodnota sigma <- 1 # odmocnina z rozptylu data <- matrix(nrow=n,ncol=opak) for (i in 1:opak){ data[,i] <- rnorm(n,mu,sigma) } # end for - mame vyrobenou matici kde kazdy sloupec je jeden vektor (nahodny vyber) o rozsahu 50 z N(0,1) # ##### Stredni hodnota normalniho rozdeleni se ZNAMYM rozptylem # # Napocitame BODOVE ODHADY (prumer kazdeho vyberu): means <- apply(data,2,mean) # dava vektor sloupcovych prumeru: # data = pouzita matice, 2 = aplikuj na sloupce (1 by bylo na radky), mean = aplikuj funkci mean (prumer). # Viz help(apply). Je to rychlejsi nez for cyklus :-). # Napocitame a vykreslime 95% intervaly spolehlivosti (INTERVALOVE ODHADY): alpha <- .05 #alpha <- .1 low.lim <- means - qnorm(1-alpha/2)*sigma/sqrt(n) # dolni mez upp.lim <- means + qnorm(1-alpha/2)*sigma/sqrt(n) # horni mez plot(means, pch=20, ylim=c(-.75,.75), xlab="Cislo vyberu", ylab="Bodovy a intervalovy odhad mu", main=paste(100*(1-alpha),"% oboustranne intervaly spolehlivosti pro mu")) # parametr pch vybira znak pro kresleni, ylim urcuje meze osy Y for (i in 1:opak){ lines(c(i,i),c(low.lim[i],upp.lim[i])) lines(c(i-.5,i+.5),c(low.lim[i],low.lim[i])) lines(c(i-.5,i+.5),c(upp.lim[i],upp.lim[i])) } abline(h=mu,lwd=2,col="red") # skutecna hodnota mu length(upp.lim[upp.limmu]) # pocet intervalu nad skutecnou hodnotou length(low.lim[low.limmu]) # pocet intervalu prekryvajicich mu opak*(1-alpha) # ocekavany pocet intervalu prekryvajicich mu # pravdepodobnost, ze interval spolehlivosti pokryje skutecnou hodnotu je # 1-alpha, proto ze 100 vyberu by v prumeru 100*(1-alpha) intervalu # spolehlivosti melo obsahovat 0 # # PRIKLAD: pomoci centralni limitni vety lze spocitat priblizne rozdeleni "podilu" (a tedy i poctu) # konfidencni intervalu prekryvajicich mu - "Uspech v jednotlivem pokusu" ma alternativni rozdeleni # s pravdepodobnosti uspechu 0.95, staci tedy spocitat stredni hodnotu a rozptyl teto nahodne veliciny... # # # # Provedte cele znovu se stejnymi daty, ale pro alpha = 0.1 nebo 0.01. Jak cekate, ze se # zmeni bodove odhady, intervaly spolehlivosti a frekvence chyb? # ##### Stredni hodnota normalniho rozdeleni s NEZNAMYM rozptylem # sds <- apply(data,2,sd) ### ### Uz jsme pouzili: rozdeleni prumeru standardizovane skutecnou smerodatnou odchylkou je normalni. ### ### Otazka: jake je rozdeleni prumeru standardizovane odhadem smerodatne odchylky? ### hist((means-0)/sigma,freq=FALSE) curve(dnorm(x),add=TRUE) hist((means-0)/sds,freq=FALSE) curve(dnorm(x),add=TRUE) # odpovida normalni hustota lepe pro male nebo velke rozsahy vyberu? ### Spocitame konfidencni intervaly s odhadnutym rozptylem low.lim2 <- means - qt(1-alpha/2,n-1)*sds/sqrt(n) upp.lim2 <- means + qt(1-alpha/2,n-1)*sds/sqrt(n) windows() # nove graficke okno plot(means, pch=20, ylim=c(-.75,.75), xlab="Cislo vyberu", ylab="Bodovy a intervalovy odhad mu", main=paste(100*(1-alpha),"% oboustranne intervaly spolehlivosti pro mu pri neznamem sigma")) for (i in 1:opak){ lines(c(i,i),c(low.lim2[i],upp.lim2[i])) lines(c(i-.5,i+.5),c(low.lim2[i],low.lim2[i])) lines(c(i-.5,i+.5),c(upp.lim2[i],upp.lim2[i])) } abline(h=mu,lwd=2,col="red") # graphics.off() # ### Jak to funguje pro data s exponencialnim rozdelenim? for (i in 1:opak){ data[,i] <- rexp(n,mu,sigma) } means <- apply(data,2,mean) ### atd. # ##### Rozptyl normalniho rozdeleni # ### za predpokladu normality ma rozptyl chi-kvadrat rozdeleni - jak vypadaji konfidencni intervaly?? # vars <- apply(data,2,var) low.lim.var <- (n-1)*vars/qchisq(1-alpha/2,n-1) upp.lim.var <- (n-1)*vars/qchisq(alpha/2,n-1) plot(vars, pch=20, ylim=c(0,2.5), xlab="Cislo vyberu", ylab="Bodovy a intervalovy odhad pro sigma^2", main=paste(100*(1-alpha),"% oboustranne intervaly spolehlivosti pro sigma^2")) for (i in 1:opak){ lines(c(i,i),c(low.lim.var[i],upp.lim.var[i])) lines(c(i-.5,i+.5),c(low.lim.var[i],low.lim.var[i])) lines(c(i-.5,i+.5),c(upp.lim.var[i],upp.lim.var[i])) } abline(h=sigma^2,lwd=2,col="red") length(upp.lim.var[upp.lim.varsigma^2]) # pocet intervalu nad skutecnou hodnotou length(low.lim.var[low.lim.varsigma^2]) # pocet intervalu prekryvajicich sigma^2 opak*(1-alpha) # ocekavany pocet intervalu prekryvajicich sigma^2 # # ### UKOL: samostatne prozkoumat vlastnosti jednostrannych konf. intervalu (horniho a dolniho odhadu) ### Pokud zbude cas: ######################################## ######## Prace s daty (odhady) ##### ######## - histogram a boxplot ##### ######## ##### ######################################## ### Ukol: maji data o gejzirech normalni rozdeleni? ######### natahneme si data ######################## data() #### vypise defaultni moznosti ############### data(faithful) #### "natahne" konkretni data ################# help(faithful) #### informace o datech ############ attach(faithful) ## funkce pro zakladni popisne statistiky summary(eruptions) summary(eruptions)[1] summary(eruptions)["Min."] min(eruptions) var(eruptions) # vyberovy rozptyl S^2 sd(eruptions) # odmocnina z vyberoveho rozptylu S sd(eruptions)**2 # S^2 sd(eruptions)^2 # S^2 jinak m.eruptions<-c(mean(eruptions),sd(eruptions)) ### vyroba vlastniho vystupu, pouzijeme typ datove struktury list (seznam) moje.summary<-list(qant=summary(eruptions),rozp=c(var(eruptions),sd(eruptions),length(eruptions))) moje.summary$rozp ############################################################# #### hustota normalniho rozdeleni ########################## x.eruptions<-seq(min(eruptions),max(eruptions), by=0.1) f.eruptions<-dnorm(x.eruptions, mean=m.eruptions[1],sd=m.eruptions[2]) #### distrib. fce normalniho rozdel. ######################## F.eruptions<-pnorm(x.eruptions, mean=m.eruptions[1],sd=m.eruptions[2]) #### empiricka distribucni fce ############################## eF.eruptions<-ecdf(eruptions) ### obrazek ################################################# pic.eruptions<-par(mfrow=c(2,3)) ### vytvori multiobrazek 2 x 3 pole plot(eruptions) ### 1. obrazek plot(waiting) ### 2. obrazek boxplot(eruptions, horizontal=FALSE) ### 3. obrazek boxplot(eruptions, waiting, notch=TRUE, names = c("erup", "wait")) ### 4. obrazek hist(eruptions, prob=TRUE) # 5. obrazek histogram (na ose y relativni cetnosti) rug(eruptions) # na osu (eruptions) prida pozorovani density(eruptions) # vypise summary jadroveho odhadu hustoty lines(density(eruptions), lty=1) # prikresli odhad do tehoz grafu, plnou carou lines(x.eruptions, f.eruptions, lty=2) # prikresli hust, N rozd. do tehoz grafu, carkovane plot(eF.eruptions) # 6. obrazek: emp. distr. fce vs Normalni distr. fce lines(x.eruptions,F.eruptions, lty=2) pic.standard <- par(mfrow=c(1,1)) qqnorm(eruptions) # 7. obrazek normalni diagram qqline(eruptions, lty=2) ################################################################### ####### zavislost histogramu na sirce binu a levem pocatecnim bode par(mfrow=c(3,3)) ### vytvori multiobrazek 3 x 3 pole hist(eruptions, main="Default: Sturges") rug(eruptions) hist(eruptions, prob=TRUE, breaks = c(1:6), main=paste("Equidistant, 5 binu")) rug(eruptions) bin <- 0.8 for (d in seq(-bin,0, length=7)){ hist(eruptions, prob=TRUE, breaks = seq(min(eruptions)+d, max(eruptions)+bin, by=bin), main=paste("Zacatek: ", round((min(eruptions)+d),digits=2), " Posun: ",round((.8/6), digits=3)), xlim=c(0.7,6)) rug(eruptions) }