title: Difference in 2D KDE produced using kde2d (R) and ksdensity2d (Matlab)
body: <p>While trying to port some code from Matlab to R I have run into a problem. The gist of the code is to produce a 2D kernel density estimate and then do some simple calculations using the estimate. In Matlab the KDE calculation was done using the function <a href="https://github.com/probml/pmtk3/blob/master/matlabTools/stats/ksdensity2d.m" rel="noreferrer">ksdensity2d.m</a>. In R the KDE calculation is done with kde2d from the MASS package. So lets say I want to calculate the KDE and just add the values (this is not what i intend to do, but it serves this purpose). In R this can be done by</p>
  
  <pre><code>    library(MASS)
      set.seed(1009)
      x &lt;- sample(seq(1000, 2000), 100, replace=TRUE)
      y &lt;- sample(seq(-12, 12), 100, replace=TRUE)
      kk &lt;- kde2d(x, y, h=c(30, 1.5), n=100, lims=c(1000, 2000, -12, 12))
      sum(kk$z)
  </code></pre>
  
  <p>which gives the answer 0.3932732. When using ksdensity2d in Matlab using the same exact data and conditions the answer is 0.3768. From looking at the code for kde2d I noticed that the bandwidth is divided by 4</p>
  
  <pre><code>    kde2d &lt;- function (x, y, h, n = 25, lims = c(range(x), range(y))) 
      {
      nx &lt;- length(x)
      if (length(y) != nx) 
       stop("data vectors must be the same length")
      if (any(!is.finite(x)) || any(!is.finite(y))) 
       stop("missing or infinite values in the data are not allowed")
      if (any(!is.finite(lims))) 
       stop("only finite values are allowed in 'lims'")
      n &lt;- rep(n, length.out = 2L)
      gx &lt;- seq.int(lims[1L], lims[2L], length.out = n[1L])
      gy &lt;- seq.int(lims[3L], lims[4L], length.out = n[2L])
      h &lt;- if (missing(h)) 
      c(bandwidth.nrd(x), bandwidth.nrd(y))
      else rep(h, length.out = 2L)
      if (any(h &lt;= 0)) 
       stop("bandwidths must be strictly positive")
      h &lt;- h/4
      ax &lt;- outer(gx, x, "-")/h[1L]
      ay &lt;- outer(gy, y, "-")/h[2L]
      z &lt;- tcrossprod(matrix(dnorm(ax), , nx), matrix(dnorm(ay), 
       , nx))/(nx * h[1L] * h[2L])
      list(x = gx, y = gy, z = z)
      }
  </code></pre>
  
  <p>A simple check to see if the difference in bandwidth is the reason for the difference in the results is then</p>
  
  <pre><code>    kk &lt;- kde2d(x, y, h=c(30, 1.5)*4, n=100, lims=c(1000, 2000, -12, 12))
      sum(kk$z)
  </code></pre>
  
  <p>which gives 0.3768013 (which is the same as the Matlab answer).</p>
  
  <p>So my question is then: Why does kde2d divide the bandwidth by four? (Or why doesn't ksdensity2d?)</p>