library(readxl)
library(knitr)
library(fitdistrplus)
## Cargando paquete requerido: MASS
## Cargando paquete requerido: survival
# 1. Carga de datos
Datos <- read_xlsx("C:/Users/LEO/Documents/Antisana/weatherdataANTISANA.csv.xlsx")

# Trabajamos con los valores originales de Precipitación
Prec_Total <- as.numeric(Datos$Precipitation)
Prec_Total <- Prec_Total[!is.na(Prec_Total)]

# --- GRÁFICA 106: TODOS LOS DATOS ORIGINALES ---
hist(Prec_Total, freq = FALSE, main="Gráfica 106. Distribución Total de Precipitación",
     xlab="mm", ylab="Densidad", col="blue")

# --- FILTRO CRÍTICO: RANGO 30 A 80 ---
Precipitacion <- Prec_Total[Prec_Total >= 30 & Prec_Total <= 80]
# ----------------------------------------------------------------

# Definimos los cortes de 30 a 86 de 8 en 8
cortes <- seq(30, 86, by = 8)

# GRÁFICA 107: Agrupación manual de 8 en 8
histograma <- hist(Precipitacion, breaks = cortes, freq = FALSE, 
                   main="Gráfica 107. Modelo de probabilidad GAMMA (Agrupación 8 en 8)",
                   xlab="mm", ylab="Densidad de probabilidad", col="blue",
                   xlim=c(30, 80))

# GRÁFICA 108: Ajuste del modelo al histograma (GAMMA)
histograma <- hist(Precipitacion, breaks = cortes, freq = FALSE, 
                   main="Gráfica 108. Ajuste de Densidad GAMMA (30-80)",
                   xlab="mm", ylab="Densidad de probabilidad", col="blue",
                   xlim=c(30, 80))
h <- length(histograma$counts)

# --- AJUSTE GAMMA ---
fit_gam <- fitdist(Precipitacion, "gamma")
shape_p <- fit_gam$estimate["shape"]
rate_p <- fit_gam$estimate["rate"]

# Curva suave GAMMA
x_seq <- seq(30, 80, length.out = 500)
lines(x_seq, dgamma(x_seq, shape = shape_p, rate = rate_p), lwd=4, col="black")

# Frecuencia simple observada (counts en los intervalos de 8mm)
Fo <- histograma$counts

# Frecuencia simple esperada usando la función pgamma
P <- c(0)
for (i in 1:h) {
  P[i] <- (pgamma(histograma$breaks[i+1], shape = shape_p, rate = rate_p) -
             pgamma(histograma$breaks[i], shape = shape_p, rate = rate_p))
}
Fe <- P * length(Precipitacion)

# Comparar tamaño real y el modelo
sum(Fe)
## [1] 60.9618
n <- length(Precipitacion)

# 2. Test de Pearson (Correlación)
plot(Fo, Fe, main="Gráfica 109: Correlación de frecuencias (Modelo Gamma)",
     xlab="Frecuencia Observada", ylab="Frecuencia esperada", col="blue3", pch=19)

abline(lm(Fe ~ Fo), col="red", lwd=2)

Correlación <- cor(Fo, Fe) * 100
Correlación
## [1] 93.03777
# TEST DE CHI-CUADRADO
grados_libertad <- (length(histograma$counts)-1)
nivel_significancia <- 0.05

x2 <- sum((Fe - Fo)^2 / Fe)
umbral_aceptacion <- qchisq(1 - nivel_significancia, grados_libertad)

# ¿Se acepta el modelo?
x2 < umbral_aceptacion
## [1] TRUE
# 3. Tabla resumen de test
Variable <- c("Precipitation (Gamma 30-80)")
tabla_resumen <- data.frame(Variable, round(Correlación, 2), round(x2, 2), round(umbral_aceptacion, 2))
colnames(tabla_resumen) <- c("Variable", "Test Pearson (%)", "Chi Cuadrado", "Umbral de aceptación")

kable(tabla_resumen, format = "markdown", caption = "Tabla. Resumen de validación para Precipitación (GAMMA)")
Tabla. Resumen de validación para Precipitación (GAMMA)
Variable Test Pearson (%) Chi Cuadrado Umbral de aceptación
Precipitation (Gamma 30-80) 93.04 7.09 12.59 
# 4. Cálculo de probabilidad
# Ejemplo: Probabilidad de lluvia mayor a 40 mm
limite_inf <- 40 
limite_sup <- 80

prob_area <- pgamma(limite_sup, shape_p, rate_p) - pgamma(limite_inf, shape_p, rate_p)

x_grafica <- seq(30, 80, length.out = 1000)
y_grafica <- dgamma(x_grafica, shape_p, rate_p)

# GRAFICAR MODELO FINAL
plot(x_grafica, y_grafica, type = "l", col = "skyblue3", lwd = 2, 
     main="Distribución GAMMA Final (Rango 30-80)",
     ylab="Densidad de Probabilidad", 
     xlab="mm",
     xlim=c(30, 80),
     xaxt="n")

# Sombreado de la probabilidad
x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
y_sombra <- dgamma(x_sombra, shape_p, rate_p)
polygon(c(x_sombra, rev(x_sombra)), c(y_sombra, rep(0, length(y_sombra))), 
        col = rgb(1, 0, 0, 0.4), border = "red")

# Eje X marcado cada 10mm
axis(1, at = seq(30, 80, by = 10), las = 2)

# Leyenda
legend("topright", 
       legend = c("Modelo Gamma", paste("Probabilidad > 40mm:", round(prob_area*100, 2), "%")), 
       fill = c(NA, rgb(1, 0, 0, 0.4)), border = c("skyblue3", "red"), bty = "n")