medde: Maximum Entropy [De]Regularized Density Estimation

Description

Density estimation based on maximum entropy methods

Usage

medde(
  x,
  v = 300,
  lambda = 0.5,
  alpha = 1,
  Dorder = 1,
  w = NULL,
  mass = 1,
  rtol = 1e-06,
  verb = 0,
  control = NULL
)

Value

An object of class "medde" with components

x: points of evaluation on the domain of the density
y: estimated function values at the evaluation points x
status: exit status from Mosek

Arguments

x: Data: either univariate or bivariate, the latter is highly experimental
v: Undata: either univariate or bivariate, univariate default is an equally spaced grid of 300 values, for bivariate data there is not (yet) a default. Making v extend well beyond the support of x is advisable to avoid weird boundary behavior of the estimated density.
lambda: total variation penalty smoothing parameter, if lambda is in [-1,0], a shape constraint is imposed. see Koenker and Mizera (2010) for further details. When Dorder = 0, the shape constraint imposes that the density is monotonically decreasing, when Dorder = 1 it imposes a concavity constraint.
alpha: Renyi entropy parameter characterizing fidelity criterion by default 1 is log-concave and 0.5 is Hellinger.
Dorder: Order of the derivative operator for the penalty default is Dorder = 1, corresponding to TV norm constraint on the first derivative, or a concavity constraint on some transform of the density. Dorder = 0 imposes a TV penalty on the function itself, or when lambda < 0 a monotonicity constraint.
w: weights associated with x,
mass: normalizing constant for fitted density,
rtol: Convergence tolerance for Mosek algorithm,
verb: Parameter controlling verbosity of solution, 0 for silent, 5 gives rather detailed iteration log.
control: Mosek control list see KWDual documentation

Author

Roger Koenker and Ivan Mizera

Details

See the references for further details. And also Mosek "Manuals". The acronym, according to the urban dictionary has a nice connection to a term used in Bahamian dialect, mostly on the Family Islands like Eleuthera and Cat Island meaning "mess with" "get involved," "get entangled," "fool around," "bother:" "I don't like to medder up with all kinda people" "Don't medder with people (chirren)" "Why you think she medderin up in their business."

This version implements a class of penalized density estimators solving: $$\min_x \phi(x_1) | A_1 x_1 - A_2 x_2 = b, 0 \le x_1, -\lambda \le x_2 \le \lambda $$ where $x$ is a vector with two component subvectors: $x_1$ is a vector of function values of the density $x_2$ is a vector of dual values, $\lambda$ is typically positive, and controls the fluctuation of the Dorder derivative of some transform of the density. When alpha = 1 this transform is simply the logarithm of the density, and Dorder = 1 yields a piecewise exponential estimate; when Dorder = 2 we obtain a variant of Silverman's (1982) estimator that shrinks the fitted density toward the Gaussian, i.e. with total variation of the second derivative of $log f$ equal to zero. See demo(Silverman) for an illustration of this case. If $\lambda$ is in $(-1,0]$ then the $x_2$ TV constraint is replaced by $x_2 \geq 0$, which for $\alpha = 1$, constrains the fitted density to be log-concave; for $\alpha = 0.5$, $-1/\sqrt f$ is constrained to be concave; and for $\alpha \le 0$, $1/f^{\alpha -1}$ is constrained to be concave. In these cases no further regularization of the smoothness of density is required as the concavity constraint acts as regularizer. As explained further in Koenker and Mizera (2010) and Han and Wellner (2016) decreasing $\alpha$ constrains the fitted density to lie in a larger class of quasi-concave densities. See demo(velo) for an illustration of these options, but be aware that more extreme $\alpha$ pose more challenges from an numerical optimization perspective. Fitting for $\alpha < 1$ employs a fidelity criterion closely related to Renyi entropy that is more suitable than likelihood for very peaked, or very heavy tailed target densities. For $\lambda < 0$ fitting for Dorder != 1 proceed at your own risk. A closely related problem is illustrated in the demo Brown which imposes a convexity constraint on $0.5 x^2 + log f(x)$. This ensures that the resulting Bayes rule, aka Tweedie formula, is monotone in $x$, as described further in Koenker and Mizera (2013).

References

Chen, Y. and R.J. Samworth, (2013) "Smoothed log-concave maximum likelihood estimation with applications", Statistica Sinica, 23, 1373--1398.

Han, Qiyang and Jon Wellner (2016) ``Approximation and estimation of s-concave densities via Renyi divergences, Annals of Statistics, 44, 1332-1359.

Koenker, R and I. Mizera, (2007) ``Density Estimation by Total Variation Regularization,'' Advances in Statistical Modeling and Inference: Essays in Honor of Kjell Doksum, V.N. Nair (ed.), 613-634.

Koenker, R and I. Mizera, (2006) ``The alter egos of the regularized maximum likelihood density estimators: deregularized maximum-entropy, Shannon, Renyi, Simpson, Gini, and stretched strings,'' Proceedings of the 7th Prague Symposium on Asymptotic Statistics.

Koenker, R and I. Mizera, (2010) ``Quasi-Concave Density Estimation'' Annals of Statistics, 38, 2998-3027.

Koenker, R and I. Mizera, (2013) ``Convex Optimization, Shape Constraints, Compound Decisions, and Empirical Bayes Rules,'' JASA, 109, 674--685.

Koenker, R and I. Mizera, (2014) ``Convex Optimization in R.'', Journal of Statistical Software, 60, 1-23.

Examples

Run this code


if (FALSE) {
#Maximum Likelihood Estimation of a Log-Concave Density
set.seed(1968)
x <- rgamma(50,10)
m <- medde(x, v = 50, lambda = -.5, verb = 5)
plot(m, type = "l", xlab = "x", ylab = "f(x)")
lines(m$x,dgamma(m$x,10),col = 2)
title("Log-concave Constraint")
}

if (FALSE) {
#Maximum Likelihood Estimation of a Gamma Density with TV constraint
set.seed(1968)
x <- rgamma(50,5)
f <- medde(x, v = 50, lambda = 0.2, verb = 5)
plot(f, type = "l", xlab = "x", ylab = "f(x)")
lines(f$x,dgamma(f$x,5),col = 2)
legend(10,.15,c("ghat","true"),lty = 1, col = 1:2)
title("Total Variation Norm Constraint")
}

Run the code above in your browser using DataLab