$JOB POLRG[30,30]
$FORTRAN POLRG
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
$FORTRAN PLOT
C PLOT 10
C ..................................................................PLOT 20
C PLOT 30
C SUBROUTINE PLOT PLOT 40
C PLOT 50
C PURPOSE PLOT 60
C PLOT SEVERAL CROSS-VARIABLES VERSUS A BASE VARIABLE PLOT 70
C PLOT 80
C USAGE PLOT 90
C CALL PLOT (NO,A,N,M,NL,NS) PLOT 100
C PLOT 110
C DESCRIPTION OF PARAMETERS PLOT 120
C NO - CHART NUMBER (3 DIGITS MAXIMUM) PLOT 130
C A - MATRIX OF DATA TO BE PLOTTED. FIRST COLUMN REPRESENTS PLOT 140
C BASE VARIABLE AND SUCCESSIVE COLUMNS ARE THE CROSS- PLOT 150
C VARIABLES (MAXIMUM IS 9). PLOT 160
C N - NUMBER OF ROWS IN MATRIX A PLOT 170
C M - NUMBER OF COLUMNS IN MATRIX A (EQUAL TO THE TOTAL PLOT 180
C NUMBER OF VARIABLES). MAXIMUM IS 10. PLOT 190
C NL - NUMBER OF LINES IN THE PLOT. IF 0 IS SPECIFIED, 50 PLOT 200
C LINES ARE USED. PLOT 210
C NS - CODE FOR SORTING THE BASE VARIABLE DATA IN ASCENDING PLOT 220
C ORDER PLOT 230
C 0 SORTING IS NOT NECESSARY (ALREADY IN ASCENDING PLOT 240
C ORDER). PLOT 250
C 1 SORTING IS NECESSARY. PLOT 260
C PLOT 270
C REMARKS PLOT 280
C NONE PLOT 290
C PLOT 300
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED PLOT 310
C NONE PLOT 320
C PLOT 330
C ..................................................................PLOT 340
C PLOT 350
SUBROUTINE PLOT(NO,A,N,M,NL,NS) PLOT 360
DIMENSION OUT(101),YPR(11),ANG(9),A(1) PLOT 370
C PLOT 380
1 FORMAT(1H1,60X,7H CHART ,I3,//) PLOT 390
2 FORMAT(1H ,F11.4,5X,101A1) PLOT 400
3 FORMAT(1H ) PLOT 410
4 FORMAT(10H 123456789) PLOT 420
5 FORMAT(10A1) PLOT 430
7 FORMAT(1H ,16X,101H. . . . . PLOT 440
1 . . . . . .) PLOT 450
8 FORMAT(1H0,9X,11F10.4) PLOT 460
C PLOT 470
C ..................................................................PLOT 480
C PLOT 490
NLL=NL PLOT 500
C PLOT 510
IF(NS) 16, 16, 10 PLOT 520
C PLOT 530
C SORT BASE VARIABLE DATA IN ASCENDING ORDER PLOT 540
C PLOT 550
10 DO 15 I=1,N PLOT 560
DO 14 J=I,N PLOT 570
IF(A(I)-A(J)) 14, 14, 11 PLOT 580
11 L=I-N PLOT 590
LL=J-N PLOT 600
DO 12 K=1,M PLOT 610
L=L+N PLOT 620
LL=LL+N PLOT 630
F=A(L) PLOT 640
A(L)=A(LL) PLOT 650
12 A(LL)=F PLOT 660
14 CONTINUE PLOT 670
15 CONTINUE PLOT 680
C PLOT 690
C TEST NLL PLOT 700
C PLOT 710
16 IF(NLL) 20, 18, 20 PLOT 720
18 NLL=50 PLOT 730
C PLOT 740
C PRINT TITLE PLOT 750
C PLOT 760
20 WRITE(6,1)NO PLOT 770
C PLOT 780
C DEVELOP BLANK AND DIGITS FOR PRINTING PLOT 790
C PLOT 800
REWIND 13 PLOT 810
WRITE (13,4) PLOT 820
REWIND 13 PLOT 830
READ (13,5) BLANK,(ANG(I),I=1,9) PLOT 840
REWIND 13 PLOT 850
C PLOT 860
C FIND SCALE FOR BASE VARIABLE PLOT 870
C PLOT 880
XSCAL=(A(N)-A(1))/(FLOAT(NLL-1)) PLOT 890
C PLOT 900
C FIND SCALE FOR CROSS-VARIABLES PLOT 910
C PLOT 920
M1=N+1 PLOT 930
YMIN=A(M1) PLOT 940
YMAX=YMIN PLOT 950
M2=M*N PLOT 960
DO 40 J=M1,M2 PLOT 970
IF(A(J)-YMIN) 28,26,26 PLOT 980
26 IF(A(J)-YMAX) 40,40,30 PLOT 990
28 YMIN=A(J) PLOT1000
GO TO 40 PLOT1010
30 YMAX=A(J) PLOT1020
40 CONTINUE PLOT1030
YSCAL=(YMAX-YMIN)/100.0 PLOT1040
C PLOT1050
C FIND BASE VARIABLE PRINT POSITION PLOT1060
C PLOT1070
XB=A(1) PLOT1080
L=1 PLOT1090
MY=M-1 PLOT1100
I=1 PLOT1110
45 F=I-1 PLOT1120
XPR=XB+F*XSCAL PLOT1130
IF(A(L)-XPR) 50,50,70 PLOT1140
C PLOT1150
C FIND CROSS-VARIABLES PLOT1160
C PLOT1170
50 DO 55 IX=1,101 PLOT1180
55 OUT(IX)=BLANK PLOT1190
DO 60 J=1,MY PLOT1200
LL=L+J*N PLOT1210
JP=((A(LL)-YMIN)/YSCAL)+1.0 PLOT1220
OUT(JP)=ANG(J) PLOT1230
60 CONTINUE PLOT1240
C PLOT1250
C PRINT LINE AND CLEAR, OR SKIP PLOT1260
C PLOT1270
WRITE(6,2)XPR,(OUT(IZ),IZ=1,101) PLOT1280
L=L+1 PLOT1290
GO TO 80 PLOT1300
70 WRITE(6,3) PLOT1310
80 I=I+1 PLOT1320
IF(I-NLL) 45, 84, 86 PLOT1330
84 XPR=A(N) PLOT1340
GO TO 50 PLOT1350
C PLOT1360
C PRINT CROSS-VARIABLES NUMBERS PLOT1370
C PLOT1380
86 WRITE(6,7) PLOT1390
YPR(1)=YMIN PLOT1400
DO 90 KN=1,9 PLOT1410
90 YPR(KN+1)=YPR(KN)+YSCAL*10.0 PLOT1420
YPR(11)=YMAX PLOT1430
WRITE(6,8)(YPR(IP),IP=1,11) PLOT1440
RETURN PLOT1450
END PLOT1460
$DECK POL.CDR
SAMPLE00015041 20
1 10 30
2 16 40
3 20 50
4 23 60
5 25 70
6 26 80
7 30 90
8 36 100
9 48 110
10 62 120
11 78 130
12 94 140
13 107 150
14 118 160
15 127 170
$EOD
.ASSIGN CDR 5
.ASSIGN LPT 6
.SET CDR POL
.EXECUTE/REL POLRG,PLOT,WES:SSP/LIB
%FIN::
.DELETE POL.CDR