medde | R Documentation |
Density estimation based on maximum entropy methods
medde(
x,
v = 300,
lambda = 0.5,
alpha = 1,
Dorder = 1,
w = NULL,
mass = 1,
rtol = 1e-06,
verb = 0,
control = NULL
)
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. |
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 |
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).
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 |
Roger Koenker and Ivan Mizera
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.
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.
## 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)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.