EVD gumbel

Patient_Admission_Data_3_ <- read_excel("Patient_Admission_Data (3).xlsx")
x <- na.omit(Patient_Admission_Data_3_$Number_Of_Days)
if(length(x) == 0){
  stop("Error: No numeric data found in Number_Of_Days column")
}
n <- length(x)  # sample size


fit_gumbel <- fgev(x, shape = 0)
loc_gumbel <- fit_gumbel$estimate["loc"]
scale_gumbel <- fit_gumbel$estimate["scale"]
k_gumbel <- 2  # loc + scale


fit_gev <- fgev(x)  # shape estimated freely
loc_gev <- fit_gev$estimate["loc"]
scale_gev <- fit_gev$estimate["scale"]
shape_gev <- fit_gev$estimate["shape"]
k_gev <- length(fit_gev$estimate)  # loc + scale + shape


dev.new(width = 10, height = 6)
hist(x,
     probability = TRUE,
     col = "lightgray",
     border = "white",
     main = "Histogram with Gumbel and General EVD Densities",
     xlab = "Number of Days",
     ylab = "Density")

x_range <- seq(min(x), max(x), length.out = 500)

# Densities
density_gumbel <- dgev(x_range, loc = loc_gumbel, scale = scale_gumbel, shape = 0)
density_gev    <- dgev(x_range, loc = loc_gev, scale = scale_gev, shape = shape_gev)

# Add curves
lines(x_range, density_gumbel, col = "blue", lwd = 2)
lines(x_range, density_gev, col = "red", lwd = 2, lty = 2)

legend("topright",
       legend = c("Gumbel (shape=0)", "General EVD"),
       col = c("blue", "red"),
       lwd = 2,
       lty = c(1,2))


dev.new(width = 6, height = 6)
p <- ppoints(n)
q_gumbel <- qgev(p, loc = loc_gumbel, scale = scale_gumbel, shape = 0)
q_gev    <- qgev(p, loc = loc_gev, scale = scale_gev, shape = shape_gev)

plot(q_gumbel, sort(x), main = "Q-Q Plot Comparison", xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", col="blue", pch=16)
points(q_gev, sort(x), col="red", pch=17)
abline(0,1,col="black", lwd=2)
legend("topleft", legend=c("Gumbel","General EVD"), col=c("blue","red"), pch=c(16,17))


AIC_gumbel <- AIC(fit_gumbel)
BIC_gumbel <- AIC_gumbel + (k_gumbel*log(n) - 2*k_gumbel)

AIC_gev <- AIC(fit_gev)
BIC_gev <- AIC_gev + (k_gev*log(n) - 2*k_gev)

cat("MODEL COMPARISON RESULTS\n")

## MODEL COMPARISON RESULTS

cat(sprintf("Gumbel (shape=0)     : AIC = %.2f, BIC = %.2f\n", AIC_gumbel, BIC_gumbel))

## Gumbel (shape=0)     : AIC = 1849.37, BIC = 1855.11

cat(sprintf("General EVD (shape~) : AIC = %.2f, BIC = %.2f\n", AIC_gev, BIC_gev))

## General EVD (shape~) : AIC = 1850.38, BIC = 1858.98