Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
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