logConDens: Compute log-concave density estimator and related quantities In logcondens: Estimate a Log-Concave Probability Density from Iid Observations

Description

Compute the log-concave and smoothed log-concave density estimator.

Usage

 1 2 logConDens(x, xgrid = NULL, smoothed = TRUE, print = FALSE, gam = NULL, xs = NULL) 

Arguments

 x Vector of independent and identically distributed numbers, not necessarily unique. xgrid Governs the generation of weights for observations. See preProcess for details. smoothed If TRUE, the smoothed version of the log-concave density estimator is also computed. print print = TRUE outputs the log-likelihood in every loop, print = FALSE does not. Make sure to tell R to output (press CTRL+W). gam Only necessary if smoothed = TRUE. The standard deviation of the normal kernel. If equal to NULL, gam is chosen such that the variances of the original sample x_1, …, x_n and \hat f_n^* coincide. xs Only necessary if smoothed = TRUE. Either provide a vector of support points where the smoothed estimator should be computed at, or leave as NULL. Then, a sufficiently width equidistant grid of points will be used.

Details

See activeSetLogCon for details on the computations.

Value

logConDens returns an object of class "dlc", a list containing the following components: xn, x, w, phi, IsKnot, L, Fhat, H, n, m, knots, mode, and sig as generated by activeSetLogCon. If smoothed = TRUE, then the returned object additionally contains f.smoothed, F.smoothed, gam, and xs as generated by evaluateLogConDens. Finally, the entry smoothed of type "logical" returnes the value of smoothed.

The methods summary.dlc and plot.dlc are used to obtain a summary and generate plots of the estimated density.

Author(s)

Kaspar Rufibach, kaspar.rufibach@gmail.com,
http://www.kasparrufibach.ch

References

Duembgen, L, Huesler, A. and Rufibach, K. (2010). Active set and EM algorithms for log-concave densities based on complete and censored data. Technical report 61, IMSV, Univ. of Bern, available at http://arxiv.org/abs/0707.4643.

Duembgen, L. and Rufibach, K. (2009). Maximum likelihood estimation of a log–concave density and its distribution function: basic properties and uniform consistency. Bernoulli, 15(1), 40–68.

Duembgen, L. and Rufibach, K. (2011). logcondens: Computations Related to Univariate Log-Concave Density Estimation. Journal of Statistical Software, 39(6), 1–28. http://www.jstatsoft.org/v39/i06

Examples

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 ## =================================================== ## Illustrate on simulated data ## =================================================== ## Set parameters n <- 50 x <- rnorm(n) res <- logConDens(x, smoothed = TRUE, print = FALSE, gam = NULL, xs = NULL) summary(res) plot(res, which = "density", legend.pos = "topright") plot(res, which = "log-density") plot(res, which = "CDF") ## Compute slopes and intercepts of the linear functions that ## compose phi slopes <- diff(res$phi) / diff(res$x) intercepts <- -slopes * res$x[-n] + res$phi[-n] ## =================================================== ## Illustrate method on reliability data ## Reproduce Fig. 2 in Duembgen & Rufibach (2009) ## =================================================== ## Set parameters data(reliability) x <- reliability n <- length(x) res <- logConDens(x, smooth = TRUE, print = TRUE) phi <- res$phi f <- exp(phi) ## smoothed log-concave PDF f.smoothed <- res$f.smoothed xs <- res$xs ## compute kernel density sig <- sd(x) h <- sig / sqrt(n) f.kernel <- rep(NA, length(xs)) for (i in 1:length(xs)){ xi <- xs[i] f.kernel[i] <- mean(dnorm(xi, mean = x, sd = h)) } ## compute normal density mu <- mean(x) f.normal <- dnorm(xs, mean = mu, sd = sig) ## =================================================== ## Plot resulting densities, i.e. reproduce Fig. 2 ## in Duembgen and Rufibach (2009) ## =================================================== plot(0, 0, type = 'n', xlim = range(xs), ylim = c(0, 6.5 * 10^-3)) rug(res$x) lines(res$x, f, col = 2) lines(xs, f.normal, col = 3) lines(xs, f.kernel, col = 4) lines(xs, f.smoothed, lwd = 3, col = 5) legend("topleft", c("log-concave", "normal", "kernel", "log-concave smoothed"), lty = 1, col = 2:5, bty = "n") ## =================================================== ## Plot log-densities ## =================================================== plot(0, 0, type = 'n', xlim = range(xs), ylim = c(-20, -5)) legend("bottomright", c("log-concave", "normal", "kernel", "log-concave smoothed"), lty = 1, col = 2:5, bty = "n") rug(res$x) lines(res\$x, phi, col = 2) lines(xs, log(f.normal), col = 3) lines(xs, log(f.kernel), col = 4) lines(xs, log(f.smoothed), lwd = 3, col = 5) ## =================================================== ## Confidence intervals at a fixed point for the density ## see help file for logConCI() ## =================================================== 

Example output

Warning messages:
1: In rgl.init(initValue, onlyNULL) : RGL: unable to open X11 display
2: 'rgl_init' failed, running with rgl.useNULL = TRUE
call: fun(...)

Estimation of a log-concave density from i.i.d. data

Number of initial observations: n = 50

Number of unique observations (or grid points): m = 50

log-likelihood: -2.34

Maximum likelihood estimate:
Mode: x[35] = 0.36
Value of log-density at mode: -0.85
Value of density at mode: 0.43

Smoothed maximum likelihood estimate:
Mode: -0.0047
Value of log-density at mode: -0.9
Value of density at mode: 0.41

Number of knots of the MLE: 5

Knots of the MLE:
x[1, 16, 21, 35, 50] =
-2.54, -0.53, -0.30,  0.36,  2.10

[1] "iter1 = 1 / L = -6.7569 / max(H) = 2.163 / #knots = 3"
[1] "iter1 = 2 / L = -6.7337 / max(H) = 1.3981 / #knots = 4"
[1] "iter1 = 3 / L = -6.7266 / max(H) = 0.4054 / #knots = 5"
[1] "iter1 = 4 / L = -6.7212 / max(H) = 0.2512 / #knots = 6"
[1] "iter1 = 5 / L = -6.719 / max(H) = 0.1588 / #knots = 7"
[1] "iter1 = 6 / L = -6.7177 / max(H) = 0.0748 / #knots = 7"
[1] "iter1 = 7 / L = -6.7176 / max(H) = 0.0363 / #knots = 8"
[1] "iter1 = 8 / L = -6.7168 / max(H) = 0.0321 / #knots = 8"
[1] "iter1 = 9 / L = -6.7166 / max(H) = 0.0207 / #knots = 9"
[1] "iter1 = 10 / L = -6.7165 / max(H) = 0.0174 / #knots = 9"
[1] "iter1 = 11 / L = -6.7164 / max(H) = 0.0074 / #knots = 9"
[1] "iter1 = 12 / L = -6.7163 / max(H) = 0.005 / #knots = 9"
[1] "iter1 = 13 / L = -6.7163 / max(H) = 2e-04 / #knots = 10"
[1] "iter1 = 14 / L = -6.7163 / max(H) = 0 / #knots = 10"
[1] "iter1 = 15 / L = -6.7163 / max(H) = 0 / #knots = 10"
[1] "10% of computation of smooth estimates done"
[1] "20% of computation of smooth estimates done"
[1] "30% of computation of smooth estimates done"
[1] "40% of computation of smooth estimates done"
[1] "50% of computation of smooth estimates done"
[1] "60% of computation of smooth estimates done"
[1] "70% of computation of smooth estimates done"
[1] "80% of computation of smooth estimates done"
[1] "90% of computation of smooth estimates done"
[1] "100% of computation of smooth estimates done"


logcondens documentation built on May 2, 2019, 6:11 a.m.