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 Wind (Viento)
Viento_Total <- as.numeric(Datos$Wind)
Viento_Total <- Viento_Total[!is.na(Viento_Total)]
# --- GRÁFICA: TODOS LOS DATOS ORIGINALES DE VIENTO ---
hist(Viento_Total, freq = FALSE, main="Distribución Total de Velocidad del Viento",
xlab="m/s", ylab="Densidad", col="blue")
# --- FILTRO CRÍTICO ---
# La distribución Gamma requiere valores estrictamente positivos (> 0)
Viento <- Viento_Total[Viento_Total > 0]
# ----------------------------------------------------------------
# GRÁFICA 107: Sin agrupación manual (Breaks automáticos de R)
histograma <- hist(Viento, freq = FALSE,
main="Gráfica 107. Modelo de probabilidad GAMMA (Viento)",
xlab="m/s", ylab="Densidad de probabilidad", col="blue")
# GRÁFICA 108: Ajuste del modelo al histograma (GAMMA)
histograma <- hist(Viento, freq = FALSE,
main="Gráfica 108. Ajuste de Densidad GAMMA (Viento)",
xlab="m/s", ylab="Densidad de probabilidad", col="blue")
h <- length(histograma$counts)
# --- AJUSTE GAMMA ---
fit_gam <- fitdist(Viento, "gamma")
shape_v <- fit_gam$estimate["shape"]
rate_v <- fit_gam$estimate["rate"]
# Curva suave GAMMA
x_seq <- seq(min(Viento), max(Viento), length.out = 500)
lines(x_seq, dgamma(x_seq, shape = shape_v, rate = rate_v), lwd=4, col="black")
# Frecuencia simple observada
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_v, rate = rate_v) -
pgamma(histograma$breaks[i], shape = shape_v, rate = rate_v))
}
Fe <- P * length(Viento)
# Comparar tamaño real y el modelo
sum(Fe)
## [1] 361.8285
n <- length(Viento)
# 2. Test de Pearson (Correlación)
plot(Fo, Fe, main="Gráfica 109: Correlación de frecuencias (Viento - 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
# 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("Wind Speed (Gamma)")
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 Viento (GAMMA)")
Tabla. Resumen de validación para Viento (GAMMA)
| Wind Speed (Gamma) |
97.49 |
19.39 |
21.03 |
# 4. Cálculo de probabilidad
# --- CAMBIO SOLICITADO: Rango 2.0 a 3.0 ---
limite_inf <- 2.0
limite_sup <- 3.0
prob_area <- pgamma(limite_sup, shape_v, rate_v) - pgamma(limite_inf, shape_v, rate_v)
x_grafica <- seq(min(Viento), max(Viento), length.out = 1000)
y_grafica <- dgamma(x_grafica, shape_v, rate_v)
# GRAFICAR MODELO FINAL
plot(x_grafica, y_grafica, type = "l", col = "skyblue3", lwd = 2,
main="Distribución GAMMA Final (Viento)",
ylab="Densidad de Probabilidad",
xlab="m/s")
# Sombreado de la probabilidad (Área entre 2.0 y 3.0)
x_sombra <- seq(limite_inf, limite_sup, length.out = 100)
y_sombra <- dgamma(x_sombra, shape_v, rate_v)
polygon(c(x_sombra, rev(x_sombra)), c(y_sombra, rep(0, length(y_sombra))),
col = rgb(1, 0, 0, 0.4), border = "red")
# Leyenda corregida
legend("topright",
legend = c("Modelo Gamma", paste("Probabilidad 2-3 m/s:", round(prob_area*100, 2), "%")),
fill = c(NA, rgb(1, 0, 0, 0.4)), border = c("skyblue3", "red"), bty = "n")
