kdem = function(x, # points to evaluate X, # observed data points s #bandwidth vector ){ ## multivariate kernel density estimator ## Make a list to store difference matrices diff.ls=list() for(i.dim in 1:ncol(X)){ diff.ls[[i.dim]] = matrix(x[,i.dim],nrow(x),nrow(X),byrow=FALSE) - matrix(X[,i.dim],nrow(x),nrow(X),byrow=TRUE) } diffmat=do.call(cbind,lapply(diff.ls,as.numeric)) Kvecx = dmvnorm(diffmat,sigma=diag(s)) Kmatx = matrix(Kvecx, nrow(x), nrow(X) ,byrow=FALSE) return( Kmatx%*%rep(1/nrow(X),nrow(X) ) ) } kdem.cv = function(X, #observed data points s #bandwidth ){ ## Calculate log likelihood of kernel # First we'll calculate the kernels everywhere diff.ls=list() for(i.dim in 1:ncol(X)){ diff.ls[[i.dim]] = matrix(X[,i.dim],nrow(X),nrow(X),byrow=FALSE) - matrix(X[,i.dim],nrow(X),nrow(X),byrow=TRUE) } diffmat=do.call(cbind,lapply(diff.ls,as.numeric)) Kvecx = dmvnorm(diffmat,sigma=diag(s)) Kmatx = matrix(Kvecx, nrow(X), nrow(X) ,byrow=FALSE) # Deleting the diagonal will remove X_i from influencing # it's own density Kmatx = Kmatx - diag(diag(Kmatx)) # Now get f^(-i)(X_i) (remember to use n-1) fvec = Kmatx%*%rep(1/(nrow(X)-1),nrow(X)) # and return sum of logs return( sum(log(fvec)) ) } kdem.cv_wrapper=function(s, X){return(-kdem.cv(X,s))} ## for running with optim ##################### Example ######################## #generate "observed" data, grid X=matrix(c(1,1,1.1, 1.1, 1.2,1), ncol=2) x=expand.grid(seq(0,2,length=20), seq(0,2,length=10)) ## CV, then evaluate opt=optim(c(.2,.2), fn=wrapper, X=X) out=kdem(x=x, X=X, s=opt$par) ## plot color.gradient <- function(x, colors=c("red","yellow","green"), colsteps=100) { return( colorRampPalette(colors) (colsteps) [ findInterval(x, seq(min(x),max(x), length.out=colsteps)) ] ) } plot(x=x[,1], y=x[,2],pch=15,cex=4, col=color.gradient(out)) points(x=X[,1], y=X[,2], pch=19)