C Output from Public domain Ratfor, version 1.01 subroutine sakj(x,z,p,iker,dens,psi,score,nx,nz,h,alpha,kappa, * xlam) C univariate kernel density-score estimator C the algorithm is basically from Silverman as adapted for Portnoy and Koenker C Annals paper on adaptive L-estimation in regression. C x--pts used for centers of kernel assumed to be sorted!!!! C z--pts at which density is calculated C p--probability associated with x's C dens--f(z), the densities at z C psi--f'(z)/f(z) the score at z C score--(log(f(z)))'', the J fn at z C nx--number of pts in x C nz--number of pts in z C iker--kernel C 0=gaussian C 1=cauchy C h--initial window size (overall)--choose zero for default C kappa--constant determining initial (default) window width C xlam--Silverman's lambda, window adjustment for each x double precision dens(nz),score(nz),psi(nz),g,h,kappa double precision z(nz),x(nx),xlam(nx),p(nx),qrange,pi double precision con1,sum,sqsum,xsd,a,fifth,hinv,half,ginv double precision xn,xker,dxker,ddxker,fact,xponen,alpha,glog,zero, * one,two,four 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 C call srtad(x,1,nx) #port sort routine now done in S interface. if(iker.eq.0) then con1= one/sqrt(2.0*pi) else if(iker.eq.1) then con1= one/pi endif C if no h is provided, calculate a default if(h.le.0.) then sum=0. sqsum=0. do 23006 i=1,nx sqsum=sqsum+x(i)*x(i)*p(i) sum=sum+x(i)*p(i) 23006 continue xsd=dsqrt(sqsum-sum*sum) c compute 'qrange' := IQR (x[i]) sum=zero c first, qrange = Q_1 [ = quantile(*, 0.25) ] do i=1,nx sum=sum+p(i) if(sum .ge. .25) then qrange = x(i) goto 23010 endif enddo 23010 continue sum=one do i=nx,1,-1 sum=sum-p(i) if(sum .le. .75) then qrange = x(i) - qrange goto 23015 endif enddo 23015 continue a=min(xsd,qrange/1.34) h=kappa*a/(xn**fifth) C see Silverman p 48 endif hinv=one/h C Stage one: compute pilot estimate of density do 23018 j=1,nx xker=0. if(iker.eq.0) then do 23022 i=1,nx xponen=(x(j)-x(i))*hinv xponen=half*xponen**2 xker=xker+p(i)*exp(-xponen)*hinv 23022 continue else if(iker.eq.1) then do 23026 i=1,nx xponen=(x(j)-x(i))*hinv xker=xker+p(i)*hinv/(1+xponen**2) 23026 continue endif xlam(j)=con1*xker 23018 continue C Stage two: Automatic window widths (Silverman p101) glog=zero do 23028 i=1,nx glog=glog+p(i)*log(xlam(i)) 23028 continue g=exp(glog) ginv=one/g do 23030 i=1,nx xlam(i)=hinv/((xlam(i)*ginv)**(-alpha)) C notice xlam no longer its own self at this pt! xlam is 1/(h*lambda(i)) C substitution of * for / thus achieved speeds things up C Stage two: new density-score estimates 23030 continue do 23032 j=1,nz xker=zero dxker=zero ddxker=zero if(iker.eq.0) then C gaussian kernel do 23036 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 23036 continue else if(iker.eq.1) then C cauchy kernel do 23040 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) 23040 continue endif dens(j)=con1*xker psi(j)=-(dxker/xker) score(j)=(dxker/xker)**2-ddxker/xker 23032 continue return end