# T² CHART DE HOTELLING CON DETECCIÓN DE OUTLIERS MULTIVARIADOS
library(MASS)
library(ggplot2)
library(MVN)
library(biotools)

set.seed(2024)

# Parámetros
n1 <- 50; n2 <- 50; p <- 3
mu_1 <- c(0.30, 0.25, 0.45)
mu_2 <- c(0.32, 0.24, 0.46)

Sigma_1 <- matrix(c(
  0.0015,  0.0008,  0.0010,
  0.0008,  0.0012, -0.0003,
  0.0010, -0.0003,  0.0020
), nrow = 3, byrow = TRUE)

Sigma_2 <- matrix(c(
  0.0018,  0.0009,  0.0011,
  0.0009,  0.0014, -0.0004,
  0.0011, -0.0004,  0.0022
), nrow = 3, byrow = TRUE)

datos_sensor1 <- mvrnorm(n = n1, mu = mu_1, Sigma = Sigma_1)
colnames(datos_sensor1) <- c("Rojo_650nm", "Verde_550nm", "NIR_850nm")
datos_sensor2 <- mvrnorm(n = n2, mu = mu_2, Sigma = Sigma_2)
colnames(datos_sensor2) <- c("Rojo_650nm", "Verde_550nm", "NIR_850nm")

# Función para T² Chart (detección de outliers)
detectar_outliers_T2 <- function(datos, alpha = 0.05) {
  n <- nrow(datos); p <- ncol(datos)
  x_bar <- colMeans(datos); S <- cov(datos); S_inv <- solve(S)
  T2_values <- numeric(n)
  for(i in 1:n) {
    diff <- datos[i, ] - x_bar
    T2_values[i] <- as.numeric(t(diff) %*% S_inv %*% diff)
  }
  F_crit <- qf(1 - alpha, df1 = p, df2 = n - p)
  UCL <- ((p * (n - 1) * (n + 1)) / (n * (n - p))) * F_crit
  outliers_idx <- which(T2_values > UCL)
  data.frame(Obs = 1:n, T2 = T2_values, UCL = UCL, Es_Outlier = T2_values > UCL)
}

# Detectar outliers en ambos sensores
res_s1 <- detectar_outliers_T2(datos_sensor1)
res_s2 <- detectar_outliers_T2(datos_sensor2)

cat(sprintf("Sensor 1: %d outliers de %d (UCL=%.4f)\n",
    sum(res_s1$Es_Outlier), n1, res_s1$UCL[1]))
cat(sprintf("Sensor 2: %d outliers de %d (UCL=%.4f)\n\n",
    sum(res_s2$Es_Outlier), n2, res_s2$UCL[1]))

# Normalidad multivariada
shapiro_1 <- mvn(data = datos_sensor1, mvnTest = "royston", multivariatePlot = "none")
shapiro_2 <- mvn(data = datos_sensor2, mvnTest = "royston", multivariatePlot = "none")
cat(sprintf("Normalidad Sensor 1: p = %.4f\n", shapiro_1$multivariateNormality$`p value`))
cat(sprintf("Normalidad Sensor 2: p = %.4f\n\n", shapiro_2$multivariateNormality$`p value`))

# Test de Box
datos_combined <- rbind(datos_sensor1, datos_sensor2)
grupo <- factor(rep(c("S1", "S2"), c(n1, n2)))
box_test <- boxM(datos_combined, grupo)
varianzas_iguales <- box_test$p.value > 0.05
cat(sprintf("Test de Box: p = %.4f → %s\n\n", box_test$p.value,
    ifelse(varianzas_iguales, "Varianzas iguales", "Varianzas desiguales")))

# T² de Hotelling dos muestras
x_bar_1 <- colMeans(datos_sensor1); x_bar_2 <- colMeans(datos_sensor2)
S_1 <- cov(datos_sensor1); S_2 <- cov(datos_sensor2)
d <- x_bar_1 - x_bar_2

S_pooled <- ((n1-1)*S_1 + (n2-1)*S_2) / (n1+n2-2)
T2_pooled <- as.numeric((n1*n2)/(n1+n2) * t(d) %*% solve(S_pooled) %*% d)
F_pooled <- ((n1+n2-p-1)/((n1+n2-2)*p)) * T2_pooled
p_val_pooled <- 1 - pf(F_pooled, p, n1+n2-p-1)
cat(sprintf("T² Pooled: F = %.4f, p = %.6f\n", F_pooled, p_val_pooled))

# T² Chart - Gráfico Sensor 1
ggplot(res_s1, aes(x = Obs, y = T2)) +
  geom_line(color = "#2E86AB", size = 1) +
  geom_point(aes(color = Es_Outlier), size = 3) +
  geom_hline(yintercept = res_s1$UCL[1], color = "red", linetype = "dashed", size = 1.5) +
  scale_color_manual(values = c("FALSE" = "#2E86AB", "TRUE" = "#e74c3c"),
                     labels = c("Normal", "Outlier")) +
  annotate("text", x = n1*0.9, y = res_s1$UCL[1]*1.1,
           label = sprintf("UCL = %.2f", res_s1$UCL[1]), color = "red", fontface = "bold") +
  labs(title = "T² Chart - Sensor 1",
       subtitle = sprintf("Outliers: %d de %d", sum(res_s1$Es_Outlier), n1),
       x = "Observación", y = "Estadístico T²", color = "Estado") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.position = "top")