Computes a Log-Concave Probability Density Estimate via an Active Set Algorithm
Given a vector of observations x_n = (x_1, …, x_n)
with not necessarily equal entries,
activeSetLogCon first computes vectors x_m = (x_1, …, x_m)
and w = (w_1, …, w_m) where w_i is the weight of each x_i s.t.
∑_{i=1}^m w_i = 1.
Then, activeSetLogCon computes a concave, piecewise
linear function \widehat φ_m on [x_1, x_m] with knots only in {x_1, …, x_m} such that
L(φ) = ∑_{i=1}^m w_i φ(x_i) - int_{-∞}^∞ exp(φ(t)) dt
is maximal. To accomplish this, an active set algorithm is used.
activeSetLogCon(x, xgrid = NULL, print = FALSE, w = NA)
x |
Vector of independent and identically distributed numbers, not necessarily unique. |
xgrid |
Governs the generation of weights for observations. See |
print |
|
w |
Optional vector of weights. If weights are provided, i.e. if |
xn |
Vector with initial observations x_1, …, x_n. |
x |
Vector of observations x_1, …, x_m that was used to estimate the density. |
w |
The vector of weights that had been used. Depends on the chosen setting for |
phi |
Vector with entries \widehat φ_m(x_i). |
IsKnot |
Vector with entries IsKnot_i = 1\{\widehat φ_m has a kink at x_i\}. |
L |
The value L(φ_m) of the log-likelihood-function L at the maximum \widehat φ_m. |
Fhat |
A vector (\widehat F_{m,i})_{i=1}^m of the same size as x with entries \widehat F_{m,i} = \int_{x_1}^{x_i} \exp(\widehat φ_m(t)) dt. |
H |
Vector (H_1, …, H_m)' where H_i is the derivative of t \to L(φ + tΔ_i) at zero and Δ_i(x) = \min(x - x_i, 0). |
n |
Number of initial observations. |
m |
Number of unique observations. |
knots |
Observations that correspond to the knots. |
mode |
Mode of the estimated density \hat f_m. |
sig |
The standard deviation of the initial sample x_1, …, x_n. |
Kaspar Rufibach, kaspar.rufibach@gmail.com,
http://www.kasparrufibach.ch
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
activeSetLogCon can be used to estimate a log-concave density. However, to generate an object of
class dlc that allows application of summary and plot we recommend to use logConDens.
The following functions are used by activeSetLogCon:
Log concave density estimation via an iterative convex minorant algorithm can be performed using
icmaLogCon.
## estimate gamma density set.seed(1977) n <- 200 x <- rgamma(n, 2, 1) res <- activeSetLogCon(x, w = rep(1 / n, n), print = FALSE) ## plot resulting functions par(mfrow = c(2, 2), mar = c(3, 2, 1, 2)) plot(res$x, exp(res$phi), type = 'l'); rug(x) plot(res$x, res$phi, type = 'l'); rug(x) plot(res$x, res$Fhat, type = 'l'); rug(x) plot(res$x, res$H, type = 'l'); rug(x) ## compute and plot function values at an arbitrary point x0 <- (res$x[100] + res$x[101]) / 2 Fx0 <- evaluateLogConDens(x0, res, which = 3)[, "CDF"] plot(res$x, res$Fhat, type = 'l'); rug(res$x) abline(v = x0, lty = 3); abline(h = Fx0, lty = 3) ## compute and plot 0.9-quantile of Fhat q <- quantilesLogConDens(0.9, res)[2] plot(res$x, res$Fhat, type = 'l'); rug(res$x) abline(h = 0.9, lty = 3); abline(v = q, lty = 3)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.