### Graf funkce G z limitniho rozdeleni K-S testu G <- function(y){ k <- 1:10000; # aproximace nekoncne sumy 1 - 2*sum((-1)^(k+1)*exp(-2*k^2*y^2)) } curve(Vectorize(G)(x),0.01,2,main="Distribucni funkce G",ylab=expression(G(x))) # vsimnete si rozdil mezi nasledujicimi vyrazy G(1:5) Vectorize(G)(1:5) #################### ### Hladina t-testu pro ruzna rozdeleni a ruzna n ### Uvazujeme normalni rozdeleni, rovnomerne a exponencialni. ### Pro Rovn. a Exp. se jedna o asymptoticky test (hladina se blizi k alpha pro n->infty) ### Pro normalni rozdeleni by mela byt hladina = alpha pro libovolne n N = c(5,10,20,50,100,200,300,500) nN=length(N) Aunif = rep(NA,nN) Aexp = rep(NA,nN) Anorm = rep(NA,nN) opak=1000 for(j in 1:nN){ n = N[j] p.rovn=numeric(opak) p.exp=numeric(opak) p.norm=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 x3=rnorm(n,mean=0,sd=1) p.norm[i]=t.test(x3,mu=0)$p.val } Aunif[j] = mean(p.rovn<=0.05) Aexp[j] = mean(p.exp<=0.05) Anorm[j] = mean(p.norm<=0.05) } plot(c(rep(N,2),mean(N)),c(Aunif,Aexp,0.05),type="n",ylab=expression(alpha),xlab="Rozsah n", main="Hladina testu",ylim=c(0,0.2)) abline(h=0.05,lty=3,lwd=2) lines(N,Aunif,lty=2,col="blue",lwd=2) points(N,Aunif,pch=16,col="blue") lines(N,Aexp,lty=2,col=2,lwd=2) points(N,Aexp,col=2,pch=16) lines(N,Anorm,lty=2,col=3,lwd=2) points(N,Anorm,col=3,pch=16) legend("topright",c("Rovn","Exp","Normal"),pch=16,lty=2,col=c(4,2,3)) #### Tentyz obrazek pro silu testu #### V pripade konzistentniho testu: sila -> 1 pro n->infty za H_1 ## Uvazujeme Rovn[0,1.2] a H_0: EX=1/2 ## a Exp(lam=1.2) a H_0: EX=1 N = c(5,10,20,50,100,200,300,500) nN=length(N) Aunif = rep(NA,nN) Aexp = rep(NA,nN) opak=1000 for(j in 1:nN){ n = N[j] 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 } Aunif[j] = mean(p.rovn<=0.05) Aexp[j] = mean(p.exp<=0.05) } plot(c(rep(N,2),mean(N)),c(Aunif,Aexp,0.05),type="n",ylab="Sila testu",xlab="Rozsah n", main="Sila testu") abline(h=0.05,lty=3,lwd=2) lines(N,Aunif,lty=2) points(N,Aunif,pch=16) lines(N,Aexp,lty=2,col=2) points(N,Aexp,col=2,pch=16) legend("bottomright",c("Rovn","Exp"),pch=16,lty=2,col=1:2,cex=1.2) ########### ### Pas spolehlivosti pro df F vs. intervaly spolehlivosti pro F(x) Iq=read.table("Iq.txt",header=T) IQ <- Iq$iq; (n <- length(IQ)); # rozsah vyberu fce <- function(x) G(x) - (1 - alpha); # Nasledujici funkce najde koren rovnice (krit <- uniroot(fce, c(0.1,10))$root) empir.df <- ecdf(IQ); F0 <- function(x) pnorm(x, mean=100, sd=15) xgrid <- seq(min(IQ)-10, max(IQ)+10, length=10000); # sit bodu dolni.pas <- pmax(empir.df(xgrid)-krit/sqrt(n),0) horni.pas <- pmin(empir.df(xgrid)+krit/sqrt(n),1) dolni.int=pmax(0,empir.df(xgrid) - pnorm(.975)*empir.df(xgrid)*(1-empir.df(xgrid))/sqrt(n)) horni.int=pmin(1,empir.df(xgrid) + pnorm(.975)*empir.df(xgrid)*(1-empir.df(xgrid))/sqrt(n)) plot(empir.df, ylab="", main="Porovnani distribucnich funkci",do.points=F,lwd=2) polygon(c(xgrid, rev(xgrid)), c(horni.pas, rev(dolni.pas)), col = "gray90", border="gray70",lty=2); polygon(c(xgrid, rev(xgrid)), c(horni.int, rev(dolni.int)), col = "lightblue", border="blue",lty=2); plot(empir.df, verticals=FALSE, add=T, lwd=2,do.points=FALSE) lines(xgrid, F0(xgrid), col="red", 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", "interval spolehlivosti"))