April 2019
 
 
by @mw

The fastest density estimation in R

There are many packages dealing with density estimation in R. They offer several advantages over manually coded algorithms, including bandwidth-selection procedures or involving some more complex features of density estimation, like derivative estimation or higher order kernels. Some of them are also coded in native C language, which should speed up the calculations and enhance memory management. Nevertheless, many of these extra features may be often unused in simple applications of density estimations. That leaves the question open: which algorithm is the fastest?

I try to look at the broadest possible set of R packages dealing with density estimation, including
np
,
kdensity
,
sm
and
GenKern
. There are some other packages which I skip here, as I wanted to make sure I estimate the density at a given point, and not across the domain, to make the numbers comparable. (As a side comment, package
GenKern
seems to have some scaling problems depending on the sample size, but given that the scaling factor is a multiple of 10, I include the package in the comparison). To finally compare the results, I write two functions to perform density estimation. The first one utilizes a standard
for
loop, and the second one offers a vectorized solution. In theory, the vectorized program should have quite good performance. To compare the solutions I use the package
microbenchmark
, which executes the codes multiple times in a random order, and compares the aggregate results.

Although there is a wide literature on density estimation, the technical implementation is quite simple. For more information and technicalities behind the method, I refer you to the set of lecture notes of Bruce Hansen, which offer a great introduction to the topic [LINK]. Given sample $\{X_i: i=1,...,n\}$, where $n$ is the sample size, and a regular kernel function $K$, which typically is a probability density function, one can estimate the density at point $x$ as

$$ f↖{∧}(x) = {1}/{n h} ∑↙{i=1}↖n K({x-X_i}/h). $$

Here, the parameter $h$ is the so-called bandwidth, which controls the amount of smoothing of the estimator. It is usually data-driven and converges to 0 as the sample size tends to infinity. For the experiment I take the same bandwidth parameter for all the estimators.

The strategy is simple. I would like to generate a sample of size 10,000 from a standard normal distribution, and estimate the density at point 1.25 using different algorithms. As emphasized above, I will utilize overall 6 different functions:

  • K.dens - a simple function to sum over the realizations of $X_i$
  • K.vdens - a vectorized version of K.dens
  • npudens from the package
    np
  • kde from the package
    kdensity
  • sm.density from the package
    sm
  • KernSec from the package
    GenKern

The plug-in code to reproduce the results is given below.

#load necessary packages
require(ks)
require(np)
require(kdensity)
require(sm)
require(GenKern)
require(microbenchmark)
 
#non-vectorized density estimates
K.dens = function(x, H, p) {
  x = as.matrix(x)
  H = as.matrix(H)
  n = length(x[,1])
  d = length(x[1,])
 
  dens = 0
  for(i in 1:n) {
    dens = dens + exp( -1/2 * t(p-x[i,]) %*% solve(H) %*% (p-x[i,]) )
  }
 
  return( dens / det(sqrt(H)) / (2*pi)^(d/2) / n )
}
 
#vectorized density estimates
K.vdens = function(x, H, p) {
  x = as.matrix(x)
  H = as.matrix(H)
  d = length(x[1,])
 
  if(d == 1) {
    u = (p - x)/sqrt(H[1,1])
    return( mean( dnorm( u ) / sqrt(H[1,1]) ) )
  }
  if(d > 1)  {
    u = sweep(sweep(x,2,p),2,sqrt(diag(H)),"/")
    return( mean( apply(sapply(1:d, FUN = function(i) dnorm( (p[i] - x[,i])/sqrt(H[i,i]) )/sqrt(H[i,i]) ), 1, FUN = prod) ) )
  }
}
 
#simulation setup
vector <- rnorm(10000, 0, 1)
point  <- c(1.25)
h      <- hns(vector)
H      <- h^2
 
#check if the values are the same
K.dens(x = vector, H = H, p = point)
K.vdens(x = vector, H = H, p = point)
npudens(bws = h, tdat = vector, edat = point)\$dens
kde(vector, h = h, eval.points = point)\$estimate
sm.density(vector, h, eval.points = point)\$estimate
KernSec(vector, xbandwidth = h, range.x=point )\$yden
#compare the performance      
ts <- microbenchmark(K.dens(x = vector, H = H, p = point),
               K.vdens(x = vector, H = H, p = point),
               npudens(bws = h, tdat = vector, edat = point)\$dens,
               kde(vector, h = h, eval.points = point)\$estimate,
               sm.density(vector, h, eval.points = point)\$estimate,
               KernSec(vector, xbandwidth = h, range.x=point )\$yden)
 
boxplot(ts)

The codes generates a plot with time performance for each of the selected functions, probably similar to the graph below.

Candlestick chart As it can be observed, the kde function from
kdensity
package is the best but the vectorized version of the simple estimator seems to be nearly as good. Interestingly, the most advanced in my opinion
np
is somewhat slower, but this is a result of an extra functionality and probably many consistency checks and verifications along the lines. It may well be that the performance of these functions differs if $X$ vector has more dimensions, but this is a topic to be investigated in the next post.

Comments

4 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> canadian pharmacy viagra canadian government approved pharmacies <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> canadian online pharmacy http://canadianmpharmacy.com/# <a href="http://bbs.jmtjmt.com/home.php?mod=space&uid=2756640&do=profile&from=space">canada viagra</a> <a href="http://www.722400.net/home.php?mod=space&uid=813377&do=profile&from=space">canadian drug store</a>

yesterday
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> canadian pharmacy viagra canada drugs <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> online pharmacies http://canadianmpharmacy.com/# <a href="http://www.google.com.uy/url?q=http://viagraswithoutdoctor.com/">canadian pharmacies that ship to us</a> <a href="http://www.google.ht/url?q=http://cialisya.com/">top rated canadian pharmacies online</a>

yesterday
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> northwest pharmacy canada canadian pharmacies without an rx <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> buy viagra online usa http://canadianbpharmacy.com/# <a href="http://www.google.vu/url?q=http://canadiannpharmacy.com/">drugs for sale</a> <a href="http://www.google.ae/url?q=http://viagraswithoutdoctor.com/">prescriptions from canada without</a>

2 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canadian pharmacies online canadian drug <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> pharmacie http://canadianbpharmacy.com/# <a href="http://www.google.nl/url?q=http://cialisvie.com/">canadian pharmacy</a> <a href="https://www.google.com.ag/url?q=http://viagraswithoutdoctor.com/">canadian pharmacies</a>

yesterday
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> pharmacy uk canadian rx <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> online pharmacies india http://canadianbpharmacy.com/# <a href="http://www.google.td/url?q=http://viagraswithoutdoctor.com/">canadian cialis</a> <a href="http://www.google.so/url?q=http://canadianbpharmacy.com/">canadian pharmaceuticals</a>

2 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> pharmacy onesource drugs for sale <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> canadian pharmacy no prescription http://canadianbpharmacy.com/# <a href="http://www.google.mw/url?q=http://canadianmpharmacy.com/">canada drugs</a> <a href="http://www.google.com.sb/url?q=http://viagraawithoutdoctor.com/">canadian prescription drugstore</a>

23 hrs ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> pharmacy canada online prescriptions online pharmacy <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> canada online pharmacy http://canadianbpharmacy.com/# <a href="http://www.google.com.sv/url?q=http://cialisya.com/">online pharmacies in usa</a> <a href="http://www.google.nu/url?q=http://canadianbpharmacy.com/">canada pharmaceuticals online</a>

2 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> buy viagra usa international pharmacy <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> canadian pharmacy cialis http://canadianbpharmacy.com/# <a href="http://www.google.cc/url?q=http://viagraswithoutdoctor.com/">canadian pharmacies that ship to us</a> <a href="http://www.google.ae/url?q=http://canadianbpharmacy.com/">top rated canadian pharmacies online</a>

19 hrs ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> online pharmacy canadian cialis <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> buy viagra now http://canadianmpharmacy.com/# <a href="http://www.google.dj/url?q=http://cialisya.com/">Northwest Pharmacy</a> <a href="http://www.google.mg/url?q=http://canadiannpharmacy.com/">canada drug pharmacy</a>

2 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> canada pharmacy online canadian pharmacies shipping to usa <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> canadian medications online http://canadianmpharmacy.com/# <a href="http://www.google.co.ug/url?q=http://viagraawithoutdoctor.com/">canada online pharmacies</a> <a href="https://www.google.com/url?q=http://canadianbpharmacy.com/">canadian pharmacies online</a>

15 hrs ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canadian pharmacy meds buy viagra usa <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> buy viagra now http://canadianbpharmacy.com/# <a href="http://www.google.ge/url?q=http://canadianbpharmacy.com/">pharmacy online</a> <a href="http://www.google.td/url?q=http://canadiannpharmacy.com/">prescriptions online</a>

2 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> trust pharmacy canada pharmacy uk <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> northwestpharmacy http://canadianmpharmacy.com/# <a href="http://www.google.com.lb/url?q=http://viagraswithoutdoctor.com/">pharmacy online</a> <a href="http://www.google.gm/url?q=http://viagraawithoutdoctor.com/">prescriptions online</a>

11 hrs ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> pharmacy near me global pharmacy canada <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> canadian pharmacies-24h http://canadianbpharmacy.com/# <a href="http://www.google.co.jp/url?q=http://cialisvie.com/">canadian pharmacy online</a> <a href="http://www.google.cat/url?q=http://canadiannpharmacy.com/">canada pharmacy online</a>

2 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> canada drug canadian pharmacy king <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> pharmacy canada http://canadianmpharmacy.com/# <a href="http://www.google.com.sb/url?q=http://cialisvie.com/">pharmacy onesource</a> <a href="http://www.google.co.kr/url?q=http://cialisya.com/">canadianpharmacyusa24h is it legal</a>

7 hrs ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canada drug pharmacy online pharmacies legitimate <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> online pharmacy canada http://canadianbpharmacy.com/# <a href="http://www.google.is/url?q=http://canadiannpharmacy.com/">online pharmacies legitimate</a> <a href="http://www.google.com.sg/url?q=http://canadiannpharmacy.com/">northwest pharmacy canada</a>

3 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> cialis canadian pharmacy canadian medications online <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> canadian pharmaceuticals online http://canadianmpharmacy.com/# <a href="http://www.google.com.fj/url?q=http://viagraswithoutdoctor.com/">online canadian pharmacy</a> <a href="http://www.google.tk/url?q=http://viagraawithoutdoctor.com/">online canadian pharmacies</a>

3 hrs ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> canadian pharmacies that ship to us pharmacy near me <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> buy viagra now http://canadianmpharmacy.com/# <a href="http://www.g.cn/url?q=http://canadianbpharmacy.com/">pharmacy online</a> <a href="http://www.google.md/url?q=http://cialisya.com/">prescriptions online</a>

4 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canadian pharmacy meds canadian pharmacy world <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> online pharmacies canada http://canadianbpharmacy.com/# <a href="http://www.google.co.kr/url?q=http://cialisya.com/">buy vistagra online safe</a> <a href="http://www.google.cc/url?q=http://canadianbpharmacy.com/">online pharmacies of canada</a>

3 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> buy vistagra online safe canada drugs <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> trust pharmacy canada http://canadianbpharmacy.com/# <a href="http://www.google.it/url?q=http://canadianbpharmacy.com/">online drug store</a> <a href="http://www.google.fi/url?q=http://canadianmpharmacy.com/">canadian pharmacy meds</a>

3 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> international pharmacy canada pharmaceuticals online <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> cialis from canada http://canadianmpharmacy.com/# <a href="http://www.google.ch/url?q=http://cialisvie.com/">international pharmacy</a> <a href="http://www.google.se/url?q=http://viagraswithoutdoctor.com/">pharmacies</a>

3 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canada drugs online canadian pharmacy king <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> cialis from canada http://canadianbpharmacy.com/# <a href="http://www.google.com.ph/url?q=http://viagraswithoutdoctor.com/">online canadian pharmacy</a> <a href="http://www.images.google.cn/url?q=http://cialisya.com/">online canadian pharmacies</a>

3 days ago
 
<a href="http://canadianmpharmacy.com/#">canadianmpharmacy.com</a> drugs for sale pharmeasy <a href=http://canadianmpharmacy.com/#>http://canadianmpharmacy.com/</a> pharmacy canada http://canadianmpharmacy.com/# <a href="http://www.google.cm/url?q=http://canadianmpharmacy.com/">northwest pharmacies online</a> <a href="http://www.google.co.ve/url?q=http://cialisya.com/">canadian rx</a>

3 days ago
 
<a href="http://canadianbpharmacy.com/#">canadianbpharmacy.com</a> canada online pharmacy canadian pharmacy no prescription <a href=http://canadianbpharmacy.com/#>http://canadianbpharmacy.com/</a> canadian prescription drugstore http://canadianbpharmacy.com/# <a href="http://www.google.gp/url?q=http://cialisvie.com/">online prescription</a> <a href="http://www.google.gy/url?q=http://canadianbpharmacy.com/">canadian pharcharmy</a>

Leave your comment



M. Wolski
Marcin Wolski, PhD
Economist
European Investment Bank
E-mail: M.Wolski (at) eib.org
Phone: +352 43 79 88708

View my LinkedIn profile View my profile
View my IDEAS/RePEc profile  IDEAS/RePEc