activeSetLogCon
first computes vectors ${\bold{x}_m} = (x_1, \ldots, x_m)$ and
${\bold{w}} = (w_1, \ldots, w_m)$ where $w_i$ is the weight of each $x_i$ s.t.
$\sum_{i=1}^m w_i = 1$.
Then, activeSetLogCon
computes a concave, piecewise
linear function $\widehat \phi_m$ on $[x_1, x_m]$ with knots only in ${x_1, \ldots, x_m}$ such that
$$L(\phi) = \sum_{i=1}^m w_i \phi(x_i) - \int_{-\infty}^\infty \exp(\phi(t)) dt$$
is maximal. In order to be able to apply the pool - adjacent - violaters algorithm, computations are performed
in the parametrization
$${\bold{\eta}}({\bold{\phi}}) = \Bigl(\phi_1, \Bigl(\eta_1 + \sum_{j=2}^i (x_i-x_{i-1})\eta_i\Bigr)_{i=2}^m \Bigr).$$
To find the maximum of $L$, a variant of the iterative convex minorant using the pool - adjacent - violaters
algorithm is used.icmaLogCon(x, xgrid = NULL, eps = 10^-8, T1 = 2000,
robustif = TRUE, print = FALSE)
preProcess
for details.robustif = TRUE
performs the robustification and Hermite interpolation procedure detailed in
Rufibach (2006, 2007), robustif = FALSE
does not. In the latter case, convergence of the algorithm
isprint = TRUE
outputs log-likelihood in every loop, print = FALSE
does not. Make sure to
tell R
to output (press CTRL+W).xgrid
.icmaLogCon
can be used to estimate a log-concave density. However, to generate an object of
class dlc
that allows application of summary
and plot
one has to
use logConDens
.
The following functions are used by icmaLogCon
:
phieta
, etaphi
, Lhat_eta
, quadDeriv
,
robust
, isoMean
.set.seed(1977)
x <- rgamma(200, 2, 1)
res <- icmaLogCon(x, T1 = 2000, robustif = TRUE, print = TRUE)
## plot resulting functions
par(mfrow = c(2, 1), mar = c(3, 2, 1, 2))
plot(x, exp(res$phi), type = 'l'); rug(x)
plot(x, res$phi, type = 'l'); rug(x)
Run the code above in your browser using DataLab