pneu <- read.csv("pneu.csv") load("sign_test.RData") save(pneu,sign.test,file="cvic_k10_3.RData") # save(sign.test,file="sign_test.RData") load(file="cvic_k10_3.RData") siteaddr <- "http://www.karlin.mff.cuni.cz/~pesta/NSTP097" datafile <- paste(siteaddr,"cvic_k10_3.RData",sep="/") load(url(datafile)) ### Sily a hladiny nobs <- 60 x <- rnorm(nobs,mean=0.5,sd=sqrt(2)) ## Jednovyb. KS test ks.test(x,y="pnorm",mean=0.5,sd=sqrt(2)) # One-sample Kolmogorov-Smirnov test # # data: x # D = 0.0801, p-value = 0.8066 # alternative hypothesis: two-sided od <- min(x)*0.9 do <- max(x)*1.1 plot(ecdf(x),xlim=c(od,do)) body <- seq(od,do,length=500) lines(body,pnorm(body,mean=0.5,sd=sqrt(2)),col="blue") ks.test(x,pnorm,mean=0.55,sd=sqrt(2)) od <- min(x)*0.9 do <- max(x)*1.1 plot(ecdf(x),xlim=c(od,do)) body <- seq(od,do,length=500) lines(body,pnorm(body,mean=0.55,sd=sqrt(2)),col="blue") ks.test(x,y="pnorm",mean=0.6,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.65,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.7,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.75,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.8,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.85,sd=sqrt(2)) ks.test(x,y="pnorm",mean=0.9,sd=sqrt(2)) # One-sample Kolmogorov-Smirnov test # # data: x # D = 0.1735, p-value = 0.04748 # alternative hypothesis: two-sided od <- min(x)*0.9 do <- max(x)*1.1 plot(ecdf(x),xlim=c(od,do)) body <- seq(od,do,length=500) lines(body,pnorm(body,mean=0.9,sd=sqrt(2.0)),col="blue") ## Jednovyberovy t-test t.test(x,mu=0.5) # One Sample t-test # # data: x # t = -0.2939, df = 59, p-value = 0.7699 # alternative hypothesis: true mean is not equal to 0.5 # 95 percent confidence interval: # 0.03747135 0.84405361 # sample estimates: # mean of x # 0.4407625 t.test(x,mu=0.55) t.test(x,mu=0.6) t.test(x,mu=0.65) t.test(x,mu=0.75) t.test(x,mu=0.8) t.test(x,mu=0.85) # # One Sample t-test # # data: x # t = -2.0305, df = 59, p-value = 0.04682 # alternative hypothesis: true mean is not equal to 0.85 # 95 percent confidence interval: # 0.03747135 0.84405361 # sample estimates: # mean of x # 0.4407625 ## Simulace sily a <- ks.test(x,y="pnorm",mean=0.55,sd=sqrt(2)) vem.ph <- function(x,test,...) { test(x,...)$p.value } vem.ph(x,ks.test,y="pnorm",mean=0.55,sd=sqrt(2)) vem.ph(x,t.test,mu=0.5) nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) ks.ph <- apply(x.vyb,2,vem.ph,ks.test,y="pnorm",mean=0.55,sd=sqrt(2)) hist(ks.ph) ks.test(ks.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: ks.ph # D = 0.0166, p-value = 0.9457 # alternative hypothesis: two-sided mean(ks.ph<0.05) # [1] 0.046 ks.ph <- apply(x.vyb,2,vem.ph,ks.test,y="pnorm",mean=0.85,sd=sqrt(2)) hist(ks.ph) ks.test(ks.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: ks.ph # D = 0.466, p-value < 2.2e-16 # alternative hypothesis: two-sided mean(ks.ph<0.05) # [1] 0.386 t.ph <- apply(x.vyb,2,vem.ph,t.test,mu=0.5) hist(t.ph) ks.test(t.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: t.ph # D = 0.0166, p-value = 0.945 # alternative hypothesis: two-sided mean(t.ph<0.05) # [1] 0.044 t.ph <- apply(x.vyb,2,vem.ph,t.test,mu=0.85) hist(t.ph) ks.test(t.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: ks.ph # D = 0.5442, p-value < 2.2e-16 # alternative hypothesis: two-sided mean(t.ph<0.05) # [1] 0.468 ## 2. Parove normalni nobs <- 60 x <- rnorm(nobs,mean=0.5,sd=sqrt(2)) y <- x+rnorm(nobs,mean=0,sd=1) mean(x) mean(y) cor(x,y) # [1] 0.8611978 ## skutecna korelace sqrt(2/3) # [1] 0.8164966 ## Parovy t-test t.test(x,y,paired=T) # Paired t-test # # data: x and y # t = 0.3692, df = 59, p-value = 0.7133 # alternative hypothesis: true difference in means is not equal to 0 # 95 percent confidence interval: # -0.9728426 1.4130681 # sample estimates: # mean of the differences # 0.2201127 y <- x+rnorm(nobs,mean=0,sd=1) t.test(x,y,paired=T) y <- x+rnorm(nobs,mean=0.2,sd=1) t.test(x,y,paired=T) # # One Sample t-test # # data: x # t = -2.0305, df = 59, p-value = 0.04682 # alternative hypothesis: true mean is not equal to 0.85 # 95 percent confidence interval: # 0.03747135 0.84405361 # sample estimates: # mean of x # 0.4407625 ## Simulace sily nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) y.vyb <- x.vyb+matrix(rnorm(nobs*nvyb,mean=0,sd=1),nrow=nobs,ncol=nvyb) oba.vyb <- rbind(x.vyb,y.vyb) vem.ph.2 <- function(data,test,...) { test(x=data[1:nobs],y=data[-c(1:nobs)],...)$p.value } t.ph <- apply(oba.vyb,2,vem.ph.2,t.test,paired=T) hist(t.ph) ks.test(t.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: t.ph # D = 0.0333, p-value = 0.2164 # alternative hypothesis: two-sided mean(t.ph<0.05) # [1] 0.065 nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) y.vyb <- x.vyb+matrix(rnorm(nobs*nvyb,mean=0.1,sd=1),nrow=nobs,ncol=nvyb) oba.vyb <- rbind(x.vyb,y.vyb) t.ph <- apply(oba.vyb,2,vem.ph.2,t.test,paired=T) hist(t.ph) ks.test(t.ph,y="punif") # One-sample Kolmogorov-Smirnov test # # data: t.ph # D = 0.1399, p-value < 2.2e-16 # alternative hypothesis: two-sided mean(t.ph<0.05) # 0.128 nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) y.vyb <- x.vyb+matrix(rnorm(nobs*nvyb,mean=0.2,sd=1),nrow=nobs,ncol=nvyb) oba.vyb <- rbind(x.vyb,y.vyb) t.ph <- apply(oba.vyb,2,vem.ph.2,t.test,paired=T) hist(t.ph) mean(t.ph<0.05) # 0.318 nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) y.vyb <- x.vyb+matrix(rnorm(nobs*nvyb,mean=0.3,sd=1),nrow=nobs,ncol=nvyb) oba.vyb <- rbind(x.vyb,y.vyb) t.ph <- apply(oba.vyb,2,vem.ph.2,t.test,paired=T) hist(t.ph) mean(t.ph<0.05) # 0.62 nvyb <- 1000 x.vyb <- matrix(rnorm(nobs*nvyb,0.5,sqrt(2)),nrow=nobs,ncol=nvyb) y.vyb <- x.vyb+matrix(rnorm(nobs*nvyb,mean=0.4,sd=1),nrow=nobs,ncol=nvyb) oba.vyb <- rbind(x.vyb,y.vyb) t.ph <- apply(oba.vyb,2,vem.ph.2,t.test,paired=T) hist(t.ph) mean(t.ph<0.05) # 0.88 ### Pneu pneu dim(pneu) pneu$met.v attach(pneu) met.v plot(pneu) cor(pneu) hist(met.v-met.d) t.test(met.v,met.d,paired=T) # Paired t-test # # data: met.v and met.d # t = 5.6503, df = 15, p-value = 4.614e-05 # alternative hypothesis: true difference in means is not equal to 0 # 95 percent confidence interval: # 2.837493 6.275007 # sample estimates: # mean of the differences # 4.55625 qt(0.975,df=15) # [1] 2.131450 sign.test(met.v,met.d) sign.test <- function (x, y = NULL, md = 0, alternative = "two.sided", conf.level = 0.95) { ## from library(BSDA) choices <- c("two.sided", "greater", "less") alt <- pmatch(alternative, choices) alternative <- choices[alt] if (length(alternative) > 1 || is.na(alternative)) stop("alternative must be one \"greater\", \"less\", \"two.sided\"") if (!missing(md)) if (length(md) != 1 || is.na(md)) stop("median must be a single number") if (!missing(conf.level)) if (length(conf.level) != 1 || is.na(conf.level) || conf.level < 0 || conf.level > 1) stop("conf.level must be a number between 0 and 1") if (is.null(y)) { dname <- paste(deparse(substitute(x))) x <- sort(x) diff <- (x - md) n <- length(x) nt <- length(x) - sum(diff == 0) s <- sum(diff > 0) estimate <- median(x) method <- c("One-sample Sign-Test") names(estimate) <- c("median of x") names(md) <- "median" names(s) <- "s" CIS <- "Conf Intervals" if (alternative == "less") { pval <- sum(dbinom(0:s, nt, 0.5)) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)][1] if (k < 1) { conf.level <- (1 - (sum(dbinom(k, n, 0.5)))) xl <- -Inf xu <- x[n] ici <- c(xl, xu) } else { ci1 <- c(-Inf, x[n - k + 1]) acl1 <- (1 - (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(-Inf, x[n - k]) acl2 <- (1 - (sum(dbinom(0:k, n, 0.5)))) xl <- -Inf xu <- (((x[n - k + 1] - x[n - k]) * (conf.level - acl2))/(acl1 - acl2)) + x[n - k] ici <- c(xl, xu) } } else if (alternative == "greater") { pval <- (1 - sum(dbinom(0:s - 1, nt, 0.5))) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)][1] if (k < 1) { conf.level <- (1 - (sum(dbinom(k, n, 0.5)))) xl <- x[1] xu <- Inf ici <- c(xl, xu) } else { ci1 <- c(x[k], Inf) acl1 <- (1 - (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(x[k + 1], Inf) acl2 <- (1 - (sum(dbinom(0:k, n, 0.5)))) xl <- (((x[k] - x[k + 1]) * (conf.level - acl2))/(acl1 - acl2)) + x[k + 1] xu <- Inf ici <- c(xl, xu) } } else { p1 <- sum(dbinom(0:s, nt, 0.5)) p2 <- (1 - sum(dbinom(0:s - 1, nt, 0.5))) pval <- min(2 * p1, 2 * p2, 1) return(rval, Confidence.Intervals) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)/2][1] if (k < 1) { conf.level <- (1 - 2 * (sum(dbinom(k, n, 0.5)))) xl <- x[1] xu <- x[n] ici <- c(xl, xu) } else { ci1 <- c(x[k], x[n - k + 1]) acl1 <- (1 - 2 * (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(x[k + 1], x[n - k]) acl2 <- (1 - 2 * (sum(dbinom(0:k, n, 0.5)))) xl <- (((x[k] - x[k + 1]) * (conf.level - acl2))/(acl1 - acl2)) + x[k + 1] xu <- (((x[n - k + 1] - x[n - k]) * (conf.level - acl2))/(acl1 - acl2)) + x[n - k] ici <- c(xl, xu) } } } else { if (length(x) != length(y)) stop("Length of x must equal length of y") xy <- sort(x - y) diff <- (xy - md) n <- length(xy) nt <- length(xy) - sum(diff == 0) s <- sum(diff > 0) dname <- paste(deparse(substitute(x)), " and ", deparse(substitute(y)), sep = "") estimate <- median(xy) method <- c("Dependent-samples Sign-Test") names(estimate) <- c("median of x-y") names(md) <- "median difference" names(s) <- "S" CIS <- "Conf Intervals" if (alternative == "less") { pval <- sum(dbinom(0:s, nt, 0.5)) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)][1] if (k < 1) { conf.level <- (1 - (sum(dbinom(k, n, 0.5)))) xl <- -Inf xu <- xy[n] ici <- c(xl, xu) } else { ci1 <- c(-Inf, xy[n - k + 1]) acl1 <- (1 - (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(-Inf, xy[n - k]) acl2 <- (1 - (sum(dbinom(0:k, n, 0.5)))) xl <- -Inf xu <- (((xy[n - k + 1] - xy[n - k]) * (conf.level - acl2))/(acl1 - acl2)) + xy[n - k] ici <- c(xl, xu) } } else if (alternative == "greater") { pval <- (1 - sum(dbinom(0:s - 1, nt, 0.5))) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)][1] if (k < 1) { conf.level <- (1 - (sum(dbinom(k, n, 0.5)))) xl <- xy[1] xu <- Inf ici <- c(xl, xu) } else { ci1 <- c(xy[k], Inf) acl1 <- (1 - (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(xy[k + 1], Inf) acl2 <- (1 - (sum(dbinom(0:k, n, 0.5)))) xl <- (((xy[k] - xy[k + 1]) * (conf.level - acl2))/(acl1 - acl2)) + xy[k + 1] xu <- Inf ici <- c(xl, xu) } } else { p1 <- sum(dbinom(0:s, nt, 0.5)) p2 <- (1 - sum(dbinom(0:s - 1, nt, 0.5))) pval <- min(2 * p1, 2 * p2, 1) loc <- c(0:n) prov <- (dbinom(loc, n, 0.5)) k <- loc[cumsum(prov) > (1 - conf.level)/2][1] if (k < 1) { conf.level <- (1 - 2 * (sum(dbinom(k, n, 0.5)))) xl <- xy[1] xu <- xy[n] ici <- c(xl, xu) } else { ci1 <- c(xy[k], xy[n - k + 1]) acl1 <- (1 - 2 * (sum(dbinom(0:k - 1, n, 0.5)))) ci2 <- c(xy[k + 1], xy[n - k]) acl2 <- (1 - 2 * (sum(dbinom(0:k, n, 0.5)))) xl <- (((xy[k] - xy[k + 1]) * (conf.level - acl2))/(acl1 - acl2)) + xy[k + 1] xu <- (((xy[n - k + 1] - xy[n - k]) * (conf.level - acl2))/(acl1 - acl2)) + xy[n - k] ici <- c(xl, xu) } } } if (k < 1) { cint <- ici attr(cint, "conf.level") <- conf.level rval <- structure(list(statistic = s, p.value = pval, estimate = estimate, null.value = md, alternative = alternative, method = method, data.name = dname, conf.int = cint)) oldClass(rval) <- "htest" return(rval) } else { result1 <- c(acl2, ci2) result2 <- c(conf.level, ici) result3 <- c(acl1, ci1) Confidence.Intervals <- round(as.matrix(rbind(result1, result2, result3)), 4) cnames <- c("Conf.Level", "L.E.pt", "U.E.pt") rnames <- c("Lower Achieved CI", "Interpolated CI", "Upper Achieved CI") dimnames(Confidence.Intervals) <- list(rnames, cnames) cint <- ici attr(cint, "conf.level") <- conf.level rval <- structure(list(statistic = s, p.value = pval, estimate = estimate, null.value = md, alternative = alternative, method = method, data.name = dname, conf.int = cint)) oldClass(rval) <- "htest" ## return(rval, Confidence.Intervals) return(rval) } }