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