Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmd05r.for
There is 1 other file named bmd05r.for in the archive. Click here to see a list.
C POLYNOMIAL REGRESSION JUNE 10, 1966
C THIS IS A SIFTED VERSION OF BMD05R ORIGINALLY WRITTEN IN
C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
DIMENSION SE(10),XX(500),YY(500),FMT(180),XXB(10),KODE(8),CONST(8)
DIMENSION X(400,6),SD(11,11),AV(11),B(10),RR(11,11)
DOUBLE PRECISION A1,A2,A3,SUM,X,SD,AV,B,RR,ONO,GSSDR,SSDR,SSAR,SU,
1A
DATA A1,A2,A3/6HPROBLM,6HFINISH,6HSPECTG/
C
918 FORMAT(45H1BMD05R - POLYNOMIAL REGRESSION - VERSION OF ,
118HJUNE 10, 1966 ,/
241H HEALTH SCIENCES COMPUTING FACILITY, UCLA//)
C
NTAPE=5
CALL USAGEB('BMD05R')
5 READ (5,900)SUM,PROB,N,NP,IPLOT,NVG,NACC,MTAPE,KVR
NPP=0
IF(SUM.EQ.A1)GO TO 4
IF(SUM.EQ.A2)GO TO 15
8 WRITE (6,919)
15 IF(6-NTAPE)16,4031,4031
16 REWIND NTAPE
GO TO 4031
4 CALL TPWD(MTAPE,NTAPE)
IF(NACC)748,748,749
748 NACC=13
749 TOL=10.0**(-NACC)
19 WRITE (6,918)
WRITE (6,902)PROB
WRITE (6,903)N
WRITE(6,9181) KVR
IF(NP*(NP-11))200,401,401
200 IF(N-501)204,402,402
204 IF(N-NP)403,403,205
205 IERROR=0
NT=NP+1
ONO=0.0D0
ON=N
ONO=ON
DO 721 I=1,500
X(I,NT)=0.0D0
721 X(I,1)=0.0D0
DO 20 I=1,NT
AV(I)=0.0D0
DO 20 J=1,NT
20 SD(I,J)=0.0D0
IF(KVR.GT.0.AND.KVR.LE.10)GO TO 21
KVR=1
WRITE (6,4000)
21 KVR=KVR*18
READ (5,916)(FMT(I),I=1,KVR)
WRITE(6,211) (FMT(I), I=1,KVR)
211 FORMAT(' THE VARIABLE FORMAT IS ',18A4,/(1H0,24X,18A4))
IF(-NVG)22,23,23
22 WRITE (6,901)
READ(5,920)SUM,M,(KODE(I),CONST(I),I=1,M)
IF(SUM.EQ.A3)GO TO 220
210 WRITE (6,921)
IERROR=-501
GO TO 23
220 DO 230 I=1,M
WRITE (6,922)KODE(I),CONST(I)
IF(KODE(I)*(KODE(I)-11))230,225,225
225 WRITE (6,923)KODE(I)
230 CONTINUE
23 READ (NTAPE,FMT)(X(I,1),X(I,NT),I=1,N)
IF(NTAPE-5)10,11,10
10 REWIND NTAPE
11 IF(NVG)25,25,24
24 CALL TRNGEN(N,NT,M,X,KODE,CONST)
IF(IERROR) 5, 25, 25
25 SUMY=0.0
SUMYY=0.0
SUMXY=0.0
SUMX=0.0
SUMXX=0.0
DO 27 I=1,N
YY(I)=X(I,NT)
XX(I)=X(I,1)
SUMX=SUMX+XX(I)
SUMY=SUMY+YY(I)
SUMXY=SUMXY+XX(I)*YY(I)
SUMXX=SUMXX+XX(I)*XX(I)
27 SUMYY=SUMYY+YY(I)*YY(I)
R=(SUMXY-SUMX*SUMY/ON)/SQRT((SUMXX-SUMX*SUMX/ON)*(SUMYY-SUMY*SUMY/
1ON))
RRR=R*R
IF(NP-1)1010,1010,1212
1212 DO 30 J=2,NP
JJ=J-1
DO 30 I=1,N
30 X(I,J)=X(I,JJ)*X(I,1)
1010 DO 36 J=1,NT
DO 35 I=1,N
35 AV(J)=AV(J)+X(I,J)
36 AV(J)=AV(J)/ONO
DO 40 J=1,NT
DO 40 I=1,N
40 X(I,J)=X(I,J)-AV(J)
DO 46 J=1,NT
DO 46 K=J,NT
DO 45 I=1,N
45 SD(J,K)=SD(J,K)+(X(I,J)*X(I,K))
46 SD(K,J)=SD(J,K)
B(1)=SD(1,NT)/SD(1,1)
A=AV(NT)-B(1)*AV(1)
SE(1)=DSQRT((SUMYY-A*SUMY-B(1)*SUMXY)/(SD(1,1)*(ON-2.0)))
WRITE (6,904)
WRITE (6,905)
WRITE (6,906)AV(1)
WRITE (6,907)AV(NT)
WRITE (6,908)A
WRITE (6,909)B(1)
WRITE (6,910)SE(1)
WRITE (6,911)R
SSAR=RRR*SD(NT,NT)
SSDR=(1.0-RRR)*SD(NT,NT)
N1=1
CALL REPORT (SSAR,SSDR,N1,N,NPP)
GSSDR=SSDR
IF(NP-1)50,50,1414
50 XXB(1)=B(1)
AXX=SNGL(A)
NXX=1
GO TO 1313
1414 DO 73 IJ=2,NP
N1=IJ
DO 52 J=1,IJ
B(J)=SD(J,NT)
DO 52 I=1,IJ
52 RR(I,J)=SD(I,J)
CALL MATINV(RR,IJ,B,D,T)
IF(T-TOL)66,191,191
66 WRITE (6,917)
SSDR=GSSDR
IF(NPP) 67, 67, 75
67 NPP=1
GO TO 75
191 DO 192 I=1,IJ
RR(I,I)=-RR(I,I)
192 B(I)=-B(I)
SSAR=0.0D0
DO 60 I=1,IJ
60 SSAR=SSAR+B(I)*SD(I,NT)
SSDR=SD(NT,NT)-SSAR
DF=IJ+1
VARI=SNGL(SSDR)/(ON-DF)
DO 62 J=1,IJ
SUMY=VARI*RR(J,J)
62 SE(J)=SQRT(SUMY)
SU=0.0D0
DO 64 I=1,IJ
64 SU=SU+B(I)*AV(I)
69 GSSDR=SSDR
A=AV(NT)-SU
WRITE (6,912)N1
WRITE (6,908)A
WRITE (6,913)
WRITE (6,915)(B(I), I=1,N1)
DO 70 I=1,N1
70 XXB(I)=B(I)
AXX=SNGL(A)
NXX=N1
WRITE (6,914)
WRITE (6,915)(SE(I), I=1,N1)
WRITE (6,904)
73 CALL REPORT (SSAR,SSDR,N1,N,NPP)
75 N1=105
WRITE (6,904)
CALL REPORT (SSAR,SSDR,N1,N,NPP)
1313 IF(IPLOT-1)5,600,5
600 CALL GRAPH (XX,YY,AXX,XXB,N,NXX)
GO TO 5
401 PRINT 4001
GO TO 4031
402 PRINT 4002
GO TO 4031
403 PRINT 4003
4031 STOP
900 FORMAT(A6,A2,I3,4I2,49X,2I2)
901 FORMAT(21H0TRANSGENERATION CARD//16H CODE CONSTANT//)
902 FORMAT(15H0PROBLEM NO. A2)
903 FORMAT(12H SAMPLE SIZEI5)
904 FORMAT(1H0)
905 FORMAT(38H0REGRESSION - ONE INDEPENDENT VARIABLE)
906 FORMAT(23H0XMEAN..... F12.5)
907 FORMAT(23H YMEAN..... F12.5)
908 FORMAT(23H0INTERCEPT (A VALUE)...F12.5)
909 FORMAT(23H REG. COEFFICIENT......F12.5)
910 FORMAT(25H STD. ERROR OF REG. COEF.F10.5)
911 FORMAT(23H0CORRELATION COEF......F12.5///)
912 FORMAT(32H0POLYNOMIAL REGRESSION OF DEGREEI3)
913 FORMAT(24H0REGRESSION COEFFICIENTS)
914 FORMAT(25H0STD. ERROR OF REG. COEF.)
915 FORMAT(10(1X,F12.5))
916 FORMAT(20A4)
917 FORMAT(71H0COMPUTATIONAL ACCURACY NOT SUFFICIENT FOR POLYNOMIALS O
1F HIGHER DEGREE)
9181 FORMAT(' NO. VF CDS.',I5)
919 FORMAT(22H0ERROR ON PROBLEM CARD)
920 FORMAT(A6,I1,8(I2,F6.0))
921 FORMAT(96H0SPECTG CARD MISPUNCHED OR OUT OF ORDER. PROGRAM WILL PR
1OCEED TO NEXT JOB, IF ANY, OR TERMINATE.)
922 FORMAT(1H I3,F12.5)
923 FORMAT(29H0ILLEGAL TRANSGENERATION CODEI3,67H SPECIFIED. PROGRAM W
1ILL CONTINUE LEAVING OUT THIS TRANSGENERATION.)
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
1IED, ASSUMED TO BE 1.)
4001 FORMAT(' DEGREE OF POLYNOMIAL INCORRECTLY SPECIFIED. IT CANNOT EX
XCEED 10')
4002 FORMAT(' NUMBER OF CASES TOO LARGE. IT CANNOT EXCEED 500.')
4003 FORMAT(' DEGREE OF POLYNOMIAL IS TOO HIGH.')
END
SUBROUTINE FORM2(SYMB,XY)
C SUBROUTINE FORM2 FOR PLOTR (IBM 360) JUNE 21, 1966
DIMENSION TEST(18)
INTEGER XY,SYMB,BLANK,TEST
DATA BLANK/2H /
DATA TEST/'B ','B ','4 ','5 ','6 ','7 ','8 ','9 ','A ','B ','C ',
1'D ','E ','F ','G ','H ','I ','/ '/
IF(XY.EQ.BLANK)GO TO 50
DO 30 I=1,17
IF(XY.NE.TEST(I))GO TO 30
C PUT IN NEXT SYMBOL OF ARRAY FOR MULTIPLE POINTS
XY=TEST(I+1)
GO TO 100
30 CONTINUE
IF(XY.EQ.TEST(18))GO TO 100
C IF OTHER THAN CHARACTERS IN ARRAY TEST PUT IN CHARACTER 2.
XY=TEST(1)
GO TO 100
C IF BLANK, PUT IN SYMBOL
50 XY=SYMB
100 RETURN
END
SUBROUTINE GRAPH(XX,YY,A,B,N,NPP)
C SUBROUTINE GRAPH FOR BMD05R JUNE 10, 1966
DIMENSION XX(500),YY(500),B(10),PP(500),FMT(104)
DIMENSION SYM(15),YYY(15)
DATA O,P/1HO,1HP/
DO 10 I=1,N
SUM=0.0
AA=1.0
DO 5 J=1,NPP
AA=AA*XX(I)
5 SUM=SUM+B(J)*AA
PP(I)=SUM+A
10 CONTINUE
WRITE (6,100)
WRITE (6,101)
DO 15 I=1,N
DIFF=YY(I)-PP(I)
15 WRITE (6,102)I,XX(I),YY(I),PP(I),DIFF
SS=10.0**10
FF=-10.0**9
XS=SS
XL=FF
DO 300 I=1,N
XS=AMIN1(XS,XX(I))
XL=AMAX1(XL,XX(I))
SS=AMIN1(SS,YY(I),PP(I))
300 FF=AMAX1(FF,YY(I),PP(I))
WRITE (6,4001)
WRITE (6,4000)
C MXY AND LYNN ARE DUMMY VARIABLES FOR SUBR SCALE.
CALL SCALE(SS,FF,100.0,MXY,S1,F1,LYNN)
302 SYM(1)=O
SYM(2)=P
X=XX(1)
YYY(1)=YY(1)
YYY(2)=PP(1)
NC=2
CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
DO 50 I=2,N
X=XX(I)
YYY(1)=YY(I)
YYY(2)=PP(I)
50 CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
NC=0
CALL PLOTR(X,XS,XL,YYY,SYM,SS,FF,NC,1)
WRITE (6,218)S1,F1
100 FORMAT(1H15X,18HTABLE OF RESIDUALS)
101 FORMAT(/6H0 NO. 4X,7HX VALUE4X,7HY VALUE4X,11HY PREDICTED4X,8HRESI
1DUAL//)
102 FORMAT(I5,F13.4,F11.4,2F13.4)
218 FORMAT(26H GRAPH SCALE EXTENDS FROM F10.4,4H TO F10.4)
4000 FORMAT(1H1,46X,38HPLOT OF OBSERVED AND PREDICTED VALUES./)
4001 FORMAT(12H0GRAPH CODES//13H O=OBSERVED Y/14H P=PREDICTED Y/21H B=O
1BSERVED=PREDICTED//)
RETURN
END
SUBROUTINE MATINV(A,N,B,D,T)
C SUBROUTINE MATINV FOR BMD05R JUNE 10, 1966
C THE FOLLOWING STATEMENT(S) HAVE BEEN MANUFACTURED BY THE TRANSLATOR---
C
DOUBLE PRECISION A , B , U , Y , V
DIMENSION A(11,11),B(10),U(10)
DIMENSION IND(10),C(10)
T=1.0
DO 50 I=1,N
IND(I)=0
50 C(I)=A(I,I)
D=1.0
DO 8 L=1,N
X=0.0
DO 4 I=1,N
IF(IND(I))4,2,4
2 IF(X-A(I,I)/C(I))3,4,4
3 K=I
X=A(I,I)/C(I)
4 CONTINUE
T=AMIN1(T,X)
11 IND(K)=1
DO 5 I=K,N
U(I)=A(I,K)
5 A(I,K)=0.D0
Y=U(K)
D=D*SNGL(Y)
DO 6 I=1,K
U(I)=A(K,I)
6 A(K,I)=0.D0
V=B(K)
B(K)=0.D0
U(K)=1.D0
DO 8 I=1,N
DO 7 J=1,I
7 A(I,J)=A(I,J)-U(I)*U(J)/Y
8 B(I)=B(I)-U(I)*V/Y
RETURN
END
SUBROUTINE PLOTR(X,ZMIN,ZMAX,Y,SYM,WMIN,WMAX,NC,NP)
C
C SUBROUTINE PLOTR (IBM 360) AUGUST 13, 1966
C
C 'PLOTR' IS A UTILITY SUBPROGRAM FOR THE BMD... PROGRAMS WHICH
C PLOTS EITHER SINGLE-LINE OR WHOLE-PAGE GRAPHS AND SETS UP
C APPROPRIATE SCALING. THE CALLING PARAMETERS ARE AS FOLLOWS -
C
C X - THE VALUE OF THE INDEPENDENT VARIABLE
C ZMIN - THE MINIMUM VALUE OF X FOR THIS PLOT
C ZMAX - THE MAXIMUM VALUE OF X FOR THIS PLOT
C Y - THE ARRAY CONTAINING THE VALUES OF UP TO 15 DEPENDENT VAR.'S
C SYM - THE ARRAY CONTAINING THE SYMBOLS TO BE PLOTTED
C WMIN - THE MINIMUM VALUE OF ALL Y'S FOR THIS PLOT
C WMAX - THE MAXIMUM VALUE OF ALL Y'S FOR THIS PLOT
C NC - THE NUMBER OF DEPENDENT VARIABLES
C NC=-1 CLOSES A SINGLE-LINE PLOT
C NC= 0 PRINTS AND CLOSES A WHOLE-PAGE PLOT
C NP - THE CONTROL VARIABLE
C NP=-1 PRINTS A SINGLE LINE
C NP=0 OR NP=1 SETS UP A WHOLE-PAGE PLOT
C
C THE PLOTTING ROUTINE MUST BE CALLED ONCE FOR EACH VALUE OF THE
C INDEPENDENT VARIABLE THAT IS TO BE PLOTTED NO MATTER WHETHER IN
C THE SINGLE-LINE OR WHOLE-PAGE MODE
C
DIMENSION Y(15),CLAB(12),GF(10),FMT(12),XY(51,101),SYM(15)
INTEGER XY,BLANKS
DATA TC,TP,BLANKS/1H.,1H+,1H /
DATA GF/ 4H 1X,,4H 2X,,4H 3X,,4H 4X,,4H 5X,,4H 6X,,
14H 7X,,4H 8X,,4H 9X,,4H 10X/
DATA FMT/'(17X',' ','5(F1','2.3,','8X)/','7X, ',' ','4(F1','2.3,',
1'8X),','F12.','3) '/
C
100 FORMAT(1H 6X5(F12.3,8X),F12.3/17X,5(F12.3,8X))
101 FORMAT(1H F12.3,1X,103A1,F12.3)
102 FORMAT (1H 13X,103A1)
1000 FORMAT(1H 14X,101A1)
1001 FORMAT(15X,20(5H+....),1H+)
C
DATA NCC/2/
C 'NCC' ON THE INITIAL ENTRY TO PLOTR IS NON-ZERO BECAUSE OF THE DATA
C STATEMENT ABOVE.
C
C 'NCC' IS 0 WHILE A PLOT IS BEING MADE. IT IS 1 OR 2 AT OTHER TIMES
C
IF(NCC) 50,48,50
C
C THE VARIABLE 'KL' CONTROLS THE FUNCTIONING OF THE OPENING AND
C CLOSING SECTIONS OF PLOTR. KL=0 INDICATES OPENING OF THE GRAPH,
C KL=1 INDICATES CLOSING.
C
50 KL=0
CALL SCALE(WMIN,WMAX,100.0,JY,YMIN,YMAX,YIJ)
YR=YMAX-YMIN
230 J=JY
IF(J*(J-10))204,201,201
C
C THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN FIXED FORMAT
C UNDER CONTROL OF KL
C
201 IF(KL)220,220,231
C
231 WRITE (6,1001)
IF(KL)250,250,220
C
220 CLAB(1)= YMIN
DO 222 I=2,11
222 CLAB(I)=CLAB(I-1)+YIJ
WRITE (6,100)(CLAB(I),I=1,11,2),(CLAB(J),J=2,10,2)
IF(KL)231,231,14
C
C THE FOLLOWING SECTION OPENS OR CLOSES A PLOT IN A VARIABLE
C FORMAT UNDER CONTROL OF KL AND JY FROM 'SCALE'
C
204 IF(J-5)205,221,207
207 J=J-5
205 JYT=5-J
221 CONTINUE
226 FMT(2)=GF(JY)
IF (KL) 225,225,227
C
225 FMT(7)=GF(JY)
TT=JY
TT=TT*YIJ/10.0
CLAB(1)= YMIN+TT
DO 223 I=2,10
223 CLAB(I)=CLAB(I-1) +YIJ
WRITE (6,FMT)(CLAB(I),I=2,10,2),(CLAB(I),I=1,9 ,2)
IF(KL)227,227,14
C
227 IF(JY-5)208,209,208
208 WRITE(6,1000)(TC,I=1,J ),(TP,(TC,I=1,4),K=1,19),TP,(TC,I=1,JYT)
IF (KL) 250,250,225
C
209 WRITE (6,1001)
IF (KL) 250,250,225
C
250 CONTINUE
NCC=0
IC=0
IF(NP)80,11,11
C
C THIS SECTION PREPARES FOR A FULL PAGE PLOT BY FILLING IN XY WITH
C BLANKS AND SETTING UP SCALING FOR THE INDEPENDENT VARIABLE 'X'
C
11 DO 1 I=1,51
DO 1 J=1,101
1 XY(I,J)=BLANKS
CALL SCALE (ZMIN,ZMAX,50.0,JX,XMIN,XMAX,XIJ)
XR=XMAX-XMIN
GO TO 48
C
C
C ENTRY TO PLOTS CAN BE USED ONLY AFTER THE CALLING PARAMETERS
C HAVE BEEN TRANSFERRED BY A CALL ON PLOTR. THE CALL ON PLOTS
C IS IDENTICAL WITH ENTRY TO PLOTR BUT IT ALLOWS THE PROGRAMMER TO
C CALL THE PLOTTING ROUTINE WITHOUT HAVING TO INCLUDE THE PARAMETERS
C
48 IF(NC)52,13,49
49 IF(NP)80,10,10
C THE FOLLOWING SECTION SETS UP A FULL PAGE BUT DOES NO PRINTING.
C THIS SECTION IS REACHED BY SPECIFYING NC POSITIVE AND NP POSITIVE.
C
10 DO 9 N=1,NC
SYMB=SYM(N)
XDIFFR=XMAX-X
IF(XDIFFR)105,106,106
105 XDIFFR=0.0
106 YDIFFR=YMAX-Y(N)
IF(YDIFFR)107,108,108
107 YDIFFR=0.0
108 L=51.0-(50.0*XDIFFR)/XR+.5
K=101.0-(100.0*YDIFFR)/YR+.5
CALL FORM2(SYMB,XY(L,K))
9 CONTINUE
GO TO 15
C
C THE FOLLOWING SECTION CONSTRUCTS AND PLOTS ONE LINE OF A MULTILINE
C GRAPH. LOCATION ALONG THE AXIS OF THE PAPER IS PRINTED AT EVERY
C STEP. THIS SECTION IS ACCESSIBLE BY SPECIFYING NC POSITIVE AND
C NP NEGATIVE.
C
80 DO 86 I=1,101
86 XY(1,I)=BLANKS
L=1
DO 95 N=1,NC
SYMB=SYM(N)
YDIFFR=YMAX-Y(N)
IF(YDIFFR)860,865,865
860 YDIFFR=0.0
865 K=101.0-(100.0*YDIFFR)/YR+.5
95 CALL FORM2(SYMB,XY(L,K))
IF(MOD(IC,5))97,96,97
96 W=TP
GO TO 98
97 W=TC
98 WRITE (6,101)X,W,(XY(1,N),N=1,101),W,X
IC=IC+1
GO TO 15
C
C THIS SECTION PLOTS OUT THE PREVIOUSLY PREPARED WHOLE PAGE GRAPH.
C IT PRINTS LOCATION ALONG THE PAPER'S AXIS EVERY FIFTH STEP. THIS
C SECTION IS ACCESSED BY SPECIFYING NC=0.
C
13 M=6-JX
LL=50+M
T=JX
IF(5-JX)131,131,135
131 T=0.0
135 RLAB=XMAX-(T*XIJ)/5.0
W=TC
K=52
DO 31 L=M,LL
K=K-1
I=MOD(L,5)
IF(I-1)2,3,2
3 W=TP
WRITE (6,101)RLAB,W,(XY(K,N),N=1,101),W,RLAB
RLAB=RLAB-XIJ
W=TC
GO TO 31
2 WRITE (6,102)W,(XY(K,N),N=1,101),W
31 CONTINUE
C
52 KL=1
GO TO 230
C
14 NCC=1
15 RETURN
END
SUBROUTINE REPORT (SSAR,SSDR,N1,N,NPP)
C SUBROUTINE REPORT FOR BMD05R JUNE 10, 1966
DIMENSION SS(10)
IF(N1-100) 10, 30, 30
10 NN=N-1
N2=NN-N1
D1=N1
D2=N2
SSARM=SSAR/D1
SSDRM=SSDR/D2
F=SSARM/SSDRM
TOTAL=SSAR+SSDR
IF(N1-1) 15, 15, 20
15 WRITE (6,900)
GO TO 21
20 WRITE (6,901)N1
21 WRITE (6,902)
WRITE (6,903)
WRITE (6,904)N1,SSAR,SSARM,F
WRITE (6,905)N2,SSDR,SSDRM
WRITE (6,906)NN, TOTAL
NPP=NPP+1
SS(NPP)=SSAR
GO TO 70
30 N1=1
WRITE (6,907)NPP
WRITE (6,921)
WRITE (6,922)
DO 60 I=1,NPP
GO TO (31,32,33,34,35,36,37,38,39,40), I
31 WRITE (6,911)N1,SS(1),SS(1)
GO TO 60
32 SSAR=SS(2)-SS(1)
IF(SSAR) 41, 42, 42
41 SSAR=0.0
42 WRITE (6,912)N1,SSAR,SSAR
GO TO 60
33 SSAR=SS(3)-SS(2)
IF(SSAR) 43, 44, 44
43 SSAR=0.0
44 WRITE (6,913)N1,SSAR,SSAR
GO TO 60
34 SSAR=SS(4)-SS(3)
IF(SSAR) 45, 46, 46
45 SSAR=0.0
46 WRITE (6,914)N1,SSAR,SSAR
GO TO 60
35 SSAR=SS(5)-SS(4)
IF(SSAR) 47, 48, 48
47 SSAR=0.0
48 WRITE (6,915)N1,SSAR,SSAR
GO TO 60
36 SSAR=SS(6)-SS(5)
IF(SSAR) 49, 50, 50
49 SSAR=0.0
50 WRITE (6,916)N1,SSAR,SSAR
GO TO 60
37 SSAR=SS(7)-SS(6)
IF(SSAR) 52, 53, 53
52 SSAR=0.0
53 WRITE (6,917)N1,SSAR,SSAR
GO TO 60
38 SSAR=SS(8)-SS(7)
IF(SSAR) 54, 55, 55
54 SSAR=0.0
55 WRITE (6,918)N1,SSAR,SSAR
GO TO 60
39 SSAR=SS(9)-SS(8)
IF(SSAR) 56, 57, 57
56 SSAR=0.0
57 WRITE (6,919)N1,SSAR,SSAR
GO TO 60
40 SSAR=SS(10)-SS(9)
IF(SSAR) 58, 59, 59
58 SSAR=0.0
59 WRITE (6,920)N1,SSAR,SSAR
60 CONTINUE
WRITE (6,905)N2,SSDR,SSDRM
WRITE (6,906)NN,TOTAL
900 FORMAT(65H0 ANALYSIS OF VARIANCE FOR SIMPLE LINEAR R
1EGRESSION)
901 FORMAT(41H0 ANALYSIS OF VARIANCE FOR I4, 19H DEGRE
1E POLYNOMIAL)
902 FORMAT(85H0 SOURCE OF VARIATION DEGREE OF S
1UM OF MEAN F)
903 FORMAT(87H FREEDOM S
1QUARES SQUARE VALUE)
904 FORMAT(37H0DUE TO REGRESSION............. I6,F19.5,2F14.5)
905 FORMAT(37H DEVIATION ABOUT REGRESSION.... I6,F19.5,F14.5)
906 FORMAT(37H TOTAL.... I6,F19.5///)
907 FORMAT(45H0 FINAL ANALYSIS OF VARIANCE FOR I4, 19H D
1EGREE POLYNOMIAL)
911 FORMAT(37H0LINEAR TERM................... I6,F19.5,F14.5)
912 FORMAT(37H QUADRATIC TERM................ I6,F19.5,F14.5)
913 FORMAT(37H CUBIC TERM.................... I6,F19.5,F14.5)
914 FORMAT(37H QUARTIC TERM.................. I6,F19.5,F14.5)
915 FORMAT(37H QUINTIC TERM.................. I6,F19.5,F14.5)
916 FORMAT(37H SEXTIC TERM................... I6,F19.5,F14.5)
917 FORMAT(37H 7TH DEGREE TERM............... I6,F19.5,F14.5)
918 FORMAT(37H 8TH DEGREE TERM............... I6,F19.5,F14.5)
919 FORMAT(37H 9TH DEGREE TERM............... I6,F19.5,F14.5)
920 FORMAT(37H 10TH DEGREE TERM.............. I6,F19.5,F14.5)
921 FORMAT(74H0 SOURCE OF VARIATION DEGREE OF S
1UM OF MEAN)
922 FORMAT(75H FREEDOM S
1QUARES SQUARE)
70 RETURN
END
C SUBROUTINE 'SCALE' CALCULATES THE SCALING FOR 'PLOTR'
C
C SUBROUTINE SCALE FOR PLOTR JUNE 21, 1966
C
SUBROUTINE SCALE(YMIN,YMAX,YINT,JY,TYMIN,TYMAX,YIJ)
DIMENSION C(10)
DATA C /1.0,1.5,2.0,3.0,4.0,5.0,7.5,10.0,15.0,20.0/
DATA TEST / 0.76293945E-05/
50 YR=YMAX-YMIN
TT=YR/YINT
J = ALOG10(TT)+TEST
E=10.0**J
TT=TT/E
I=0
IF(TT-1.0+TEST)205,201,201
205 TT=TT*10.0
E=E/10.0
201 I=I+1
IF(9-I)1,2,2
1 E=E*10.0
I=1
2 IF(TT-C(I))233,202,201
233 YIJ=C(I)*E
GO TO 203
202 Y=YMIN/C(I)
J=Y
T=J
IF(0.0001-ABS(T-Y))204,233,233
204 YIJ=C(I+1)*E
203 X=((YMAX+YMIN)/YIJ-YINT )/2.0+.00001
K=X
IF(K)235,240,240
235 Y=K
IF(X-Y)236,240,236
236 K=K-1
240 TYMIN=K
TYMIN=YIJ*TYMIN
TYMAX=TYMIN+YINT*YIJ
IF (YMAX-TYMAX-TEST)11,11,201
11 YIJJ=C(I)*E
XT=((YMAX+YMIN)/YIJJ-YINT)/2.0+.00001
KT=XT
IF (KT) 1235,1240,1240
1235 YT=KT
IF (XT.NE.YT) KT=KT-1
1240 TYMINT=KT
TYMINT=YIJJ*TYMINT
TYMAXT=TYMINT+YINT*YIJJ
IF (YMAX-TYMAXT.GT.TEST) GO TO 10
TYMIN=TYMINT
TYMAX=TYMAXT
YIJ=YIJJ
K=KT
10 TT=YINT/10.0
JY=TT+.000001
YIJ=YINT*(YIJ/10.0)
J=TYMIN/ YIJ
IF (K)242,241,241
242 J=J-1
241 J=J*JY+JY-K
JY=J
RETURN
END
SUBROUTINE TPWD(NT1,NT2)
C SUBROUTINE TPWD FOR BMD05R JUNE 10, 1966
IF(NT1)40,10,12
10 NT1=5
12 IF(NT1-NT2)14,19,14
14 IF(NT2.EQ.5)GO TO 18
17 REWIND NT2
19 IF(NT1-5)18,24,18
18 IF(NT1-6)22,40,22
22 REWIND NT1
24 NT2=NT1
28 RETURN
40 WRITE (6,49)
49 FORMAT(25H ERROR ON TAPE ASSIGNMENT)
STOP
END
SUBROUTINE TRNGEN(N,NT,M,X,LCODE,CONST)
C SUBROUTINE TRNGEN FOR BMD05R JUNE 10, 1966
DIMENSION LCODE(8),CONST(8)
DOUBLE PRECISION X(400,6)
ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
4 DO 210 I=1,M
FM=CONST(I)
JUMP=LCODE(I)
DO 85 J=1,N
D=X(J,NT)
GO TO(10,20,30,40,50,60,70,80,90,110),JUMP
10 IF(D)99,11,12
11 D2=0.0
GO TO 8
12 D2=SQRT(D)
GO TO 8
20 IF(D)99,21,22
21 D2=1.0
GO TO 8
22 D2=SQRT(D)+SQRT(D+1.0)
GO TO 8
30 IF(D)99,99,31
31 D2=ALOG10(D)
GO TO 8
40 D2=EXP(D)
GO TO 8
50 IF(-D)52,11,99
52 IF(D-1.0)53,54,99
54 D2=3.14159265/2.0
GO TO 8
53 D2=ASN(SQRT(D))
GO TO 8
60 SAMP=N
A=D/(SAMP+1.0)
B=A+1.0/(SAMP+1.0)
IF(A)99,61,62
61 IF(-B)64,11,99
64 D2=ASN(SQRT(B))
GO TO 8
62 IF(B)99,65,66
65 D2=ASN(SQRT(A))
GO TO 8
66 D2=ASN(SQRT(A))+ASN(SQRT(B))
GO TO 8
70 IF(D)71,99,71
71 D2=1.0/D
GO TO 8
80 D2=D+FM
GO TO 8
90 D2=D*FM
GO TO 8
99 WRITE (6,105)J,JUMP
GO TO 85
110 IF(D)111,11,111
111 D2=D**FM
8 X(J,NT)=D2
85 CONTINUE
210 CONTINUE
1000 RETURN
105 FORMAT(18H0THE VALUE OF CASEI4,83H OF THE DEPENDENT VARIABLE VIOLA
1TED THE RESTRICTION FOR TRANSGENERATION OF THE TYPE, I3,1H./52H TH
2E PROGRAM CONTINUED LEAVING THIS VALUE UNCHANGED.)
END