C-----------------------------------------------------
C
C  Program name:  YATES
C
C  Date:          October 1987
C
C  Revised:       August 1988
C
C-----------------------------------------------------
C
      REAL*8      C(128), EFF(128), PN
      INTEGER     R, N, CH1, CH2
      CHARACTER*7 TC(128), TE(128)
C
      WRITE(*,100)
  100 FORMAT(////,25X,'YATES ALGORITHM',////)
C
C  Enter data:
C
  101 WRITE(*,103)
  103 FORMAT(8(/),1X,'CHOOSE OPTION FOR ENTERING DATA:  ',
     +////,5X,' (1) ENTER DATA FROM KEYBOARD',
     +/,5X,' (2) FETCH STORED DATA FILE ',
     +////,1X,'ENTER NUMBER OF SELECTION:  ',/)
C
      READ(*,*) CH1
      IF (CH1.GE.1.AND.CH1.LE.2) GOTO 107
C
      WRITE(*,105)
  105 FORMAT(//,1X,'INVALID SELECTION.  PLEASE TRY ',
     +'AGAIN')
      GOTO 101
C
  107 GOTO (110,120), CH1
C
  110 CALL KEYBD(R,TC,TE,C)
      GOTO 125
C
  120 CALL FETCH(R,TC,TE,C)
C
C  Perform Yates algorithm and display results
C
  125 N=2**R
      DO 127 I=1,N
        EFF(I) = C(I)
  127 CONTINUE
C
      CALL YATES(R,EFF)
      WRITE(*,130)
  130 FORMAT(///,1X,'ESTIMATED EFFECTS CALCULATED',
     +' BY YATES ALGORITHM:',/)
      DO 140 I=1,N
        WRITE(*,135) TE(I),EFF(I)
  135   FORMAT(1X,A7,1X,'=',1X,F9.4)
  140 CONTINUE
      CALL PAUSE
C
C  Perform further analyses
C
  300 WRITE (*,301)
  301 FORMAT(8(/),1X,'CHOOSE OPTION FOR FURTHER ',
     +'ANALYSIS:',
     +//,10X,'1.  CONSTRUCT A NORMAL PLOT OF FITTED ',
     +'EFFECTS',
     +/,10X,'2.  CONSTRUCT A HALF-NORMAL PLOT OF ',
     +'ABSOLUTE FITTED EFFECTS',
     +/,10X,'3.  USE REVERSE YATES ALGORITHM',
     +/,10X,'4.  SAVE INPUT CELL MEANS IN DISK FILE',
     +/,10X,'5.  EXIT FROM THE PROGRAM',
     +//,1X,'ENTER NUMBER OF SELECTION:  ',/)
      READ(*,*) CH2
      IF(CH2.GE.1.AND.CH2.LE.5) GOTO 305
C
      WRITE(*,304)
  304 FORMAT(///,1X,'INVALID SELECTION, PLEASE TRY ',
     +'AGAIN')
      GOTO 300
C
  305 GOTO (310,310,330,340,350),CH2
C
  310 CALL PLOTS(CH2,N,EFF)
      GOTO 300
C
  330 CALL RYATES(N,R,EFF,TC,TE,C)
      GOTO 300
C
C
  340 CALL SAVE(R,C)
      GOTO 300
C
  350 WRITE(*,355)
  355 FORMAT(////)
      STOP
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE KEYBD(R,TC,TE,C)
C
      INTEGER     R,N
      REAL*8      C(128)
      CHARACTER*7 TC(128), TE(128)
C
  102 WRITE(*,104)
  104 FORMAT(1X,'ENTER NUMBER OF FACTORS ',
     +'( UP TO 7 ):  ',/)
      READ(*,*) R
C
C  Check for error in entering number of factors
C
      IF (R.GE.1.AND.R.LE.7) GOTO 110
      WRITE(*,105)
  105 FORMAT(///,1X,'NUMBER OF FACTORS MUST BE ',
     +'BETWEEN 1 AND 7',/,
     +1X,'PLEASE TRY AGAIN')
      GOTO 102
C
C  Create character strings representing treatment
C  combinations and treatment effects.
C
  110 CALL STRNGS(R,TC,TE)
C
C  Read cell means from keyboard.
C
      N=2**R
      WRITE(*,115)
  115 FORMAT(///,1X,'ENTER CELL MEAN FOR SPECIFIED',
     +' TREATMENT COMBINATION:',/)
      DO 125 I=1,N
        WRITE(*,120) TC(I)
  120   FORMAT(1X,A7,' = ',/)
        READ(*,*) C(I)
  125 CONTINUE
C
      RETURN
      END
C
C
C-----------------------------------------------------------
C
      SUBROUTINE FETCH(R,TC,TE,C)
C
      INTEGER     R, N
      REAL*8      C(128)
      CHARACTER*7 TC(128), TE(128)
      CHARACTER*8 FILENAME
C
      WRITE(*,61)
   61 FORMAT(//,' PLEASE ENTER DATA FILE NAME:  ',/)
      READ(*,63) FILENAME
   63 FORMAT(A8)
      IF (FILENAME.EQ.' ') GO TO 200
      OPEN (1,FILE=FILENAME,STATUS='OLD')
         READ(1,192) R
  192    FORMAT(I3)
         READ(1,*)   C
         CLOSE(1)
C
      CALL STRNGS(R,TC,TE)
      N=2**R
      CALL PRDATA(N,TC,C)
C
  200 RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE PRDATA(N,TC,C)
C
      INTEGER     N
      REAL*8      C(128)
      CHARACTER*7 TC(128)
C
      WRITE(*,100)
  100 FORMAT(///,1X,'CELL MEANS:',/)
      DO 140 I=1,N
        WRITE(*,120) TC(I),C(I)
  120   FORMAT(1X,A7,1X,'=',1X,F9.4)
  140 CONTINUE
      CALL PAUSE
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE SAVE(R,C)
C
      REAL*8      C(128)
      INTEGER     R
      CHARACTER*8 FILENAME
C
   70 WRITE(*,80)
   80 FORMAT(//,' PLEASE ENTER DATA FILE NAME (UP TO 8 CHARS):  ',/)

      READ(*,81) FILENAME
   81 FORMAT(A8)
      IF (FILENAME.EQ.' ') GO TO 70

      OPEN (1,FILE=FILENAME,STATUS='NEW')
        WRITE(1,92) R
   92   FORMAT(I3)
        WRITE(1,*)  C
        CLOSE(1)
C
      RETURN
      END
C
C
C----------------------------------------------------
C
      SUBROUTINE STRNGS(R,TC,TE)
C
      INTEGER R,NUM,BIN(7),DECIML
      CHARACTER*7 TC(128),TE(128),STRNGC,STRNGE
      CHARACTER*1 ARRAYC(7),ARRAYE(7),LTRSC(7),
     +LTRSE(7)
C
      EQUIVALENCE (STRNGC,ARRAYC),(STRNGE,ARRAYE)
C
      DATA LTRSC/'a','b','c','d','e','f','g'/
      DATA LTRSE/'A','B','C','D','E','F','G'/
C
C  Create the character strings that represent
C  treatment combinations and treatment effects by
C  using the binary representation of a number to
C  determine which letters appear in the string.
C
      NUM=2**R
      TC(1)='(1)'
      TE(1)='I'
      DO 200 I=1,NUM-1
C
C  First, convert a decimal number into binary
C
        DECIML=I
        DO 210 J=1,R
          N=2**(R-J)
          L=R-J+1
          IF (DECIML.GE.N) THEN
            BIN(L)=1
            DECIML=DECIML-N
          ELSE
            BIN(L)=0
          ENDIF
  210   CONTINUE
C
C  Now build the character strings associated with
C  the binary number.  First blank out the character
C  string,
C
        DO 220 J=1,7
          ARRAYC(J)=' '
          ARRAYE(J)=' '
  220   CONTINUE
C
C  then fill in the appropriate letters.
C
        K=1
        DO 230 J=1,R
          IF (BIN(J).EQ.0) GOTO 230
          ARRAYC(K)=LTRSC(J)
          ARRAYE(K)=LTRSE(J)
          K=K+1
  230   CONTINUE
C
C  Finally, put the strings into the arrays to be
C  used in the program.
C
        TC(I+1)=STRNGC
        TE(I+1)=STRNGE
  200 CONTINUE
C
      RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE YATES(R,EFF)
C
      REAL*8 EFF(128),STO(128)
      INTEGER I1,R,N
C
      N=2**R
      I1=2**(R-1)
C
      DO 830 I=1,R
        DO 810 J=1,I1
          STO(J)=EFF(2*J-1) + EFF(2*J)
          STO(J+I1)=EFF(2*J) - EFF(2*J-1)
  810   CONTINUE
        DO 820 J=1,N
          EFF(J)=STO(J)
  820   CONTINUE
  830 CONTINUE
C
      DO 840 I=1,N
        EFF(I)=EFF(I)/N
  840 CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE PLOTS(TYPE,NUM,EFF)
C
      INTEGER TYPE,NUM,CH2,NN
      REAL*8 EFF(128),P(128)
C
  404 WRITE(*,405)
  405 FORMAT(///,1X,'CHOOSE OPTION REGARDING GRAND ',
     +'MEAN:',
     +//,10X,'1.  INCLUDE GRAND MEAN IN THE PLOT',
     +/,10X,'2.  DO NOT INCLUDE GRAND MEAN IN THE ',
     +'PLOT',
     +//,1X,'ENTER NUMBER OF SELECTION:  ',/)
      READ(*,*) CH2
      IF (CH2.LT.1.OR.CH2.GT.2) GOTO 430
C
      GOTO (410,420),CH2
C
  410 CALL NPLOT(TYPE,NUM,EFF)
      GOTO 490
C
  420 DO 425 I=2,NUM
        P(I-1)=EFF(I)
  425 CONTINUE
      NN=NUM-1
      CALL NPLOT(TYPE,NN,P)
      GOTO 490
C
  430 WRITE(*,435)
  435 FORMAT(///,1X,'INVALID SELECTION, PLEASE TRY ',
     +'AGAIN')
      GOTO 404
C
  490 RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE NPLOT(TYPE,NUM,VALS)
C
      INTEGER TYPE,NUM
      REAL*8 VALS(128),AVAL(128),SORTV(128),
     +SORTA(128),ARG,PN,Y(128)
C
      GOTO (510,530),TYPE
C
C  Contruct a normal plot
C
  510 DO 520 I=1,NUM
        ARG=(DBLE(I)-.50D0)/DBLE(NUM)
        Y(I)=PN(ARG)
  520 CONTINUE
      CALL SORT(NUM,VALS,SORTV)
      WRITE(*,525)
  525 FORMAT(///,20X,'PLOT OF Y=NORMAL QUANTILE ',
     +'VS. X=FITTED EFFECT')
      CALL PLTTR(NUM,SORTV,Y)
      GOTO 560
C
C  Construct a half-normal plot
C
  530 DO 540 I=1,NUM
        ARG=(DBLE(I)-.50D0)/DBLE(NUM)
        Y(I)=PN((ARG/2.0D0) + .50D0)
        AVAL(I)=DABS(VALS(I))
  540 CONTINUE
      CALL SORT(NUM,AVAL,SORTA)
      WRITE(*,545)
  545 FORMAT(///,13X,'PLOT OF Y=HALF-NORMAL ',
     +'QUANTILE VS. X=ABSOLUTE FITTED EFFECT')
      CALL PLTTR(NUM,SORTA,Y)
C
  560 RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE RYATES(N,R,EFF,TC,TE,C)
C
      INTEGER N,R
      REAL*8 EFF(128),EFFR(128),Y(128),T(128),
     +RES(128),C(128),PN,ARG,SORTR(128)
      CHARACTER*7 TC(128),TE(128)
C
C  Use subroutine INPUT to specify values for effects,
C  either the estimated value from Yates algorithm or
C  zero.  Estimated values from Yates algorithm are in
C  EFF, revised values to be used in reverse Yates
C  algorithm are put in EFFR.
C
      CALL INPUT(N,EFF,TE,EFFR)
C
C  Calculate fitted cell means by reverse Yates
C  algorithm, calculate residuals, and display results
C
      CALL YATES(R,EFFR)
      DO 610 I=1,N
        T(I)=EFFR(N-I+1)*(2**R)
  610 CONTINUE
      WRITE(*,615)
  615 FORMAT(///,12X,'TC',9X,'CELL MEAN',5X,
     +'FITTED CELL MEAN',5X,'RESIDUAL',/)
      DO 625 I=1,N
        EFFR(I)=T(I)
        RES(I)=C(I)-EFFR(I)
        WRITE(*,620) TC(I),C(I),EFFR(I),RES(I)
  620   FORMAT(11X,A7,5X,F9.4,8X,F9.4,9X,F9.4)
  625 CONTINUE
C
      CALL PAUSE
C
C  Display residual plot and normal plot
C
      WRITE(*,630)
  630 FORMAT(///,20X,'PLOT OF Y=RESIDUAL VS. ',
     +'X=FITTED CELL MEAN')
      CALL PLTTR(N,EFFR,RES)
C
      CALL SORT(N,RES,SORTR)
      DO 635 I=1,N
        ARG=(DBLE(I)-.50D0)/DBLE(N)
        Y(I)=PN(ARG)
  635 CONTINUE
      WRITE(*,640)
  640 FORMAT(///,20X,'PLOT OF Y=NORMAL QUANTILE ',
     +'VS. X=RESIDUAL')
      CALL PLTTR(N,SORTR,Y)
C
  645 RETURN
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE INPUT(N,ORIG,TE,REV)
C
      INTEGER N
      CHARACTER*7 TE(128)
      CHARACTER*1 CH
      REAL*8 ORIG(128),REV(128)
C
      WRITE(*,700)
  700 FORMAT(///,1X,'TO BEGIN THE REVERSE YATES ',
     +'ALGORITHM, SPECIFY WHICH EFFECTS SHOULD BE',
     +/,1X,'USED AS CALCULATED BY THE YATES ',
     +'ALGORITHM, AND WHICH SHOULD BE SET TO ZERO.',/)
C
      DO 730 I=1,N
  705   WRITE(*,710) TE(I),ORIG(I)
  710   FORMAT(1X,A7,' = ',F9.4,10X,
     +'SET TO ZERO?  (ENTER Y OR N)  ',/)
        READ(*,715) CH
  715   FORMAT(A1)
        IF (CH.NE.'N'.AND.CH.NE.'n'.AND.
     +CH.NE.'Y'.AND.CH.NE.'y') GOTO 720
        IF (CH.EQ.'N'.OR.CH.EQ.'n') REV(N-I+1)=ORIG(I)
        IF (CH.EQ.'Y'.OR.CH.EQ.'y') REV(N-I+1)=0.0D0
        GOTO 730
C
C  Note that the order of the effect values is
C  reversed
C
  720   WRITE(*,725)
  725   FORMAT(/,1X,'INVALID SELECTION, PLEASE TRY ',
     +'AGAIN')
        GOTO 705
C
  730 CONTINUE
C
      WRITE(*,735)
  735 FORMAT(///,1X,'EFFECT VALUES TO BE USED',
     +' BY THE REVERSE YATES ALGORITHM:',/)
      DO 750 I=1,N
        WRITE(*,740) TE(I),REV(N-I+1)
  740   FORMAT(1X,A7,1X,'=',1X,F9.4)
  750 CONTINUE
C
      CALL PAUSE
C
      RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE SORT(N,OBS,SOBS)
C
C  Sort algorithm due to Loeser, Communications of
C  the ACM, Vol. 17, No. 3, page 143.
C
      REAL*8 OBS(N),SOBS(N),S
C
      DO 190 J=1,N
        SOBS(J)=OBS(J)
  190 CONTINUE
C
      I=1
  101 IF (I-N) 102,102,103
  102 I=I+I
      GO TO 101
  103 M=I-1
  104 M=M/2
      IF (M) 110,110,105
  105 K=N-M
      DO 109 J=1,K
        I=J+M
  106   I=I-M
        IF (I) 109,109,107
  107   L=I+M
        IF (SOBS(L)-SOBS(I)) 108,108,109
  108   S=SOBS(I)
        SOBS(I)=SOBS(L)
        SOBS(L)=S
        GO TO 106
  109 CONTINUE
      GO TO 104
  110 RETURN
      END
C
C-----------------------------------------------------
C
      DOUBLE PRECISION FUNCTION PN(P)
C
C  Algorithm AS 111, Applied Statistics (1977),
C  vol. 26, pp 118-121.
C
      REAL*8 A0,A1,A2,A3,B1,B2,B3,B4,C0,C1,C2,C3,D1,
     +D2,P,Q,R
      A0=2.50662823884D0
      A1=-18.61500062529D0
      A2=41.39119773534D0
      A3=-25.44106049637D0
      B1=-8.47351093090D0
      B2=23.08336743743D0
      B3=-21.06224101826D0
      B4=3.13082909833D0
      C0=-2.78718931138D0
      C1=-2.29796479138D0
      C2=4.85014127135D0
      C3=2.32121276858D0
      D1=3.54388924762D0
      D2=1.63706781897D0
      Q=P-.50D0
      IF (DABS(Q) .GT. .42D0) GO TO 100
      R=Q*Q
      PN=Q*(((A3*R + A2)*R + A1)*R + A0)/
     +((((B4*R + B3)*R + B2)*R + B1)*R + 1.0D0)
      RETURN
C
  100 R=P
      IF (Q .GT. 0.0D0) R=1.0D0-P
      R=DSQRT(-DLOG(R))
      PN=(((C3*R + C2)*R +C1)*R + C0)/
     +((D2*R + D1)*R + 1.0D0)
      IF (Q .LT. 0.0D0) PN=-PN
      RETURN
C
      END
C
C
C-----------------------------------------------------
C
      SUBROUTINE PLTTR(N,HORIZ,VERT)
C
      INTEGER I,J,K,COUNT,PRNTFLG(75)
      REAL*8 VERT(128),HORIZ(128),UP,LO,MX,MN,X(128),
     +DX,DY,XMIN,XMAX
      CHARACTER*1 PRNTLN(68),PRNTCHR(6),GO
      DATA PRNTCHR/' ','*','*','*','*','*'/
C
C  Find largest and smallest horizontal coordinates
C  and determine the horizontal scaling factor, DX.
C
      MN=9999999999.9D0
      MX=-9999999999.9D0
      DO 900 I=1,N
       IF (HORIZ(I) .GT. MX) MX=HORIZ(I)
       IF (HORIZ(I) .LT. MN) MN=HORIZ(I)
  900 CONTINUE
      IF (MX.NE.MN) GOTO 902
      MX=MX+1.0D0
      MN=MN-1.0D0
  902 XMIN=MN
      XMAX=MX
      DX=(MX-MN)/64.0D0
C
C  Do the same for the vertical scale.
C
      MN=9999999999.9D0
      MX=-9999999999.9D0
      DO 905 I=1,N
       IF (VERT(I) .GT. MX) MX=VERT(I)
       IF (VERT(I) .LT. MN) MN=VERT(I)
  905 CONTINUE
      IF (MX.NE.MN) GOTO 908
      MX=MX+1.0D0
      MN=MN-1.0D0
  908 DY=(MX-MN)/15.0D0
C
C  Draw the plot a line at a time.
C
      WRITE (*,910)
  910 FORMAT(/,'   Y',5X,'|',/,9X,'|')
      LO=MX
      COUNT=3
C
      DO 940 I=1,17
        DO 915 J=1,65
         PRNTFLG(J)=1
  915   CONTINUE
        UP=LO
        LO=UP-DY
        DO 917 K=1,N
          IF ((VERT(K) .GT. LO) .AND.
     +        (VERT(K) .LE. UP)) THEN
            J=INT(SNGL((HORIZ(K)-XMIN)/DX))+1
            PRNTFLG(J)=PRNTFLG(J)+1
            IF (PRNTFLG(J) .GT. 6) PRNTFLG(J)=6
          END IF
  917   CONTINUE
        IF (COUNT.EQ.3) THEN
          WRITE (*,920) UP,
     +                  (PRNTCHR(PRNTFLG(J)),J=1,65)
  920     FORMAT (1X,F8.2,'+ ',65A1)
          COUNT=0
        ELSE
          WRITE (*,925) (PRNTCHR(PRNTFLG(J)),J=1,65)
  925     FORMAT (9X,'| ',65A1)
          COUNT=COUNT+1
        ENDIF
  940 CONTINUE
C
C  Put the labels under the horizontal axis.
C
      WRITE (*,945)
  945 FORMAT (10X,'-+-------+-------+-------+----',
     +'---+-------+-------+-------+-------+----')
      DX=(XMAX-XMIN)/8.0D0
      DO 950 I=1,9
        X(I)=XMIN+DBLE(I-1)*DX
  950 CONTINUE
C
      WRITE (*,960) (X(I),I=1,9)
  960 FORMAT(7X,9(F6.1,2X),'X')
C
      CALL PAUSE
C
      RETURN
      END
C
C-----------------------------------------------------
C
      SUBROUTINE PAUSE
C
      CHARACTER*1 GO
C
C  Pause to let the user view a display.
C
      WRITE(*,70)
   70 FORMAT(/,1X,'STRIKE CARRIAGE RETURN ',
     +'WHEN YOU WISH TO CONTINUE  ',/)
      READ(*,75) GO
   75 FORMAT(A1)
C
      RETURN
      END
C
C----------------------------------------------------