activeSetLogCon.mode: Computes the Modally-Constrained Log-Concave Probability Density Maximum Likelihood Estimate via an Active Set Algorithm.

Description

This is an adapted version of activeSetLogCon from the logcondens package for computing the MLE of a log-concave density with known location of mode.

Given a vector of observations ${\bold{x}_n} = (x_1, \ldots, x_n)$ with potentially distinct or nondistinct entries, activeSetLogCon.mode 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$. The vector $\bold{x}_m$ contains the fixed location of the mode, mode. Then, activeSetLogCon.mode computes a concave, piecewise linear function $\widehat \phi_m^0$ on $[x_1, x_m]$ with $p$ knots only in ${x_1, \ldots, x_m}$ and with mode value, mode, such that $$L(\phi) = \sum_{i=1}^m w_i \phi(x_i) - \int_{-\infty}^\infty \exp(\phi(t)) dt$$ is maximal. To accomplish this, an active set algorithm is used.

Usage

activeSetLogCon.mode(x, xgrid = NULL, mode=x[1], print = FALSE, w
= NA, prec=10^-10)

Arguments

Vector of independent and identically distributed numbers, not necessarily unique.

xgrid

Governs the generation of weights for observations. See preProcess for details.

mode

This is the constrained value for the location of the mode.

print = TRUE outputs the log-likelihood in every loop, print = FALSE does not. Make sure to tell R to output (press CTRL+W).

Optional vector of weights. If weights are provided, i.e., if w != NA, then xgrid is ignored.

prec

Governs precision of various subfunctions, e.g., the Newton-Raphson procedure.

Value

xnVector with initial observations $x_1, \ldots, x_n$.
xVector of observations $x_1, \ldots, x_m$ that was used to estimate the density, i.e., points that include all possible knots of the estimate. Note that $\bold{x}$ always includes the mode value mode, since that point is a possible knot! Note also that this $x$ is not identical to the $x$ passed in (xn is identical). This vector is referred to as 'z' in Doss (2013).
wThe vector of weights that had been used. Depends on the chosen setting for xgrid. Of the same length as x. The weight corresponding to the mode will be 0 if the mode is not a data point, and otherwise will be nonzero.
LThe value $L(\widehat {\bold{\phi}}_m^0)$ of the log-likelihood-function $L$ at the maximum $\widehat {\bold{\phi}}_m^0$.
MINumeric vector of length 2 giving the endpoints of the modal interval.
IsKnotVector with entries IsKnot$_i = 1{\widehat{\phi}_m^0$ has a kink at $x_i$}.
IsMICAnalogous to IsKnot; stands for "Is Modally Inactive Constraint," i.e., denotes whether the modal constraints are active or inactive. It is a numeric vector of length 2, corresponding to whether the mode is a left-knot or a right-knot. Just as with IsKnot, a 1 denotes an inactive constraint and a 0 denotes an active one. Thus a 0 indicates that the constraint that the estimate be equal in value at the mode and the nearest knot to the left or to the right, respectively, is active. Note also that if max(IsMIC)==1 then the corresponding index in IsKnot is a 1 (i.e., IsKnot[dlcMode$idx] ==1 ).
constrknots[constr] is equal to MI; that is, constr is a numeric (integral) vector of length two with values in $1, \ldots, p$ indicating which of the $p$ knots are the left and right of the modal interval.
knotsknots equals x[IsKnot>0], gives the values of the points that are knots.
phiVector with entries $\widehat \phi_m(x_i)$, $i=1,\ldots,m$. Named "phi" not "phihat" for backwards compatibility.
fhatVector with entries $\widehat{f}_m^0(x_i) = e^{\widehat{\phi}_m^0(x_i)}$, $i=1,\ldots, m$.
FhatA vector $(\widehat F_{m,i}^0)_{i=1}^m$ of the same size as ${\bold{x}}$ with entries
$$\widehat F_{m,i}^0 = \int_{x_1}^{x_i} \exp(\widehat \phi_m^0(t)) dt.$$
HNumeric vector $(H_1, \ldots, H_{m})'$ where $H_i$ is the derivative of $$t \to L(\phi + t\Delta_i)$$ at zero and $\Delta_i(x) = \min(x - x_i, 0)$ if $x_i$ is less than dlcMode$val or $\Delta_i(x) = \min(x_i - x, 0)$ if $x_i$ is greater than dlcMode$val. If $x_i$ is the mode (i.e., equals dlcMode$val) $H_i$, is set to $0$. The corresponding values for the mode are accessed via H.m.
Note that in the unconstrained problems the derivatives in the directions $\min(x_i - x, 0)$ and $\min(x-x_i,0)$ are equal, but in the constrained problem these derivatives are not equal.
H.mVector $(H.m_1, H.m_2)'$ where $H.m_1$ is the derivative of $$t \to L(\phi + t\Delta_i)$$ at zero and $\Delta_1(x) = \min(x - a, 0)$ and $\Delta_2(x) = \min(a-x, 0)$, where $a$ is the mode.
nNumber of initial observations, i.e., length of xn.
m1Number of unique observations. This count excludes the mode if the mode is not a data point (or if xgrid is not NULL then excludes the mode if it is not in the output of preProcess).
mNumber of points used to compute the estimator, i.e., unique observations as well as the mode, i.e., length of x. So is either $m_1 + 1$ or $m_1$ depending on whether dlcMode$isx is FALSE or TRUE, respectively.
dlcModeA list, of class "dlc.mode", with components $val, $idx, and $isx. dlcMode$val gives the constrained mode value, dlcMode$idx gives the corresponding index in x, and dlcMode$isx is TRUE or FALSE depending on whether the value is or is not equal to an element of the vector preProcess(x, xgrid)$x (where x is the argument passed in, not the value returned).
Note, when the mode is not an x value, w[dlcMode$idx] == 0. This can often be used in place of an explicit check via $isx as to whether the mode is or is not an x value.
sigThe standard deviation of the initial sample $x_1, \ldots, x_n$.
phi.fAll outputs named "name.f" are functions corresponding to name. So, phi.f(x) equals $\widehat{\phi}^0_m(x)$.
fhat.fIs a function such that fhat.f(x) equals $\widehat{f}_m^0(x)$.
Fhat.fIs a function such that Fhat.f(x) equals $\widehat{F}_m^0(x)$.
EL.fEL.f(l,u) $= \int_{l}^{u} \widehat{F}_m^0(t) dt$
Note that this is not analogous to H or H.m, which are derivatives of the log likelihood and so have subtracted an integral of the empirical cdf.
ER.fER.f(l,u) $= \int_{l}^{u} (1-\widehat{F}_m^0(t)) dt$
E.fEquals EL.f. Included so as to be compatible (i.e., follow inheritance principles) with activeSetLogCon, which returns an E.f variable.
phiPLNumeric vector of length $m$ with values $(\widehat{\phi}_m^0)'(x_i-)$
phiPRNumeric vector of length $m$ with values $(\widehat{\phi}_m^0)'(x_i+)$
phiPL.fIs a function such that phiPL.f(x) equals $(\widehat{\phi}_m^0)'(x-)$.
phiPR.fIs a function such that phiPR.f(x) equals $(\widehat{\phi}_m^0)'(x+)$.

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

Doss, C. R. (2013). Shape-Constrained Inference for Concave-Transformed Densities and their Modes. PhD thesis, Department of Statistics, University of Washington, in preparation.

Doss, C. R. and Wellner, J. A. (2013). Inference for the mode of a log-concave density. Technical Report, University of Washington, in preparation.

Examples

Run this code

## estimate gamma density

set.seed(1977)
n <- 200
x <- rgamma(n, 2, 1)
TRUEMODE <- 1; ## (2-1)*1
res <- activeSetLogCon.mode(x, mode=TRUEMODE, w = rep(1 / n, n), print = FALSE)

## plot resulting functions
par(mfrow = c(2, 2), mar = c(3, 2, 1, 2))
plot(res$x, res$fhat, type = 'l'); rug(res$xn)
plot(res$x, res$phi, type = 'l'); rug(res$xn)
plot(res$x, res$Fhat, type = 'l'); rug(res$xn)
plot(res$x, res$H, type = 'l'); rug(res$xn)

## Or can use the ".f" functions 
xpts <- seq(from=0, to=9, by=.01)
par(mfrow = c(2, 2), mar = c(3, 2, 1, 2))
plot(xpts, res$fhat.f(xpts), type = 'l'); rug(res$xn) 
plot(xpts, res$phi.f(xpts), type = 'l'); rug(res$xn) 
## these are not analogous to res$H.
plot(xpts, res$EL.f(upper=xpts), type = 'l'); rug(res$xn) 
plot(xpts, res$ER.f(lower=xpts), type = 'l'); rug(res$xn)


## 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
alpha <- .1
q <- quantilesLogConDens(1-alpha, res)[2]
plot(res$x, res$Fhat, type = 'l'); rug(res$x)
abline(h = 1-alpha, lty = 3); abline(v = q, lty = 3)

Run the code above in your browser using DataLab