Dilit-d Dzeknajf

Implementiramo delete-d jackknife

dd_jack <- function(x, h, d) {
  n <- ifelse(is.vector(x), length(x), nrow(x))
  r <- n - d
  N <- choose(n, d)
  
  leave_d_out <- apply(combn(1:n, d), 2, function(s){
    h(subset(x, !(1:n %in% s)))
  })
  
  r/d*(N-1)*var(leave_d_out)/N
}

i paralelnu verziju da bude koliko toliko brze

pardd_jack <- function(x, h, d) {
  library(parallel)
  cl <- makeCluster(detectCores(), type = "FORK")
  
  n <- ifelse(is.vector(x), length(x), nrow(x))
  r <- n - d
  N <- choose(n, d)
  
  leave_d_out <- parApply(cl, combn(1:n, d), 2, function(s){
    h(subset(x, !(1:n %in% s)))
  })
  
  ret <- r/d*(N-1)*var(leave_d_out)/N
  stopCluster(cl)
  ret
}

Testiramo da je delete-1 jednak obicnom…

x <- rnorm(20)
dd_jack(x, mean, 1)
V_jack(x, mean)

Odlicno!

Implementirajmo prosto slucajno uzorkovanje za ubrzanje dd jackknife ocena

srs_dd_jack <- function(x, h, d, m) {
  n <- ifelse(is.vector(x), length(x), nrow(x))
  r <- n - d
  # each column is a sample of size r, there are m columns
  samples <- replicate(m, sample(x, r))
  h_values <- apply(samples, 2, h)
  r/d*(m-1)/m*var(h_values)
}

x <- rnorm(20)
pardd_jack(x, mean, 2)
srs_dd_jack(x, mean, 2, 50)

Dobra je aproksimacija.

Pogledajmo kako se menja greska aproksimacije ddjack za ocenu ocekivanja i disperzije normalne raspodele, za uzorak obima 20, vrednosti za d = 3 i vrednosti m iz skupa {10, 20, 30, 50, 70, 90, 120, 160, 200, 275, 400}

mvals <- c(10, 20, 30, 50, 70, 90, 120, 160, 200, 275, 400)
set.seed(825)
x <- rnorm(20)
ddj_m <- pardd_jack(x, mean, 3)
approx_m <- replicate(1000, sapply(mvals, srs_dd_jack, x = x, h = mean, d = 3))
approx_m <- rowMeans(approx_m)

ddj_v <- pardd_jack(x, var, 3)
approx_v <- replicate(1000, sapply(mvals, srs_dd_jack, x = x, h = var, d = 3))
approx_v <- rowMeans(approx_v)

library(latex2exp)
ddj_m
mm <- rbind(as.integer(mvals),
            approx_m,
            abs(approx_m - ddj_m)*10^3,
            abs(approx_m - ddj_m) / ddj_m)

plot(mvals, abs(approx_m - ddj_m)/ddj_m, xaxt = "n", type = "l",
     xlab = "vrednosti m", ylab = "relativna greška")
axis(side = 1, at = mvals, labels = mvals)

ddj_v

vv <- rbind(as.integer(mvals),
            approx_v,
            abs(approx_v - ddj_v)*10^3,
            abs(approx_v - ddj_v) / ddj_v)

plot(mvals, abs(approx_v - ddj_v)/ddj_v, xaxt = "n", type = "l",
     xlab = "vrednosti m", ylab = "relativna greška")
axis(side = 1, at = mvals, labels = mvals)

res <- rbind(mm, vv)
xtable::xtable(res, digits = 3)

Uporedimo klasican primer medijane

set.seed(406)
n <- 100
N <- 1000
library(parallel)
cl <- makeCluster(detectCores(), type = "FORK")
normal <- matrix(rnorm(n*N), nrow = N)
estimates <- list(med=median)
norm_results <- sapply(estimates, function(f) {
  res <- parApply(cl, normal, 1, function(X){
    est <- f(X)
    d1 <- dd_jack(X, f, 1)
    d2 <- dd_jack(X, f, 2)
    d3 <- srs_dd_jack(X, f, 3, 10000)
    d25 <- srs_dd_jack(X, f, 25, 10000)
    d50 <- srs_dd_jack(X, f, 50, 10000)
    c(est = est, d1 = d1, d2 = d2, d3 = d3, d25 = d25, d50 = d50)
  })
  c(est = mean(res['est',]),
    var = var(res['est',]),
    d1 = mean(res['d1',]),
    d2 = mean(res['d2',]),
    d3 = mean(res['d3',]),
    d25 = mean(res['d25',]),
    d50 = mean(res['d50',]))
})
stopCluster(cl)

round(norm_results, digits = 3)

      med
est 0.003
var 0.015
d1  0.031
d2  0.023
d3  0.026
d25 0.018
d50 0.016