Skip to content

Instantly share code, notes, and snippets.

@xvzftube

xvzftube/ash.R

Created September 3, 2023 01:00
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save xvzftube/3bc379b87ae8347cf9c887c3b72ecdde to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save xvzftube/3bc379b87ae8347cf9c887c3b72ecdde to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. xvzftube created this gist Sep 3, 2023.
    28 changes: 28 additions & 0 deletions ash.R
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,28 @@
    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")