Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmodaa.bas
There are 2 other files named cmodaa.bas in the archive. Click here to see a list.
00020REM**********************************************************
00030REM CMODAA CMODAA CMODAA CMODAA CMODAA CMODAA
00040REM**********************************************************
00050REM
00060REM THIS MODULE ASSISTS THE USER IN SPECIFYING A BETA DISTRIBUTION
00070REM AS A PRIOR DISTRIBUTION ON PI.
00080REM
00100REM
00110REM****************************************************
00120 FILES RFILE1,RFILE2,RFILE3
00160RESTORE#1
00161 INPUT# 1,I1,I2,I3
00170SCRATCH#1
00171 PRINT # 1,36,I2,I3
00180S9=0
00190V7=0
00195V2=1
00200S8=0
00210GOSUB 5160
00220 DIM S(4,5)
00230PRINT L$
00240PRINT " PRIOR DISTRIBUTION - BETA PASCAL MODEL"
00250PRINT
00260PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A BETA DISTRIBUTION"
00270PRINT "TO YOUR PRIOR BELIEFS ABOUT PI. WE BEGIN BY ASKING YOU TO"
00280PRINT "SPECIFY THE 25TH, 50TH AND 75TH PERCENTILES OF YOUR PRIOR"
00290PRINT "DISTRIBUTION."
00300PRINT
00310PRINT "SPECIFY 25TH. YOUR BETTING ODDS ARE 3 TO 1 THAT PI IS"
00320PRINT "GREATER THAN THIS VALUE.";
00330GOSUB 9000
00340IF O1>0 THEN 370
00350PRINT "REENTER. 25TH MUST GREATER THAN 0 AND LESS THAN 1."
00360GOTO 330
00370IF O1 >= 1 THEN 350
00380Q1=O1
00390PRINT
00400PRINT "SPECIFY 50TH. YOUR BETTING ODDS ARE EVEN THAT PI IS"
00410PRINT "GREATER THAN THIS VALUE.";
00420GOSUB 9000
00430IF O1>Q1 THEN 470
00440PRINT "REENTER. 50TH MUST BE > THAN 25TH BUT NOT GREATER THAN .97."
00450GOTO 420
00460GOTO 3500
00470IF O1>.97 THEN 440
00480Q2=O1
00490PRINT
00500PRINT "SPECIFY 75TH. YOUR BETTING ODDS ARE 1 TO 3 THAT PI IS"
00510PRINT "GREATER THAN THIS VALUE.";
00520GOSUB 9000
00530IF O1>Q2 THEN 560
00540PRINT "REENTER. 75TH MUST BE GREATER THAN 50TH AND LESS THAN 1."
00550GOTO 520
00560IF O1 >= 1 THEN 540
00570Q3=O1
00580REM*****************************************************************
00590REM
00600REM Q1 25TH PERCENTILE
00610REM Q2 50TH PERCENTILE
00620REM Q3 75TH PERCENTILE
00630PRINT
00640PRINT "POSSIBLE APPROXIMATE DISTRIBUTIONS ARE BEING COMPUTED."
00650X=SQR(Q2*(1-Q1))-SQR(Q1*(1-Q2))
00660X=X*X
00670C=.114/X
00680X7=1/X
00690A=C*Q2+1/3
00700B=C*(1-Q2)+1/3
00710K5=1
00720GOSUB 1610
00730REM
00740REM CHECK TO SEE IT THE FIT IS RIGHT ON
00750REM
00760IF ABS((S(1,2)-Q1))<.01 THEN 780
00770GOTO 850
00780IF ABS((S(1,3)-Q2))<.01 THEN 800
00790GOTO 850
00800IF ABS((S(1,4)-Q3)) >= .01 THEN 850
00810GOTO 2510
00820REM
00830REM FIT IS NOT RIGHT ON SO PRINT 4 APPROXIMATIONS
00840REM
00850PRINT L$
00860PRINT "HERE ARE SOME OF THE PERCENTILES OF FOUR BETA DISTRIBUTIONS"
00870PRINT "THAT HAVE BEEN FITTED TO YOUR PERCENTILE SPECIFICATIONS."
00880PRINT
00890PRINT " 10TH 25TH 50TH 75TH 90TH"
00900GOSUB 2480
00910X=SQR(Q2*(1-Q3))-SQR(Q3*(1-Q2))
00920X=X*X
00930C=.114/X
00940A=C*Q2+1/3
00950B=C*(1-Q2)+1/3
00960K5=2
00970GOSUB 1610
00980C2=.057*(1/X+X7)
00990A=C2*Q2+1/3
01000B=C2*(1-Q2)+1/3
01010K5=3
01020GOSUB 1610
01030A=5*P(1,1)+3*P(2,1)+P(3,1)
01040A=A/9
01050B=5*P(1,2)+3*P(2,2)+P(3,2)
01060B=B/9
01070K5=4
01080GOSUB 1610
01090PRINT
01100GOSUB 4390
01110GOSUB 9000
01120IF O1=0 THEN 1140
01130GOTO 1450
01140PRINT
01150PRINT "RESPECIFY PERCENTILES."
01160PRINT
01170PRINT "25TH";
01180GOSUB 9000
01190IF O1>0 THEN 1220
01200PRINT "REENTER. 25TH MUST BE GREATER THAN 0 AND LESS THAN 1."
01210GOTO 1180
01220IF O1 >= 1 THEN 1200
01230Q1=O1
01240PRINT "50TH";
01250GOSUB 9000
01260IF O1>Q1 THEN 1290
01270PRINT "REENTER. 50TH MUST BE GREATER THAN 25TH AND LESS THAN 1."
01280GOTO 1250
01290IF O1 >= 1 THEN 1270
01300Q2=O1
01310PRINT "75TH";
01320GOSUB 9000
01330IF O1>Q2 THEN 1360
01340PRINT "REENTER. 75TH MUST BE GREATER THAN 50TH AND LESS THAN 1."
01350GOTO 1310
01360IF O1 >= 1 THEN 1340
01370Q3=O1
01380GOTO 650
01390REM*************************************************************
01400REM
01410REM USER DOES NOT WANT TO RESPECIFY PERCENTILES
01420REM PICKS ONE OF FOUR APPROXIMATE DISTRIBUTIONS AS
01430REM TENTATIVE FIT
01440REM
01450IF O1=1 THEN 1520
01460IF O1=2 THEN 1520
01470IF O1=3 THEN 1520
01480IF O1=4 THEN 1520
01490PRINT "REENTER. MUST BE 0 OR NUMBER OF A DISTRIBUTION."
01500GOSUB 9000
01510GOTO 1120
01520K5=O1
01530A=P(K5,1)
01540B=P(K5,2)
01550GOTO 3830
01560REM**********************************************************
01570REM
01580REM IF EITHER PARAMETER IS NOT GREATER THAN 1.15
01590REM ADJUSTMENTS ARE MADE TO BOTH SO THAT BOTH ARE
01600REM
01610IF A>B THEN 2180
01620IF A>1.15 THEN 2230
01630IF K5=4 THEN 1830
01640IF K5=3 THEN 1750
01650REM A IS LESS THAN 1.15
01660REM
01670REM SET A=1.15*(K5-1)*.4 AND LET Q2 BE AN ESTIMATE OF MEAN
01680REM
01690B=B*(1.15+(K5-1)*.4)/A
01700A=1.15+(K5-1)*.4
01710GOTO 2230
01720REM
01730REM SET A=1.15 AND LET Q2 BE ESTIMATE OF MODE
01740REM
01750B=.15/(Q2)+.85
01760A=1.15
01770GOTO 2230
01780REM
01790REM SET M0=A+B WHERE B IS ESTIMATED TWICE-ONCE ASSUMING
01800REM Q2=MEAN AND AGAIN ASSUMING Q2=MODE KEEP UPPING M0
01810REM UNTIL BOTH PARAMETERS ARE GREATER THAN 1.15
01820REM
01830A=1.15
01840M0=A/Q2
01850M0=.5*M0+.5*A/Q2+1
01860B=M0-A
01870IF B>1.15 THEN 2230
01880A=A+.1
01890GOTO 1840
01900REM******************************************************
01910REM
01920REM B IS LESS THAN 1.15
01930REM
01940IF K5=3 THEN 2050
01950IF K5=4 THEN 2150
01960REM
01970REM SET B=1.15 AND LET Q2 BE ESTIMATE OF MEAN
01980A=A*(1.15+(K5-1)*.4)/B
01990B=1.15+(K5-1)*.4
02000GOTO 2230
02010REM
02020REM SET M0=A+B WHERE IT IS ESTIMATED TWICE-ONCE UNDER
02030REM ASSUMPTION Q2=MEAN AND AGAIN ASSUMING Q2=MODE AVERAGE
02040REM
02050B=1.15
02060M0=B/(1-Q2)
02070M0=.5*M0+.5*((2*Q2-B-1)/(Q2-1))
02080A=M0-B
02090IF A>1.15 THEN 2230
02100B=B+.1
02110GOTO 2050
02120REM
02130REM SET B=1.15 AND LET Q2=MODE
02140REM
02150B=1.15
02160A=(1-.85*Q2)/(1-Q2)
02170GOTO 2230
02180IF B<1.15 THEN 1940
02190REM *********************************************************
02200REM
02210REM COMPUTE PERCENTILES
02220REM
02230P1=.1
02240P(K5,1)=A
02250P(K5,2)=B
02260GOSUB 4120
02270S(K5,1)=J2
02280P1=.9
02290GOSUB 4150
02300S(K5,5)=J2
02310P1=.75
02320E1=S(K5,1)
02330A7=.1
02340GOSUB 4180
02350S(K5,4)=J2
02360P1=.25
02370E1=S(K5,1)
02380A7=.1
02390GOSUB 4180
02400S(K5,2)=J2
02410P1=.5
02420E2=S(K5,4)
02430A8=.75
02440GOSUB 4180
02450S(K5,3)=J2
02460IF S9=1 THEN 2540
02470IF K5=1 THEN 2500
02480PRINT USING 2490,K5,S(K5,1),S(K5,2),S(K5,3),S(K5,4),S(K5,5)
02490: ## ##.## ##.## ##.## ##.## ##.##
02500RETURN
02510PRINT L$
02520PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE BETA"
02530PRINT "DISTRIBUTION FITTED TO YOUR PERCENTILE SPECIFICATIONS."
02540PRINT
02550M=A+B
02560: HYPOTHETICAL SAMPLE SIZE (M)######.##
02570: STANDARD DEVIATION ##.##
02580: 10TH PERCENTILE ##.##
02590: 25TH PERCENTILE ##.##
02600: 50TH (MEDIAN) ##.##
02610: 75TH PERCENTILE ##.##
02620: 90TH PERCENTILE ##.##
02630: 50% HDR ##.## -##.##
02640: 75% HDR ##.## -##.##
02650: MODE ##.##
02660: 95% HDR ##.## -##.##
02670: MEAN ##.##
02680: PARAMETER A ######.##
02690: PARAMETER B ######.##
02700J5=.5
02710J8=S(K5,4)/S(K5,2)
02720GOSUB 7000
02730IF V7=1 THEN 2760
02740PRINT USING 2560,A+B
02750IF V7=0 THEN 2770
02760PRINT USING 2650,(A-1)/(M-2)
02770PRINT USING 2580,S(K5,1)
02780PRINT USING 2590,S(K5,2)
02790PRINT USING 2600,S(K5,3)
02800PRINT USING 2610,S(K5,4)
02810PRINT USING 2620,S(K5,5)
02820IF V7=1 THEN 2850
02830IF J2<.99 THEN 2840
02835J2=.99
02840GOTO 2870
02850J1=H(1)
02860J2=H(2)
02870PRINT USING 2630,J1,J2
02880IF V7=1 THEN 2950
02890H(1)=J1
02900H(2)=J2
02910J5=.75
02920GOSUB 7000
02930IF J2<.99 THEN 2940
02935J2=.99
02940GOTO 2970
02950J1=H(3)
02960J2=H(4)
02970PRINT USING 2640,J1,J2
02980IF V7=1 THEN 3050
02990H(3)=J1
03000H(4)=J2
03010J5=.95
03020GOSUB 7000
03030IF J2<.99 THEN 3040
03035J2=.99
03040GOTO 3070
03050J1=H(5)
03060J2=H(6)
03070PRINT USING 2660,J1,J2
03080IF V7=1 THEN 3880
03090H(5)=J1
03100H(6)=J2
03110IF S9=1 THEN 3175
03120PRINT "IF YOU DO NOT FEEL THAT THE HYPOTHETICAL SAMPLE SIZE (M)"
03130PRINT "REFLECTS YOUR PRIOR INFORMATION ABOUT PI YOU CAN SPECIFY A"
03140PRINT "DIFFERENT VALUE FOR M. THIS WILL NOT AFFECT THE MEDIAN BUT"
03150PRINT "WILL CHANGE THE HDRS AND OTHER PERCENTILES. A LARGER M WILL"
03160PRINT "RESULT IN SHORTER INTERVALS, AND A SMALLER M IN LONGER ONES."
03170PRINT
03175IF V2=0 THEN 3300
03176V2=0
03180M2=2
03190M7=M2/S(K5,3)
03200M8=(M2-1)/S(K5,3)+2
03210IF M7>M8 THEN 3230
03220M7=M8
03230IF M7-M2>2 THEN 3260
03240M2=M2+1
03250GOTO 3190
03260M7=INT(M7+1)
03270IF M7<M THEN 3290
03280M7=M
03290REM
03300:IF YOU WANT TO CHANGE M TYPE NEW VALUE (AT LEAST###.##).
03310PRINT USING 3300,M7
03320PRINT "IF YOU DO NOT WANT TO CHANGE M TYPE '0'";
03330GOSUB 9000
03340IF O1=0 THEN 3410
03350IF O1 >= M7 THEN 3380
03360PRINT "REENTER. INPUT MUST BE 0 OR AN ACCEPTABLE M VALUE."
03370GOTO 3330
03380M=O1
03381IF M<25 THEN 3390
03382A=S(K5,3)*M
03383B=M-A
03384GOTO 3750
03390O1=S(K5,3)
03400GOTO 3560
03410PRINT
03420IF S8=0 THEN 3440
03430GOTO 3470
03440PRINT "TO CHANGE THE CENTERING OF THE DISTRIBUTION, SPECIFY"
03450S8=1
03460PRINT "A DIFFERENT MEDIAN. THIS WILL NOT AFFECT THE VALUE OF M."
03470REM
03480PRINT "IF YOU WANT TO CHANGE MEDIAN TYPE NEW VALUE ELSE '0'.";
03490GOSUB 9000
03500IF O1=0 THEN 4040
03510IF O1>0 THEN 3540
03520PRINT "REENTER. INPUT MUST BE 0, OR VALUE BETWEEN 0 AND 1."
03530GOTO 3490
03540V2=1
03542IF O1<1 THEN 3560
03550GOTO 3520
03560L4=O1*M
03570L3=O1*(M-2)+1
03575IF L4>L3 THEN 3580
03576A=L3
03577L3=L4
03578L4=A
03580A=.5*(L3+L4)
03582B=M-A
03585GOTO 3630
03590L4=A
03600L3=A0
03605GOTO 3630
03610L4=A0
03620L3=A
03624A=.5*(L3+L4)
03630S9=1
03640P1=.5
03650B=M-A
03660GOSUB 4120
03670IF ABS(J2-O1)<.004 THEN 3750
03680IF J2>O1 THEN 3720
03690L3=A
03700A=.5*(L4+A)
03710GOTO 3650
03720L4=A
03730A=.5*(L3+A)
03740GOTO 3650
03750IF A<1.15 THEN 3770
03760IF B>1.15 THEN 3830
03770PRINT "THE MEDIAN YOU SPECIFIED IS NOT CONSISTENT WITH THE M VALUE"
03780PRINT "YOU ARE USING. EXTREME VALUES (HIGH OR LOW) REQUIRE LARGER"
03790PRINT "M VALUES. YOU SHOULD EITHER MODERATE YOUR MEDIAN ESTIMATE"
03800PRINT "OR INCREASE THE M. WHEN YOU WANT TO CONTINUE TYPE '1'."
03810GOSUB 9000
03820GOTO 1530
03830PRINT L$
03840PRINT "HERE ARE SOME CHARACTERISTICS OF THE BETA DISTRIBUTION YOU"
03850PRINT "ARE NOW CONSIDERING."
03860IF S9=0 THEN 2540
03870GOTO 1610
03880REM
03890PRINT "THIS COMPLETES THE FITTING OF A PRIOR DISTRIBUTION ON PI."
03900PRINT "IF YOU ARE NOT GOING DIRECTLY TO THE POSTERIOR ANALYSIS YOU"
03910PRINT "SHOULD RECORD THE PARAMETERS."
03920PRINT
03930PRINT "IF YOU WANT TO GO TO THE POSTERIOR ANALYSIS TYPE '1'."
03940PRINT "IF YOU DO NOT TYPE '0'.";
03950GOSUB 9000
03960IF O1=1 THEN 4020
03970IF O1=0 THEN 4010
03980PRINT
03990PRINT "REENTER. INPUT MUST BE 0 OR 1."
04000GOTO 3950
04010CHAIN "RSTRT"
04020SCRATCH#2
04021 PRINT # 2,A,B
04030CHAIN "CMODAC"
04040V7=1
04050PRINT L$
04060PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE PRIOR"
04070PRINT "DISTRIBUTION FITTED TO YOUR PRIOR BELIEFS ABOUT PI."
04080PRINT
04090PRINT USING 2680,A
04100PRINT USING 2690,B
04110GOTO 2560
04120E1=0
04130A7=0
04140GOSUB 5220
04150E2=1
04160A8=.999999
04170A6=P1
04180S5=(A6-A7)/(A8-A7)
04190IF S5<.85 THEN 4210
04200S5=.85
04210IF S5>.15 THEN 4230
04220S5=.15
04230J2=E1+S5*(E2-E1)
04240J1=0
04250GOSUB 5000
04260IF ABS(P-P1)<.005 THEN 4360
04270IF J2>.99 THEN 4370
04280IF P>P1 THEN 4320
04290E1=J2
04300A7=P
04310GOTO 4180
04320A8=P
04330E2=J2
04340GOTO 4180
04350GOTO 4360
04360IF J2<.99 THEN 4380
04370J2=.99
04380RETURN
04390PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE"
04400PRINT "WHICH MOST CLOSELY CORRESPONDS TO YOUR PRIOR BELIEFS."
04410PRINT "YOU CAN EITHER TENTATIVELY ACCEPT THIS DISTRIBUTION OR"
04420PRINT "RESPECIFY THE PERCENTILES."
04430PRINT
04440PRINT "IF YOU WANT ONE OF THESE DISTRIBUTIONS TYPE ITS NUMBER."
04450PRINT "IF YOU WANT TO RESPECIFY THE PERCENTILES TYPE '0'."
04460RETURN
05000REM ***************************************************
05005REM BETA CDF ROUTINE
05010REM INPUT A B J2
05015REM OUTPUT P
05020REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220
05025 DIM W(16),O(16)
05030GOSUB 5340
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A>1 THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=0 THEN 5115
05105IF D9=1 THEN 5115
05110P=P+W(I1)*EXP(LOG(D9)*(A-1)+LOG(1-D9)*(B-1))
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=0 THEN 5155
05135A=C6
05140B=C7
05145P=1-P
05150J2=1-J2
05155RETURN
05160FOR I1=1 TO 16
05165READ W(I1),O(I1)
05170NEXT I1
05175DATA 2.71525E-02,-.989401
05180DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
05185DATA .124629,-.755404,.149596,-.617876
05190DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
05195DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
05200DATA .149596,.617876,.124629,.755404
05205DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
05210DATA .989401
05215RETURN
05220G9=A+B
05225GOSUB 5850
05230F0=G0
05235G9=A
05240GOSUB 5850
05245F0=F0-G0
05250G9=B
05255GOSUB 5850
05260F0=F0-G0
05265IF A+B>85 THEN 5275
05270F0=EXP(F0)
05275RETURN
05280W1=(B*J2)**(1/3)
05285W2=(A*(1-J2))**(1/3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)**(1/3)
05305W2=(A*(1-J1))**(1/3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05340REM 2/16/76 CHANGED TO ALL LOG
05345D2=F0+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350IF D2<-80 THEN 5370
05355IF D2>85 THEN 5380
05360D2=EXP(D2)
05365RETURN
05370D2=1.E-37
05375RETURN
05380D2=1.E+37
05385RETURN
05390REM
05395REM END OF BETA CDF ROUTINE
05400REM *******************************************************
05850REM ****************************************************
05852REM LOG GAMMA ROUTINE
05853REM INPUT G9
05854REM OUTPUT G0
05860G5=G9
05863IF G9 <= 1.E+30 THEN 5872
05866G0=1.E+38
05869RETURN
05872IF G9>1.E-09 THEN 5881
05875G0=0
05878RETURN
05881IF G9<1.E+10 THEN 5890
05884G0=G9*(LOG(G9)-1)
05887RETURN
05890G6=1
05893IF 18<G5 THEN 5905
05896G6=G6*G5
05899G5=G5+1
05902GOTO 5893
05905R8=1/G5**2
05908G0=(G5-.5)*LOG(G5)-G5+.918939-LOG(G6)
05911C1=8.33333E-02
05914C2=2.77778E-03
05917C3=7.93651E-04
05920C4=5.95238E-04
05923G0=G0+1/G5*(C1-(R8*(C2+(R8*(C3-(R8*(C4)))))))
05926RETURN
05927REM END OF LOG GAMMA ROUTINE
05928REM ****************************************************
07000REM**************************************************************
07002REM BETA HDR ROUTINE
07004REM INPUT A B J5
07006REM OUTPUT J1 J2
07008J8=2
07010GOSUB 5220
07012REM
07014J7=1
07016GOSUB 7082
07018IF J8<-99 THEN 7028
07020P9=0
07022GOSUB 5000
07024P0=P
07026IF ABS(P0-J5)>.001 THEN 7030
07028RETURN
07030IF P0>J5 THEN 7038
07031J7=J8
07032J8=J8*J8
07034P9=P0
07036GOTO 7016
07038J9=J7
07040J0=J8
07042J8=(J5-P9)/(P0-P9)
07044IF J8>.15 THEN 7050
07046J8=.15
07048GOTO 7054
07050IF J8<.85 THEN 7054
07052J8=.85
07054J8=J8*(J0-J9)+J9
07056GOSUB 7082
07058IF J8<-99 THEN 7066
07060GOSUB 5000
07062J3=P
07064IF ABS(J3-J5)>.001 THEN 7068
07066RETURN
07068IF J3>J5 THEN 7076
07070J9=J8
07072P9=J3
07074GOTO 7042
07076J0=J8
07078P0=J3
07080GOTO 7042
07082J=(A-1)/(B-1)
07084J1=(J8**J-1)/(J8**(J+1)-1)
07086J2=J8*J1
07088IF J2<1 THEN 7094
07090J8=J8*.95
07092GOTO 7084
07094RETURN
07096REM END OF BETA HDR ROUTINE
07098REM*************************************************
08000REM **********************************************************
08001REM ROUTINE CALCULATES THE CDF FOR NORMAL DISTRIBUTION
08002REM INPUT Y3
08003REM OUTPUT P
08004REM
08005Y4=ABS(Y3)
08010X1=X
08015X=Y3
08020T=1/(1+.231642*Y4)
08025D=.398942*EXP(-X*X/2)
08030C1=1.33027
08035C2=1.82126
08040C3=1.78148
08045C4=.356564
08050C5=.319382
08055P=1-D*T*((((C1*T-C2)*T+C3)*T-C4)*T+C5)
08060IF X >= 0 THEN 8070
08065P=1-P
08070X=X1
08075RETURN
08076REM
08077REM END OF NORMAL CDF ROUTINE
08078REM **********************************************************
09000REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09035REM*************END ROUTINE
09999 END