gnp <- function(n)
  {
    fn <- ecdf(rnorm(n))
    plot(fn,xlim=c(-2.5,2.5),cex=0.3,
         main=paste("n =",deparse(substitute(n))))
    x <- seq(-2.5,2.5,length=500)
    lines(x,pnorm(x),col="blue")
    eps <- min(diff(knots(fn)))/2
    maxdif.u <- max(abs(fn(knots(fn)-eps)-pnorm(knots(fn))))
    maxdif.l <- max(abs(fn(knots(fn)+eps)-pnorm(knots(fn))))
    return(max(maxdif.u,maxdif.l))
  }

plotCI <- function (x, y = NULL, uiw, liw = uiw, ui, li, err = "y", ylim = NULL, 
    xlim = NULL, type = "p", col = par("col"), barcol = col, 
    pt.bg = par("bg"), sfrac = 0.01, gap = 1, lwd = par("lwd"), 
    lty = par("lty"), labels = FALSE, add = FALSE, xlab, ylab, 
    minbar, maxbar, ...) 
{
  ## from gplots library
  invalid <- function (x) 
    {
      if (missing(x) || is.null(x) || length(x) == 0) 
        return(TRUE)
      if (is.list(x)) 
        return(all(sapply(x, invalid)))
      else if (is.vector(x)) 
        return(all(is.na(x)))
      else return(FALSE)
    }
  
  
    if (is.list(x)) {
        y <- x$y
        x <- x$x
    }
    if (invalid(xlab)) 
        xlab <- deparse(substitute(x))
    if (invalid(ylab)) {
        if (is.null(y)) {
            xlab <- ""
            ylab <- deparse(substitute(x))
        }
        else ylab <- deparse(substitute(y))
    }
    if (is.null(y)) {
        if (is.null(x)) 
            stop("both x and y NULL")
        y <- as.numeric(x)
        x <- seq(along = x)
    }
    if (err == "y") 
        z <- y
    else z <- x
    if (invalid(ui)) 
        ui <- z + uiw
    if (invalid(li)) 
        li <- z - liw
    if (!invalid(minbar)) 
        li <- ifelse(li < minbar, minbar, li)
    if (!invalid(maxbar)) 
        ui <- ifelse(ui > maxbar, maxbar, ui)
    if (err == "y") {
        if (is.null(ylim)) 
            ylim <- range(c(y, ui, li), na.rm = TRUE)
        if (is.null(xlim) && !is.R()) 
            xlim <- range(x, na.rm = TRUE)
    }
    else if (err == "x") {
        if (is.null(xlim)) 
            xlim <- range(c(x, ui, li), na.rm = TRUE)
        if (is.null(ylim) && !is.R()) 
            ylim <- range(x, na.rm = TRUE)
    }
    if (!add) {
        if (invalid(labels) || labels == FALSE) 
            plot(x, y, ylim = ylim, xlim = xlim, col = col, xlab = xlab, 
                ylab = ylab, ...)
        else {
            plot(x, y, ylim = ylim, xlim = xlim, col = col, type = "n", 
                xlab = xlab, ylab = ylab, ...)
            text(x, y, label = labels, col = col, ...)
        }
    }
    if (is.R()) 
        myarrows <- function(...) arrows(...)
    else myarrows <- function(x1, y1, x2, y2, angle, code, length, 
        ...) {
        segments(x1, y1, x2, y2, open = TRUE, ...)
        if (code == 1) 
            segments(x1 - length/2, y1, x1 + length/2, y1, ...)
        else segments(x2 - length/2, y2, x2 + length/2, y2, ...)
    }
    if (err == "y") {
        if (gap != FALSE) 
            gap <- strheight("O") * gap
        smidge <- par("fin")[1] * sfrac
        if (!is.null(li)) 
            myarrows(x, li, x, pmax(y - gap, li), col = barcol, 
                lwd = lwd, lty = lty, angle = 90, length = smidge, 
                code = 1)
        if (!is.null(ui)) 
            myarrows(x, ui, x, pmin(y + gap, ui), col = barcol, 
                lwd = lwd, lty = lty, angle = 90, length = smidge, 
                code = 1)
    }
    else {
        if (gap != FALSE) 
            gap <- strwidth("O") * gap
        smidge <- par("fin")[2] * sfrac
        if (!is.null(li)) 
            myarrows(li, y, pmax(x - gap, li), y, col = barcol, 
                lwd = lwd, lty = lty, angle = 90, length = smidge, 
                code = 1)
        if (!is.null(ui)) 
            myarrows(ui, y, pmin(x + gap, ui), y, col = barcol, 
                lwd = lwd, lty = lty, angle = 90, length = smidge, 
                code = 1)
    }
    points(x, y, col = col, lwd = lwd, bg = pt.bg, type = type, 
        ...)
    invisible(list(x = x, y = y))
}

ci.asym <- function(x,alpha=0.05)
  {
    prum <- mean(x)
    sig2 <- var(x)
    n <- length(x)
    skok <- qt(1-alpha/2,length(x)-1)*sqrt(sig2/n)
    c(prum-skok,prum,prum+skok)
  }

pneu <- read.csv(file="pneu.csv")

save(gnp,plotCI,ci.asym,pneu,file="cvic_k10_2.RData")


plotCI(c(0.5,1,0.75,3), uiw=0.4, err = "y")


siteaddr <- "http://www.karlin.mff.cuni.cz/~pesta/NSTP097"
datafile <- paste(siteaddr,"cvic_k10_2.RData",sep="/")
load(url(datafile))


x=rnorm(25)
edf=ecdf(x)
edf(-0.25)
plot(edf)

x=rbeta(35,0.5,0.5)
edf=ecdf(x)
edf(0.25)
plot(edf)

par(mfrow=c(2,2))
gnp(10)
gnp(50)
gnp(500)
gnp(2000)
par(mfrow=c(1,1))



## intervaly spolehlivosti pro stredni hodnotu

smp <- rnorm(20,2,1)
prum <- mean(smp)

ci.asym(smp)

nobs <- 20
nvyb <- 100
data.mat <- matrix(rnorm(nobs*nvyb,2,1),nrow=nobs,ncol=nvyb)

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<2 & 2<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=2,col="red")



nobs <- 20
nvyb <- 1000
data.mat <- matrix(rnorm(nobs*nvyb,2,1),nrow=nobs,ncol=nvyb)

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<2 & 2<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=2,col="red")



nobs <- 100
nvyb <- 1000
data.mat <- matrix(rnorm(nobs*nvyb,2,1),nrow=nobs,ncol=nvyb)

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<2 & 2<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=2,col="red")



nobs <- 1000
nvyb <- 1000
data.mat <- matrix(rnorm(nobs*nvyb,2,1),nrow=nobs,ncol=nvyb)

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<2 & 2<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))



co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=2,col="red")




nobs <- 20
nvyb <- 1000
data.mat <- matrix(rexp(nobs*nvyb,5),nrow=nobs,ncol=nvyb)
skut.sh <- 0.2

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]< skut.sh & skut.sh<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=skut.sh,col="red")




nobs <- 100
nvyb <- 1000
data.mat <- matrix(rexp(nobs*nvyb,5),nrow=nobs,ncol=nvyb)
skut.sh <- 0.2

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<skut.sh & skut.sh<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=skut.sh,col="red")



nobs <- 1000
nvyb <- 1000
data.mat <- matrix(rexp(nobs*nvyb,5),nrow=nobs,ncol=nvyb)
skut.sh <- 0.2

vs.ci <- apply(data.mat,2,ci.asym)

(pokryti <- sum(vs.ci[1,]<skut.sh & skut.sh<vs.ci[3,])/nvyb)
(prum.sirka <- mean(vs.ci[3,]-vs.ci[1,]))

co <- 1:100
plotCI(vs.ci[2,co],uiw=(vs.ci[3,co]-vs.ci[1,co])/2,gap=0.15,sfrac=0.002,
       ylab="Int. spol. pro str. hodnotu",xlab="Vyber")
abline(h=skut.sh,col="red")






pneu

dim(pneu)

pneu$met.v

attach(pneu)

met.v


mean(pneu)
var(pneu)
n <- nrow(pneu)
s3 <- 1/(n-1)*sum((met.v-mean(met.v))^2)
s3
cor(pneu)

a <- ecdf(met.v)
b <- ecdf(met.d)

plot(a)
lines(b,lty=2,col="blue")


ci.asym(met.v)
ci.asym(met.d)