# [R] kernel density in 2D

John Fieberg John.Fieberg at dnr.state.mn.us
Fri Aug 1 22:39:07 CEST 2003

```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)
}

```