var/da/da_biascorr_airmass/regress_one.f90

   1   SUBROUTINE REGRESS_ONE(NX,XBAR,YBAR,XCOV,YCOV,XYCOV,NXEIGS,CCOF,CCOF0,RESERR,RESCOV,XVEC)
   2
   3   USE RAD_BIAS
   4   IMPLICIT NONE
   5
   6 ! ROUTINE          CALCF2  SUBROUTINE  *****  |EYRE.TOVFFG@
   7
   8 ! PURPOSE          EIGENVECTOR REGRESSION ROUTINE
   9
  10 ! VERSION          3.00,190491,J.R.EYRE
  11
  12 ! DESCRIPTION      ROUTINE TO CALCULATE EIGENVECTOR REGRESSION COEFFICIENT
  13 !                  Y AGAINST VECTOR X FROM COVARIANCE MATRICES XYCOV AND
  14 !                  XCOV, WHICH CONTAINS THE COVARIANCE MATRICES OF Y AND X
  15
  16 ! ARGUMENTS        XBAR   R*8  MEAN VECTOR OF PREDICTORS, LENGTH NX
  17 !                  YBAR   R*8  MEAN VALUE OF PREDICTANT
  18 !                  XCOV   R*8  COVARIANCE MATRIX (NX*NX) OF PREDICTORS
  19 !                  XYCOV  R*8  CROSS-COVARIANCE MATRIX (NX) OF PREDICTORS/PREDICTANT
  20 !                  NX     I*4  NUMBER OF PREDICTORS
  21 !                  NXEIGS  I*4  NUMBER OF EIGENVECTORS OF PREDICTOR USED IN
  22 !                              THE REGRESSION EQUATION
  23 !                  CCOF   R*8  OUTPUT MATRIX (NX) OF REGRESSION COEFF
  24 !                              SUCH THAT: Y = CCOF0 + CCOF.X
  25 !                  CCOF0  R*8  OFFSET IN ABOVE EQUATION
  26 !                  RESERR R*8  RESIDUAL ERROR IN PREDICANT
  27 !                  RESCOV R*8  RESIDUAL VARIANCE IN PREDICTANT
  28 !                  XVEC   R*8  OUTPUT MATRIX (NX*NX) OF PREDICTOR EIGENVECTORS
  29
  30 ! SUBPROGRAMS      F02ABF, F01ABF (nag)
  31
  32 ! FILES            UNIT(6)  DIAGNOSTIC OUTPUT
  33
  34 ! 14 JAN 88. MADISON VERSION CHANGED BACK FOR HOMER.
  35 ! 13 MAY 97. REWRITTEN IN F90 (B.Harris)
  36
  37   INTEGER, PARAMETER :: JWORK=4900
  38   INTEGER, PARAMETER :: JMPRED=70
  39   INTEGER            :: NXEIG
  40
  41   INTEGER, INTENT(IN) :: NX, NXEIGS
  42   INTEGER NXE
  43
  44   REAL(KIND=LONG),    INTENT(IN) :: XBAR(NX), YBAR, XCOV(NX,NX), XYCOV(NX), YCOV
  45
  46   REAL(KIND=LONG), INTENT(OUT) :: CCOF(NX), CCOF0
  47   REAL(KIND=LONG), INTENT(OUT) :: RESERR, RESCOV
  48   REAL(KIND=LONG), INTENT(OUT) :: XVEC(NX,NX)
  49
  50   REAL(KIND=LONG) :: XDEV(NX), XXBAR(NX)
  51   REAL(KIND=LONG) :: YVEC, XVAL(NX), YVAL, DET, D(NX)
  52   REAL(KIND=LONG) :: XW1(NX,NX), XW2(NX,NX), YW1, XYW1(NX), YXW1(NX), XINV(NX,NX), XYW2(NX)
  53   REAL(KIND=LONG) :: W1(JWORK)
  54
  55   INTEGER :: IFAIL = 0
  56   INTEGER :: IPIV(JMPRED)
  57   INTEGER :: NWORK = JWORK
  58
  59   INTEGER :: J, JJ, I, II
  60   REAL    :: RR
  61 !  REAL    :: offdiag
  62   REAL(KIND=LONG)  :: offdiag(NX)
  63   NXE=1
  64   NXEIG=NX
  65
  66   CCOF=0.
  67   XVEC=0.
  68      D=0.
  69
  70 ! CHECK INPUT PARAMETERS
  71   IF ((NX < 1) .OR. (NX > JMPRED)) STOP 701
  72   IF ((NXEIG < 1) .OR. (NXEIG > NX)) STOP 702
  73
  74
  75 ! PREDICTOR = XCOV, SIZE = NX*NX
  76 ! PREDICTANT = YCOV
  77 ! CROSS-COVARIANCE = XYCOV, SIZE = NX
  78
  79   DO I=1, NX
  80     XDEV(I) = SQRT(XCOV(I,I))
  81   ENDDO
  82
  83   PRINT *, 'PREDICTOR STD DEV AND MEAN NEW  F02FAF'
  84   PRINT 99, XDEV(1:NX)
  85   PRINT 99, XBAR(1:NX)
  86
  87 ! *** CALCULATE EIGENVECTOR REGRESSION COEFFICIENTS
  88
  89 !     CALCULATE EIGENVECTORS AND EIGENVALUES OF XCOV:
  90 !     XVEC AND XVAL.
  91
  92 ! PRINT *, 'PREDICTOR COVARIANCE MATRIX XCOV'
  93 ! DO J=1, NX
  94 !   PRINT 99, XCOV(1:NX,J)
  95 ! ENDDO
  96 ! DO J=1, NX
  97 !   DO I=1, NX
  98 !     XW1(I,J) = XCOV(I,J)/SQRT(XCOV(I,I)*XCOV(J,J))
  99 !   ENDDO
 100 ! ENDDO
 101 ! PRINT *, 'PREDICTOR CORRELATION MATRIX'
 102 ! DO J=1, NX
 103 !   PRINT 99, XW1(1:NX,J)
 104 ! ENDDO
 105
 106 ! CALL F02ABF(XCOV,NX,NX,XVAL,XVEC,NX,W1,IFAIL)
 107   XVEC=XCOV
 108 !  CALL F02FAF('V','L',NX,XVEC,NX,XVAL,W1,JWORK,IFAIL)
 109 ! tranform matrix to tridiagonal form
 110   call tred2(XVEC,NX,NX,XVAL,offdiag)
 111 ! QL decomposition
 112   call tqli(XVAL,offdiag,NX,NX,XVEC)
 113
 114   IF (IFAIL /= 0) THEN
 115     PRINT *, 'PROBLEM WITH PREDICTOR EIGENVALUE CALCULATION'
 116     STOP
 117   ENDIF
 118
 119   PRINT 89
 120   89 FORMAT (/1X,'XVAL = EIGENVALUES OF PREDICTOR in one and two zero')
 121   PRINT 99, XVAL(1:NX)
 122   99 FORMAT (1X,10F12.4)
 123
 124   PRINT *, 'EIGENVECTORS OF PREDICTOR BY COLUMN'
 125    XVEC(1:NX,1:(NXEIG-NXEIGS))=0
 126   DO I=1, NX
 127     PRINT 99, XVEC(I,1:NX)
 128
 129   ENDDO
 130
 131   XXBAR = MATMUL(XBAR,XVEC)
 132   PRINT *, 'MEAN VALUE OF PREDICTORS IN XVEC SPACE'
 133   PRINT 99, XXBAR(1:NX)
 134
 135   XYW1 = MATMUL(XYCOV,XVEC)
 136
 137 ! CALCULATE REGRESSION MATRIX OF EIGENVECTOR COEFFICIENTS:
 138
 139   D(NXE:NXEIG) = XYW1(NXE:NXEIG)/XVAL(NXE:NXEIG)
 140
 141 ! CALCULATE REGRESSION MATRIX OF Y AGAINST X :  CCOF :
 142 ! CCOF = XVEC * D
 143
 144   CCOF = MATMUL(XVEC(:,NXE:NXEIG),D(NXE:NXEIG))
 145
 146 ! ***     CALCULATE OFFSET VECTOR, CCOF0 :
 147 !         Y - YBAR  =  CCOF . (X-XBAR)
 148 !         THERFORE  Y = CCOF.X + CCOF0, WHERE CCOF0 = YBAR -CCOF.XBAR
 149
 150   CCOF0 = YBAR - SUM(CCOF(NXE:NX)*XBAR(NXE:NX))
 151
 152 ! PRINT 86
 153 ! 86 FORMAT (/1X,'PREDICTIVE MATRIX, CCOF')
 154 ! DO I=1, NX
 155 !   PRINT 99, CCOF(I)
 156 ! ENDDO
 157
 158 ! PRINT *, 'CCOF IN XVEC SPACE'
 159 ! DO I=1, NX
 160 !   PRINT 99, D(I)
 161 ! ENDDO
 162   XYW1=0
 163   XYW1(NXE:NXEIG) = MATMUL(XYCOV,XVEC(:,NXE:NXEIG))
 164   XYW1(NXE:NXEIG) = XYW1(1:NXEIG) * XYW1(NXE:NXEIG)
 165   XYW1(NXE:NXEIG) = XYW1(1:NXEIG)/XVAL(NXE:NXEIG)
 166
 167   PRINT *, 'EXPLANATION OF VARIANCE'
 168   PRINT 99, YCOV
 169   PRINT *, ' '
 170   PRINT 99, 1.0
 171   DO I=1, NXEIG
 172     PRINT 99, XYW1(I)/YCOV
 173   ENDDO
 174   PRINT 99, (YCOV - SUM(XYW1(1:NXEIG)))/YCOV
 175   PRINT *, ' '
 176   PRINT 99, YCOV - SUM(XYW1(1:NXEIG))
 177
 178 ! PRINT 85
 179 ! 85 FORMAT (/1X,'OFFSET VECTOR, CCOF0')
 180 ! PRINT 99, CCOF0
 181
 182 ! *** CALCULATE RESIDUAL COVARIANCE MATRIX
 183 !     = YCOV + CCOF.XCOV.CCOF(TR) - 2.XYCOF.CCOF(TR)
 184
 185   YXW1 = CCOF
 186   XYW1 = MATMUL(CCOF,XCOV)
 187   XYW1 = XYW1 - 2*XYCOV
 188   RESCOV = YCOV + DOT_PRODUCT(XYW1,YXW1)
 189
 190 ! RESCOV NOW CONTAINS RESIDUAL COVARIANCE(CORRELATION) MATRIX
 191   RR = RESCOV
 192   IF (RR < 0.0) THEN
 193     PRINT 51, I, RR
 194     51 FORMAT (/1X,'RES ERR COV FOR ELEMENT',I5,' = ',E12.4/1X,'SET TO ZERO')
 195     RR=0.0
 196   ENDIF
 197   RESERR = SQRT(RR)
 198
 199 ! PRINT 84
 200 ! 84 FORMAT (/1X,'RESIDUAL ERROR VECTOR, RESERR')
 201 ! PRINT 99, RESERR
 202
 203
 204   RETURN
 205   END SUBROUTINE REGRESS_ONE