[R] kernel density in 2D

[John Fieberg]

Has anyone written code in R to do kernel density estimation in 2-D with
sphered data?  In other words:

1.  Transform the data to have unit covariance matrix.
2.  Use one of the 2-D kernel estimators in R (e.g., kde2d, bkde2D,
sm.density...) to obtain fhat(x).
3.  Transform back to the original scale.

I have data for which the two dimensions are highly correlated and my
thinking was that this approach might be worth exploring.  I attempted
to write a function to do this (below), but think there must be
something wrong w/ the code I wrote - as things look fine until I  do
the back-transformation.  Any help would be greatly appreciated!

sphere.k<-function(x){
# Function applies the bkde2D kernel density function to the sphered
data #    and then transforms back
#  Inputs:  x data matrix with 2 columns
  
  x<-as.matrix(x)
 
#  Sphere data using cholesky decomposition
#    bn = S^(-1/2)
  bn<-chol(ginv(cov(x)))

#  mu = mean of data
mu<-apply(x,2,mean)

# Transform
  newdata<-(x-outer(rep(1,nrow(x)),mu))%*%t(bn)

# Estimate bandwidth using eq 4.14 in Silverman (1986)
  hsp<-0.96*(nrow(x))^(-1/6)

 est1<-bkde2D(newdata, bandwidth=c(hsp,hsp), gridsize=c(50,50))

# Translate grid back
  ynew<-cbind(est1$x1, est1$x2)
  temp2<-ynew%*%ginv(t(bn))+outer(rep(1,nrow(ynew)),apply(x,2,mean))

# re-order the data for producing contour plots (if necessary)
  tempx<-order(ynew[,1])
  tempy<-order(ynew[,2])
  est1$x1<-temp2[tempx,1]
  est1$x2<-temp2[tempy,2]
  est1$fhat<-est1$fhat[tempx, tempy]*det(bn)
  return(est1)
}