I’m proposing to add the Grenander estimator to scipy.stats as a nonparametric maximum likelihood estimator (MLE) of a probability density under a monotonicity constraint.
Whereas scipy.stats.ecdf provides the MLE of a distribution function, there is no corresponding MLE for a density without additional structure. Imposing monotonicity yields a well-defined and classical solution, called the Grenander estimator, which this proposal makes available in SciPy.
The estimator can equivalently be viewed as a data-adaptive histogram with decreasing heights and variable bin widths. However, unlike histograms or kernel density estimators, it requires no bandwidth or tuning parameters.
To my knowledge, there is currently no implementation of the Grenander estimator in a standard Python library. Given its close connections to ecdf and isotonic_regression and its lack of tuning parameters, scipy.stats seems like the natural place to make it available.
A draft implementation and tests are available in this PR: gh-24244. The implementation follows the standard characterization as the left derivative of the least concave majorant of the empirical CDF and is built directly on scipy.optimize.isotonic_regression, connecting naturally to existing SciPy functionality.