SUBROUTINE HHDD(X,N,XIQR,K,HHD) C----------------------------------------------------------------------- C VERSION: January 1992 C C JOACHIM ENGEL C C PURPOSE: C C GENERAL SUBROUTINE FOR AUTOMATIC BANDWIDTH SELECTION WITH C HEIDELBERG METHOD for estimation of k-th derivative C C COMPUTATION OF OPTIMAL BANDWIDTH USING THE ASYMPTOTIC FORMULA C WITH ITERATIVE REFINEMENT, FOR A Normal Kernel C ITERATION METHOD, iterations C Bandwidth inflation for estimating functional: C A = CO*H*N**2/35 or correspondigly for derivatives C (A BANDWIDTH FOR ESTIMATING FUNCTIONAL) C C PARAMETERS C C INPUT X(N) DATA (MUST BE SORTED) C INPUT N LENGTH OF X C INPUT XIQR INTERQUARTILE RANGE C INPUT K ORDER OF DERIVATIVE (0, 1 or 2) C C OUTPUT HHD ESTIMATED OPTIMAL BANDWIDTH C C C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(N) pi=3.141592645d0 rt2pi = dsqrt(2.d0*pi) xn=dble(n) m=k+2 l=m+2 C- C- C------- COMPUTE CONSTANTS FOR LATER USE c xiqr=x(nint(.75*xn))-x(nint(.25*xn)+1) c ITER=3*m-l+1 iter=3 EXT=1./(2.*m+1) EXINF=2./(2*m+1)/(l+m+1.) c c if to estimate density itself (k=0) c if (k.eq.0) then c c- **estimate inflation constant c h2=(.920*xiqr)/(xn**(1.d0/7.d0)) h3=(.912*xiqr)/(xn**(1.d0/9.d0)) s2 = 0.d0 s3 = 0.d0 do 20 i = 1, n-1 do 30 j = i+1, n d2 = ((x(i) - x(j))/h2)**2 d3 = ((x(i) - x(j))/h3)**2 e2 = dexp(-0.5d0*d2) e3 = dexp(-0.5d0*d3) s2 = s2 + (d2**2 - 6.d0*d2 + 3.d0)*e2 s3 = s3 + (d3**3 - 15.d0*d3**2 + 45.d0*d3 - 15.d0)*e3 30 continue 20 continue rhat2 = (2.d0*s2)/((xn**2)*(h2**5)*rt2pi) rhat2 = rhat2 + 3.d0/(rt2pi*xn*(h2**5)) rhat3 = (-2.d0*s3)/((xn**2)*(h3**7)*rt2pi) rhat3 = rhat3 + 15.d0/(rt2pi*xn*(h3**7)) chat = 1.459058157d0*(rhat2**(1.d0/5.d0))/(rhat3**(1.d0/7.d0)) c write(*,*) 'rhat2=',rhat2, ' rhat3=',rhat3,' chat=',chat c c CO1=chat*(XN**EXINF) C- C------- COMPUTE CONSTANT OF ASYMPTOTIC FORMULA CONST=1.d0/(2.d0*dsqrt(pi)) a=1.132795764/rhat3**(1.d0/7.d0)*xn**(-1.d0/7.d0) C- C------ LOOP OVER ITERATIONS DO 100 IT=1,ITER C- C------- ESTIMATE Functional c s=0.d0 do 40 i = 1, n-1 do 50 j = i+1, n d2 = ((x(i) - x(j))/a)**2 e2 = dexp(-0.5d0*d2) s = s + (d2**2 - 6.d0*d2 + 3.d0)*e2 50 continue 40 continue r2 = (2.d0*s)/((xn**2)*(a**5)*rt2pi) r2 = r2 + 3.d0/(rt2pi*xn*(a**5)) C- C- C- C------- ESTIMATE BANDWIDTH BY ASYMPTOTIC FORMULA HHD=(CONST/(R2*XN))**(.2d0) A=HHD*CO1 100 continue c write(*,*) 'Final rhat2=',r2 endif c c if to estimate first derivative (k=1) c if (k.eq.1) then c c- **estimate inflation constant c c h3=(.912*xiqr)/(xn**(1.d0/9.d0)) h4=(.914*xiqr)/(xn**(1.d0/11.d0)) s3 = 0.d0 s4 = 0.d0 do 21 i = 1, n-1 do 31 j = i+1, n d3 = ((x(i) - x(j))/h3)**2 d4 = ((x(i) - x(j))/h4)**2 e3 = dexp(-0.5d0*d3) e4 = dexp(-0.5d0*d4) s3 = s3 + (d3**3 - 15.d0*d3**2 + 45.d0*d3 - 15.d0)*e3 s4=s4+ (d4**4-28.d0*d4**3+210.d0*d4**2-420.d0*d4 +105.d0)*e4 31 continue 21 continue rhat3 = (-2.d0*s3)/((xn**2)*(h3**7)*rt2pi) rhat3 = rhat3 + 15.d0/(rt2pi*xn*(h3**7)) rhat4 = (2.d0*s4)/((xn**2)*(h4**9)*rt2pi) rhat4 = rhat4 + 105.d0/(rt2pi*xn*(h4**9)) chat = 1.489839163d0*(rhat3**(1.d0/7.d0))/(rhat4**(1.d0/9.d0)) c c CO1=chat*(XN**EXINF) C- C------- COMPUTE CONSTANT OF ASYMPTOTIC FORMULA CONST=3.d0/4.d0/dsqrt(pi) a=1.317592936/rhat4**(1.d0/9.d0)*xn**(-1.d0/9.d0) C- C------ LOOP OVER ITERATIONS DO 101 IT=1,ITER C- C------- ESTIMATE Functional c s=0.d0 do 41 i = 1, n-1 do 51 j = i+1, n d3 = ((x(i) - x(j))/a)**2 e3 = dexp(-0.5d0*d3) s = s + (d3**3 -15.d0*d3**2 + 45.d0*d3 -15.d0)*e3 51 continue 41 continue r3 = (-2.d0*s)/((xn**2)*(a**7)*rt2pi) r3 = r3 + 15.d0/(rt2pi*xn*(a**7)) C- C- C- C------- ESTIMATE BANDWIDTH BY ASYMPTOTIC FORMULA HHD=(CONST/(R3*XN))**(1.d0/7.d0) A=HHD*CO1 101 continue c write(*,*) 'rhat3=',r3 c endif if (k.eq.2) then c c- **estimate inflation constant c h4=(.914*xiqr)/(xn**(1.d0/11.d0)) h5=(.919*xiqr)/(xn**(1.d0/13.d0)) s4 = 0.d0 s5 = 0.d0 do 22 i = 1, n-1 do 32 j = i+1, n d4 = ((x(i) - x(j))/h4)**2 d5 = ((x(i) - x(j))/h5)**2 e4 = dexp(-0.5d0*d4) e5 = dexp(-0.5d0*d5) s4 = s4+ (d4**4 - 28.d0*d4**3 + 210.d0*d4**2 - . 420.d0*d4 +105.d0)*e4 s5 = s5 + (d5**5 - 45.d0*d5**4 + 630.d0*d5**3 - . 3150.d0*d5**2 +4725.d0*d5-945.d0)*e5 32 continue 22 continue rhat4 = (2.d0*s4)/((xn**2)*(h4**9)*rt2pi) rhat4 = rhat4 + 105.d0/(rt2pi*xn*(h4**9)) rhat5 = (-2.d0*s5)/((xn**2)*(h5**11)*rt2pi) rhat5 = rhat5 + 945.d0/(rt2pi*xn*(h5**11)) chat = 1.679304212d0*(rhat4**(1.d0/9.d0))/(rhat5**(1.d0/11.d0)) c write(*,*) 'rhat4=',rhat4, ' rhat5=',rhat5,' chat=',chat c CO1=chat*(XN**EXINF) C- C------- COMPUTE CONSTANT OF ASYMPTOTIC FORMULA CONST=15.d0/8.d0/dsqrt(pi) a=1.495649885/rhat5**(1.d0/11.d0)*xn**(-1.d0/11.d0) C- C------ LOOP OVER ITERATIONS DO 102 IT=1,ITER C- C------- ESTIMATE Functional c s=0.d0 do 42 i = 1, n-1 do 52 j = i+1, n d4 = ((x(i) - x(j))/a)**2 e4 = dexp(-0.5d0*d4) s = s + (d4**4 -28.d0*d4**3 + 210.d0*d4**2 - . 420.d0*d4 + 105.d0)*e4 52 continue 42 continue r4 = (2.d0*s)/((xn**2)*(a**9)*rt2pi) r4 = r4 + 105.d0/(rt2pi*xn*(a**9)) C- C- C- C------- ESTIMATE BANDWIDTH BY ASYMPTOTIC FORMULA HHD=(CONST/(R4*XN))**(1.d0/9.d0) A=HHD*CO1 102 continue c write(*,*) 'rhat4=' ,r4 c endif C- C- RETURN END