C             POLYNOMIAL REGRESSION                   JUNE 10, 1966
C        THIS IS A SIFTED VERSION OF BMD05R ORIGINALLY WRITTEN IN
C        FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C        AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
      DIMENSION SE(10),XX(500),YY(500),FMT(180),XXB(10),KODE(8),CONST(8)
      DIMENSION X(400,6),SD(11,11),AV(11),B(10),RR(11,11)
      DOUBLE PRECISION A1,A2,A3,SUM,X,SD,AV,B,RR,ONO,GSSDR,SSDR,SSAR,SU,
     1A
      DATA A1,A2,A3/6HPROBLM,6HFINISH,6HSPECTG/
C
  918 FORMAT(45H1BMD05R - POLYNOMIAL REGRESSION - VERSION OF ,
     118HJUNE 10, 1966     ,/
     241H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
C
      NTAPE=5
	CALL USAGEB('BMD05R')
    5 READ (5,900)SUM,PROB,N,NP,IPLOT,NVG,NACC,MTAPE,KVR
      NPP=0
      IF(SUM.EQ.A1)GO TO 4
      IF(SUM.EQ.A2)GO TO 15
 8    WRITE (6,919)
 15   IF(6-NTAPE)16,4031,4031
 16   REWIND NTAPE
      GO TO 4031
 4    CALL TPWD(MTAPE,NTAPE)
      IF(NACC)748,748,749
 748  NACC=13
 749  TOL=10.0**(-NACC)
 19   WRITE (6,918)
      WRITE (6,902)PROB
      WRITE (6,903)N
      WRITE(6,9181) KVR
      IF(NP*(NP-11))200,401,401
  200 IF(N-501)204,402,402
  204 IF(N-NP)403,403,205
 205  IERROR=0
      NT=NP+1
      ONO=0.0D0
      ON=N
      ONO=ON
      DO 721 I=1,500
      X(I,NT)=0.0D0
 721  X(I,1)=0.0D0
      DO 20 I=1,NT
      AV(I)=0.0D0
      DO 20 J=1,NT
 20   SD(I,J)=0.0D0
      IF(KVR.GT.0.AND.KVR.LE.10)GO TO 21
      KVR=1
      WRITE (6,4000)
 21   KVR=KVR*18
      READ (5,916)(FMT(I),I=1,KVR)
      WRITE(6,211) (FMT(I), I=1,KVR)
 211  FORMAT(' THE VARIABLE FORMAT IS  ',18A4,/(1H0,24X,18A4))
      IF(-NVG)22,23,23
 22   WRITE (6,901)
      READ(5,920)SUM,M,(KODE(I),CONST(I),I=1,M)
      IF(SUM.EQ.A3)GO TO 220
 210  WRITE (6,921)
      IERROR=-501
      GO TO 23
 220  DO 230 I=1,M
      WRITE (6,922)KODE(I),CONST(I)
      IF(KODE(I)*(KODE(I)-11))230,225,225
 225  WRITE (6,923)KODE(I)
 230  CONTINUE
 23   READ (NTAPE,FMT)(X(I,1),X(I,NT),I=1,N)
      IF(NTAPE-5)10,11,10
 10   REWIND NTAPE
 11   IF(NVG)25,25,24
 24   CALL TRNGEN(N,NT,M,X,KODE,CONST)
      IF(IERROR) 5, 25, 25
 25   SUMY=0.0
      SUMYY=0.0
      SUMXY=0.0
      SUMX=0.0
      SUMXX=0.0
      DO 27 I=1,N
      YY(I)=X(I,NT)
      XX(I)=X(I,1)
      SUMX=SUMX+XX(I)
      SUMY=SUMY+YY(I)
      SUMXY=SUMXY+XX(I)*YY(I)
      SUMXX=SUMXX+XX(I)*XX(I)
 27   SUMYY=SUMYY+YY(I)*YY(I)
      R=(SUMXY-SUMX*SUMY/ON)/SQRT((SUMXX-SUMX*SUMX/ON)*(SUMYY-SUMY*SUMY/
     1ON))
      RRR=R*R
      IF(NP-1)1010,1010,1212
 1212 DO 30 J=2,NP
      JJ=J-1
      DO 30 I=1,N
 30   X(I,J)=X(I,JJ)*X(I,1)
 1010 DO 36 J=1,NT
      DO 35 I=1,N
 35   AV(J)=AV(J)+X(I,J)
 36   AV(J)=AV(J)/ONO
      DO 40 J=1,NT
      DO 40 I=1,N
 40   X(I,J)=X(I,J)-AV(J)
      DO 46 J=1,NT
      DO 46 K=J,NT
      DO 45 I=1,N
 45   SD(J,K)=SD(J,K)+(X(I,J)*X(I,K))
 46   SD(K,J)=SD(J,K)
      B(1)=SD(1,NT)/SD(1,1)
      A=AV(NT)-B(1)*AV(1)
      SE(1)=DSQRT((SUMYY-A*SUMY-B(1)*SUMXY)/(SD(1,1)*(ON-2.0)))
      WRITE (6,904)
      WRITE (6,905)
      WRITE (6,906)AV(1)
      WRITE (6,907)AV(NT)
      WRITE (6,908)A
      WRITE (6,909)B(1)
      WRITE (6,910)SE(1)
      WRITE (6,911)R
      SSAR=RRR*SD(NT,NT)
      SSDR=(1.0-RRR)*SD(NT,NT)
      N1=1
      CALL REPORT (SSAR,SSDR,N1,N,NPP)
      GSSDR=SSDR
      IF(NP-1)50,50,1414
   50 XXB(1)=B(1)
      AXX=SNGL(A)
      NXX=1
      GO TO 1313
 1414 DO 73 IJ=2,NP
      N1=IJ
      DO 52 J=1,IJ
       B(J)=SD(J,NT)
       DO 52 I=1,IJ
 52    RR(I,J)=SD(I,J)
      CALL MATINV(RR,IJ,B,D,T)
      IF(T-TOL)66,191,191
 66   WRITE (6,917)
      SSDR=GSSDR
      IF(NPP) 67, 67, 75
 67   NPP=1
      GO TO 75
 191  DO 192 I=1,IJ
      RR(I,I)=-RR(I,I)
 192  B(I)=-B(I)
      SSAR=0.0D0
      DO 60 I=1,IJ
 60   SSAR=SSAR+B(I)*SD(I,NT)
      SSDR=SD(NT,NT)-SSAR
      DF=IJ+1
      VARI=SNGL(SSDR)/(ON-DF)
      DO 62 J=1,IJ
      SUMY=VARI*RR(J,J)
 62   SE(J)=SQRT(SUMY)
      SU=0.0D0
      DO 64 I=1,IJ
 64   SU=SU+B(I)*AV(I)
 69   GSSDR=SSDR
      A=AV(NT)-SU
      WRITE (6,912)N1
      WRITE (6,908)A
      WRITE (6,913)
      WRITE (6,915)(B(I), I=1,N1)
      DO 70 I=1,N1
   70 XXB(I)=B(I)
      AXX=SNGL(A)
      NXX=N1
      WRITE (6,914)
      WRITE (6,915)(SE(I), I=1,N1)
      WRITE (6,904)
 73   CALL REPORT (SSAR,SSDR,N1,N,NPP)
 75   N1=105
      WRITE (6,904)
      CALL REPORT (SSAR,SSDR,N1,N,NPP)
 1313 IF(IPLOT-1)5,600,5
 600  CALL GRAPH (XX,YY,AXX,XXB,N,NXX)
      GO TO 5
  401 PRINT 4001
      GO TO 4031
  402 PRINT 4002
      GO TO 4031
  403 PRINT 4003
 4031 STOP
 900  FORMAT(A6,A2,I3,4I2,49X,2I2)
 901  FORMAT(21H0TRANSGENERATION CARD//16H CODE   CONSTANT//)
 902  FORMAT(15H0PROBLEM NO.   A2)
 903  FORMAT(12H SAMPLE SIZEI5)
 904  FORMAT(1H0)
 905  FORMAT(38H0REGRESSION - ONE INDEPENDENT VARIABLE)
 906  FORMAT(23H0XMEAN.....            F12.5)
 907  FORMAT(23H YMEAN.....            F12.5)
 908  FORMAT(23H0INTERCEPT (A VALUE)...F12.5)
 909  FORMAT(23H REG. COEFFICIENT......F12.5)
 910  FORMAT(25H STD. ERROR OF REG. COEF.F10.5)
 911  FORMAT(23H0CORRELATION COEF......F12.5///)
 912  FORMAT(32H0POLYNOMIAL REGRESSION OF DEGREEI3)
 913  FORMAT(24H0REGRESSION COEFFICIENTS)
 914  FORMAT(25H0STD. ERROR OF REG. COEF.)
  915 FORMAT(10(1X,F12.5))
  916 FORMAT(20A4)
 917  FORMAT(71H0COMPUTATIONAL ACCURACY NOT SUFFICIENT FOR POLYNOMIALS O
     1F HIGHER DEGREE)
 9181 FORMAT(' NO. VF CDS.',I5)
 919  FORMAT(22H0ERROR ON PROBLEM CARD)
 920  FORMAT(A6,I1,8(I2,F6.0))
 921  FORMAT(96H0SPECTG CARD MISPUNCHED OR OUT OF ORDER. PROGRAM WILL PR
     1OCEED TO NEXT JOB, IF ANY, OR TERMINATE.)
 922  FORMAT(1H I3,F12.5)
 923  FORMAT(29H0ILLEGAL TRANSGENERATION CODEI3,67H SPECIFIED. PROGRAM W
     1ILL CONTINUE LEAVING OUT THIS TRANSGENERATION.)
 4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
     1IED, ASSUMED TO BE 1.)
 4001 FORMAT(' DEGREE OF POLYNOMIAL INCORRECTLY SPECIFIED.  IT CANNOT EX
     XCEED 10')
 4002 FORMAT(' NUMBER OF CASES TOO LARGE.  IT CANNOT EXCEED 500.')
 4003 FORMAT(' DEGREE OF POLYNOMIAL IS TOO HIGH.')
      END
      SUBROUTINE FORM2(SYMB,XY)
C        SUBROUTINE FORM2 FOR PLOTR (IBM 360)         JUNE 21, 1966
      DIMENSION TEST(18)
      INTEGER XY,SYMB,BLANK,TEST
      DATA BLANK/2H  /
      DATA TEST/'B ','B ','4 ','5 ','6 ','7 ','8 ','9 ','A ','B ','C ',
     1'D ','E ','F ','G ','H ','I ','/ '/
      IF(XY.EQ.BLANK)GO TO 50
      DO 30 I=1,17
      IF(XY.NE.TEST(I))GO TO 30
C        PUT IN NEXT SYMBOL OF ARRAY FOR MULTIPLE POINTS
      XY=TEST(I+1)
      GO TO 100
   30 CONTINUE
      IF(XY.EQ.TEST(18))GO TO 100
C        IF OTHER THAN CHARACTERS IN ARRAY TEST PUT IN CHARACTER 2.
      XY=TEST(1)
      GO TO 100
C        IF BLANK, PUT IN SYMBOL
   50 XY=SYMB
  100 RETURN
      END
      SUBROUTINE GRAPH(XX,YY,A,B,N,NPP)
C             SUBROUTINE GRAPH FOR BMD05R             JUNE 10, 1966
      DIMENSION XX(500),YY(500),B(10),PP(500),FMT(104)
      DIMENSION SYM(15),YYY(15)
      DATA O,P/1HO,1HP/
      DO 10 I=1,N
      SUM=0.0
      AA=1.0
      DO 5 J=1,NPP
      AA=AA*XX(I)
    5 SUM=SUM+B(J)*AA
      PP(I)=SUM+A
   10 CONTINUE
      WRITE (6,100)
      WRITE (6,101)
      DO 15 I=1,N
      DIFF=YY(I)-PP(I)
   15 WRITE (6,102)I,XX(I),YY(I),PP(I),DIFF
      SS=10.0**10
      FF=-10.0**9
      XS=SS
      XL=FF
      DO 300 I=1,N
      XS=AMIN1(XS,XX(I))
      XL=AMAX1(XL,XX(I))
      SS=AMIN1(SS,YY(I),PP(I))
  300 FF=AMAX1(FF,YY(I),PP(I))
      WRITE (6,4001)
      WRITE (6,4000)
C     MXY AND LYNN ARE DUMMY VARIABLES FOR SUBR SCALE.
      CALL SCALE(SS,FF,100.0,MXY,S1,F1,LYNN)
  302 SYM(1)=O
      SYM(2)=P
      X=XX(1)
      YYY(1)=YY(1)
      YYY(2)=PP(1)
      NC=2
      CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
      DO 50 I=2,N
      X=XX(I)
      YYY(1)=YY(I)
      YYY(2)=PP(I)
   50 CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
      NC=0
      CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
      WRITE (6,218)S1,F1
  100 FORMAT(1H15X,18HTABLE OF RESIDUALS)
  101 FORMAT(/6H0 NO. 4X,7HX VALUE4X,7HY VALUE4X,11HY PREDICTED4X,8HRESI
     1DUAL//)
  102 FORMAT(I5,F13.4,F11.4,2F13.4)
  218 FORMAT(26H GRAPH SCALE EXTENDS FROM F10.4,4H TO F10.4)
 4000 FORMAT(1H1,46X,38HPLOT OF OBSERVED AND PREDICTED VALUES./)
 4001 FORMAT(12H0GRAPH CODES//13H O=OBSERVED Y/14H P=PREDICTED Y/21H B=O
     1BSERVED=PREDICTED//)
      RETURN
      END
      SUBROUTINE MATINV(A,N,B,D,T)
C             SUBROUTINE MATINV FOR BMD05R            JUNE 10, 1966
C THE FOLLOWING STATEMENT(S) HAVE BEEN MANUFACTURED BY THE TRANSLATOR---
C
      DOUBLE PRECISION  A      , B      , U      , Y      , V
      DIMENSION A(11,11),B(10),U(10)
      DIMENSION IND(10),C(10)
      T=1.0
      DO 50 I=1,N
      IND(I)=0
   50 C(I)=A(I,I)
      D=1.0
      DO 8 L=1,N
      X=0.0
      DO 4 I=1,N
      IF(IND(I))4,2,4
    2 IF(X-A(I,I)/C(I))3,4,4
    3 K=I
      X=A(I,I)/C(I)
    4 CONTINUE
      T=AMIN1(T,X)
   11 IND(K)=1
      DO 5 I=K,N
      U(I)=A(I,K)
    5 A(I,K)=0.D0
      Y=U(K)
      D=D*SNGL(Y)
      DO 6 I=1,K
      U(I)=A(K,I)
    6 A(K,I)=0.D0
      V=B(K)
      B(K)=0.D0
      U(K)=1.D0
      DO 8 I=1,N
      DO 7 J=1,I
    7 A(I,J)=A(I,J)-U(I)*U(J)/Y
    8 B(I)=B(I)-U(I)*V/Y
      RETURN
      END
      SUBROUTINE PLOTR(X,ZMIN,ZMAX,Y,SYM,WMIN,WMAX,NC,NP)
C
C             SUBROUTINE PLOTR (IBM 360)              AUGUST 13, 1966
C
C     'PLOTR' IS A UTILITY SUBPROGRAM FOR THE BMD... PROGRAMS WHICH
C     PLOTS EITHER SINGLE-LINE OR WHOLE-PAGE GRAPHS AND SETS UP
C     APPROPRIATE SCALING.  THE CALLING PARAMETERS ARE AS FOLLOWS -
C
C     X - THE VALUE OF THE INDEPENDENT VARIABLE
C     ZMIN - THE MINIMUM VALUE OF X FOR THIS PLOT
C     ZMAX - THE MAXIMUM VALUE OF X FOR THIS PLOT
C     Y - THE ARRAY CONTAINING THE VALUES OF UP TO 15 DEPENDENT VAR.'S
C     SYM - THE ARRAY CONTAINING THE SYMBOLS TO BE PLOTTED
C     WMIN - THE MINIMUM VALUE OF ALL Y'S FOR THIS PLOT
C     WMAX - THE MAXIMUM VALUE OF ALL Y'S FOR THIS PLOT
C     NC - THE NUMBER OF DEPENDENT VARIABLES
C               NC=-1 CLOSES A SINGLE-LINE PLOT
C               NC= 0 PRINTS AND CLOSES A WHOLE-PAGE PLOT
C     NP - THE CONTROL VARIABLE
C               NP=-1 PRINTS A SINGLE LINE
C               NP=0 OR NP=1 SETS UP A WHOLE-PAGE PLOT
C
C     THE PLOTTING ROUTINE MUST BE CALLED ONCE FOR EACH VALUE OF THE
C     INDEPENDENT VARIABLE THAT IS TO BE PLOTTED NO MATTER WHETHER IN
C     THE SINGLE-LINE OR WHOLE-PAGE MODE
C
      DIMENSION Y(15),CLAB(12),GF(10),FMT(12),XY(51,101),SYM(15)
      INTEGER XY,BLANKS
      DATA TC,TP,BLANKS/1H.,1H+,1H /
      DATA GF/            4H 1X,,4H 2X,,4H 3X,,4H 4X,,4H 5X,,4H 6X,,
     14H 7X,,4H 8X,,4H 9X,,4H 10X/
      DATA FMT/'(17X',' ','5(F1','2.3,','8X)/','7X, ',' ','4(F1','2.3,',
     1'8X),','F12.','3)  '/
C
 100   FORMAT(1H 6X5(F12.3,8X),F12.3/17X,5(F12.3,8X))
 101  FORMAT(1H F12.3,1X,103A1,F12.3)
  102 FORMAT (1H 13X,103A1)
 1000 FORMAT(1H  14X,101A1)
 1001 FORMAT(15X,20(5H+....),1H+)
C
      DATA NCC/2/
C    'NCC' ON THE INITIAL ENTRY TO PLOTR IS NON-ZERO BECAUSE OF THE DATA
C    STATEMENT ABOVE.
C
C    'NCC' IS 0 WHILE A PLOT IS BEING MADE.  IT IS 1 OR 2 AT OTHER TIMES
C
      IF(NCC) 50,48,50
C
C    THE VARIABLE 'KL' CONTROLS THE FUNCTIONING OF THE OPENING AND
C    CLOSING  SECTIONS OF PLOTR.  KL=0 INDICATES OPENING OF THE GRAPH,
C    KL=1 INDICATES CLOSING.
C
   50 KL=0
      CALL SCALE(WMIN,WMAX,100.0,JY,YMIN,YMAX,YIJ)
      YR=YMAX-YMIN
  230 J=JY
      IF(J*(J-10))204,201,201
C
C     THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN FIXED FORMAT
C     UNDER CONTROL OF KL
C
  201 IF(KL)220,220,231
C
  231 WRITE (6,1001)
      IF(KL)250,250,220
C
  220 CLAB(1)= YMIN
      DO 222 I=2,11
  222 CLAB(I)=CLAB(I-1)+YIJ
      WRITE (6,100)(CLAB(I),I=1,11,2),(CLAB(J),J=2,10,2)
      IF(KL)231,231,14
C
C     THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN A VARIABLE
C     FORMAT UNDER CONTROL OF KL AND JY FROM 'SCALE'
C
  204 IF(J-5)205,221,207
  207 J=J-5
  205 JYT=5-J
  221 CONTINUE
  226 FMT(2)=GF(JY)
      IF (KL) 225,225,227
C
  225 FMT(7)=GF(JY)
      TT=JY
      TT=TT*YIJ/10.0
      CLAB(1)= YMIN+TT
      DO 223 I=2,10
  223 CLAB(I)=CLAB(I-1) +YIJ
      WRITE (6,FMT)(CLAB(I),I=2,10,2),(CLAB(I),I=1,9 ,2)
      IF(KL)227,227,14
C
  227 IF(JY-5)208,209,208
  208 WRITE(6,1000)(TC,I=1,J    ),(TP,(TC,I=1,4),K=1,19),TP,(TC,I=1,JYT)
      IF (KL) 250,250,225
C
  209 WRITE (6,1001)
      IF (KL) 250,250,225
C
  250 CONTINUE
      NCC=0
      IC=0
      IF(NP)80,11,11
C
C    THIS SECTION PREPARES FOR A FULL PAGE PLOT BY FILLING IN XY WITH
C     BLANKS AND SETTING UP SCALING FOR THE INDEPENDENT VARIABLE 'X'
C
   11 DO 1 I=1,51
       DO 1 J=1,101
   1  XY(I,J)=BLANKS
      CALL SCALE (ZMIN,ZMAX,50.0,JX,XMIN,XMAX,XIJ)
      XR=XMAX-XMIN
      GO TO 48
C
C
C     ENTRY TO PLOTS CAN BE USED ONLY AFTER THE CALLING PARAMETERS
C     HAVE BEEN TRANSFERRED BY A CALL ON PLOTR.  THE CALL ON PLOTS
C     IS IDENTICAL WITH ENTRY TO PLOTR BUT IT ALLOWS THE PROGRAMMER TO
C     CALL THE PLOTTING ROUTINE WITHOUT HAVING TO INCLUDE THE PARAMETERS
C
   48 IF(NC)52,13,49
   49 IF(NP)80,10,10
C    THE FOLLOWING SECTION SETS UP A FULL PAGE BUT DOES NO PRINTING.
C    THIS SECTION IS REACHED BY SPECIFYING NC POSITIVE AND NP POSITIVE.
C
   10 DO 9 N=1,NC
      SYMB=SYM(N)
      XDIFFR=XMAX-X
      IF(XDIFFR)105,106,106
  105 XDIFFR=0.0
  106 YDIFFR=YMAX-Y(N)
      IF(YDIFFR)107,108,108
  107 YDIFFR=0.0
  108 L=51.0-(50.0*XDIFFR)/XR+.5
      K=101.0-(100.0*YDIFFR)/YR+.5
      CALL FORM2(SYMB,XY(L,K))
 9     CONTINUE
      GO TO 15
C
C     THE FOLLOWING SECTION CONSTRUCTS AND PLOTS ONE LINE OF A MULTILINE
C     GRAPH.  LOCATION ALONG THE AXIS OF THE PAPER IS PRINTED AT EVERY
C     STEP.  THIS SECTION IS ACCESSIBLE BY SPECIFYING NC POSITIVE AND
C     NP NEGATIVE.
C
   80 DO 86 I=1,101
   86 XY(1,I)=BLANKS
       L=1
      DO 95 N=1,NC
      SYMB=SYM(N)
      YDIFFR=YMAX-Y(N)
      IF(YDIFFR)860,865,865
  860 YDIFFR=0.0
  865 K=101.0-(100.0*YDIFFR)/YR+.5
   95 CALL FORM2(SYMB,XY(L,K))
      IF(MOD(IC,5))97,96,97
   96 W=TP
      GO TO 98
   97 W=TC
   98 WRITE (6,101)X,W,(XY(1,N),N=1,101),W,X
      IC=IC+1
      GO TO 15
C
C    THIS SECTION PLOTS OUT THE PREVIOUSLY PREPARED WHOLE PAGE GRAPH.
C    IT PRINTS LOCATION ALONG THE PAPER'S AXIS EVERY FIFTH STEP.  THIS
C    SECTION IS ACCESSED BY SPECIFYING NC=0.
C
   13 M=6-JX
      LL=50+M
      T=JX
      IF(5-JX)131,131,135
  131 T=0.0
  135 RLAB=XMAX-(T*XIJ)/5.0
      W=TC
      K=52
      DO 31 L=M,LL
      K=K-1
      I=MOD(L,5)
      IF(I-1)2,3,2
    3 W=TP
      WRITE (6,101)RLAB,W,(XY(K,N),N=1,101),W,RLAB
      RLAB=RLAB-XIJ
      W=TC
      GO TO 31
    2 WRITE (6,102)W,(XY(K,N),N=1,101),W
   31 CONTINUE
C
   52 KL=1
      GO TO 230
C
   14 NCC=1
   15 RETURN
      END
      SUBROUTINE REPORT (SSAR,SSDR,N1,N,NPP)
C             SUBROUTINE REPORT FOR BMD05R            JUNE 10, 1966
      DIMENSION SS(10)
      IF(N1-100) 10, 30, 30
 10   NN=N-1
      N2=NN-N1
      D1=N1
      D2=N2
      SSARM=SSAR/D1
      SSDRM=SSDR/D2
      F=SSARM/SSDRM
      TOTAL=SSAR+SSDR
      IF(N1-1) 15, 15, 20
 15   WRITE (6,900)
      GO TO 21
 20   WRITE (6,901)N1
 21   WRITE (6,902)
      WRITE (6,903)
      WRITE (6,904)N1,SSAR,SSARM,F
      WRITE (6,905)N2,SSDR,SSDRM
      WRITE (6,906)NN, TOTAL
      NPP=NPP+1
      SS(NPP)=SSAR
      GO TO 70
 30   N1=1
      WRITE (6,907)NPP
      WRITE (6,921)
      WRITE (6,922)
      DO 60 I=1,NPP
      GO TO (31,32,33,34,35,36,37,38,39,40), I
 31   WRITE (6,911)N1,SS(1),SS(1)
      GO TO 60
 32   SSAR=SS(2)-SS(1)
      IF(SSAR) 41, 42, 42
 41   SSAR=0.0
 42   WRITE (6,912)N1,SSAR,SSAR
      GO TO 60
 33   SSAR=SS(3)-SS(2)
      IF(SSAR) 43, 44, 44
 43   SSAR=0.0
 44   WRITE (6,913)N1,SSAR,SSAR
      GO TO 60
 34   SSAR=SS(4)-SS(3)
      IF(SSAR) 45, 46, 46
 45   SSAR=0.0
 46   WRITE (6,914)N1,SSAR,SSAR
      GO TO 60
 35   SSAR=SS(5)-SS(4)
      IF(SSAR) 47, 48, 48
 47   SSAR=0.0
 48   WRITE (6,915)N1,SSAR,SSAR
      GO TO 60
 36   SSAR=SS(6)-SS(5)
      IF(SSAR) 49, 50, 50
 49   SSAR=0.0
 50   WRITE (6,916)N1,SSAR,SSAR
      GO TO 60
 37   SSAR=SS(7)-SS(6)
      IF(SSAR) 52, 53, 53
 52   SSAR=0.0
 53   WRITE (6,917)N1,SSAR,SSAR
      GO TO 60
 38   SSAR=SS(8)-SS(7)
      IF(SSAR) 54, 55, 55
 54   SSAR=0.0
 55   WRITE (6,918)N1,SSAR,SSAR
      GO TO 60
 39   SSAR=SS(9)-SS(8)
      IF(SSAR) 56, 57, 57
 56   SSAR=0.0
 57   WRITE (6,919)N1,SSAR,SSAR
      GO TO 60
 40   SSAR=SS(10)-SS(9)
      IF(SSAR) 58, 59, 59
 58   SSAR=0.0
 59   WRITE (6,920)N1,SSAR,SSAR
 60   CONTINUE
      WRITE (6,905)N2,SSDR,SSDRM
      WRITE (6,906)NN,TOTAL
 900  FORMAT(65H0               ANALYSIS OF VARIANCE FOR SIMPLE LINEAR R
     1EGRESSION)
 901  FORMAT(41H0                ANALYSIS OF VARIANCE FOR I4, 19H  DEGRE
     1E POLYNOMIAL)
 902  FORMAT(85H0   SOURCE OF VARIATION              DEGREE OF         S
     1UM OF         MEAN          F)
 903  FORMAT(87H                                      FREEDOM          S
     1QUARES       SQUARE       VALUE)
 904  FORMAT(37H0DUE TO REGRESSION.............      I6,F19.5,2F14.5)
 905  FORMAT(37H DEVIATION ABOUT REGRESSION....      I6,F19.5,F14.5)
 906  FORMAT(37H                      TOTAL....      I6,F19.5///)
 907  FORMAT(45H0              FINAL ANALYSIS OF VARIANCE FOR I4, 19H  D
     1EGREE POLYNOMIAL)
 911  FORMAT(37H0LINEAR TERM...................      I6,F19.5,F14.5)
 912  FORMAT(37H QUADRATIC TERM................      I6,F19.5,F14.5)
 913  FORMAT(37H CUBIC TERM....................      I6,F19.5,F14.5)
 914  FORMAT(37H QUARTIC TERM..................      I6,F19.5,F14.5)
 915  FORMAT(37H QUINTIC TERM..................      I6,F19.5,F14.5)
 916  FORMAT(37H SEXTIC TERM...................      I6,F19.5,F14.5)
 917  FORMAT(37H 7TH DEGREE TERM...............      I6,F19.5,F14.5)
 918  FORMAT(37H 8TH DEGREE TERM...............      I6,F19.5,F14.5)
 919  FORMAT(37H 9TH DEGREE TERM...............      I6,F19.5,F14.5)
 920  FORMAT(37H 10TH DEGREE TERM..............      I6,F19.5,F14.5)
 921  FORMAT(74H0   SOURCE OF VARIATION              DEGREE OF         S
     1UM OF         MEAN)
 922  FORMAT(75H                                      FREEDOM          S
     1QUARES       SQUARE)
 70   RETURN
      END
C     SUBROUTINE 'SCALE' CALCULATES THE SCALING FOR 'PLOTR'
C
C             SUBROUTINE SCALE FOR PLOTR              JUNE 21, 1966
C
      SUBROUTINE SCALE(YMIN,YMAX,YINT,JY,TYMIN,TYMAX,YIJ)
      DIMENSION C(10)
      DATA C           /1.0,1.5,2.0,3.0,4.0,5.0,7.5,10.0,15.0,20.0/
        DATA TEST / 0.76293945E-05/
   50 YR=YMAX-YMIN
      TT=YR/YINT
      J = ALOG10(TT)+TEST
      E=10.0**J
      TT=TT/E
      I=0
      IF(TT-1.0+TEST)205,201,201
  205 TT=TT*10.0
      E=E/10.0
 201  I=I+1
      IF(9-I)1,2,2
    1 E=E*10.0
      I=1
    2 IF(TT-C(I))233,202,201
  233 YIJ=C(I)*E
      GO TO 203
  202 Y=YMIN/C(I)
      J=Y
      T=J
      IF(0.0001-ABS(T-Y))204,233,233
  204 YIJ=C(I+1)*E
  203 X=((YMAX+YMIN)/YIJ-YINT )/2.0+.00001
      K=X
      IF(K)235,240,240
  235 Y=K
      IF(X-Y)236,240,236
  236 K=K-1
  240 TYMIN=K
      TYMIN=YIJ*TYMIN
      TYMAX=TYMIN+YINT*YIJ
      IF (YMAX-TYMAX-TEST)11,11,201
   11 YIJJ=C(I)*E
      XT=((YMAX+YMIN)/YIJJ-YINT)/2.0+.00001
      KT=XT
      IF (KT) 1235,1240,1240
 1235 YT=KT
      IF (XT.NE.YT) KT=KT-1
 1240 TYMINT=KT
      TYMINT=YIJJ*TYMINT
      TYMAXT=TYMINT+YINT*YIJJ
      IF (YMAX-TYMAXT.GT.TEST) GO TO 10
      TYMIN=TYMINT
      TYMAX=TYMAXT
      YIJ=YIJJ
      K=KT
   10 TT=YINT/10.0
      JY=TT+.000001
      YIJ=YINT*(YIJ/10.0)
      J=TYMIN/ YIJ
      IF (K)242,241,241
  242 J=J-1
  241 J=J*JY+JY-K
      JY=J
      RETURN
      END
      SUBROUTINE TPWD(NT1,NT2)
C             SUBROUTINE TPWD FOR BMD05R              JUNE 10, 1966
      IF(NT1)40,10,12
 10   NT1=5
 12   IF(NT1-NT2)14,19,14
   14 IF(NT2.EQ.5)GO TO 18
 17   REWIND NT2
   19 IF(NT1-5)18,24,18
 18   IF(NT1-6)22,40,22
 22   REWIND NT1
 24   NT2=NT1
 28   RETURN
 40   WRITE (6,49)
 49   FORMAT(25H ERROR ON TAPE ASSIGNMENT)
      STOP
      END
      SUBROUTINE TRNGEN(N,NT,M,X,LCODE,CONST)
C             SUBROUTINE TRNGEN FOR BMD05R            JUNE 10, 1966
      DIMENSION           LCODE(8),CONST(8)
      DOUBLE PRECISION X(400,6)
      ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
 4    DO 210 I=1,M
      FM=CONST(I)
      JUMP=LCODE(I)
      DO 85 J=1,N
      D=X(J,NT)
      GO TO(10,20,30,40,50,60,70,80,90,110),JUMP
   10 IF(D)99,11,12
   11 D2=0.0
      GO TO 8
   12 D2=SQRT(D)
      GO TO 8
   20 IF(D)99,21,22
   21 D2=1.0
      GO TO 8
   22 D2=SQRT(D)+SQRT(D+1.0)
      GO TO 8
   30 IF(D)99,99,31
   31 D2=ALOG10(D)
      GO TO 8
   40 D2=EXP(D)
      GO TO 8
   50 IF(-D)52,11,99
   52 IF(D-1.0)53,54,99
   54 D2=3.14159265/2.0
      GO TO 8
   53 D2=ASN(SQRT(D))
      GO TO 8
   60 SAMP=N
      A=D/(SAMP+1.0)
      B=A+1.0/(SAMP+1.0)
      IF(A)99,61,62
   61 IF(-B)64,11,99
   64 D2=ASN(SQRT(B))
      GO TO 8
   62 IF(B)99,65,66
   65 D2=ASN(SQRT(A))
      GO TO 8
   66 D2=ASN(SQRT(A))+ASN(SQRT(B))
      GO TO 8
   70 IF(D)71,99,71
   71 D2=1.0/D
      GO TO 8
   80 D2=D+FM
      GO TO 8
   90 D2=D*FM
      GO TO 8
 99   WRITE (6,105)J,JUMP
      GO TO 85
  110 IF(D)111,11,111
  111 D2=D**FM
    8 X(J,NT)=D2
   85 CONTINUE
  210 CONTINUE
 1000 RETURN
 105  FORMAT(18H0THE VALUE OF CASEI4,83H OF THE DEPENDENT VARIABLE VIOLA
     1TED THE RESTRICTION FOR TRANSGENERATION OF THE TYPE, I3,1H./52H TH
     2E PROGRAM CONTINUED LEAVING THIS VALUE UNCHANGED.)
      END





























