### MATEMATICKA STATISTIKA 1 ### R cviceni c. 2 (cviceni c. 9) ### Kolmogorovovuv-Smirnovuv test, t-test ### ### -------------------------------------------------- ##### Nastaveni pracovniho adresare ##### ================================================= setwd("H:/nmsa331") rm(list=ls()) ## vycisteni ### alpha pro testy a intervaly spolehlivosti alpha <- 0.05 Hosi=read.table("Hosi.txt",header=TRUE) head(Hosi) summary(Hosi) # Opet se budeme zabyvat porodni hmotnosti # Udelame si tentokrat podvyber o rozsahu 100 set.seed(27112018); # nastavuje generator (pseudo)-nahodnych cisel hmot100 <- sample(Hosi$por.hmot, 100) n=length(hmot100) mean(hmot100) mean(Hosi$por.hmot) # zkoumani normality rozdeleni: qqnorm(hmot100) qqline(hmot100) hist(hmot100,prob=T) curve(dnorm(x,mean=mean(hmot100),sd=sd(hmot100)),min(hmot100),max(hmot100),col="red",add=T,lwd=2) #provedeni t-testu: t.test(hmot100,mu=3300) # testova statistika (tstat=sqrt(n)*(mean(hmot100)-3300)/sd(hmot100)) # kriticka hodnota qt(1 - alpha/2, df=n-1) # p-hodnota 2*(1-pt(abs(tstat), df=n-1)) # interval spolehlivosti: # dolni mez a dolni mez mean(hmot100)-qt(1 - alpha/2, df=n-1)/sqrt(n)*sd(hmot100) mean(hmot100)+qt(1 - alpha/2, df=n-1)/sqrt(n)*sd(hmot100) # pristup k jednotlivym polozkam t-testu: tt=t.test(hmot100,mu=3300) names(tt) tt$stat tt$p.val tt$conf.int tt$alter # jina spolehlivost t.test(hmot100, mu = 3300, conf.level = 0.99) # 1-stranny test t.test(hmot100,mu=3500,alternative="less") #p-hodnota tstat=(mean(hmot100)-3500)*sqrt(n)/sd(hmot100) pt(tstat,df=n-1) # interval spolehlivosti zalozeny na CLV: # dolni mez a dolni mez mean(hmot100)-qnorm(1 - alpha/2)/sqrt(n)*sd(hmot100) mean(hmot100)+qnorm(1 - alpha/2)/sqrt(n)*sd(hmot100) #### nyni porovname vysledky pro cela data a nas podvyber t.test(hmot100,mu=3300) t.test(Hosi$por.hmot,mu=3300) t.test(hmot100,mu=3500,alternative="less") t.test(Hosi$por.hmot,mu=3500,alternative="less") ############################# Iq=read.table("Iq.txt",header=T) head(Iq) summary(Iq) dim(Iq) # Casto se uvadi, napr. # # https://cs.wikipedia.org/wiki/Inteligen%C4%8Dn%C3%AD_kvocient # ze IQ v populaci ma priblizne normalni rozdeleni se stredni hodnotou 100 # a smerodatnou odchylkou 15. Otestujte, zda jsou nase data o zacich v souladu # s timto tvrzenim. # # Pro jednoduchost si ulozime data do promenne IQ IQ <- Iq$iq; (n <- length(IQ)); # rozsah vyberu ks.test(IQ, y="pnorm", mean=100, sd=15) #### # manualni vypocet # abychom nemuseli porad opakovat parametru, tak si F_0 ulozime do F_0 F0 <- function(x) pnorm(x, mean=100, sd=15) IQ.sorted=sort(IQ) Knplus <- max((1:n)/n - F0(IQ.sorted)) Knminus <- max(F0(IQ.sorted) - (0:(n-1))/n) (Kn <- max(Knplus, Knminus)) # distribucni funkce limitniho rozdeleni G <- function(y){ k <- 1:10000; # aproximace nekoncne sumy 1 - 2*sum((-1)^(k+1)*exp(-2*k^2*y^2)) } # p-hodnota 1 - G(sqrt(n)*Kn) # viz p-hodnota ve vystupu fce ks.test() # kdybychom chteli kritickou hodnotu: # (1-alfa)-kvantil rozdeleni G (tj. kriticka hodnota testu) se da spocist nasledovne: alfa <- 0.05; fce <- function(x) G(x) - (1 - alfa); # Nasledujici funkce najde koren rovnice (krit <- uniroot(fce, c(0.1,10))$root) # coz porovnavame s sqrt(n)*Kn ####### # # # Graficka ilustrace # Empiricka distribucni funkce empir.df <- ecdf(IQ); xgrid <- seq(min(IQ)-10, max(IQ)+10, length=10000); # sit bodu # Dolni pas spolehlivosti pro skutecnou d.f. dolni.pas <- pmax(empir.df(xgrid)-krit/sqrt(n),0) # Horni pas spolehlivosti pro skutecnou d.f. horni.pas <- pmin(empir.df(xgrid)+krit/sqrt(n),1) plot(empir.df, ylab="", main="Porovnani distribucnich funkci",do.points=F) # Vyplneni sedou barvou polygon(c(xgrid, rev(xgrid)), c(horni.pas, rev(dolni.pas)), col = "gray90", border="gray70",lty=2); # Empiricka distribucni funkce jeste jednou (protoze predchozi prikaz ji prekryl sedym pasem) plot(empir.df, verticals=FALSE, add=T, lwd=2,do.points=FALSE) # Predpokladana distribucni funkce za nulove hypotezy lines(xgrid, F0(xgrid), col="blue", lwd=2); # Bod nejvetsiho rozdilu k0 <- which.max(F0(IQ.sorted) - (0:(n-1))/n); abline(v=IQ.sorted[k0], col="red"); legend("bottomright", col=c("black", "grey70", "blue"), lty=c(1,2,1), legend=c("empiricka d.f.", "pas spolehlivosti", "d.f. za nulove hypotezy")) # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # Pozor. Funkce ks.test() nedava spravnou p-hodnotu, pokud parametry rozdeleni # nezname, ale odhadovali bychom je z dat. # V nasledujici simulacni studii budeme generovat nahodny vyber o rozsahu n=111 z normalniho # rozdeleni N(100, 15^2). Budeme porovnavat skutecnou hladinu testu # v pripadech, ze tyto parametry jsou dany a v pripadech, ze bychom si je # odhadovali z dat. nopak <- 10000; # pocet nahodnych vyberu pval.zname <- numeric(nopak); # sem si budu ukladat p-hodnoty v pripade znamych parametru pval.nezname <- numeric(nopak); # sem si budu ukladat p-hodnoty v pripade neznamych parametru n<-100 for(i in 1:nopak){ Y <- rnorm(n, mean=100, sd=15); # vygenerovany vyber z N(100,15^2) pval.zname[i] <- ks.test(Y, y="pnorm", mean=100, sd=15)$p.value; pval.nezname[i] <- ks.test(Y, y="pnorm", mean=mean(Y), sd=sd(Y))$p.value; } # Odhad skutecne hladiny testu pri znamych parametrech mean(pval.zname <= 0.05) # Odhad skutecne hladiny testu pri neznamych parametrech mean(pval.nezname <= 0.05) ##### # K-S test pro diskretni rozdeleni nopak=10000 pval=numeric(nopak) n=100 # rozsah vyberu N=5 # parametry rozdeleni Bi(N,p) p=1/3 F0=function(x) pbinom(x, size=N, prob=p) set.seed(123) for(i in 1:nopak){ x = rbinom(n,size=N,prob=p) xi=unique(c(x,0:N)) Kn=max(abs(ecdf(x)(xi)-F0(xi))) pval[i]=1-G(sqrt(n)*Kn) } mean(pval<=0.05) ##### # # T-Test: Simulace za (nulove) hypotezy # A. Rovnomerne rozdeleni a exp. rozdeleni za nulove hypotezy n=100 opak=1000 p.rovn=numeric(opak) p.exp=numeric(opak) for(i in 1:opak){ x1=runif(n) p.rovn[i]=t.test(x1,mu=1/2)$p.val x2=rexp(n,rate=1) p.exp[i]=t.test(x2,mu=1)$p.val } mean(p.rovn<=0.05) mean(p.exp<=0.05) # B. Rovnomerne rozdeleni a exp. rozdeleni za alternativy hypotezy n=100 opak=1000 p.rovn=numeric(opak) p.exp=numeric(opak) for(i in 1:opak){ x1=runif(n,0,1.2) p.rovn[i]=t.test(x1,mu=1/2)$p.val x2=rexp(n,rate=1.2) p.exp[i]=t.test(x2,mu=1)$p.val } mean(p.rovn<=0.05) mean(p.exp<=0.05) ##