C*********************************************************************** C THIS SUBROUTINE CALCULATES AN ESTIMATE OF F(0) AND ITS VARIANCE C USING A FOURIER SERIES. IT NEEDS AS INPUT THE UNGROUPED DATA C VECTOR - X,THE MAXIMUM VALUE OF THE DATA - WIDTH, THE NUMBER C OF DATA POINTS - CNT, AND WHETHER THE NUMBER OF PARAMETERS HAS C BEEN SPECIFIED (SET) AND HOW MANY IF IT IS (NPSET) - ALL ARE C IN COMMON. C C SUBROUTINES CALLED: HEADER,NAME C*********************************************************************** SUBROUTINE FSEST (IND) INCLUDE 'PARMTR.INC' C*********************************************************************** C DECLARATIONS C*********************************************************************** INTEGER CNT LOGICAL YES, WARN, MSET, STRT, SET C*********************************************************************** C COMMON STATEMENTS C*********************************************************************** COMMON /PDOPT/ STRT, SET, PSTRT(10,MAXPAR), NPSET(10) COMMON /SOLN/ FZERO, VARF, FMAX, FMIN, FZ(100) COMMON /NUM/ XL(MAXLIN), WIDTH, N, CNT, CONV(3), VARN, IDF, WARN COMMON /MEASUR/ X(MAXOBJ), AA(MAXOBJ,4) DOUBLE PRECISION PAR, VCMAT, G, XLL COMMON /DPAR/ PAR(MAXPR2), VCMAT(MAXPAR,MAXPAR), G(MAXPAR), 1 XLL, NPAR, INDEX C*********************************************************************** C INITIALIZE VARIABLES AND WRITE OUT HEADER STATEMENTS C*********************************************************************** NP=MAXPAR MSET=SET XN=FLOAT(CNT) YES=.FALSE. FZERO=1./WIDTH NPAR=0 IF (NPSET(IND).EQ.0) MSET=.FALSE. IF (MSET) NPAR=NPSET(IND) CALL HEADER (1) CALL NAME (INDEX) WRITE (6,130) I=0 STOP=SQRT(2./(XN+1.))/WIDTH DO 10 K=1,NP*2 10 PAR(K)=0.D0 C*********************************************************************** C CALCULATE PARAMETER ESTIMATES C*********************************************************************** 20 I=I+1 DO 30 K=1,CNT 30 PAR(I)=PAR(I)+COS((FLOAT(I)*3.14159*X(K))/WIDTH) PAR(I)=2.*PAR(I)/(XN*WIDTH) C*********************************************************************** C THE LOGICAL VARIABLE YES IS SET TO TRUE IF THE NUMBER OF PARAMETER C IN THE MODEL HAS BEEN CALCULATED. C THEN PAR(NPAR+1),...,PAR(2*NPAR) ARE CALCULATED TO GET VARIANCES. C*********************************************************************** IF (YES.AND.(I.LT.2*NPAR)) GO TO 20 IF (YES.AND.(I.GE.2*NPAR)) GO TO 60 FZERO=FZERO+PAR(I) WRITE (6,100) I,PAR(I),FZERO,STOP C*********************************************************************** C IF THE NUMBER OF PARAMETERS HAS BEEN SET THEN SEE IF THAT NUMBER C HAS BEEN REACHED. C*********************************************************************** IF (.NOT.MSET) GO TO 40 IF (I.EQ.NPAR) YES=.TRUE. WRITE (6,150) GO TO 20 C*********************************************************************** C SEE IF THE ABSOLUTE VALUE OF THE PARAMETER ESTIMATE IS GREATER C THAN THE STOPPING RULE VALUE TO DETERMINE IF THE PARAMETER SHOULD C BE INCLUDED IN THE MODEL. C*********************************************************************** 40 IF (DABS(PAR(I)).GE.STOP) GO TO 50 YES=.TRUE. NPAR=I-1 FZERO=FZERO-PAR(I) WRITE (6,150) WRITE (6,110) NPAR GO TO 20 C*********************************************************************** C FOR THE STOPPING RULE APPROACH NEVER LET THE NUMBER OF PARAMETERS C EXCEED 6. C*********************************************************************** 50 IF (I.LT.6) GO TO 20 YES=.TRUE. NPAR=I WRITE (6,150) WRITE (6,120) GO TO 20 C*********************************************************************** C NEXT CALCULATE THE COVARIANCE MATRIX OF THE COEFFICIENTS IN THE C SERIES AND PAR VARIANCE FOR F(0). C*********************************************************************** 60 VARF=0.0 CALL HEADER (1) CALL NAME (INDEX) IF ((NPAR.EQ.0).OR.(XN.EQ.1.)) RETURN DO 90 I=1,NPAR IF (I.EQ.NPAR) GO TO 80 DO 70 J=I+1,NPAR VCMAT(I,J)=((PAR(I+J)+PAR(J-I))/DBLE(WIDTH) 1 -PAR(I)*PAR(J))/DBLE(XN-1.) VCMAT(J,I)=VCMAT(I,J) 70 VARF=VARF+2.*SNGL(VCMAT(I,J)) 80 I2=I*2 VCMAT(I,I)=((PAR(I2)+2./DBLE(WIDTH))/DBLE(WIDTH) 1 -PAR(I)*PAR(I))/DBLE(XN-1.) 90 VARF=VARF+SNGL(VCMAT(I,I)) C*********************************************************************** C FORMAT STATEMENTS C*********************************************************************** RETURN C 100 FORMAT (4X,'A(',I2.2,')',4X,G12.5,6X,G12.5,7X,G15.8) 110 FORMAT ('0The absolute value of the parameter is less', 1' than the stopping rule value.'/ 2' Thus ',I2,' term(s) will be used in the model.') 120 FORMAT ('0The stopping rule may have allowed more than 6', 1' terms in the model.'/ 2' This porgram will only use at the most 6 terms unless', 3' the user specifies a'/ 4' set number of terms greater than 6.') 130 FORMAT (//'0The ungrouped version of the Fourier series', 1' estimator determines the'/ 2' number of terms by a stopping rule which minimizes', 3' the mean integrated'/ 4' squared error. Below the value for each parameter,', 5' the cumulative value'/ 6' of F(0) and the stopping rule value are given.'/ 7'0Parameter Estimate Cumulative F(0)', 5 4X,'Stopping Rule Value'/1X,66(1H-)) 150 FORMAT (1X,66(1H-)) END