Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-06 - 43,50422/cmod9.bas
There are 2 other files named cmod9.bas in the archive. Click here to see a list.
00011C0=0
00020 REM CMOD9
00100 FILES RFILE1,RFILE2,RFILE3,,,,,,RF9
00140RESTORE#1
00141  INPUT#  1,I1,I2,I3
00150SCRATCH#1
00151  PRINT #  1,9,I2,I3
00160RESTORE#2
00161  INPUT#  2,A8,B8,T1,T2
00170  DIM X(25),N(25),V(25),G(25),H(25),A(100),Y(5)
00180PRINT L$
00190PRINT "       M-GROUP PROPORTIONS POSTERIOR DISTRIBUTIONS"
00200PRINT
00210PRINT "THIS MODULE OBTAINS POSTERIOR DISTRIBUTIONS FOR THE PROBABILITIES"
00220PRINT "OF SUCCESS (PI) IN M DIFFERENT GROUPS FOR THE CASE IN WHICH SAMPLE"
00230PRINT "SIZES ARE NOT THE SAME FOR ALL GROUPS.  A MODEL USING THE VARIANCE-"
00240PRINT "STABILIZING (ARC-SINE) TRANSFORMATION ON THE SAMPLE PROPORTIONS"
00250PRINT "IS EMPLOYED, BUT ALL RESULTS ARE GIVEN IN TERMS OF PI VALUES."
00290IF A8+B8>.1 THEN 600
00300PRINT
00310PRINT "INPUT PARAMETERS (A AND B) OF YOUR FITTED PRIOR DISTRIBUTION"
00320PRINT "FOR A TYPICAL GROUP PROBABILITY OF SUCCESS (PI)."
00330PRINT "A = ";
00340GOSUB 9000
00350IFO1<.1THEN360
00351IFO1<=9999THEN380
00360PRINT "REENTER.  INPUT MUST BE BETWEEN .1 AND 9999."
00370GOTO 330
00380A8=O1
00390PRINT "B = ";
00400 GOSUB 9000
00410IF O1<.1THEN420
00411IFO1<=9999THEN440
00420PRINT "REENTER.  INPUT MUST BE BETWEEN .1 AND 9999."
00430GOTO 390
00440B8=O1
00450PRINT
00460PRINT "NOW INPUT PRIOR NUMBER OF GROUPS (T1 AND T2) ASSOCIATED WITH"
00470PRINT "YOUR KNOWLEDGE CONCERNING THE MEAN AND VARIANCE OF THE PI VALUES."
00480PRINT "T1 (# OF GROUPS FOR MEAN) = ";
00490GOSUB 9000
00500IFO1<1THEN510
00501IFO1<=100THEN530
00510PRINT "REENTER.  INPUT MUST BE BETWEEN 1 AND 100."
00520GOTO 480
00530T1=O1
00540PRINT "T2 (# OF GROUPS FOR VARIANCE) = ";
00550GOSUB 9000
00560IF O1<5THEN570
00561IFO1<=100THEN590
00570PRINT "REENTER.  INPUT MUST BE BETWEEN 5 AND 100."
00580GOTO 540
00590T2=O1
00600PRINT
00601PRINT "AT PRESENT, PRIOR PARAMETER VALUES ARE"
00602PRINT "A =";A8
00603PRINT "B =";B8
00604PRINT "T1=";T1
00605PRINT "T2=";T2
00606PRINT "IF THESE ARE SATISFACTORY, TYPE '1'.  TO INPUT NEW VALUES, TYPE '0'"
00620GOSUB 9000
00621IF O1=C0 THEN 300
00622IF O1=1 THEN 625
00623PRINT "REENTER.  INPUT MUST BE 0 OR 1."
00624GOTO 620
00625SCRATCH#2
00626  PRINT #  2,A8,B8,T1,T2
00630PRINT L$
00640PRINT "INPUT THE ACTUAL NUMBER OF GROUPS (M) IN YOUR ANALYSIS."
00650PRINT "(PROGRAM IS SET UP TO HANDLE VALUES OF M FROM 2 TO 25.)"
00660PRINT "M = ";
00670GOSUB 9000
00680IFO1<2THEN690
00681IFO1<=25THEN710
00690PRINT "REENTER.  INPUT MUST BE BETWEEN 2 AND 25."
00700GOTO 660
00710M0=O1
00720PRINT
00730PRINT "INPUT DATA BY GROUP:  NUMBER OF SUCCESSES FOLLOWED BY A COMMA"
00740PRINT "AND THEN THE NUMBER OF OBSERVATIONS FOR THE GROUP."
00745PRINT "NOTE:  KEEP TRACK OF GROUP NUMBER, SINCE THIS IS THE WAY GROUPS"
00746PRINT "WILL BE REFERRED TO IN SUBSEQUENT ANALYSES."
00747GOSUB 750
00748GOTO 780
00750PRINT
00760PRINT "GROUP  X , N"
00770PRINT "------------"
00775RETURN
00780FOR O=1 TO M0
00781GOSUB 790
00782NEXT O
00783GOTO 1040
00790REM
00800PRINTO
00810GOSUB 9050
00820IFO1<>INT(O1)THEN830
00821IFO2<>INT(O2)THEN830
00822GOTO840
00830GOTO 980
00840IFO1<C0THEN1000
00841IFO1>O2THEN1000
00850IF O2<1 THEN 1020
00860X(O)=O1
00870N(O)=O2
00880V(O)=1/(4*O2+2)
00890N4=1/(O2+1)
00900N5=O1*N4
00910V0=SQR(N5)
00920GOSUB 4010
00930V5=V4
00940V0=SQR(N5+N4)
00950GOSUB 4010
00960G(O)=(V4+V5)/2
00970RETURN
00980PRINT "REENTER BOTH VALUES. THEY MUST BOTH BE INTEGERS."
00990GOTO 800
01000PRINT "REENTER BOTH VALUES.  X MUST BE BETWEEN 0 AND N."
01010GOTO 800
01020PRINT "REENTER BOTH VALUES.  N MUST BE AT LEAST 1."
01040PRINT
01050:YOU HAVE NOW ENTERED DATA FOR ALL ## GROUPS.  IF YOU WISH
01060PRINT  USING 1050,M0
01070PRINT "TO CONTINUE, TYPE '100'.  TO REENTER ALL YOUR DATA, TYPE '0'."
01072PRINT "TO REENTER DATA FOR A SPECIFIC GROUP, TYPE THE GROUP NUMBER."
01080GOSUB 9000
01090IF O1=100 THEN 1140
01100IF O1=C0 THEN 630
01105IFO1<>INT(O1)THEN1110
01106IFO1<1THEN1110
01107IFO1<=M0THEN1121
01110PRINT "REENTER.  INPUT MUST BE 100, 0, OR AN INTEGER BETWEEN 1 AND";M0
01120GOTO 1080
01121O=O1
01122PRINT "OLD VALUES"
01123GOSUB 750
01125:  ## ###,###
01126PRINT  USING 1125,O1,X(O1),N(O1)
01127PRINT
01128PRINT "NEW VALUES"
01129GOSUB 790
01130PRINT
01131GOTO 1070
01140V0=SQR(A8/(A8+B8))
01150GOSUB 4010
01160T3=V4
01180N0=T2-1
01200G1=.25*T1/(T1+1)*(N0-2)/(A8+B8+1)
01250FOR O=1 TO M0
01260H(O)=G(O)
01270NEXT O
01280FOR I9=1 TO 100
01290H1=C0
01300FOR O=1 TO M0
01310H1=H1+H(O)
01320NEXT O
01330H1=H1/M0
01340H2=C0
01350FOR O=1 TO M0
01360H3=H(O)-H1
01370H2=H2+H3*H3
01380NEXT O
01390U=(T1*T3+M0*H1)/(T1+M0)
01400P=(H2+G1+T1*M0*(H1-T3)*(H1-T3)/(T1+M0))/(N0+M0)
01410H5=C0
01420FOR O=1 TO M0
01430H3=(P*G(O)+V(O)*U)/(P+V(O))
01440H4=H3-H(O)
01450H5=H5+H4*H4
01460H(O)=H3
01470NEXT O
01480H6=H5/H2
01490IF H6<1.E-08 THEN 1560
01500NEXT I9
01510PRINT
01520PRINT "THERE ARE CONVERGENCE PROBLEMS WHICH WILL MAKE THE JOINT"
01530PRINT "ESTIMATES OF THE PI VALUES SUSPECT.  THIS WILL NOT AFFECT"
01540PRINT "OTHER RESULTS AND THE ANALYSIS WILL CONTINUE."
01550PRINT
01560REM
01570A9=1
01580GOSUB 1600
01590GOTO 1660
01600PRINT L$
01610PRINT "           POSTERIOR JOINT MODAL ESTIMATES"
01620PRINT
01630PRINT "GROUP         X           N         X/N        ESTIMATE*"
01640PRINT "--------------------------------------------------------"
01650RETURN
01660FOR O=1 TO M0
01670C2=.5/N(O)
01680C3=SIN(H(O))
01690P0=C3*C3
01700:####      #####       #####      ##.###        ##.###
01710PRINT  USING 1700,O,X(O),N(O),X(O)/N(O),P0
01720IF A9<10 THEN 1760
01730GOSUB 1800
01740IF O=M0 THEN 1790
01750GOSUB 1600
01760A9=A9+1
01770NEXT O
01780GOSUB 1800
01790GOTO 1890
01800PRINT "----------------------------------------------------------"
01810PRINT "* THESE ESTIMATES ARE BASED ON THE JOINT MODE OF THE"
01820PRINT "  TRANSFORMED PI VALUES.  (NOTE:  THEY ARE NOT IDENTICAL"
01830PRINT "  TO THE JOINT MODE OF THE PI VALUES THEMSELVES.)"
01840PRINT
01850PRINT "WHEN YOU WANT TO CONTINUE TYPE '1'.";
01860GOSUB 9000
01870A9=C0
01880RETURN
01890REM
01900F8=C0
01910X0=.00001
01920X7=X0
01930GOSUB 4600
01940F7=I0
01950FOR X7=X0+.005 TO 10 STEP .005
01960GOSUB 4600
01970IF F7>I0 THEN 2000
01980F7=I0
01990NEXT X7
02000X9=X7-.0025
02010F6=SGN(4-F7)
02020F7=ABS(4-F7)
02030F7=F6*INT(F7)
02040F8=1
02050PRINT "BEFORE GOING ON, SOME LENGTHY COMPUTATIONS ARE REQUIRED."
02060PRINT "PLEASE BE PATIENT."
02070GOSUB 9900
02080PRINT L$
02090PRINT "YOU MAY NOW OBTAIN INFORMATION ABOUT THE MARGINAL DISTRIBUTION"
02100PRINT "OF PI FOR ANY GROUP.  WHICH GROUP WOULD YOU LIKE TO CONSIDER"
02110PRINT "FIRST";
02120GOSUB 9000
02130IFO1<>INT(O1)THEN2140
02131IFO1<1THEN2140
02132IFO1<=M0THEN2160
02140PRINT "REENTER.  INPUT MUST BE AN INTEGER BETWEEN 1 AND";M0
02150GOTO 2120
02160K2=O1
02170F8=2
02180PRINT
02190PRINT "MORE LENGTHY COMPUTATIONS ARE REQUIRED TO OBTAIN THE MEAN"
02195PRINT "OF THE MARGINAL POSTERIOR DISTRIBUTION FOR THIS GROUP."
02200PRINT
02201PRINT "                                  ESTIMATES"
02202PRINT "GROUP    X     N       X/N     JOINT   MARGINAL*"
02203PRINT "------------------------------------------------"
02204:  ##   ###   ###     ##.###     ##.###    ##.### 
02205 O9=SIN(H(K2))*SIN(H(K2))
02210GOSUB 9900
02220Y2=Y(2)/Y(1)
02230P9=SIN(Y2)*SIN(Y2)
02240PRINT  USING 2204,K2,X(K2),N(K2),X(K2)/N(K2),O9,P9 
02260PRINT "------------------------------------------------"
02270PRINT "* THIS ESTIMATE IS BASED ON THE MEAN OF THE TRANSFORMED"
02280PRINT "  PI VALUE.  (NOTE:  IT IS NOT IDENTICAL TO THE MEAN"
02285PRINT "  OF THE PI VALUE ITSELF.)"
02290PRINT
02300PRINT "IF YOU WOULD LIKE TO SEE SOME APPROXIMATE PERCENTILES,"
02310PRINT "AND/OR APPROXIMATE PROBABILITIES THAT PI EXCEEDS VALUES"
02320PRINT "WHICH YOU SPECIFY, TYPE '1'.  TO SKIP THIS PART OF THE"
02330PRINT "ANALYSIS, TYPE '0'.  (THESE PROBABILITIES AND PERCENTILES"
02340PRINT "ARE BASED ON A NORMAL APPROXIMATION TO THE MARGINAL POST-"
02350PRINT "ERIOR DISTRIBUTION OF THE TRANSFORMED PI VALUES.)"
02360GOSUB 9000
02370IF O1=C0 THEN 3420
02380IF O1=1 THEN 2410
02390PRINT "REENTER.  INPUT MUST BE 0 OR 1."
02400GOTO 2360
02410PRINT "MORE LENGTHY COMPUTATIONS ARE REQUIRED, THIS TIME TO"
02412PRINT "OBTAIN THE STANDARD DEVIATION OF THE MARGINAL DISTRIBUTION"
02413PRINT "FOR THIS GROUP.  WITH THE MEAN AND STANDARD DEVIATION"
02414PRINT "EVALUATED, HOWEVER, THE APPROXIMATE PERCENTILES AND"
02415PRINT "PROBABILITIES CAN BE CALCULATED VERY RAPIDLY."
02416F8=3
02420GOSUB 9900
02424PRINT
02425PRINT "HALFWAY THERE|"
02430F8=4
02440GOSUB 9900
02450Y3=SQR((Y(3)+Y(4))/Y(1))
02470P(1)=-1.2816
02480P(2)=-.6745
02490P(3)=C0
02500P(4)=.6745
02510P(5)=1.2816
02520A9=C0
02530GOSUB 2550
02540 GOTO 2621
02550PRINT L$
02560PRINT "HERE ARE SOME APPROXIMATE PERCENTILES FOR THE POSTERIOR"
02570PRINT "MARGINAL DISTRIBUTION OF PI FOR GROUP";K2
02580PRINT
02590PRINT "  X    N     X/N     10TH     25TH     50TH     75TH     90TH"
02600PRINT "-------------------------------------------------------------"
02610RETURN
02619:##.###
02620:###
02621SCRATCH#9
02622 PRINT#9USING2620,X(K2) 
02623RESTORE#9 
02624INPUT#9,Z$ 
02625PRINTZ$;"  ";
02626SCRATCH#9
02627PRINT#9USING2620,N(K2)
02628RESTORE#9
02629INPUT#9,Z$
02630PRINTZ$;"  "; 
02631SCRATCH#9
02632PRINT#9USING2619,X(K2)/N(K2)
02633RESTORE#9
02634INPUT#9,Z$
02635PRINTZ$;
02640FOR K5=1 TO 5
02650P7=P(K5)*Y3+Y2
02660IF P7>0 THEN 2690
02670P7=C0
02680GOTO 2740
02690IF P7<1.57079 THEN 2720
02700P7=1
02710GOTO 2740
02720P7=SIN(P7)
02730P7=P7*P7
02740SCRATCH#9
02741PRINT#9USING2619,P7
02742RESTORE#9
02743INPUT#9,Z$
02744PRINT"   ";Z$;
02760NEXT K5
02770PRINT
02780PRINT "-------------------------------------------------------------"
02790PRINT "IF YOU WANT THE PROBABILITIES THAT PI EXCEEDS"
02800PRINT "CERTAIN VALUES WHICH YOU SPECIFY (UP TO 5 AT A TIME)"
02810PRINT "TYPE THE NUMBER OF VALUES YOU WANT TO SPECIFY, ELSE '0'."
02860GOSUB 9000
02870IF O1=C0 THEN 3420
02880IFO1<>INT(O1)THEN2890
02881IFO1<C0THEN2890
02882IFO1<=5THEN2910
02890PRINT "REENTER.  INPUT MUST BE AN INTEGER BETWEEN 0 AND 5."
02900GOTO 2860
02910K6=O1
02920GOTO 2950
02930PRINT "REENTER.  VALUE MUST BE BETWEEN 0 AND 1."
02940GOTO 2990
02950PRINT
02960FOR K5=1 TO K6
02970:'CCCC ##
02980Z$="IMAGE"
02981PRINTUSING2970,Z$,K5
02990GOSUB 9000
03000IF O1 <= 0 THEN 2930
03010IF O1 >= 1 THEN 2930
03020P(K5)=O1
03030NEXT K5
03040A9=C0
03050GOSUB 3070
03060GOTO 3200
03070PRINT L$
03080PRINT "HERE IS THE APPROXIMATE PROBABILITY THAT PI FOR GROUP";K2
03090PRINT "EXCEEDS PI(0)."
03100PRINT
03102PRINT "                     VALUES OF PI(0)"
03103PRINT "                     ---------------------------------------------"
03120PRINT"  X    N    X/N"; 
03130FOR K5=1 TO K6
03133 SCRATCH #9
03140:##.###
03150PRINT#9USING 3140,P(K5) 
03152RESTORE#9
03154INPUT#9,Z$
03156PRINT"    ";Z$;
03160NEXT K5
03170PRINT
03180PRINT "------------------------------------------------------------------"
03190RETURN
03199:##.###
03200:###
03201 SCRATCH#9 
03202PRINT#9USING 3200,X(K2) 
03203 RESTORE#9 
03204 INPUT#9,Z$ 
03205 PRINTZ$; 
03206SCRATCH#9
03207 PRINT #9USING3200,N(K2)
03208RESTORE#9
03209INPUT#9,Z$
03210PRINT"  ";Z$;"  ";
03211SCRATCH#9
03212PRINT#9USING3199,X(K2)/N(K2)
03213RESTORE#9
03214INPUT#9,Z$
03215PRINTZ$;
03220FOR K5=1 TO K6
03230V0=SQR(P(K5))
03240GOSUB 4010
03250G6=V4
03260W0=(G6-Y2)/Y3
03270GOSUB 4130
03280U0=1-W4
03290:      ##.###
03295 SCRATCH #9
03320PRINT#9USING 3290,U0 
03322 RESTORE#9
03324 INPUT #9,Z$
03325 PRINT"      ";Z$;
03330NEXT K5
03340PRINT
03350PRINT "------------------------------------------------------------------"
03360PRINT "HOW MANY MORE VALUES DO YOU WANT TO SPECIFY (0-5) ";
03370GOTO 2860
03420PRINT
03430PRINT "IF YOU WISH TO CONSIDER ANOTHER GROUP, TYPE ITS"
03435PRINT "NUMBER.  IF YOU ARE THROUGH WITH THIS ANALYSIS,"
03436PRINT "TYPE 0."
03440GOSUB 9000
03450IF O1=C0 THEN 3520
03460IFO1<>INT(O1)THEN3470
03461IFO1<1THEN3470
03462IFO1<=M0THEN3490
03470PRINT "REENTER.  INPUT MUST BE AN INTEGER BETWEEN 0 AND";M0
03480GOTO 3440
03490PRINT L$
03500GOTO 2160
03520CHAIN "RSTRT"
04000 REM
04010V1=ATN(1)*4
04020IF ABS(V0)<1.E-10 THEN 4050
04030V2=ATN(SQR(1-V0*V0)/V0)+.5*V1*(1-SGN(V0))
04040GOTO 4060
04050V2=V1/2-V0
04060IF ABS(V0)<ABS(V1/2-V2) THEN 4090
04070V3=SGN(V0)*ABS(V0)
04080GOTO 4100
04090V3=SGN(V0)*ABS(V1/2-V2)
04100V4=V3
04110RETURN
04120REM
04130W1=ABS(W0)
04140IF W1 >= 6 THEN 4220
04150W2=.398942*EXP(-W1*W1/2)
04160W3=1/(1+.231642*W1)
04170W4=((((1.33027*W3-1.82126)*W3+1.78148)*W3-.356564)*W3+.319381)
04180W4=1-W4*W3*W2
04190IF W0>0 THEN 4210
04200W4=1-W4
04210IF ABS(W1)<6 THEN 4230
04220W4=.5+SGN(W0)*.5
04230RETURN
04300PRINT
04310PRINT "YOU MAY NOW OBTAIN THE POSTERIOR PROBABILITY THAT PI FOR GROUP"
04320:## IS LESS THAN ANY SPECIFIED PI ZERO BETWEEN .01 AND .99 .
04330PRINT  USING 4320,K2
04340PRINT "(TYPE '0' TO SKIP THIS PART OF THE ANALYSIS.)"
04350PRINT "WHAT VALUE OF PI ZERO WOULD YOU LIKE";
04360GOSUB 9000
04370IF O1=C0 THEN 3420
04380IFO1<.01THEN4390
04381IFO1<=.99THEN4410
04390PRINT "REENTER.  INPUT MUST BE 0 OR A VALUE BETWEEN .01 AND .99 ."
04400GOTO 4350
04410V0=SQR(O1)
04420GOSUB 4000
04430G0=V4
04440F8=5
04450GOSUB 9900
04460PRINT
04470PRINT
04480PRINT "COMPUTING CUMULATIVE POSTERIOR PROBABILITY"
04490PRINT
04500:PROB(PI<.##) = .##
04510PRINT  USING 4500,P9,Y(5)/Y(1)
04520PRINT
04530PRINT "IF OTHER PROB. IS TO BE EVALUATED, ENTER 1; OTHERWISE,ENTER 0"
04540GOSUB 9000
04550GOTO 3420
04560GOTO 4300
04600REM
04610W6=C0
04620G9=C0
04630FOR I=1 TO M0
04640V9=X7/(X7+V(I))
04650W6=W6+V9
04660G9=G9+G(I)*V9
04670NEXT I
04680G9=G9/W6
04690G4=G1+(T1*W6/(T1+W6))*(G9-T3)*(G9-T3)
04700P1=C0
04710FOR I=1 TO M0
04720G4=G4+X7/(X7+V(I))*(G(I)-G9)*(G(I)-G9)
04730P1=P1+LOG(X7+V(I))
04740NEXT I
04750G4=G4/X7
04760G4=G4+LOG(T1+W6)+(N0+2)*LOG(X7)
04770I0=-(G4+P1)/2
04780IF F8=C0 THEN 4970
04790F9=1
04800IF F8=1 THEN 4950
04810W5=X7/(X7+V(K2))
04820G5=G(K2)
04830F9=W5*G5+(1-W5)*(T1*T3+W6*G9)/(T1+W6)
04840IF F8=2 THEN 4950
04850IF F8=4 THEN 4940
04860V6=X7*(1-W5)*(1+(1-W5)/(T1+W6))
04870IF F8=5 THEN 4900
04880F9=V6
04890GOTO 4950
04900W0=(G0-F9)/SQR(V6)
04910GOSUB 4120
04920F9=W4
04930GOTO 4950
04940F9=(F9-Y2)*(F9-Y2)
04950IF I0+F7+LOG(F9)<-80 THEN 4980
04960Y1=EXP(I0+F7+LOG(F9))
04970RETURN
04980Y1=C0
04990RETURN
09000REM
09005INPUT O1
09015IF O1=-9999 THEN 9025
09020RETURN
09025CHAIN "RSTRT"
09050REM
09055INPUT O1,O2
09065IF O1=-9999 THEN 9080
09070IF O2=-9999 THEN 9080
09075RETURN
09080CHAIN "RSTRT"
09600REM   ROMBERG INTEGRATION ROUTINE
09615E0=.00005
09625X7=X0
09630GOSUB 9885
09635A(1)=Y1
09640X7=X8
09645GOSUB 9885
09650A(1)=(Y1+A(1))*.5
09655H0=X8-X0
09660IF ABS(H0)<.000001 THEN 9875
09665H1=H0
09670E1=E0/ABS(H0)
09675D3=C0
09680P0=1
09685J0=1
09690FOR I9=2 TO 100
09695Y0=A(1)
09700D2=D3
09705H2=H1
09710H1=.5*H1
09715P0=.5*P0
09720X7=X0+H1
09725S0=C0
09730FOR J=1 TO J0
09735GOSUB 9885
09740S0=S0+Y1
09745X7=X7+H2
09750NEXT J
09755A(I9)=.5*A(I9-1)+P0*S0
09760Q0=1
09765J9=I9-1
09770FOR J=1 TO J9
09775I8=I9-J
09780Q0=4*Q0
09785A(I8)=A(I8+1)+(A(I8+1)-A(I8))/(Q0-1)
09790NEXT J
09795D3=ABS(Y0-A(1))
09800IF (I9-5)<C0 THEN 9815
09805IF (D3-E1) <= 0 THEN 9875
09810IF (D3-D2) >= 0 THEN 9835
09815J0=J0+J0
09820NEXT I9
09825PRINT "100 ITERATIONS WITHOUT REACHING CONVERGENCE."
09830GOTO 9840
09835PRINT "CONVERGENCE PROBLEMS, PROBABLY DUE TO ROUNDING ERROR"
09840PRINT "CRITERION = ";E0,"MINIMUM ATTAINED =";D2*H0
09845PRINT "ANALYSIS WILL CONTINUE.  IF, HOWEVER, DISCREPANCY BETWEEN"
09850PRINT "THE ABOVE TWO VALUES IS GREAT, SUBSEQUENT RESULTS MAY BE IN"
09855PRINT "ERROR."
09860IF I9>99.5 THEN 9875
09865Y1=Y0*H0
09870RETURN
09875Y1=H0*A(1)
09880RETURN
09885GOSUB 4600
09890RETURN
09900REM
09905X0=.00001
09910I5=.5
09915Y(F8)=C0
09920X8=X0+X9*I5
09925GOSUB 9600
09930Y(F8)=Y(F8)+Y1
09935IF Y(F8) <= 0 THEN 790
09940IF Y1/Y(F8) <= .00005 THEN 9960
09945X0=X8
09950I5=2*I5
09955GOTO 9920
09960RETURN
09999END