Why does this simple weighted quantile differ from Hmisc::wtd.quantile? Which method is to be preferred? It just struck me today that we should be able to use the weights method of stats::density.default to roll our own simple weighted quantile function (I'll call it my_wtd_q) with base functionality:
my_wtd_q = function(x, w, prob, n = 4096) {
  with(
    density(x, weights = w/sum(w), n = n), 
    x[which.max(cumsum(y*(x[2L] - x[1L])) >= prob)]
  )
}

The idea being that if we have a weighted PDF, we can construct the weighted CDF, and thereby the weighted quantiles.
However, the implied results look pretty different from those of what I took as the "canonical" function for weighted quantiles in R, namely Hmisc::wtd.quantile:
set.seed(3049)
p = seq(0, 1, length.out = 100)
png('~/Desktop/wtd_quantile.png', width = 1920, height = 1920, res = 100)
par(mfrow = c(2, 2), mar = c(0, 0, 0, 0), oma = c(5.1, 4.1, 4.1, 2.1))
for (nn in 10^(2:5)) {
  x = rnorm(nn)
  w = rchisq(nn, df = ceiling(abs(x^3)))
    
  add_x = nn %in% c(1000, 10000)
  add_y = nn %in% c(10, 1000)

  my_wtd_y = sapply(p, my_wtd_q, x = x, w = w)
  hmisc_y = sapply(p, Hmisc::wtd.quantile, x = x, w = w)
  unwtd_y = sapply(p, quantile, x = x)

  matplot(
    p, cbind(my_wtd_y, hmisc_y, unwtd_y),
    xaxt = if (!add_x) 'n', xlab = '',
    yaxt = if (!add_y) 'n', ylab = '',
    type = 'l', lty = 1L, lwd = 2L, las = 1L, main = ''
  )
  title(line = -1, sprintf('n = %s', prettyNum(nn, big.mark = ',')))
  if (add_x) mtext(side = 1L, 'Quantile', line = 3)
  if (add_y) mtext(side = 2L, 'Inverse CDF', line = 3)
  legend(
    'topleft', col = 1:3, lwd = 2L,
    legend = c('Simple Weighted', 'Hmisc::wtd.quantile', 'Unweighted')
  )
}
title(
  'Comparison of Weighted Quantile Methods\nVarious Sample Sizes', 
  outer = TRUE
)
dev.off()


It certainly seems like the two approaches are asymptotically equivalent, but I'm curious about the source of the divergence at smaller sample sizes (Perhaps my method is making stronger smoothness assumptions? Or they're using different implied assumptions about interpolation?)
Can we say for sure one method is "better" than the other?
 A: You could also try:
quantreg::rq(x~1, weights = w, tau = p)
use of density(...) is contra-indicated since it presumes smoothness and
because it inevitably inherits the tail behavior of the kernel in the extremes
as illustrated clearly in your plot for n = 100.
A: The difference between your method and wtd.quantile (and also survey::svyquantile) is the use of smoothing to make a PDF before constructing the CDF
There's no need to smooth: we know that the unsmoothed CDF converges uniformly to the true CDF (under no assumptions on the distribution and weak assumptions on the weights).
On the other hand, it is known in the unweighted setting that a bootstrap for quantiles works better with a bit of kernel smoothing if the true density is differentiable (discusssed here). So it makes sense that you might want to smooth for weighted quantiles.
On the other other hand, you only get the (asymptotic) improvement with the right amount of smoothing. The linked paper notes "The optimal bandwidth is of the same order, though with a different constant, as that which minimises the mean integrated squared error of [the density]"
And, of course, if the distribution doesn't have a density, smoothing has got to make the estimation less accurate.
