📄 kern_reg.f90
字号:
IF(MOD(i-k,2) == 0) ww=ww+a3(i,k)
xx=xx+ww*sw(k)
END DO
sw(i)=xx
END DO
sw(1)=a3(1,1)*sw(1)
ELSE
DO i=iord,2,-1
DO k=1,i-1
sw(i)=sw(i)+c(i,k)*sw(k)
END DO
END DO
END IF
RETURN
END SUBROUTINE lreg
SUBROUTINE freg(sw, nue, kord, iboun, y, c, icall, a)
!------------------------------------------------------------------
! SHORT-VERSION: MAY, 1995
! PURPOSE:
! FINAL COMPUTATION OF A SMOOTHED VALUE VIA LEGENDRE POLYNOMIALS
! PARAMETERS :
! ************** INPUT *******************
! SW(KORD+1): SUM OF DATA WEIGHTS FOR LEGENDRE POLYNOM.
! NUE : ORDER OF DERIVATIVE
! KORD : ORDER OF KERNEL
! IBOUN : 0: INTERIOR KERNEL, ELSE BOUNDARY KERNEL
! C(KORD+1) : SEQUENCE OF POLYN. COEFF. FOR BOUND. KERNEL
! ICALL : PARAMETER USED TO INITIALISE COMPUTATION
! : OF A MATRIX
! ************** WORK ******************
! A(7,7) : MATRIX OF COEFFICIENTS
! ************** OUTPUT ******************
! Y : COMPUTED ESTIMATE
!--------------------------------------------------------------------
REAL (dp), INTENT(IN) :: sw(7)
INTEGER, INTENT(IN) :: nue
INTEGER, INTENT(IN) :: kord
INTEGER, INTENT(IN) :: iboun
REAL (dp), INTENT(OUT) :: y
REAL (dp), INTENT(IN) :: c(7)
INTEGER, INTENT(IN OUT) :: icall
REAL (dp), INTENT(OUT) :: a(7,7)
INTEGER :: i, j
REAL (dp) :: ww
!-
!------- DEFINITION OF LEGENDRE COEFFICIENTS FOR BOUNDARY
IF(icall == 0 .AND. iboun /= 0) THEN
a(2,2)=2./3.
a(1,3)=.6
a(3,3)=.4
a(2,4)=4./7.
a(4,4)=8./35.
a(1,5)=27./63.
a(3,5)=28./63.
a(5,5)=8./63.
a(2,6)=110./231.
a(4,6)=72./231.
a(6,6)=16./231.
a(1,7)=143./429.
a(3,7)=182./429.
a(5,7)=88./429.
a(7,7)=16./429.
icall=1
END IF
IF(iboun /= 0) THEN
!-
!------- COMPUTATION OF THE SMOOTHED VALUE AT BOUNDARY
y=c(1)*sw(1)+c(2)*a(2,2)*sw(2)
DO j=3,kord+1
ww=a(j,j)*sw(j)
DO i=j-2,1,-2
ww=ww+a(i,j)*sw(i)
END DO
y=y+c(j)*ww
END DO
ELSE
!-
!------- COMPUTATION OF THE SMOOTHED VALUE AT INTERIOR
IF(nue == 0) THEN
IF(kord == 2) y=-.1*sw(3)+.6*sw(1)
IF(kord == 4) y=(sw(5)-4.*sw(3)+9.*sw(1))/12.
IF(kord == 6) y=-7.2115379E-02*sw(7)+.25961537*sw(5) &
-.4375*sw(3)+.75*sw(1)
END IF
IF(nue == 1) THEN
IF(kord == 3) y=(3.*sw(4)-10.*sw(2))/14.
IF(kord == 5) y=(-15.*sw(6)+48.*sw(4)-55.*sw(2))/44.
END IF
IF(nue == 2) THEN
IF(kord == 4) y=(-5.*sw(5)+14.*sw(3)-9.*sw(1))/6.
IF(kord == 6) y=2.01923*sw(7)-5.76923*sw(5)+5.25*sw(3) -1.5*sw(1)
END IF
IF(nue == 3) y=4.772727*sw(6)-12.272727*sw(4)+7.5*sw(2)
IF(nue == 4) y=-36.34615*sw(7)+88.84615*sw(5)-52.5*sw(3)
END IF
!-
RETURN
END SUBROUTINE freg
SUBROUTINE kerncl(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
! SHORT-VERSION JANUARY 1995
! KERNEL SMOOTHING, CONVENTIONAL ALGORITHM,GENERAL
! DESIGN, LOCAL BANDWIDTH ALLOWED
! PARAMETERS :
! INPUT T(N) INPUT GRID
! 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)
! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE
! INPUT NY 0: GLOBAL, 1: LOCAL BANDWIDTH INPUTED IN Y
! 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 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:n)
REAL (dp), INTENT(IN) :: tt(:)
INTEGER, INTENT(IN) :: m
REAL (dp), INTENT(IN OUT) :: y(:)
REAL (dp) :: c(7), c1(7)
INTEGER :: ist, i, iboun, iord
REAL (dp) :: bb, bmax, wid, s0, s1, sn, bmin
!-
!------ COMPUTE KERNEL COEFFICIENTS FOR INTERIOR AND SOME CONSTANTS
CALL coffi(nue, kord, c)
iord=kord+1
bb=b
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=0.1D0*bmin
ist=1
!-
!------- LOOP OVER OUTPUT GRID
DO i=1,m
IF(ny /= 0) bb=y(i)
IF(bb > bmax) bb=bmax
IF(bb < bmin) bb=bmin
wid=bb
s1=tt(i)-bb
iboun=0
!-
!------- COMPUTE LEFT BOUNDARY KERNEL
IF(s1 < s(0)) THEN
s1=s(0)
wid=bb+bb+s(0)-tt(i)
CALL coffb(nue,kord,(tt(i)-s(0))/wid,1,c1)
iboun=1
END IF
!-
!------- COMPUTE RIGHT BOUNDARY KERNEL
IF(tt(i)+bb > s(n)) THEN
s1=s(n)-(bb+bb)
wid=tt(i)-s1
CALL coffb(nue,kord,(s(n)-tt(i))/wid,-1,c1)
iboun=-1
END IF
!-
!------ SEARCH FIRST S-POINT OF SMOOTHING INTERVAL
2 IF(s(ist) <= s1) THEN
ist=ist+1
GO TO 2
END IF
3 IF(s(ist-1) > s1) THEN
ist=ist-1
GO TO 3
END IF
!-
!------- IF BANDWIDTH IS TOO SMALL NO SMOOTHING
IF(s(ist) >= tt(i)+wid .OR. ist == n) THEN
y(i)=x(ist)
IF(nue > 0) y(i)=0.
ELSE
!-
!----- COMPUTE SMOOTHED DATA AT TT(I)
IF (iboun /= 0) THEN
CALL smo(s, x, n, tt(i), wid, nue, iord, iboun, ist, s1, c1, y(i))
ELSE
CALL smo(s, x, n, tt(i), wid, nue, iord, iboun, ist, s1, c, y(i))
END IF
END IF
END DO
!-
RETURN
END SUBROUTINE kerncl
SUBROUTINE smo(s, x, n, tau, wid, nue, iord, iboun, ist, s1, c, y)
!-----------------------------------------------------------------------
! SHORT-VERSION JANUARY 1995
! PERFORMS ONE SMOOTHING STEP
! PARAMETERS :
! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE
! INPUT X(N) DATA
! INPUT N LENGTH OF X
! INPUT TAU POINT WHERE FUNCTION IS ESTIMATED
! INPUT WID ONE SIDED BANDWIDTH
! INPUT NUE ORDER OF DERIVATIVE (0-4)
! INPUT IORD ORDER OF KERNEL POLYNOMIAL
! INPUT IBOUN TYPE OF BOUNDARY
! INPUT IST INDEX OF FIRST POINT OF SMOOTHING INTERVAL
! INPUT S1 LEFT BOUNDARY OF SMOOTHING INTERVAL
! INPUT C(7) KERNEL COEFFICIENTS
! OUTPUT Y SMOOTHED VALUE AT TAU
! WORK WO(7) WORK ARRAY
!-----------------------------------------------------------------------
REAL (dp), INTENT(IN) :: s(0:)
REAL (dp), INTENT(IN) :: x(:)
INTEGER, INTENT(IN) :: n
REAL (dp), INTENT(IN) :: tau
REAL (dp), INTENT(IN) :: wid
INTEGER, INTENT(IN) :: nue
INTEGER, INTENT(IN) :: iord
INTEGER, INTENT(IN) :: iboun
INTEGER, INTENT(IN) :: ist
REAL (dp), INTENT(IN) :: s1
REAL (dp), INTENT(IN) :: c(7)
REAL (dp), INTENT(OUT) :: y
REAL (dp) :: wo(7)
INTEGER :: jend, ibeg, incr, i, j
REAL (dp) :: yy, yyy, w, widnue
!-
y=0.
jend=0
ibeg=2
IF(iboun /= 0 .OR. (nue /= 1 .AND. nue /= 3)) ibeg=1
incr=2
IF(iboun /= 0) incr=1
!-
!------ COMPUTE INITIAL KERNEL VALUES
IF(iboun > 0) THEN
yy=(tau-s1)/wid
wo(ibeg)=yy
DO i=ibeg,iord-incr,incr
wo(i+incr)=wo(i)*yy
END DO
ELSE
DO i=ibeg,iord,incr
wo(i)=1.
END DO
END IF
!-
!------ LOOP OVER SMOOTHING INTERVAL
DO j=ist,n
yy=(tau-s(j))/wid
IF(yy < -1.) THEN
yy=-1.
jend=1
END IF
yyy=yy
IF(iboun == 0) THEN
yy=yy*yy
IF(nue == 1 .OR. nue == 3) yyy=yy
END IF
!-
!------ LOOP FOR COMPUTING WEIGHTS
w=0.
DO i=ibeg,iord,incr
w=w+c(i)*(wo(i)-yyy)
wo(i)=yyy
yyy=yyy*yy
END DO
y=y+w*x(j)
IF(jend == 1) EXIT
END DO
!-
!------- NORMALIZING FOR NUE>0
IF(nue > 0) THEN
widnue=wid**nue
y=y/widnue
END IF
!-
RETURN
END SUBROUTINE smo
SUBROUTINE kerncp(t, x, n, b, nue, kord, ny, s, tt, m, y)
!-----------------------------------------------------------------------
! SHORT-VERSION JANUARY 1997
! KERNEL SMOOTHING, CONVENTIONAL ALGORITHM,GENERAL
! DESIGN, LOCAL BANDWIDTH ALLOWED, WITHOUT BOUNDARY KERNELS,
! JUST NORMALIZING
! PARAMETERS :
! INPUT T(N) INPUT GRID
! 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)
! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE
! INPUT NY 0: GLOBAL, 1: LOCAL BANDWIDTH INPUTED IN Y
! 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 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:n)
REAL (dp), INTENT(IN) :: tt(:)
INTEGER, INTENT(IN) :: m
REAL (dp), INTENT(IN OUT) :: y(:)
REAL (dp) :: c(7), c1(7)
INTEGER :: ist, i, iboun, iord
REAL (dp) :: bb, bmax, wid, s0, s1, sn, bmin
!-
!------ COMPUTE KERNEL COEFFICIENTS FOR INTERIOR AND SOME CONSTANTS
CALL coffi(nue, kord, c)
iord=kord+1
bb=b
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=0.1D0*bmin
ist=1
!-
!------- LOOP OVER OUTPUT GRID
DO i=1,m
IF(ny /= 0) bb=y(i)
IF(bb > bmax) bb=bmax
IF(bb < bmin) bb=bmin
wid=bb
s1=tt(i)-bb
iboun=0
!-
!------- COMPUTE LEFT BOUNDARY KERNEL
IF(s1 < s(0)) THEN
s1=s(0)
wid=bb+bb+s(0)-tt(i)
CALL coffb(nue,kord,(tt(i)-s(0))/wid,1,c1)
iboun=1
END IF
!-
!------- COMPUTE RIGHT BOUNDARY KERNEL
IF(tt(i)+bb > s(n)) THEN
s1=s(n)-(bb+bb)
wid=tt(i)-s1
iboun=-1
END IF
!-
!------ SEARCH FIRST S-POINT OF SMOOTHING INTERVAL
2 IF(s(ist) <= s1) THEN
ist=ist+1
GO TO 2
END IF
3 IF(s(ist-1) > s1) THEN
ist=ist-1
GO TO 3
END IF
!-
!------- IF BANDWIDTH IS TOO SMALL NO SMOOTHING
IF(s(ist) >= tt(i)+wid .OR. ist == n) THEN
y(i)=x(ist)
IF(nue > 0) y(i)=0.
ELSE
!-
!----- COMPUTE SMOOTHED DATA AT TT(I)
CALL smop(s, x, n, tt(i), wid, nue, iord, iboun, ist, s1, c, y(i))
END IF
END DO
!-
RETURN
END SUBROUTINE kerncp
SUBROUTINE smop(s, x, n, tau, wid, nue, iord, iboun, ist, s1, c, y)
!-----------------------------------------------------------------------
! SHORT-VERSION JANUARY 1997
! PERFORMS ONE SMOOTHING STEP
! PARAMETERS :
! INPUT S(0:N) HALF POINT INTERPOLATION SEQUENCE
! INPUT X(N) DATA
! INPUT N LENGTH OF X
! INPUT TAU POINT WHERE FUNCTION IS ESTIMATED
! INPUT WID ONE SIDED BANDWIDTH
! INPUT NUE ORDER OF DERIVATIVE (0-4)
! INPUT IORD ORDER OF KERNEL POLYNOMIAL
! I⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -