# Nonparametric Estimation of the Underlying Distribution of Binned Continuous
Data

Ejike R. Ugba & Jan Gertheiss† To whom correspondence should be addressed:
<EMAIL_ADDRESS>of Mathematics and Statistics, Helmut Schmidt
University, Hamburg, Germany

Abstract

Estimating cumulative distribution functions (CDF) and probability density
functions (PDF) is a fundamental practice in applied statistics. However,
challenges often arise when dealing with data arranged in grouped intervals.
This paper discusses a suitable and highly flexible non-parametric density
estimation approach for binned distributions, based on cubic monotonicity-
preserving splines—known as cubic spline interpolation. Results from
simulation studies demonstrate that this approach outperforms many widely used
heuristic methods. Additionally, the application of this method to a dataset
of train delays in Germany and micro census data on distance and travel time
to work yields both meaningful but also some questionable results.

Keywords: Non-parametric, Interpolated CDF, Monotone spline, Cut-points,
Binning, Grouped data

##  1 Introduction

Nonparametric density estimation has gained much attention in diverse fields
of empirical research. Unlike its parametric counterpart with many restrictive
assumptions, this approach is independent of distributional assumptions. The
literature is inundated with many valuable materials on the nonparametric
estimation of density functions. Some benchmark references include Silverman
(1986), Eubank (1988), Loader (1999b), and Wasserman (2006). The cardinal
point of most nonparametric estimations hinges on the estimation of
probability density functions (PDFs) and cumulative distribution functions
(CDFs), considered crucial cornerstones of statistical estimation and
inference. However, when information on a supposedly continuous attribute of
interest is captured in groups of intervals rather than the original
measurement scale, estimating the CDFs and PDFs of the underlying continuous
variable is often very problematic. A common situation in practice is in
surveys where data on income, for instance, is only available in terms of
intervals or where an age bracket is used instead of capturing the actual age
of respondents. Although some mathematical heuristics may be employed to
obtain some descriptive statistics of interest, such, of course, fail to tell
the full story regarding the true distribution being estimated.

In what follows, we assume that the said grouped data $z$ in terms of
intervals arises from an underlying continuous variable $y^{*}$ such that

$z=j,\;\text{if}\ \tau_{j-1}<y^{*}\leq\tau_{j},$ (1)

with $j=1,\dots,r$ and the thresholds $\tau_{j}$ being known. We could use $z$
to estimate the entire distribution of $y^{*}$, given that we, at least, know
the empirical cumulative distribution function (eCDF) of $y^{*}$ at
$\tau_{j}$, $j=1,\dots,r-1$. In the simplest case, we can just use linear
interpolation to estimate the entire cdf. This is exactly what is accomplished
by drawing a relative cumulative frequency graph. However, linear
interpolation may induce some systematic bias because it implicitly assumes
that the underlying data is uniformly distributed within each interval defined
by the cutpoints. This also means that the underlying pdf is piecewise
constant with potential jumps at the cutpoints, a rather unrealistic
assumption in many applications.

Alternatively, one may consider (monotone) splines of higher order, such as
cubic splines, to obtain a smoother and differentiable function. An example of
this approach is found in a recent paper by von Hippel et al. (2017), where a
similar technique was applied when analyzing binned income data of US
counties. Our study will rigorously evaluate the performance of monotone
spline interpolated cumulative distribution functions (CDFs). This evaluation
will be conducted through an extensive simulation study, encompassing an
underlying CDF with multiple distributional shapes and support.

The subsequent sections of this study are organized as follows: A concise
overview of the related literature is presented in Section 2. We delve into
the heuristic approach and the kernel/eCDF estimation method in Sections 3 and
4, respectively. This is followed by an in-depth presentation of the
interpolated CDF in Section 5. In Section 6, we conduct a comprehensive
simulation study to assess the performance of the monotone spline interpolated
CDF. Real data examples are then explored in Section 7. Finally, Section 8
offers a conclusive discussion of our findings.

##  2 Previous Work

A readily available nonparametric approach to density estimation is the so-
called kernel density estimation, the main ideas of which can be traced back
to works such as Rosenblatt (1956), Parzen (1962), Nadaraya (1964), Yamato
(1973), and Azzalini (1981). These early contributions employed kernel
functions to estimate the distribution of univariate, independent, and
identically distributed (iid) random variables. Liu and Yang (2008)
investigated the asymptotic properties of the integrated smooth kernel
estimator for multivariate and weakly dependent data. In contrast, Funke and
Palmes (2016) explored estimating non-negative valued random variables via
convolution power kernels. Despite these advancements, selecting an
appropriate bandwidth for kernel-type estimation remains an issue in many
empirical applications.

An underexplored avenue of nonparametric density estimation involves using
polynomial curves and spline functions to estimate the cumulative distribution
functions (CDFs) and probability density functions (PDFs) of measurement
attributes. Kim (1996), for instance, discusses the use of Bezier curves, a
variant of X-splines, for nonparametric density estimation. Kim et al. (1999)
demonstrated that such an application of Bezier curves competes favorably with
results obtained from kernel-based estimators and shares the same asymptotic
properties. Bezier curves were also employed by Kim et al. (2003) for
smoothing the Kaplan-Meier estimator and by Bae et al. (2005) for smoothing a
bivariate Kaplan-Meier estimator. Kim and Park (2012) discussed selection
techniques for Bezier points, while Bae et al. (2014) discussed the
nonparametric estimation of the distribution function using the Bezier curve.
However, the subjectivity in determining the shape parameter and the number of
adequate Bezier points in various situations remains a concern.

Monotone spline estimators were employed by Xue and Wang (2009) to estimate
the CDFs of random variables in simulation and real data studies. Huang et al.
(2016) applied spline and kernel estimators to estimate the interval-
aggregated data’s PDF. Furthermore, Charles et al. (2010) used cubic monotone
control smoothing splines to estimate finite interval-defined CDFs. Rizzi et
al. (2015) demonstrated how binned data can be ungrouped via a penalized
composite link model, assuming the original distribution meets a specific
smoothness criterion. Jargowsky and Wheeler (2018) proposed a mean-constrained
integration over brackets (MCIB) approach to improve existing binned data
estimation methods. In a recent paper, Ghosh et al. (2018) proposed a method
for estimating the latent density of ordinal data using a combination of
Bernstein Polynomials and information on latent distribution provided by cut-
points. Their estimation approach is based on the Anderson-Darling statistic
rather than the widely used maximum likelihood estimation and yields good
results. However, such an approach is plagued with a loss of efficiency
arising from its inability to constrict the shape of the latent density to
unimodal when such a situation arises. Some other handful attempts have also
been made in the past towards estimating the distribution of binned data using
both parametric and nonparametric approaches, which in many instances fail to
perform adequately given the discontinuous structure of such data. As hinted
in Section 1, the method used by von Hippel et al. (2017) seems promising and
will be explored further in this study.

##  3 Heuristic Estimator

A simple and widely used method of binned data estimation is the midpoint
estimator (ME), which assigns cases to all class midpoints, enabling easy
computation of some sample or population characteristics of interest. Consider
the threshold model in (1), where the distinct bin frequencies,
$\eta_{1},\eta_{2},\dots,\eta_{r}$ (as shown in Table 1) are subdivisions of
the total observation $n$. Assuming $\tau_{j-1}$ and $\tau_{j}$ are finite and
known, a straightforward estimation approach is to assign corresponding cases
to the respective bin midpoints with the midpoint estimator

$\text{ME}_{j}=(\tau_{j-1}+\tau_{j})/2,\quad j=1,\dots,r,$ (2)

(see, e.g., Heitjan, 1989) and then compute some basic descriptive statistics.
The bin consistency property of ME and its other desirable properties (see von
Hippel et al., 2015) makes it widely acceptable for most grouped data
analysis. A different dimension to it is with open-ended grouped data where
the last class limit, $\tau_{r}$, is considered finite but unknown (e.g., in
most grouped income data), a common practice is to assume the particular bin
in question follows a Pareto distribution with shape parameter $\gamma>1$
(Henson, 1967), and with its midpoint given by
$\mu_{r}=\tau_{r-1}\gamma/(\gamma-1)$, the so-called Pareto midpoint estimator
(PME) (see Henson, 1967). von Hippel et al. (2015) consider the harmonic
rather than the arithmetic mean of the Pareto distribution a more robust
choice, given instead the following robust Pareto mean estimator (RPME):
$\mu^{\prime}_{r}=\tau_{r-1}(1+1/\gamma)$, where the estimate of $\gamma$ is
given by

$\hat{\gamma}=\frac{\ln\\{(\eta_{r-1}+\eta_{r})/\eta_{r}\\}}{\ln(\tau_{r-1}/\tau_{r-2})}$
(3)

Table 1: Frequency distribution of a binned continuous variable ($y^{*}$) with $L_{b}$, $U_{b}$ and $\eta$ the lower and upper class limits, and the frequencies, respectively. $L_{b}$ $U_{b}$ |  $\eta$
---|---
$\tau_{0}\quad\ \ \ \leq y^{*}<\quad\tau_{1}$ |  $\eta_{1}$
$\tau_{1}\quad\ \ \ \leq y^{*}<\quad\tau_{2}$ |  $\eta_{2}$
$\vdots$ |  $\vdots$
$\tau_{r-2}\quad\leq y^{*}<\quad\tau_{r-1}$ |  $\eta_{r-1}$
$\tau_{r-1}\quad\leq y^{*}<\quad\tau_{r}$ |  $\eta_{r}$

However, while such an approach does provide a quick remedy to the problem at
hand, it goes without saying that such an imposed parametric assumption on a
single data point could yield unreliable population estimates and misleading
conclusions. Consequently, a few parametric approaches have been proposed. For
instance, the generalized beta (GB) family has been proposed (McDonald and
Ransom, 2008) and applied (see, von Hippel et al., 2015) in several binned
data arising from large-scale survey studies. To improve the utility of the GB
family, a multimodel GB estimator (MGBE) was further developed in von Hippel
et al. (2015). However, besides not yielding correct fit to data in most
empirical applications, the GB family dependence on iterative methods slows
things down, making it difficult to use in large-scale studies (von Hippel et
al., 2017).

In addition to the familiar measures of central tendency and spread, in most
instances, it is usually of interest to obtain different fractiles of binned
data, including quartiles, deciles, percentiles, and n-tiles in general. Thus,
with binned data, such quantities are generally obtained either with linear
CDF interpolation of observed classes (a relative frequency distribution
graph) or the widely used interpolation formula, both yielding the same
result. For instance, the quantiles of a given binned data could be obtained
as follows

$Q(q)=\tau_{L}+\frac{nq-C}{\eta}W,\quad q\in(0,1),$ (4)

where $q$ is the desired quantile (a value between 0 and 1) and $n$ is the
total number of observations. Assuming $\xi$ denotes the class containing $q$,
$\tau_{L}$ is the lower boundary of $\xi$, $C$ is the cumulative frequency of
the class preceding $\xi$, $\eta$ is the frequency of $\xi$, and $W$ is the
width of $\xi$. The following special quantile: $Q_{.25}$, $Q_{.50}$ and
$Q_{.75}$ respectively denote the first quartile $(Q_{1})$, the median
$(Q_{2})$, and the third quartile $(Q_{3})$ of a given distribution, with the
interquartile range given by

$IQR=Q_{3}-Q_{1}.$ (5)

These and similar measures are often employed empirically to study several
population characteristics of interest. Our goal in this study is to examine
an alternative non-parametric CDF and density estimation approach with which
such measures could be optimally estimated and to compare the same with the
status quo.

In the next section, we present a formal definition of the CDF, how it relates
to the PDF and other advanced approaches to their estimations.

##  4 Kernel and eCDF Estimators

Assuming $X_{1},X_{2},\dots X_{n}$ are independent identically distributed
real-valued random variables with a continuous PDF $f(\tau)$, and
$F_{X}(\tau)=P(X\leq\tau)$ denotes the cumulative distribution function (CDF)
of the same, where the right-hand side represents the probability of $X$
taking on a value less than or equal to $\tau$. $F_{X}(\tau)$ represents a
right-continuous monotonically increasing function (Park, 2018) that uniquely
describes $f(\tau)$. For a continuous $X$, the following relationship between
$F_{X}(\tau)$ and $f(\tau)$ applies:

$F_{X}(\tau)=\int_{-\infty}^{\tau}f_{X}(t)\,dt.$ (6)

This relationship implies that it is possible to obtain $F_{X}(\tau)$ from
$f(\tau)$ if known and vice versa. As earlier hinted in Section 2, several
non-parametric approaches to CDF and density estimations already exist in the
literature. A typical and common one is the so-called kernel density estimator
(Rosenblatt, 1956; Parzen, 1962). In an ideal continuous case, where
$X_{1},X_{2},\dots X_{n}$ are supposed to be observed, the kernel estimator of
the CDF is expressed as follows,

$\hat{F_{h}}(\tau)=\frac{1}{nh}\sum_{i=1}^{n}\mathbb{\kappa}\left(\frac{\tau-
X_{i}}{h}\right),$ (7)

where $\kappa(x)=\int_{-\infty}^{x}K(t)\,dt$ is an integral over a kernel
function $K$, with the bandwidth denoted by $h$. While the choice of $K$ is of
secondary importance, that of $h$ plays a crucial role in the shape of the
estimator given in (7). Two widely used selection approaches for $h$ include
cross-validation and plug-in methods (Jones et al., 1996; Loader, 1999a;
Henderson and Parmeter, 2015). While the first method offers a data-driven
solution to the estimation problem, the second approach seeks to minimize the
asymptotic mean integrated squared error while substituting the unknown
density in the optimization with a pilot estimate. These notwithstanding, the
selection of $h$ and extension to non-continuous data still pose serious
concerns in practice. For instance, trying to apply the second approach to the
group midpoints of the grouped data could result in unreliable estimates since
too much probability mass is falsely allocated to the midpoints and too little
to the unobserved $X_{i}$ (Scott and Sheather, 1985; Wang and Wertelecki,
2013). Unfortunately, choosing a large $h$ to circumvent such problems incurs
a further loss of information about the underlying true distribution (Reyes et
al., 2019), thus offering no helpful solution to the problem. Furthermore,
Minoiu and Reddy (2014) analyzed the performance of kernel density methods
applied to grouped data to estimate poverty using Monte Carlo simulations and
household surveys, with findings indicating that such technique gives rise to
biases in poverty estimates.

A different frequently used estimator of $F_{X}(\tau)$ is the empirical
cumulative distribution function (eCDF) denoted by $\hat{F_{n}}(\tau)$. Given
$n$ number of observations, $\hat{F_{n}}(\tau)$ assigns a relative frequency
(i.e., estimated probability) of $1/n$ to an ordered sequence of the observed
data, resulting in a step function that increases by $1/n$ at each datum. This
may be formally expressed as follows:

$\hat{F_{n}}(\tau)=\frac{1}{n}\sum_{j=1}^{n}\mathbb{1}\\{X_{j}\leq\tau\\}\quad\forall\tau\in\rm
I\\!R$ (8)

Figure 1:  The empirical cumulative distribution function (eCDF, unbinned) of
an iid normally distributed random sample $(n=10000)$, and the eCDF after the
observations are binned, considering a small, medium, and large number of bins
$B$.

where $\mathbb{1}\\{\cdot\\}$ is the indicator function,
$\lim_{\tau\to-\infty}F_{X}(\tau)=0$ and $\lim_{\tau\to\infty}F_{X}(\tau)=1$.
However, similar to the previously discussed estimator, estimating an
underlying CDF from eCDF makes absolute sense with ungrouped continuous random
variables, given a sufficiently large number of observations. Grouped data
requires having a large number of bins to make any meaningful approximation of
an underlying CDF. One has that.

$F_{X}(\tau)=\lim_{B\to\infty}F_{X_{B}}(\tau),$ (9)

where $F_{X_{B}}(\tau)$ is the eCDF of a binned data and $B$ the number of
bins. In practice, a relatively small number of bins ($B$ < 10 or a bit more)
is typically used to capture various population characteristics of interest,
such that eCDFs obtained under such settings, irrespective of the total number
of frequencies, are mere step functions with large jumps that are not
differentiable. This is further illustrated in Figure 1, showing the eCDF of
an unbinned random sample from the normal distribution and the binned version
of the same with $B\in\\{7,15,50\\}$. As observed, with a large number of
bins, $F_{X_{B}}(\tau)$ approximates $F_{X}(\tau)$. However, as a large $B$ is
practically unattainable for single or multiple study characteristics, one at
least needs a smooth differentiable function, irrespective of $B$, with which
the underlying distribution could be approximated. Thus, we will next consider
interpolating the CDFs with splines of higher order.

##  5 Monotone Spline Interpolated CDF

A basic and prevalent form of CDF interpolation for binned data is the
relative cumulative frequency plot (RCFP), which maps the cutpoints $\tau_{j}$
of binned data to the realized cumulative frequencies, connecting them with
straight lines (see, e.g., von Hippel et al., 2017). However, RCFP produces a
function with rough edges at the respective connecting nodes and, as such,
does not represent a smooth function. At best, it offers a crude estimation of
an underlying CDF. A robust and more accurate estimate of $F_{X}(\tau)$ from
binned data could be achieved using appropriate monotone interpolating cubic
splines (MICS) or similar functions (von Hippel et al., 2017). While a more
elaborate simulation study follows in the next section, to buttress the point
a bit more, Figure 2 shows a typical estimated CDF of the standard normal
(upper row) and the gamma (lower row) distributions for $r=\\{3,5\mbox{ and
}10\\}$ using MICS and RCFP. While MICS provides a smooth differentiable curve
that passes through the nodes, RCFP connects the points with straight lines
that are somewhat discontinuous at the nodes, even though both shapes seem to
align together given a large $r$. Moreover, the indicated median (dashed grey
lines) shows how both lines coincide at this point given the normal
distribution but differ for the gamma distribution.

Figure 2:  Estimated CDF of the standard normal (upper row) and the gamma
(lower row) distributions for $r=\\{3,5\mbox{ and }10\\}$ using monotone
interpolating cubic splines (MICS) and relative cumulative frequency plots
(RCFP). Dashed grey horizontal and vertical lines indicate the second quartile
on MICS compared to RCFP.

A whole lot of literature exists on interpolation approaches, with some
benchmark references being papers by Fritsch and Butland (1980), Fritsch and
Carlson (1980), Lehmann et al. (1999) and Hagan and West (2006).

Consider a mesh of data points $\\{\tau_{1},\tau_{2},\dots,\tau_{r}\\}$ and
corresponding values $\\{F_{1},F_{2},\dots,F_{r}\\}$ for a generic but unknown
function $F(\tau)$. Given the following definitions of the mesh sizes and
gradients:

$h_{j}=\tau_{j+1}-\tau_{j}$ (10)

$m_{j}=\frac{F_{j+1}-F_{j}}{h_{j}},$ (11)

with $j=1,\dots,r-1$, the cubic interpolating spline function of (1) could be
defined by piecewise cubic polynomials that pass through consecutive points
and may be expressed as follows:

$S(\tau)=a_{j}+b_{j}(\tau-\tau_{j})+c_{j}(\tau-\tau_{j})^{2}+d_{j}(\tau-\tau_{j})^{3},$
(12)

where $\tau\in[\tau_{j},\tau_{j+1}]$ and $j=1,\dots,r-1$,
$a_{j}=F(\tau_{j})\equiv F_{j}$, $b_{j}=F^{\prime}(\tau_{j})$, and the other
two parameters defined as follows (see, Hagan and West, 2006):

$c_{j}=\frac{3m_{j}-b_{j+1}-2b_{j}}{h_{j}}$ (13)

$d_{j}=\frac{b_{j+1}+b_{j}-2m_{j}}{h_{j}^{2}}$ (14)

The cubic interpolating spline function has the following boundary conditions:
$S^{\prime}(\tau_{0})=0$ and $S^{\prime}(\tau_{r})=1$, where the pair
$(\tau_{0},\tau_{r})$ correspond to the minimum and maximum values of a given
distribution, which can also be $-\infty$ and $\infty$, respectively.

As previously noted, $F_{X}(\tau)$, which we seek to estimate, is a right-
continuous monotone increasing function and, as such, would require being
estimated with a corresponding monotone increasing function. This ensures that
the geometry of the underlying curve is preserved.

An excellent example is Hyman’s filter of the Fritsch-Butland algorithm
(Hyman, 1983; Fritsch and Butland, 1980; Butland, 1980; Fritsch and Carlson,
1980). It offers a monotonicity-preserving algorithm that could be implemented
in principle on almost any cubic-spline-like interpolation method. Hyman
(1983) first defines the following.

$b_{1}=\frac{(2h_{1}+h_{2})m_{1}-h_{1}m_{2}}{h_{1}+h_{2}},$ (15)

$b_{r}=\frac{(2h_{r-1}+h_{r-2})m_{r-1}-h_{r-1}m_{r-2}}{h_{r-1}+h_{r-2}},$ (16)

followed by the popular Hyman’s monotonicity constraint which ensures that no
spurious extrema are introduced in the interpolated function. The constrain is
given as follows (Hyman, 1983):

$b_{j}=\left\\{\begin{array}[]{ c l
}\mbox{min}(\mbox{max}(0,b_{j}),3\mbox{min}(m_{j-1},m_{j}))&\quad\textrm{if
}m_{j-1},m_{j}>0\\\
\mbox{max}(\mbox{min}(0,b_{j}),3\mbox{max}(m_{j-1},m_{j}))&\quad\textrm{if
}m_{j-1},m_{j}<0\end{array},\right.$ (17)

where the former choice is made when the curve is increasing and the latter
when the curve is decreasing. However, only the first constraint applies here
since the current study focuses on an increasing $F$. See Hagan and West
(2006) for further discussion of the properties of (17) and alternative
approaches.

Thus, given a set of binned data as illustrated in Table 1, where
$\tau_{0},\tau_{1},\dots,\tau_{r}$, are the thresholds (including minimum and
maximum) with which the continuous attribute was discretized, and
$(0,c_{1},\dots,c_{j-1},1)$ the corresponding cumulative relative frequencies,
the underlying CDF could be mirrored by monotonicity-preserving spline
interpolation of these $r$ pairs of points (see, e.g., von Hippel et al.,
2017). Obtaining the underlying CDF of the binned data this way also makes it
possible to realize the underlying PDF, considering that the derivative of the
spline interpolator is an estimator of the density (Lii and Rosenblatt, 1975).

##  6 Simulation Study

The MICS approach was examined in a simulation study against results from the
true theoretical distribution and compared to related approaches. Four
distinct underlying distributions (shapes) were considered with a relatively
small and large number of observations ($n=100\mbox{ and }1000$). These
include the normal: $y^{*}_{i1}\sim\mathcal{N}(3,\,1)$, gamma:
$y^{*}_{i2}\sim\mathcal{G}(1,\,2)$, generalized extreme value:
$y^{*}_{i3}\sim\mathcal{G}v(1,\,2,\,0)$ and triangular:
$y^{*}_{i4}\sim\mathcal{T}(0,\,1,\,0.5)$ distributions. The finite ranges of
observations $(\tau_{0},\tau_{r})$ considered for the various distributions
are respectively $(0,10)$, $(0,8)$, $(-4,35)$ and $(0,1)$.

The generated continuous distributions were discretized using $r=6$
equidistant cut-points and subsequently compressed to frequency tables,
yielding a binned dataset, $y^{b}_{ik}$, $k=1,\dots 4$, with $r+1$ classes.
The entire process was replicated a thousand times, resulting in a thousand
pairs of datasets $\\{(y^{*}_{il},y^{b}_{ik}),\ i=1,\dots,n,\ k=1,\dots,4\\}$
for the four specified distributions.

With the datasets generated, we proceeded with the estimation of the CDFs and
PDFs of the individual distributions (i) directly from the continuous sampling
distribution $y^{*}_{ik}$ and (ii) from the binned data $y^{b}_{ik}$, using
the MICS, heuristic, eCDF, and kernel estimation approaches. These were
subsequently compared to the theoretical distributions obtained with prior
specified parameters. To simplify things, we shall pay particular attention to
estimating quantiles, means, and standard deviations from the various
distributions. The estimated quantiles with each method include .25, .50 and
.75. In other words, quantiles obtained from the continuous sampling
distributions and those from the binned data are compared with the true
theoretical quantiles, the same with means and standard deviations.

Figure 3:  Pairwise differences of quantile estimates $(Q_{1},Q_{2},Q_{3})$ of
binned data generated from an underlying continuous data with sample size
$n=1000$. D denotes estimates from the sampling distribution, and S, H, E, and
K are estimates from the spline, heuristic, eCDF, and kernel estimators,
respectively. Dotted horizontal red lines represent the mark of over and
underestimation of the underlying measure. Figure 4:  Pairwise differences of
mean $\Delta(\mu)$ and standard deviation $\Delta(\sigma)$ estimates of binned
data generated from an underlying continuous data with sample size $n=1000$. D
denotes estimates from the sampling distribution, and S, H, E, and K are
estimates from the spline, heuristic, eCDF, and kernel estimators,
respectively. Dotted horizontal red lines represent the mark of over and
underestimation of the underlying measure.

For easy presentation of the results, the pairwise differences of the
estimated population quantiles $\Delta(Q)$, means $\Delta(\mu)$, and standard
deviations $\Delta(\sigma)$, with reference to the true theoretical values,
were obtained from a large sample size $(n=1000)$ and a relatively small
sample size $(n=100)$. Thus, assuming $Q_{q}$, $\mu$ and $\sigma$ respectively
denote the theoretical quantiles, means, and standard deviations;
$\hat{Q}_{q}$, $\hat{\mu}$, and $\hat{\sigma}$ respectively denote the
estimated quantiles, means, and standard deviations from the sampling
distribution or from the binned data, the paired differences are obtained as
follows,

$\Delta(Q_{q})=Q_{q}-\hat{Q}_{q},\quad q\in(.25,\ .50,\ .75)$ (18)

$\Delta(\mu)=\mu-\hat{\mu}$ (19)

$\Delta(\sigma)=\sigma-\hat{\sigma},$ (20)

Figure 5:  Pairwise differences of quantile estimates $(Q_{1},Q_{2},Q_{3})$ of
binned data generated from an underlying continuous data with sample size
$n=100$. D denotes estimates from the sampling distribution, and S, H, E, and
K are estimates from the spline, heuristic, eCDF, and kernel estimators,
respectively. Dotted horizontal red lines represent the mark of over and
underestimation of the underlying measure. Figure 6:  Pairwise differences of
mean $\Delta(\mu)$ and standard deviation $\Delta(\sigma)$ estimates of binned
data generated from an underlying continuous data with sample size $n=100$. D
denotes estimates from the sampling distribution, and S, H, E, and K are
estimates from the MICS, heuristic, eCDF, and kernel estimators, respectively.
Dotted horizontal red lines represent the mark of over and underestimation of
the underlying measure.

The calculated pairwise differences of the various metrics for the different
distributions are shown in Figures 3 and 4, where D, for instance, denotes the
pairwise difference of estimates from the sampling distribution and the
theoretical values. Similarly, S, H, E, and K denote the pairwise differences
of estimates from the MICS, heuristic, eCDF, and kernel estimators,
respectively. Dotted red lines in the plots indicate if the different methods
tend to overestimate or underestimate the actual values.

As observed in Figure 3, in addition to approximating the actual theoretical
quantities, the estimated quantiles with the MICS approach (S) outperform
those from the popular and widely used approaches. They are also very close to
those from the empirical distribution (D). It adapts very well to both mild
and strong departures from normality (captured by the different
distributions). In all four settings, the MICS approach further recovers much
information about the spread of the underlying data that was lost due to
binning.

Similar patterns were also observed with the mean and standard deviation
estimates reported in Figure 4. Besides yielding mean and standard deviation
values that well approximate those from both the theoretical and empirical
distributions (D), it stays consistent irrespective of changes in the shapes
of the underlying distributions. The performance of the estimators was further
assessed using a relatively small sample size $(n=100)$. Results are given in
Figures 5 and 6, where we see almost a replica of the previous results, except
for larger deviations from the theoretical values due to sampling errors.
Meaning, that the MICS approach also performs very well, compared to competing
methods, for smaller sample sizes. In general, the MICS’s improvement over the
heuristic and the other approaches is more pronounced with skewed
distributions.

##  7 Real World Data Applications

###  7.1 Train Punctuality

For a real-world application of the discussed approaches, we consider the data
on Deutsche Bahn (DB) delays (Deutsche Bahn, 2022). In the year 2022, only
around two-thirds of all long-distance trains (ICE and IC) reportedly reached
their destination on time (65.2%). This has been attributed to poor and scarce
infrastructure, intensive construction activity, and (from the second quarter)
a rapidly growing traffic volume in long-distance and local transport, which
all put rail operations under pressure. According to official records, long-
distance transport went 10% points below the previous year (75.2%) in
punctuality. Table 2 shows the proportion of punctual stops in relation to all
en route and final stops. A stop is considered on time by Deutsche Bahn if the
arrival time is less than 6 or 16 minutes behind schedule, respectively. The
corresponding aggregated data in terms of relative frequencies are published
by DB on a monthly basis (Deutsche Bahn, 2022). The punctuality statistics in
Table 2 were obtained from more than 800,000 journeys by DB passenger trains
(PT) in a month, all stops of more than 20,000 monthly trips in long-distance
traffic (LT), and approximately 780,000 monthly trips in local traffic (RT) –
including all S-Bahn trains.

Table 2:  Aggregated delays of all DB passenger trains (PT), long-distance trains (LT), and regional trains (RT) as published by Deutsche Bahn; % counts displayed. | Deutsche Bahn Delays in 2022
---|---
| PT | LT | RT
Month |  $\leq$ 5:59 Min |  $\leq$ 15:59 Min |  $\leq$ 5:59 Min |  $\leq$ 15:59 Min |  $\leq$ 5:59 Min |  $\leq$ 15:59 Min
Jan | 95.1 | 99.0 | 80.9 | 92.2 | 95.5 | 99.2
Feb | 93.9 | 98.6 | 76.0 | 88.7 | 94.4 | 98.8
Mar | 93.3 | 98.5 | 71.1 | 85.9 | 93.9 | 98.8
Apr | 93.2 | 98.5 | 69.1 | 85.2 | 93.9 | 98.8
May | 91.5 | 97.8 | 62.7 | 79.4 | 92.3 | 98.3
Jun | 87.7 | 96.7 | 58.0 | 77.3 | 88.5 | 97.3
Jul | 88.7 | 97.2 | 59.9 | 79.1 | 89.6 | 97.7
Aug | 88.6 | 97.1 | 56.8 | 76.8 | 89.5 | 97.7
Sep | 90.6 | 97.8 | 62.8 | 81.7 | 91.4 | 98.3
Oct | 89.8 | 97.7 | 63.2 | 82.0 | 90.6 | 98.1
Nov | 89.5 | 97.7 | 61.1 | 80.6 | 90.3 | 98.1
Dec | 89.8 | 97.5 | 61.3 | 79.6 | 90.7 | 98.0
Figure 7:  Estimated cumulative distributions of the proportion of punctual stops in relation to all en route and final stops of the three classes of Deutsche Bahn (PT, LT, and RT) for the months of January and June. Figure 8:  Estimated population characteristics of delays in the months of January and June of Deutsche Bahn passenger trains (PT), long-distance trains (LT), and Regional trains (RT). Where S, H, E, and K are spline, heuristic, eCDF, and kernel estimates respectively. Table 3:  Estimated population characteristics of the monthly and averaged delays of Deutsche Bahn long-distance trains (LT). The interquartile range (IQR) and the mean from the MICS (S), heuristic (H), eCDF (E), and kernel (K) estimators are reported. | Spline | Heuristic | eCDF | Kernel
---|---|---|---|---
Month | IQR | $\mu$ | IQR | $\mu$ | IQR | $\mu$ | IQR | $\mu$
Jan | 2.97 | 13.31 | 3.09 | 19.78 | 0.00 | 17.02 | 2.63 | 10.31
Feb | 3.54 | 18.19 | 3.29 | 26.05 | 0.00 | 22.62 | 0.44 | 13.43
March | 4.62 | 22.15 | 5.88 | 31.16 | 5.07 | 27.01 | 0.22 | 15.96
April | 5.02 | 23.19 | 6.86 | 32.51 | 5.08 | 28.13 | 5.04 | 16.64
May | 7.03 | 31.16 | 10.37 | 42.72 | 5.08 | 36.43 | 5.04 | 21.73
June | 8.34 | 34.20 | 11.65 | 46.66 | 5.10 | 39.42 | 32.03 | 23.68
July | 7.35 | 31.73 | 10.78 | 43.49 | 5.10 | 36.99 | 5.05 | 22.11
Aug | 8.73 | 34.93 | 11.90 | 47.60 | 5.10 | 40.13 | 62.10 | 24.15
Sep | 6.35 | 28.16 | 9.46 | 38.92 | 5.09 | 33.36 | 5.04 | 19.82
Oct | 6.24 | 27.74 | 9.30 | 38.38 | 5.09 | 32.93 | 5.04 | 19.56
Nov | 6.77 | 29.70 | 10.08 | 40.90 | 5.10 | 34.94 | 5.04 | 20.81
Dec | 7.07 | 30.99 | 10.45 | 42.53 | 5.09 | 36.25 | 5.04 | 21.63
Average | 6.17 | 27.12 | 8.59 | 37.56 | 4.24 | 32.10 | 11.06 | 19.15

Given the first month (January), for instance, it is seen that 95.1% and 99.0%
of all passenger traffic trains respectively showed up not more than 5 and 15
minutes after schedule. The long-distance trains had 80.9% and 92.2%, while
the regional trains had 95.5% and 99.2% respectively. However, with just two
pairs of data points each not much could be said about the actual distribution
of the trains’ punctuality. For a better estimate of the underlying
distribution, one at least considers two more extreme cases (data points),
namely, 0% and 100% delays. The resulting four pairs of points could then be
viewed as real points on the unobserved CDFs of train punctuality and could
help to retrace the shape of the underlying distribution. While it makes sense
in this scenario to assume that 0% delays correspond to 0 minutes, one problem
is how to fix the upper limit, $\tau_{r}$, at which delays correspond to 100%.
Per experience, delays may take up to 3 hours (180 mins). We would consider
this the train delays’ upper limit for estimating the underlying train delay
CDFs. The obtained CDFs (via MICS) for the months of January and June for each
of the three train types (PT, LT, and RT) are shown in Figure 7, with the last
points at (300, 1) an additional pseudo observation used to stabilize the
curve at the endpoints. On the one hand, the shapes of the CDFs of PT and RT
look quite the same in January and June, suggesting a similar delay time
frame. On the other hand, some differences are seen with the CDFs of the long-
distance trains (LT), and, particularly for June delays, the estimated CDF
looks odd. Interpolating straight lines would (presumably) produce more
realistic results here.

Figure 9:  Estimated population characteristics of delays of Deutsche Bahn
long-distance train (LT) in the month of July using different topmost upper
limits ($\tau_{r}$). Where S, H, E, and K are MICS, heuristic, eCDF, and
kernel estimates respectively.

Also, while we could consider $\tau_{r}=180$ mins as the train delays upper
limit, it could also be interesting to investigate what difference it makes
using other upper limits on the different estimation methods. This is a well-
known problem in several practical applications of binned data (see, e.g.,
Jargowsky and Wheeler, 2018; Ho and Reardon, 2012; Atkinson et al., 2011).
Thus, we investigate the effect of varying $\tau_{r}$ on the estimated
population characteristics and consider the different estimation approaches.
For this, we obtained the quantiles ($Q1$, $Q2$ and $Q3$) and means $(\mu)$
for the long-distance train (LT) given $\tau_{r}=\\{30,\ 60,\ 120,\ 180\\}$.
Results are given in Figure 9. As observed, quantile estimates (Q1, Q2, and
Q3) from MICS (S), heuristic (H), and eCDF (E) approaches remain invariant to
increasing $\tau_{r}$. On the other hand, the kernel estimates exhibited large
variations with the median and third quartile, as indicated by up and down
movements in the plots. The estimated means, on the other hand, increase
monotonically with each rise in $\tau_{r}$. This is not surprising considering
that even under normal circumstances, means are generally known to be affected
by extreme values. In this instance, making decisions based on quantile
estimates would be a good idea, rather than relying on the means.

Continuing further with our fixed $\tau_{r}$ setting, the means and IQRs from
the different approaches were subsequently obtained and compared in Figure 8,
while Table 3 presents both the monthly and averaged estimates for the long-
distance train (LT). Results show some form of dissimilarities between the
different methods. As with some other estimates, the MICS estimates capture a
realistic situation across the months, with values rising from January,
climaxing in June till August (summer), dropping afterward, and picking up
again in December. Furthermore, estimates from E and K are smaller than those
from S and H, mainly due to less variability captured by the former,
particularly with E which is basically a step function with constant values
recorded at different quantile steps. The long-distance train (LT) out of the
three train categories had the largest estimated mean delays and widest spread
(in terms of IQR). On average (i.e., annually), the estimated mean delays in
minutes and IQR in parenthesis for 2022 using MICS are 27.12 (6.17); 37.56
(8.59), 32.10 (4.24), and 19.15 (11.06) for the heuristic, eCDF and kernel
estimators, respectively.

###  7.2 Distance-Time to Work

Table 4:  Upper segment: % frequencies of the distance in kilometers traveled by all workers (All), self-employed workers (SE), and the employees (EM) to work in Germany. Lower segment: the amount of time in minutes spent en route to the workplace. Actual observations are given in groups with the lower (LL) and upper class limits (UL) also indicated. |  | Commuters
---|---|---
LL | UL | All-workers (%) | Self-employed (%) | Employees (%)
Distance (km) |  |  |
0 | 5 | 27.5 | 47.7 | 26.3
5 | 10 | 22.6 | 20.0 | 22.7
10 | 25 | 30.1 | 20.4 | 30.8
25 | 50 | 14.6 | 7.7 | 15.1
50 | 100 | 5.2 | 4.2 | 5.2
Time (min)
0 | 10 | 21.4 | 41.6 | 20.2
10 | 30 | 50.9 | 41.7 | 51.5
30 | 60 | 22.7 | 13.1 | 23.2
60 | 90 | 5.0 | 3.7 | 5.1
Figure 10:  Estimated cumulative distributions of distance in kilometers to
work and time in minutes to work of the self-employed (SE), employees (EM),
and all workers (All) in Germany.

The next real data example deals with a segment of micro-census data on the
distance and time required to get to work by workers in Germany in 2020.
According to Destatis (2022), the term ‘micro-census’ means “small population
census” and represents the largest annual household survey of official
statistics in Germany. The survey has been carried out jointly by the federal
and state statistical offices since 1957. Hence, with roughly 810,000 people
in around 370,000 private households and communal accommodations, about 1% of
Germany’s population was surveyed about their working and living conditions.
The statistics office had the federal territory divided into areas with
approximately the same number of apartments, giving all districts/cities equal
chances of being selected for the survey. Thus, selection for the survey
depends on a respondent’s living apartment being selected. Table 4 shows the %
frequencies of distance in kilometers and the time in minutes traveled to work
by all workers (All), the self-employed (SE), and the employees (EM).

Similar to the previous example, the distribution of distance and time to work
for all workers, the self-employed, and the employees were estimated using
MICS and the other estimation methods. The estimated CDFs are shown in Figure
10, with the SE curve steeper than EM in both distance and time CDFs. The
means, first, second, and third quartiles from the estimated distribution were
obtained using the MICS and other estimation methods (see Figure 11). While
MICS and heuristic results do not differ so much, both differ from the other
two approaches and are more and more pronounced with the quantile estimates.
With the exception of the estimated mean from the eCDF, the estimated mean
from the other three is nearly uniformly shaped. In summary, the obtained
results suggest that employees in Germany spend more time and distance to work
than self-employed individuals.

Figure 11:  Estimated population characteristics of the distance in kilometers
(top row) and time in minutes (bottom row) of employees (EM), self-employed
(SE), and all workers (All) in Germany. Where S, H, E, and K are spline,
heuristic, eCDF, and kernel estimates, respectively.

##  8 Discussion

Nonparametric density estimation provides an alternative means of obtaining
crucial estimates of different population characteristics in situations where
the parametric approach fails. In this study, we have discussed how such could
be used to estimate the distribution of binned data. Rather than handling
binned data in an aggregated form, leading to the loss of valuable
information, we would rather estimate the underlying (unobserved) CDF with
monotone-preserving cubic splines, often leading to more accurate and reliable
estimates. This approach frequently outperforms other widely used approaches,
as shown in our simulation studies and real data examples.

However, as observed in our study, a typical open question while estimating
the population characteristics of binned data is how to handle open-ended
classes, a common problem in various large-scale surveys. In situations with
no standard benchmark limits to fall on, this could present a daunting
challenge to making reliable inferences from binned data. As noted earlier in
Section 1, common suggestions, such as those in studies of income
distributions, include assuming that the quantity in the upper tail of the
distribution follows a Pareto distribution (von Hippel et al., 2015), which we
consider a costly assumption. Jargowsky and Wheeler (2018) proposed an
improvement on this approach by introducing the mean-constrained integration
over brackets (MCIB) approach, which retains the same Pareto distribution
assumption for the upper tail but employs a different approach for estimating
the shape parameter of the specified distribution. Alternatively, Rizzi et al.
(2015) demonstrated how binned data can be ungrouped using the penalized
composite link model. This approach can handle wide, open-ended intervals and
incorporate extra information about how the tail area could be occupied,
assuming such information is available. However, estimation with this approach
is based on an iteratively reweighted least-squares algorithm, with the
selection of a smoothing parameter through the minimization of AIC, which is
inherently parametric. Another proposal by von Hippel et al. (2017) involves
constraining the parameter shaping the top bin to the grand mean if known,
such that the mean of the fitted distribution matches the grand mean. When the
grand mean is unknown, an ad hoc estimate is obtained by temporarily setting
the upper bound in question to twice the lower bound, calculating the mean of
a step PDF fit to all bins, unbinding the top bin, and proceeding as though
the mean were known.

Notwithstanding these approaches, the inherent subjectivity associated with
them leaves the initial problem not entirely solved. In our study, we
investigated the effect of varying upper limits in our Deutsche Bahn delay
study, for instance, and observed that some population characteristics are not
invariant to the choice of missing $\tau_{r}$. The mean, in particular, is
most affected. However, we did note that quantile estimates remain invariant
to the choice of $\tau_{r}$ and should be reported in scenarios where the mean
estimate could be influenced by extreme values.

In conclusion, we consider the MICS density estimation presented in this study
a promising non-parametric approach for handling binned data. The certainty of
having the interpolating curve pass through all data points is assured, which
is preferable to fitting inflexible curves that fail to capture the complex
variation of actual data. Moreover, the MICS approach easily adapts to
different underlying distributional shapes, as demonstrated in the simulation
study, and is not affected by whether bin sizes are equally or unequally
spaced.

## References

  * Atkinson et al. (2011) Atkinson, A. B., T. Piketty, and E. Saez (2011).  Top incomes in the long run of history.  Journal of Economic Literature 49(1), 3–71.
  * Azzalini (1981) Azzalini, A. (1981).  A note on the estimation of a distribution function and quantiles by a kernel method.  Biometrika 68, 326–328.
  * Bae et al. (2005) Bae, W., H. Choi, B. Park, and C. Kim (2005).  Smoothing techniques for the bivariate kaplan-meier esstimator.  Communications in Statistics - Theory and Methods 34, 1659–1674.
  * Bae et al. (2014) Bae, W., R. Kim, and C. Kim (2014).  Nonparametric estimation of distribution function using bezier curve.  Communications for Statistical Applications and Methods 21, 105–114.
  * Butland (1980) Butland, J. (1980).  A method of interpolating reasonable-shaped curves through any data.  In Proceedings of Computer Graphics, pp. 409–422. Online Publications Ltd.
  * Charles et al. (2010) Charles, J., S. Sun, and C. Martin (2010).  Cumulative distribution estimation via control theoretic smoothing splines.  In Three Decades of Progress in Control Sciences, pp. 95–104.
  * Destatis (2022) Destatis (2022).  Berufspendler.  https://www.destatis.de/DE/Themen/Arbeit/Arbeitsmarkt/Erwerbstaetigkeit/Tabellen/pendler1.html.  Online; accessed 29 June 2022.
  * Deutsche Bahn (2022) Deutsche Bahn (2022).  Erläuterung Pünktlichkeitswerte Dezember 2022.  https://www.deutschebahn.com/de/konzern/konzernprofil/zahlen_fakten/puenktlichkeitswerte-6878476.  Online; accessed 29 January 2023.
  * Eubank (1988) Eubank, R. (1988).  Spline Smoothing and Nonparametric Regression.  New York: Marcel Dekker.
  * Fritsch and Butland (1980) Fritsch, F. and J. Butland (1980).  An Improved Monotone Piecewise Cubic Interpolation Algorithm.  Lawrence Livermore National Laboratory preprint UCRL-85104.
  * Fritsch and Carlson (1980) Fritsch, F. and R. Carlson (1980).  Monotone piecewise cubic interpolation.  SIAM J. Numer. Anal 17, 238–246.
  * Funke and Palmes (2016) Funke, B. and C. Palmes (2016).  A note on estimating cumulative distribution functions by the use of convolution power kernels.  In Statistics and Probability Letters, Volume 121, pp. 90–98.
  * Ghosh et al. (2018) Ghosh, S., C. Burns, D. Prager, L. Zhang, and G. Hui (2018).  On nonparametric estimation of the latent distribution for ordinal data.  Computational Statistics & Data Analysis 119, 86–98.
  * Hagan and West (2006) Hagan, P. and G. West (2006).  Interpolation methods for curve construction.  Applied Mathematical Finance 13(2), 89 – 129.
  * Heitjan (1989) Heitjan, D. (1989).  Inference from grouped continuous data: A review: Rejoinder.  Statistical Science 4, 182–183.
  * Henderson and Parmeter (2015) Henderson, D. and C. Parmeter (2015).  Applied Nonparametric Econometrics.  New York: Cambridge University Press.
  * Henson (1967) Henson, M. (1967).  Trends in the Income of Families and Persons in the United States.  pp. 1947–1964. Washington, D.C: U.S. Department of Commerce, Bureau of the Census.
  * Ho and Reardon (2012) Ho, A. D. and S. F. Reardon (2012).  Estimating achievement gaps from test scores reported in ordinal ’proficiency’ categories.  Journal of Educational and Behavioral Statistics 37(4), 489–517.
  * Huang et al. (2016) Huang, J., X. Wang, X. Wu, and L. Zhou (2016).  Estimation of a probability density function using interval aggregated data.  Journal of Statistical Computation and Simulation 86, 3093–3105.
  * Hyman (1983) Hyman, J. (1983).  Accurate monotonicity preserving cubic interpolation.  SIAM J. Sci. Stat. Comput 4(4), 645 – 654.
  * Jargowsky and Wheeler (2018) Jargowsky, P. A. and C. A. Wheeler (2018).  Estimating income statistics from grouped data: Mean-constrained integration over brackets.  Sociological Methodology 48(1), 337–374.
  * Jones et al. (1996) Jones, M., J. Marron, and S. Sheather (1996).  A brief survey of bandwidth selection for density estimation.  Journal of the American Statistical Association 91(433), 401–407.
  * Kim (1996) Kim, C. (1996).  Nonparametric density estimation via the bezier curve.  In ASA Proceedings of the Section on Statistical Graphics, pp. 25–28.
  * Kim et al. (1999) Kim, C., W. Kim, C. Hong, B. Park, and M. Jeong (1999).  Smoothing techniques via the bezier curve.  Communications in Statistics - Theory and Methods 28, 1577–1596.
  * Kim et al. (2003) Kim, C., W. Kim, B. U. Park, and J. Lim (2003).  Bezier curve smoothing of the kaplan-meier estimator.  Annals of the Institute of Statistical Mathematics 55, 359–367.
  * Kim and Park (2012) Kim, C. and H. J. Park (2012).  On the selection of bezier points in bezier curve smoothing.  Korean Journal of Applied Statistics 25, 1049–1058.
  * Lehmann et al. (1999) Lehmann, T., C. Gonner, and K. Spitzer (1999).  Survey: Interpolation methods in medical image processing.  IEEE Transactions on Medical Imaging 18(11), 1049.
  * Lii and Rosenblatt (1975) Lii, K. and M. Rosenblatt (1975).  Asymptotic behavior of a spline estimate of a density function.  Comput. Math. Appl. 1, 223 – 235.
  * Liu and Yang (2008) Liu, R. and L. Yang (2008).  Kernel estimation of multivariate cumulative distribution function.  Journal of Nonparametric Statistics 20(8), 661–677.
  * Loader (1999a) Loader, C. (1999a).  Bandwidth selection: classical or plug-in?  Annals of Statistics 27(2), 415 – 438.
  * Loader (1999b) Loader, C. (1999b).  Local Regression and Likelihood.  London: Springer.
  * McDonald and Ransom (2008) McDonald, J. and M. Ransom (2008).  The generalized beta distribution as a model for the distribution of income: Estimation of related measures of inequality.  In Modeling Income Distributions and Lorenz Curves, pp. 147–166. New York: Springer.
  * Minoiu and Reddy (2014) Minoiu, C. and S. Reddy (2014).  Kernel density estimation on grouped data: The case of poverty assessment.  The Journal of Economic Inequality 12(2), 163–189.
  * Nadaraya (1964) Nadaraya, E. (1964).  Some new estimators for distribution functions.  Theory of Probability and its Applications 9, 497–500.
  * Park (2018) Park, K. I. (2018).  Fundamentals of Probability and Stochastic Processes with Applications to Communications.  Springer.
  * Parzen (1962) Parzen, E. (1962).  On estimation of a probability density function and mode.  The Annals of Mathematical Statistics 33(3), 1065–1076.
  * Reyes et al. (2019) Reyes, M., M. Francisco-Fernández, R. Cao, and D. Barreiro-Ures (2019).  Kernel distribution estimation for grouped data.  statistics and operations research transactions 43(2), 259–288.
  * Rizzi et al. (2015) Rizzi, S., J. Gampe, and P. Eilers (2015).  Efficient estimation of smooth distributions from coarsely grouped data.  Am J Epidemiol. 182(2), 138–147.
  * Rosenblatt (1956) Rosenblatt, M. (1956).  Remarks on some nonparametric estimates of a density function.  The Annals of Mathematical Statistics 27(3), 832 – 837\.
  * Scott and Sheather (1985) Scott, D. and S. Sheather (1985).  Kernel density estimation with binned data.  Commun Stat Theory Methods 14, 1353–1359.
  * Silverman (1986) Silverman, B. (1986).  Density Estimation for Statistics and Data Analysis.  York: Chapman and Hall.
  * von Hippel et al. (2017) von Hippel, P., D. Hunter, and M. Drown (2017).  Better estimates from binned income data: Interpolated cdfs and mean-matching.  Sociological Science 4, 641–655.
  * von Hippel et al. (2015) von Hippel, P., S. Scarpino, and I. Holas (2015).  Robust estimation of inequality from binned incomes.  Sociological Methodology 46, 212–51.
  * Wang and Wertelecki (2013) Wang, B. and W. Wertelecki (2013).  Density estimation for data with rounding errors.  Computational Statistics & Data Analysis 65, 4 – 12.
  * Wasserman (2006) Wasserman, L. (2006).  All of Nonparametric Statistics.  London: Springer.
  * Xue and Wang (2009) Xue, L. and J. Wang (2009).  Distribution function estimation by constrained polynomial spline regression.  Journal of Nonparametric Statistics 22(4), 443–457.
  * Yamato (1973) Yamato, H. (1973).  Uniform convergence of an estimator of a distribution function.  Bulletin of Mathematical Statistics 15, 69–78.