wilcox_R: What the package does (short line)

akerdmul <-
function(x,pts=NA,hval=NA,aval=.5,fr=.8,pr=FALSE,plotit=TRUE,theta=50,
phi=25,expand=.5,scale=FALSE,xlab="X",ylab="Y",zlab="",ticktype="simple"){
#
# Compute adaptive kernel density estimate
# for multivariate data
# (See Silverman, 1986)
#
#  Use expected frequency as initial estimate of the density
#
# hval is the span used by the kernel density estimator
# fr is the span used by the expected frequency curve
# pr=T, returns density estimates at pts
# ticktype="detailed" will create ticks as done in two-dimensional plot
#
library(MASS)
library(akima)
if(is.na(pts[1]))pts<-x
if(ncol(x)!=ncol(pts))stop("Number of columns for x and pts do not match")
if(!is.matrix(x))stop("Data should be stored in a matrix")
fhat <- rdplot(x,pyhat=TRUE,plotit=FALSE,fr=fr)
n<-nrow(x)
d<-ncol(x)
pi<-gamma(.5)^2
cd<-c(2,pi)
if(d==2)A<-1.77
if(d==3)A<-2.78
if(d>2){
for(j in 3:d)cd[j]<-2*pi*cd[j-2]/n  # p. 76
}
if(d>3)A<-(8*d*(d+2)*(d+4)*(2*sqrt(pi))^d)/((2*d+1)*cd[d])  # p. 87
if(is.na(hval))hval<-A*(1/n)^(1/(d+4))  # Silverman, p. 86
svec<-NA
for(j in 1:d){
sig<-sqrt(var(x[,j]))
temp<-idealf(x[,j])
iqr<-(temp$qu-temp$ql)/1.34
A<-min(c(sig,iqr))
x[,j]<-x[,j]/A
svec[j]<-A
}
hval<-hval*sqrt(mean(svec^2)) # Silverman, p. 87
# Now do adaptive; see Silverman, 1986, p. 101
gm<-exp(mean(log(fhat[fhat>0])))
alam<-(fhat/gm)^(0-aval)
dhat<-NA
nn<-nrow(pts)
for(j in 1:nn){
temp1<-t(t(x)-pts[j,])/(hval*alam)
temp1<-temp1^2
temp1<-apply(temp1,1,FUN="sum")
temp<-.5*(d+2)*(1-temp1)/cd[d]
epan<-ifelse(temp1<1,temp,0) # Epanechnikov kernel, p. 76
dhat[j]<-mean(epan/(alam*hval)^d)
}
if(plotit && d==2){
fitr<-dhat
iout<-c(1:length(fitr))
nm1<-length(fitr)-1
for(i in 1:nm1){
ip1<-i+1
for(k in ip1:length(fitr))if(sum(x[i,]==x[k,])==2)iout[k]<-0
}
fitr<-fitr[iout>=1]
mkeep<-x[iout>=1,]
fit<-interp(mkeep[,1],mkeep[,2],fitr)
persp(fit,theta=theta,phi=phi,xlab=xlab,ylab=ylab,zlab=zlab,expand=expand,
scale=scale,ticktype=ticktype)
}
m<-"Done"
if(pr)m<-dhat
m
}