Ajuste de Distribuciones a Datos Climáticos

Introducción

En este documento, ajustamos diferentes distribuciones a un conjunto de datos climáticos utilizando la librería fitdistrplus en R. Luego comparamos las distribuciones utilizando el criterio de información de Akaike (AIC) y visualizamos los resultados.

Dashboard Datos Climaticos

https://ayrgthon.github.io/dcp.html

Cargar los datos

## Loading required package: MASS

## Loading required package: survival

# Seleccionar columnas relevantes de datos climáticos
climatic_data <- Datos_aves_CV[, c("t_max", "solar.radiation", "speed.wind", "t_min", "vapr", "prec", "tavg")]

head(climatic_data)

##   t_max solar.radiation speed.wind t_min vapr prec tavg
## 1  32.4           17862        2.6  23.1 2.88  213 27.8
## 2  32.3           16833        2.6  23.0 2.87  177 27.6
## 3  32.2           17810        2.5  22.7 2.85  276 27.5
## 4  32.2           16813        2.6  22.9 2.86  178 27.5
## 5  32.3           17822        2.6  22.8 2.86  328 27.5
## 6  32.3           17650        2.6  22.8 2.87  192 27.6

library(fitdistrplus)

fit_and_summarize <- function(X) {
  # Función interna para ajuste con manejo de errores
  safe_fit <- function(dist_name) {
    tryCatch({
      fitdist(X, dist_name)
    }, error = function(e) {
      NULL
    })
  }
  
  fw <- safe_fit("weibull")
  fg <- safe_fit("gamma")
  fln <- safe_fit("lnorm")
  fex <- safe_fit("exp")
  fca <- safe_fit("cauchy")
  flo <- safe_fit("logis")
  
  fits <- list(Weibull = fw, Gamma = fg, Lognormal = fln, Exponencial = fex, Cauchy = fca, Logistica = flo)
  fits <- fits[!sapply(fits, is.null)]
  
  aic_values <- sapply(fits, function(fit) fit$aic)
  df_aic <- data.frame(Distribucion = names(aic_values), AIC = aic_values)
  df_aic_ordenado <- df_aic[order(df_aic$AIC), ]
  
  parameters <- lapply(fits, function(fit) fit$estimate)
  
  list(
    summary = df_aic_ordenado,
    parameters = parameters
  )
}

results <- lapply(climatic_data, fit_and_summarize)

## <simpleError in optim(par = vstart, fn = fnobj, fix.arg = fix.arg, obs = data,     gr = gradient, ddistnam = ddistname, hessian = TRUE, method = meth,     lower = lower, upper = upper, ...): non-finite finite-difference value [1]>

Resultados para temperatura maxima

# Mostrar resultados para 't_max'
print(results$t_max$summary)

##             Distribucion      AIC
## Cauchy            Cauchy 21521.28
## Weibull          Weibull 21927.63
## Logistica      Logistica 22575.88
## Gamma              Gamma 24390.66
## Lognormal      Lognormal 24645.11
## Exponencial  Exponencial 55773.55

print(results$t_max$parameters)

## $Weibull
##    shape    scale 
## 27.38498 32.67238 
## 
## $Gamma
##     shape      rate 
## 351.58635  10.98412 
## 
## $Lognormal
##    meanlog      sdlog 
## 3.46458873 0.05445697 
## 
## $Exponencial
##       rate 
## 0.03124138 
## 
## $Cauchy
##   location      scale 
## 32.1646065  0.4641085 
## 
## $Logistica
##   location      scale 
## 32.1837099  0.7864144

Weibull: La distribución Weibull tiene un parámetro de forma (shape) de 27.38498 y un parámetro de escala (scale) de 32.67238 Esto sugiere que la variable t_max sigue una distribución con una forma muy pronunciada y centrada alrededor de 32.67238.

Cauchy: La distribución Cauchy tiene una ubicación (location) de 32.1646065 y una escala (scale) de 0.4641085. La Cauchy, conocida por sus colas largas, sugiere que t_max puede tener valores extremos con más frecuencia que otras distribuciones.

Logística: La distribución Logística tiene una ubicación (location) de 32.1837099 y una escala (scale) de 0.7864144, indicando que t_max sigue una distribución simétrica alrededor de 32.18, con una dispersión moderada.

# Configurar el layout para gráficos
par(mfrow=c(2,2))
plot.legend <- c("Weibull", "Cauchy", "Logistica")
X <- climatic_data$t_max

# Ajustar distribuciones a 't_max'
fw <- fitdist(X, "weibull")
fln <- fitdist(X, "lnorm")
fg <- fitdist(X, "gamma")

# Comparación de densidades
denscomp(list(fw, fln, fg), legendtext=plot.legend)

# Gráfico QQ
qqcomp(list(fw, fln, fg), legendtext=plot.legend)

# Comparación de CDF
cdfcomp(list(fw, fln, fg), legendtext=plot.legend)

# Comparación de PP
ppcomp(list(fw, fln, fg), legendtext=plot.legend)