kern_reg.f90

The subroutines glkern.f and lokern.f use an efficient and fast algorithm for automatically adapti

CALL kernel(t, x, n, g1, kord, kord+2, nyg, s, w1(2:,2), m1-1, w1(2:,3))

DO  j=2,m1
  w1(j,3)=w1(j,3)*w1(j,3)
END DO
DO  j=1,m
  y(j)=g2
  IF(tt(j) < s(0)+g1) THEN
    y(j)=g2*(1.0 + ((tt(j)-g1-s(0))/g1)**2)
    y(j)=MIN(y(j), bmax)
  ELSE IF(tt(j) > s(n)-g1) THEN
    y(j)=g2*(1.0 + ((tt(j)-s(n)+g1)/g1)**2)
    y(j)=MIN(y(j), bmax)
  END IF
END DO

CALL kernp(w1(2:,2), w1(2:,3), m1-1, g2, nuev, kordv, nyl, w1(1:,1), wm, m, y)
!-
!-------- 22. FINISH
DO  j=1,m
  xh=bmin**(2*kord+1)*ABS(y(j))*vi/const
  xxh=const*ABS(ban(j))/vi/bmax**(2*kord+1)
  IF(ban(j) < xh) THEN
    ban(j)=bmin
  ELSE
    IF(y(j) < xxh) THEN
      ban(j)=bmax
    ELSE
      ban(j)=(const*ban(j)/y(j)/vi)**ex
    END IF
  END IF
END DO
!-
!-------- 23. COMPUTE SMOOTHED FUNCTION WITH LOCAL PLUG-IN BANDWIDTH
230 y(1:m)=ban(1:m)
CALL kernel(t, x, n, b, nue, kord, nyl, s, tt, m, y)

RETURN
END SUBROUTINE lokern



!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
!C     SEVERAL KERNEL SMOOTHING SUBROUTINES WHICH ARE USED BY GLKERN.F
!C     AND LOKERN.F, VERSION JANUARY 1997
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C     THIS FILE CONTAINS:
!C
!C     SUBROUTINE RESEST(T,X,N,RES,SNR,SIGMA2)
!C                FOR VARIANCE ESTIMATION
!C
!C     SUBROUTINE KERNEL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                DRIVER SUBROUTINE FOR KERNEL REGRESSION ESTIMATION
!C                CALLS FAST OR CONVENTIAL KERNEL ROUTINE
!C
!C     SUBROUTINE KERNP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                DRIVER SUBROUTINE FOR KERNEL REGRESSION ESTIMATION
!C                WITHOUT USE OF BOUNDARY KERNELS
!C---------------------------------------------------------------------
!C     SUBROUTINE KERNFA(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                FAST ALGORITHM FOR KERNEL ESTIMATION
!C     SUBROUTINE DREG(SW,A1,A2,IORD,X,SL,SR,T,B,IFLOP)
!C                USED BY SUBROUTINE KERNFA,KERNFP
!C     SUBROUTINE LREG(SW,A3,IORD,D,DOLD,Q,C)
!C                USED BY SUBROUTINE KERNFA,KERNFP
!C     SUBROUTINE FREG(SW,NUE,KORD,IBOUN,Y,C,ICALL,A)
!C                USED BY SUBROUTINE KERNFA,KERNFP
!C     SUBROUTINE KERNFP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                FAST ALGORITHM FOR KERNP ESTIMATION WITHOUT BOUNDARY
!C---------------------------------------------------------------------
!C     SUBROUTINE KERNCL(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                CONVENTIONAL ALGORITHM FOR KERNEL ESTIMATION
!C     SUBROUTINE SMO(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y)
!C                SINGLE ESTIMATION STEP, USED BY KERNCL
!C     SUBROUTINE KERNCP(T,X,N,B,NUE,KORD,NY,S,TT,M,Y)
!C                CONVENTIONAL ALGORITHM WITHOUT BOUNDARY KERNELS
!C     SUBROUTINE SMOP(S,X,N,TAU,WID,NUE,IORD,IBOUN,IST,S1,C,Y)
!C                SINGLE ESTIMATION STEP, USED BY KERNCP
!C---------------------------------------------------------------------
!C     SUBROUTINE COFFI(NUE,KORD,C)
!C                KERNEL COEFFICIENT OF POLYNOMIAL KERNELS USED BY
!C                KERNCL,KERNCP AND KERNFP
!C     SUBROUTINE COFFB(NUE,KORD,Q,IBOUN,C)
!C                KERNEL COEFFICIENT OF POLYNOMIAL BOUNDARY KERNELS
!C                USED BY KERNFA AND KERNCL
!C---------------------------------------------------------------------
!C     SUBROUTINE COFF(X,N,FA)
!C                SIMPLE SUBROUTINE FOR ARRAY INITIALIZATION
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


SUBROUTINE resest(t, x, n, res, snr, sigma2)
!--------------------------------------------------------------------
!    VERSION: JUNE, 1996

!    PURPOSE:

!    COMPUTES ONE-LEAVE-OUT RESIDUALS FOR NONPARAMETRIC ESTIMATION
!    OF RESIDUAL VARIANCE (LOCAL LINEAR APPROXIMATION FOLLOWED BY
!    REWEIGHTING)

!  PARAMETERS:

!  INPUT   T(N)      ABSCISSAE (ORDERED: T(I)<=T(I+1))
!  INPUT   X(N)      DATA
!  INPUT   N         LENGTH OF DATA ( >2 )
!  OUTPUT  RES(N)    RESIDUALS AT T(1),...,T(N)
!  OUTPUT  SNR       EXPLAINED VARIANCE OF THE TRUE CURVE
!  OUTPUT  SIGMA2    ESTIMATION OF SIGMA**2 (RESIDUAL VARIANCE)

!--------------------------------------------------------------------

REAL (dp), INTENT(IN)   :: t(:)
REAL (dp), INTENT(IN)   :: x(:)
INTEGER, INTENT(IN)     :: n
REAL (dp), INTENT(OUT)  :: res(:)
REAL (dp), INTENT(OUT)  :: snr
REAL (dp), INTENT(OUT)  :: sigma2

! Local variables

INTEGER    :: i
REAL (dp)  :: dn, g1, g2, ex, ex2, sx, tt

!-
sigma2=0.
ex=x(1)*(t(2)-t(1))
ex2=x(1)*ex
DO  i=2,n-1
  tt=t(i+1)-t(i-1)
  IF(tt /= 0.) g1=(t(i+1)-t(i))/tt
  IF(tt == 0.) g1=.5
  g2=1.-g1
  res(i)=(x(i)-g1*x(i-1)-g2*x(i+1))/SQRT(1.+g1*g1+g2*g2)
  sigma2=sigma2+res(i)*res(i)
  sx=x(i)*tt
  ex=ex+sx
  ex2=ex2+x(i)*sx
END DO
tt=t(3)-t(2)
IF(tt /= 0.) g1=(t(1)-t(2))/tt
IF(tt == 0.) g1=.5
g2=1.-g1
res(1)=(x(1)-g1*x(3)-g2*x(2))/SQRT(1.+g1*g1+g2*g2)
tt=t(n-1)-t(n-2)
IF(tt /= 0.) g1=(t(n-1)-t(n))/tt
IF(tt == 0.) g1=.5
g2=1.-g1
res(n)=(x(n)-g1*x(n-2)-g2*x(n-1))/SQRT(1.+g1*g1+g2*g2)
sigma2=(sigma2+res(1)*res(1)+res(n)*res(n))/n
!-
sx=x(n)*(t(n)-t(n-1))
dn=2.*(t(n)-t(1))
ex=(ex+sx)/dn
ex2=(ex2+x(n)*sx)/dn
IF(ex2 == 0) snr=0.
IF(ex2 > 0) snr=1-sigma2/(ex2-ex*ex)
RETURN
END SUBROUTINE resest



SUBROUTINE kernel(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
!       SHORT-VERSION MAY, 1995

!       DRIVER SUBROUTINE FOR KERNEL SMOOTHING, CHOOSES BETWEEN
!       STANDARD AND O(N) ALGORITHM

!  PARAMETERS :

!  INPUT    T(N)         INPUT GRID (REGRESSION DESIGN)
!  INPUT    X(N)         DATA, GIVEN ON T(N)
!  INPUT    N            LENGTH OF X
!  INPUT    B            ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH)
!  INPUT    NUE          ORDER OF DERIVATIVE (0-4)
!  INPUT    KORD         ORDER OF KERNEL (<=6); DEFAULT IS KORD=NUE+2
!  INPUT    NY           0: GLOBAL BANDWIDTH (DEFAULT)
!                        1: VARIABLE BANDWIDTHS, GIVEN IN Y AS INPUT
!  INPUT    S(0:N)       INTERPOLATION SEQUENCE
!  INPUT    TT(M)        OUTPUT GRID.  MUST BE PART OF INPUT GRID FOR IEQ=0
!  INPUT    M            NUMBER OF POINTS WHERE FUNCTION IS ESTIMATED,
!                         OR  LENGTH OF TT. DEFAULT IS M=400
!  INPUT    Y(M)         BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0
!  OUTPUT   Y(M)         ESTIMATED REGRESSION FUNCTION

!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)      :: t(:)
REAL (dp), INTENT(IN)      :: x(:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: b
INTEGER, INTENT(IN)        :: nue
INTEGER, INTENT(IN)        :: kord
INTEGER, INTENT(IN)        :: ny
REAL (dp), INTENT(IN)      :: s(0:)
REAL (dp), INTENT(IN)      :: tt(:)
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN OUT)  :: y(:)

! Local variable

REAL (dp) :: chan

!-
!------  COMPUTING CHANGE POINT
chan=(5. + kord)*MAX(1.0, SQRT(REAL(n)/REAL(m)))
!------
IF(b*(n-1)/(t(n)-t(1)) < chan) THEN
  CALL kerncl(t, x, n, b, nue, kord, ny, s, tt, m, y)
ELSE
  CALL kernfa(t, x, n, b, nue, kord, ny, s, tt, m, y)
END IF
!-
RETURN
END SUBROUTINE kernel



SUBROUTINE kernp(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
!       SHORT-VERSION JANUARY, 1997

!       DRIVER SUBROUTINE FOR KERNEL SMOOTHING, CHOOSES BETWEEN
!       STANDARD AND O(N) ALGORITHM WITHOUT USING BOUNDARY KERNELS

!  PARAMETERS :

!  INPUT    T(N)         INPUT GRID (REGRESSION DESIGN)
!  INPUT    X(N)         DATA, GIVEN ON T(N)
!  INPUT    N            LENGTH OF X
!  INPUT    B            ONE SIDED BANDWIDTH (FOR NY=1 MEAN BANDWIDTH)
!  INPUT    NUE          ORDER OF DERIVATIVE (0-4)
!  INPUT    KORD         ORDER OF KERNEL (<=6); DEFAULT IS KORD=NUE+2
!  INPUT    NY           0: GLOBAL BANDWIDTH (DEFAULT)
!                        1: VARIABLE BANDWIDTHS, GIVEN IN Y AS INPUT
!  INPUT    S(0:N)       INTERPOLATION SEQUENCE
!  INPUT    TT(M)        OUTPUT GRID. MUST BE PART OF INPUT GRID FOR IEQ=0
!  INPUT    M            NUMBER OF POINTS WHERE FUNCTION IS ESTIMATED,
!                         OR  LENGTH OF TT.  DEFAULT IS M=400
!  INPUT    Y(M)         BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0
!  OUTPUT   Y(M)         ESTIMATED REGRESSION FUNCTION

!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)      :: t(:)
REAL (dp), INTENT(IN)      :: x(:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: b
INTEGER, INTENT(IN)        :: nue
INTEGER, INTENT(IN)        :: kord
INTEGER, INTENT(IN)        :: ny
REAL (dp), INTENT(IN)      :: s(0:)
REAL (dp), INTENT(IN)      :: tt(:)
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN OUT)  :: y(:)

! Local variable

REAL (dp) :: chan

!-
!------  COMPUTING CHANGE POINT
chan=(5.+kord)*MAX(1., SQRT(REAL(n)/REAL(m)))
!------
IF(b*(n-1)/(t(n)-t(1)) < chan) THEN
  CALL kerncp(t, x, n, b, nue, kord, ny, s, tt, m, y)
ELSE
  CALL kernfp(t, x, n, b, nue, kord, ny, s, tt, m, y)
END IF
!-
RETURN
END SUBROUTINE kernp



SUBROUTINE kernfa(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
!       SHORT-VERSION: MAY, 1995

!       PURPOSE:

!       COMPUTATION OF KERNEL ESTIMATE USING O(N) ALGORITHM BASED ON
!       LEGENDRE POLYNOMIALS, GENERAL SPACED DESIGN AND LOCAL
!       BANDWIDTH ALLOWED. (NEW INITIALISATIONS OF THE LEGENDRE SUMS
!       FOR NUMERICAL REASONS)

!  PARAMETERS :

!  INPUT    T(N)         INPUT GRID
!  INPUT    X(N)         DATA, GIVEN ON T(N)
!  INPUT    N            LENGTH OF X
!  INPUT    B            ONE SIDED BANDWIDTH
!  INPUT    NUE          ORDER OF DERIVATIVE (0-4)
!  INPUT    KORD         ORDER OF KERNEL (<=6)
!  INPUT    NY           0, GLOBAL BANDWIDTH; 1, LOCAL BANDWIDTH IN Y
!  INPUT    S(0:N)       HALF POINT INTERPOLATION SEQUENCE
!  INPUT    TT(M)        OUTPUT GRID
!  INPUT    M            NUMBER OF POINTS TO ESTIMATE
!  INPUT    Y(M)         BANDWITH SEQUENCE FOR NY=1, DUMMY FOR NY=0
!  OUTPUT   Y(M)         ESTIMATED FUNCTION

!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)      :: t(:)
REAL (dp), INTENT(IN)      :: x(:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: b
INTEGER, INTENT(IN)        :: nue
INTEGER, INTENT(IN)        :: kord
INTEGER, INTENT(IN)        :: ny
REAL (dp), INTENT(IN)      :: s(0:n)
REAL (dp), INTENT(IN)      :: tt(:)
INTEGER, INTENT(IN)        :: m
REAL (dp), INTENT(IN OUT)  :: y(:)


INTEGER   :: j, k, iord, init, icall, i, iboun
INTEGER   :: jl, jr, jnr, jnl

REAL (dp) :: c(7), sw(7), xf(7), dold
REAL (dp) :: a(7,7), a1(7), a2(7), a3(7,7), cm(7,6)
REAL (dp) :: bmin, bmax, bb, s0, sn, wwl, wwr, wid, wr, wido
!-
!------ COMPUTE CONSTANTS FOR LATER USE
s0=1.5*t(1)-0.5*t(2)
sn=1.5*t(n)-0.5*t(n-1)
bmin=(sn-s0)*.6D0/DBLE(n)*DBLE(kord-1)
bmax=(s(n)-s(0))*.5
IF(kord == 2) bmin=bmin*0.1D0
iord=kord+1
DO  k=3,iord
  a1(k)=DBLE(2*k-1)/DBLE(k)
  a2(k)=DBLE(1-k)/DBLE(k)
END DO
!-
init=0
icall=0
dold=0.d0
!-
!------ SMOOTHING LOOP
DO  i=1,m
  bb=b
  IF (ny == 1) bb=y(i)
  IF(bb < bmin) bb=bmin
  IF(bb > bmax) bb=bmax
  iboun=0
!-
!------ COMPUTE LEFT BOUNDARY KERNEL
  IF(tt(i) < s(0)+bb) THEN
    wwl=s(0)
    wwr=s(0)+bb+bb
    wid=wwr-tt(i)
    iboun=1
    CALL coffb(nue,kord,(tt(i)-s(0))/wid,iboun,c)
  END IF
!-
!------ COMPUTE RIGHT BOUNDARY KERNEL
  IF(tt(i)+bb > s(n)) THEN
    wwl=s(n)-(bb+bb)
    wwr=s(n)
    wid=tt(i)-wwl
    iboun=-1
    CALL coffb(nue,kord,(s(n)-tt(i))/wid,iboun,c)
  END IF
!-
!------ NO BOUNDARY
  IF(iboun == 0) THEN
    wid=bb
    wwl=tt(i)-bb
    wwr=tt(i)+bb
  END IF
!-
!------ INITIALISATION FOR INIT=0
  IF(init == 0) THEN
    sw(1:iord)=0.
    jl=1
    DO  j=1,n
      IF(s(j-1) < wwl) THEN
        jl=j+1
      ELSE
        IF(s(j) > wwr) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
      END IF
    END DO
    jr=j-1
    wr=wwr
    init=1
    GO TO 6666
  ELSE
    init=init+1
  END IF
!-
!------ COMPARE OLD SUM WITH NEW SMOOTHING INTERVALL TT(I)-B,TT(I)+B
  IF(s(jr-1) >= wwl) THEN
    jnr=jr
    jnl=jl
    IF(s(jr) > wwr) THEN
      DO  j=jr,jl,-1
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i-1),wido,-1)
        jnr=j-1
        IF(s(jnr) <= wwr) EXIT
      END DO
    END IF
    IF(s(jl-1) < wwl) THEN
      DO  j=jl,jr
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i-1),wido,-1)
        jnl=j+1
        IF(s(j) >= wwl) EXIT
      END DO
    END IF
!-
!------ UPDATING OF SW
    CALL lreg(sw,a3,iord,(tt(i)-tt(i-1))/wid,dold,wido/wid,cm)
    IF(jnr == jr) THEN
      DO  j=jr+1,n
        IF(s(j) > wwr) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
        jnr=j
      END DO
    END IF
    jr=jnr
    IF(jl == jnl) THEN
      DO   j=jl-1,1,-1
        IF(s(j-1) < wwl) EXIT
        CALL dreg(sw,a1,a2,iord,x(j),s(j-1),s(j),tt(i),wid,1)
        jnl=j
      END DO
    END IF
    jl=jnl
  ELSE
!-
!------ NEW INITIALISATION OF SW
    sw(1:iord)=0.
    DO  j=jr,n