Nothing
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,con1
double precision z(nz),x(nx),xlam(nx),p(nx),qrange,pi,four
double precision con1,sum,sqsum,xsd,a,fifth,hinv,half,g,ginv
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
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.