Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod4.bas
There are 2 other files named cmod4.bas in the archive. Click here to see a list.
00020REM ***************************************************************
00030REM CMOD4 CMOD4 CMOD4 CMOD4 CMOD4 CMOD4
00040REM***************************************************************
00050REM
00060REM PRIOR FOR STANDARD DEVIATION OF TWO PARAMETER NORMAL
00070REM
00080REM************************************************************
00090S8=0
00100FILES RFILE1,RFILE2,RFILE3,RF4,RF5,RF6,RF7,RF8,RF9
00140RESTORE#1
00141 INPUT# 1,I1,I2,I3
00150SCRATCH#1
00151 PRINT # 1,4,I2,I3
00155S8=0
00160DATA 1,0,0,0,1,0,0,0,1,3,2,1
00170 DIM J(4),S(4,5),W(4,3),V(76),C(72,2)
00180MAT READ W
00190REM******************************************************************
00200REM N1=DEGREES OF FREEDOM
00210REM l1 =scale factor
00220REM******************************************************************
00230D9=0
00240S9=0
00250MAT READ C
00260FOR K5=1 TO 72
00270K6=INT(C(K5,1)*.1)
00280V(K5+4)=K6/1000
00290C(K5,1)=C(K5,1)-10*K6
00300NEXT K5
00310PRINT L$
00320PRINT " PRIOR DISTRIBUTION ON THE STANDARD DEVIATION"
00330PRINT
00340PRINT "THIS MODULE WILL ASSIST YOU IN FITTING AN INVERSE CHI DISTRIBUTION"
00350PRINT "TO YOUR PRIOR BELIEFS ABOUT THE STANDARD DEVIATION OF A NORMAL"
00360PRINT "DISTRIBUTION"
00370PRINT
00380PRINT "WE BEGIN BY ASKING YOU TO SPECIFY THE 25TH, 50TH AND"
00390PRINT "75TH PERCENTILES OF YOUR PRIOR DISTRIBUTION."
00400PRINT
00410PRINT "SPECIFY 25TH PERCENTILE. YOUR BETTING ODDS ARE 3 TO 1 THAT THE "
00420PRINT "STANDARD DEVIATION IS GREATER THAN THIS VALUE. INPUT 25TH.";
00430GOSUB 9000
00440Q1=O1
00450IF Q1>0 THEN 480
00460PRINT "STANDARD DEVIATION MUST BE POSITIVE."
00470GOTO 410
00480PRINT
00490PRINT "SPECIFY 50TH PERCENTILE. YOUR BETTING ODDS ARE EVEN THAT THE"
00500PRINT "STANDARD DEVIATION IS LESS THAN THIS VALUE. INPUT 50TH.";
00510GOSUB 9000
00520Q2=O1
00530IF Q2>Q1 THEN 580
00540PRINT L$
00550PRINT "THE 50TH PERCENTILE MUST BE LARGER THAN THE 25TH,AND THE"
00560PRINT "75TH MUST BE LARGER THAN THE 50TH. PLEASE RESPECIFY."
00570GOTO 400
00580Q8=Q2
00590PRINT
00600PRINT "SPECIFY 75TH PERCENTILE. YOUR BETTING ODDS ARE 1 TO 3 THAT THE"
00610PRINT "STANDARD DEVIATION IS GREATER THAN THIS VALUE. INPUT 75TH.";
00620GOSUB 9000
00630Q3=O1
00640IF Q3>Q2 THEN 660
00650GOTO 540
00660PRINT
00670PRINT"FOUR POSSIBLE APPROXIMATE PRIOR DISTRIBUTIONS ARE NOW BEING"
00680PRINT"FOR YOUR CONSIDERATION."
00690FOR K5=1 TO 4
00700R9=1.E+30
00710FOR K6=1 TO 72
00720R=W(K5,1)*(Q2/Q1-C(K6,1))**2
00730R=R+W(K5,2)*(Q3/Q1-C(K6,2))**2
00740R=R+W(K5,3)*(Q2/Q3-C(K6,1)/C(K6,2))**2
00750IF R>R9 THEN 920
00760R9=R
00770IF K6 <= 16 THEN 840
00780IF K6 <= 28 THEN 880
00790IF K6 <= 32 THEN 820
00800V(K5)=K6-22
00810GOTO 890
00820V(K5)=8+(K6-28)*.5
00830GOTO 890
00840V(K5)=2+(K6-1)*.2
00850IF V(K5) <> 2 THEN 870
00860V(K5)=2.01
00870GOTO 890
00880V(K5)=5+(K6-16)*.25
00890J(K5)=K6
00900NEXT K6
00910GOTO 930
00920IF K5=4 THEN 900
00930IF K5=1 THEN 1000
00940IF K5=2 THEN 1090
00950IF K5=3 THEN 1220
00960IF K5=4 THEN 1310
00970NEXT K5
00980GOTO 1480
00990REM ASSUME 25TH AND 50TH ARE MORE ACCURATELY ESTIMATED
01000IF J(1) <> 72 THEN 1060
01010K=(Q2-Q1)/.64
01020GOSUB 2310
01030J(1)=L1
01040V(1)=N1
01050GOTO 1350
01060J(1)=Q1/V(J(1)+4)
01070GOTO 1350
01080REM ASSUME 25TH AND 75TH ARE MOST ACCURATELY ESTIMATED
01090IF J(2) <> 72 THEN 1180
01100Q6=Q2
01110Q2=Q1+(.64/1.34)*(Q3-Q1)
01120K=(Q2-Q1)/.64
01130GOSUB 2310
01140J(2)=L1
01150V(2)=N1
01160Q2=Q6
01170GOTO 1350
01180J(2)=Q3/(C(J(2),2)*V(J(2)+4))
01190GOTO 1350
01200REM ASSUME 50TH AND 75TH ARE MOST ACCURATELY ESTIMATED
01210Q8=J(3)
01220IF J(3) <> 72 THEN 1280
01230K=(Q3-Q2)/.695
01240GOSUB 2310
01250V(3)=N1
01260J(3)=L1
01270GOTO 1350
01280J(3)=Q2/(C(J(3),1)*V(J(3)+4))
01290GOTO 1350
01300REM ASSUME ALL ARE EQUALLY WELL ESTIMATED
01310IF Q8 <> J(4) THEN 1340
01320J(4)=Q1/V(J(1)+4)
01330GOTO 1350
01340J(4)=.5*Q2/(C(J(4),1)*V(J(4)+4))+.5*Q3/(C(J(4),2)*V(J(4)+4))
01350N1=V(K5)
01360L1=J(K5)
01370GOSUB 1700
01380IF K5 <> 1 THEN 1450
01390PRINT L$
01400PRINT "HERE ARE THE PERCENTILES OF FOUR INVERSE CHI DISTRIBUTIONS"
01410PRINT "FITTED TO YOUR PERCENTILE SPECIFICATIONS."
01420PRINT
01430REM
01440PRINT" 10TH 25TH 50TH 75TH 90TH"
01450 PRINT USING 1460,K5,S(K5,1),S(K5,2),S(K5,3),S(K5,4),S(K5,5)
01460: ##. #####.## #####.## #####.## #####.## #####.##
01470GOTO 970
01480PRINT
01490PRINT "COMPARE THE PERCENTILES OF THESE DISTRIBUTIONS AND DECIDE"
01500PRINT "WHICH MOST CLOSELY CORRESPONDS TO YOUR PRIOR BELIEFS. YOU"
01510PRINT "CAN EITHER TENTATIVELY ACCEPT THIS DISTRIBUTION OR SPECIFY"
01520PRINT "NEW VALUES FOR THE PERCENTILES."
01530PRINT
01540PRINT "IF YOU WANT ONE OF THE DISTRIBUTION TYPE ITS NUMBER.";
01550PRINT "IF YOU WANT TO RESPECIFY THE PERCENTILES TYPE '0'."
01560GOSUB 9000
01570IF O1=0 THEN 1990
01580GOTO 1620
01590PRINT "REENTER. INPUT MUST BE 0 OR THE NUMBER OF A DISTRIBUTION.";
01600GOSUB 9000
01610GOTO 1570
01620IF O1=1 THEN 1680
01630IF O1=2 THEN 1680
01640IF O1=3 THEN 1680
01650IF O1=4 THEN 1680
01660PRINT
01670GOTO 1590
01680K5=O1
01690GOTO 2440
01700E1=0
01710E2=.5*L1
01720P0=.1
01730GOSUB 1910
01740E1=X3
01750E2=L1*.65
01760P0=.25
01770GOSUB 1910
01780P0=.5
01790E1=X3
01800E2=L1*.87
01810GOSUB 1910
01820E1=X3
01830E2=1.4*L1
01840P0=.75
01850GOSUB 1910
01860E1=X3
01870P0=.9
01880E2=2.2*L1
01890GOSUB 1910
01900RETURN
01910GOSUB 3870
01920S(K5,INT(P0/.25+1.5))=X3
01930RETURN
01940REM **** END OF ROUTINE TO FIND PERCENTILE POINTS *******************
01950PRINT
01960PRINT "RESPECIFY YOUR PERCENTILES."
01970PRINT
01980GOTO 2000
01990PRINT
02000PRINT "INPUT 25TH";
02010GOSUB 9000
02020C1=O1
02030PRINT
02040IF C1>0 THEN 2070
02050PRINT "REENTER. STANDARD DEVIATION MUST BE POSITIVE."
02060GOTO 2000
02070PRINT "INPUT 50TH";
02080GOSUB 9000
02090C2=O1
02100PRINT
02110IF C2>C1 THEN 2150
02120PRINT "REENTER. 50TH MUST BE LARGER THAN 25TH."
02130GOTO 2070
02140GOTO 2000
02150PRINT "INPUT 75TH";
02160GOSUB 9000
02170C3=O1
02180IF C3>C2 THEN 2210
02190PRINT "REENTER. 75TH MUST BE LARGER THAN 50TH."
02200GOTO 2150
02210Q1=C1
02220Q2=C2
02230Q3=C3
02240GOTO 660
02250:FOUR POSSIBLE APPROXIMATE PRIOR DISTRIBUTIONS ARE NOW BEING
02260:COMPUTED FOR YOUR CONSIDERATION.
02270REM *****************************************************************
02280REM ***** ROUTINE TO SOLVE FOR INVERSE CHI PARAMETERS ASSUMING"
02290REM ***** THE VARIANCE AND MEDIAN ARE KNOWN *******
02300REM ******************************************************************
02310K=K**2
02320B=-22*K/3-Q2*Q2
02330A=2*K
02340C=20*K/3+2*Q2*Q2/3
02350F=(-B+SQR(B**2-4*A*C))/(2*A)
02360IF F<500 THEN 2380
02370F=500
02380L=Q2**2*(F-2/3)
02390L=SQR(L)
02400N1=F
02410L1=L
02420RETURN
02430REM ****** SWITCH FROM PERCENTILE TO HYPOTHETICAL SAMPLE SIZE *******
02440M6=S(K5,3)
02450K7=K5
02460N1=V(K5)
02470L1=J(K5)
02480GOTO 2920
02490GOSUB 7300
02500Q5=L1*J1
02510Q6=L1*J2
02520RETURN
02530PRINT "IF YOU DO NOT FEEL THAT THIS HYPOTHETICAL SAMPLE SIZE ( M )"
02540PRINT "REFLECTS YOUR PRIOR INFORMATION ABOUT THE STANDARD DEVIATION"
02550PRINT "YOU CAN SPECIFY A DIFFERENT ONE. A DIFFERENT M WILL NOT AFFECT"
02560PRINT "THE MEDIAN, BUT WILL CHANGE THE HDRS AND OTHER PERCENTILES. A"
02570PRINT "LARGER M WILL SHORTEN THE HDR INTERVALS, AND A SMALLER M WILL"
02580PRINT "LENGTHEN THEM."
02590PRINT
02600REM
02610PRINT "IF YOU WANT TO CHANGE M TYPE THE NEW VALUE (GREATER THAN 3)."
02620PRINT "IF YOU DO NOT WANT TO CHANGE M TYPE '0'."
02630REM
02640GOSUB 9000
02650IF O1=0 THEN 3290
02660IF O1>2000 THEN 2700
02670IF O1>3 THEN 2720
02680PRINT "REENTER. M MUST BE GREATER THAN 3."
02690GOTO 2630
02700PRINT "REENTER. M MUST NOT BE GREATER THAN 2000."
02710GOTO 2630
02720N1=O1-1
02730S9=1
02740GOTO 3450
02750E1=0
02760E2=.65*L1
02770P0=.25
02780GOSUB 3870
02790S(K7,2)=X3
02800IF S9=1 THEN 2860
02810E1=X3
02820E2=.87*L1
02830P0=.5
02840GOSUB 3870
02850S(K7,3)=X3
02860E2=1.4*L1
02870E1=X3
02880P0=.75
02890GOSUB 3870
02900S(K7,4)=X3
02910K5=K7
02920PRINT L$
02930PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE INVERSE CHI"
02940PRINT "DISTRIBUTION YOU ARE NOW CONSIDERING."
02950PRINT
02960 PRINT USING 2980,N1+1
02970M9=L1/SQR(N1+1)
02980: HYPOTHETICAL SAMPLE SIZE(M)####.##
02990IF D9=0 THEN 3060
03000 PRINT USING 3020,N1
03010 PRINT USING 3030,L1
03020: DEGREES OF FREEDOM ###########.##
03030: SCALE PARAMETER ###########.##
03040: MODE ###########.##
03050 PRINT USING 3040,M9
03060 PRINT USING 3070,S(K5,2)
03070: 25TH PERCENTILE ###########.##
03080: 75TH PERCENTILE ###########.##
03090: 50TH (MEDIAN) #######.##
03100 PRINT USING 3090,S(K5,3)
03110 PRINT USING 3080,S(K5,4)
03120: 50% HDR #######.## TO#######.##
03130: 75% HDR #######.## TO#######.##
03140: 95% HDR #######.## TO#######.##
03150G=N1
03160J5=.5
03170GOSUB 2490
03180 PRINT USING 3120,Q5,Q6
03190J5=.75
03200GOSUB 2490
03210 PRINT USING 3130,Q5,Q6
03220J5=.95
03230GOSUB 2490
03240 PRINT USING 3140,Q5,Q6
03250IF D9=1 THEN 3700
03260PRINT
03270IF S9=0 THEN 2530
03280GOTO 3530
03290PRINT
03300IF S8=1 THEN 3350
03310PRINT "YOU CAN CHANGE THE CENTERING OF THE DISTRIBUTION BY"
03320PRINT "SPECIFYING A DIFFERENT MEDIAN. THIS WILL NOT AFFECT"
03330PRINT "THE HYPOTHETICAL SAMPLE SIZE."
03340PRINT
03350PRINT "IF YOU WANT TO CHANGE THE MEDIAN TYPE THE NEW VALUE."
03360PRINT "IF YOU DO NOT TYPE '0'."
03370GOSUB 9000
03380IF O1=0 THEN 3650
03390M6=O1
03400IF O1>0 THEN 3440
03410PRINT "REENTER. INPUT MUST BE '0' OR POSITIVE VALUE FOR MEDIAN.";
03420GOSUB 9000
03430GOTO 3380
03440S8=1
03450L1=1
03460P0=.5
03470E1=0
03480E2=.86
03490GOSUB 3870
03500L1=M6/X3
03510S(K5,3)=L1*X3
03520GOTO 2750
03530PRINT "IF YOU WANT TO CHANGE M TYPE NEW VALUE ELSE '0'.";
03540GOSUB 9000
03550IF O1=0 THEN 3290
03560IF O1>9999 THEN 3600
03570IF O1>3 THEN 3620
03580PRINT "REENTER. M MUST BE GREATER THAN 3."
03590GOTO 3540
03600PRINT "REENTER. M MUST BE LESS THAN 9999."
03610GOTO 3540
03620N1=O1-1
03630GOTO 3450
03640D9=1
03650PRINT L$
03660D9=1
03670PRINT "HERE ARE SOME OF THE CHARACTERISTICS OF THE INVERSE CHI"
03680PRINT "DISTRIBUTION FITTED TO YOUR PRIOR BELIEFS ABOUT SIGMA."
03690GOTO 2950
03700PRINT
03710PRINT "THIS COMPLETES THE SPECIFICATION OF A PRIOR DISTRIBUTION"
03720PRINT "ON SIGMA. IF YOU DO NOT WANT TO FIT A PRIOR DISTRIBTION"
03730PRINT "ON THE MEAN YOU SHOULD RECORD THE PARAMETERS OF YOUR"
03740PRINT "PRIOR DISTRIBUTION ON SIGMA (DEGREES AND SCALE)."
03750PRINT
03760PRINT "IF YOU WANT TO SPECIFY THE PRIOR ON THE MEAN TYPE '1'."
03770PRINT "IF YOU DO NOT TYPE '0'."
03780GOSUB 9000
03790IF O1=0 THEN 3860
03800IF O1=1 THEN 3840
03810PRINT "REENTER. MUST BE '0' OR '1'."
03820GOSUB 9000
03830GOTO 3780
03840SCRATCH#2
03841 PRINT # 2,N1,L1,0,0,S(K5,3),0,0,0,0
03850CHAIN "CMOD5"
03860CHAIN "RSTRT"
03870REM*********************************************************
03880IF E1 <> 0 THEN 3910
03890P8=0
03900GOTO 3940
03910X3=E1
03920GOSUB 4000
03930P8=P
03940X3=E2
03950GOSUB 4000
03960P9=P
03970X3=.5*(E1+E2)
03980GOSUB 4000
03990GOTO 4060
04000X=L1*L1/X3/X3
04010G=N1
04020GOSUB 5775
04030GOSUB 5500
04040P=1-P
04050RETURN
04060IF ABS(P-P0)<.001 THEN 4220
04070IF P<P0 THEN 4190
04080E2=X3
04090P9=P
04100P=(P0-P8)/(P9-P8)
04110IF P<.85 THEN 4140
04120P=.85
04130GOTO 4160
04140IF P>.15 THEN 4160
04150P=.15
04160X3=P*(E2-E1)+E1
04170GOSUB 4000
04180GOTO 4060
04190E1=X3
04200P8=P
04210GOTO 4100
04220RETURN
04230DATA 6011.41,2.1944
04240DATA 5731.39,2.0971
04250DATA 5491.36,2.0152
04260DATA 5281.35,1.9474
04270DATA 5091.33,1.8907
04280DATA 4931.32,1.8392
04290DATA 4781.31,1.7965
04300DATA 4641.29,1.7605
04310DATA 4521.29,1.7292
04320DATA 4411.27,1.7003
04330DATA 4301.27,1.6737
04340DATA 4211.26,1.6487
04350DATA 4121.25,1.627
04360DATA 4031.25,1.6116
04370DATA 3961.24,1.5896
04380DATA 3881.23,1.5736
04390DATA 3801.23,1.553
04400DATA 3721.22,1.5354
04410DATA 3641.22,1.5211
04420DATA 3571.21,1.5041
04430DATA 3501.21,1.4922
04440DATA 3441.2,1.4809
04450DATA 3381.2,1.4675
04460DATA 3321.19,1.4581
04470DATA 3271.19,1.4478
04480DATA 3221.19,1.4364
04490DATA 3171.18,1.4277
04500DATA 3121.18,1.4197
04510DATA 3041.17,1.4029
04520DATA 2961.17,1.3888
04530DATA 2891.16,1.3767
04540DATA 2821.16,1.3651
04550DATA 2701.15,1.3454
04560DATA 2591.14,1.3264
04570DATA 2501.14,1.3125
04580DATA 2411.13,1.297
04590DATA 2341.13,1.2857
04600DATA 2271.12,1.275
04610DATA 2201.12,1.2652
04620DATA 2151.12,1.2577
04630DATA 2091.11,1.2491
04640DATA 2041.11,1.242
04650DATA 2001.11,1.2354
04660DATA 1961.1,1.2279
04670DATA 1911.1,1.2239
04680DATA 1881.1,1.2179
04690DATA 1841.1,1.213
04700DATA 1811.1,1.2081
04710DATA 1781.09,1.2036
04720DATA 1751.09,1.2001
04730DATA 1721.09,1.1963
04740DATA 1691.09,1.1922
04750DATA 1661.09,1.1886
04760DATA 1641.09,1.1856
04770DATA 1621.09,1.1822
04780DATA 1591.08,1.1791
04790DATA 1571.08,1.1765
04800DATA 1551.08,1.1733
04810DATA 1531.08,1.1717
04820DATA 1511.08,1.1681
04830DATA 1491.08,1.1666
04840DATA 1481.08,1.1641
04850DATA 1461.08,1.1618
04860DATA 1441.08,1.1603
04870DATA 1431.07,1.1578
04880DATA 1411.07,1.1557
04890DATA 1401.07,1.1543
04900DATA 1381.07,1.1519
04910DATA 1371.07,1.1504
04920DATA 1351.07,1.1482
04930DATA 1341.07,1.147
04940DATA 1331.06,1.1452
05500REM ********************************************************
05501REM CHI-SQUARE CDF ROUTINE-LOWER TAIL
05502REM INPUT G X
05503REM OUTPUT P
05504REM PRIOR GOSUB 5775
05505REM
05506REM
05510IF G>30 THEN 5795
05515T0=.5*G-1
05520X=X*.5
05525T=81.4983
05530P=5.575E-35*(T+X)**T0
05535T=69.9622
05540P=P+4.0883E-30*(T+X)**T0
05545T=61.0585
05550P=P+2.45182E-26*(T+X)**T0
05555T=53.6086
05560P=P+3.60577E-23*(T+X)**T0
05565T=47.1531
05570P=P+2.01052E-20*(T+X)**T0
05575T=41.4517
05580P=P+5.35019E-18*(T+X)**T0
05585T=36.3584
05590P=P+7.8198E-16*(T+X)**T0
05595T=31.776
05600P=P+6.89418E-14*(T+X)**T0
05605T=27.6359
05610P=P+3.91774E-12*(T+X)**T0
05615T=23.8873
05620P=P+1.50701E-10*(T+X)**T0
05625T=20.4915
05630P=P+4.07286E-09*(T+X)**T0
05635T=17.418
05640P=P+7.96081E-08*(T+X)**T0
05645T=14.6427
05650P=P+1.15132E-06*(T+X)**T0
05655T=12.1461
05660P=P+1.25447E-05*(T+X)**T0
05665T=9.9121
05670P=P+1.04461E-04*(T+X)**T0
05675T=7.92754
05680P=P+6.72163E-04*(T+X)**T0
05685T=6.18154
05690P=P+3.36935E-03*(T+X)**T0
05695T=4.66508
05700P=P+.013226*(T+X)**T0
05705T=3.37077
05710P=P+4.07325E-02*(T+X)**T0
05715T=2.29256
05720P=P+9.81663E-02*(T+X)**T0
05725T=1.4256
05730P=P+.183323*(T+X)**T0
05735T=.766097
05740P=P+.258807*(T+X)**T0
05745T=.311239
05750P=P+.258774*(T+X)**T0
05755T=5.90199E-02
05760P=P+.142812*(T+X)**T0
05761IF X<80 THEN 5765
05762P=1
05763RETURN
05765P=1-(T1*P*EXP(-X))
05770RETURN
05775G9=.5*G
05778IF G>30 THEN 5790
05780GOSUB 5850
05785T1=1/EXP(G0)
05790RETURN
05795X2=((X/G)**(1/3)-(1-2/9/G))/SQR(2/9/G)
05800Y3=X2
05805GOSUB 8000
05810RETURN
05812REM
05813REM END OF CHI-SQUARE CDF ROUTINE
05815REM *********************************************************
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 ****************************************************
07300REM *******************************************************
07301REM INVERSE CHI HIGHEST DENSITY REGION ROUTINE
07302REM INPUTS G J5
07303REM OUTPUTS J1 J2
07304REM
07305IF G>30 THEN 7307
07306GOSUB 5775
07307J8=1
07308GOSUB 7337
07309X=1/(J1*J1)
07310GOSUB 5500
07311J3=P
07312X=1/(J2*J2)
07313GOSUB 5500
07314J3=J3-P
07315IF ABS(J3-J5)>.0001 THEN 7317
07316RETURN
07317IF J3>J5 THEN 7320
07318J8=J8+1
07319GOTO 7308
07320J9=J8-1
07321J0=J8
07322J8=(J0+J9)/2
07323GOSUB 7337
07324X=1/(J1*J1)
07325GOSUB 5500
07326J3=P
07327X=1/(J2*J2)
07328GOSUB 5500
07329J3=J3-P
07330IF ABS(J3-J5)>.0001 THEN 7332
07331RETURN
07332IF J3>J5 THEN 7335
07333J9=J8
07334GOTO 7322
07335J0=J8
07336GOTO 7322
07337J=J8*(1+EXP(-2*J8/(G+1)))/(1-EXP(-2*J8/(G+1)))
07338J1=1/SQR(J+J8)
07339J2=1/SQR(J-J8)
07340RETURN
07341REM
07342REM END OF ROUTINE FOR INVERSE CHI HDR'S
07343REM**************************************************
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)
08021IF X*X/2<80 THEN 8025
08022D=0
08023GOTO 8030
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