File: POLRG.FT of Tape: Various/ETH/eth11-1
(Source file text)
C SAMPLE MAIN PROGRAM FOR POLYNOMIAL REGRESSION - POLRG
C THE FOLLOWING ROUTINES ARE USED: GDATA,ORDER,MINV,MULTR
C AND AN OPTIONAL PLOTTING SUBROUTINE COULD BE ADDED
C
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE
C PRODUCT OF N*(M+1) WHERE N IS THE NUMBER OF OBSERVATIONS AND
C M IS THE HIGHEST DEGREE POLYNOMIAL SPECIFIED.
DIMENSION X(75)
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE
C PRODUCT OF M*M.
DIMENSION DI(16)
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO
C (M+2)*(M+1)/2.
DIMENSION D(15)
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO M.
DIMENSION B(4),SB(4),T(4),E(4)
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO
C (M+1).
DIMENSION XBAR(5),STD(5),COE(5),SUMSQ(5),ISAVE(5)
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO 10.
DIMENSION ANS(10)
C THE FOLLOWING DIMENSION WILL BE USED IF THE PLOT OF OBSERVED
C DATA AND ESTIMATES IS DESIRED. THE SIZE OF THE DIMENSION, IN
C THIS CASE, MUST BE GREATER THAN OR EQUAL TO N*3. OTHERWISE,
C THE SIZE OF THE DIMENSION MAY BE SET TO 1.
DIMENSION P(45)
COMMON IOUT,IN
C
C OUTPUT CHANNEL = IOUT, INPUT CHANNEL = IN
1 FORMAT(A4,A2,I5,I2,I1)
2 FORMAT(2F6.0)
3 FORMAT(//27H POLYNOMIAL REGRESSION.....,A4,A2)
4 FORMAT(/23H NUMBER OF OBSERVATION,I6)
5 FORMAT(/32H POLYNOMIAL REGRESSION OF DEGREE,I3)
6 FORMAT(/12H INTERCEPT,F15.5)
7 FORMAT(/26H REGRESSION COEFFICIENTS/(6F12.5))
8 FORMAT(/8X,24HANALYSIS OF VARIANCE FOR,I4,19H DEGREE POLYNOMI
1AL/)
9 FORMAT(/' SOURCE OF VARIATION DEGREE'6X'SUM OF'8X'MEAN'
110X'F IMPROVEMENT'/24X'OF'8X'SQUARES'7X'SQUARE'6X
2'VALUE IN TERMS OF'/21X'FREEDOM'39X'SUM OF SQUARES')
10 FORMAT(/' DUE TO REGRESSION',I6,4X,4F13.5)
11 FORMAT(' DEVIATION ABOUT'/8X'REGRESSION',I6,4X,2F13.5)
12 FORMAT(' TOTAL',12X,I6,4X,F13.5///)
13 FORMAT(/17H NO IMPROVEMENT )
14 FORMAT(//27X,18HTABLE OF RESIDUALS//16H OBSERVATION NO.,5X
1,7HX VALUE,7X,7HY VALUE,7X,10HY ESTIMATE,7X,8HRESIDUAL/)
15 FORMAT(/,3X,I6,F18.5,F14.5,F17.5,F15.5)
C
IN=1
IOUT=2
C READ PARAMETER CARD
100 READ(IN,1) PR,PR1,N,M,NPLOT
C PR=PROBLEM NUMBER (MAY BE ALPHAMERIC)
C PR1=PROBLEM NUMBER (CONTINUED)
C N=NUMBER OF OBSERVATIONS
C M=HIGHEST DEGREE POLYNOMIAL SPECIFIED
C NPLOT=OPTION CODE FOR PLOTTING
C 0 IF PLOT IS NOT DESIRED
C 1 IF PLOT IS DESIRED
IF(N.EQ.0) STOP
C PRINT PROBLEM NUMBER AND N.
WRITE(IOUT,3) PR,PR1
WRITE(IOUT,4) N
C READ INPUT DATA
L=N*M
DO 110 I=1,N
J=L+I
C X(I) IS THE INDEPENDENT VARIABLE, AND X(J) IS THE DEPENDENT
C VARIABLE.
110 READ(IN,2) X(I),X(J)
CALL GDATA(N,M,X,XBAR,STD,D,SUMSQ)
MM=M+1
SUM=0.0
NT=N-1
DO 200 I=1,M
ISAVE(I)=I
C FORM SUBSET OF CORRELATION COEFFICIENT MATRIX
CALL ORDER(MM,D,MM,I,ISAVE,DI,E)
C INVERT THE SUBMATRIX OF CORRELATION COEFFICIENTS
CALL MINV(DI,I,DET,B,T)
CALL MULTR(N,I,XBAR,STD,SUMSQ,DI,E,ISAVE,B,SB,T,ANS)
C PRINT THE RESULT OF CALCULATION
WRITE(IOUT,5) I
IF(ANS(7)) 140,130,130
130 SUMIP=ANS(4)-SUM
IF(SUMIP) 140,140,150
140 WRITE(IOUT,13)
GO TO 210
150 WRITE(IOUT,6) ANS(1)
WRITE(IOUT,7) (B(J),J=1,I)
WRITE(IOUT,8) I
WRITE(IOUT,9)
SUM=ANS(4)
WRITE(IOUT,10) I,ANS(4),ANS(6),ANS(10),SUMIP
NI=ANS(8)
WRITE(IOUT,11) NI,ANS(7),ANS(9)
WRITE(IOUT,12) NT,SUMSQ(MM)
C SAVE COEFFICIENTS FOR CALCULATION OF Y ESTIMATES
COE(1)=ANS(1)
DO 160 J=1,I
160 COE(J+1)=B(J)
LA=I
200 CONTINUE
C TEST WHETHER PLOT IS DESIRED
210 IF(NPLOT) 100,100,220
C CALCULATE ESTIMATES
220 NP3=N+N
DO 230 I=1,N
NP3=NP3+1
P(NP3)=COE(1)
L=I
DO 230 J=1,LA
P(NP3)=P(NP3)+X(L)*COE(J+1)
230 L=L+N
C COPY OBSERVED DATA
N2=N
L=N*M
DO 240 I=1,N
P(I)=X(I)
N2=N2+1
L=L+1
240 P(N2)=X(L)
C PRINT TABLE OF RESIDUALS
WRITE(IOUT,3) PR,PR1
WRITE(IOUT,5) LA
WRITE(IOUT,14)
NP2=N
NP3=N+N
DO 250 I=1,N
NP2=NP2+1
NP3=NP3+1
RESID=P(NP2)-P(NP3)
250 WRITE(IOUT,15) I,P(I),P(NP2),P(NP3),RESID
CALL PLOT(LA,P,N,3,0,80,1)
GO TO 100
END