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?


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: enter image description here

| cite | improve this answer | |
  • $\begingroup$ +1 This is a general answer because most kernels are developed to have mathematically tractable expressions for their PDFs and CDFs. $\endgroup$ – whuber Jan 21 '16 at 15:18

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:


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  

Sure enough, it works:

enter image description here

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

| cite | improve this answer | |
  • 3
    $\begingroup$ I believe you can use rbeta to speed this up considerably. The parameters are $3$ and $3$. Just double the result and subtract $1$. Look at this: hist(rbeta(1e6, 3, 3)*2 - 1, freq=FALSE) (about $0.2$ seconds). $\endgroup$ – whuber Jan 21 '16 at 17:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.