Incorrect computation of acceleration parameter in BCA intervals

[bwiernik]

Describe the bug
Currently, the BCa bootstrap intervals computed by bcai() compute the acceleration parameter from the bootstrap estimates and the mean estimate.

bayestestR/R/bci.R

Line 277 in 603b91d

U <- (sims - 1) * (mean(x, na.rm = TRUE) - x)

This isn't consistent with the definition given by https://projecteuclid.org/journals/statistical-science/volume-11/issue-3/Bootstrap-confidence-intervals/10.1214/ss/1032280214.full (equations 6.5, 6.6, 6.7), which computes this parameter based on influence values/leave-one-out/jackknife values for the data points (not bootstrap replicates).

By comparison, boot::bca.ci() computes the a parameter using influence functions or jackknife procedures.

bca.ci <- function (boot.out, conf = 0.95, index = 1, t0 = NULL, t = NULL, 
  L = NULL, h = function(t) t, hdot = function(t) 1, hinv = function(t) t, 
  ...) 
{
  t.o <- t
  if (is.null(t) || is.null(t0)) {
    t <- boot.out$t[, index]
    t0 <- boot.out$t0[index]
  }
  t <- t[is.finite(t)]
  w <- qnorm(sum(t < t0)/length(t))
  if (!is.finite(w)) 
    stop("estimated adjustment 'w' is infinite")
  alpha <- (1 + c(-conf, conf))/2
  zalpha <- qnorm(alpha)
  if (is.null(L)) 
    L <- empinf(boot.out, index = index, t = t.o, ...)
  a <- sum(L^3)/(6 * sum(L^2)^1.5)
  if (!is.finite(a)) 
    stop("estimated adjustment 'a' is NA")
  adj.alpha <- pnorm(w + (w + zalpha)/(1 - a * (w + zalpha)))
  qq <- norm.inter(t, adj.alpha)
  cbind(conf, matrix(qq[, 1L], ncol = 2L), matrix(hinv(h(qq[, 
    2L])), ncol = 2L))
}

Thus, the results of bcai() are currently not correct, as they are not correctly computing the acceleration parameter.

cf. https://stats.stackexchange.com/questions/590257/using-the-bootstrap-estimates-to-compute-the-acceleration-parameter-for-bca-boot