Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmod22.bas
There are 2 other files named cmod22.bas in the archive. Click here to see a list.
00012REM 10/4/76
00015FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00020REM*****************************************************************
00025REM CMOD22 CMOD22 CMOD22 CMOD22 CMOD22 CMOD22
00030REM******************************************************************
00035REM CONSENSUS PREPOSTERIOR--GEORGE WOODWORTH
00040REM*******************************************************************
00045REM SHORTEST CONSENSUS CREDIBILITY REGION - CALLING PROGRAM
00050REM INPUT M0,M1,M2,M3,M4,M5,M6,M7,M8,M9
00055REM OUTPUT El OR V1,V2
00075RESTORE#1
00076 INPUT# 1,I1,I2,I3
00080SCRATCH#1
00081 PRINT # 1,53,I2,I3
00085PRINT L$
00090PRINT " CONSENSUS PREPOSTERIOR"
00095PRINT " "
00100PRINT "IF YOU DON'T YET HAVE THE PARAMETERS OF THE PRIOR MARGINAL DISTRI-"
00105PRINT "BUTIONS ON THE MEAN AND STANDARD DEVIATION FOR BOTH INVESTIGATOR"
00110PRINT "AND ADVERSARY, USE THE TWO PARAMETER NORMAL MODEL OF THE CADA"
00115PRINT "MONITOR TO OBTAIN THESE. TO EXIT THIS MODULE TO OBTAIN THESE"
00120PRINT "PARAMETERS, TYPE A '0'. OTHERWISE, TYPE '1' AND CONTINUE. ";
00125REM PRIOR D OF F FOR STD DEV = MS-1
00130REM PRIOR SCALE PARAM = SQR(MS)*SIGMA-TILDE
00135REM PRIOR MEAN FOR POP MEAN = W-DOT
00140REM PRIOR SCALE PARAM = (MS/MM)*SIGMA-TILDE**2
00145REM SEE NOVICK AND JACKSON P 224.
00150GOSUB 9000
00155C0=O1
00160IF C0=0 THEN 770
00165PRINT " "
00170PRINT " "
00175GOSUB 3515
00180GOSUB 4825
00185PRINT L$
00190PRINT "INPUT THE PARAMETERS OF THE INVESTIGATOR'S PRIOR MARGINAL DISTRI-"
00195PRINT "BUTIONS ON THE MEAN AND THE STANDARD DEVIATION OF THE NORMAL"
00200PRINT "DISTRIBUTION."
00205PRINT "FROM THE PRIOR ON THE STANDARD DEVIATION:"
00210PRINT "INPUT MODAL ESTIMATE OF STANDARD DEVIATION.";
00215GOSUB 9000
00220IF O1>0 THEN 235
00225PRINT "STANDARD DEVIATION MUST BE POSITIVE. PLEASE RESPECIFY.";
00230GOTO 215
00235PRINT "INPUT DEGREES OF FREEDOM";
00240GOSUB 9000
00245N1=O1
00250IF N1 >= 6 THEN 265
00255PRINT "DEGREES OF FREEDOM MUST BE AT LEAST 6. PLEASE RESPECIFY.";
00260GOTO 240
00265T8=O1*SQR(N1+1)
00270PRINT " "
00275PRINT "FROM THE PRIOR ON THE MEAN:"
00280PRINT "INPUT POINT ESTIMATE OF THE MEAN";
00285GOSUB 9000
00290N5=O1
00295PRINT "INPUT STANDARD DEVIATION OF T DISTRIBUTION ON MEAN.";
00300GOSUB 9000
00305IF O1>0 THEN 320
00310PRINT "STANDARD DEVIATION MUST BE POSITIVE. PLEASE RESPECIFY.";
00315GOTO 300
00320N3=O1*O1*(N1-2)
00325PRINT L$
00330PRINT "INPUT THE PARAMETERS OF THE ADVERSARY'S PRIOR MARGINAL DISTRIBU-"
00335PRINT "TIONS ON THE MEAN AND THE STANDARD DEVIATION OF THE NORMAL"
00340PRINT "DISTRIBUTION."
00345PRINT "FROM THE PRIOR ON THE STANDARD DEVIATION:"
00350PRINT "INPUT MODAL ESTIMATE OF STANDARD DEVIATION."
00355GOSUB 9000
00360IF O1>0 THEN 375
00365PRINT "STANDARD DEVIATION MUST BE POSITIVE. PLEASE RESPECIFY.";
00370GOTO 355
00375PRINT "INPUT DEGREES OF FREEDOM";
00380GOSUB 9000
00385N2=O1
00390IF N2 >= 6 THEN 405
00395PRINT "DEGREES OF FREEDOM MUST BE AT LEAST 6. PLEASE RESPECIFY.";
00400GOTO 380
00405T9=O1*SQR(N2+1)
00410PRINT ""
00415PRINT "FROM THE PRIOR ON THE MEAN:"
00420PRINT "INPUT POINT ESTIMATE OF THE MEAN."
00425GOSUB 9000
00430N6=O1
00435PRINT "INPUT STANDARD DEVIATION OF T DISTRIBUTION ON MEAN."
00440GOSUB 9000
00445IF O1>0 THEN 460
00450PRINT "STANDARD DEVIATION MUST BE POSITIVE. PLEASE RESPECIFY.";
00455GOTO 440
00460N4=O1*O1*(N2-2)
00465M3=N1+1
00470M4=N2+1
00475M5=N5
00480M6=N6
00485M7=T8/SQR(M3)
00490M8=T9/SQR(M4)
00495M1=T8*T8/N3
00500M2=T9*T9/N4
00505PRINT L$
00510PRINT "INPUT P FOR THE P% SMALLEST CONSENSUS CREDIBILITY INTERVAL-SCCI"
00515GOSUB 9000
00520M0=O1
00525M0=M0/100
00530IF M0<.01 THEN 540
00535IF M0 <= .99 THEN 555
00540PRINT "MUST BE PERCENT BE IN THE INTERVAL 1 THROUGH 99."
00545PRINT "PLEASE RESPECIFY."
00550GOTO 515
00555PRINT "IF YOU WANT THE ENDPOINTS OF THE PRIOR S C C I "
00560PRINT "INPUT ZERO FOR N. IF YOU WANT THE EXPECTED LENGTH OF THE"
00565PRINT "POSTERIOR CONSENSUS HIGHEST DENSITY REGION, INPUT THE NUMBER OF"
00570PRINT "SAMPLE OBSERVATIONS FOR N."
00575PRINT "N=";
00580GOSUB 9000
00585M9=INT(O1)
00590PRINT ""
00595N1=M3-1
00600N2=M4-1
00605N3=M7*M7*M3/M1
00610N4=M8*M8*M4/M2
00615N5=M5
00620N6=M6
00625N8=SQR(N3/N1)
00630N9=SQR(N4/N2)
00635IF M9>2 THEN 705
00640IF M9=1 THEN 690
00645IF M9=2 THEN 690
00650IF M9<0 THEN 680
00655GOSUB 1615
00660GOSUB 2235
00665PRINT "LOWER ENDPOINT OF SCCI IS ";V1;" UPPER ENDPOINT IS ";V2
00670PRINT "LENGTH OF INTERVAL IS ";V2-V1
00675GOTO 750
00680PRINT "N MUST BE EQUAL TO ZERO OR GREATER THAN 2 AND LESS THAN 10000."
00685GOTO 555
00690PRINT "THIS PROGRAM CANNOT CALCULATE EXPECTED LENGTH FOR N=1 OR 2."
00695PRINT "USE N=0 OR N>2."
00700GOTO 555
00705IF M9>9999 THEN 680
00709B2=(M3-1)/2
00710A=(M9-1)/2
00715PRINT L$
00720PRINT "PLEASE BE PATIENT SINCE REACHING A CONSENSUS OFTEN TAKES A WHILE."
00725Z3=M9+M3-2
00730GOSUB 1055
00735GOSUB 4265
00740PRINT " "
00745PRINT "EXPECTED LENGTH OF THE REGION IS ";L4
00750PRINT "IF YOU WISH TO EVALUATE WITH ANOTHER SAMPLE SIZE TYPE '1' ELSE '0'"
00755GOSUB 9000
00760C0=O1
00765IF C0=1 THEN 570
00770CHAIN "CMOD2A"
00775REM END OF CALLING PROGRAM
01055REM ****************************************************************
01065REM SETUP ROUTINE FOR INTEGRATION
01075REM INPUT A B2 M0 M1 M2 M3 M9 N1 N2 Z3
01085REM OUTPUT G1 G2 G3 N1 N2 Y6 Y7 Y8 Y9
01095REM OUTPUT T U V Z ARRAYS
01105 DIM T(3),U(2),V(2),X(3),Y(2),Z(3)
01115Z(3)=.774597
01125Z(2)=0
01135Z(1)=-Z(3)
01145X(3)=5/9
01155X(2)=8/9
01165X(1)=X(3)
01175Y(2)=.57735
01185Y(1)=-Y(2)
01195V(2)=1
01205V(1)=1
01215B=B2-.5
01225G9=B
01235GOSUB 5305
01245G8=G0
01255G9=B2
01265GOSUB 5305
01275G8=G8-G0
01285G9=A+B2
01295GOSUB 5305
01305G8=G8+G0
01315G9=A+B
01325GOSUB 5305
01335G8=G8-G0
01345G3=EXP(G8)
01355FOR I4=1 TO 2
01365S0=(Y(I4)+1)/2
01375GOSUB 4945
01385GOSUB 1775
01395U(I4)=R0
01405NEXT I4
01415FOR I7=1 TO 3
01425Q2=(Z(I7)+1)/2
01435V0=Z3-1
01445GOSUB 3145
01455T(I7)=Y
01465NEXT I7
01475G9=(Z3+1)/2
01485GOSUB 5305
01495G8=G0
01505G9=(Z3-1)/2
01515GOSUB 5305
01525G8=G8+G0
01535G9=Z3/2
01545GOSUB 5305
01555G8=G8-2*G0
01565G3=EXP(G8)*SQR((Z3-1)/Z3)*G3
01575G1=M9*M1/(M9+M1)
01585G2=M9*M2/(M9+M2)
01595N1=M3+M9-1
01605N2=M4+M9-1
01615Q2=(1+M0)/2
01625V0=N1
01635GOSUB 3145
01645Y6=Y
01655V0=N2
01665GOSUB 3145
01675Y7=Y
01685Q2=M0
01695V0=N1
01705GOSUB 3145
01715Y8=Y
01725V0=N2
01735GOSUB 3145
01745Y9=Y
01755RETURN
01765REM END OF SETUP ROUTINE
01775REM *****************************************************************
01785REM BETA PERCENTILE ROUTINE
01795REM INPUT A B S0
01805REM OUTPUT R0
01815GOSUB 3515
01825GOSUB 4945
01835GOSUB 3515
01845Q=1-S0
01855GOSUB 3615
01865R3=V
01875R4=(R3*R3-3)/6
01885R5=1/(2*A-1)
01895R6=1/(2*B-1)
01905R7=2/(R5+R6)
01915R1=R3*SQR(R7+R4)/R7-(R6-R5)*(R4+5/6-2/(3*R7))
01925R1=A/(A+B*EXP(2*R1))
01935J2=R1
01945GOSUB 4495
01955S1=P
01965W6=S1-S0
01975IF ABS(W6)>.0001 THEN 2005
01985R0=R1
01995RETURN
02005R2=R1-W6*SQR((A*B)/(A+B+1))/(A+B)
02015IF R2 >= 0 THEN 2035
02025R2=R1/2
02035IF R2 <= 1 THEN 2055
02045R2=(1+R1)/2
02055FOR I3=1 TO 10
02065J2=R2
02075GOSUB 4495
02085S2=P
02095W6=S2-S0
02105R9=(S2-S1)/(R2-R1)
02115IF ABS(W6)>.0001 THEN 2145
02125R0=R2
02135RETURN
02145R1=R2
02155S1=S2
02165R2=R2-W6/R9
02175NEXT I3
02185PRINT "BETI DID NOT CONVERGE IN 10 ITERATIONS"
02195PRINT "A=",A,"B=",B,"DESIRED TAIL AREA=",S0
02205PRINT "RESULTING TAIL AREA AFTER 10 ITERATIONS IS ",S1
02215RETURN
02225REM END OF BETA PERCENTILE ROUTINE
02235REM *****************************************************************
02245REM CONSENSUS HIGHEST DENSITY REGION ROUTINE
02255REM FINDS ENDPOINTS OF REGION
02265REM INPUT N1 N2 N3 N4 N5 N6 N8 N9 Y6 Y7 Y8 Y9
02275REM OUTPUT V1 V2
02285N8=SQR(N3/N1)
02295N9=SQR(N4/N2)
02305W5=N8
02306IF N8<N9 THEN 2308
02307W5=N9
02308W5=W5*.001
02315X2=N5-Y6*N8
02325X3=N5+Y6*N8
02335X4=N6-Y7*N9
02345X5=N6+Y7*N9
02355X6=N5-Y8*N8
02365X7=N6-Y9*N9
02375IF X7>X6 THEN 2545
02385X0=X4
02395GOSUB 2865
02405D5=D4
02415IF D5>0 THEN 2455
02425V1=X4
02435V2=X5
02445RETURN
02455X8=X4
02465X9=X7-W5
02475X0=X9
02485GOSUB 2865
02495D7=D4
02505IF D7<0 THEN 2805
02515V1=X9
02525V2=N7+D7/2
02535RETURN
02545X0=X2
02555GOSUB 2865
02565D5=D4
02575IF D5<0 THEN 2615
02585V1=X2
02595V2=X3
02605RETURN
02615IF X6>X4 THEN 2715
02625X8=X2
02635X9=X6-W5
02645X0=X9
02655GOSUB 2865
02665D7=D4
02675IF D7>0 THEN 2805
02685V1=X9
02695V2=N7-D7/2
02705RETURN
02715X0=X4
02725GOSUB 2865
02735D7=D4
02745IF D7>0 THEN 2785
02755V1=X4
02765V2=X5
02775RETURN
02785X8=X2
02795X9=X4
02805K4=10
02815GOSUB 3715
02825V1=L3
02835V2=N7
02845RETURN
02855REM END OF CONSENSUS HIGHEST DENSITY REGION ROUTINE
02865REM *****************************************************************
02875REM D ROUTINE
02885REM DIFFERENCE BETWEEN PERCENTILE POINTS OF TWO
02895REM T DISTRIBUTIONS
02905REM INPUT M0 N1 N2 N5 N6 N8 N9 X0
02915REM OUTPUT D4
02925T2=(X0-N5)/N8
02935T3=(X0-N6)/N9
02945G=N1
02955J2=T2
02965GOSUB 5595
02975T4=M0+P
02985G=N2
02995J2=T3
03005GOSUB 5595
03015T5=M0+P
03025Q2=T4
03035V0=N1
03045GOSUB 3145
03055T6=N5+Y*N8
03065Q2=T5
03075V0=N2
03085GOSUB 3145
03095T7=N6+Y*N9
03105D4=T6-T7
03115N7=(T6+T7)/2
03125RETURN
03135REM END OF D ROUTINE
03145REM *********************************************************************
03155REM T PERCENTILE ROUTINE
03165REM INPUT Q2 V0
03175REM OUTPUT Y
03185Q=Q2
03195E1=.001
03205GOSUB 3615
03215Y1=V*SQR(V0/(V0-2))
03225G=V0
03235J2=Y1
03245N=V0
03255GOSUB 5845
03265GOSUB 5595
03275Q3=P
03285W4=Q3-Q2
03295IF ABS(W4)>E1 THEN 3325
03305Y=Y1
03315RETURN
03325Y2=Y1-3*W4
03335FOR I2=1 TO 10
03345J2=Y2
03355GOSUB 5595
03365Q4=P
03375W4=Q4-Q2
03385Z1=(Q4-Q3)/(Y2-Y1)
03395IF ABS(W4)>E1 THEN 3425
03405Y=Y2
03415RETURN
03425Y1=Y2
03435Q3=Q4
03445Y2=Y2-W4/Z1
03455NEXT I2
03465PRINT "JCTINV DID NOT CONVERGE IN 10 ITERATIONS"
03475PRINT "DEGREES OF FREEDOM=",V0,"DESIRED TAIL AREA=",Q2
03485PRINT "RESULTING TAIL AREA AFTER 10 ITERATIONS IS",Y
03495RETURN
03505REM END OF T PERCENTILE ROUTINE
03515REM *********************************************************************
03525REM NORMAL PERCENTILE ROUTINE - INITIALIZATION
03535REM INPUT Q
03545REM OUTPUT V
03555REM GOSUB 2900 REQUIRED ONCE PRIOR TO GOSUB 2925
03565K0=2.30753
03575K1=.27061
03585K2=.99229
03595K3=.04481
03605RETURN
03615REM NORMAL PERCENTILE ROUTINE - PERCENTILES
03625Q1=Q
03635S=-1
03645IF Q<.5 THEN 3675
03655Q1=1-Q
03665S=1
03675T1=SQR(LOG(1/(Q1*Q1)))
03685V=S*(T1-(K0+K1*T1)/(1+T1*(K2+K3*T1)))
03695RETURN
03705REM END OF NORMAL PERCENTILE ROUTINE
03715REM ********************************************************************
03725REM ROOTFINDER ROUTINE
03735REM INPUT D5 D7 E1 N8 N9 X8 X9
03745REM OUTPUT L3
03755E=N8
03756IF N8<N9 THEN 3758
03757E=N9
03758E=E*100*E1
03765L0=X8
03775L2=X9
03785L1=(L0*D7-L2*D5)/(D7-D5)
03795X0=L1
03805GOSUB 2865
03815D6=D4
03825FOR K5=1 TO K4
03835L3=D6*D7*L0/((D5-D6)*(D5-D7))
03845L3=L3+D5*D7*L1/((D6-D5)*(D6-D7))
03855L3=L3+D5*D6*L2/((D7-D5)*(D7-D6))
03865IF L3>L0 THEN 3895
03875L3=(L0+L1)/2
03885GOTO 3915
03895IF L3<L2 THEN 3915
03905L3=(L1+L2)/2
03915X0=L3
03925GOSUB 2865
03935D8=D4
03945IF ABS(D8)>E THEN 3965
03955RETURN
03965IF D5*D6<0 THEN 4085
03975L0=L1
03985D5=D6
03995IF L3>L1 THEN 4055
04005L1=(L0+L2)/2
04015X0=L1
04025GOSUB 2865
04035D6=D4
04045GOTO 4185
04055L1=L3
04065D6=D8
04075GOTO 4185
04085L2=L1
04095D7=D6
04105IF L3<L1 THEN 4165
04115L1=(L0+L2)/2
04125X0=L1
04135GOSUB 2865
04145D6=D4
04155GOTO 4185
04165L1=L3
04175D6=D8
04185NEXT K5
04195PRINT "NO ROOT FOUND IN ";K4;" ITERATIONS"
04205PRINT" PARAAMETER VALUES"
04215PRINT "NUI=";N1;" NUA=";N2;" KAPPAI=";N3;" KAPPAA=";N4;
04225PRINT " ZETAI=";N5;" ZETAA=";N6;" P=";M0
04235PRINT "VALUE OF FUNCTION D(x) IS ";D8;" AT X= ";L3
04245RETURN
04255REM END OF ROOTFINDER ROUTINE
04265REM *****************************************************************
04275REM INTEGRATION FOR EXPECTED LENGTH OF CHDR
04285REM INPUT G1 G2 M0 THRU M9 N1 THRU N6
04295REM INPUT T U V Z ARRAYS
04305REM OUTPUT L4
04315L4=0
04325FOR I8=1 TO 3
04335L5=0
04345FOR I9=1 TO 2
04355L6=M5+T(I8)*M7*SQR(M3/((1-U(I9))*G1*(N1-1)))
04365L7=M3*(M7*M7)*U(I9)/(1-U(I9))
04375N5=(M1*M5+M9*L6)/(M1+M9)
04385N6=(M2*M6+M9*L6)/(M2+M9)
04395N3=(M3*(M7*M7)+L7+G1*((L6-M5)**2))/(M1+M9)
04405N4=(M4*(M8*M8)+L7+G2*((L6-M6)**2))/(M2+M9)
04415GOSUB 2235
04425L5=L5+(V2-V1)*SQR(1-U(I9))*V(I9)/2
04435NEXT I9
04445L4=L4+L5*X(I8)/(2*SQR(1+T(I8)*T(I8)/(N1-1)))
04455NEXT I8
04465L4=L4*G3
04475RETURN
04485REM END OF INTEGRATION ROUTINE
04495REM ***************************************************
04505REM BETA CDF ROUTINE
04515REM INPUT A B J2
04525REM OUTPUT P
04535REM GOSUB'S TO BE CALLED PRIOR 5135 AND 5195
04545 DIM W(16),O(16)
04555GOSUB 5185
04565J1=0
04575IF A+B>85 THEN 5065
04585P=0
04595C6=0
04605IF A>1 THEN 4665
04615C6=A
04625C7=B
04635A=C7
04645B=C6
04655J2=1-J2
04665D0=(J2-J1)*.5
04675D1=(J1+J2)*.5
04685FOR I1=1 TO 16
04695D9=D0*O(I1)+D1
04705IF D9=0 THEN 4735
04715IF D9=1 THEN 4735
04725P=P+W(I1)*EXP(LOG(D9)*(A-1)+LOG(1-D9)*(B-1))
04735NEXT I1
04745P=P*F0
04755P=P*D0
04765IF C6=0 THEN 4815
04775A=C6
04785B=C7
04795P=1-P
04805J2=1-J2
04815RETURN
04825FOR I1=1 TO 16
04835READ W(I1),O(I1)
04845NEXT I1
04855DATA 2.71525E-02,-.989401
04865DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
04875DATA .124629,-.755404,.149596,-.617876
04885DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
04895DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
04905DATA .149596,.617876,.124629,.755404
04915DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
04925DATA .989401
04935RETURN
04945G9=A+B
04955GOSUB 5305
04965F0=G0
04975G9=A
04985GOSUB 5305
04995F0=F0-G0
05005G9=B
05015GOSUB 5305
05025F0=F0-G0
05035IF A+B>85 THEN 5055
05045F0=EXP(F0)
05055RETURN
05065W1=(B*J2)**(1/3)
05075W2=(A*(1-J2))**(1/3)
05085GOSUB 5155
05095I1=P
05105W1=(B*J1)**(1/3)
05115W2=(A*(1-J1))**(1/3)
05125GOSUB 5155
05135P=I1-P
05145RETURN
05155Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05165GOSUB 5955
05175RETURN
05185IF A+B>85 THEN 5215
05195D2=F0*J2**(A-1)*(1-J2)**(B-1)
05205RETURN
05215D2=F0+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05225IF D2>85 THEN 5255
05235D2=EXP(D2)
05245RETURN
05255D2=1.E+37
05265RETURN
05275REM
05285REM END OF BETA CDF ROUTINE
05295REM *******************************************************
05305REM ****************************************************
05315REM LOG GAMMA ROUTINE
05325REM INPUT G9
05335REM OUTPUT G0
05345G5=G9
05355IF G9 <= 1.E+30 THEN 5385
05365G0=1.E+38
05375RETURN
05385IF G9>1.E-09 THEN 5415
05395G0=0
05405RETURN
05415IF G9<1.E+10 THEN 5445
05425G0=G9*(LOG(G9)-1)
05435RETURN
05445G6=1
05455IF 18<G5 THEN 5495
05465G6=G6*G5
05475G5=G5+1
05485GOTO 5455
05495R8=1/G5**2
05505G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
05515C1=8.33333E-02
05525C2=2.77778E-03
05535C3=7.93651E-04
05545C4=5.95238E-04
05555G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
05565RETURN
05575REM END OF LOG GAMMA ROUTINE
05585REM ****************************************************
05595REM ****************************************************************
05605REM STUDENT'S T CDF ROUTINE
05615REM INPUT G J2
05625REM OUTPUT P
05635REM PRIOR GOSUB 5135
05645Y3=J2
05655IF G <= 120 THEN 5685
05665GOSUB 5955
05675GOTO 5835
05685IF G=1 THEN 5915
05695J1=0
05705N=G
05715P=0
05725D0=(J2-J1)*.5
05735D1=(J1+J2)*.5
05745FOR I1=1 TO 16
05755D9=D0*O(I1)+D1
05765IF D9=0 THEN 5795
05775IF D9=1 THEN 5795
05785P=P+W(I1)*(EXP(-(N+1)/2*LOG(1+D9*D9/N)))
05795NEXT I1
05805P=P*F0
05815P=P*D0
05825P=P+.5
05835RETURN
05845G9=(N+1)/2
05855GOSUB 5305
05865F0=G0
05875G9=N/2
05885GOSUB 5305
05895F0=EXP(F0-G0)/SQR(3.14159*N)
05905RETURN
05915REM FOLLOWING FOR NU=1
05925P=.5+1/3.14159*ATN(Y3)
05935RETURN
05945REM END STUDENT'S T CDF ROUTINE
05955REM **********************************************************
05965REM ROUTINE CALCULATES THE CDF FOR NORMAL DISTRIBUTION
05975REM INPUT Y3
05985REM OUTPUT P
05995REM
06005Y4=ABS(Y3)
06015X=Y3
06025T=1/(1+.231642*Y4)
06035D=.398942*EXP(-X*X/2)
06045C1=1.33027
06055C2=1.82126
06065C3=1.78148
06075C4=.356564
06085C5=.319382
06095P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
06105IF X >= 0 THEN 6125
06115P=1-P
06125RETURN
06135REM
06145REM END OF NORMAL CDF ROUTINE
06155REM **********************************************************
09000REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09035REM*************END ROUTINE
09999END