ash_plots_ggplot.R

library(MASS)
library(tidyverse)

waiting <- geyser$waiting

n <- length(waiting)
print(n)
m <- 16
a <- min(waiting)
b <- max(waiting)
h <- 7.27037
step <- 0.1
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, step))
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")
# use ggplot to plot the density and a histogram of the data
df_waiting <- tibble(waiting = waiting)
p <- tibble(vec_z = vec_z, f_ash = f_ash) |>
    ggplot(aes(x = vec_z, y = f_ash)) +
    geom_histogram(
        data = df_waiting, aes(x = waiting, y = after_stat(density)),
        bins = 20, color = "black", fill = "white"
    ) +
    geom_line(color = "tomato") +
    theme_minimal()

ggsave("ash_hist.jpg", p, width = 6, height = 4, dpi = 300)
# Plot the CDF
# Compute the estimated CDF
cdf_ash <- cumsum(f_ash * step) # 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
p <- tibble(vec_x = vec_x, ecdf_values = ecdf_values, cdf_ash = cdf_ash) |>
    ggplot(aes(x = vec_x, y = ecdf_values)) +
    geom_step() +
    geom_step(aes(x = vec_z, y = cdf_ash), color = "tomato") +
    # geom_point(aes(x = vec_z, y = cdf_ash), color = "black") +
    theme_minimal()

ggsave("ash.jpg", p, width = 6, height = 4, dpi = 300)