Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - decus/20-0137/bstat/bstat.for
There is 1 other file named bstat.for in the archive. Click here to see a list.
C	WESTERN MICHIGAN UNIVERSITY
C	BSTAT.F4 (FILE NAME ON LIBRARY DECTAPE)
C	BSTAT, 1.1.2 (CALLING NAME, SUBLST #)
C	BASIC STATISTICS (MEANS, ST. DEVS., VAR., SUM, SUM SQ.,
C	 COEF. OF VAR., ST. ERROR OF MEAN, 3RD AND 4TH MOM., MAX,
C	 MIN, RANGE, ONE WAY ANOVA, T AND F TESTS)
C	PROGRAMMED BY SAM ANEMA PATTERNED AFTER A PROGRAM GIVEN
C	 BY WAYNE STATE UNIVERSITY, PUT ON SYSTEM BEFORE 11/30/71
C	STATISTICAL CONSULTANT - DR. M. STOLINE
C	FORWMU PROGS. USED:  TTYPTY, EXISTS, PRINTS, DEVCHG
C	APLIB PROGS. USED:  GETFOR, IO, FISHER
C	LIBRARY DECTAPE PROGS. USED:  USAGE.MAC
C	INTERNAL SUBR. USED:  SETUP, FUNC
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
c	Spurious comma removed from ENCODE statement.
c		Paul T. Robinson, Wesleyan Univ, 18 Oct 1980
c
C---------------IFMT(96) - NO MORE THAN 6 FORMAT CARDS
C---------------IDENT(16) - 80 CHARS. TO IDENTIFY OUTPUT.
C---------------XDIF(50,50) - NO MORE THAN 50 VARS. ENTERING DATA BY VARS.
C---------------100 COVERS LIMITATIONS FOR BOTH METHODS OF INPUT
      DIMENSION IFMT(96),IDENT(16),XDIF(50,50)
      DIMENSION SUMC(100),SUMF(100), BIG(100), SMA(100)
      DIMENSION AMEAN(100),SD(100)
      DIMENSION N(100),V(100),SUM(100),SUMSQ(100)
C---------------JNAM - SPACE FOR FILENAME AND EXT
	DOUBLE PRECISION JNAM
      FIRST=.TRUE.
      INP=2
      IOUT=3
      IND=-4
      NOUTD=-1
      WRITE(NOUTD,2)
2     FORMAT(' ',31X,'WMU'//23X,'BASIC STATISTICS PROGRAM'///)
C	CALL USAGE('BSTAT')
C---------------TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICODER)
C---------------IDV1 IS RETURNED.  OTHER ARGS. ARE INPUT.  1 CAUSES
C--------------- OUTPUT? TO PRINT
      CALL IO(IOUT,IDV1,NOUTD,IND,1,ICODER)
703	CLOSE(UNIT=1)
	CALL JOBNUM(JNUM)
	JNUM=100+JNUM
	ENCODE(10,1709,JNAM)JNUM
1709	FORMAT('JOB',I3,'.DAT')
	CALL DELETE(JNAM)
C---------------IDV2 IS RETURNED.  OTHER ARGS. ARE INPUT.
C--------------- 0 CAUSES INPUT? TO PRINT
      CALL IO(INP,IDV2,NOUTD,IND,0,ICODER)
	IF(IDV2.EQ.'TTY')OPEN(UNIT=1,DEVICE='DSK',FILE=JNAM,
     1ACCESS='SEQOUT',RECORDSIZE=80,MODE='ASCII')
      TOT=0.
      SSTOT=0.
      WRITE(NOUTD,705)
705   FORMAT(' TYPE 1 FOR ANALYSIS BY SAMPLE, 2 FOR ANALYSIS BY VARIABLE
     1'/)
      READ(IND,20)NHOW
      IF(NHOW.EQ.1)GO TO 720
      VR='VAR.'
709   WRITE(NOUTD,10)
10    FORMAT(' ','HOW MANY VARIABLES?'/)
      READ(IND,20)NVAR
20    FORMAT(I)
      IF(NVAR.GT.50)GO TO 706
      DO 45 I=1,NVAR
      N(I)=0
      SUM(I)=0
      SUMSQ(I)=0.
      SUMC(I)=0.
      SUMF(I)=0.
      BIG(I)=-1.E30
      SMA(I)=1.E30
      DO 45 J=1,NVAR
      XDIF(I,J)=0.0
45    CONTINUE
      NTOT=0
      GO TO 707
706   WRITE(NOUTD,708)
708   FORMAT(' ANALYSIS BY VARIABLE MAY NOT BE DONE ON MORE THAN 50 VARI
     1ABLES'/)
      GO TO 709
707   WRITE(NOUTD,21)
21    FORMAT(' ENTER IDENTIFICATION'/)
      READ(IND,701)(IDENT(I),I=1,16)
701   FORMAT(16A5)
C---------------IFMT, ISTD ARE RETURNED. OTHER ARGS. ARE INPUT.
C--------------- 96=NO. OF OBJ. TIME FORMAT WORDS (6 LINES)
C--------------- 2 MEANS F-TYPE FORMAT ONLY.
      CALL GETFOR(IND,NOUTD,IFMT,ISTD,96,2)
      IF(ISTD.EQ.1)IFMT(1)='(10F)'
      IDC=INP
      IF(IDV2.NE.'TTY')GO TO 805
      IF(ISTD.EQ.1)WRITE(NOUTD,740)
      IF(ISTD.NE.1)WRITE(NOUTD,741)
      CALL SETUP(100000,NOUTD,IND,JNAM)
      IDC=1
	CLOSE(UNIT=1)
	OPEN(UNIT=1,DEVICE='DSK',FILE=JNAM,
     1ACCESS='SEQIN',RECORDSIZE=80,MODE='ASCII')
C---------------READ DATA AND CALC. FOR INPUT BY VARS.
805	READ(IDC,IFMT,END=899)(V(I),I=1,NVAR)
807   DO 809 I=1,NVAR
      N(I)=N(I)+1
      SUM(I)=SUM(I)+V(I)
      SUMSQ(I)=SUMSQ(I)+V(I)**2
      SUMC(I)=SUMC(I)+V(I)**3
      SUMF(I)=SUMF(I)+V(I)**4
      NTOT=NTOT+1
      TOT=TOT+V(I)
      SSTOT=SSTOT+V(I)*V(I)
      IF(SMA(I).GT.V(I))SMA(I)=V(I)
      IF(BIG(I).LT.V(I))BIG(I)=V(I)
      DO 710 J=I,NVAR
710   XDIF(I,J)=XDIF(I,J)+(V(I)-V(J))**2
809   CONTINUE
      GO TO 805
720   VR='SAMP'
      WRITE(NOUTD,21)
      READ(IND,701)(IDENT(I),I=1,16)
      WRITE(NOUTD,721)
721   FORMAT(' HOW MANY SAMPLES?'/)
      READ(IND,20)NVAR
      DO 451 I=1,NVAR
      N(I)=0
      SUM(I)=0
      SUMSQ(I)=0.
      SUMC(I)=0.
      SUMF(I)=0.
      BIG(I)=-1.E30
      SMA(I)=1.E30
      DO 451 J=1,NVAR
      XDIF(I,J)=0.0
451    CONTINUE
      NTOT=0
801   WRITE(NOUTD,30)
30    FORMAT(' WHAT ARE THE SAMPLE SIZES?'/)
      READ(IND,40)(N(J),J=1,NVAR)
40    FORMAT(10I)
      LINES=0
      DO 972 J=1,NVAR
      FLIN=FLOAT(N(J))/10.0+.999
972   LINES=LINES+IFIX(FLIN)
      CALL GETFOR(IND,NOUTD,IFMT,ISTD,96,2)
      IF(ISTD.EQ.1)IFMT(1)='(10F)'
      IDC=INP
      IF(IDV2.NE.'TTY')GO TO 739
      IF(ISTD.EQ.1)WRITE(NOUTD,740)
740   FORMAT(' ENTER DATA (AT MOST 10 PER LINE)'/)
      IF(ISTD.NE.1)WRITE(NOUTD,741)
741   FORMAT(' ENTER DATA'/)
C---------------ALLOW USER TO MAKE CORRECTIONS IF INPUT IS BY TTY
      CALL SETUP(LINES,NOUTD,IND,JNAM)
	KINC=LINES/NVAR
	KX=1-KINC
	CLOSE(UNIT=1)
	OPEN(UNIT=1,DEVICE='DSK',FILE=JNAM,
     1ACCESS='SEQIN',RECORDSIZE=80,MODE='ASCII')
      IDC=1
C---------------READ DATA AND CALC. FOR INPUT BY SAMPLE.
739   DO 100 I=1,NVAR
      ITR=0
      NCY=0
      IF(N(I)-100)60,60,70
60    JEND=N(I)
      GO TO 80
70    JEND=100
80	READ(IDC,IFMT)(V(J),J=1,JEND)
      DO 90 J=1,JEND
      SUM(I)=SUM(I)+V(J)
      SUMSQ(I)=SUMSQ(I)+V(J)*V(J)
      SUMC(I)=SUMC(I)+V(J)**3
      SUMF(I)=SUMF(I)+V(J)**4
      TOT=TOT+V(J)
      SSTOT=SSTOT+V(J)*V(J)
      IF(ITR)110,120,110
120   BIG(I)=V(J)
      SMA(I)=V(J)
      ITR=1
110   IF(BIG(I)-V(J))140,130,130
140   BIG(I)=V(J)
130   IF(SMA(I)-V(J))90,90,160
160   SMA(I)=V(J)
90    CONTINUE
      NCY=NCY+JEND
      IF(NCY-N(I))170,180,180
170   NT=N(I)-NCY
      IF(NT-100)190,190,200
190   JEND=NT
      GO TO 80
200   JEND=100
      GO TO 80
180   NTOT=NTOT+N(I)
100   CONTINUE
C---------------BOTH METHODS OF INPUT GO THIS WAY
899   WRITE(IOUT,23)(IDENT(I),I=1,16)
23     FORMAT(/' ',16A5)
      DO 195 I=1,NVAR
      XN=N(I)
      AMEAN(I)=SUM(I)/XN
195   SD(I)=SQRT((SUMSQ(I)-XN*(AMEAN(I)**2))/(XN-1.))
182   WRITE(NOUTD,210)
210   FORMAT(' ','WOULD YOU LIKE OPTION 1?'/)
      READ(IND,220)NOPT1
220   FORMAT(A3)
      IF(NOPT1-'YES')240,230,240
C       OPTION 1
230   WRITE(IOUT,250)VR
250   FORMAT(///30X,'OPTION 1'//' ',1X,A4,3X,'SIZE',11X,
     1'MEAN',10X,'VARIANCE',6X,'STANDARD DEVIATION'//)
      DO 300 I=1,NVAR
      VAR=SD(I)**2
300   WRITE(IOUT,260)I,N(I),AMEAN(I),VAR,SD(I)
260   FORMAT(' ',I4,I7,F19.5,F18.5,F19.5/)
      WRITE(IOUT,330)VR
330   FORMAT(/' ',28X,'SUM OF SQUARES',3X,'COEFFICIENT',5X,'STANDARD',
     1/' ',1X,A4,3X,'SIZE',8X,'SUM',4X,'OF OBSERVATIONS',2X,
     2'OF VARIATION',2X,'ERROR OF MEAN'//)
      DO 310 I=1,NVAR
      XN=N(I)
      SEM=SD(I)/SQRT(XN)
      CV=(SD(I)/AMEAN(I))*100.
310   WRITE(IOUT,320)I,N(I),SUM(I),SUMSQ(I),CV,SEM
320   FORMAT(' ',I4,I7,2X,2G,F13.3,F12.4/)
240   WRITE(NOUTD,241)
241   FORMAT(/' ','WOULD YOU LIKE OPTION 2?'/)
      READ(IND,220)NOPT1
      IF(NOPT1-'YES')370,342,370
342   WRITE(IOUT,305)VR
305   FORMAT(///30X,'OPTION 2'//' ',A4,9X,'3RD MOMENT',
     14X,'SKEWNESS',10X,'4TH MOMENT',4X,'KURTOSIS'//)
      DO 350 I=1,NVAR
      XN=N(I)
      XMOM3=SUMC(I)/XN-(3.*AMEAN(I)*SUMSQ(I))/XN+2.*AMEAN(I)**3
      XMOM4=SUMF(I)/XN-(4.*AMEAN(I)*SUMC(I))/XN+(6.*(AMEAN(I)**2)
     1*SUMSQ(I))/XN-3.*(AMEAN(I)**4)
      SKEWN=XMOM3/((XN-1.)*SD(I)**2/XN)**1.5
      XKUR=XMOM4/((XN-1.)*SD(I)**2/XN)**2
      WRITE(IOUT,340)I,XMOM3,SKEWN,XMOM4,XKUR
340   FORMAT(' ',I3,F20.5,F12.3,F20.5,F12.3/)
350   CONTINUE
      WRITE(IOUT,371)VR
371   FORMAT(/' ',A4,6X,'MAXIMUM VALUE',7X,'MINIMUM VALUE',
     113X,'RANGE'//)
      DO 372 I=1,NVAR
      RANGE=BIG(I)-SMA(I)
372   WRITE(IOUT,373)I,BIG(I),SMA(I),RANGE
373   FORMAT(' ',I3,3F20.5/)
C---------------BYPASS OPTION 3 IF INPUT BY VAR.
370   IF(NHOW.EQ.2)GO TO 490
      WRITE(NOUTD,401)
401   FORMAT(' ','WOULD YOU LIKE OPTION 3?'/)
      READ(IND,220)NOPT1
      IF(NOPT1-'YES')490,400,490
400   WRITE(IOUT,410)
410   FORMAT(///30X,'OPTION 3')
      UP=0.
      ALG=0.
      XNTOT=NTOT
      C=0.
      SST=0.
      SSB=0.
      DO 430 I=1,NVAR
      XN=N(I)
      C=C+SUM(I)
      SST=SST+SUMSQ(I)
      UP=UP+(1./(XN-1))
      ALG=ALG+(XN-1)*ALOG(SD(I)**2)
430   SSB=SSB+(SUM(I)**2/XN)
      SST=SST-(C**2)/XNTOT
      SSB=SSB-(C**2)/XNTOT
      SSW=SST-SSB
      NUM=NVAR-1
      NOM=NTOT-NVAR
      NOT=NTOT-1
      XNUM=NUM
      XNOM=NOM
      BMS=SSB/XNUM
      WMS=SSW/XNOM
      F=BMS/WMS
      PROB=FISHER(NUM,NOM,F)
      CHI=(1./((UP-1./XNOM)/(3.*XNUM)+1.))*(XNOM*ALOG(WMS)-ALG)
      WRITE(IOUT,420)CHI,NUM
420   FORMAT(//27X,'BARTLETTS TEST'/10X,'CHI SQUARE = ',F8.3,10X,
     1'D.F. = ',I4)
      WRITE(IOUT,470)
470   FORMAT(//20X,'ONE WAY ANALYSIS OF VARIANCE'//2X,'SOURCE',5X,
     1'SUM OF SQ.',5X,'D.F.',4X,'MEAN SQ.',8X,'F',4X,'SIGNIFICANT AT'/)
      WRITE(IOUT,450)SSB,NUM,BMS,F,PROB,SSW,NOM,WMS,SST,NOT
450   FORMAT(2X,'BETWEEN',F16.5,I7,F15.5,F8.2,F10.3/'  WITHIN',F17.5,I7,
     1F15.5/3X,'TOTAL',F17.5,I7)
C---------------BOTH METHODS OF INPUT GO THIS WAY
490   WRITE(NOUTD,502)
502   FORMAT(//' ','WOULD YOU LIKE OPTION 4?'/)
      READ(IND,220)NOPT1
      IF(NOPT1-'YES')703,501,703
501   WRITE(IOUT,500)
500   FORMAT(///30X,'OPTION 4'//24X,'STUDENT T AND F TESTS'//)
C---------------NHOW=2 MEANS INPUT BY VAR.
      IF(NHOW.EQ.2)GO TO 725
      WRITE(NOUTD,520)
520   FORMAT(' ','ARE POPULATION VARIANCES ASSUMED EQUAL?'/)
      READ(IND,530)YES
530   FORMAT(A3)
      IF(YES.EQ.'YES')WRITE(IOUT,1601)
1601  FORMAT(' FOR THE FOLLOWING T-VALUES THE POPULATION VARIANCES ARE 
     1ASSUMED'/' TO BE EQUAL.'/)
      IF(YES.EQ.'NO')WRITE(IOUT,602)
602   FORMAT(' FOR THE FOLLOWING T-VALUES THE POPULATION VARIANCES ARE
     1 ASSUMED'/' TO BE UNEQUAL.'/)
      GO TO 726
C---------------HERE IF INPUT BY VAR.
725   WRITE(IOUT,727)
727   FORMAT(' THE FOLLOWING T-TESTS ARE CORRELATED T-TESTS'/)
726   IF(FIRST)WRITE(IOUT,1603)
1603  FORMAT(' THE PROBABILITY ASSOCIATED WITH EACH T-VALUE IS OBTAINED
     1 BY ASSUMING'/' A TWO-TAILED TEST.  A ONE-TAILED TEST MAY BE DERI
     2VED BY HALVING THE '/' PROBABILITY VALUE GIVEN.'/)
      FIRST=.FALSE.
C---------------NHOW=2 MEANS INPUT BY VAR.
      IF(NHOW.EQ.2)GO TO 730
      WRITE(IOUT,525)
525   FORMAT(50X,'DF   DF'/' SAMP VS. SAMP',2X,'T-VALUE',3X,'DF',2X,
     1'PROB.',5X,'F-VALUE',2X,'NUM',2X,'DEN'/)
      GO TO 742
730   WRITE(IOUT,731)
731   FORMAT(' VAR VS. VAR',4X,'T-VALUE',3X,'DF',2X,'PROB.')
742   NVAR1=NVAR-1
      DO 599 I=1,NVAR1
      I1=I+1
      DO 599 J=I1,NVAR
      XI=N(I)
      XJ=N(J)
      IF(NHOW.EQ.2)GO TO 728
C---------------YES MEANS POP. VAR. ARE ASSUMED EQ.
      IF(YES-'YES')540,550,540
550   T=(AMEAN(I)-AMEAN(J))/SQRT(((SD(I)**2*(XI-1.)+SD(J)**2*
     1(XJ-1.))/(XI+XJ-2.))*(XI+XJ)/(XI*XJ))
      NDF=N(I)+N(J)-2
      GO TO 591
540   T=(AMEAN(I)-AMEAN(J))/SQRT(SD(I)**2/XI+SD(J)**2/XJ)
      NDF=((SD(I)**2/XI+SD(J)**2/XJ)**2)/((SD(I)**2/XI)**2/(XI+1.)
     1+(SD(J)**2/XJ)**2/(XJ+1.))-2.0
591   NFN=N(I)-1
      NFD=N(J)-1
      F=(SD(I)/SD(J))**2
      PROB=FISHER(1,NDF,T**2)
C---------------OUTPUT IF INPUT IS BY SAMPLE
      WRITE(IOUT,593)I,J,T,NDF,PROB,F,NFN,NFD
      GO TO 599
728   DEN=SQRT((XDIF(I,J)-XI*(AMEAN(I)-AMEAN(J))**2)/(XI*(XI-1.0)))
      IF(DEN.EQ.0.0)GO TO 751
      T=(AMEAN(I)-AMEAN(J))/DEN
      GO TO 752
751   WRITE(IOUT,753)
753   FORMAT(/' THE SD OF THE DIFFERENCES BETWEEN VARS. ',I4,' AND ',I4 
     1,' IS ZERO')
      GO TO 599
752   NDF=N(I)-1
      PROB=FISHER(1,NDF,T**2)
C---------------OUTPUT IF INPUT IS BY VAR.
      WRITE(IOUT,729)I,J,T,NDF,PROB
729   FORMAT(/1X,I3,I8,F10.3,I6,F7.3)
599   CONTINUE
593   FORMAT(/1X,I3,I8,F10.3,I6,F7.3,F11.3,I6,I5)
      GO TO 703
      END
C---------------ALLOW TO MAKE CORRECTIONS IF INPUT IS BY TTY
C---------------ALL ARGS. ARE INPUT.
      SUBROUTINE SETUP(LINES,NOUTD,IND,JNAM)
      DIMENSION X(16)
	DOUBLE PRECISION JNAM
      DO 3 J=1,LINES
1     READ(IND,10,END=20)(X(I),I=1,16)
10    FORMAT(16A5)
3     WRITE(1,10)(X(I),I=1,16)
20    CONTINUE
200   WRITE(NOUTD,21)
21    FORMAT(' OK? '/)
      READ(IND,22)HUH
22    FORMAT(A3)
      IF(HUH.EQ.'NO')GO TO 25
      IF(HUH.EQ.'YES')GO TO 26
      IF(HUH.NE.'OK')GO TO 200
      GO TO 26
25      K=1
	CLOSE(UNIT=1)
	OPEN(UNIT=1,DEVICE='DSK',FILE=JNAM,
     1ACCESS='RANDOM',RECORDSIZE=80,MODE='ASCII')
40    READ(1#K,10,END=41)(X(I),I=1,16)
      DO 27 I=1,16
      IF(X(I).NE.' ')GO TO 28
27    CONTINUE
28    DO 29 J=I,16
      IF(X(J).EQ.' ')GO TO 30
29    CONTINUE
      J=17
30    WRITE(NOUTD,31)K,(X(I),I=1,J-1)
31    FORMAT(' ',I3,2X,16A5)
      K=K+1
	IF(K.LE.LINES)GO TO 40
41    WRITE(NOUTD,45)
45    FORMAT(//)
101   WRITE(NOUTD,42)
42    FORMAT('+CORRECTIONS? ',$)
      READ(IND,43)K,(X(I),I=1,16)
43    FORMAT(I,16A5)
      IF(K.EQ.0)GO TO 100
      WRITE(1#K,10)(X(I),I=1,16)
      GO TO 101
100	CONTINUE
26    RETURN
      END