Density Estimation with rbbnp

This article is a guide to kernel density estimation with bias_bound_density().

Estimating a density

bias_bound_density() takes a data vector and returns a density estimate with bias-aware pointwise confidence intervals:

bias_bound_density(
  x,                    # Data vector
  eval = NULL,          # Evaluation points (automatic if NULL)
  bw = "cv",            # Positive number, "cv", or "silverman"
  conf_level = 0.95,
  kernel = "schennach",
  control = NULL        # Named list of advanced settings
)

A basic fit, on data drawn from the 2-fold uniform distribution:

# Generate data from the two-fold-uniform distribution
set.seed(42)
x_data <- runif(1000) + runif(1000)

# Estimate density
fit <- bias_bound_density(x_data, bw = 0.08, kernel = "schennach")

# Display results
fit
#> Bias-Bounded Density Estimation
#> 
#> Call:
#> bias_bound_density(x = x_data, bw = 0.08, kernel = "schennach")
#> 
#> Sample size: n = 1000
#> Bandwidth:   h = 0.0800 (user-specified) 
#> Kernel:      Schennach2004 
#> 
#> Bias bound parameters:
#>   A = 12.5357, r = 3.0000
#>   bias bound b1x = 0.0139
#> 
#> Evaluation points: 100 (range: [-0.1746, 2.1780])
#> Confidence level: 95%
#> 
#> Use summary() for detailed statistics
#> Use plot() to visualize results

The result is an S3 object of class bbnp_density. Use coef() for the key parameters and summary() for a fuller report:

# Key parameters
coef(fit)
#>        A        r        h 
#> 12.53574  3.00000  0.08000

# Detailed summary
summary(fit)
#> Summary: Bias-Bounded Density Estimation
#> ============================================================
#> 
#> Call:
#> bias_bound_density(x = x_data, bw = 0.08, kernel = "schennach")
#> 
#> Sample Information:
#>   Sample size (n):  1000
#>   Bandwidth (h):    0.0800
#>   Kernel function:  Schennach2004
#> 
#> Bias Bound Parameters:
#>   A (amplitude):    12.5357
#>   r (decay rate):   3.0000
#>   b1x (bias bound): 0.0139
#>   Xi interval:      [2.4074, 4.5544]
#> 
#> Range Estimation:
#>   Density estimates:
#>    min Q1.25% median   mean Q3.75%    max 
#> 0.0000 0.1111 0.4273 0.4217 0.7125 0.8798 
#> 
#>   Standard errors:
#>    min   mean    max 
#> 0.0000 0.0341 0.0563

Visualizing the fit

The default plot shows the estimate together with the bias band and the confidence interval:

plot(fit)

Element	Description
Estimate	Kernel density estimate \(\hat{f}(x)\)
Bias bound	Bias range: \([\hat{f} - \bar{b}, \hat{f} + \bar{b}]\)
95% CI	\([\hat{f} - \bar{b} - z\hat{\sigma}, \hat{f} + \bar{b} + z\hat{\sigma}]\)

The type = "ft" plot shows how the smoothness parameters are detected from the empirical Fourier transform:

plot(fit, type = "ft")

The legend labels each line:

Empirical |phi|: the empirical Fourier transform \(|\hat{\phi}(\xi)|\)
Fitted envelope: the estimated envelope \(A|\xi|^{-r}\), drawn over the selected window
Shaded band: the selected frequency range \([\underline{\xi}, \bar{\xi}_n]\)

The slope \(r\) measures smoothness: a larger \(r\) means a smoother function.

Choosing the bandwidth

Set bw to a positive number for a fixed bandwidth, or to "cv" or "silverman" for automatic selection. Cross-validation is the default and minimizes the leave-one-out least-squares density objective. The fast path uses linear binning and a zero-padded FFT; it is checked against the direct objective in the package tests:

fit_cv <- bias_bound_density(x_data, bw = "cv", kernel = "normal")
cat("CV bandwidth:", coef(fit_cv)["h"], "\n")
#> CV bandwidth: 0.1280561

plot(fit_cv) + ggtitle("Cross-Validation Bandwidth")

Silverman’s rule of thumb is a faster alternative. The normal-kernel choice is the stats::bw.nrd0() rule; the package multiplies it by 2.34 for Epanechnikov and uses a heuristic multiplier of 1.2 for the two infinite-order kernels:

fit_silv <- bias_bound_density(
  x_data,
  bw = "silverman",
  kernel = "normal"
)
cat("Silverman bandwidth:", coef(fit_silv)["h"], "\n")
#> Silverman bandwidth: 0.09390781

plot(fit_silv) + ggtitle("Silverman's Rule Bandwidth")

Kernel functions

Four kernels are available. The formal coverage theorem assumes an infinite-order kernel whose Fourier transform equals one near the origin, so "schennach" is the recommended theorem-compatible default. The finite-order kernels are useful comparisons, but a fit with them should not be described as covered by the same theorem:

Kernel	Order	Theory status
schennach	Infinite	Recommended theorem-compatible default
sinc	Infinite	Infinite-order comparison
normal	2	Finite-order comparison
epanechnikov	2	Finite-order comparison

The same data fitted with each kernel:

fit_sch <- bias_bound_density(x_data, kernel = "schennach")
fit_sinc <- bias_bound_density(x_data, kernel = "sinc")
fit_norm <- bias_bound_density(x_data, kernel = "normal")
fit_epan <- bias_bound_density(x_data, kernel = "epanechnikov")

grid.arrange(
  plot(fit_sch) + ggtitle("Schennach2004 (infinite order)"),
  plot(fit_sinc) + ggtitle("Sinc (infinite order)"),
  plot(fit_norm) + ggtitle("Normal (2nd order)"),
  plot(fit_epan) + ggtitle("Epanechnikov (2nd order)"),
  ncol = 2
)

The default "schennach" kernel is the package’s infinite-order kernel.

Bias and the bandwidth tradeoff

Changing the bandwidth moves the two pieces of the interval in opposite directions. A larger bandwidth lowers the variance (the sampling band narrows) but raises the bias (the bias band widens); a smaller bandwidth does the reverse. The plots below show the same data at the optimal bandwidth and at half of it:

fit_opt <- bias_bound_density(x_data, bw = "cv", kernel = "sinc")
h_opt <- unname(coef(fit_opt)["h"])

# Bias-bound interval at the optimal bandwidth
result_opt <- fit_opt

# ... and at a smaller (undersmoothed) bandwidth
result_under <- bias_bound_density(x_data, bw = h_opt * 0.5, kernel = "sinc")

grid.arrange(
  plot(result_opt) + ggtitle(paste0("Optimal bandwidth (h = ", round(h_opt, 3), ")")),
  plot(result_under) + ggtitle(paste0("Undersmoothed (h = ", round(h_opt/2, 3), ")")),
  ncol = 1
)

At the optimal bandwidth the bias band is the larger component, so the total interval here is wider than at the undersmoothed bandwidth. That is expected. Under the theorem’s smoothness, window, and kernel conditions, the bias allowance supports inference without shrinking the MSE-oriented bandwidth, and a single sample can still produce a wide band. The Theoretical Background article explains why a naive interval can under-cover when smoothing bias is non-negligible and why undersmoothing, the classical alternative, pays a variance cost.

More options

Advanced settings live in control. For example, set the Fourier frequency window manually with frequency_range:

# Default (automatic)
fit_auto <- bias_bound_density(x_data, bw = 0.08)

# Custom range
fit_custom <- bias_bound_density(
  x_data,
  bw = 0.08,
  control = list(frequency_range = c(2, 8))
)

grid.arrange(
  plot(fit_auto, type = "ft") + ggtitle("Automatic Frequency Range"),
  plot(fit_custom, type = "ft") + ggtitle("Custom Range [2, 8]"),
  ncol = 1
)

The estimate, evaluation grid, and confidence intervals are all available on the fitted object:

# Confidence intervals as matrix
ci <- confint(fit)
head(ci, 10)
#>             lower      upper
#>  [1,] 0.000000000 0.01391246
#>  [2,] 0.000000000 0.01391246
#>  [3,] 0.000000000 0.01391246
#>  [4,] 0.000000000 0.01391246
#>  [5,] 0.000000000 0.01391246
#>  [6,] 0.000000000 0.02323744
#>  [7,] 0.000000000 0.03945872
#>  [8,] 0.000000000 0.05674960
#>  [9,] 0.001988962 0.07633667
#> [10,] 0.014643115 0.09834454

# Evaluation points
x_points <- fit$x
head(x_points)
#> [1] -0.1746134 -0.1508492 -0.1270850 -0.1033208 -0.0795566 -0.0557924

# Density estimates
f_hat <- fit$density
head(f_hat)
#> [1] 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.002945549