akj.r
``````subroutine sakj(x,z,p,iker,dens,psi,score,nx,nz,h,alpha,kappa,xlam)
# univariate kernel density-score estimator
# the algorithm is basically from Silverman as adapted for Portnoy and Koenker
# Annals paper on adaptive L-estimation in regression.

# x--pts used for centers of kernel assumed to be sorted!!!!
# z--pts at which density is calculated
# p--probability associated with x's
# dens--f(z), the densities at z
# psi--f'(z)/f(z) the score at z
# score--(log(f(z)))'', the J fn at z
# nx--number of pts in x
# nz--number of pts in z
# iker--kernel
#          0=gaussian
#	   1=cauchy
# h--initial window size (overall)--choose zero for default
# kappa--constant determining initial (default) window width
# xlam--Silverman's lambda, window adjustment for each x

double precision dens(nz),score(nz),psi(nz),h,kappa
double precision z(nz),x(nx),xlam(nx),p(nx),qrange,pi
double precision con1,sum,sqsum,xsd,a,fifth,hinv,half
double precision xn,xker,dxker,ddxker,fact,xponen,alpha,glog,zero,one,two
parameter( zero = 0.d0)
parameter( one = 1.d0)
parameter( two = 2.d0)
parameter( four = 4.d0)
parameter( half = 0.5d0)
parameter( fifth = 0.2d0)
parameter( pi = 3.141593d0)
xn=nx

# call srtad(x,1,nx) # port sort routine now done in S interface.
if(iker==0) con1=one/sqrt(2.0*pi)
else if(iker==1) con1=one/pi

# if no h is provided, calculate a default
if(h<=0.) {
sum=0.
sqsum=0.
do i=1,nx
{ sqsum=sqsum+x(i)*x(i)*p(i)
sum=sum+x(i)*p(i)
}
xsd=dsqrt(sqsum-sum*sum)
sum=zero
for(i=1;i<nx;i=i+1) {
sum=sum+p(i)
if(sum<.25) next
else { qrange=x(i);break }
}
sum=one
for(i=nx;i>0;i=i-1) {
sum=sum-p(i)
if(sum>.75) next
else { qrange=x(i)-qrange;break }
}
a=min(xsd,qrange/1.34)
h=kappa*a/(xn**fifth) # see Silverman p 48
}
hinv=one/h
# Stage one:  compute pilot estimate of density
do j=1,nx {
xker=0.
if(iker==0) {
do i=1,nx {
xponen=(x(j)-x(i))*hinv
xponen=half*xponen**2
xker=xker+p(i)*exp(-xponen)*hinv
}
}
else if(iker==1) {
do i=1,nx {
xponen=(x(j)-x(i))*hinv
xker=xker+p(i)*hinv/(1+xponen**2)
}
}
xlam(j)=con1*xker
}
# Stage two:  Automatic window widths (Silverman p101)
glog=zero
do i=1,nx
{ glog=glog+p(i)*log(xlam(i))}
g=exp(glog)
ginv=one/g
do i=1,nx {
xlam(i)=hinv/((xlam(i)*ginv)**(-alpha))
}
# notice xlam no longer its own self at this pt! xlam is 1/(h*lambda(i))
# substitution of * for / thus achieved speeds things up

# Stage two:  new density-score estimates
do j=1,nz {
xker=zero
dxker=zero
ddxker=zero
if(iker==0) {  # gaussian kernel
do i=1,nx {
xponen=(z(j)-x(i))*xlam(i)
fact=exp(-half*xponen*xponen)*xlam(i)
xker=xker+p(i)*fact
dxker=dxker-p(i)*fact*xponen*xlam(i)
ddxker=ddxker-
p(i)*fact*(one - xponen**2)*xlam(i)**2
}
}
else if(iker==1) {  # cauchy kernel
do i=1,nx {
xponen=(z(j)-x(i))*xlam(i)
fact=xlam(i)/(one+xponen**2)
xker=xker+p(i)*fact
dxker=dxker-p(i)*two*xponen*fact**2
ddxker=ddxker-
p(i)*two*(fact**2)*(xlam(i)-
four*(xponen**2)*fact)
}
}
dens(j)=con1*xker
psi(j)=-(dxker/xker)
score(j)=(dxker/xker)**2-ddxker/xker
}
return
end
``````