variacni.rada <- function(X, row.names){
  X  <- as.numeric(X)
  U  <- sort(unique(X))
  nu <- length(U)
  nj <- rep(0, nu)
  for(j in 1 : nu){
    nj[j] <- sum(X == U[j])
  }
  
  n  <- sum(nj)
  pj <- nj/n
  Nj <- cumsum(nj)
  Fj <- cumsum(pj)
  
  variacni.rada <- data.frame(nj = nj, pj = pj, Nj = Nj, Fj = Fj)
  row.names(variacni.rada) <- row.names
  variacni.rada
}

rel.barplot <- function(N, col = 1:length(N), border = 'black', names = 1:length(N), main = '', xlab = '', 
                        ylab = 'relativni cetnost', xlim = c(0.2, 2), ylim = c(-0.03, 1.03), 
                        density = 60, cex = 1, mtext = '', a = 0, axes = axes){
  n       <- sum(N)
  l       <- length(N)
  
  barplot(matrix(N / n, l, 1), col = col, border = border, density = density,
          main = main, xlim = xlim, ylim = ylim, ylab = ylab, axes = T, las = 1)
  legend('topright',legend = rev(names), fill = rev(col), density = density, bty = 'n')
  mtext(xlab, 1, line = 1)
  
  stred <- 0.7
  cn    <- cumsum(N) / n
  cn2   <- (N / n) / 2
  vyska <- c(cn2[1], cn[1 : (length(cn) - 1)] + cn2[2 : length(cn2)])
  vyska <- vyska + a
  text(stred, vyska, paste(N, '; ',round(N / n * 100, 2),'%',sep = ''), cex = cex)
}


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))
}