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.
x-grid
point where the KDE is highest to estimate the population mode. $\endgroup$ – BruceET Apr 18 '19 at 18:18