Given a vector of observations \({\bold{x}_n} = (x_1, \ldots, x_n)\)
with not necessarily equal entries,
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)
Vector of independent and identically distributed numbers, not necessarily equal.
Governs the generation of weights for observations. See preProcess
for details.
An arbitrary real number, typically small. Iterations are halted if the directional derivative of \({\bold{\eta}} \to L({\bold{\eta}})\) in the direction of the new candidate is \(\le \varepsilon\).
Maximal number of iterations to perform.
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
is no longer guaranteed.
print = TRUE
outputs log-likelihood in every loop, print = FALSE
does not. Make sure to
tell R
to output (press CTRL+W).
Vector of observations \(x_1, \ldots, x_m\) that was used to estimate the density.
The vector of weights that had been used. Depends on the chosen setting for xgrid
.
Vector with entries \(\widehat f_m(x_i).\)
Vector with initial observations \(x_1, \ldots, x_n\).
The value \(L(\widehat \phi_m)\) of the log-likelihood-function \(L\) at the maximum \(\widehat \phi_m.\)
Number of iterations performed.
The standard deviation of the initial sample \(x_1, \ldots, x_n\).
Rufibach K. (2006) Log-concave Density Estimation and Bump Hunting for i.i.d. Observations. PhD Thesis, University of Bern, Switzerland and Georg-August University of Goettingen, Germany, 2006. Available at http://www.zb.unibe.ch/download/eldiss/06rufibach_k.pdf.
Rufibach, K. (2007) Computing maximum likelihood estimators of a log-concave density function. J. Stat. Comput. Simul. 77, 561--574.
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
:
# NOT RUN {
set.seed(1977)
x <- rgamma(200, 2, 1)
# }
# NOT RUN {
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