Bootstrap your own statistic

Runs the ported Numerics bootstrap against resampling, fitting and statistic functions you write. This is the class upstream exposes as four delegates – ResampleFunction, FitFunction, StatisticFunction and JackknifeFunction – so any quantity you can compute from a fitted parameter set can be given a confidence interval, not just the built-in distribution quantiles bootstrap_analysis() covers.

Usage

bootstrap_custom(
  data,
  resample,
  fit = NULL,
  statistic,
  jackknife = NULL,
  replicates = 1000,
  alpha = 0.1,
  ci_method = "Percentile",
  seed = 12345,
  parameters = NULL,
  inner_replicates = NULL,
  max_retries = NULL,
  run_type = "regular",
  fit_with_covariance = NULL,
  original_covariance = NULL,
  pivotal_links = NULL,
  pivotal_invalid_draw_policy = "drop",
  regularize_pivotal_covariances = TRUE,
  pivotal_z_limit = NULL,
  add_pivotal_jitter = FALSE,
  pivotal_jitter_scale = 0.01
)

Arguments

data: the original sample: a non-empty numeric vector.
resample, statistic: the two always-required functions, with the signatures above.
fit: the fitting function, required by run_type = "regular" and unused – and so refused – by run_type = "pivotal", which fits through fit_with_covariance instead.
jackknife: the leave-one-out function, required by ci_method = "BCa" and unused by every other method. NULL by default.
replicates: number of bootstrap replicates.
alpha: the interval's total tail probability: 0.1 gives a 90% interval, alpha / 2 in each tail.
ci_method: one of "Percentile" (the default), "BiasCorrected", "Normal", "BootstrapT" or "BCa". "Normal" and "BootstrapT" work on the ported cube-root transform of the statistic; "BootstrapT" runs the studentized workflow, which nests inner_replicates further resample-and-fit pairs inside every replicate.
seed: PRNG seed; 12345 is the C# default.
parameters: the original parameter vector the replicates are compared against. NULL, the default, uses fit(data), which is what the bootstrap ordinarily means by it.
inner_replicates: inner replicates for ci_method = "BootstrapT", ignored by every other method. NULL leaves the ported default (300) in force.
max_retries: the maximum number of times a single failed replicate is retried before it is counted in failed_replicates. NULL leaves the ported default (MaxRetries, 20) in force. Each retry is another crossing into R, so lowering it caps the worst case rather than changing the typical one.
run_type: "regular" (the default) or "pivotal"; see the section above.
fit_with_covariance: the covariance-aware fitting function run_type = "pivotal" fits through, with the signature above. Required by that run type and refused by the other.
original_covariance: the parent fit's covariance matrix, one row and one column per parameter. NULL, the default, takes it from fit_with_covariance(data). Pivotal only.
pivotal_links: one entry per parameter, each a link function name ("Identity", "Log", "Logit", "Probit", "ComplementaryLogLog", "YeoJohnson" or "FisherZ") or NULL for the identity, given as a list. The standardization happens in link space, so "Log" on a scale parameter keeps every reinflated draw positive. NULL, the default, is the identity throughout. Pivotal only.
pivotal_invalid_draw_policy: what to do with a draw the transform could not produce (a non-finite standardized vector, or one outside pivotal_z_limit): "drop" (the default, leaving it out of the ensemble), "use_raw" (keeping the untransformed fit) or "use_parent" (substituting the original fit). Pivotal only.
regularize_pivotal_covariances: whether each link-space covariance is made symmetric positive definite before it is factored. TRUE by default, as upstream. Pivotal only.
pivotal_z_limit: an absolute limit on every component of the standardized vector, beyond which the draw is invalid and the policy above applies. NULL, the default, is no limit. Pivotal only.
add_pivotal_jitter, pivotal_jitter_scale: whether to add Gaussian jitter to the standardized vector before the limit is checked, and its base standard deviation (the applied scale is this divided by the square root of the parameter count). FALSE and 0.01 by default, as upstream. Pivotal only.

Value

A list with, per statistic, estimate (the statistic of parameters, not a bootstrap average), lower, upper, standard_error, mean (the mean over valid replicates, so mean - estimate is the bias estimate) and valid_count; the same first three for the fitted parameters as parameter_estimate, parameter_lower and parameter_upper; and replicates, failed_replicates, alpha, ci_method and run_type. A pivotal run adds pivotal_diagnostics (a list of the six replicate counts) and the raw block: raw_estimate, raw_lower, raw_upper, raw_standard_error, raw_mean, raw_valid_count, raw_parameter_estimate, raw_parameter_lower and raw_parameter_upper.

The four functions

Getting an argument order wrong is the likeliest mistake here, and the C++ side cannot tell a swapped pair from a deliberate one, so each signature is given exactly:

resample(data, parameters, rng): Returns one bootstrap sample. data is the original sample, parameters is the current parameter vector, and rng is a handle on THIS replicate's generator – draw from it with rng_uniform() and rng_integers(), not with stats::runif() or base::sample(), or the seeded run stops being reproducible and stops agreeing with Python. The returned sample need not be the same length as data.
fit(data): Returns the parameter vector fitted to data: one number per parameter, the same count every time.
statistic(parameters): Returns the numbers to put intervals on, computed from a fitted parameter vector: one or more, the same count every time.
jackknife(data, index): Returns data with observation index left out. index counts from 0, matching the ported delegate, so the R spelling is data[-(index + 1)] – the naive data[-index] is wrong for every value of index, not just one: at index = 0 it is data[-0], which R evaluates to numeric(0), the EMPTY vector, and that is refused by name; at every later index it is the right LENGTH but drops the wrong observation (one off), which nothing can catch. Only the "BCa" method uses it; every other method ignores it.
fit_with_covariance(data): Returns list(parameters = , covariance = ): the parameter vector fitted to data and its covariance matrix, one row and one column per parameter. Only run_type = "pivotal" uses it, and that run type uses it INSTEAD of fit.

The pivotal run type

run_type = "pivotal" is upstream's other bootstrap mode. Rather than treating each resampled fit as a draw from the sampling distribution, it standardizes the fit against the original one through the resample's OWN covariance and reinflates it through the original's, so a replicate fitted on an unusually flat likelihood contributes an appropriately smaller step. It therefore needs a covariance with every fit, which is what fit_with_covariance supplies in place of fit.

The original fit is the parent every draw is compared against. Left alone it is fit_with_covariance(data); parameters replaces its parameter vector and original_covariance its covariance, either independently of the other.

The result gains pivotal_diagnostics – the six replicate counts the run kept, from requested_replicates down to retained_pivotal_replicates – and a second interval block, raw_lower/raw_upper and companions, which is the plain percentile interval of the RAW covariance-aware fits before the transform. Comparing the two is the point of reporting both. Only ci_method = "Percentile" exists after a pivotal run, and asking for another is refused before the first replicate rather than after all of them.

How many times your functions are called

replicates calls of each of resample, fit and statistic, plus one extra statistic call to learn how many values it returns and one fit call when parameters is not supplied. A failed replicate is retried up to max_retries times (20 by default), and "BCa" adds one jackknife + fit + statistic per observation. "BootstrapT" is the expensive one: it multiplies the resample and fit counts by inner_replicates, so the ported defaults (10,000 x 300) would be three million crossings back into R. Start small. A pivotal run calls fit_with_covariance in place of fit, once per replicate plus once up front for the parent fit, and statistic once per retained draw plus once per raw fit – so about twice as often as a regular run of the same length.

Reproducing a run in Python

The draws come from the core's seeded Mersenne Twister, so an identical bootstrap_custom() call in Python resamples the identical observations. The numbers your functions compute from them are your own R code, though, and R and Python do not guarantee identical rounding for the same formula: arithmetic (+ - * /) is IEEE-deterministic and does reproduce, while log, exp and friends come from each platform's math library. Note also that R's sum() and mean() accumulate in extended precision where Python does not, so an explicit loop is the portable spelling.

Examples

# \donttest{
x <- c(4.1, 5.2, 4.8, 5.5, 4.9, 5.1, 5.3, 4.7)
# An ordinary iid bootstrap of the mean. rng_integers() draws on [0, n), counting from 0, so
# the index is shifted by one for R.
res <- bootstrap_custom(
  data = x,
  resample = function(data, parameters, rng) {
    data[rng_integers(rng, length(data), 0, length(data)) + 1L]
  },
  fit = function(data) {
    acc <- 0
    for (xi in data) acc <- acc + xi
    acc / length(data)
  },
  statistic = function(parameters) parameters,
  replicates = 500, seed = 12345
)
c(res$lower, res$estimate, res$upper)
#> [1] 4.699375 4.950000 5.163125

# The pivotal run type. It fits through `fit_with_covariance` instead of `fit`: a Normal
# location-scale MLE here, whose covariance is diag(s2 / n, s2 / 2n) in closed form. The
# "Log" link standardizes the scale parameter in log space, keeping every draw positive.
fit_with_cov <- function(data) {
  n <- length(data)
  mu <- sum(data) / n
  s2 <- sum((data - mu)^2) / n
  list(parameters = c(mu, sqrt(s2)),
       covariance = matrix(c(s2 / n, 0, 0, s2 / (2 * n)), nrow = 2L, ncol = 2L))
}
piv <- bootstrap_custom(
  data = x,
  resample = function(data, parameters, rng) {
    data[rng_integers(rng, length(data), 0, length(data)) + 1L]
  },
  statistic = function(parameters) parameters,
  fit_with_covariance = fit_with_cov,
  run_type = "pivotal",
  pivotal_links = list(NULL, "Log"),
  replicates = 200, seed = 12345
)
# Two interval blocks: the pivotal ensemble, and the raw fits it was built from.
rbind(pivotal = c(piv$lower[1], piv$upper[1]), raw = c(piv$raw_lower[1], piv$raw_upper[1]))
#>            [,1]     [,2]
#> pivotal 4.56776 5.140545
#> raw     4.75000 5.175000
unlist(piv$pivotal_diagnostics)
#>        requested_replicates     rejected_raw_replicates 
#>                         200                           0 
#>       failed_raw_replicates     accepted_raw_replicates 
#>                           0                         200 
#>  invalid_pivotal_replicates retained_pivotal_replicates 
#>                           0                         200 
# }