regress_one.f90
References to this file elsewhere.
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