ash cdf

ash_cdf.R

library(MASS)
waiting <- geyser$waiting

n <- length(waiting)
print(n)
m <- 3
a <- min(waiting)
b <- max(waiting)
h <- 7.27037
delta <- h / m

# get the bin counts
br <- seq(a - delta * m, b + 2 * delta * m, delta)
histg <- hist(waiting, breaks = br, plot = FALSE)
nk <- histg$counts
K <- abs((1 - m):(m - 1))

fhat <- function(x) {
    i <- max(which(x > br))
    k <- (i - m + 1):(i + m - 1)
    vk <- nk[k]
    sum((1 - K / m) * vk) / (n * h)
}

vec_z <- as.matrix(seq(a, b + h, .1))
f_ash <- apply(vec_z, 1, fhat)

br2 <- seq(a, b + h, h)
# hist(waiting, breaks = br2, freq = FALSE, main = "ASH", ylim = c(0, max(f_ash)))
# lines(vec_z, f_ash, xlab = "waiting", ylab = "density", col = "tomato")

cdf_ash <- cumsum(f_ash * 0.1) # Assuming a uniform step of 0.1 in vec_z

# Plot the CDF
# Compute the estimated CDF
cdf_ash <- cumsum(f_ash * 0.1) # Assuming a uniform step of 0.1 in vec_z

# Compute the real ECDF of your data
ecdf_data <- ecdf(waiting)

# Create a vector of values at which to evaluate the real ECDF
vec_x <- seq(a, b, length.out = length(vec_z))

# Evaluate the real ECDF at those values
ecdf_values <- ecdf_data(vec_x)

# Plot the estimated CDF and the real ECDF
plot(vec_z, cdf_ash, type = "l", col = "blue", xlab = "waiting", ylab = "CDF", main = "Estimated CDF vs. Real ECDF")
lines(vec_x, ecdf_values, col = "red") # Plot the real ECDF in red
legend("topright", legend = c("Estimated CDF", "Real ECDF"), col = c("blue", "red"), lty = 1)