C                                            *** STAT PACK ***
C     ROUTINE TO CALCULATE POINT BISERIAL CORRELATIONS
C     CALLING SEQUENCE: CALL PTBIS(NV,NC,MV,MC,DATA,STD,IO,NAMES)
C     WHERE NV - NUMBER OF VARIABLES USED
C           NC - NUMBER OF OBSERVATIONS USED
C           MV - NUMBER OF VARIABLES DIMENSIONED FOR
C           MC - NUMBER OF OBSERVATIONS DIMENSIONED FOR
C           DATA - MATRIX CONTAINING DATA
C           STD - VECTOR CONTAINING VARIABLE STANDARD DEVIATIONS
C           IO - VECTOR DIMENSIONED AT LEAST NO
C           NAMES- VECTOR CONTAINING VARIABLE NAMES
C
C     PROGRAM MAKES USE OF SINGLETON SORT TECHNIQUE FOR DICHOTOMY
C      ORIGINAL ROUTINE WRITTEN BY SAM ANEMA.  SUGESTED FOR 
C     INCLUSION IN STP BY SAM ANEMA (FACULTY COORDINATOR - WMU)
C
      SUBROUTINE PTBIS(NV,NC,MV,MC,DATA,STD,IO,NAMES)
      DIMENSION IV(51),IU(16),IL(16),NAMES(1),STD(1),INPUT(15)
      DIMENSION NUMBER(3),IO(1),DATA(MC,MV)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON /EXTRA/HEDR(70),NSZ
      EQUIVALENCE (IDV,IV(51))
      ALL=0
      IB=0
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT(' WHICH VARIABLES? ',$)
      CALL ALPHA(IV,50,NN,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 1
      IF(IHELP.NE.1) GO TO 200
      WRITE(IDLG,3)
3     FORMAT('0ENTER THE VARIABLES TO BE USED IN THE POINT-BISERIAL'/
     1' CORRELATION.  EITHER VARIABLE NAMES(IF THEY HAVE BEEN'/
     2' SUPPLIED) OR VARIABLE NUMBERS MAY BE USED.  UP TO 50'/
     3' VARIABLES MAY BE ENTERED.  THE DICHOTOMOUS VARIABLE'/
     4' MAY BE ENTERED WITH THE OTHERS, HOWEVER IT IS NOT NECESSARY')
      WRITE(IDLG,6)
      GO TO 1
200   DO 201 I=1,NN
      IF(IV(I).GE.1) GO TO 201
      ALL=1
201   CONTINUE
C
C     DICHOTOMOUS VARIABLE
C
4     IF(ICC.NE.2) WRITE(IDLG,5)
5     FORMAT(' WHICH IS THE DICHOTOMOUS VARIABLE? ',$)
      CALL ALPHA(IV(51),1,N,IRET,IHELP,IERR,NAMES,NV)
      IF(IRET.EQ.1) RETURN
      IF(IERR.EQ.1) GO TO 4
      IF(IHELP.NE.1) GO TO 7
      WRITE(IDLG,6)
6     FORMAT('0THE DICHOTOMOUS VARIABLE WILL BE CORRELATED WITH ALL'/
     1' THE OTHER VARIABLES SPECIFIED USING A POINT-BISERIAL'/
     2' CORRELATION.  FIRST THE VARIABLE WILL BE CHECKED TO SEE'/
     3' IF ALL THE OBSERVATIONS ARE EQUAL TO 1 OF 2 VALUES, IN'/
     4' WHICH CASE THE ANALYSIS WILL BE PERFORMED.  IF A THIRD OR'/
     5' MORE VALUE EXISTS FOR THE DICHOTOMOUS VARIABLE, IT WILL BE'/
     1' NECESSARY TO ENTER A BREAK POINT WHEN INSTRUCTED.'/)
      GO TO 4
C
C     CHECK DICHOTOMOUS VARIABLE
C
7     IF(IV(51).GT.0) GO TO 9
      WRITE(IDLG,8)
8     FORMAT(' IT IS ILLEGAL TO SPECIFY ALL VARIABLES AS DICHOTOMOUS')
      GO TO 4
C
C     SORT ON DICHOTOMOUS VARIABLE SINGLETON NETHOD BUT BY INDEX
C
9     DO 15 I=1,NC
15    IO(I)=I
      M=1
      II=1
      J=NC
21    IF(II.GE.J) GO TO 28
22    K=II
      IJ=(J+II)/2
      T=DATA(IO(IJ),IDV)
      IF(DATA(IO(II),IDV).LE.T) GO TO 23
      ISAV=IO(IJ)
      IO(IJ)=IO(II)
      IO(II)=ISAV
      T=DATA(IO(IJ),IDV)
23    LL=J
      IF(DATA(IO(J),IDV).GE.T) GO TO 25
      ISAV=IO(IJ)
      IO(IJ)=IO(J)
      IO(J)=ISAV
      T=DATA(IO(IJ),IDV)
      IF(DATA(IO(II),IDV).LE.T) GO TO 25
      ISAV=IO(IJ)
      IO(IJ)=IO(II)
      IO(II)=ISAV
      T=DATA(IO(IJ),IDV)
      GO TO 25
24    ISAV=IO(LL)
      IO(LL)=IO(K)
      IO(K)=ISAV
25    LL=LL-1
      IF(DATA(IO(LL),IDV).GT.T) GO TO 25
      TT=DATA(IO(LL),IDV)
26    K=K+1
      IF(DATA(IO(K),IDV).LT.T) GO TO 26
      IF(K.LE.LL) GO TO 24
      IF((LL-II).LE.(J-K)) GO TO 27
      IL(M)=II
      IU(M)=LL
      II=K
      M=M+1
      GO TO 29
27    IL(M)=K
      IU(M)=J
      J=LL
      M=M+1
      GO TO 29
28    M=M-1
      IF(M.EQ.0) GO TO 40
      II=IL(M)
      J=IU(M)
29    IF((J-II).GE.11) GO TO 22
      IF(II.EQ.1) GO TO 21
      II=II-1
30    II=II+1
      IF(II.EQ.J) GO TO 28
      ISAV=IO(II+1)
      T=DATA(IO(II+1),IDV)
      IF(DATA(IO(II),IDV).LE.T) GO TO 30
      K=II
31    IO(K+1)=IO(K)
      K=K-1
      IF(T.LT.DATA(IO(K),IDV)) GO TO 31
      IO(K+1)=ISAV
      GO TO 30
C
C     SORT DONE CHECK FOR DICHOTOMY
C
40    SAV=DATA(IO(1),IDV)
      K=1
      DO 41 I=2,NC
      IF(SAV.EQ.DATA(IO(I),IDV)) GO TO 41
      K=K+1
      IF(K.GT.2) GO TO 45
      NL=I
      SAV=DATA(IO(I),IDV)
41    CONTINUE
      IF(K.GT.1) GO TO 60
      WRITE(IDLG,42) NAMES(IDV)
42    FORMAT(' ALL OBSERVATIONS OF THE VARIABLE: ',A5,' ARE EQUAL')
      RETURN
C
C     NOT A 2 VALUE VARIABLE ASK FOR BREAK POINT
C
45    IF(ICC.NE.2) WRITE(IDLG,46) NAMES(IDV)
46    FORMAT(' WHAT IS THE BREAKPOINT FOR VARIABLE: ',A5,'? ',$)
      IB=1
      READ(ICC,47,END=45) INPUT
47    FORMAT(15A1)
      IF(INPUT(1).EQ.'!') RETURN
      IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.
     1(INPUT(3).NE.'L').OR.(INPUT(4).NE.'P')) GO TO 49
      WRITE(IDLG,48)
48    FORMAT(' TO ENTER THE BREAKPOINT FOR THE DICHOTOMOUS VARIABLE'/
     1' CHOSE A VALUE SOMEWHERE IN THE RANGE OF THE VARIABLE AND'/
     2' ENTER IT.  THE VARIABLE WILL THEN BE BROKEN INTO 2 GROUPS'/
     3' BASED ON THIS VALUE; 1 GROUP WILL CONTAIN ALL OBSERVATIONS'/
     4' WHERE THE VALUE IS LESS THAN OR EQUAL TO THE BREAKPOINT VALUE'/
     5' AND THE OTHER GROUP WILL CONTAIN ALL OBSERVATIONS WHERE'/
     6' THE VALUE IS GREATER THAN THE BREAKPOINT.  THE MEDIAN MAY'/
     7' ALSO BE SPECIFIED AS THE BREAKPOINT BY TYPING "MEDIAN"'/)
      GO TO 45
49    IF((INPUT(1).NE.'M').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'D')
     1.OR.(INPUT(4).NE.'I').OR.(INPUT(5).NE.'A').OR.(INPUT(6).NE.'N'))
     3 GO TO 51
      BREAK=DATA(IO((NC+1)/2),IDV)
      IF((NC.AND.1).EQ.0) BREAK=(DATA(IO(NC/2),IDV)+DATA(IO(NC/2+1),
     1IDV))/2.
C
C     CHECK TO SEE THAT MEDIAN IS LESS THAN LAST VALUE
C
      IF(BREAK.NE.DATA(IO(NC),IDV)) GO TO 54
      DO 56 I=NC/2,1,-1
      IF(BREAK.EQ.DATA(IO(I),IDV)) GO TO 56
      BREAK=DATA(IO(I),IDV)
      GO TO 54
56    CONTINUE
      PAUSE 'PROBLEM PT BISERIAL'
51    ENCODE(15,47,NUMBER) INPUT
      DECODE(15,52,NUMBER) BREAK
52    FORMAT(G)
      IF((BREAK.GE.DATA(IO(1),IDV)).AND.(BREAK.LT.DATA(IO(NC),IDV)) )
     1GO TO 54
      WRITE(IDLG,53)
53    FORMAT(' BREAKPOINT WAS NOT IN RANGE')
      GO TO 45
54    DO 55 I=1,NC
      IF(DATA(IO(I),IDV).LE.BREAK) GO TO 55
      NL=I
      GO TO 60
55    CONTINUE
      PAUSE 'BREAK POINT NOT FOUND'
C
C     DONE WITH BREAK VARIABLE NOW DO ANALYSIS
C
60    IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(J),J=1,NSZ)
5566  FORMAT('1',70A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      LINES=10
      WRITE(IOUT,62) NAMES(IDV)
62    FORMAT('0',10X,'**** POINT BISERIAL CORRELATION ****'/
     110X,' VARIABLE:',A5,' IS THE DICHOTOMOUS VARIABLE')
      IF(IB.EQ.1) WRITE(IOUT,63) BREAK
63    FORMAT(3X,' THE BREAKPOINT USED TO SPLIT THE VARIABLE IS',
     1G15.7)
      IF(IB.EQ.1) LINES=LINES+1
      NH=NC-NL+1
      NB=NL-1
      WRITE(IOUT,64) NB,NH
64    FORMAT(1X,' THE LOWER GROUP SIZE IS ',I4,
     1' AND THE UPPER GROUP SIZE IS ',I4)
      WRITE(IOUT,65) NAMES(IDV)
65    FORMAT('0',56X,'POINT-BISERIAL'/13X,'MEAN OF',7X,'MEAN OF',
     16X,'STANDARD',7X,'CORRELATION WITH'/' VARIABLE   ',
     2'LOW GROUP     HIGH GROUP    DEVIATION',8X,'VARIABLE:',A5)
      IF(ALL.EQ.1) NN=NV
      DO 61 I=1,NN
      IF(ALL.EQ.1) IVAR=I
      IF(ALL.NE.1) IVAR=IV(I)
      IF(IVAR.EQ.IDV) GO TO 61
      XBL=0
      DO 66 J=1,NL-1
66    XBL=XBL+DATA(IO(J),IVAR)
      XBH=0
      DO 67 J=NL,NC
67    XBH=XBH+DATA(IO(J),IVAR)
      R=0
      XMP=XBH/NH
      XMQ=XBL/NB
      IF(IOUT.NE.21) GO TO 71
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 71
      CALL PRNTHD
      WRITE(IOUT,62) NAMES(IDV)
      LINES=11
      IF(IB.EQ.1) WRITE(IOUT,63) BREAK
      IF(IB.EQ.1) LINES=LINES+1
      WRITE(IOUT,64) NB,NH
      WRITE(IOUT,65) NAMES(IDV)
71    IF(STD(IVAR).EQ.0) GO TO 69
      R=(XMP-XMQ)*SQRT(FLOAT(NH*NB)/FLOAT(NC**2))/STD(IVAR)
      WRITE(IOUT,68) NAMES(IVAR),XMQ,XMP,STD(IVAR),R
68    FORMAT(3X,A5,2X,G14.5,1X,G14.5,1X,G14.5,2X,F10.6)
      GO TO 61
69    WRITE(IOUT,70) NAMES(IVAR),XMQ,XMP,STD(IVAR)
70    FORMAT(3X,A5,2X,G14.5,1X,G14.5,1X,G14.5,2X,'CANNOT CALCULATE')
61    CONTINUE
      RETURN
      END
C                                   *** STAT PACK ***
C     SUBROUTINE TO CREATE BAR GRAPHS (HORIZ HISTOGRAMS)
C     CALLING SEQUENCE: CALL BARGR(NV,NC,MV,MC,DATA,AV,NAMES)
C     WHERE NV - NUMBER OF VARIABLES USED
C           NC - NUMBER OF OBSERVATIONS USED
C           MV - MAXIMUM NUMBER OF VARIABLES
C           MC - MAXIMUM NUMBER OF OBSERVATIONS
C           DATA - MATRIX CONTAINING DATA
C           AV - EXTRA VECTOR (AT LEAST NC LONG)
C           NAMES - VECTOR CONTAINING VARIABLE NAMES
C
C     THIS ROUTINE WAS WRITTEN IN RESPONSE TO A NEED EXPRESSED BY
C     MIKE KEENAN AND TOM BRAYTON OF THE MANAGEMENT DEPARTMENT. 
C     THEY BELIEVED THAT THE HIST COMMAND WAS NOT VERSITALE ENOUGH
C     AND NEEDED REVISION.  TOO MANY PEOPLE ARE USE TO HIST SO THE
C     DECISION WAS MADE TO CREATE THE BARGR COMMAND TO MEET THEIR
C     NEEDS
C
      SUBROUTINE BARGR(NV,NC,MV,MC,DATA,AV,NAMES)
      DIMENSION DATA(MC,MV),HV(50),BV(50),PCNT(50),IFRQ(50)
      DIMENSION AV(1),NAMES(1),IN(40),IV(40)
      DIMENSION NET(2),IU(16),IL(16),ICHAR(80),CHAR(4)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON /EXTRA/ HEDR(70),NSZ
      NGRP=50
      AUTO=0
      NOUT=20
      IF(IOUT.EQ.21) NOUT=80
      DASH='-'
      PLUS='+'
      ARROW='^'
      BLANK=' '
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT(' WHICH VARIABLES? ',$)
      CALL ALPHA(IV,40,NN,IRET,IHELP,IERR,NAMES,NV)
      IF(IERR.EQ.1) GOTO 1
      IF(IRET.EQ.1) RETURN
      IF(IHELP.NE.1) GO TO 5
      WRITE(IDLG,3)
3     FORMAT(' ENTER THE VARIABLE NUMBERS OR NAMES OF VARIABLES (IF '/
     1' THEY HAVE BEEN DEFINED) TO BE GRAPHED.  IF ALL VARIABLES ARE'/
     2' TO BE EXAMINED TYPE "*".  YOU WILL  BE ASKED TO SPECIFY'/
     3' RANGES')
      WRITE(IDLG,4)
4     FORMAT('0RANGES ARE USED TO DETERMINE WHICH OBSERVATIONS FIT'/
     1' INTO EACH BAR OF THE GRAPH.  RANGES ARE ENTERED ONE PER LINE;'/
     2' EACH RANGE CONSISTING OF A MINIMUM AND A MAXIMUM VALUE'/
     3' SEPARATED BY A COMMA.  WHEN THE LAST RANGE HAS BEEN ENTERED'/
     4' TYPE A <CR>, OR A ^Z(CONTROL Z).  IF YOU DESIRE THE RANGES'/
     5' TO BE AUTOMATICALLY CALCULATED ANSWER "AUTO" TO THE INQUIRY'/
     6' FOR RANGES.  IF THERE ARE FEWER INDIVIDUAL VALUES THAN'/
     7' RANGES, EACH INDIVIDUAL VALUE IS PLACED IN ITS OWN RANGE,'/
     8' OTHERWISE RANGES OF EQUAL MAGNITUDE ARE ASSUMED.'/
     9' IF FEWER RANGES ARE DESIRED THEY MAY BE SPECIFIED'/
     1' BY FOLLOWING THE "AUTO" WITH A "/" AND THE NUMBER OF RANGES'/)
      GO TO 1
5     IF(ICC.NE.2) WRITE(IDLG,6)
6     FORMAT(' ENTER RANGES 1 PER LINE'/)
      J=1
7     IF(ICC.NE.2) WRITE(IDLG,8)
8     FORMAT('+? ',$)
      READ(ICC,9,END=24) IN
9     FORMAT(40A1)
      IF(IN(1).EQ.'!') RETURN
      IF((IN(1).NE.'H').OR.(IN(2).NE.'E').OR.(IN(3).NE.'L').OR.
     1(IN(4).NE.'P')) GO TO 10
      WRITE(IDLG,4)
      GO TO 7
10    IF((IN(1).NE.'A').OR.(IN(2).NE.'U').OR.(IN(3).NE.'T').OR.
     1(IN(4).NE.'O')) GO TO 20
      AUTO=1
      IF(IN(5).NE.'/') GO TO 25
      DO 11 I=1,2
11    NET(I)=' '
      J=1
      I=5
12    I=I+1
      IF((IN(I).LT.'0').OR.(IN(I).GT.'9')) GO TO 15
      IF(J.GT.2) GO TO 13
      NET(J)=IN(I)
      J=J+1
      GO TO 12
13    WRITE(IDLG,14)
14    FORMAT(' MAXIMUM OF 50 GROUPS FOR BARGRAPHS')
      GO TO 5
15    IF(NET(1).EQ.' ') GO TO 7
      IF(NET(2).NE.' ') GO TO 16
      NET(2)=NET(1)
      NET(1)=' '
16    ENCODE(2,9,I) NET
      DECODE(2,17,I) NGRP
17    FORMAT(I2)
      IF(NGRP.GT.50) GO TO 13
      IF(NGRP.LT.1) GO TO 7
      GO TO 25
20    IF((IN(1).EQ.' ').AND.(IN(2).EQ.' ')) GO TO 24
      IF(J.GT.50) GO TO 13
      REREAD 21,BV(J),HV(J)
21    FORMAT(2F)
      IF(BV(J).LT.HV(J)) GO TO 31
      C=BV(J)
      BV(J)=HV(J)
      HV(J)=C
31    IF(J.LT.2) GO TO 22
      DO 32 K=1,J-1
      IF(HV(J).LT.BV(K)) GO TO 32
      IF(BV(J).GT.HV(K)) GO TO 32
      WRITE(IDLG,33) K
33    FORMAT(' THIS RANGE OVERLAPS RANGE #',I2,' - PLEASE REENTER'/)
      GO TO 7
32    CONTINUE
22    J=J+1
      GO TO 7
C
C     OK USER HAS SUBMITTED ALL HIS INFO NOW RUN
C
24    NGRP=J-1
      IF(NGRP.LT.1) GO TO 5
25    ALL=0
      DO 26 I=1,NN
      IF(IV(I).LT.0) ALL=1
26    CONTINUE
      IF(ALL.EQ.1) NN=NV
      DO 27 I=1,NN
      IANLY=I
      IF(ALL.NE.1) IANLY=IV(I)
      IF(AUTO.EQ.1) GO TO 39
C
C     IF NOT AUTO - RANGES WERE ENTERED BY HAND
C
      DO 28 K=1,NGRP
28    IFRQ(K)=0
      DO 29 J=1,NC
      DO 30 K=1,NGRP
      IF(DATA(J,IANLY).LT.BV(K)) GO TO 30
      IF(DATA(J,IANLY).GT.HV(K)) GO TO 30
      IFRQ(K)=IFRQ(K)+1
      GO TO 29
30    CONTINUE
29    CONTINUE
      NWORK=NGRP
      GO TO 100
C
C     AUTO WAS USED  SORT AND FORM FREQ ON INDIVIDUAL VALUES 
C     OTHERWISE ON RANGES  PARTITIONING SORT - ACM
C
39    DO 40 J=1,NC
40    AV(J)=DATA(J,IANLY)
      M=1
      II=1
      J=NC
41    IF(II.GE.J) GO TO 48
42    K=II
      IJ=(J+II)/2
      T=AV(IJ)
      IF(AV(II).LE.T) GO TO 43
      AV(IJ)=AV(II)
      AV(II)=T
      T=AV(IJ)
43    LL=J
      IF(AV(J).GE.T) GO TO 45
      AV(IJ)=AV(J)
      AV(J)=T
      T=AV(IJ)
      IF(AV(II).LE.T) GO TO 45
      AV(IJ)=AV(II)
      AV(II)=T
      T=AV(IJ)
      GO TO 45
44    AV(LL)=AV(K)
      AV(K)=TT
45    LL=LL-1
      IF(AV(LL).GT.T) GO TO 45
      TT=AV(LL)
46    K=K+1
      IF(AV(K).LT.T) GO TO 46
      IF(K.LE.LL) GO TO 44
      IF((LL-II).LE.(J-K)) GO TO 47
      IL(M)=II
      IU(M)=LL
      II=K
      M=M+1
      GO TO 49
47    IL(M)=K
      IU(M)=J
      J=LL
      M=M+1
      GO TO 49
48    M=M-1
      IF(M.EQ.0) GO TO 59
      II=IL(M)
      J=IU(M)
49    IF((J-II).GE.11) GO TO 42
      IF(II.EQ.1) GO TO 41
      II=II-1
50    II=II+1
      IF(II.EQ.J) GO TO 48
      T=AV(II+1)
      IF(AV(II).LE.T) GO TO 50
      K=II
51    AV(K+1)=AV(K)
      K=K-1
      IF(T.LT.AV(K)) GO TO 51
      AV(K+1)=T
      GO TO 50
C
C     NOW TRY FOR FREQUENCIES ON INDIVIDUAL VALUES
C
59    NWORK=1
      K=1
      IFRQ(K)=1
      HV(K)=AV(1)
      BV(K)=AV(1)
      IF(NC.LT.2) GO TO 100
      DO 60 J=2,NC
      IF(AV(J).EQ.HV(K)) GO TO 60
      K=K+1
      IF(K.GT.NGRP) GO TO 69
      HV(K)=AV(J)
      BV(K)=AV(J)
      IFRQ(K)=0
60    IFRQ(K)=IFRQ(K)+1
      NWORK=K
      GO TO 100
C
C     DIDN'T WORK WITH INDIVIDUAL VALUE BREAK INTO GROUPS OF EQUAL SIZE
C
69    AINC=(AV(NC)-AV(1))/NGRP
      BV(1)=AV(1)
      HV(1)=AV(1)+AINC-AINC/10000.
      IF(NGRP.LT.2) GO TO 71
      DO 70 J=2,NGRP
      BV(J)=BV(J-1)+AINC
70    HV(J)=BV(J)+AINC-AINC/10000.
      HV(NGRP)=AV(NC)
71    DO 72 K=1,NGRP
72    IFRQ(K)=0
      DO 73 J=1,NC
      DO 74 K=1,NGRP
      IF(AV(J).LT.BV(K)) GO TO 74
      IF(AV(J).GT.HV(K)) GO TO 74
      IFRQ(K)=IFRQ(K)+1
      GO TO 73
74    CONTINUE
73    CONTINUE
      NWORK=NGRP
C
C     FREQUENCIES ARE CALCULATED NOW FORM PERCENT AND CUM. PERCENT
C
100   TOTAL=0
      IATOT=100
      NATOT=100
      DO 101 J=1,NWORK
      PCNT(J)=0
101   TOTAL=TOTAL+IFRQ(J)
      IF(TOTAL.EQ.0) GO TO 120
      PCTMAX=0
      DO 102 J=1,NWORK
      PCNT(J)=(IFRQ(J)/TOTAL)*100.
      IF(PCNT(J).GT.PCTMAX) PCTMAX=PCNT(J)
102   CONTINUE
109   IF(PCTMAX.GT.NATOT) GO TO 120
      IATOT=NATOT
      NATOT=NATOT-20
      IF(NATOT.LE.0) NATOT=NATOT+18.
      IF(NATOT.EQ.18) NATOT=10.
      GO TO 109
120   IF(IOUT.NE.21) WRITE(IOUT,5566) (HEDR(J),J=1,NSZ)
5566  FORMAT('1',72A1)
      IF(IOUT.EQ.21) CALL PRNTHD
      WRITE(IOUT,103) NAMES(IANLY)
103   FORMAT('0',10X,'***** BAR GRAPH FOR VARIABLE: ',A5,' *****')
      NDASH=4
      IF(IOUT.EQ.21) NDASH=19
      WRITE(IOUT,104)((DASH,K=1,NDASH),PLUS,L=1,4)
104   FORMAT('0',2X,'RANGE OF VALUES',7X,'FREQ',2X,'PCENT',
     11X,'+',80A1)
      ATOT=IATOT
      AINC=ATOT/NOUT
      DO 105 J=1,NWORK
      L=PCNT(J)/AINC+.5
      IF(L.LT.1) GO TO 108
      DO 107 K=1,L
107   ICHAR(K)='X'
108   IF(L.LT.1) WRITE(IOUT,111) BV(J),HV(J),IFRQ(J),PCNT(J)
      IF(L.GE.1) WRITE(IOUT,111) BV(J),HV(J),IFRQ(J),PCNT(J)
     1,(ICHAR(K),K=1,L)
111   FORMAT(1X,G10.4,'- ',G10.4,2X,I4,2X,F5.1,1X,'I',80A1)
105   CONTINUE
      ATOT=IATOT
      CHAR(1)=ATOT/4.
      CHAR(2)=2*(ATOT/4)
      CHAR(3)=3*(ATOT/4)
      CHAR(4)=IATOT
      WRITE(IOUT,116)((DASH,K=1,NDASH),PLUS,L=1,4)
116   FORMAT(25X,'----',8X,'+',80A1)
      ITOTAL=TOTAL
      WRITE(IOUT,117)ITOTAL,((BLANK,K=1,NDASH),ARROW,L=1,4)
117   FORMAT(25X,I4,9X,80A1)
      IF(IOUT.EQ.21) WRITE(IOUT,118)(CHAR(L),L=1,4)
118   FORMAT(38X,4(15X,F5.1))
      IF(IOUT.NE.21) WRITE(IOUT,119)(CHAR(L),L=1,4)
119   FORMAT(38X,5F5.1)
      WRITE(IOUT,121)
121   FORMAT(54X,'PRECENTAGE')
27    CONTINUE
      RETURN
      END
C                                                  **** STAT PACK ****
C
C     ROUTINE TO OUTPUT SIZE OF OVERLAYS IN DECIMLE WORDS.
C     CALLING SEQUENCE: CALL SIZZ
C
C     TO OFTEN IN PAST THE SIZE OF INDIVIDUAL OVERLAYS WAS LEFT TO 
C     GUESS.  NOW SIZES OF OVERLAYS ARE BECOMING MORE STANDARD.  NORM
C     GRANT (W. M. U.) IS RESPONSIBLE FOR THE SIZE MACRO ROUTINE.
C
      SUBROUTINE SIZZ
      DO 6 I=1,16
      L=I
      CALL SIZE (L,LOC)
      TYPE 7,I,LOC
7     FORMAT(1X,'OVERLAY',I3,3X,I9)
6     CONTINUE
      RETURN
      END
      SUBROUTINE MABNK(NV,NC,MV,MC,DATA,NAMES)
      DIMENSION XDAT(125),DATA(MC,MV),NAMES(1),INPUT(10)
      DIMENSION DATED(2),NNS(18,6)
      COMMON /DEV/ ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      DOUBLE PRECISION FILN
      EQUIVALENCE (XDAT,NNS)
      IDSK=1
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT(' BANK NAME? ',$)
      READ(ICC,3,END=60) INPUT
      IF(INPUT(1).EQ.'!') RETURN
3     FORMAT(10A1)
      IF(INPUT(1).NE.' ') GO TO 24
      WRITE(IDLG,23)
23    FORMAT(' FOR ADDITIONAL INFORMATION TYPE HELP')
      GO TO 1
24    IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'L').
     1OR.(INPUT(4).NE.'P')) GO TO 22
      WRITE(IDLG,21)
21    FORMAT(' THE "MABNK" COMMAND IS USED TO CREATE A DATA BANK.'/
     1' DATA WILL BE STORED IN A BINARY FILE - DO NOT TRY TO PRINT'/
     2' THE BANK FILE!  NAMES WILL BE STORED ALONG WITH THE DATA.'/
     3' IF NAMES HAVE NOT BEEN ASSIGNED TO VARIABLES, V1,V2,V2,ETC.'/
     4' WILL BE USED. THE BANK NAME MAY BE UP TO 6 CHARACTERS WITHOUT'/
     5' AN EXTENSION, ".BNK" WILL BE ADDED.  TO ACCESS A STORED DATA'/
     6' BANK USE THE "ACBNK" COMMAND'/)
      GO TO 1
22    I=1
4     IF(INPUT(I).EQ.' ') GO TO 10
      IF(INPUT(I).EQ.'.') GO TO 8
C     CHECK FOR RIGHT BRACKET
      IF(INPUT(I).NE."555004020100) GO TO 6
      WRITE(IDLG,5)
5     FORMAT(' NO PROJECT-PROGRAMMER NUMBER NECESSARY '/
     1' - YOUR NUMBER ASSUMED.'/)
      GO TO 10
6     I=I+1
      IF(I.LE.6) GO TO 4
      WRITE(IDLG,7)
7     FORMAT(' MAXIMUM OF 6 CHARACTERS FOR BANK NAME')
      GO TO 1
8     WRITE(IDLG,9)
9     FORMAT(' NO EXTENSION NECESSARY ".BNK" WILL BE ADDED'/)
10    INPUT(I)='.'
      INPUT(I+1)='B'
      INPUT(I+2)='N'
      INPUT(I+3)='K'
      I=I+4
      IF(I.EQ.10) GO TO 12
      DO 11 J=I,10
11    INPUT(J)=' '
12    ENCODE(10,3,FILN) INPUT
      IPRO="055
13    IF(ICC.NE.2) WRITE(IDLG,14)
14    FORMAT(' WHAT PROTECTION? ',$)
      READ(ICC,3,END=60) INPUT
      IF(INPUT(1).EQ.'!') RETURN
      IF((INPUT(1).NE.'H').OR.(INPUT(2).NE.'E').OR.(INPUT(3).NE.'L').
     1OR.(INPUT(4).NE.'P')) GO TO 16
      WRITE(IDLG,15)
15    FORMAT(' ENTER THE 3 DIGIT PROTECTION CODE YOU DESIRE'/
     1' IF A RETURN IS GIVEN THE PROTECTION WILL BE ASSUMED 055')
      GO TO 13
16    IF(INPUT(1).EQ.' ') GO TO 20
      DO 17 J=1,3
      IF((INPUT(J).GE.'0').AND.(INPUT(J).LE.'7')) GO TO 17
      WRITE(IDLG,18)
18    FORMAT(' ALL DIGITS OF THE PROTECTION MUST BE BETWEEN 0 AND 7')
      GO TO 13
17    CONTINUE
      ENCODE(3,3,ISAV)(INPUT(I),I=1,3)
      DECODE(3,19,ISAV) IPRO
19    FORMAT(O3)
20    OPEN(UNIT=IDSK,FILE=FILN,ACCESS='RANDOM',DEVICE='DSK',
     1RECORD SIZE=126,MODE='BINARY',PROTECTION=IPRO)
      CALL GETPPN(IPROJ,IPROG)
      CALL DATE(DATED)
      I=1
      J=0
      VER='V2   '
      IREC=1
      WRITE(IDSK#IREC)NV,NC,I,DATED,IPROJ,IPROG,VER,(J,K=1,117)
C
C     PUT DATA OUT
C
      DO 30 I=1,NV
      DO 31 J=1,NC,125
      JEND=J+124
      IF(JEND.GT.NC) JEND=NC
      DO 40 K=J,JEND
40    XDAT(K-J+1)=DATA(K,I)
      IF(JEND.EQ.(J+124)) GO TO 32
      LAST=JEND-J+2
      DO 41 K=LAST,125
41    XDAT(K)=-9999E-20
32    IREC=IREC+1
31    WRITE(IDSK#IREC) XDAT
30    CONTINUE
      DO 50 I=1,NV,6
      IEND=I+5
      IF(IEND.GT.NV) IEND=NV
      DO 51 K=I,IEND
      LNAM=NAMES(K)
      DECODE(5,3,LNAM) (INPUT(J),J=1,5)
      IF((INPUT(1).LE.'Z').AND.(INPUT(1).GE.'A')) GO TO 56
      INPUT(1)='V'
      ENCODE(4,52,LNAM) K
52    FORMAT(I4)
      DECODE(4,3,LNAM)(INPUT(J),J=2,5)
54    IF(INPUT(2).NE.' ') GO TO 55
      DO 53 J=3,5
53    INPUT(J-1)=INPUT(J)
      INPUT(5)=' '
      GO TO 54
55    ENCODE(5,3,LNAM) (INPUT(J),J=1,5)
56    NNS(1,K-I+1)=LNAM
      NNS(10,K-I+1)=0
      ENCODE(40,57,NNS(2,K-I+1)) K,FILN
57    FORMAT('CREATED IN STP - VAR ',I3,' BANK ',2A5)
51    CONTINUE
      IREC=IREC+1
50    WRITE(IDSK#IREC) XDAT
      END FILE (IDSK)
      CALL RELEAS (IDSK)
60    RETURN
      END
C                                          *** STAT PACK ***
C     SUBROUTINE FOR CALCULATING PROBABILITIES ASSOC. WITH T
C     F, AND CHI SQUARES.
C     CALLING SEQUENCE: CALL PROB
C
C     PROGRAM MAKES USE OF FUNCTION FISHER.  USER HAS OPTION
C     OF ENTERING STRING COMMANDS OR TYPE OF PROBABILTIY
C     DESIRED, ALLOWING THE MACHINE TO ASK FOR THE INDIVIDUAL
C     VALUES IT NEEDS.  FOR CHI SQUARE PROBABILITIES THE
C     SUBROUTINES CHIPRB AND CUNO WERE ADDED.
C
      SUBROUTINE PROB
      COMMON /DEV/ICC,IDATA,IOUT,IDLG,IDSK
      COMMON /PRNT/ LINPP,ICOPS,RUNPRG
      COMMON/EXTRA/HEDR(70),NSZ
      IF(IOUT.EQ.21) CALL PRNTHD
      LINES=2
      IF(ICC.NE.2) WRITE(IDLG,60)
60    FORMAT(' "?" INDICATES PROGRAM IS WAITING FOR  INSTRUCT.')
      IF(IOUT.NE.21) WRITE(IOUT,5566)(HEDR(I),I=1,NSZ)
5566  FORMAT('1',70A1)
1     IF(ICC.NE.2) WRITE(IDLG,2)
2     FORMAT('0? ',$)
      READ(ICC,3,END=55) HELP
3     FORMAT(A4,20X)
      IF(HELP.NE.'HELP') GO TO 5
      WRITE(IDLG,4)
4     FORMAT('0PROBABILITIES ASSOCIATED WITH T TESTS, F TESTS',
     1', AND CHI SQUARES'/' ARE GIVEN IN THIS SECTION.  EACH MAY BE',
     2' OBTAINED IN TWO DIFFERENT'/' WAYS: EITHER BY USE OF A ',
     3'COMMAND STRING (WHERE ALL PERTINENT DATA IS'/' GIVEN)',
     4', OR BY SPECIFYING THE TYPE OF PROBABILITY DESIRED, AND'/
     5' ALLOWING THE PROGRAM TO ASK FOR THE NECESSARY INFORMATION',
     6'.  THE'/' THREE COMMAND STRINGS ARE:'/'0CHI SQR - CHI<CHI.',
     7' SQ. VALUE>,<DEGREES OF FREEDOM>'/' T TESTS - T',
     8'<T SCORE VALUE>,<DEGREES OF FREEDOM>'/' F TEST  - F<F SCORE',
     9' VALUE>,<NUM. DEG. OF FREED.>,<DEN. DEG. OF FREED.>'/
     1' WHERE <> INDICATES THE NUMERIC VALUE SPECIFIED IS TO BE',
     2' INSERTED.'/'0THE OTHER METHOD CONSISTS OF GIVING THE CODE',
     3' (T,F, OR CHI) AND'/' ALLOWING THE PROGRAM TO ASK FOR',
     4' THE NECESSARY VALUES.  A "?"'/' INDICATES THE PROGRAM IS ',
     5'WAITING FOR A COMMAND.  TO EXIT FROM THE'/' PROBABILITY',
     6' SECTION TYPE "EXIT".')
      GO TO 1
5     IF(HELP.EQ.'EXIT') RETURN
      IF(HELP.EQ.'     ')RETURN
      IF(HELP.EQ.'!') RETURN
      REREAD 6,HELP,CHI,DOF
6     FORMAT(A3,2F)
      IF(HELP.NE.'CHI') GO TO 20
C     CHI SQR
      NDG=DOF
      IF(NDG.LE.0) GO TO 9
7     CALL CHIPRB(NDG,CHI,PROB,IERR)
      IF(IERR.EQ.1) WRITE(IDLG,107)
107   FORMAT(' NO PROBABILITY CALCULATED -'
     1,' DF OR CHI SQ. OUTSIDE LIMITS')
      IF(IERR.EQ.1) GO TO 1
      IF(IOUT.NE.21) GO TO 15
      LINES=LINES+1
      IF(LINES.LE.LINPP) GO TO 15
      CALL PRNTHD
      LINES=3
15    WRITE(IOUT,8)CHI,NDG,PROB
8     FORMAT(' THE PROB. FOR A CHI SQ OF',F7.3,' WITH',I5,' DEGREES',
     1' OF FREEDOM IS',F6.3)
      GO TO 1
9     IF(CHI.GT.0) GO TO 10
      IF(ICC.NE.2) WRITE(IDLG,11)
11    FORMAT('+WHAT IS THE VALUE OF THE CHI. SQ.? ',$)
      READ (ICC,12) CHI
12    FORMAT(F)
10    IF(ICC.NE.2) WRITE(IDLG,13)
13    FORMAT('+HOW MANY DEGREES OF FREEDOM? ',$)
      READ (ICC,12) DOF
      NDG=DOF
      DOF=NDG
      IF(NDG.GT.0) GO TO 7
      WRITE(IDLG,14)
14    FORMAT(' THE NUMBER OF DEGREES OF FREEDOM MUST BE GREATER THAN',
     1' 0.'/)
      GO TO 10
C     T TEST
20    REREAD 21,HELP,T,DOF,DOF2
21    FORMAT(A1,3F)
      IF(HELP.NE.'T') GO TO 30
      NDG=DOF
      IF(NDG.LE.0) GO TO 24
22    TSQ=T**2
      PROB2=FISHER(1,NDG,TSQ)
      PROB=.5*PROB2
      IF(IOUT.NE.21) GO TO 28
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 28
      CALL PRNTHD
      LINES=4
28    WRITE(IOUT,23) T,NDG,PROB,PROB2
23    FORMAT(' THE PROB. FOR A T OF',F7.3,' WITH',I5,' DEGREE',
     1'S OF FREEDOM IS:',/10X,'ONE TAILED',F7.4,'; TWO TAILED',F7.4)
      GO TO 1
24    IF(T.GT.0) GO TO 26
      IF(ICC.NE.2) WRITE(IDLG,25)
25    FORMAT('+WHAT IS THE VALUE OF THE T? ',$)
      READ(ICC,12) T
26    IF(ICC.NE.2) WRITE(IDLG,13)
      READ(ICC,12)DOF
      NDG=DOF
      DOF=NDG
      IF(NDG.GT.0) GO TO 22
      WRITE(IDLG,14)
      GO TO 26
C     F TEST
30    IF(HELP.NE.'F') GO TO 50
      F=T
      NDG=DOF
      DOF=NDG
      NDG2=DOF2
      DOF2=NDG2
      IF(NDG.LE.0) GO TO 33
      IF(NDG2.LE.0) GO TO 37
31    PROB=FISHER(NDG,NDG2,F)
      IF(IOUT.NE.21) GO TO 39
      LINES=LINES+2
      IF(LINES.LE.LINPP) GO TO 39
      CALL PRNTHD
      LINES=4
39    WRITE(IOUT,32) F,NDG,NDG2,PROB
32    FORMAT(' PROB FOR  AN F OF',F7.3,' WITH',I5,' DEGREES OF',
     1' FREEDOM IN THE NUMERATOR'/10X,'AND',I5,' DEGREES OF FREEDOM',
     2' IN THE DENOMINATOR IS',F7.4)
      GO TO 1
33    IF(F.GT.0) GO TO 35
      IF(ICC.NE.2) WRITE(IDLG,34)
34    FORMAT('+WHAT IS THE VALUE OF THE F? ',$)
      READ(ICC,12)F
35    IF(ICC.NE.2) WRITE(IDLG,36)
36    FORMAT('+HOW MANY DEGREES OF FREEDOM IN THE NUMERATOR? ',$)
      READ (ICC,12) DOF
      NDG=DOF
      DOF=NDG
      IF(NDG.GT.0) GO TO 37
      WRITE(IDLG,14)
      GO TO 35
37    IF(NDG2.GT.0) GO TO 31
      IF(ICC.NE.2) WRITE(IDLG,38)
38    FORMAT('+HOW MANY DEGREES OF FREEDOM IN THE DENOMINATOR? ',$)
      READ(ICC,12) DOF2
      NDG2=DOF2
      DOF2=NDG2
      IF(NDG2.GT.0) GO TO 31
      WRITE(IDLG,14)
      GO TO 37
50    WRITE(IDLG,51)
51    FORMAT(' THE COMMAND JUST GIVEN DOES NOT EXIST')
      GO TO 1
55    RETURN
      END
C                                    *** STAT PACK ***
C     SUBROUTINE USED FOR DETERMINING CHI SQUARE PROBABILITIES
C     CALLING SEQUENCE: CALL CHIPRB(K,X,Y,IERR)
C     WHERE K - NUMBER OF DEGREES OF FREEDOM
C           X - CHI SQUARE VALUE
C           Y - PROBABILITY ASSOCIATED WITH CHI SQUARE
C           IERR - ERROR WAS ENCOUNTERED WHEN ATTEMPTING TO CALCULATE
C                  PROBABILITY
C
C     ROUTINE WAS WRITTEN BY CHARLES NAGY OF WESTERN, AND ADAPTED
C     TO STP BECAUSE THE FISHER ROUTINE WOULD NOT GIVE ENOUGH 
C     ACCURACY FOR THE PSYCHOLOGISTS.  CALLS SUBROUTINE CUNO.
C
      SUBROUTINE CHIPRB(K,X,Y,IERR)
      DIMENSION F(25)
      F(1)=.5
      F(2)=.598706326
      F(3)=.691462461
      F(4)=.773372648
      F(5)=.841344746
      F(6)=.894350226
      F(7)=.933192799
      F(8)=.959940843
      F(9)=.977249868
      F(10)=.987775527
      F(11)=.993790335
      F(12)=.997020237
      F(13)=.998650102
      F(14)=.999422975
      F(15)=.999767371
      F(16)=.999911583
      F(17)=.999968329
      F(18)=.999989311
      F(19)=.999996602
      F(20)=.999998983
      F(21)=.999999713
      F(22)=.999999924
      F(23)=.999999981
      F(24)=.999999996
      F(25)=.999999999
      IERR=0
      Y=0
      IF((K.LE.0).OR.(K.GT.100)) IERR=1
      IF(X.GT.141) IERR=1
      IF(IERR.EQ.1) RETURN
      IF(X.LE.0) GO TO 13
      IF(K.GE.4) GO TO 4
      GO TO (1,2,3),K
1     P=SQRT(X)
      CALL CUNO(P,S,IERR,F)
      IF(IERR.EQ.1) RETURN
      Y=2.*S-1
      GO TO 13
2      Y=1.-(1./EXP(X/2.))
      GO TO 13
3     P=SQRT(X)
      CALL CUNO(P,S,IERR,F)
      IF(IERR.EQ.1) RETURN
      Y=(2.*S-1)-P/(1.25331414*EXP(X/2.))
      GO TO 13
4     M=K/2
      IF(K.EQ.2*M) GO TO 6
      P=SQRT(X)
      CALL CUNO(P,S,IERR,F)
      IF(IERR.EQ.1) RETURN
      Y=2.*S-1.
      S=X/2.
      C=1./(.62665707*P*EXP(S))
      P=S
      T=.5
      GO TO 7
6     C=1./EXP(X/2.)
      Y=1.-C
      S=0
      P=1
      T=0
7     DO 8 I=1,M-1
      T=T+1
      P=P*(X/(T*2.))
8     S=S+P
      Y=Y-C*S
13    Y=1.-Y
      END
C                                      *** STAT PACK ***
C     SUBROUTINE USED IN FINDING PROB FOR CHI SQUARE
C     CALLING SEQUENCE: CALL CUNO(X,Y,IERR,F)
C     ORIGINALLY WRITTEN BY CHARLES NAGY OF WMU.
C
      SUBROUTINE CUNO(X,Y,IERR,F)
      DIMENSION F(1)
      W=X
      IF(W.LT.0) W=-W
      IF(W.LE.6.125) GO TO 2
1      Z=1
      GO TO 7
2     K=INT(4.*W)+1
      A=.25*(K-1.)
      IF(W-A) 10,3,4
3     Z=F(K)
      GO TO 7
4     IF(W-(A+.125))6,6,5
5     K=K+1
      A=A+.25
6     H=W-A
      ASQ=A*A
      C1=((-ASQ+10.)*ASQ-15.)*A
      C2=(6.*ASQ-36.)*ASQ+18.
      C3=(-30.*ASQ+90.)*A
      C4=120.*(ASQ-1.)
      C5=-360.*A
      C6=(((((C1*H+C2)*H+C3)*H+C4)*H+C5)*H+720.)*H
      Z=F(K)+C6/(720.*SQRT(6.28318531*EXP(ASQ)))
7     Y=Z
      IF(X.LT.0.) Y=1.-Z
      RETURN
10    IERR=1
      RETURN
      END

