kernf77.txt

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

      DOUBLE PRECISION B0,B1,B2,Y0(M),Y1(M),Y2(M)
      DOUBLE PRECISION WN(0:N,5),W1(M1,3)

      DATA TL,TU/1.0,0.0/,S/N*0.0,-1.0/,SIG/-1.0/
      OPEN(81,FILE='raw.dat')
      DATA X / 1.9666, 1.9000, 1.6449, 1.4275, 2.4000, 1.8487,
     . 1.3383, 1.9514, 1.5722, 1.1978, 1.3109, 1.6940, 1.0816,
     . 1.1070, 1.2650, 1.1143,  .5717,  .9492,  .3116,  .7942,
     .  .5670,  .3555,  .3243,  .3242,  .7647,  .4167,  .7583,
     .  .9462, 1.4337, 1.8820, 2.3054, 2.1040, 3.8500, 3.7052,
     . 4.3047, 4.4505, 4.3425, 4.3516, 3.7578, 4.0112, 3.3286,
     . 2.9024, 2.4331, 1.5544, 1.1867,  .6870, -.3630,  .0436,
     . -.8197, -.2883,-1.0559, -.9795,-1.7802,-2.2206,-1.3922,
     .-1.6511,-1.7694,-1.1309,-2.1912,-1.5785,-2.6189,-1.8125,
     .-2.6155,-1.4585,-2.0951,-2.1428,-2.4827,-2.4171,-2.6610,
     .-2.6509,-2.6475,-3.1040,-2.7631,-2.9486,-3.1323/
      DO 10 I=1,N
        T(I)=(I-0.5)/DBLE(N)
10    WRITE(81,'(2F9.5)') T(I),X(I)
      CLOSE(81)
      DO 20 J=1,M
20      TT(J)=DBLE(J-1)*(T(N)-T(1))/DBLE(M-1)+T(1)

      NUE=0 
      KORD=NUE+2
      CALL GLKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,B0,Y0)
      NUE=1 
      KORD=NUE+2
      S(0)=1.0
      S(N)=0.0
      CALL GLKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,B1,Y1)
      NUE=2 
      KORD=NUE+2
      S(0)=1.0
      S(N)=0.0
      CALL GLKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,B2,Y2)
      PRINT *,'Plugin bandwidths used in glkern:'
      PRINT *,'regression function:',B0
      PRINT *,'first derivative:   ',B1
      PRINT *,'second derivative:  ',B2
      OPEN(82,FILE='est.dat')
      WRITE(82,'(4G14.6)') (TT(J),Y0(J),Y1(J),Y2(J),J=1,M)       
      CLOSE(82)
999   STOP
      END
The output of the code will give

 Plugin bandwidths used in glkern:
 regression function: 4.845240675004461E-02
 first derivative:    8.135588838823485E-02
 second derivative:   .1190205218068572
The estimated functions are written in a file named 'est.dat'. The first lines 
of this data-file should contain 
   .666667E-02   1.88514      -.412435      -38.6205
   .996656E-02   1.87886      -.564513      -38.5011
   .132664E-01   1.87257      -.716592      -38.3817
   .165663E-01   1.86628      -.868670      -38.2623
   .198662E-01   1.86000      -1.02075      -38.1429
This file may be handled by some graphics software. An example of an figure 
created with gnuplot can be find below or alternatively as ps-file.
 
2. A second example calls glkern and lokern in order to smooth a data set of 250 
points which are read from an input file 'raw.dat'. Here, we suppose that they 
stem from a random sample but are already ordered in the first variable. As 
before, the output will be written in a file 'est.dat'.

      PROGRAM
      INTEGER N,M,IHOM,NUE,KORD,IRND,ISMO,M1,I,J
      PARAMETER(N=250,M=300,IHOM=0,NUE=0,KORD=NUE+2)
      PARAMETER(IRND=0,ISMO=0,M1=400)
      DOUBLE PRECISION T(N),X(N),TT(M),TL,TU,S(0:N),SIG,B
      DOUBLE PRECISION YG(M),YL(M),BAN(M),WN(0:N,5),W1(M1,3),WM(M)

      DATA TL,TU/1.0,0.0/,S/N*0.0,-1.0/,SIG/-1.0/
      OPEN(81,FILE='raw.dat')
      DO 10 I=1,N
10    READ(81,*) T(I),X(I)
      CLOSE(81)
      DO 20 J=1,M
20      TT(J)=DBLE(J-1)*(T(N)-T(1))/DBLE(M-1)+T(1)

      CALL GLKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,B,YG)
      CALL LOKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,WM,BAN,YL)

      OPEN(82,FILE='est.dat')
      WRITE(82,'(4G14.6)') (TT(J),YG(J),YL(J),BAN(J),J=1,M)       
      CLOSE(82)
999   STOP
      END
3. The third example calls lokern in order to smooth a data set of 250 points 
with an input bandwidth array. The data and the bandwidth array with an 
associated output grid is read from an input file 'raw.dat' and 
'band.dat',respectively, whereas the output will be written in a file 'est.dat' 
as before. Here, we assume a data set from (ordered) random design points with 
heteroscedastic variance.

      PROGRAM
      INTEGER N,M,IHOM,NUE,KORD,IRND,ISMO,M1,I,J
      PARAMETER(N=250,M=300,IHOM=1,NUE=0,KORD=NUE+2)
      PARAMETER(IRND=0,ISMO=1,M1=400)
      DOUBLE PRECISION T(N),X(N),TT(M),TL,TU,S(0:N),SIG
      DOUBLE PRECISION Y(M),BAN(M),WN(0:N,5),W1(M1,3),WM(M)

      DATA TL,TU/1.0,0.0/,S/N*0.0,-1.0/
      OPEN(81,FILE='raw.dat')
      DO 10 I=1,N
10    READ(81,*) T(I),X(I)
      CLOSE(81)
      OPEN(82,FILE='band.dat')
      DO 20 J=1,M
20    READ(82,*) TT(J),BAN(J)
      CLOSE(82)
      
      CALL LOKERN(T,X,N,TT,M,IHOM,NUE,KORD,IRND,ISMO,M1,TL,TU,S,
     .   SIG,WN,W1,WM,BAN,Y)

      OPEN(82,FILE='est.dat')
      WRITE(82,'(2G14.6)') (TT(J),Y(J),J=1,M)       
      CLOSE(82)
999   STOP
      END
Goto top, general remarks, computational remarks, recommendations, statistical 
remarks, and references.



Statistical Remarks
The proposed algorithms for kernel regression estimation performed quite well in 
many simulations and in asymptotic theory. Nevertheless, there can be several 
reasons why the smoothing results may behave worse in special examples. Some of 
them are mentioned below. Also computing errors or instabilities may be 
observed, please inform us about any suspect error. A critical look at the 
results is always necessary. 
    The method is sensitive to correlations among the error variables. If some 
    kind of positive correlation occurs, the typical result will be an 
    undersmoothed regression function. There exist some proposals on how to 
    adapt some kinds of dependencies but they are not included in the algorithms 
    used by glkern.f and lokern.f up to now. 
    Typically, any derivative estimator has severe problems in the boundary 
    regions. 
    In situations where a very good kernel regression estimator exists, e.g. for 
    very large sample size, small residual variance, very smooth regression 
    functions, estimating optimal bandwidths with the plugin method seems to be 
    more difficult. Here, we observed sometimes clearly suboptimal results for 
    the automatic bandwidths. Nevertheless, the fits might still be good in 
    absolute error. 
    The method attempts to optimize bandwidth choice with respect to an 
    integrated squared error criterion. This may be significantly different from 
    other criteria, e.g. average squared error for non-equidistant design. 
Of course, we are interested in improving our method. Please let us know if you 
have problems in applying the software and if you observe unsatisfactory 
results. 
Goto top, general remarks, computational remarks, recommendations, examples, and 
references.



References
On the global iterative plugin bandwidth estimator:
T. Gasser, A. Kneip, and W. K鰄ler (1991). A flexible and fast method for 
automatic smoothing. Journal of the American Statistical Association, 86, 
643-652. 
On the local plugin bandwidth estimator:
M. Brockmann, T. Gasser, and E. Herrmann (1993). Locally adaptive bandwidth 
choice for kernel regression estimators. Journal of the American Statistical 
Association, 88, 1302-1309. 
On nonparametric variance estimation:
T. Gasser, L. Sroka, and C. Jennen-Steinmetz (1986). Residual and residual 
pattern in nonlinear regression. Biometrika, 73, 625-633. 
On adapting heteroscedasticity:
E. Herrmann (1997). Local bandwidth choice in kernel regression estimation. 
Journal of Graphical and Computational Statistics, 6, 35-54. 
On the fast algorithm for kernel regression estimator:
T. Gasser and A. Kneip (1989) discussion of Buja, A., Hastie, T. and Tibshirani, 
R.: Linear smoothers and additive models, The Annals of Statistics, 17, 532-535. 

B. Seifert, M. Brockmann, J. Engel, and T. Gasser (1994). Fast algorithms for 
nonparametric curve estimation. J. Computational and Graphical Statistics 3, 
192-213. 
On the special kernel estimator for random design:
E. Herrmann (1996). On the convolution type kernel regression estimator. 
Preprint 1833, FB Mathematik, Technische Universit鋞 Darmstadt (available from 
the preprint server of the mathematical department of the technical university 
of Darmstadt )
Goto top, general remarks, computational remarks, recommendations, examples, and 
statistical remarks.


[ University of Zurich | Biostatistics | Software Page ]

 Comments to Eva Herrmann: [e-mail]

 Last update: 10-December-98 / mz