Generating samples from kernels other than Epanechnikov's

Question

Generating samples from kernel densities with Gaussian kernel is very simple, the same with uniform, or triangular kernels. There is also a smart algorithm for generating samples from Epanchenikov kernel, however is it possible to generate samples from other kernels (e.g. biweight, cosine, optcosine) in more efficient way than by using than accept-reject algorithm?

Answer 1

I suggest the inverse cdf method. For the biweight kernel, $$ \begin{align}F(x)&=\int_{-1}^x \frac{15}{16}(1-u^2)^2\text{d}u\\ &= \int_{-1}^x \frac{15}{16}(1-2u^2+u^4)\text{d}u\\ &=\frac{15}{16}\left\{ x+1-2\frac{x^3+1}{3}+\frac{x^5+1}{5}\right\}\\\end{align} $$ which can be inverted numerically, while, for the cosine kernel, $$ \begin{align}F(x)&=\int_{-1}^x \frac{\pi}{4}\cos(\pi u/2)\text{d}u\\ &= \frac{1}{2}\{\sin(\pi x/2)+1\}\end{align}$$ which produces $$F^{-1}(u)=2\arcsin(2u-1)/\pi$$ Note that the biweight cdf looks very similar to the logistic cdf $$G(x)=\dfrac{1}{1+\exp\{-4x^2+3\text{sign}(x)x^2/2\}}$$which helps with the inversion:

Answer 2

My answer is maybe as much a question related to Xi'an's +1 answer.

I took up the idea of numerically inverting the biweight cdf, using this code:

library(MASS)

inverse  <-  function (f, lower = -1, upper = 1) {
  function (y) uniroot((function (x) f(x) - y), lower = lower, upper = upper)[1]
}

biweight_inverse <- inverse(function (x) (abs(x)<=1)*15/16*(x+1-2*(x^3+1)/3+(x^5+1)/5))

reps <- 1e6
biweight_realizations <- matrix(NA,reps)
for (i in 1:reps){
  biweight_realizations[i] <- biweight_inverse(runif(1))$root  
}
truehist(biweight_realizations)

Sure enough, it works:

I however suspect that OP wants to avoid rejection sampling for efficiency considerations. Now, my code is awfully slow in view of the need to perform numerical inversion.

So I wonder whether (a) the inverse cdf method is really a good idea when the quantile function needs to be found numerically or (b) my code is just bad. (I am sure the loop could be avoided, but I doubt this is what slows down the simulation.)