data <- read.delim('16-anova-head.txt')
udaje <- data[data$sex == 'f' , c('head.L', 'bizyg.W')]
udaje <- na.omit(udaje)
n <- dim(udaje)[1]
(tab <- data.frame(n, min_X = min(udaje$head.L), max_X = max(udaje$head.L), 
                   min_Y = min(udaje$bizyg.W), max_Y = max(udaje$bizyg.W)))

MVN::mvn(udaje, mvnTest = 'mardia')$multivariateNormality
Z <- MASS::kde2d(udaje$head.L, udaje$bizyg.W, n = 50, 
                 lims = c(165, 205, 115, 155))
par(mar = c(0, 4, 0, 4), family = 'Times')
GA::persp3D(Z$x, Z$y, Z$z, xlab = 'delka hlavy (v mm)', ticktype = 'simple',
            ylab = 'sirka tvare (v mm)', zlab = "relativni cetnost", 
            col.palette = heat.colors, border = 'darkred',
            theta = 5, phi = 30)

par(mar = c(4, 4, 1, 1), family = 'Times')
elipsa <- car::dataEllipse(udaje$head.L, udaje$bizyg.W, level = 0.90, 
                           xlim = c(165, 205), ylim = c(115, 155), xlab = '', 
                           ylab = 'sirka tvare (v mm)', main = '', pch = 21,  
                           col = 'red3', bg = 'yellow', lwd = 2, las = 1)
mtext('delka hlavy (v mm)', side = 1, line = 2.1)

library(sp) # funkce point.in.polygon()
# lezi bod uvnitr (TRUE; 1) nebo vne (FALSE; 0) elipsy?
rozhodnuti <- as.logical(point.in.polygon(udaje$head.L, udaje$bizyg.W, elipsa[,1], elipsa[,2])) 
table(rozhodnuti) # 10 vne, 90 uvnitr
10 / n # 10% hodnot lezi mimo elipsu spolehlivosti