################################################################## ################################################################## ######### Intervaly spolehlivosti: ################################################################## # 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 a vykreslime 95% intervaly spolehlivosti: 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). 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 # # Provedte znovu se stejnymi daty, ale pro alpha = 0.1. Jak cekate, ze se zmeni bodove odhady, intervaly spolehlivosti a frekvence chyb? # ##### Stredni hodnota normalniho rozdeleni s neznamym rozptylem # sds <- apply(data,2,sd) low.lim2 <- means - qt(1-alpha/2,n-1)*sds/sqrt(n) upp.lim2 <- means + qt(1-alpha/2,n-1)*sds/sqrt(n) X11() # 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") # ##### Rozptyl normalniho rozdeleni # 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 # # # ################################################################## ######### Testovani hypotez: ################################################################## # # Predpokladejme, ze za platnosti H_0 ma testova statistika U rozdeleni N(0,1), # at pro jednoduchost mu_0 = 0 a testujeme proti alternative mu_1 > 0 # grid <- seq(from=-3.5, to=7.5, by=.01) plot(grid, dnorm(grid), type="l", lwd=2, col="blue", main="Hustota statistiky U pri H_0 a 2 alternativach", ylab="Hustoty", xlab="U(x)") # hustota N(0,1) abline(v=qnorm(.95), lwd=3, col="red") #vertikalni cara s kritickou hodnotou legend(x=-3.5, y=.4, legend=c("krit. hodnota","mu = 0 (Ho)","mu = 2.5","mu = 4"), lwd=c(3,2,1,1), col=c("red","blue","brown","green"), bty="n") # # Nez budete pokracovat, uvedomte si, co je chyba I. druhu a co je zde kriticky obor testu !!! # # Alternativy: lines(grid, dnorm(grid,mean=2.5,sd=1), type="l", lwd=1, col="brown") # # Nez budete pokracovat, uvedomte si, co je zde sila testu a co je chyba II. druhu !!! # lines(grid, dnorm(grid,mean=4,sd=1), type="l", lwd=1, col="green") # # Porovnejte chyby II. druhu (sily testu) pro obe alternativy - dal se takovy vysledek ocekavat? # pnorm(qnorm(.95)-2.5) # pravdepodobnost chyby 2. druhu pro mu=2.5 pnorm(qnorm(.95),mean=2.5,sd=1) # ekvivalentni zpusob vypoctu pnorm(qnorm(.95)-4) # pravdepodobnost chyby 2. druhu pro mu=4 pnorm(qnorm(.95),mean=4,sd=1) # ekvivalentni zpusob vypoctu # # # ---------------------------------------------------------------- # # Jednostranne testy s ruznou hladinou spolehlivosti, jmenovite alfa = 1%, 5% a 10%: # # Jak se "seradi" jednotlive kriticke hodnoty? # Ktery test (tj. na jake hladine vyznamnosti) bude zamitat H_0 nejcasteji? Proc? # grid <- seq(from=-3.5, to=3.5, by=.01) plot(grid, dnorm(grid,mean=0,sd=1), type="l", lwd=2, col="blue", main=c("Kriticke hodnoty c pro jednostr. test, hladiny alpha = 10, 5 a 1 %","pro normalni rozdeleni se znamym rozptylem"), xlab="U(x)", ylab="Hustota N(0,1)") legend(x=-3.5, y=.4, legend=c(paste("c = ", round(qnorm(.90),3),", hladina 10%"), paste("c = ", round(qnorm(.95),3),", hladina 5%"), paste("c = ", round(qnorm(.99),3),", hladina 1%"), "Hustota N(0,1)"), lty=c(2,1,3,1), lwd=2, col=c("red","red","red","blue"), bty="n") abline(v=qnorm(.90), lwd=2, lty=2, col="red") abline(v=qnorm(.95), lwd=2, lty=1, col="red") abline(v=qnorm(.99), lwd=2, lty=3, col="red") # # ---------------------------------------------------------------- # Jednostranny versus oboustranny test: # # Chyba I. druhu je u obou testu stejna, tj. 5%. Porovnejte chybu II. druhu (silu testu) pro alternativu mu = 3 !! Jak by to dopadlo napr. pri alternative mu = -2? # grid <- seq(from=-3.5, to=3.5, by=.01) plot(grid, dnorm(grid), type="l", lwd=2, col="blue", main="Kriticke hodnoty c na 5% pro jednostr. a oboustr. test", ylab="Hustota N(0,1)", xlab="U(x) s vyznacenymi krit. hodnotami") abline(v=qnorm(.95), lwd=2, col="red", lty=1) abline(v=qnorm(1-.05/2), lwd=2, col="red", lty=3) lines(c(-qnorm(1-.05/2),-qnorm(1-.05/2)), c(-.1,.25), lwd=2, col="red", lty=3) lines(grid, dnorm(grid,3,1), type="l", lwd=1, col="green") legend(x=-3.5, y=.36, legend=c(paste("c = ", round(qnorm(.95),3)," pro jednostranny 5% test"), paste("c = -+", round(qnorm(1-.05/2),3)," pro oboustranny 5% test"), "N(3,1) pro alternativu mu=3", "Hustota U(x) za H_0"), lty=c(1,3,1,1), lwd=c(2,2,1,2), col=c("red","red","green","blue"), bty="n") # # ----------------------------------------------------------------- # # Oboustranny test H0: mu=mu_0=0 proti H1:mu<>0): # Nagenerujeme 300 pozorovani z N(0.2,1) # data1 <- rnorm(300, mean=0.2, sd=1) # #library(ctest) # natahneme si prislusnou knihovnu (ctest jako Classical Tests), # v novejsich verzich uz je soucasti knihovny stats a neni treba ji natahovat #library(help=ctest) # muzeme se podivat, jake funkce (testy) mame k dispozici # # Tady mame generovana data, takze zname rozptyl -> meli bychom pouzit test zalozeny na normalnim rozdeleni. Ten zde nenajdeme (prakticky se nepouziva), pouzijeme t-test. Jak na to? help(t.test) # help(t.test,package="ctest") # takhle se muselo volat drive t.test(x=data1, alternative="two.sided", mu=0, conf.level=.95) # !!! Zalezi na konkretnim vyberu, soused muze mit jiny vysledek !!! # p-hodnota je nejmensi hladina, na ktere bychom jeste hypotezu zamitli # # Zkuste generovat data s jinou stredni hodnotou a zkousejte citlivost testu na rozdil od mu=0 # A 300 pozorovani je asi trochu moc, dejte rekneme jenom 20 nebo 25 pozorovani # # Zkusime dva vybery data1 <- rnorm(30, 0.5, 1) data2 <- rnorm(20, 0.6, 2) # muzeme pouzit dvouvyberovy t-test? t.test(x=data1,y=data2,alternative="two.sided",var.equal=TRUE,conf.level=.95) # testujme rovnost rozptylu pomoci F-testu help(var.test) var.test(data1,data2,alternative="two.sided") # q() # Koncime ...