C PLRG 10
C ..................................................................PLRG 20
C PLRG 30
C SAMPLE MAIN PROGRAM FOR POLYNOMIAL REGRESSION - POLRG PLRG 40
C PLRG 50
C PURPOSE PLRG 60
C (1) READ THE PROBLEM PARAMETER CARD FOR A POLYNOMIAL REGRES-PLRG 70
C SION, (2) CALL SUBROUTINES TO PERFORM THE ANALYSIS, (3) PLRG 80
C PRINT THE REGRESSION COEFFICIENTS AND ANALYSIS OF VARIANCE PLRG 90
C TABLE FOR POLYNOMIALS OF SUCCESSIVELY INCREASING DEGREES, PLRG 100
C AND (4) OPTIONALLY PRINT THE TABLE OF RESIDUALS AND A PLOT PLRG 110
C OF Y VALUES AND Y ESTIMATES. PLRG 120
C PLRG 130
C REMARKS PLRG 140
C THE NUMBER OF OBSERVATIONS, N, MUST BE GREATER THAN M+1, PLRG 150
C WHERE M IS THE HIGHEST DEGREE POLYNOMIAL SPECIFIED. PLRG 160
C IF THERE IS NO REDUCTION IN THE RESIDUAL SUM OF SQUARES PLRG 170
C BETWEEN TWO SUCCESSIVE DEGREES OF THE POLYNOMIALS, THE PLRG 180
C PROGRAM TERMINATES THE PROBLEM BEFORE COMPLETING THE ANALY- PLRG 190
C SIS FOR THE HIGHEST DEGREE POLYNOMIAL SPECIFIED. PLRG 200
C PLRG 210
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED PLRG 220
C GDATA PLRG 230
C ORDER PLRG 240
C MINV PLRG 250
C MULTR PLRG 260
C PLOT (A SPECIAL PLOT SUBROUTINE PROVIDED FOR THE SAMPLE PLRG 270
C PROGRAM.) PLRG 280
C PLRG 290
C METHOD PLRG 300
C REFER TO B. OSTLE, 'STATISTICS IN RESEARCH', THE IOWA STATE PLRG 310
C COLLEGE PRESS', 1954, CHAPTER 6. PLRG 320
C PLRG 330
C ..................................................................PLRG 340
C PLRG 350
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE PLRG 360
C PRODUCT OF N*(M+1), WHERE N IS THE NUMBER OF OBSERVATIONS AND M PLRG 370
C IS THE HIGHEST DEGREE POLYNOMIAL SPECIFIED.. PLRG 380
C PLRG 390
DIMENSION X(1100) PLRG 400
C PLRG 410
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE PLRG 420
C PRODUCT OF M*M.. PLRG 430
C PLRG 440
DIMENSION DI(100) PLRG 450
C PLRG 460
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO PLRG 470
C (M+2)*(M+1)/2.. PLRG 480
C PLRG 490
DIMENSION D(66) PLRG 500
C PLRG 510
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO M.. PLRG 520
C PLRG 530
DIMENSION B(10),E(10),SB(10),T(10) PLRG 540
C PLRG 550
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO (M+1).. PLRG 560
C PLRG 570
DIMENSION XBAR(11),STD(11),COE(11),SUMSQ(11),ISAVE(11) PLRG 580
C PLRG 590
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO 10.. PLRG 600
C PLRG 610
DIMENSION ANS(10) PLRG 620
C PLRG 630
C THE FOLLOWING DIMENSION WILL BE USED IF THE PLOT OF OBSERVED DATA PLRG 640
C AND ESTIMATES IS DESIRED. THE SIZE OF THE DIMENSION, IN THIS PLRG 650
C CASE, MUST BE GREATER THAN OR EQUAL TO N*3. OTHERWISE, THE SIZE PLRG 660
C OF DIMENSION MAY BE SET TO 1. PLRG 670
C PLRG 680
DIMENSION P(300) PLRG 690
C PLRG 700
C ..................................................................PLRG 710
C PLRG 720
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE PLRG 730
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION PLRG 740
C STATEMENT WHICH FOLLOWS. PLRG 750
C PLRG 760
C DOUBLE PRECISION X,XBAR,STD,D,SUMSQ,DI,E,B,SB,T,ANS,DET,COE PLRG 770
C PLRG 780
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS PLRG 790
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS PLRG 800
C ROUTINE. PLRG 810
C PLRG 820
C ...............................................................PLRG 830
C PLRG 840
1 FORMAT(A4,A2,I5,I2,I1) PLRG 850
2 FORMAT(2F6.0) PLRG 860
3 FORMAT(27H1POLYNOMIAL REGRESSION.....,A4,A2/) PLRG 870
4 FORMAT(23H0NUMBER OF OBSERVATIONS,I6//) PLRG 880
5 FORMAT(32H0POLYNOMIAL REGRESSION OF DEGREE,I3) PLRG 890
6 FORMAT(12H0 INTERCEPT,E20.7) PLRG 900
7 FORMAT(26H0 REGRESSION COEFFICIENTS/(6E20.7)) PLRG 910
8 FORMAT(1H0/24X,24HANALYSIS OF VARIANCE FOR,I4,19H DEGREE POLYNOMIPLRG 920
1AL/) PLRG 930
9 FORMAT(1H0,5X,19HSOURCE OF VARIATION,7X,9HDEGREE OF,7X,6HSUM OF,9XPLRG 940
1,4HMEAN,10X,1HF,9X,20HIMPROVEMENT IN TERMS/33X,7HFREEDOM,8X,7HSQUAPLRG 950
2RES,7X,6HSQUARE,7X,5HVALUE,8X,17HOF SUM OF SQUARES) PLRG 960
10 FORMAT(20H0 DUE TO REGRESSION,12X,I6,F17.5,F14.5,F13.5,F20.5) PLRG 970
11 FORMAT(32H DEVIATION ABOUT REGRESSION ,I6,F17.5,F14.5) PLRG 980
12 FORMAT(8X,5HTOTAL,19X,I6,F17.5///) PLRG 990
13 FORMAT(17H0 NO IMPROVEMENT) PLRG1000
14 FORMAT(1H0//27X,18HTABLE OF RESIDUALS//16H OBSERVATION NO.,5X,7HX PLRG1010
1VALUE,7X,7HY VALUE,7X,10HY ESTIMATE,7X,8HRESIDUAL/) PLRG1020
15 FORMAT(1H0,3X,I6,F18.5,F14.5,F17.5,F15.5) PLRG1030
C PLRG1040
C ..................................................................PLRG1050
C PLRG1060
C READ PROBLEM PARAMETER CARD PLRG1070
C PLRG1080
100 READ (5,1,END=999) PR,PR1,N,M,NPLOT PLRG1090
C PLRG1100
C PR....PROBLEM NUMBER (MAY BE ALPHAMERIC) PLRG1110
C PR1...PROBLEM NUMBER (CONTINUED) PLRG1120
C N.....NUMBER OF OBSERVATIONS PLRG1130
C M.....HIGHEST DEGREE POLYNOMIAL SPECIFIED PLRG1140
C NPLOT.OPTION CODE FOR PLOTTING PLRG1150
C 0 IF PLOT IS NOT DESIRED. PLRG1160
C 1 IF PLOT IS DESIRED. PLRG1170
C PLRG1180
C PRINT PROBLEM NUMBER AND N. PLRG1190
C PLRG1200
WRITE (6,3) PR,PR1 PLRG1210
WRITE (6,4) N PLRG1220
C PLRG1230
C READ INPUT DATA PLRG1240
C PLRG1250
L=N*M PLRG1260
DO 110 I=1,N PLRG1270
J=L+I PLRG1280
C PLRG1290
C X(I) IS THE INDEPENDENT VARIABLE, AND X(J) IS THE DEPENDENT PLRG1300
C VARIABLE. PLRG1310
C PLRG1320
110 READ (5,2) X(I),X(J) PLRG1330
C PLRG1340
CALL GDATA (N,M,X,XBAR,STD,D,SUMSQ) PLRG1350
C PLRG1360
MM=M+1 PLRG1370
SUM=0.0 PLRG1380
NT=N-1 PLRG1390
C PLRG1400
DO 200 I=1,M PLRG1410
ISAVE(I)=I PLRG1420
C PLRG1430
C FORM SUBSET OF CORRELATION COEFFICIENT MATRIX PLRG1440
C PLRG1450
CALL ORDER (MM,D,MM,I,ISAVE,DI,E) PLRG1460
C PLRG1470
C INVERT THE SUBMATRIX OF CORRELATION COEFFICIENTS PLRG1480
C PLRG1490
CALL MINV (DI,I,DET,B,T) PLRG1500
C PLRG1510
CALL MULTR (N,I,XBAR,STD,SUMSQ,DI,E,ISAVE,B,SB,T,ANS) PLRG1520
C PLRG1530
C PRINT THE RESULT OF CALCULATION PLRG1540
C PLRG1550
WRITE (6,5) I PLRG1560
IF(ANS(7)) 140,130,130 PLRG1570
130 SUMIP=ANS(4)-SUM PLRG1580
IF(SUMIP) 140, 140, 150 PLRG1590
140 WRITE (6,13) PLRG1600
GO TO 210 PLRG1610
150 WRITE (6,6) ANS(1) PLRG1620
WRITE (6,7) (B(J),J=1,I) PLRG1630
WRITE (6,8) I PLRG1640
WRITE (6,9) PLRG1650
SUM=ANS(4) PLRG1660
WRITE (6,10) I,ANS(4),ANS(6),ANS(10),SUMIP PLRG1670
NI=ANS(8) PLRG1680
WRITE (6,11) NI,ANS(7),ANS(9) PLRG1690
WRITE (6,12) NT,SUMSQ(MM) PLRG1700
C PLRG1710
C SAVE COEFFICIENTS FOR CALCULATION OF Y ESTIMATES PLRG1720
C PLRG1730
COE(1)=ANS(1) PLRG1740
DO 160 J=1,I PLRG1750
160 COE(J+1)=B(J) PLRG1760
LA=I PLRG1770
200 CONTINUE PLRG1780
C PLRG1790
C TEST WHETHER PLOT IS DESIRED PLRG1800
C PLRG1810
210 IF(NPLOT) 100, 100, 220 PLRG1820
C PLRG1830
C CALCULATE ESTIMATES PLRG1840
C PLRG1850
220 NP3=N+N PLRG1860
DO 230 I=1,N PLRG1870
NP3=NP3+1 PLRG1880
P(NP3)=COE(1) PLRG1890
L=I PLRG1900
DO 230 J=1,LA PLRG1910
P(NP3)=P(NP3)+X(L)*COE(J+1) PLRG1920
230 L=L+N PLRG1930
C PLRG1940
C COPY OBSERVED DATA PLRG1950
C PLRG1960
N2=N PLRG1970
L=N*M PLRG1980
DO 240 I=1,N PLRG1990
P(I)=X(I) PLRG2000
N2=N2+1 PLRG2010
L=L+1 PLRG2020
240 P(N2)=X(L) PLRG2030
C PLRG2040
C PRINT TABLE OF RESIDUALS PLRG2050
C PLRG2060
WRITE (6,3) PR,PR1 PLRG2070
WRITE (6,5) LA PLRG2080
WRITE (6,14) PLRG2090
NP2=N PLRG2100
NP3=N+N PLRG2110
DO 250 I=1,N PLRG2120
NP2=NP2+1 PLRG2130
NP3=NP3+1 PLRG2140
RESID=P(NP2)-P(NP3) PLRG2150
250 WRITE (6,15) I,P(I),P(NP2),P(NP3),RESID PLRG2160
C PLRG2170
CALL PLOT (LA,P,N,3,0,1) PLRG2180
C PLRG2190
GO TO 100 PLRG2200
999 STOP
END PLRG2210