Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod58.bas
There are 2 other files named cmod58.bas in the archive. Click here to see a list.
00020 DIM W(4,3),C(70,3),S(4,7),A(5)
00030M1=1
00040REM ***************************************************************
00045REM CMOD58 CMOD58 CMOD58 CMOD58 CMOD58
00050M3=1
00060M2=0
00070REM ACTIVE ACTIVE ACTIVE ACTIVE ACTIVE ACTIVE
00080REM THIS ROUTINE USES A CHI SQUARE TO FIT A GAMMA, HENCE THE DF
00090REM USE BYE THIS PROGRAM ARE 2Z WHERE Z IS THE USERS SAMPLE SIZE
00100REM AN THE SCALE OF THE USER IS EQUAL TO 1/2(S(K5,7)), WHERE S(K5,7)
00110REM IS THE SCALE PARA FOR THE CHI SQUARE DIST
00120REM***************************************************************
00130REM
00140REM PRIOR FOR POISSON
00150REM
00160PRINT" 10TH 25TH 50TH 75TH 90TH"
00170A1=2.45262
00180B1=-2.47084
00190A2=2.01371
00200B2=-2.0064
00210Q7=0
00220S8=0
00230O3=0
00240 FILES RFILE1,RFILE2,RFILE3
00280RESTORE#1
00281 INPUT# 1,I1,I2,I3
00290SCRATCH#1
00291 PRINT # 1,58,I2,I3
00300MAT READ W
00310DATA 1,0,0,0,1,0,0,0,1,1,1,1
00320MAT READ A
00330DATA .1,.25,.5,.75,.9
00340MAT READ C
00350GOTO 450
00360REM******************************************************************
00370REM N1=DEGREES OF FREEDOM
00380REM l1 =scale factor
00390REM******************************************************************
00400PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT"
00410PRINT "1. PRIOR FOR BINOMINAL APPROXMATION."
00420PRINT "2. PRIOR FOR POISSON PROCESS."
00430GOSUB 9000
00440IF O1=1 THEN 460
00450M1=2
00460D9=0
00470S9=0
00480PRINT L$
00490IF M1=1 THEN 520
00500PRINT " PRIOR DISTRIBUTION FOR THE GAMMA POISSON."
00510GOTO 540
00520PRINT " PRIOR DISTRIBUTION FOR THE GAMMA POISSON BINOMINAL"
00530PRINT " APPROXIMATION."
00540PRINT
00550PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A GAMMA DISTRIBUTION TO"
00560PRINT "YOUR PRIOR BELIEFS ABOUT THE INTENSITY PARAMETER OF A POISSON"
00570PRINT "DISTRIBUTION."
00580PRINT
00590PRINT "WE BEGIN BY ASKING YOU TO SPECIFY THE 25TH, 50TH AND"
00600PRINT "75TH PERCENTILES OF YOUR PRIOR DISTRIBUTION."
00610PRINT
00620PRINT "SPECIFY 25TH PERCENTILE. YOUR BETTING ODDS ARE 3 TO 1 THE"
00630PRINT "INTENSITY IS GREATER THAN THIS VALUE IN A GIVEN INTERVAL."
00640PRINT "INPUT 25TH.";
00650GOSUB 9000
00660Q1=O1
00670IF Q1>0 THEN 700
00680PRINT "X MUST BE POSITIVE."
00690GOTO 620
00700PRINT
00710PRINT "SPECIFY 50TH PERCENTILE. YOUR BETTING ODDS ARE EVEN THE"
00720PRINT "INTENSITY IS LESS THAN THIS VALUE IN A GIVEN INTERVAL."
00730PRINT "INPUT 50TH.";
00740GOSUB 9000
00750Q2=O1
00760IF Q2>Q1 THEN 810
00770PRINT L$
00780PRINT "THE 50TH PERCENTILE MUST BE LARGER THAN THE 25TH,AND THE"
00790PRINT "75TH MUST BE LARGER THAN THE 50TH. PLEASE RESPECIFY."
00800GOTO 610
00810Q8=Q2
00820PRINT
00830PRINT "SPECIFY 75TH PERCENTILE. YOUR BETTING ODDS ARE 1 TO 3 THE"
00840PRINT "INTENSITY IS GREATER THAN THIS VALUE IN A GIVEN INTERVAL."
00850PRINT "INPUT 75TH.";
00860GOSUB 9000
00870Q3=O1
00880IF Q3>Q2 THEN 900
00890GOTO 770
00900PRINT L$
00910PRINT"FOUR POSSIBLE APPOXIMATE PRIOR DISTRBUTIONS ARE NOW BEING"
00920PRINT"COMPUTED FOR YOUR CONSIDERATION."
00930FOR K5=1 TO 4
00940O2=0
00950R9=1.E+30
00960FOR K6=3 TO 70
00970R=W(K5,1)*(Q1/Q2-C(K6,1))**2
00980R=R+W(K5,2)*(Q1/Q3-C(K6,2))**2
00990R=R+W(K5,3)*(Q2/Q3-C(K6,3))**2
01000IF R>R9 THEN 1050
01010R9=R
01020J(K5)=K6
01030NEXT K6
01040GOTO 1080
01050G=K6-1
01060S(K5,6)=G
01070R=0
01080IF K5=1 THEN 1140
01090IF K5=2 THEN 1290
01100IF K5=3 THEN 1380
01110IF K5=4 THEN 1490
01120GOTO 1740
01130REM ASSUME 25TH AND 50TH ARE MORE ACCURATELY ESTIMATED
01140GOSUB 1640
01150IF J(1) <> 70 THEN 1200
01160U1=Q2/Q1
01170V1=A1+B1*U1
01180G=1/V1/V1
01190S(K5,6)=G
01200GOSUB 3750
01210S1=Q1/S(K5,2)
01220S2=Q2/S(K5,3)
01230S=(S1+S2)/2
01240O2=1
01250GOSUB 3750
01260O2=0
01270S(K5,7)=S
01280GOTO 1710
01290IF J(2) <> 70 THEN 1320
01300G=1373.71*(Q1/Q3)-1021.13
01310S(K5,6)=G
01320GOSUB 3750
01330S1=Q1/S(K5,2)
01340S2=Q3/S(K5,4)
01350GOTO 1230
01360REM ASSUME 50TH AND 75TH ARE MOST ACCURATELY ESTIMATED
01370Q8=J(3)
01380IF J(3) <> 70 THEN 1430
01390U1=Q3/Q2
01400V1=A2+B2*U1
01410G=1/V1/V1
01420S(K5,6)=G
01430GOSUB 3750
01440S1=Q2/S(K5,3)
01450S2=Q3/S(K5,4)
01460GOTO 1230
01470GOTO 1630
01480REM ASSUME ALL ARE EQUALLY WELL ESTIMATED
01490IF S(3,6) <> 70 THEN 1570
01500U1=Q2/Q1
01510U2=Q3/Q2
01520V1=A1+B1*U1
01530V2=A2+B2*U2
01540V3=(V1+V2)/2
01550G=1/V3/V3
01560S(K5,6)=G
01570GOSUB 3750
01580S1=Q1/S(K5,2)
01590S2=Q2/S(K5,3)
01600S3=Q3/S(K5,4)
01610S=(S1+S2+S3)/3
01620GOTO 1240
01630REM
01640PRINT
01650PRINT "HERE ARE PERCENTILES OF FOUR GAMMA DISTRIBUTIONS FITTED"
01660PRINT "TO YOUR PERCENTILE SPECIFICATIONS."
01670PRINT
01680PRINT USING 1690
01690: 10TH 25TH 50TH 75TH 90TH
01700RETURN
01710PRINT USING 1720,K5,S(K5,1),S(K5,2),S(K5,3),S(K5,4),S(K5,5)
01720: # #####.## #####.## #####.## #####.## #####.##
01730NEXT K5
01740PRINT
01750PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE"
01760PRINT "WHICH MOST CLOSELY CORRESPONDS TO YOUR PRIOR BELIEFS. YOU"
01770PRINT "CAN EITHER TENTATIVELY ACCEPT THIS DISTRIBUTION OR SPECIFY"
01780PRINT "NEW VALUES FOR THE PERCENTILES."
01790PRINT
01800PRINT "IF YOU WANT ONE OF THE DISTRIBUTION TYPE ITS NUMBER.";
01810PRINT "IF YOU WANT TO RESPECIFY THE PERCENTILES TYPE '0'.";
01820GOSUB 9000
01830IF O1=0 THEN 620
01840GOTO 1880
01850PRINT "REENTER. INPUT MUST BE 0 OR THE NUMBER OF A DISTRIBUTION.";
01860GOSUB 9000
01870GOTO 1830
01880IF O1=1 THEN 1950
01890IF O1=2 THEN 1950
01900IF O1=3 THEN 1950
01910IF O1=4 THEN 1950
01920PRINT
01930K5=O1
01940GOTO 1850
01950K5=O1
01960IF M1=2 THEN 4350
01970GOSUB 2510
01980GOTO 2020
01990:FOUR POSSIBLE APPROXIMATE PRIOR DISTRIBUTIONS ARE NOW BEING
02000:COMPUTED FOR YOUR CONSIDERATION.
02010RETURN
02020PRINT "IF YOU DO NOT FEEL THAT THIS HYPOTHETICAL SAMPLE SIZE ( M )"
02030PRINT "REFLECTS YOUR PRIOR INFORMATION ABOUT THE STANDARD DEVIATION"
02040PRINT "YOU CAN SPECIFY A DIFFERENT ONE. A DIFFERENT M WILL NOT AFFECT"
02050PRINT "THE MEDIAN, BUT WILL CHANGE THE HDRS AND OTHER PERCENTILES. A"
02060PRINT "LARGER M WILL SHORTEN THE HDR INTERVALS, AND A SMALLER M WILL"
02070PRINT "LENGTHEN THEM."
02080PRINT
02090REM
02100PRINT "IF YOU WANT TO CHANGE M TYPE THE NEW VALUE (GREATER THAN 1.5)."
02110PRINT "IF YOU DO NOT WANT TO CHANGE M TYPE '0'."
02120REM
02130S9=1
02140GOSUB 9000
02150IF O1=0 THEN 2870
02160IF O1>1000 THEN 2200
02170IF O1 >= 1.5 THEN 2220
02180PRINT "REENTER. M MUST BE GREATER THAN 1.5."
02190GOTO 2120
02200PRINT "REENTER. M MUST NOT BE GREATER THAN 1000."
02210GOTO 2120
02220O2=0
02230P0=.5
02240G=2*O1
02250GOSUB 5000
02260S=S(K5,3)/J2
02270S(K5,7)=S
02280S(K5,6)=G
02290O2=1
02300GOSUB 3750
02310O2=0
02320GOSUB 2510
02330GOTO 2080
02340E1=0
02350E2=.65*L1
02360P0=.25
02370GOSUB 3390
02380S(K7,2)=X3
02390IF S9=1 THEN 2450
02400E1=X3
02410E2=.87*L1
02420P0=.5
02430GOSUB 3390
02440S(K7,3)=X3
02450E2=1.4*L1
02460E1=X3
02470P0=.75
02480GOSUB 3390
02490S(K7,4)=X3
02500K5=K7
02510PRINT L$
02520PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE GAMMA"
02530PRINT "DISTRIBUTION YOU ARE NOW CONSIDERING."
02540PRINT
02550PRINT USING 2570,S(K5,6)/2
02560M9=S(K5,7)*(S(K5,6)-2)
02570: DEGREES OF FREEDOM ###########.##
02580IF D9=0 THEN 2650
02590GOTO 2610
02600PRINT USING 2620,S(K5,6)/2
02610PRINT USING 2630,1/(2*(S(K5,7)))
02620: DEGREES OF FREEDOM ###########.##
02630: SCALE PARAMETER ###########.##
02640: MODE ###########.##
02650PRINT USING 2660,S(K5,2)
02660: 25TH PERCENTILE ###########.##
02670: 75TH PERCENTILE ###########.##
02680: 50TH (MEDIAN) #######.##
02690PRINT USING 2680,S(K5,3)
02700PRINT USING 2670,S(K5,4)
02710: 50% HDR ######.## -######.##
02720: 75% HDR ######.## -######.##
02730: 95% HDR ######.## -######.##
02740G=S(K5,6)
02750J5=.5
02760GOSUB 6620
02770PRINT USING 2710,J1*S(K5,7),J2*S(K5,7)
02780J5=.75
02790GOSUB 6620
02800PRINT USING 2720,J1*S(K5,7),J2*S(K5,7)
02810J5=.95
02820GOSUB 6620
02830PRINT USING 2730,S(K5,7)*J1,J2*S(K5,7)
02840IF D9=1 THEN 3310
02850PRINT
02860RETURN
02870PRINT
02880IF S8=1 THEN 2930
02890PRINT "YOU CAN CHANGE THE CENTERING OF THE DISTRIBUTION BY"
02900PRINT "SPECIFYING A DIFFERENT MEDIAN. THIS WILL NOT AFFECT"
02910PRINT "THE HYPOTHETICAL SAMPLE SIZE."
02920PRINT
02930PRINT "IF YOU WANT TO CHANGE THE MEDIAN TYPE THE NEW VALUE."
02940PRINT "IF YOU DO NOT TYPE '0'."
02950GOSUB 9000
02960IF O1=0 THEN 3100
02970M6=O1
02980IF O1>0 THEN 3020
02990PRINT "REENTER. INPUT MUST BE '0' OR POSITIVE VALUE FOR MEDIAN.";
03000GOSUB 9000
03010GOTO 2960
03020S8=1
03030S=O1/(S(K5,3)/S(K5,7))
03040S(K5,7)=S
03050O2=1
03060GOSUB 3750
03070O2=0
03080GOSUB 2510
03090GOTO 2930
03100PRINT L$
03110D9=1
03120PRINT
03130GOSUB 4350
03140PRINT "IF THE NUMBER OF OBSERVATIONS IMPLIED BY THIS DISTRIBUTION"
03150PRINT "ARE NOT CONGRUENT WITH YOUR PRIOR BELIVES YOU CAN CHANGE"
03160PRINT "THE NUMBER OF OBSERVATIONS. THIS WILL EFFECT THE PERCENTILES"
03170PRINT "BUT NOT THE NUMBER OF INTERVALS. IF YOU DO NOT WANT TO"
03180PRINT "CHANGE THE NUMBER OF OBSERVATIONS TYPE '0'.";
03190GOSUB 9000
03200IF O1=0 THEN 3250
03210G=O1*2
03220S(K5,6)=G
03230GOSUB 3750
03240GOTO 3120
03250PRINT
03260M1=0
03270PRINT L$
03280PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE GAMMA"
03290PRINT "DISTRIBUTION YOU ARE NOW CONSIDERING."
03300GOSUB 2600
03310PRINT
03320PRINT "THIS COMPLETES THE SPECIFICATION OF A PRIOR DISTRIBUTION"
03330PRINT "ON PI. YOU WILL NEED THE SCALE PARAMETER AND DEGRESS OF"
03340PRINT "FOR THE POSTERIOR ANALYSIS. WHEN YOU ARE READY TO"
03350PRINT "CONTINUE TYPE '1'.";
03360GOSUB 9000
03370IF O1 <> 1 THEN 3350
03380CHAIN "RSTRT"
03390REM*********************************************************
03400IF E1 <> 0 THEN 3430
03410P8=0
03420GOTO 3460
03430X3=E1
03440GOSUB 3520
03450P8=P
03460X3=E2
03470GOSUB 3520
03480P9=P
03490X3=.5*(E1+E2)
03500GOSUB 3520
03510GOTO 3580
03520X=L1*L1/X3/X3
03530G=N1
03540REMgosub5775
03550GOSUB 5500
03560P=1-P
03570RETURN
03580IF ABS(P-P0)<.001 THEN 3740
03590IF P<P0 THEN 3710
03600E2=X3
03610P9=P
03620P=(P0-P8)/(P9-P8)
03630IF P<.85 THEN 3660
03640P=.85
03650GOTO 3680
03660IF P>.15 THEN 3680
03670P=.15
03680X3=P*(E2-E1)+E1
03690GOSUB 3520
03700GOTO 3580
03710E1=X3
03720P8=P
03730GOTO 3620
03740RETURN
03750FOR I=1 TO 5
03760P0=A(I)
03770GOSUB 5000
03780IF O2 <> 1 THEN 3810
03790S(K5,I)=J2*S
03800GOTO 3820
03810S(K5,I)=J2
03820NEXT I
03830RETURN
03840PRINT
03850PRINT "IF THE NUMBER OF INTERVALS IMPLIED BY THIS DISTRIBUTION ARE"
03860PRINT "ARE NOT CONGRUENT WITH YOUR PRIOR BELIVES YOU CAN CHANGE THE"
03870O2=1
03880PRINT "NUMBER OF INTERVALS BY TYPING THE DESIRED NUMBER OF INTERVALS"
03890PRINT "IF YOU DO NOT WISH TO CHANGE THE NUMBER OF INTERVALS TYPE '0'.";
03900M3=0
03910GOSUB 9000
03920IF O1=0 THEN 3120
03930IF O1>0 THEN 3960
03940PRINT "INTERVALS MUST BE GREATER THEN 0, RENTER INTERVALS"
03950GOTO 3890
03960S=1/(2*O1)
03970IF M2=1 THEN 4050
03980M2=1
03990PRINT " "
04000PRINT "A CHANGE EFFECTS THE OBSERTVATIONS AND PERCENTILES AS FOLLOWS:"
04010PRINT "A)MORE INTERVALS OBSERVATIONS CONSTANT PERCENTILES LOWER."
04020PRINT "B)MORE INTERVALS PERCENTILES CONSTANT OBSERVATIONS HIGHER."
04030PRINT "C)LESS INTERVALS OBSERVATIONS CONSTANT PERCENTILES LOWER."
04040PRINT "D)LESS INTERVALS PERCENTILES CONSTANT OBSERVATIONS HIGHER"
04050PRINT
04060PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
04070PRINT "1. OBSERVATIONS CONSTANT"
04080PRINT "2. PERCENTILES CONSTANT";
04090GOSUB 9000
04100IF O1=2 THEN 4180
04110IF O1 <> 1 THEN 4060
04120G=S(K5,6)
04130S(K5,7)=S
04140GOSUB 3750
04150GOSUB 4350
04160M2=1
04170GOTO 3850
04180G=S(K5,3)/S+2/3
04190IF G>2 THEN 4280
04200PRINT "A POSSION PROCESS OF THIS TYPE IMPLIES LESS THEN 1 SUCCESS"
04210PRINT USING 4220,S
04220:IN #### INTERVALS. BY DEFINTION THE PROCESS IS DEFINED.
04230PRINT "BY THE TIME TILL SUCESS OCCUR. IF YOU WISH THE SUCCESS"
04240PRINT "TO BE SET TO ONE TYPE '1', ELSE '0'";
04250GOSUB 9000
04260IF O1 <> 1 THEN 4350
04270G=3
04280S(K5,7)=S
04290IF G<2000 THEN 4310
04300G=2000
04310S(K5,6)=G
04320GOSUB 4350
04330GOTO 3840
04340REM
04350PRINT L$
04360PRINT "HERE ARE SOME OF THE CHARACTERICS OF THE DISTRIBUTIONS"
04370PRINT "YOU ARE NOW CONSIDERING."
04380PRINT
04390PRINT USING 4400,S(K5,6)/2
04400: HYPOTHETICAL OBSERVATIONS######.##
04410PRINT USING 4420,1/(2*S(K5,7))
04420: HYPOTHETICAL INTERVALS#########.##
04430GOSUB 2650
04440IF M3=1 THEN 3840
04450RETURN
05000REM***********************************************************
05002REM CHI-SQUARE PERCENTILE
05004REM INPUTS P0
05006REM OUTPUTS J2
05008J=P0
05010P9=P0
05012P2=P0
05014GOSUB 5066
05016X=J2
05018GOSUB 5500
05020IF ABS(P-P9)<.0001 THEN 5064
05022P4=P
05024E7=J2
05026P2=P2-.009
05028GOSUB 5066
05030E1=J2
05032E2=E7
05034X=E1
05036GOSUB 5500
05038P3=P
05040X3=(P9-P3)/(P4-P3)
05042X3=X3*(E2-E1)+E1
05044X=X3
05046GOSUB 5500
05048IF ABS(P9-P)<.0001 THEN 5064
05050IF P>P9 THEN 5058
05052E1=X3
05054P3=P
05056GOTO 5040
05058E2=X3
05060P4=P
05062GOTO 5040
05064RETURN
05066P3=P2
05068IF P2 <= .5 THEN 5072
05070P2=1-P2
05072A1=SQR(LOG(1/P2/P2))
05074A2=2.51552+.802853*A1+.010328*A1*A1
05076A2=A2/(1+1.43279*A1+.189269*A1*A1+.001308*A1*A1*A1)
05078J2=A1-A2
05080IF P3 <= .5 THEN 5086
05082P2=P3
05084GOTO 5088
05086J2=-J2
05088J2=G*(1-2/9/G+J2*SQR(2/9/G)) ** 3
05090RETURN
05092REM
05094REM********************* END OF CSQ PERCENTILE ********************
05345REM ********************************************************
05350REM CHI-SQUARE DENSITY FUNCTION
05355REM ********************************************************
05360G9=.5*G
05365GOSUB 5860
05370F0=G0
05375F0=.5*(G-2)*LOG(J2)-.5*J2/S0/S0
05380D2=F0
05385RETURN
05500REM ********************************************************
05510REM CHI-SQUARE CDF ROUTINE-LOWER TAIL
05520REM INPUT G X
05530REM OUTPUT P
05540REM*************************
05550REM PEIZER-PRATT-APPROXIMATION
05560REM PROB(CHISQ<X) WITH DF=G
05570REM*************************
05580X2=X/2
05590T=X2-G*.5+1/3-.04/G
05591IF X2>0 THEN 5600
05592X2=.0001
05600Y=(G*.5-.5)/X2
05610IF Y<.00001 THEN 5680
05620IF ABS(Y-1)<.00001 THEN 5700
05630G6=(1-Y*Y+2*Y*LOG(Y))/(1-Y)**2
05640IF G6>-1 THEN 5710
05650Y3=T
05660GOTO 5720
05670GOTO 5710
05680G6=1
05690GOTO 5710
05700G6=0
05710Y3=T*SQR((1+G6)/X2)
05720GOSUB 8000
05740RETURN
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 ****************************************************
06600REM HDR ROUTINE FOR CHI-SQUARE
06605REM INPUTS G,J5
06610REM OUTPUTS J1,J2
06615REM*************************************************************
06620J8=6
06625GOSUB 6770
06630X=J1
06635GOSUB 5500
06640J3=P
06645X=J2
06650GOSUB 5500
06655J3=P-J3
06660IF ABS(J3-J5)>.00001 THEN 6670
06665RETURN
06670IF J3>J5 THEN 6685
06675J8=J8+5
06680GOTO 6625
06685J9=J8-5
06690J0=J8
06695J8=(J0+J9)/2
06700GOSUB 6770
06705X=J1
06710GOSUB 5500
06715J3=P
06720X=J2
06725GOSUB 5500
06730J3=P-J3
06735IF ABS(J3-J5)>.00001 THEN 6745
06740RETURN
06745IF J3>J5 THEN 6760
06750J9=J8
06755GOTO 6695
06760J0=J8
06765GOTO 6695
06770J=J8*(EXP(2*J8/(G-2))+1)/(EXP(2*J8/(G-2))-1)
06775J1=J-J8
06776IF J1>0 THEN 6780
06777J1=0
06780J2=J+J8
06785RETURN
06790REM
06795REM END OF CHID-SQUARE HDR ROUTINE
06800REM************************************************************
07000DATA .604875,.221152,.365616
07001DATA .417785,.208926,.50008
07002DATA .5131,.295542,.575994
07003DATA .57297,.357183,.62339
07004DATA .614846,.403747,.656664
07005DATA .646022,.440746,.682247
07006DATA .670508,.470813,.702174
07007DATA .690499,.496166,.718562
07008DATA .706955,.517977,.732688
07009DATA .72112,.536926,.744572
07010DATA .733521,.553554,.754653
07011DATA .74408,.568399,.763895
07012DATA .753802,.581965,.77204
07013DATA .761989,.594056,.779612
07014DATA .769689,.604903,.785906
07015DATA .776515,.614808,.791753
07016DATA .783131,.624412,.797328
07017DATA .788631,.632834,.802446
07018DATA .794293,.641007,.807016
07019DATA .798748,.648356,.811716
07020DATA .803904,.65569,.815632
07021DATA .807971,.662127,.819493
07022DATA .811755,.667953,.82285
07023DATA .815545,.674144,.826617
07024DATA .819322,.67935,.829161
07025DATA .822513,.684764,.832526
07026DATA .825623,.689512,.835141
07027DATA .828536,.694323,.838012
07028DATA .83172,.699103,.840551
07029DATA .834176,.703491,.843337
07030DATA .836875,.707275,.845139
07031DATA .83925,.71146,.847734
07032DATA .84186,.714982,.849289
07033DATA .844066,.718803,.851596
07034DATA .846497,.722458,.853468
07035DATA .848492,.725841,.855449
07036DATA .85065,.728931,.85691
07037DATA .85234,.732293,.859156
07038DATA .85418,.734736,.860166
07039DATA .856125,.73795,.861966
07040DATA .857561,.740817,.863865
07041DATA .859136,.743314,.865188
07042DATA .860804,.745969,.866596
07043DATA .861874,.748785,.868787
07044DATA .863724,.751111,.86962
07045DATA .86501,.753011,.870523
07046DATA .866328,.755632,.872223
07047DATA .86778,.757816,.873281
07048DATA .869266,.760053,.874362
07049DATA .870141,.761824,.875518
07050DATA .87179,.764343,.876751
07051DATA .872764,.766295,.87801
07052DATA .87381,.768382,.879346
07053DATA .874929,.769872,.879925
07054DATA .876071,.772134,.88136
07055DATA .877287,.773786,.882022
07056DATA .877775,.775564,.883557
07057DATA .879084,.777378,.884304
07058DATA .88047,.778547,.88424
07059DATA .881054,.780518,.885891
07060DATA .882537,.78262,.886785
07061DATA .883254,.783263,.886793
07062DATA .883983,.784719,.887708
07063DATA .884726,.787013,.889556
07064DATA .885533,.788622,.890563
07065DATA .886353,.790261,.891587
07066DATA .88724,.791188,.89174
07067DATA .88814,.79297,.892843
07068DATA .889056,.793935,.893008
07069DATA .890039,.795864,.89419
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
09999END