variacni_rada<-function(X, nazvy){
  X <- as.numeric(X)
  
  nj <- rep(0,max(X))
  for(j in min(X):max(X)){
    nj[j]<-sum(X==j)
  }
  
  n  <- sum(nj)
  pj <- nj/n
  Nj <- cumsum(nj)
  Fj <- cumsum(pj)
  
  (variacni.rada <- data.frame(nj=nj, Nj=Nj, pj=pj, Fj=Fj))
  row.names(variacni.rada) <- nazvy
  variacni.rada
}


dotplot<-function(X, Y, main='Dotplot', xlab='X', ylab='Y', xlim=c(min(X),max(X)), 
                  ylim=c(min(Y),max(Y)), col='black', pch=21, bg='white', cex=1, lwd=1){
  rand <- rnorm(length(X),0,0.03)
  X2 <- X+rand
  plot(X2,Y,type='p',main=main,xlab=xlab,ylab=ylab,
       xlim=xlim,ylim=ylim,col=col,pch=pch,bg=bg,lwd=lwd)
}

norm2       <- function(x, y, mu1, mu2, sigma1, sigma2){
  rho       <- 0
  Sigma     <- matrix(c(sigma1^2, sigma2*sigma1*rho, sigma1*sigma2*rho, sigma2^2), 2, 2, byrow=T)
  xy        <- c(x[1] - mu1, y[1] - mu2)
  konstanta <- 1/(2*pi*sqrt(sigma1*sigma2*(1-rho^2)))
  hustota   <- konstanta*exp(-1/2*t(xy)%*%solve(Sigma)%*%xy)
  return(hustota)
}

#Scheffé
Scheffe <- function(X, group, names, alpha){
ID <- group
r <- length(unique(ID))
n <- length(X)

Xi. <- Mi. <- ni <- NULL
for(i in 1:r){
  Xi.[i] <- sum(X[ID==i])
  ni[i]  <- length(X[ID==i])
  Mi.[i] <- sum(X[ID==i])/ni[i]
}
X.. <- sum(X)
M.. <- mean(X)

SA <- sum(ni*(Mi.-M..)^2)
fA <- r-1
ST <- sum((X-M..)^2)
SE <- ST-SA
fE <- n-r
Fa <- (SA/fA)/(SE/fE)

Scheffe.R <- matrix(NA, r, r)
Scheffe.L <- matrix(NA, r, r)

Sh <- sqrt(SE/fE)
for(k in 1:r){
  for(j in 1:r){
    Scheffe.R[k,j] <- Sh*sqrt((r-1)*(1/ni[k]+1/ni[j])*qf(1-alpha,r-1,n-r))
    Scheffe.L[k,j] <- abs(Mi.[k]-Mi.[j])
  }
}
Scheffe.R <- data.frame(Scheffe.R, row.names=names)
names(Scheffe.R) <- names
Scheffe.L <- data.frame(Scheffe.L, row.names=names)
names(Scheffe.L) <- names
return(list(R=Scheffe.R, L=Scheffe.L))
}