medde {REBayes}

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
)

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

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).

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

Author(s)

Roger Koenker and Ivan Mizera

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.

See Also

This function is based on an earlier function of the same name in the deprecated package MeddeR that was based on an R-Matlab interface. A plotting method is available, or medde estimates can be added to plots with the usual lines(meddefit, ... invocation. For log concave estimates there is also a quantile function qmedde and a random number generation function rmedde, eventually there should be corresponding functionality for other alphas.

Examples

## Not run: 
#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")

## End(Not run)

## Not run: 
#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")

## End(Not run)