# T² DE HOTELLING PARA ESPECTROSCOPÍA DE REFLECTANCIA (UNA MUESTRA)
library(MASS)
library(ggplot2)

set.seed(2024)

# Parámetros
n <- 20
p <- 3
mu_0 <- c(0.30, 0.25, 0.45)
names(mu_0) <- c("Rojo_650nm", "Verde_550nm", "NIR_850nm")

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

# Generar datos de normal trivariada
datos <- mvrnorm(n = n, mu = mu_0, Sigma = Sigma)
colnames(datos) <- names(mu_0)

# Estadísticos muestrales
x_bar <- colMeans(datos)
S <- cov(datos)
d <- x_bar - mu_0
S_inv <- solve(S)

# Estadístico T²
T2 <- as.numeric(n * t(d) %*% S_inv %*% d)

# Transformación a F
F_stat <- ((n - p) / (p * (n - 1))) * T2
gl1 <- p
gl2 <- n - p
alpha <- 0.05
F_crit <- qf(1 - alpha, df1 = gl1, df2 = gl2)
p_valor <- 1 - pf(F_stat, df1 = gl1, df2 = gl2)
decision <- ifelse(F_stat > F_crit, "RECHAZAR H0", "NO RECHAZAR H0")

cat(sprintf("T² = %.4f\n", T2))
cat(sprintf("F(%d, %d) = %.4f\n", gl1, gl2, F_stat))
cat(sprintf("F crítico = %.4f\n", F_crit))
cat(sprintf("p-valor = %.6f\n", p_valor))
cat(sprintf("Decisión: %s\n\n", decision))

# Tests univariados por banda
for(i in 1:3) {
  t_stat <- sqrt(n) * d[i] / sqrt(diag(S)[i])
  p_uni <- 2 * (1 - pt(abs(t_stat), df = n - 1))
  cat(sprintf("%s: t = %.4f, p = %.4f\n", names(mu_0)[i], t_stat, p_uni))
}

# Matriz de correlaciones
R <- cov2cor(S)
cat("\nMatriz de correlaciones:\n")
print(round(R, 4))

# Visualización
df_compare <- data.frame(
  Banda = rep(c("Rojo\n650 nm", "Verde\n550 nm", "NIR\n850 nm"), 2),
  Tipo = rep(c("Teórico", "Observado"), each = 3),
  Reflectancia = c(mu_0, x_bar),
  SE = c(0, 0, 0, sqrt(diag(S)/n))
)

ggplot(df_compare, aes(x = Banda, y = Reflectancia, fill = Tipo)) +
  geom_col(position = position_dodge(width = 0.8), alpha = 0.8, width = 0.7) +
  geom_errorbar(aes(ymin = Reflectancia - 1.96*SE, ymax = Reflectancia + 1.96*SE),
                position = position_dodge(width = 0.8), width = 0.3, size = 1) +
  scale_fill_manual(values = c("Teórico" = "#2E86AB", "Observado" = "#A23B72")) +
  labs(title = "Reflectancia Espectral: Estándar vs Observado",
       subtitle = sprintf("Test T² de Hotelling: %s (p = %.4f)", decision, p_valor),
       x = "Banda Espectral", y = "Reflectancia", fill = "") +
  theme_minimal(base_size = 14) +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
        plot.subtitle = element_text(hjust = 0.5, size = 12),
        legend.position = "top")