1
Density Estimation

1.1 Introduction

Use of sampled observations to approximate distributions has a long history. An important milestone was Pearson (1895, 1902a, 1902b), who noted that the limiting case of the hypergeometric series can be written as in the equation below and who introduced the Pearsonian system of probability densities. This is a broad class given as a solution to the differential equation

The different families of densities (Type I-VI) are found by solving this differential equation under varying conditions on the constants. It turns out that the constants are then expressible in terms of the first four moments of the probability density function (pdf)  f, so that they can be estimated given a set of observations using the method of moments; see Kendall and Stuart (1963).

If the unknown pdf  f is known to belong to a known parametric family of density functions satisfying suitable regularity conditions, then the maximum likelihood (MLE; Fisher 1912, 1997) can be used to estimate the parameters of the density, thereby estimating the density itself. This method has very powerful statistical properties, and continues to be perhaps the most popular method of estimation in statistics. Often, the MLE is the solution to an estimating equation, as is also the case for the method of least squares. These procedures then come under the general framework of M-estimation. Two other related approaches that use ranks of the observations are the so-called L-estimation and R-estimation, where the statistics are respectively linear combinations of the order statistics or of their ranks. These estimation methods are covered in many standard textbooks. Some examples are Rao (1973, chapters 4 and 5), Serfling (1986, chapters 7, 8, and 9), and Sen and Srivastava (1990).

1.1.1 Orthogonal polynomials

Yet another approach worth mentioning here is the use of Orthogonal polynomials (see Szego 2003). In this method, the unknown density is approximated by a sum of weighted linear combinations of a set of basis functions. encov (1962) provides a general description whereas other reviews are in Wegman (1972) and Silverman (1980). Additional background information and further references can be found in Beran et al. (2013, Chapter 3) and Kendall and Stuart (1963, Chapter 6). The essential idea behind the use of Orthogonal polynomials is as follows (see Rosenblatt 1971):

Suppose that the pdf

belongs to the space of all square integrable functions with respect to the weight function G, i.e.,

holds, where denotes the real line. Also, let {Gl(x)} be a complete and orthonormal sequence of functions in . Then  f admits an expansion

which converges to  f in , where al is defined as

This formula immediately suggests an unbiased estimator of the coefficient al using sampled observations, followed by a substitution in the expansion for  f.

As an example, we take a brief look at the Gram-Charlier series representation followed by a further extension due to Schwartz (1967). The Gram-Charlier series of Type A is based on Hermite polynomials Hl and the standard normal pdf ?. Note that, for Edgeworth expansion based methods, one would consider the Fourier transform of the product Hl(x)?(x) and move on to an expansion that uses the cumulant generating function (see Kendall and Stuart 1963).

First of all consider the pdf  f such that it can be expressed as

For conditions under which this is valid, see two theorems due to Cramér quoted in Kendall and Stuart (1963, pp. 161-162) as well as some historical notes in Cramér (1972).

Here ? is the standard normal pdf, i.e.,

and Hl is the Hermite polynomial of degree l, i.e.,

Using the orthogonality property of the Hermite polynomials, i.e.,

we have

In other words, the coefficients cl are

Due to previous detailed work by Chebyshev, the Hermite polynomials are also known as the Chebyshev-Hermite polynomials. In fact, contributions of Laplace are also known. See Sansone (2004) and Szego (2003) for additional information.

The above formula for cl implies that these coefficients may be estimated from a given set of observations X1, ., Xn from  f as sample means of Hermite polynomials, i.e.,

Since with increasing l, estimation of higher order moments are involved, this method however is not optimal. From a statistical view-point, one option is to consider a finite sum.

To this end, Schwartz (1967) considers a pdf  f that is square integrable (or simply bounded) and seeks to give an approximation of the form

where Mn is a sequence of integers depending on the sample size n, dl, n are estimated from observed data, and Gl are Hermite functions

The Hermite functions Gl(x) form a complete orthonormal set over the real line. Examples of Hermite polynomials and Hermite functions are in Figure 1.1 and Figure 1.2. Moreover, due to a theorem of Cramér (see Schwartz 1967), |Gl(x)| is bounded above by a constant that does not depend on x or l. Since  f is square integrable,  f can be expanded (orthogonal series expansion) as

where

Schwartz (1967) proposes the estimator

where Mn 8 and Mn = o(n) as n 8 and the coefficients are estimators based on the sample means of Hermite functions

Under some conditions on the rth derivative (r = 2) of  f(x)?(x), Schwartz (1967) derives asymptotic properties of his estimator including the rate of convergence to zero of the mean integrated squared error (MISE).

Figure 1.1 Rescaled Hermite polynomials of degree l for l = 0, 1, 2 and the corresponding Hermite functions (right) Gl(x). These functions are related via the relation , where , where Hl is the Hermite polynomial of degree l.

Figure 1.2 Rescaled Hermite polynomials of degree l for l = 3, 4, 5 and the corresponding Hermite functions (right) Gl(x). These functions are related via the relation , where , where Hl is the Hermite polynomial of degree l.

There are various textbooks and review papers that give excellent overvews of nonparametric density estimation techniques. Basic developments and related information can, for instance, be found in Watson and Leadbetter (1964a, 1964b), Shapiro (1969), Rosenblatt (1971), Bickel and Rosenblatt (1973), Rice and Rosenblatt (1976), Tapia and Thompson (1978), Wegman (1982), Silverman (1986), Hart (1990), Jones and Sheather (1991), Müller and Wang (1994), Devroye (1987), Müller (1997), Loader (1999), and Heidenreich et al. (2013). Various textbooks have addressed applied aspects and included various theoretical results on general kernel smoothing methods. Some examples are Bowman and Azzalini (1997), Wand and Jones (1995), Simonoff (1996), Scott (1992), Thompson and Tapia (1987), and others.

In this chapter, we focus on a selection of ideas for density estimation with independently and identically distributed (iid) observations, restricting ourselves to continuous random variables. We start with the univariate case and the multivariate case is mentioned in the sequel.

Let X, X1, X2, ., Xn be iid real-valued univariate continuous random variables with an absolutely continuous cumulative distribution function

where  f(x) denotes the probability density function (pdf). The pdf  f will be assumed to be a three times continuously differentiable function with finite derivatives. Further conditions will be added in the sequel.

The problem is a nonparametric estimation of , using X1, X2, ., Xn.

1.2 Histograms

The most widely used nonparametric density estimation method is the histogram, especially for univariate random variables. The idea has a long history and the name "histogram" seems to have been used for the first time by Karl Pearson 1895). Basic information on the use of the histogram as a graphical tool to display frequency distributions can be found in any elementary statistical textbook.

Construction of a histogram proceeds as follows. We consider the univariate case. Let

be a partition of the real line...