# Usiamo a titolo di esempio i dati in Esempio 
# Caricamento dati ----
dati <- read.table("Esempio.txt", header = TRUE, sep = ",")
head(dati)
str(dati)

# Salviamo nell'oggetto x l'altezza 
x <- dati$HEIGHT
n <- length(x)


# Indici di tendenza centrale ----

## Media aritmetica ----

M <- mean(x)
M

# Verifica proprietà
M >= min(x) & M <= max(x) # a: min <= M <= max
round(sum(x - M), 10)     # b: scarti nulli

# c: M minimizza la somma degli scarti quadratici
a       <- sort(c(seq(min(x), max(x), 0.5), mean(x)))
varTest <- sapply(a, function(a) mean((x - a)^2))
min.a   <- which.min(varTest) # individua la posizione dell'osservazione corrispondente al minimo

# Grafico 
plot(a, varTest, type = "b",
     xlab = "a", ylab = expression(sum((x[i]-a)^2/n)),
     main = "La media minimizza la somma degli scarti quadratici")
points(a[min.a], varTest[min.a], pch = 20, col = "red", cex = 1.5)


## Altre medie analitiche ----

# Media geometrica
Mg <- prod(x)^(1/n)
Mg

# Media armonica
Ma <- n / sum(1/x)
Ma

# Media di potenza (r = 2)
r  <- 2
M2 <- (sum(x^r) / n)^(1/r)
M2

# Verifica ordinamento: Mar <= Mg <= M <= M2
Ma <= Mg & Mg <= M & M <= M2

# Altra proprietà 
M2^2 >= M^2

## Quantili e medie di posizione ----

# Le tre convenzioni per i quantili empirici
cbind(sort(x),
      CONV1 = 1:n / n, 
      CONV2 = 0:(n-1) / n,
      CONV3 = (0.5:(n-0.5)) / n)


## Mediana ----

median(x)
quantile(x, probs = 0.5)

# Calcolo manuale
xs <- sort(x)
ifelse(n %% 2 == 0, (xs[n/2] + xs[n/2+1])/2, xs[(n+1)/2])

# Me minimizza la somma degli scarti assoluti
a       <- sort(c(seq(min(x), max(x), 0.5), median(x)))
AbsTest <- sapply(a, function(a) sum(abs(x - a)))
min.a   <- which.min(AbsTest)

# Grafico
plot(a, AbsTest, type = "b",
     xlab = "a", ylab = expression(sum(abs(x[i]-a))),
     main = "La mediana minimizza la somma degli scarti assoluti")
points(a[min.a], AbsTest[min.a], pch = 20, col = "red", cex = 1.5)


## Quartili, decili e quantili ----

Q1 <- quantile(x, probs = 0.25); Q1
Q3 <- quantile(x, probs = 0.75); Q3

quantile(x, probs = seq(0, 1, 0.1))          # decili
quantile(x, probs = c(0.025, 0.975))         # percentili 2.5 e 97.5


## Five numbers summary ----

fivenum(x)
summary(x)

# NA rimossi automaticamente
## Media sfrondata ----

mean(x, trim = 0.1)       # esclude 10% per lato
mean(sort(x)[10:83])      # equivalente manuale
mean(x)                   # confronto con media ordinaria


# Indici di dispersione ----

## Varianza e SQM ----
var(x)
sqrt(var(x)) 
sd(x)

# Per fumatori (SMOKES==1) e non fumatori (SMOKES==2)
var(x[dati$SMOKES == 1])
sd(x[dati$SMOKES == 1])

var(x[dati$SMOKES == 2])
sd(x[dati$SMOKES == 2])

mean(x[dati$SMOKES == 1])
mean(x[dati$SMOKES == 2])


## Scarto interquartile, MAD, coefficiente di variazione ----

SI  <- c(IQR(x[dati$SMOKES == 1]),
         IQR(x[dati$SMOKES == 2]))

MAD <- c(mad(x[dati$SMOKES == 1], constant = 1),
         mad(x[dati$SMOKES == 2], constant = 1))

# Verifica calcolo MAD
mad(x, constant = 1)
1 * median(abs(x - median(x)))

# Coefficiente di variazione
CV  <- c(sd(x[dati$SMOKES == 1]) / mean(x[dati$SMOKES == 1]),
         sd(x[dati$SMOKES == 2]) / mean(x[dati$SMOKES == 2]))

rbind(SI, MAD, CV)


# Simmetria e curtosi ----

## Decili come strumento esplorativo ----

quantile(x,           probs = seq(0, 1, 0.1), digits = 2)
quantile(dati$WEIGHT, probs = seq(0, 1, 0.1), digits = 2)


## Indice di Galton e indice K ----

Qh <- fivenum(x,           na.rm = TRUE)
Qw <- fivenum(dati$WEIGHT, na.rm = TRUE)

G_HEIGHT <- (Qh[4] + Qh[2] - 2*Qh[3]) / (Qh[4] - Qh[2]); G_HEIGHT
G_WEIGHT <- (Qw[4] + Qw[2] - 2*Qw[3]) / (Qw[4] - Qw[2]); G_WEIGHT

Dh <- quantile(x,           probs = seq(0, 1, 0.1))
Dw <- quantile(dati$WEIGHT, probs = seq(0, 1, 0.1))

K_HEIGHT <- (Dh[10] + Dh[2] - 2*Dh[6]) / (Dh[10] - Dh[2]); K_HEIGHT
K_WEIGHT <- (Dw[10] + Dw[2] - 2*Dw[6]) / (Dw[10] - Dw[2]); K_WEIGHT


## Indice di Pearson e indice di simmetria gamma ----

3*(mean(x) - median(x)) / sd(x)
3*(mean(dati$WEIGHT) - median(dati$WEIGHT)) / sd(dati$WEIGHT)

library(moments)
skewness(x)
skewness(dati$WEIGHT)


## Curtosi ----

kurtosis(x, na.rm = TRUE)


# Moda ----

table(dati$ACTIVITY)

# Classe modale per variabile continua
x_clas <- cut(x, breaks = c(61, 65, 69, 72, 75), include.lowest = TRUE)
table(x_clas)

library(MASS)
data(geyser)

# Esempio di distribuzione bimodale
hist(geyser$waiting,
     breaks = c(43, seq(48, 108, by = 5)),
     prob   = TRUE,
     main   = "Tempi di attesa fra eruzioni (geyser)",
     xlab   = "Minuti", col = "lightblue", border = "white")


# Densità kernel ----

xx <- c(148, 155, 162, 168, 172, 174, 175, 178, 186)
xr <- seq(135, 200, 0.5)

## Nucleo rettangolare ----

rett <- function(x, xx, h = 1) {
  n <- length(xx)
  r <- 0
  for (i in 1:n) r <- r + (abs(x - xx[i]) <= h/2) / (n*h)
  r
}

plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
     pch = 4, cex = 1.5, main = "Nucleo rettangolare (h = 10)",
     xlab = "x", ylab = "densità")
abline(h = 0)

for (i in seq_along(xx)) {
  yc <- jitter(1/180)
  segments(xx[i]-10, yc, xx[i]+10, yc, col = i+1, lwd = 3)
  segments(xx[i]-10, 0,  xx[i]-10, yc, col = i+1, lwd = 2, lty = 2)
  segments(xx[i]+10, 0,  xx[i]+10, yc, col = i+1, lwd = 2, lty = 2)
}
points(xr, rett(xr, xx, 20), lwd = 2.5, type = "s")

## Nucleo gaussiano ----

nuclg <- function(x, xx, h = 1)
  sapply(x, function(xi) sum(dnorm(xi, xx, h) / length(xx)))

plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
     pch = 4, cex = 1.5, main = "Nucleo gaussiano (h = 4)",
     xlab = "x", ylab = "densità")
abline(h = 0)
for (i in seq_along(xx))
  curve(dnorm(x, xx[i], 4) / length(xx), col = i+1, lwd = 2, lty = 2, add = TRUE)
lines(xr, nuclg(xr, xx, 4), lwd = 2.5)


## Effetto del parametro di lisciamento ----

par(mfrow = c(1, 3))
for (h in c(2.5, 4, 6)) {
  plot(xx, rep(0, 9), ylim = c(0, 0.065), xlim = c(135, 200),
       pch = 4, cex = 1.5, main = paste("h =", h), xlab = "x", ylab = "densità")
  abline(h = 0)
  for (i in seq_along(xx))
    curve(dnorm(x, xx[i], h)/length(xx), col = i+1, lwd = 2, lty = 2, add = TRUE)
  lines(xr, nuclg(xr, xx, h), lwd = 2.5)
}
par(mfrow = c(1, 1))


## density() in R ----

den <- density(geyser$waiting)
plot(den, main = paste("Tempi di attesa eruzioni — bw =", round(den$bw, 2)))

par(mfrow = c(1, 3))
plot(density(geyser$waiting, bw =  1), main = "bw = 1 (troppo piccolo)")
plot(density(geyser$waiting, bw =  4), main = "bw = 4 (adeguato)")
plot(density(geyser$waiting, bw = 10), main = "bw = 10 (troppo grande)")
par(mfrow = c(1, 1))


# Trasformazione delle variabili ----

## Standardizzazione ----

z <- (x - mean(x)) / sd(x)
round(mean(z), 10)   # → 0
round(var(z),  10)   # → 1


## Trasformazione logaritmica ----

# La media del log NON è il log della media (disuguaglianza di Jensen)
log(mean(dati$WEIGHT))
mean(log(dati$WEIGHT))