# Estimating the population median from a kernel density estimator

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?

• Use the sample median of the data you used to get the KDE to estimate the median of the population. Use the sample mean of the data to estimate the population mean. Use the x-grid point where the KDE is highest to estimate the population mode. – BruceET Apr 18 '19 at 18:18
• @BruceET Thanks. I suppose those are unbiased estimates, but I'm thinking of a situation where I want to or need to derive the estimates for the mean and median from the kernel density estimate. – colorlace Apr 18 '19 at 19:11
• Do you know what software generated the your KDE? While playing around with the default KDE in R, I have learned some things that may be helpful. If your KDEs are from R, I'll give code for some productive experiments. If not, I'll try to discuss findings in generic terms. Or if you're done with this topic, please let me know that. – BruceET Apr 19 '19 at 4:23
• R code for your productive experiments would be much appreciated! – colorlace Apr 19 '19 at 16:43
• OK. Will resurrect that and hope to post in 4-5 hrs. – BruceET Apr 19 '19 at 18:52

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")

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.