## radially symmetric kernel (Gussian kernel)
## data, x  - where to evaluate kernel function
## b - bandwidth matrix
RadSym <- function(data, x, b) {
  u <- t(solve(b, t(data) - x))
  exp(-rowSums(u^2)/2) / (2*pi)^(ncol(u)/2)
}

## multivariate extension of Scott's bandwidth rule
Scott <- function(data)
  t(chol(cov(data))) * nrow(data) ^ (-1/(ncol(data)+4))

## compute KDE at x given data
mvkde <- function(x, data, bandwidth=Scott, kernel=RadSym) {
  # bandwidth may be a function or matrix
  if(is.function(bandwidth))
    bandwidth <- bandwidth(data)
  mean(kernel(data, x, bandwidth))
}

## compute KDE at (matrix) x given data
smvkde <- function(x, ...)
  apply(x, 1, mvkde, ...)

## scale data to between range
scale_rng <- function(x, rng=c(0,1)) {
  rngx <- range(x, na.rm=TRUE)
  z <- (x-rngx[1])/diff(rngx)
  z*diff(rng) + rng[1]
}
## Univariate example with 'airquality' data
## compute univariate KDE and plot 
data("airquality")
aq <- subset(airquality, !is.na(Solar.R), select="Solar.R")
hist(airquality$Solar.R, xlab='Solar.R', breaks=20, freq=F,
     main='Histogram of Solar.R')
for(mu in unique(aq$Solar.R)) {
  curve(scale_rng(dnorm(x, mean=mu, sd=5.7), rng=c(0,0.0005)),
        from=0, to=350, add=TRUE, n=300)
}

## compute KDE on a grid of Solar.R values
dens.Solar.R <- seq(min(aq$Solar.R),max(aq$Solar.R),length.out=100)
dens.grid <- expand.grid(Solar.R=dens.Solar.R)
dens.vals <- smvkde(dens.grid, data=aq)
lines(dens.Solar.R, scale_rng(dens.vals, c(0,0.006)))

## Bivariate example with 'airquality' data
## compute bivariate KDE and plot contours
data("airquality")
aq <- subset(airquality, !is.na(Ozone) & !is.na(Solar.R),
             select=c("Ozone", "Solar.R"))

## compute KDE on a grid of Ozone and Solar.R values
dens.Ozone <- seq(min(aq$Ozone),max(aq$Ozone),length.out=100)
dens.Solar.R <- seq(min(aq$Solar.R),max(aq$Solar.R),length.out=100)
dens.grid <- expand.grid(Ozone=dens.Ozone, Solar.R=dens.Solar.R)
dens.vals <- smvkde(dens.grid, data=aq)

## arrange density values into matrix for easy plotting
dens.mtrx <- matrix(dens.vals, 100, 100)
contour(x=dens.Ozone, y=dens.Solar.R, z=dens.mtrx,
        xlab="Ozone", ylab="Solar.R")
points(aq$Ozone, aq$Solar.R, pch=20)