Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/polrg.cdk
There are 2 other files named polrg.cdk in the archive. Click here to see a list.
$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