ash

ash.R

library(MASS)
waiting <- geyser$waiting
n <- length(waiting)
m <- 14
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")