Estimating the population median from a kernel density estimator

Question

I have a 1-d kernel density estimate in the form of two vectors:
x_grid is a vector of x-values at which the density function was sampled
density is a vector of corresponding density values (estimated by a kernel density estimate) at each x value in the x_grid.

If I want the median value of this density function ($m$ where $\int_{-\infty}^{m}{f(x) =0.5}$), then I can approximate this by taking the integral (using a numerical technique such as Trapezoidal rule) from $-\infty$ to $x_0, x_1, x_2, ...$ until I get (close enough to) $0.5$.

Essentially, I'm just summing up the areas of trapezoids under the curve until I reach a cumulative sum of $0.5$.

How can I do something similar, but for the mean/expected value of the density function?

Answer 1

Here are some experiments using the default kernel density estimator in R. It is designed in a way that makes estimation of the median practical for reasonably large samples.

I have chosen to illustrate with a sample $V_1, V_2, \dots, V_{500}$ of size $n = 500$ from $\mathsf{Gamma}(\text{shape}=5, \text{rate}=0.1),$ a right-skewed distribution. Thus the population mean is $\mu = 5/0.1 = 50,$ the population mode is $\delta = (5-1)/0.1 = 40,$ and the population median (from R) is $\eta = 46.71.$ It is feasible to get reasonable estimates of all three of these measures of centrality from the default KDE in R. [If needed, you can find information about gamma distributions on Wikipedia.]

qgamma(.5, 5, .1)
[1] 46.70909

Exploring the output from the KDE: The function density(v) produces two vectors: one with evenly spaced $x$-values and the other with corresponding heights ($y$-values).

set.seed(418)  # for reproducibility
v = rgamma(500, 5, .1)
density(v)

Call:
        density.default(x = v)

Data: v (500 obs.);     Bandwidth 'bw' = 5.43

       x                 y            
 Min.   : -7.166   Min.   :1.675e-06  
 1st Qu.: 33.912   1st Qu.:2.219e-04  
 Median : 74.990   Median :2.524e-03  
 Mean   : 74.990   Mean   :6.080e-03  
 3rd Qu.:116.068   3rd Qu.:1.230e-02  
 Max.   :157.146   Max.   :2.018e-02  

Here is a density histogram of the sample, along with the KDE.

lbl = "Density of GAMMA(5,.1) with KDE for n = 500"
hist(v, prob=T, col="skyblue2", main=lbl);  rug(v)
  lines(density(v), lwd=2)

It is convenient to give names to the two vectors. Also, we verify that the points along the horizontal grid are equally spaced (we call the common distance wd), which makes numerical integration easy, allowing us to see that the area under the KDE-curve is essentially $1.$

ht = density(v)$y;  x = density(v)$x
dx = diff(density(v)$x); summary(dx)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3215  0.3215  0.3215  0.3215  0.3215  0.3215 
wd = mean(dx) # only if all 'dx's equal
sum(ht*wd)
[1] 1.00097

Sample, population and KDE means: The sample mean $(49.68)$ is computed from the data, the population mean $(50)$ by formula from the distribution parameters, and the KDE-mean $(49.73)$ by numerical integration using the common width wd, heights ht, and $x$-values from the KDE.

mean(v); 5/.1
[1] 49.68172
[1] 50

 sum(x*ht*wd)
[1] 49.72998

Population and KDE-modes: By formula the population mode is $4/.1 = 40.$ The modal interval of the histogram is $(40,50].$ And we find the maximum of the KDE curve to be about $43.96.$ (With a different seed, the sample might have produced a KDE-mode that does not fall in the modal histogram interval.)

4/.1
[1] 40
mean(x[ht==max(ht)])
[1] 43.96074

Sample, population, and KDE-medians: The sample median $(46.94)$ is found directly from the sample with median, the population median $(46.71)$ by numerical integration with qgamma, and the KDE-median $(47.50)$ by numerical integration (up to the median) from the KDE output. All three values are near $47.$

median(v);  qgamma(.5, 5, .1)
[1] 46.93721
[1] 46.70909
kde.med = min(x[cumsum(ht*wd)>=.5]);  kde.med
[1] 47.49778

Empirical, population, and KDE CDFs: The ECDF of a sample of size $n$ has jumps of size $1/n$ at sorted sample values (with jumps of $k/n$ in discrete or rounded data $k$ observations tied at the same value).

The top figure below compares the ECDF with the population CDF curve; solid lines show the location of the population median, and dotted lines the location of the sample median. Similarly the bottom figure compares the KDE-CDF and the population CDF (and the respective medians). In case it is of interest, my R code for the figure is shown following the figure.

par(mfrow=c(2,1))  # allow for two graphs per panel
 hdr1 = "ECDF of Sample (Dots) and CDF of GAMMA(5,.1)"
 plot(ecdf(v), main=hdr1)
    lines(c(-10,median(v), median(v)), c(.5, .5, 0), lwd=2, lty="dashed")
  curve(pgamma(x, 5, .1), add=T, lwd=3, col="green2")
    lines(c(-10, qgamma(.5, 5, .1), qgamma(.5, 5, .1)), c(.5, .5, 0), col="red")
 hdr2 = "CDF from KDE and CDF of GAMMA(.5, .1)"
 plot(x, kde.cdf, pch=20, main=hdr2)
     lines(c(-10, kde.med, kde.med), c(.5,.5,0), lwd=2, lty="dashed")
   curve(pgamma(x, 5, .1), add=T, lwd=2, col="green2")
     lines(c(-10, qgamma(.5, 5, .1), qgamma(.5, 5, .1)), c(.5, .5, 0), col="red")
par(mfrow=c(1,1))  # return to single figure plotting

Notes: (1) For your KDE, you need to determine whether the horizontal grid points are equally spaced as in this example with the default KDE in R. If not you will have to use variable rectangle widths in numerical integration. (2) I did not restrict the KDE to positive values. It is a good idea to restrict the KDE to match the support of the distribution involved. (3) Notice that the summaries of $x$ and $y$-values given in the first section of this Answer are not directly applicable to estimating population mean, median, and mode. (4) Disclaimer: Much of what I know about the particular details of R's default KDE was learned by doing this investigation. People on this site who know more about the design of the density function in R, may have important comments to add.