# T² DE HOTELLING PARA DOS SENSORES INDEPENDIENTES
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")

x_bar_1 <- colMeans(datos_sensor1)
x_bar_2 <- colMeans(datos_sensor2)
S_1 <- cov(datos_sensor1)
S_2 <- cov(datos_sensor2)

# 1. Test de normalidad multivariada (Royston)
shapiro_1 <- mvn(data = datos_sensor1, mvnTest = "royston", multivariatePlot = "none")
shapiro_2 <- mvn(data = datos_sensor2, mvnTest = "royston", multivariatePlot = "none")
normalidad_1 <- shapiro_1$multivariateNormality$`p value` > 0.05
normalidad_2 <- shapiro_2$multivariateNormality$`p value` > 0.05
cat(sprintf("Normalidad Sensor 1: p = %.4f → %s\n",
    shapiro_1$multivariateNormality$`p value`,
    ifelse(normalidad_1, "NO RECHAZAR", "RECHAZAR")))
cat(sprintf("Normalidad Sensor 2: p = %.4f → %s\n\n",
    shapiro_2$multivariateNormality$`p value`,
    ifelse(normalidad_2, "NO RECHAZAR", "RECHAZAR")))

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

# 3. T² de Hotelling
d <- x_bar_1 - x_bar_2
alpha <- 0.05

# Caso varianzas iguales (Pooled)
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
gl1_p <- p; gl2_p <- n1 + n2 - p - 1
p_valor_pooled <- 1 - pf(F_pooled, gl1_p, gl2_p)

# Caso varianzas desiguales (Nel-Van der Merwe)
S_combined <- (S_1/n1) + (S_2/n2)
T2_unpooled <- as.numeric(t(d) %*% solve(S_combined) %*% d)
A <- S_1/n1; B <- S_2/n2
nu_hat <- sum(diag((A+B) %*% solve(A+B)))^2 /
  (sum(diag(A %*% solve(A+B)))^2/(n1-1) + sum(diag(B %*% solve(A+B)))^2/(n2-1))
F_unpooled <- T2_unpooled / p
p_valor_unpooled <- 1 - pf(F_unpooled, p, nu_hat)

cat("Varianzas Iguales (Pooled):\n")
cat(sprintf("  T² = %.4f, F(%d,%d) = %.4f, p = %.6f\n", T2_pooled, gl1_p, gl2_p, F_pooled, p_valor_pooled))
cat("\nVarianzas Desiguales (Nel-Van der Merwe):\n")
cat(sprintf("  T² = %.4f, F(%d,%.2f) = %.4f, p = %.6f\n", T2_unpooled, p, nu_hat, F_unpooled, p_valor_unpooled))

# Decisión final según Test de Box
if(varianzas_iguales) {
  cat(sprintf("\nRecomendación: Usar Pooled → %s\n",
      ifelse(F_pooled > qf(0.95, gl1_p, gl2_p), "RECHAZAR H0", "NO RECHAZAR H0")))
} else {
  cat(sprintf("\nRecomendación: Usar Nel-Van der Merwe → %s\n",
      ifelse(F_unpooled > qf(0.95, p, nu_hat), "RECHAZAR H0", "NO RECHAZAR H0")))
}

# Visualización
estadisticos_comp <- data.frame(
  Banda = rep(c("Rojo (650 nm)", "Verde (550 nm)", "NIR (850 nm)"), 2),
  Sensor = rep(c("Sensor 1", "Sensor 2"), each = 3),
  Media = c(x_bar_1, x_bar_2),
  SE = c(sqrt(diag(S_1)/n1), sqrt(diag(S_2)/n2))
)
estadisticos_comp$Banda <- factor(estadisticos_comp$Banda,
  levels = c("Rojo (650 nm)", "Verde (550 nm)", "NIR (850 nm)"))

ggplot(estadisticos_comp, aes(x = Banda, y = Media, fill = Sensor)) +
  geom_col(position = position_dodge(0.8), alpha = 0.8, width = 0.7) +
  geom_errorbar(aes(ymin = Media - 1.96*SE, ymax = Media + 1.96*SE),
                position = position_dodge(0.8), width = 0.3, size = 1) +
  scale_fill_manual(values = c("Sensor 1" = "#2E86AB", "Sensor 2" = "#A23B72")) +
  labs(title = "Comparación de Reflectancia entre Sensores",
       x = "Banda Espectral", y = "Reflectancia Media", fill = "") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        legend.position = "top")