Trailing-Edge
-
PDP-10 Archives
-
-
There are no other files named in the archive.
CBMD01R SIMPLE LINEAR REGRESSION FEBRUARY 4, 1965
C THIS IS A SIFTED VERSION OF BMD01R ORIGINALLY WRITTEN IN
C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
COMMON Y , X , FADJ , LTABLE
COMMON NNF , NNS , DIXON , XX , XY , YY
COMMON TXX , TXY , TYY , NDF1 , NDF2 , NDF3
COMMON NDF4 , NDF5 , TXXT , TXYT , TYYT , EXXT
COMMON EXYT , EYYT , SE , SEMS , SXXT , SXYT
COMMON SYYT , STPE , ADJT , ADJTMS , BM , FXX
COMMON BW , FYY , BT , FB
DIMENSION SY(600),SYY(600),NREP(600)
1,SX(600),SXX(600),SXY(600),Y(599),X(599),
2NS(499),NF(500),JCOMB(500),KOM(21),SSX(500),SSY(500),
3SSXX(500),SSYY(500),SSXY(500),NSREP(500),
4 SUMX(600),SUMXX(600),SUMY(600),SUMYY(600),SUMXY(600),IREPS(
5600),MAT(180)
DIMENSION NCX(8),CONX(8),NCY(8),CONY(8)
DOUBLE PRECISION TODE,COB,A123,B123,C123,D123,E123,F123,IA1,TRGN,
1ID,IE,IF1
EQUIVALENCE (D123,ID),(E123,IE),(F123,IF1)
DATA A123,B123,C123,D123,E123,F123/6HPROBLM,6HFINISH,6HSPECTG,
16HSAMSIZ,6HSUBSET,6HCOMSEL/
C
416 FORMAT('1BMD01R - SIMPLE LINEAR REGRESSION - REVISED ',
1'SEPTEMBER 2, 1968'/
240H HEALTH SCIENCES COMPUTING FACILITY,UCLA// )
C
100 FORMAT(2A6,I3,I1,I3,I2,I1,F6.0,40X,2I2)
101 FORMAT(A6,22I3)
102 FORMAT(12I6)
105 FORMAT(9(I2,F6.2))
107 FORMAT(1H0,16X,18HCOMSEL CARD NUMBERI3,64H MISPUNCHED OR OUT OF OR
1DER. PROGRAM PROCEEDS WITHOUT THIS CARD.)
108 FORMAT(A6,I2,21I3)
210 FORMAT(1H I4,I4,I3,5XF6.2,5XF6.2)
211 FORMAT(21H0SUBSET SPECIFICATION/)
212 FORMAT(17H NO START END/)
213 FORMAT(I5,I6,I6)
214 FORMAT(15H0COMBINATION NOI4//17H SUBSETS INCLUDED)
215 FORMAT(1H 20I5)
1002 FORMAT(18A4)
NTAPE=5
CALL USAGEB('BMD01R')
888 READ (5,100)TODE,COB,K,LTRANS,LSUB,LCOMB,LTABLE,DIXON,MTAPE,MFT
IF( TODE.EQ.B123) GO TO 401
400 IF( TODE.EQ.A123) GO TO 403
4020 TODE =A123
402 WRITE (6,404)TODE
404 FORMAT(47H0CONTROL CARDS INCORRECTLY ORDERED OR PUNCHED, A6,12H CA
1RD ERROR.)
401 IF(NTAPE-5)478,478,477
477 REWIND NTAPE
478 STOP
403 CALL TPWD(MTAPE,NTAPE)
CALL VFCHCK(MFT)
IERROR=0
405 IF(LSUB-501)406,4020,4020
406 IF(K*(K-600))3,4020,4020
3 NCARDS=(K+21)/22
L2=0
DO 4 J=1,NCARDS
L1=L2+1
L2=L2+22
IF(K-L2)33,34,34
33 L2=K
34 READ (5,101)IA1,(NREP(I),I=L1,L2)
IF(ID.EQ.IA1)GO TO 4
36 IERROR=-21
4 CONTINUE
IF(IERROR)37,38,38
37 TODE=D123
GO TO 402
38 WRITE (6,416)
WRITE (6,415)COB,K,LTRANS,MFT
415 FORMAT(14H PROBLEM CODE A6,/
115H NO. OF GROUPS I5,/
217H TRANSGENERATION I2,/
326H NO. OF VAR. FMT. CARD(S) I2)
NTY=0
NTX=0
KX=2*K
DO 222 I=1,K
SX(I)=0.0
SY(I)=0.0
SYY(I)=0.0
SXX(I)=0.0
SXY(I)=0.0
222 CONTINUE
MFT=MFT*18
READ (5,1002)(MAT(I),I=1,MFT)
WRITE(6,2221) (MAT(I), I=1,MFT)
2221 FORMAT(' THE VARIABLE FORMAT IS ',25A4,/(1H ,24X,25A4))
408 LTRANS=LTRANS+1
IF(LTRANS*(LTRANS-5))409,4020,4020
409 GO TO (410,411,412,413),LTRANS
411 READ (5,414)TRGN,NTX,(NCX(I),CONX(I),I=1,NTX)
1410 IF (TRGN .EQ.C123) GO TO 4105
GO TO 4025
412 READ (5,414)TRGN,NTY,(NCY(I),CONY(I),I=1,NTY)
GO TO 1410
413 READ (5,414)TRGN,NTX,(NCX(I),CONX(I),I=1,NTX)
IF( TRGN.EQ.C123) GO TO 412
GO TO 4025
C DIMENSION NCX(8),CONX(8),NCY(8),CONY(8)
414 FORMAT(A6,I1,8(I2,F6.0))
4025 TODE=C123
GO TO 402
4105 IF(-NTX)4026,4117,4117
4026 DO 4107 I=1,NTX
IF(NCX(I)*(NCX(I)-11))4107,4106,4106
4106 WRITE (6,4000)NCX(I)
4000 FORMAT(1H0,12X,27HTRANSGENERATION CODE NUMBER I3,36H IS ILLEGAL. P
1ROGRAM CANNOT PROCEED.)
IERROR=-22
4107 CONTINUE
4117 IF(-NTY)4108,4109,4109
4108 DO 4185 I=1,NTY
IF(NCY(I)*(NCY(I)-11))4185,4182,4182
4182 WRITE (6,4000)NCY(I)
IERROR=-22
4185 CONTINUE
4109 IF(-IERROR)410,410,401
410 DO 40 I=1,K
IERROR=I
N=NREP(I)
IF(-N)420,37,37
420 READ (NTAPE,MAT)(Y(L),X(L),L=1,N)
GO TO (13,7,10,11),LTRANS
7 CALL TRANS1(X,N,NTX,NCX,CONX,IERROR)
GO TO 6
10 CALL TRANS1(Y,N,NTY,NCY,CONY,IERROR)
GO TO 6
11 CALL TRANS1(X,N,NTX,NCX,CONX,IERROR)
GO TO 10
6 IF(IERROR)4025,13,13
13 DO 35 J=1,N
SX(I)=SX(I)+X(J)
SXX(I)=SXX(I)+X(J)*X(J)
SY(I)=SY(I)+Y(J)
SYY(I)=SYY(I)+Y(J)*Y(J)
SXY(I)=SXY(I)+X(J)*Y(J)
35 CONTINUE
40 CONTINUE
306 CALL COVAR(SX,SXX,SY,SYY,SXY,NREP,K)
WRITE (6,416)
CALL WRITE
IF(-LSUB)42,888,4020
42 NCARDS=(LSUB+10)/11
L2=0
DO 47 J=1,NCARDS
L1=L2+1
L2=L2+11
IF(LSUB-L2)43,44,44
43 L2=LSUB
44 READ (5,101)IA1,(NS(I),NF(I),I=L1,L2)
IF(IA1.EQ.IE)GO TO 47
434 IERROR=-22
47 CONTINUE
IF(IERROR)475,445,445
475 TODE=E123
GO TO 402
445 DO 45 I=1,LSUB
SSX(I)=0.0
SSXX(I)=0.0
SSY(I)=0.0
SSYY(I)=0.0
SSXY(I)=0.0
NSREP(I)=0
45 CONTINUE
DO 50 J=1,LSUB
NNS=NS(J)
NNF=NF(J)
M=0
DO 46 I=NNS,NNF
M=M+1
SSX(J)=SSX(J)+SX(I)
SSXX(J)=SSXX(J)+SXX(I)
SSY(J)=SSY(J)+SY(I)
SSYY(J)=SSYY(J)+SYY(I)
SSXY(J)=SSXY(J)+SXY(I)
NSREP(J)=NSREP(J)+NREP(I)
SUMX(M)=SX(I)
SUMXX(M)=SXX(I)
SUMY(M)=SY(I)
SUMYY(M)=SYY(I)
SUMXY(M)=SXY(I)
IREPS(M)=NREP(I)
46 CONTINUE
IF(LTABLE)4020,50,500
500 WRITE (6,416)
WRITE (6,211)
WRITE (6,212)
WRITE (6,213)J,NNS,NNF
NT=NNF-NNS+1
CALL COVAR(SUMX,SUMXX,SUMY,SUMYY,SUMXY,IREPS,NT)
CALL WRITE
50 CONTINUE
WRITE (6,211)
WRITE (6,212)
DO 51 I=1,LSUB
51 WRITE (6,213)I,NS(I),NF(I)
LTABLE=9999
IF(LCOMB)4020,60,52
60 CALL COVAR(SSX,SSXX,SSY,SSYY,SSXY,NSREP,LSUB)
CALL WRITE
GO TO 888
52 JOBCOM=0
DO 88 IJK =1,LCOMB
READ (5,108)IA1,NOCOM,(KOM(I),I=1,21)
JOBCOM=JOBCOM+1
IF(IA1.EQ.IF1)GO TO 525
523 WRITE (6,107)IJK
GO TO 88
525 IF(NOCOM-21)56,56,2344
2344 IF(NOCOM-99)2345,888,2345
2345 WRITE (6,2346)JOBCOM
2346 FORMAT( 52HMORE THAN 21 SUBSETS WERE CALLED FOR COMBINATION NO.I4)
GO TO 88
56 DO 54 J=1,NOCOM
SY(J)=0.0
SYY(J)=0.0
SX(J)=0.0
SXX(J)=0.0
SXY(J)=0.0
54 NREP(J)=0.0
JJ=0
DO 200 K=1,NOCOM
I=KOM(K)
IF(I-LSUB)53,53,950
950 WRITE (6,951)JOBCOM
951 FORMAT(115HTHE NUMBER OF A REQUESTED SUBSET IS GREATER THAN ANY EX
1ISTING SUBSET. THIS SUBSET WAS REQUESTED IN COMBINATION NO.I4)
GO TO 88
53 JJ=JJ+1
SY(JJ)=SY(JJ)+SSY(I)
SYY(JJ)=SYY(JJ)+SSYY(I)
SX(JJ)=SX(JJ)+SSX(I)
SXX(JJ)=SXX(JJ)+SSXX(I)
SXY(JJ)=SXY(JJ)+SSXY(I)
NREP(JJ)=NREP(JJ)+NSREP(I)
200 CONTINUE
CALL COVAR(SX,SXX,SY,SYY,SXY,NREP,NOCOM)
WRITE (6,416)
WRITE (6,214)JOBCOM
WRITE (6,215)(KOM(I),I=1,NOCOM)
CALL WRITE
88 CONTINUE
GO TO 888
END
SUBROUTINE COVAR(SUMX,SUMXX,SUMY,SUMYY,SUMXY,IREPS,NT)
C COVAR SUBROUTINE FOR BMD01R 9/4/63
COMMON Y , X , FADJ , LTABLE
COMMON NNF , NNS , DIXON , XX , XY , YY
COMMON TXX , TXY , TYY , NDF1 , NDF2 , NDF3
COMMON NDF4 , NDF5 , TXXT , TXYT , TYYT , EXXT
COMMON EXYT , EYYT , SE , SEMS , SXXT , SXYT
COMMON SYYT , STPE , ADJT , ADJTMS , BM , FXX
COMMON BW , FYY , BT , FB
DIMENSION SUMX(600),SUMY(600),SUMXX(600),SUMYY(600),SUMXY(600
1),IREPS(600), X(599),Y(599)
C
X(1)=0.0
XX=0.0
XY=0.0
YY=0.0
TXX=0.0
TXY=0.0
TYY=0.0
Y(1)=0.0
MREPS=0
DO 3 I=1,NT
X(1)=X(1)+SUMX(I)
XX=XX+SUMXX(I)
MREPS=MREPS+IREPS(I)
XY=XY+SUMXY(I)
Y(1)=Y(1)+SUMY(I)
YY=YY+SUMYY(I)
D1VR=IREPS(I)
IF(D1VR .EQ. 0.0) GO TO 3
TXX=TXX+SUMX(I)**2/D1VR
TXY=TXY+(SUMX(I)*SUMY(I))/D1VR
TYY=TYY+SUMY(I)**2/D1VR
3 CONTINUE
REPS=MREPS
SXXT=XX-X(1)**2/REPS
SXYT=XY-(X(1)*Y(1))/REPS
SYYT=YY-Y(1)**2/REPS
TXXT=TXX-X(1)**2/REPS
TXYT=TXY-X(1)*Y(1)/REPS
TYYT=TYY-Y(1)**2/REPS
EXXT=SXXT-TXXT
EXYT=SXYT-TXYT
EYYT=SYYT-TYYT
IF(LTABLE-9999)30,20,30
30 IF(LTABLE-99)20,25,20
25 NT=NNF-NNS+1
20 NDF1=NT-1
NDF2=MREPS-NT
NDF3=MREPS-1
NDF4=MREPS-NT-1
NDF5=MREPS-2
NDF6=1
DF1=NDF1
DF2=NDF2
DF3=NDF3
DF4=NDF4
DF5=NDF5
SE = EYYT
IF(EXXT .NE. 0)
1SE=EYYT-EXYT**2/EXXT
SEMS = 0.0
IF(DF4 .NE. 0.0)
1SEMS=SE/DF4
STPE = SYYT
IF(SXXT .NE. 0)
1STPE=SYYT-SXYT**2/SXXT
ADJT=STPE-SE
IF(NDF1)7,7,8
7 ADJTMS=0.0
GO TO 9
8 ADJTMS = 0.0
IF(DF1 .NE. 0.0)
1ADJTMS=ADJT/DF1
9 IF(NDF1)10,10,11
10 BM=0.0
GO TO 12
11 BM = 0.0
IF(TXXT .NE. 0.0)
1BM=TXYT/TXXT
12 BW = 0.0
IF(EXXT .NE. 0.0)
1BW=EXYT/EXXT
BT = 0.0
IF(SXXT .NE. 0.0)
1BT=SXYT/SXXT
IF(DF1)13,13,14
13 FXX=0.0
FYY=0.0
FADJ=0.0
GO TO 15
14 FXX = 0.0
IF(DF1.NE.0.0 .AND. EXXT.NE.0.0)
1FXX=(TXXT/DF1)/(EXXT/DF2)
FYY = 0.0
IF(DF1.NE.0.0 .AND. EYYT.NE.0.0)
1FYY=(TYYT/DF1)/(EYYT/DF2)
FADJ = 0.0
IF(SEMS .NE. 0.0)
1FADJ=ADJTMS/SEMS
15 FB = 0.0
IF(EXXT.NE.0.0 .AND. SEMS.NE.0.0)
1FB=(EXYT**2/EXXT)/SEMS
RETURN
END
SUBROUTINE TPWD(NT1,NT2)
CTPWD SUBROUTINE TPWD FOR BMD01R VERSION OF SEPT. 26, 1963
DIMENSION MAT(20)
C
IF(NT1)40,10,12
10 NT1=5
12 IF(NT1-NT2)14,19,14
14 IF(NT2-5)15,19,17
15 REWIND NT2
GO TO 19
17 REWIND NT2
19 IF(NT1-5)18,24,18
18 IF(NT1-6)22,40,22
40 WRITE (6,49)
STOP
49 FORMAT(25H ERROR ON TAPE ASSIGNMENT)
22 REWIND NT1
24 NT2=NT1
28 RETURN
END
SUBROUTINE TRANS1(X,NX,NTG,NCODE,CONS,IE)
CTRANS1 SUBROUTINE TRANS1 FOR BMD01R MAY 5, 1964
DIMENSION X(599),NCODE(8),CONS(8)
ASN(Q000FL)=ATAN(Q000FL/SQRT(1.0-Q000FL**2))
C
DO 100 I=1,NTG
JUMP=NCODE(I)
C=CONS(I)
DO 90 J=1,NX
T=X(J)
GO TO (1,2,3,4,5,6,7,8,9,10),JUMP
1 IF(T)99,11,12
11 X(J)=0.0
GO TO 90
12 X(J)=SQRT(T)
GO TO 90
2 IF(T) 99,13,14
13 X(J)=1.0
GO TO 90
14 X(J)=SQRT(T)+SQRT(T+1.0)
GO TO 90
3 IF(T) 99,99,15
15 X(J)=ALOG10(T)
GO TO 90
4 X(J)=EXP(T)
GO TO 90
5 IF(T) 99,16,17
16 X(J)=0.0
GO TO 90
17 IF(T-1.0) 18,19,99
18 X(J)=ASN(SQRT(T))
GO TO 90
19 X(J)=3.14159/2.0
GO TO 90
6 A=NX
A=T/A
IF(A) 99,20,21
20 X(J)=ASN(SQRT(1.0/FLOAT(NX+1)))
GO TO 90
21 IF(A-1.0) 22,22,99
22 A=NX+1
B=T/A
X(J)=ASN(SQRT(B))+ASN(SQRT(B+1.0/A))
GO TO 90
7 IF(T)24,99,24
24 X(J)=1.0/T
GO TO 90
8 X(J)=T+C
GO TO 90
9 X(J)=T*C
GO TO 90
10 IF(T) 99,25,26
25 X(J)=0.0
GO TO 90
26 X(J)=T**C
90 CONTINUE
100 CONTINUE
GO TO 200
99 WRITE (6,4010)J,IE,JUMP
GO TO 90
4010 FORMAT(18H THE VALUE OF CASE I4,19H IN TREATMENT GROUP I4,58H VIOL
1ATED THE RESTRICTIONS FOR TRANSGENERATION OF THE TYPE I3,1H./51H T
2HE PROGRAM CONTINUED LEAVING THE VALUE UNCHANGED.)
200 RETURN
END
SUBROUTINE VFCHCK(NVF)
CVFCHCK SUBROUTINE TO CHECK FOR PROPER NUMBER OF VARIABLE FORMAT CRDS
C
IF(NVF)10,10,20
20 IF(NVF-10)50,50,10
10 WRITE (6,4000)
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
1IED, ASSUMED TO BE 1.)
NVF=1
50 RETURN
END
SUBROUTINE WRITE
C WRITE SUBROUTINE FOR BMD01R 9/4/63
COMMON Y , X , FADJ , LTABLE
COMMON NNF , NNS , DIXON , XX , XY , YY
COMMON TXX , TXY , TYY , NDF1 , NDF2 , NDF3
COMMON NDF4 , NDF5 , TXXT , TXYT , TYYT , EXXT
COMMON EXYT , EYYT , SE , SEMS , SXXT , SXYT
COMMON SYYT , STPE , ADJT , ADJTMS , BM , FXX
COMMON BW , FYY , BT , FB
DIMENSION X(599),Y(599)
C
301 FORMAT(1H039X28HANALYSIS OF COVARIANCE TABLE//)
302 FORMAT(12H0 SOURCE OF4X7HDEGREES13X28HSUMS OF SQUARES AND PRODUCT
1S10X1H*9X26HDEVIATION ABOUT REGRESSION)
303 FORMAT(23H VARIATION FREEDOM51X1H*)
304 FORMAT(1H 32X2HXX14X2HXY14X2HYY7X1H*9X4HY*Y*6X2HDF8X12HMEAN SQ Y*Y
1*//)
305 FORMAT(11H MEANS6XI5,F17.4,2F16.4,4H *//)
306 FORMAT(11H WITHIN6XI5,F17.4,2F16.4,4H *F16.4,I5,F19.4//)
307 FORMAT(11H TOTAL6XI5,F17.4,2F16.4,4H *F16.4,I5//)
308 FORMAT(1H 12X62HDIFFERENCE FOR TESTING AMONG ADJUSTED MEANS......
1...........*F16.4,I5,F19.4////)
309 FORMAT(28H REGRESSION COEFFICIENTS46X12HF RATIOS FOR16X2HDF//)
310 FORMAT(14H MEANSF15.4,45X7HX MEANSF17.4,2H (I3,1H,I5,1H)/)
311 FORMAT(14H WITHINF15.4,45X7HY MEANSF17.4,2H (I3,1H,I5,1H)/)
312 FORMAT(14H TOTALF15.4,45X8HY* MEANSF16.4,2H (I3,1H,I5,1H)/
1)
313 FORMAT(1H 73X9HB WITHIN F15.4,2H (I3,1H,I5,1H))
DF5=NDF5
DF3=NDF3
NDF6=1
WRITE (6,301)
WRITE (6,302)
WRITE (6,303)
WRITE (6,304)
WRITE (6,305)NDF1,TXXT,TXYT,TYYT
WRITE (6,306)NDF2,EXXT,EXYT,EYYT,SE,NDF4,SEMS
WRITE (6,307)NDF3,SXXT,SXYT,SYYT,STPE,NDF5
WRITE (6,308)ADJT,NDF1,ADJTMS
WRITE (6,309)
WRITE (6,310)BM,FXX,NDF1,NDF2
WRITE (6,311)BW,FYY,NDF1,NDF2
WRITE (6,312)BT,FADJ,NDF1,NDF4
WRITE (6,313)FB,NDF6,NDF4
IF(DIXON) 401,400,401
401 WRITE (6,402)
402 FORMAT(/21H BIOASSAY INFORMATION//)
XLAM = 0.0
IF(DF5.NE.0.0 .AND. BT.NE.0.0)
1XLAM=SQRT(STPE/DF5)/BT
WRITE (6,403)XLAM
403 FORMAT(8H LAMBDA=F12.4/)
IF(SXXT.EQ.0.0) GO TO 1001
IF((SXYT**2/SXXT-(STPE/DF5)*DIXON**2).EQ.0.0) GO TO 1002
SL12=1./(2.*DF5+.5)+(STPE/DF5)/(SXYT**2/SXXT-(STPE/DF5)*DIXON**2)
GO TO 2000
1001 SL12=1.0/(2.0*DF5+0.5)-(1.0/DIXON**2)
GO TO 2000
1002 SL12 = 0.0
IF((2.0*DF5+0.5) .NE. 0.0)
1SL12=1.0/(2.0*DF5+0.5)
2000 IF(SL12.LE.0.0) GO TO 450
SL12=XLAM*SQRT(SL12)
WRITE (6,404)SL12
404 FORMAT('0CONFIDENCE LIMIT FOR LAMBDA=',F12.4/)
410 SM23=2.0*XLAM/SQRT(DF3+1.0)
IF(SL12.LE.0.0) GO TO 400
NNX=4.0*(XLAM+DIXON*SL12)**2/SM23**2+0.5
400 RETURN
450 WRITE(6, 460)
460 FORMAT('0THE APPROXIMATION USED IN THIS PROGRAM IS INAPPROPRIATE F
1OR THIS PROBLEM.'/'0THE CONFIDENCE LEVEL FOR LAMBDA IS NOT COMPUTE
2D.')
GO TO 410
END