Compute bootstrap variability bands

Source: R/uncertainty_quantification.R

Compute variability bands for the optimized trend filtering estimate via the sample quantiles of the bootstrap ensemble generated by bootstrap_trendfilter().

vbands(obj, level = 0.95)

Arguments

obj

A 'bootstrap_trendfilter' object.
level

The level of the pointwise variability bands. Defaults to level = 0.95.

Value

A tibble with column set c("x","lower_band","upper_band").

References

  1. Politsch et al. (2020a). Trend filtering – I. A modern statistical tool for time-domain astronomy and astronomical spectroscopy. MNRAS, 492(3), p. 4005-4018. [Publisher] [arXiv].

  2. Politsch et al. (2020b). Trend Filtering – II. Denoising astronomical signals with varying degrees of smoothness. MNRAS, 492(3), p. 4019-4032. [Publisher] [arXiv].

See also

bootstrap_trendfilter()

Examples

# Example 1: Phase-folded light curve of an eclipsing binary star system
#
# The apparent brightness over time of a star system that has two suns
# that regularly eclipse one another from our vantage point on Earth. Here,
# the time series is stacked according to the orbital period of the binary
# system, with the primary eclipse occuring at `phase = 0` and the input
# domain ranging from -0.5 to 0.5.

data("eclipsing_binary")
head(eclipsing_binary)

#> # A tibble: 6 × 3
#>    phase  flux std_err
#>    <dbl> <dbl>   <dbl>
#> 1 -0.499 0.938  0.0102
#> 2 -0.498 0.930  0.0102
#> 3 -0.496 0.944  0.0102
#> 4 -0.494 0.963  0.0102
#> 5 -0.494 0.934  0.0102
#> 6 -0.494 0.949  0.0102

x <- eclipsing_binary$phase
y <- eclipsing_binary$flux
weights <- 1 / eclipsing_binary$std_err^2

cv_tf <- cv_trendfilter(
  x = x,
  y = y,
  weights = weights,
  max_iter = 1e4,
  obj_tol = 1e-6
)

boot_tf <- bootstrap_trendfilter(
  obj = cv_tf,
  algorithm = "nonparametric",
  edf = cv_tf$edf_min["MAE"]
)
bands <- vbands(boot_tf)


# Example 2: The "Lyman-alpha forest" in the spectrum of a distant quasar

data("quasar_spectrum")
head(quasar_spectrum)

#> # A tibble: 6 × 3
#>   log10_wavelength   flux weights
#>              <dbl>  <dbl>   <dbl>
#> 1             3.55  0.424  0.0417
#> 2             3.55 -2.11   0.125
#> 3             3.55 -3.78   0.128
#> 4             3.55  3.60   0.121
#> 5             3.55 -4.89   0.134
#> 6             3.55  1.76   0.132 

x <- quasar_spectrum$log10_wavelength
y <- quasar_spectrum$flux
weights <- quasar_spectrum$weights

sure_tf <- sure_trendfilter(x, y, weights)
boot_tf <- bootstrap_trendfilter(
  obj = sure_tf,
  algorithm = "parametric",
  edf = sure_tf$edf_min
)
bands <- vbands(boot_tf)