logcondens (version 2.1.5)

logConCI: Compute pointwise confidence interval for a density assuming log-concavity

Description

Compute approximate confidence interval for the true log-concave density, on a grid of points. Two main approaches are implemented: In the first, the confidence interval at a fixed point is based on the pointwise asymptotic theory for the log-concave maximum likelihood estimator (MLE) developed in Balabdaoui, Rufibach, and Wellner (2009). In the second, the confidence interval is estimated via the boostrap.

Usage

logConCI(res, xx0, conf.level = c(0.8, 0.9, 0.95, 0.99)[3], 
    type = c("DR", "ks", "nrd", "ECDFboot", "NPMLboot")[2], 
    htype = c("hscv", "hlscv", "hpi", "hns")[4], BB = 500)

Arguments

res

An object of class dlc, usually a result of a call to logConDens.

xx0

Vector of grid points at which to calculate the confidence interval.

conf.level

Confidence level for the confidence interval(s). The default is 95%.

type

Vector of strings indicating type of confidence interval to compute. When type = ks is chosen, then htype should also be specified. The default is type = ks.

htype

Vector of strings indicating bandwidth selection method if type = ks. The default is htype = hns.

BB

number of iterations in the bootstrap if type = NPMLboot or type = ECDFboot. The default is BB = 500.

Value

The function returns a list containing the following elements:

fhat

MLE evaluated at grid points.

up_DR

Upper confidence interval limit when type = DR.

lo_DR

Lower confidence interval limit when type = DR.

up_ks_hscv

Upper confidence interval limit when type = ks and htype = hscv.

lo_ks_hscv

Lower confidence interval limit when type = ks and htype = hscv.

up_ks_hlscv

Upper confidence interval limit when type = ks and htype = hlscv.

lo_ks_hlscv

Lower confidence interval limit when type = ks and htype = hlscv.

up_ks_hpi

Upper confidence interval limit when type = ks and htype = hpi.

lo_ks_hpi

Lower confidence interval limit when type = ks and htype = hpi.

up_ks_hns

Upper confidence interval limit when type = ks and htype = hns.

lo_ks_hns

Lower confidence interval limit when type = ks and htype = hns.

up_nrd

Upper confidence interval limit when type = nrd.

lo_nrd

Lower confidence interval limit when type = nrd.

up_npml

Upper confidence interval limit when type = NPMLboot.

lo_npml

Lower confidence interval limit when boot = NPMLboot.

up_ecdf

Upper confidence interval limit when boot = ECDFboot.

lo_ecdf

Lower confidence interval limit when boot = ECDFboot.

Details

In Balabdaoui et al. (2009) it is shown that (if the true density is strictly log-concave) the limiting distribution of the MLE of a log-concave density \(\widehat f_n\) at a point \(x\) is

$$n^{2/5}(\widehat f_n(x)-f(x)) \to c_2(x) \bar{C}(0).$$

The nuisance parameter \(c_2(x)\) depends on the true density \(f\) and the second derivative of its logarithm. The limiting process \(\bar{C}(0)\) is found as the second derivative at zero of a particular operator (called the "envelope") of an integrated Brownian motion plus \(t^4\).

Three of the confidence intervals are based on inverting the above limit using estimated quantiles of \(\bar{C}(0)\), and estimating the nuisance parameter \(c_2(x)\). The options for the function logConCI provide different ways to estimate this nuisance parameter. If type = "DR", \(c_2(x)\) is estimated using derivatives of the smoothed MLE as calculated by the function logConDens (this method does not perform well in simulations and is therefore not recommended). If type="ks", \(c_2(x)\) is estimated using kernel density estimates of the true density and its first and second derivatives. This is done using the R package ks, and, with this option, a bandwidth selection method htype must also be chosen. The choices in htype correspond to the various options for bandwidth selection available in ks. If type = "nrd", the second derivative of the logarithm of the true density in \(c_2(x)\) is estimated assuming a normal reference distribution.

Two of the confidence intervals are based on the bootstrap. For type = "ECDFboot" confidence intervals based on re-sampling from the empirical cumulative distribution function are computed. For type = "NPMLboot" confidence intervals based on re-sampling from the nonparametric maximum likelihood estimate of log-concave density are computed. Bootstrap confidence intervals take a few minutes to compute!

The default option is type = "ks" with htype = "hns". Currently available confidence levels are 80%, 90%, 95% and 99%, with a default of 95%.

Azadbakhsh et al. (2014) provides an empirical study of the relative performance of the various approaches available in this function.

References

Azadbakhsh, M., Jankowski, H. and Gao, X. (2014). Computing confidence intervals for log-concave densities. Comput. Statist. Data Anal., to appear.

Baladbaoui, F., Rufibach, K. and Wellner, J. (2009) Limit distribution theory for maximum likelihood estimation of a log-concave density. Ann. Statist., 37(3), 1299--1331.

Tarn Duong (2012). ks: Kernel smoothing. R package version 1.8.10. http://CRAN.R-project.org/package=ks

Examples

# NOT RUN {
# }
# NOT RUN {
## ===================================================
## Confidence intervals at a fixed point for the density
## ===================================================
data(reliability)
x.rel <- sort(reliability)

# calculate 95<!-- % confidence interval(s) and plot the result -->
grid <- seq(min(x.rel), max(x.rel), length.out = 200)
res <- logConDens(x.rel)
ci  <- logConCI(res, grid, type = c("nrd", "ECDFboot"))	

par(las = 1, mar = c(2.5, 3.5, 0.5, 0.5))
hist(x.rel, n = 25, col = gray(0.9), main = "", freq = FALSE, 
    xlab = "", ylab = "", ylim = c(0, 0.0065), border = gray(0.5))
lines(grid, ci$fhat, col = "black", lwd = 2)
lines(grid, ci$lo_nrd, col = "red", lwd = 2, lty = 2)
lines(grid, ci$up_nrd, col = "red", lwd = 2, lty = 2)
lines(grid, ci$lo_ecdf, col = "blue", lwd = 2, lty = 2)
lines(grid, ci$up_ecdf, col = "blue", lwd = 2, lty = 2)
legend("topleft", col = c("black", "blue", "red"), lwd = 2, lty = c(1, 2, 2), legend = 
c("log-concave NPMLE", "CI for type = nrd", "CI for type = ECDFboot"), bty = "n")
# }