Package {tiltdens}

tiltdens: Tilted and Data-Sharpened Nonparametric Density Estimation

Description

High-order nonparametric density estimators built by perturbing a conventional kernel estimator, either by re-weighting the observations ("tilting") or by moving them ("data sharpening"). The perturbation is chosen so that the estimator inherits the fast convergence rate of an infinite-order kernel estimator, such as the sinc or trapezoidal flat-top estimator, while remaining a proper non-negative density without the oscillatory tails those estimators suffer from. Two criteria are provided: minimising the L2 distance to an infinite-order comparator, following Doosti and Hall (2016) doi:10.1111/rssb.12112, and minimising a cross-validation criterion that needs no comparator and is much faster, following Doosti, Hall and Mateu (2018) doi:10.1016/j.jspi.2017.12.003.

Author(s)

Maintainer: Hassan Doosti hassan.doosti@mq.edu.au [copyright holder]

See Also

Useful links:

• https://github.com/DoostiH/tiltdens

• Report bugs at https://github.com/DoostiH/tiltdens/issues

Bandwidth selectors

Description

Three rules, for three different jobs.

Usage

bw_comparator_cv(
  x,
  comparator = c("sinc", "trapezoid"),
  q_step = 0.05,
  q_max = NULL,
  q_max_factor = 2,
  q_grid = NULL
)

bw_flattop(x, c_thresh = 2, kn = NULL, t_max = NULL, n_grid = 4000)

bw_nrd_robust(x, kernel = "gaussian")

bw_canonical_factor(kernel = "gaussian")

bw_convert(bw, from, to)

Arguments

Details

'bw_comparator_cv()' chooses the bandwidth of an infinite-order comparator estimator by minimising least-squares cross-validation,

CV(h) = \int \check f(x|h)^2 dx - \frac{2}{n}\sum_i \check f_{-i}(x_i|h),

equation (2.7) of Doosti and Hall (2016). It is the default throughout the package, and it is the rule the bandwidths used for the published simulations were computed with.

'bw_flattop()' implements the empirical rule of Politis (2003): take \hat q to be the smallest frequency beyond which the modulus of the empirical characteristic function stays below c\sqrt{\log_{10}(n)/n} over a stretch of length 'kn', and set h = 1/\hat q. It is also used internally to cap the cross-validation search.

'bw_nrd_robust()' is a rule of thumb for the *conventional* estimator: (4/3)^{1/5} n^{-1/5} times a robust scale, rescaled for the kernel in use (see below).

'bw_canonical_factor()' and 'bw_convert()' handle the fact that bandwidths mean different things for different kernels.

Value

A single positive number. 'bw_comparator_cv()' attaches an attribute '"cv"', a data frame of the frequency grid and criterion, which is worth plotting: the criterion is multimodal, and seeing that is more informative than the single number.

Why the cross-validation search is bounded

Least-squares cross-validation for an infinite-order kernel is multimodal in 'h', and its deepest minimum is often a spurious one at a very small bandwidth. Past the frequency at which the empirical characteristic function decays into sampling noise, the criterion is chasing noise, so its global minimum badly undersmooths. The search is therefore capped at 'q_max_factor' times the frequency identified by 'bw_flattop()'.

Canonical bandwidths

A bandwidth is only meaningful relative to the kernel it scales: 'bw = 0.5' smooths far less with the Epanechnikov kernel, supported on [-1,1], than with the Gaussian, which has unit variance. Marron and Nolan (1988) resolve this with the canonical factor

\delta_K = \left(R(K)/\mu_2(K)^2\right)^{1/5},

where R(K) = \int K^2 and \mu_2(K) = \int u^2K(u)\,du. Bandwidths in units of \delta_K are comparable across kernels: the same number produces the same amount of smoothing, and the asymptotically optimal bandwidths coincide.

'bw_canonical_factor()' returns \delta_K. 'bw_convert()' maps a bandwidth from one kernel to the equivalent for another, multiplying by the ratio of their factors. The estimators use this automatically: when you do not supply 'bw', the bandwidth chosen for the comparator is rescaled to the kernel in use, so switching kernels changes the shape of the fit without changing how much it is smoothed.

References

Doosti, H. and Hall, P. (2016). Making a non-parametric density estimator more attractive, and more accurate, by data perturbation. *Journal of the Royal Statistical Society B* **78**, 445-462.

Politis, D. N. (2003). Adaptive bandwidth choice. *Journal of Nonparametric Statistics* **15**, 517-533.

Marron, J. S. and Nolan, D. (1988). Canonical kernels for density estimation. *Statistics and Probability Letters* **7**, 195-199.

Examples

set.seed(1)
x <- c(rnorm(50, -1.5), rnorm(50, 1.5))

bw_comparator_cv(x)
bw_flattop(x)
bw_nrd_robust(x)

## The criterion is multimodal; the search is capped to avoid the
## spurious minima at small bandwidths.
cv <- attr(bw_comparator_cv(x), "cv")
plot(cv$q, cv$cv, type = "l", xlab = "frequency 1/h", ylab = "CV")

## Bandwidths are not comparable across kernels until rescaled.
bw_canonical_factor("gaussian")
bw_canonical_factor("epanechnikov")
bw_convert(0.5, from = "gaussian", to = "epanechnikov")

Infinite-order kernel density estimators

Description

Evaluate the two comparator estimators used by the tilting and data sharpening methods. Both are built from kernels of unlimited order: they converge very fast when the underlying density is smooth, but they can take negative values and they oscillate in the tails. Removing those defects without giving up the convergence rate is the point of [tilt_density()].

Usage

sinc_density(x, bw = NULL, n = 512, from = NULL, to = NULL)

trapezoid_density(x, bw = NULL, n = 512, from = NULL, to = NULL)

Arguments

Details

'sinc_density()' uses L(u) = \sin(u)/(\pi u), whose Fourier transform is flat on [-1, 1] and zero outside.

'trapezoid_density()' uses L(u) = (\cos u - \cos 2u)/(\pi u^2), whose Fourier transform is 1 on [-1, 1] and tapers linearly to 0 at |z| = 2 (Politis 2003; Politis and Romano 1999).

Value

An object of class '"density"', so that 'plot()', 'lines()' and 'print()' work as they do for [stats::density()].

References

Politis, D. N. (2003). Adaptive bandwidth choice. *Journal of Nonparametric Statistics* **15**, 517-533.

Examples

set.seed(1)
x <- c(rnorm(50, -1.5), rnorm(50, 1.5))
plot(sinc_density(x))
abline(h = 0, col = "grey")   # the estimate dips below zero

Integrated squared error against a known density

Description

Approximates \int (\hat f - f)^2 by the trapezoidal rule over the grid of the fit. Useful for simulation studies where the truth is known.

Usage

ise(object, true_density)

Arguments

Value

A single number.

Examples

set.seed(1)
fit <- tilt_density(rnorm(80), m = 3)
ise(fit, dnorm)

Kernels for the perturbed estimator

Description

The estimator \hat f(x) = \sum_i p_i K_h(x - x_i) needs a kernel 'K'. The papers use the standard normal throughout, noting in Section 4.1 of Doosti and Hall (2016) that it "simplified numerical work"; the theory in Section 3.1 is in fact stated for kernels satisfying condition (3.1), which holds when 'K' is a k-fold convolution of a Laplace density. Both families are available here, along with the usual compactly supported kernels.

Usage

tilt_kernels()

tilt_kernel(name = "gaussian")

Arguments

Details

A custom kernel may be supplied as a list with elements 'dens', 'ft' and 'support'. 'dens' and 'ft' must be vectorised, 'ft' must satisfy 'ft(0) == 1', and 'support' is the half-width of the support or 'Inf'. The self-convolution is then computed numerically.

Value

'tilt_kernel()' returns a list with components 'name', 'dens' (the kernel), 'ft' (its Fourier transform), 'self_conv' (the self-convolution at bandwidth 1) and 'support' ('Inf' for the unbounded kernels). 'tilt_kernels()' returns the available names.

Available kernels

'"gaussian"'

The standard normal density. The default, and the one used for every numerical result in both papers.

'"laplace"'

\tfrac12 e^{-|u|}. The k = 1 case of condition (3.1) of the 2016 paper.

'"laplace2"'

\tfrac14 (1+|u|)e^{-|u|}, the twofold convolution of a Laplace density and the k = 2 case of condition (3.1). Given as an example in the paper.

'"laplace3"', '"laplace4"'

Higher-order members of the same family, progressively smoother.

'"epanechnikov"'

\tfrac34 (1-u^2) on [-1,1]. The minimum-variance choice under the usual asymptotics.

'"biweight"'

\tfrac{15}{16}(1-u^2)^2 on [-1,1].

'"triweight"'

\tfrac{35}{32}(1-u^2)^3 on [-1,1].

'"triangular"'

1-|u| on [-1,1].

Why any kernel keeps the problem convex

The quadratic form is A_{ij} = (K\ast K)(x_i - x_j), whose Fourier transform is \phi_K^2 \ge 0. A function with a non-negative Fourier transform has a positive semidefinite Gram matrix, so 'A' is positive semidefinite whatever 'K' is, and weight selection remains a convex quadratic program. Equivalently, p^\top A p = \int \hat f^2 \ge 0 by construction.

A caution about bandwidths

Bandwidths are on each kernel's own scale, so they are not comparable across kernels: 'bw = 0.5' means something different for the Gaussian, which has unit variance, than for the Epanechnikov, which is supported on [-1,1]. If you switch kernels, let the bandwidth be chosen afresh rather than carrying a number over.

Examples

tilt_kernels()

set.seed(1)
x <- c(rnorm(50, -1.5), rnorm(50, 1.5))

## The Laplace-convolution kernel of the paper's condition (3.1).
fit <- tilt_density(x, m = 3, kernel = "laplace2")
fit

Contour plot of a bivariate fit

Description

Contour plot of a bivariate fit

Usage

## S3 method for class 'tiltdens_md'
plot(x, n = 60, points = TRUE, ...)

Arguments

Value

'x', invisibly.

Evaluate a fitted density at new points

Description

Evaluate a fitted density at new points

Usage

## S3 method for class 'tiltdens'
predict(object, newdata = NULL, ...)

## S3 method for class 'tiltdens_md'
predict(object, newdata = NULL, ...)

Arguments

Value

A numeric vector of density values.

Examples

set.seed(1)
fit <- tilt_density(rnorm(60), m = 3)
predict(fit, newdata = c(-1, 0, 1))

Print a tilted or sharpened density fit

Description

Print a tilted or sharpened density fit

Usage

## S3 method for class 'tiltdens'
print(x, digits = getOption("digits") - 2L, ...)

Arguments

Value

'x', invisibly.

Print a multivariate fit

Description

Print a multivariate fit

Usage

## S3 method for class 'tiltdens_md'
print(x, digits = getOption("digits") - 2L, ...)

Arguments

Value

'x', invisibly.

Real-coded genetic algorithm for small box-constrained problems

Description

Minimises 'fn' over the box '[lower, upper]'. This is a tidied, dependency-free version of the algorithm used for the data sharpening results of Doosti and Hall (2016).

Usage

real_ga(
  fn,
  lower,
  upper,
  max_iterations = 130L,
  population = 50L,
  crossover_rate = 0.8,
  mutation_rate = 0.4,
  mutation_step = 0.02,
  blend_gamma = 0.3,
  selection_pressure = 8,
  trace = FALSE
)

Arguments

Value

A list with 'par', 'value' and 'history'.

Data-sharpened density estimation

Description

Finds shifts 'q' for the estimator

\hat f(x \mid h, q) = \frac{1}{nh}\sum_i K\!\left(\frac{x - x_i - q_i}{h}\right)

by minimising its L2 distance to an infinite-order comparator. This is the second form of data perturbation in Doosti and Hall (2016): the observations are moved rather than re-weighted.

Usage

sharpen_density(
  x,
  m = Inf,
  comparator = c("sinc", "trapezoid"),
  bw = NULL,
  comparator_bw = NULL,
  kernel = "gaussian",
  breaks = c("equal", "modal"),
  min_block = NULL,
  max_shift = 1,
  polish = TRUE,
  control = list(),
  n = 512,
  from = NULL,
  to = NULL
)

Arguments

Details

The objective is not convex, so the result depends on the random seed. Set one before calling if you need reproducibility. With 'm = Inf' the search is over 'length(x)' variables and is slow; the paper reports run times of a few minutes per sample of size 100. In the published simulations tilting was usually at least as accurate and far cheaper, so [tilt_density()] and [tilt_density_cv()] are the better default.

Value

An object of class '"tiltdens"'; see [tilt_density()]. The 'shifts' and 'sharpened' components hold 'q' and 'x + q'.

References

Doosti, H. and Hall, P. (2016). Making a non-parametric density estimator more attractive, and more accurate, by data perturbation. *Journal of the Royal Statistical Society B* **78**, 445-462.

See Also

[tilt_density()], [tilt_density_cv()]

Examples

set.seed(1)
x <- c(rnorm(15, -1.5), rnorm(15, 1.5))

## A deliberately tiny search budget, so that this runs in about a second.
## It is far too small for real use; see below for a realistic call.
fit <- sharpen_density(x, m = 3, polish = FALSE,
                       control = list(max_iterations = 10, population = 10))
plot(fit)
round(unique(fit$shifts), 3)


## Realistic settings. Expect this to take a few seconds, and much longer
## with m = Inf, where the search ranges over one shift per observation.
x <- c(rnorm(50, -1.5), rnorm(50, 1.5))
fit <- sharpen_density(x, m = 3)
fit

Summarise a tilted or sharpened density fit

Description

Summarise a tilted or sharpened density fit

Usage

## S3 method for class 'tiltdens'
summary(object, ...)

## S3 method for class 'summary.tiltdens'
print(x, digits = getOption("digits") - 2L, ...)

Arguments

Value

An object of class '"summary.tiltdens"': a list with the method, sample size, bandwidth, kernel, the number of distinct weights and their range, the criterion value, and the minimum of the fitted density. It has a 'print' method.

Examples

set.seed(1)
summary(tilt_density_cv(c(rnorm(30, -1.5), rnorm(30, 1.5))))

Cross-validation criterion for a tilted estimator

Description

Evaluates

CV(h, p) = \int \hat f(x \mid h,p)^2 dx - \frac{2}{n}\sum_i \hat f_{-i}(x_i \mid h,p),

equation (2.1) of Doosti, Hall and Mateu (2018), where \hat f_{-i} is the estimator with the term in x_i deleted. As in the paper, 'p' is *not* re-standardised after deletion; that step turns out to be unnecessary and skipping it saves computation.

Usage

tilt_cv(x, bw, weights = NULL, kernel = "gaussian")

Arguments

Details

The first term is p^\top A p and the second is linear in 'p', so this criterion has the same quadratic form as the 2016 distance criterion. The only difference is the linear term: the 2016 method compares against a comparator estimator, the 2018 method against the leave-one-out fit. That is why one solver serves both.

Value

A single number, with attributes '"A"' and '"C"' holding the quadratic and linear parts, which [tilt_density_cv()] reuses.

References

Doosti, H., Hall, P. and Mateu, J. (2018). Nonparametric tilted density function estimation: a cross-validation criterion. *Journal of Statistical Planning and Inference* **197**, 51-68.

Examples

set.seed(1)
x <- rnorm(60)
tilt_cv(x, bw = 0.4)

Tilted density estimation by distance to a comparator

Description

Chooses non-negative weights 'p' summing to one for the estimator

\hat f(x \mid h, p) = \sum_i p_i K_h(x - x_i)

by minimising its L2 distance to an infinite-order comparator estimator. This is the method of Doosti and Hall (2016): it inherits the fast convergence rate of the comparator while staying a proper, non-negative density with no oscillatory tails.

Usage

tilt_density(
  x,
  m = Inf,
  comparator = c("sinc", "trapezoid"),
  bw = NULL,
  comparator_bw = NULL,
  bw_grid = NULL,
  kernel = "gaussian",
  breaks = c("optimal", "equal", "modal"),
  min_block = NULL,
  constraint = "none",
  constraint_range = NULL,
  constraint_grid_n = 100L,
  n = 512,
  from = NULL,
  to = NULL
)

Arguments

Details

The problem is a convex quadratic program over the probability simplex, so it has a unique solution.

Value

An object of class '"tiltdens"', which inherits from '"density"', so 'plot()', 'lines()' and 'points()' work directly. Components beyond those of [stats::density()] include 'weights', 'distance2' (the attained squared L2 distance), 'comparator', 'comparator_bw', 'n_groups' and 'group_of'.

Choosing 'm'

'm = Inf' lets every observation take its own weight, which is the most flexible option and the one the paper calls T_n. 'm = 3' restricts the weights to three distinct values (T_3); it is faster, and in the paper's simulations it was often more accurate, because fewer free parameters means less overfitting.

When 'm' is finite the block boundaries matter as much as the number of blocks, and 'breaks' controls them. The default '"optimal"' treats the boundaries as part of the optimisation, which is what Section 4.1 of the 2016 paper specifies. On bimodal data it typically places them at the trough between the modes without being told to, which is the same answer '"modal"' arrives at by construction from a pilot estimate.

References

Doosti, H. and Hall, P. (2016). Making a non-parametric density estimator more attractive, and more accurate, by data perturbation. *Journal of the Royal Statistical Society B* **78**, 445-462.

See Also

[tilt_density_cv()] for the faster cross-validation variant, [sharpen_density()] for perturbing the observations instead of their weights.

Examples

set.seed(1)
x <- c(rnorm(25, -1.5), rnorm(25, 1.5))

fit <- tilt_density(x, m = 3)
fit
plot(fit)
lines(sinc_density(x), col = "red", lty = 2)
abline(h = 0, col = "grey")

## The comparator dips below zero; the tilted estimate cannot.
min(sinc_density(x)$y)
min(fit$y)


## Choosing the block boundaries, rather than fixing them, lowers the
## criterion the method is minimising. Searching them is the slow option:
## it examines every admissible pair, so its cost grows with n^2.
c(equal   = tilt_density(x, m = 3, breaks = "equal")$distance2,
  modal   = tilt_density(x, m = 3, breaks = "modal")$distance2,
  optimal = tilt_density(x, m = 3, breaks = "optimal")$distance2)

Tilted density estimation by cross-validation

Description

Chooses the bandwidth and the tilt weights together by minimising [tilt_cv()], rather than by comparing against an infinite-order estimator. This is method (III) of Doosti, Hall and Mateu (2018). It needs no comparator and is much cheaper than [tilt_density()], which is the point: the 2016 construction is accurate but slow, and this makes it practical.

Usage

tilt_density_cv(
  x,
  max_groups = 3,
  bw_grid = NULL,
  rho = 0.8,
  kernel = "gaussian",
  breaks = c("equal", "modal", "optimal"),
  min_block = NULL,
  constraint = "none",
  constraint_range = NULL,
  constraint_grid_n = 100L,
  n = 512,
  from = NULL,
  to = NULL
)

Arguments

Details

The sequential scheme of the paper is used. Start with uniform weights and choose the bandwidth by ordinary least-squares cross-validation. At step 'r', minimise 'CV(h, p)' jointly over 'h' and over weights taking 'r + 1' distinct values, subject to 'h >= rho * h_previous'; the constraint reflects the optimal bandwidth growing as the bias shrinks. Stop at 'max_groups' blocks.

Value

An object of class '"tiltdens"'; see [tilt_density()]. The 'trace' component records the number of blocks, bandwidth and criterion value at each step, which is useful for seeing whether extra blocks bought anything.

References

Doosti, H., Hall, P. and Mateu, J. (2018). Nonparametric tilted density function estimation: a cross-validation criterion. *Journal of Statistical Planning and Inference* **197**, 51-68.

See Also

[tilt_density()], [tilt_cv()]

Examples

set.seed(1)
x <- c(rnorm(50, -1.5), rnorm(50, 1.5))

fit <- tilt_density_cv(x)
fit
fit$trace
plot(fit)