Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/bmd/bmd03s.for
There is 1 other file named bmd03s.for in the archive. Click here to see a list.
C BIOLOGICAL ASSAY - PROBIT ANALYSIS JUNE 14, 1966
C THIS IS A SIFTED VERSION OF BMD03S ORIGINALLY WRITTEN IN
C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE
C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION.
C EXTERNAL DIFF
DOUBLE PRECISION AAA, CODE,Q000HL,Q001HL
DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
1, SMP(1001), P(1001), X(1001), SMY(1001)
2, A(3), DC(9), F(3)
3, CD(10), DV(10), D(10)
4, TEMP(1000), RX(1000)
5,DC0(9),CD0(9)
COMMON CODE , R , SN , ICODE , DOSE , CONST
COMMON CG , M , K , IPR , SMP , P
COMMON SMY , X , A , DC , F
ER = 0.0
CALL USAGEB('BMD03S')
C READ IN DATA
998 CONTINUE
CALL READIN
SAVC = A(3)
DO 2099 JK=1,K
A(3) = SAVC
IK = 1+8*(JK-1)
C TRANSFORM DOSE.
WRITE (6,9998)JK,K
CALL TRNGEN( DOSE, CONST(IK), M, ICODE(IK), 1, IERROR, X)
IF(IERROR)2098, 70,2098
70 CONTINUE
WRITE (6,9993)(X(I),I=1,M)
73 CONTINUE
IF(A(2))2000,2000,2001
C GET LEAST SQUARES FIT
2000 CONTINUE
I1 = 1
DO 2210 I = 1,M
IF(1.0-SMP(I)) 2210,2210,2206
2206 IF( SMP(I) -0.0) 2210,2210,2208
2208 TEMP(I1) = SMY(I)
RX(I1) = X(I)
M2 = I1
I1 = I1 +1
2210 CONTINUE
CALL LINREG ( RX, TEMP, M2, A(2), A(1))
2001 CONTINUE
CALL MNVAR ( X, M, XAV, XVAR)
WRITE (6,9988)A(1),A(2),A(3)
IF (A(2))54,54,2201
54 WRITE (6,46)
GO TO 2098
2201 CONTINUE
DO 7002 I = 1,M
7002 X(I) = X(I) - XAV
A(1) = A(1) + A(2) *XAV
C SOLVE FOR MAXIMUM LIKLIHOOD BY NEWTONS METHOD
EPS = .00001
NPR = IPR
C FIRST TIME GET BEST APPROX TO ALPHA, BETA
CG1 = CG
CG = 1.0
DO 2376 NK = 1,2
IF (CG) 204, 2041, 2041
204 NV = 3
GO TO 2042
2041 NV = 2
2042 CONTINUE
NIT = 50
CALL NEWT
1(DC,F,NV,NIT,TEMP,EPS,ER,NPR,A)
CG = CG1
2444 IF(ER)2098,2043,2044
2044 CONTINUE
2043 CONTINUE
IF (CG1)2376,2377,2377
2376 CONTINUE
2377 CONTINUE
CG = -1.0
CALL DIFF
CG = CG1
HAMU = (5.0- A(1))/A(2)
ALPHA = A(1)-A(2)*XAV
ALPHA1 = A(1)
A(1) = A(1) - A(2) *XAV
HAMU0 = (5.0-A(1))/A(2)
HASI = 1.0/A(2)
C COMPUTR COVARIANCE MATRICES
IF (CG1) 305, 207, 207
C TWO DIMENSION CASE
207 CONTINUE
N1 = 2
NN = 4
DATA Q000HL/6H. /
AAA=Q000HL
IDF = M-2
DC(4) = DC(5)
DC(3) = DC(2)
DV(1) = -DC(1)/HASI**2
DV(4) = -(HAMU*(HAMU*DC(1) + 2.0*DC(2)) + DC(5))/HASI**4 + 2.0*(HA
1MU*F(1)+F(2))/HASI**3
DV(2) = (HAMU*DC(1) +DC(2))/HASI**3
DV(3) = DV(2)
DC0(1) =-DC(1)
DC0(2) =-DC(1)*XAV-DC(2)
DC0(4) = -DC(4)-2.0*DC(2)*XAV -DC(1)*XAV**2
DC0(3) = DC0(2)
GO TO 306
C THREE DIMENSION CASE
305 CONTINUE
N1 = 3
NN = 9
DATA Q001HL/6H, C. /
AAA=Q001HL
IDF = M-3
IF(SN(M+1)) 3051,3051,3052
3052 CONTINUE
IDF = IDF + 1
3051 CONTINUE
DV(1) = -DC(1)/HASI**2
DV(5) = -(HAMU*(HAMU*DC(1)+2.0*DC(2))+DC(5))/HASI**4 + 2.0*(HAMU*F
1(1)+F(2))/HASI**3
DV(9) = -DC(9)
DV(2) = (HAMU*DC(1)+DC(2))/HASI**3
DV(4) = DV(2)
DV(3) = -DC(3)/HASI
DV(7) = DV(3)
DV(6) = (HAMU*DC(3) +DC(6))/HASI**2
DV(8) = DV(6)
DC0(1) = -DC(1)
DC0(2) = - DC(1)*XAV-DC(4)
DC0(3) = -DC(3)
DC0(5) = -DC(1)*XAV**2 -2.0*DC(4)*XAV -DC(5)
DC0(6) = -DC(3)*XAV-DC(6)
DC0(9) = -DC(9)
DC0(4) = DC0(2)
DC0(7) = DC0(3)
DC0(8) = DC0(6)
306 CONTINUE
DO 3061 I=1,NN
3061 DC(I) = -DC(I)
CALL INVERT( DC, CD, N1, DET,TEMP(20) )
CALL INVERT( DV, D, N1, DET,TEMP(20))
CALL INVERT (DC0,CD0,N1,DET,TEMP(20))
C COMPUTE CHI SQUARE
M1 = M
IF (CG1)308,3081,3081
308 CONTINUE
M1 = M + 1
3081 CONTINUE
CHI = 0.0
DO 307 I = 1,M
RI = R(I)
PP = A(3) + P(I)*(1.0-A(3))
CHI = CHI + (RI-SN(I)*PP)**2 / (PP-PP**2) /SN(I)
307 CONTINUE
C PRINT RESULTS
WRITE (6,9885)
WRITE (6,9980)ALPHA1, F(1), ALPHA, A(2), F(2),A(2), A(3), F(3), A(
13), HAMU0, HASI
WRITE (6,9981)AAA ,AAA
DO 3071 I= 1,N1
I1 = 1 + N1*(I-1)
I3 = N1*I
GO TO (3071,3170,3171),N1
3170 WRITE (6,9883)(CD(IJ),IJ=I1,I3),(CD0(IJ),IJ=I1,I3)
GO TO 3071
3171 WRITE (6,9884)(CD(IJ),IJ=I1,I3),(CD0(IJ),IJ=I1,I3)
3071 CONTINUE
WRITE (6,9982)AAA
DO 3072 I= 1,N1
I1 = 1 + N1*(I-1)
I3 = N1*I
3072 WRITE (6,9983)( D(IJ),IJ=I1,I3)
WRITE (6,9984)CHI,IDF
2098 CONTINUE
WRITE (6,9989)JK,K
A(2) = 0.0
2099 CONTINUE
WRITE (6,38)CODE
GO TO 998
C NORMAL OUTPUT FORMAT
9883 FORMAT(3X, 2F12.4, 35X, 1H*, 2F12.4)
9884 FORMAT(3X,3F12.4,23X, 1H*, 3F12.4 )
9885 FORMAT (67H0THE FINAL VALUES OF THE PARAMETERS AND THE PARTIAL DER
1IVATIVES ARE /
2 5X,40HFOR PROBIT LINE (ALPHA1+BETA*(X-XBAR)) 17X,1H*,5X,
3 35HFOR PROBIT LINE (ALPHA+BETA*(X)) )
9980 FORMAT(
15X, 7HALPHA1= F8.3, 10X,6HDL/DA=3X,E10.4 ,13X,1H*,5X,7HALPHA =F8.3
2/5X,7HBETA =F8.3, 10X,6HDL/DB= 3X,E10.4, 13X, 1H*, 5X,7HBETA =F8
3.3/ 5X, 1HC 5X,1H=, F8.3, 10X, 6HDL/DC= 3X, E10.4,13X,1H*, 5X,
41HC,5X, 1H= F8.3 /
5 120H0THE ESTIMATES MU-HAT AND SIGMA-HAT OF THE MEAN AND STANDARD
6 DEVIATION OF THE TRANSFORMED DOSE TOLERANCE ARE AS FOLLOWS. /
7 5X, 10H MU-HAT= 3X, F9.4 /
8 5X, 10HSIGMA-HAT= 3X, F9.4 )
9981 FORMAT(39H0THE COVARIANCE MATRICES ARE AS FOLLOWS /
1, 5X, 16HFOR ALPHA1, BETA A6,35X, 1H*,5X,15HFOR ALPHA, BETA A6)
9982 FORMAT ( 5X, 13HFOR MU, SIGMA A6 )
9983 FORMAT (3X, 9F12.4)
9984 FORMAT (12H0CHI-SQUARE=3X, F12.4, 6H, WITHI5, 20H DEGREES OF FREED
1OM. )
9985 FORMAT (50H THE MEAN AND VARIANCE OF THE TRANSFORMED DOSE ARE /
1 7H MEAN= F8.4, 7H VAR= F8.4)
9988 FORMAT (54H0THE INITIAL VALUES AT START OF ITERATION ARE, ALPHA=
1 F9.3, 3X,5HBETA= F9.3, 3X, 2HC= F9.3)
9993 FORMAT (3X,9F12.2 )
9998 FORMAT (///
1 45H0TRANSFORMED DOSE FOR TRANSFORMATION CARD NO I3,3H OF
2 ,I3)
9989 FORMAT (46H END FOR OUTPUT USING TRANSFORMATION CARD NO I3, 3H OF
1, I3 //)
38 FORMAT(1H05X,26HEND OF OUTPUT FOR PROBLEM A6 )
C ERROR OUTPUT FORMAT
46 FORMAT (48H0THE DATA GIVEN IS UNWORKABLE, BETA IS NEGATIVE. )
END
C SUBROUTINE DERIV2 FOR BMD03S JUNE 14, 1966
C DERIVATIVES WITH RESPECT TO A AND B.
SUBROUTINE DERIV2 ( X, YY, Z, P, R, SN, M, C, F, D)
DIMENSION X(1),Z(1),YY(1),P(1),SN(1)
1, F(1), D(1), R(1)
CQ = 1.0-C
DO 10 I=1,4
D(I) = 0.0
10 CONTINUE
F(1) = 0.0
F(2) = 0.0
12 CONTINUE
PCC = 0.0
DO 100I=1,M
PP = C+P(I)*CQ
Q = 1.0- PP
PQ = PP*Q
XLP = (R(I) - PP*SN(I)) /PQ
XLPP = -R(I)/PP**2 -(SN(I)-R(I))/Q**2
PA = CQ*Z(I)
PB = PA*X(I)
PC = 1.0-P(I)
PAA = -PA*YY(I)
PAB = -PB*YY(I)
PBB = PAB*X(I)
PAC = -Z(I)
PBC = PAC*X(I)
F(1) = F(1) + XLP*PA
F(2) = F(2) + XLP*PB
D(1) = D(1) + XLPP*PA*PA + XLP*PAA
D(2) = D(2) +XLPP*PA*PB + XLP*PAB
D(4) = D(4) + XLPP*PB*PB + XLP*PBB
100 CONTINUE
D(3) = D(2)
RETURN
END
C SUBROUTINE DERIV3 FOR BMD03S JUNE 14, 1966
C DERIVATIVES WITH RESPECT TO A, B, AND C.
SUBROUTINE DERIV3 ( X, YY, Z, P, R, SN, M, C, F, D)
DIMENSION X(1),Z(1),YY(1),P(1),SN(1)
1, F(1), D(1), R(1)
CQ = 1.0-C
DO 10 I=1,3
F(I) = 0.0
D(I) = 0.0
10 CONTINUE
DO 12 I=4,9
D(I) = 0.0
12 CONTINUE
PCC = 0.0
M1 = M
IF (C) 14,13,14
14 M1 = M+1
13 CONTINUE
DO 100I=1,M1
PP = C+P(I)*CQ
Q = 1.0- PP
PQ = PP*Q
XLP = (R(I) - PP*SN(I)) /PQ
XLPP = -R(I)/PP**2 -(SN(I)-R(I))/Q**2
PA = CQ*Z(I)
PB = PA*X(I)
PC = 1.0-P(I)
PAA = -PA*YY(I)
PAB = -PB*YY(I)
PBB = PAB*X(I)
PAC = -Z(I)
PBC = PAC*X(I)
F(1) = F(1) + XLP*PA
F(2) = F(2) + XLP*PB
F(3) = F(3) + XLP*PC
D(1) = D(1) + XLPP*PA*PA + XLP*PAA
D(2) = D(2) +XLPP*PA*PB + XLP*PAB
D(5) = D(5) + XLPP*PB*PB + XLP*PBB
D(3) = D(3) + XLPP*PA*PC + XLP*PAC
D(6) = D(6) + XLPP*PB*PC + XLP*PBC
D(9) = D(9) + XLPP*PC*PC + XLP*PCC
100 CONTINUE
D(4) = D(2)
D(8) = D(6)
D(7) = D(3)
900 RETURN
END
C SUBROUTINE DIFF FOR BMD03S JUNE 14, 1966
C DIFFERENTIAL COEFICIENTS AND FUNCTION AUXILIARY SUBROUTINE FOR NEWT
SUBROUTINE DIFF
DOUBLE PRECISION CODE
DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
1, SMP(1001), P(1001), Z(1001), X(1001), YY(1001),SMY(1001),Y(1001)
2, A( 3), DC(9), F(3)
COMMON CODE , R , SN , ICODE , DOSE , CONST
COMMON CG , M , K , IPR , SMP , P
COMMON SMY , X , A , DC , F
C PREDICT PROBIT YY AND PROBABILITY P.
DO 200 I=1,M
Y(I) = A(1) + A(2)*X(I)
YY(I)=Y(I) -5.0
CALL PHI2(YY(I), P(I), Z(I))
200 CONTINUE
M1 = M+1
P(M1) = 0.0
YY(M1)=0.0
X(M1) = 0.0
Z(M1)=0.0
2049 IF (CG)300,2050,2050
C C CONSTANT
2050 CONTINUE
CALL DERIV2 ( X, YY, Z, P, R, SN, M, A(3), F, DC)
RETURN
C OPTIMIZE C
300 CONTINUE
CALL DERIV3 ( X, YY, Z, P, R, SN, M, A(3), F, DC)
900 RETURN
END
C SUBROUTINE INVERT FOR BMD03S JUNE 14, 1966
SUBROUTINE INVERT ( G, R, N, DET, ER)
DIMENSION G(1), R(1)
ER = 0.0
1 K = N
GO TO (10,20,30),K
10 R(1)= 1.0/G(1)
RETURN
20 CONTINUE
DET = G(1)*G(4) - G(3)*G(2)
IF (DET)80,99,80
99 ER = 1.0
199 RETURN
80 CONTINUE
R(1) = G(4)/DET
R(2) = -G(2)/DET
R(3) = -G(3)/DET
R(4) = G(1)/DET
RETURN
30 CONTINUE
R(1) = (G(5)*G(9) - G(6)*G(8))
R(4) = -(G(4)*G(9) - G(6)*G(7))
R(7) = (G(4)*G(8) - G(5)*G(7))
R(2) = -(G(2)*G(9) - G(3)*G(8))
R(5) = (G(1)*G(9) - G(3)*G(7))
R(8) = -(G(1)*G(8) - G(2)*G(7))
R(3) = (G(2)*G(6) - G(3)*G(5))
R(6) = -(G(1)*G(6) - G(3)*G(4))
R(9) = (G(1)*G(5) - G(2)*G(4))
DET = R(1)*G(1) + R(4)*G(2) + R(7)*G(3)
35 IF(DET)84,99,84
84 CONTINUE
DO 32 I = 1, 9
32 R(I) = R(I)/DET
88 RETURN
9999 FORMAT (1H E12.4 )
END
C SUBROUTINE LINREG FOR BMD03S JUNE 14, 1966
C REGRESS Y AS A FUNCTION OF X. Y=AX + B.
SUBROUTINE LINREG ( X, Y, M, A, B)
DIMENSION X(1), Y(1)
XM =M
SX = 0.0
SY = 0.0
SXX= 0.0
SXY= 0.0
DO 10 I=1,M
SXY = SXY +X(I) * Y(I)
SX = SX+ X(I)
SY= SY+Y(I)
SXX = SXX+ X(I)**2
10 CONTINUE
A = (XM*SXY -SX*SY)/(XM*SXX-SX**2)
B = (SY -A*SX)/XM
1000 RETURN
END
C SUBROUTINE MNVAR FOR BMD03S JUNE 14, 1966
C COMPUTES MEAN AND VARIANCE
SUBROUTINE MNVAR ( X, N, XB, XV)
DIMENSION X(1)
XN1 = N-1
XN = N
XV=0.0
XB=0.0
DO 10 I=1,N
XB = XB+X(I)
10 CONTINUE
XB = XB/XN
DO 20 I=1,N
XV = XV + (X(I) - XB)**2
20 CONTINUE
XV = XV/XN1
RETURN
END
C SUBROUTINE NEWT FOR BMD03S JUNE 14, 1966
SUBROUTINE NEWT
1(D,F,N,NI,TEMP,EPS,ER,NPR,A)
DIMENSION D(1),F(1),TEMP(1)
1, A(1)
ER=0.0
NM = N+1
NF = N
NN = N**2
IPR = NI+1
IF (NPR)8,8,5
5 IPR = 1
8 CONTINUE
DO 200 K =1,NI
CALL DIFF
CALL INVERT( D, TEMP(NM), NF,DET, TEMP(20) )
9 DEL = 0.0
IF(DET)10,99,10
10 CONTINUE
DO 35 I = 1, N
TEMP(I) = 0.0
IJ = I
DO 30 J = 1, N
IJ = IJ+N
TEMP(I) = TEMP(I) - F(J)*TEMP(IJ)
30 CONTINUE
DEL = DEL + TEMP(I)**2
A(I) = A(I) + TEMP(I)
35 CONTINUE
IF(K-IPR)18,15,18
15 CONTINUE
IPR = IPR + NPR
WRITE (6,800)K
WRITE (6,805)(F(J), J=1,N)
WRITE (6,801)
WRITE (6,805)(A(J), J=1,N)
WRITE (6,802)
DO 16J=1,N
JL = N*(J-1)+1
JU = JL + N-1
WRITE (6,805)(D(IJ),IJ=JL,JU)
16 CONTINUE
18 CONTINUE
IF (DEL - EPS) 80, 80, 190
190 CONTINUE
200 CONTINUE
WRITE (6,806)NI
ER = 1.0
GO TO 80
99 CONTINUE
WRITE (6,807)
ER = -1.0
NI = K
80 RETURN
800 FORMAT (16H0ITERATION COUNT I5, 20H. VALUE OF FUNCTION. )
801 FORMAT ( 27H APPROXIMATION TO SOLUTION.)
805 FORMAT ( 3H , 7E16.6)
806 FORMAT ( 43H0CRITERION FOR CONVERGENCE IS NOT MET AFTER I5,
1 11H ITERATIONS )
802 FORMAT (34H DIFFERENTIAL COEFFICIENT MATRIX. )
807 FORMAT (78H0MATRIX OF DIFFERENTIAL COEFICIENTS IS SINGULAR, PROBLE
1M CANNOT BE CONTINUED. )
END
C FUNCTION PHINV FOR BMD03S JUNE 14, 1966
C USES HASTINGS APPROXIMATIONS
FUNCTION PHINV (P)
IF(P-.5)10,100,12
100 PHINV = 0.0
RETURN
12 CONTINUE
Q = 1.0-P
T = -1.0
GO TO 13
10 Q = P
T = 1.0
13 CONTINUE
IF(Q)121,121,122
121 CONTINUE
PHINV = -T*999.0
RETURN
122 CONTINUE
E=SQRT(ALOG(1.0/Q**2))
X = -E+ ((.010328*E+.802853)*E +2.515517) /(((.001308*E+.189269)*E
1 + 1.432788)*E+1.0)
PHINV=T*X
RETURN
END
C SUBROUTINE PHI2 FOR BMD03S JUNE 14, 1966
C COMPUTES PHI(X) AND PHI-PRIME(X)
SUBROUTINE PHI2( X, PHI, PP)
SQ2 = .141421356E+01
Y = ABS(X)/SQ2
YY = X*X/2.0
PP = EXP(-YY)*.39894228
P1 = PP*2.0*SQ2
ET = 1.0/ (1.0+.3275911*Y)
PH = 1.0-((((.94064607*ET-.128782245E+01)*ET+1.25969513)*ET-.25212
18668)*ET+.225836846)*ET*P1
IF(X) 1, 4, 3
1 PH = -PH
3 PHI = .5 + PH/2.0
RETURN
4 PHI = 0.0
RETURN
END
C SUBROUTINE READIN FOR BMD03S JUNE 14, 1966
SUBROUTINE READIN
DOUBLE PRECISION PPP,AAA,TRG,AA,PRBLM,TODE,CODE,FMT
1,Q003HL,Q004HL,Q005HL,Q006HL,Q007HL,Q008HL
DIMENSION DOSE(1001), R(1001), SN(1001), ICODE(64), CONST(64)
1, SMP(1001), P(1001), X(1001), SMY(1001)
2, A(3), DC(9), F(3)
3, AA(3), FMT(120)
COMMON CODE , R , SN , ICODE , DOSE , CONST
COMMON CG , M , K , IBG , SMP , P
COMMON SMY , X , A , DC , F
C
11 FORMAT(53H1BMD03S BIOLOGICAL ASSAY, PROBIT ANALYSIS - REVISED ,2X
1,'FEBRUARY 11, 1969'/
242H HEALTH SCIENCES COMPUTING FACILITY, UCLA.//)
C
DATA PPP,AAA,TRG/6HPROBLM,6HFINISH,6HSPECTG/
DATA Q003HL/6HOPTIMI/
DATA Q004HL/6HZE C. /
DATA Q005HL/6H /
DATA Q006HL/6HC HELD/
DATA Q007HL/6H CONST/
DATA Q008HL/6HANT. /
1 NTAP=5
2 ER = 0.0
READ (5,10)PRBLM, CODE, IC, CG, M, R(1), SN(1), A(1), A(2), IBG, K
1, IPR , MTAP, IVF
IF(PRBLM.NE.PPP) GO TO 990
49 CONTINUE
IF(M*(M-1001))495,993,993
495 M1=M+1
R(M1) = R(1)
SN(M1) = AMAX1(SN(1),0.0)
50 CONTINUE
K=MIN0(K,8)
IF(K) 991,991,51
51 DO 501 I=1,K
I2 = 1+8*(I-1)
I3 = I2+7
READ(5,19)TODE,K1,(ICODE(I1),CONST(I1),I1=I2,I3)
IF(TODE.NE.TRG)GO TO 994
501 CONTINUE
IF(-K)502,887,887
502 IF(IVF.GT.0.AND.IVF.LE.10)GO TO 52
IVF=1
WRITE (6,4000)
52 IVF=IVF*12
READ (5,47)(FMT(I),I=1, IVF )
CALL TPWD(MTAP,NTAP)
53 READ (NTAP,FMT)(DOSE(I),I=1,M)
READ (NTAP,FMT)(R(I), I=1,M)
READ (NTAP,FMT)(SN(I), I =1,M)
C GET SAMPLE P ( = R/S )
C GET INVERSE NORMAL OF SAMPLE P.
DO 60 I=1,M
IF(SN(I)) 54,54,540
540 CONTINUE
SMP(I)=R(I)/SN(I)
IF( SMP(I)*(SMP(I)-1.0))58,58,54
58 CONTINUE
60 CONTINUE
C CHECK OPTION FOR C
IF(ABS(CG)-1.0) 261,28,28
261 IF(CG) 29,262,262
262 IF(IC)263,30,263
263 IF(CG) 29,28,29
C WE GUESS C
28 IF (SN(M1)) 281,281,282
281 CONTINUE
SN(M1) = 0.0
R(M1) = 0.0
SMALN = SN(1)
SMALP = SMP(1)
GO TO 2821
282 CONTINUE
SMALP = R(M1) /SN(M1)
SMALN = SN(M1)
2821 CONTINUE
DO 301 I=1,M
IF(SMP(I) - SMALP) 300,301,301
300 SMALN = SN(I)
SMALP = SMP(I)
301 CONTINUE
A(3) = SMALP
IF(A(3)) 54,302,303
302 A(3) = .5/SMALN
303 CONTINUE
GO TO 291
C USER GUESS C
29 A(3) = ABS(CG)
291 CONTINUE
AA(1)=Q003HL
AA(2)=Q004HL
AA(3)=Q005HL
CG = -1.0
GO TO 32
C HOLD C CONSTANT
30 AA(1)=Q006HL
AA(2)=Q007HL
AA(3)=Q008HL
A(3) = CG
32 CONTINUE
IF (SN(M1))611,611,610
611 SMP(M1) = 0.0
SMY(M1) = 5.0+ PHINV(0.0)
GO TO 612
610 CONTINUE
SMP(M1) = R(M1)/SN(M1)
SMY(M1) = 5.0+ PHINV(SMP(M1))
612 CONTINUE
QC = 1.0- A(3)
DO 6110 I = 1,M
SMP(I) = ( SMP(I) - A(3))/QC
SMY(I) = 5.0+ PHINV(SMP(I))
6110 CONTINUE
C PRINT INPUT DATA
210 CONTINUE
WRITE (6,11)
WRITE (6,9990)CODE
WRITE (6,9991)SN(M1), R(M1), AA
IF (IPR)25,21,25
21 CONTINUE
WRITE (6,9992)
WRITE (6,9993)(DOSE(I) ,I=1,M)
WRITE (6,9994)
WRITE (6,9993)(R(I),I=1,M)
WRITE (6,9996)
WRITE (6,9993)(SN(I), I=1,M)
IF(ER)41,25,41
41 WRITE (6,9985)
WRITE (6,38)CODE
GO TO 1
25 CONTINUE
WRITE (6,9981)
WRITE (6,9993)(SMY(I),I= 1,M)
WRITE (6,9982)
WRITE (6,9993)(SMP(I),I= 1,M)
RETURN
990 CONTINUE
IF(PRBLM.EQ.AAA)GO TO 999
GO TO 992
C DATA ERROR EXIT
54 CONTINUE
995 ER = 1.0
GO TO 210
C LAST PROBLEM COMPLETED.
999 CONTINUE
WRITE (6,9995)
887 IF(NTAP.EQ.5)GO TO 889
REWIND NTAP
889 STOP
C PROBLEM CARD ERROR
992 CONTINUE
WRITE (6,11)
WRITE (6,9997)
GO TO 887
991 CONTINUE
WRITE (6,11)
WRITE (6,9998)CODE,K
GO TO 887
993 WRITE (6,11)
WRITE (6,9999)
GO TO 887
994 WRITE (6,11)
WRITE (6,9989)KODE
K=-1001
GO TO 501
10 FORMAT( A6, A6, I1, F5.5, I4, 4F6.0, 18X, I2,2I1,2I2)
19 FORMAT ( A6, I1, 8( I2,F6.0))
47 FORMAT (12A6)
4000 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
1IED, ASSUMED TO BE 1.)
9989 FORMAT(48H0CONTROL CARD ERROR. SPECTG CARD EXPECTED BUT A A6,35H C
1ARD WAS FOUND. EXECUTION DELETED.)
9990 FORMAT (8H PROBLEM 1X,A6, 1H. )
9991 FORMAT (35H CONTROL GROUP DATA. SAMPLE SIZE N= F6.0, 13H, RESPONSE
1 R= , F6.0, 1H.
2 3X, 3A6)
9992 FORMAT (20H UNTRANSFORMED DOSE. )
9993 FORMAT (3X,9F12.2 )
9994 FORMAT (10H RESPONSE. )
38 FORMAT(1H05X,26HEND OF OUTPUT FOR PROBLEM A6 )
9981 FORMAT (14H SAMPLE PROBIT )
9982 FORMAT (19H SAMPLE PROBABILITY )
9985 FORMAT (33H0DATA INCORRECT,PROBLEM DELETED. )
9995 FORMAT (52H0NO NEW PROBLEMS, PROGRAM TERMINATED BY FINISH CARD. )
9996 FORMAT (13H SAMPLE SIZE. )
9997 FORMAT (66H0NON PROBLEM CARD ENCOUNTERED CHECK DATA DECK. EXECUTIO
1N DELETED. )
9998 FORMAT(75H THE NUMBER OF SPECIAL TRANSGENERATION CARDS IS ILLEGAL.
1 PROBLEM CARD CODE= 3X, A6, 3X, 2HT=, 3X, I3, 1H. )
9999 FORMAT(54H0ILLEGAL NUMBER OF DOSES SPECIFIED. EXECUTION DELETED.)
END
C SUBROUTINE TPWD FOR BMD03S JUNE 14, 1966
SUBROUTINE TPWD(NT1,NT2)
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)
STOP
49 FORMAT(25H ERROR ON TAPE ASSIGNMENT)
END
C SUBROUTINE TRNGEN FOR BMD03S JUNE 14, 1966
SUBROUTINE TRNGEN (DOSE,CONST,M,ICODE,K,IERROR,X)
DIMENSION DOSE(1000),CONST(10),ICODE(10),X(1000)
ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
IERROR = 0
DO 2 J2 = 1,M
2 X(J2) = DOSE(J2)
IF (ICODE(1)) 1000,1000,5
5 CONTINUE
DO 210 JK = 1,8
JUMP = ICODE(JK)
FM = CONST(JK)
IF(JUMP*(JUMP-11)) 200,1000,1000
200 DO 8 J=1,M
D = X(J)
GO TO (10,20,30,40,50,60,70,80,90,110), JUMP
C 1 SQRT(DOSE)
10 IF(D)99,11,12
11 TD = 0.0
GO TO 8
12 TD = SQRT(D)
GO TO 8
C 2 SQRT(DOSE) + SQRT(DOSE + 1)
20 IF( D) 99,21,22
21 TD = 1.0
GO TO 8
22 TD = SQRT(D) + SQRT(D+1.0)
GO TO 8
C 3 LOG10( DOSE)
30 IF(D) 99,99,31
31 TD =ALOG10(D)
GO TO 8
C 4 EXP(DOSE).
40 TD = EXP(D)
GO TO 8
C 5 ARCSINE(DOSE)
50 IF(-D)52,11,99
52 IF(D-1.0) 53,54,99
54 TD = .157079632E+01
GO TO 8
53 CONTINUE
TD = ASN(A)
GO TO 8
C 6 ARCSIN(SQRT(DOSE/(M+1))) + ARCSIN(SQRT(A+10(M+1)))
60 SAMP = M
A=D/(SAMP+1.0)
B=A+1.0/(SAMP+1.0)
IF (A)99,61,61
61 TD = ASN(SQRT(A))
62 IF (B) 99,65,65
65 TD = TD + ASN(SQRT(B))
GO TO 8
C 7 1/ DOSE
70 IF(D) 71,99,71
71 TD = 1.0/D
GO TO 8
C 8 DOSE + CONST
80 TD = D + FM
GO TO 8
C 9 DOSE X CONST
90 TD = D*FM
GO TO 8
C 10 DOSETO POWER OF CONST
110 TD = D**FM
8 X(J) = TD
210 CONTINUE
1000 RETURN
99 IERROR =-999
WRITE (6,9980)JUMP,DOSE(J),D
GO TO 1000
9980 FORMAT(43H0TRANSFORMATION CANNOT BE PERFORMED. CODE = I4,
1 6H DOSE= F12.4, 5H ARG= F12.4)
END