Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmoda.bas
There are 2 other files named cmoda.bas in the archive. Click here to see a list.
00020REM *******************************************************
00030REM CMODA CMODA CMODA CMODA CMODA CMODA
00040REM *******************************************************
00050REM
00060REM POSTERIOR ANALYSIS FOR M-GROUP PROPORTIONS (EQUAL N)
00070REM USING ARC-SINE TRANSFORMATION
00080REM
00090REM******************************************************
00100 FILES RFILE1,RFILE2,RFILE3,,,,,,RF9
00140RESTORE#1
00141 INPUT# 1,I1,I2,I3
00150SCRATCH#1
00151 PRINT # 1,10,I2,I3
00160RESTORE#2
00161 INPUT# 2,A8,B8,T1,T2
00170 DIM F(26),G(26),H(26),S(26)
00180 DIM X(26),J(10)
00190PRINT L$
00200PRINT " M-GROUP PROPORTIONS POSTERIOR DISTRIBUTIONS"
00210PRINT
00220PRINT "THIS MODULE OBTAINS POSTERIOR DISTRIBUTIONS FOR THE PROBABILITIES"
00230PRINT "OF SUCCESS (PI) IN M DIFFERENT GROUPS FOR THE CASE IN WHICH SAMPLE"
00240PRINT "SIZES FOR ALL GROUPS ARE THE SAME. A MODEL USING THE VARIANCE-"
00250PRINT "STABILIZING (ARC-SINE) TRANSFORMATION ON THE SAMPLE PROPORTIONS"
00260PRINT "IS EMPLOYED, BUT ALL RESULTS ARE GIVEN IN TERMS OF PI VALUES."
00270REM *******************************************************************
00280REM CHECK TO SEE IF PRIOR PARAMETERS HAVE BEEN SPECIFIED. IF NOT,
00290REM GET VALUES FOR A8,B8,T1,T2
00300IF A8+B8>.1 THEN 600
00310PRINT
00320PRINT "INPUT PARAMETERS (A AND B) OF YOUR FITTED PRIOR DISTRIBUTION"
00330PRINT "FOR A TYPICAL GROUP PROBABILITY OF SUCCESS (PI)."
00340PRINT "A = ";
00350GOSUB 9000
00360 IF O1<.1 THEN 370
00362IF O1<=9999 THEN 390
00370PRINT "REENTER. INPUT MUST BE BETWEEN .1 AND 9999."
00380GOTO 340
00390A8=O1
00400PRINT "B = ";
00410GOSUB 9000
00420 IF O1<.1 THEN 430
00422 IF O1<=9999 THEN 450
00430PRINT "REENTER. INPUT MUST BE BETWEEN .1 AND 9999."
00440GOTO 400
00450B8=O1
00460PRINT
00470PRINT "NOW INPUT PRIOR NUMBER OF GROUPS (T1 AND T2) ASSOCIATED WITH"
00480PRINT "YOUR KNOWLEDGE CONCERNING THE MEAN AND VARIANCE OF THE PI VALUES."
00490PRINT "T1 (# OF GROUPS FOR MEAN) = ";
00500GOSUB 9000
00510 IF O1<1 THEN 520
00515 IF O1<=100 THEN 540
00520PRINT "REENTER. INPUT MUST BE BETWEEN 1 AND 100"
00530GOTO 490
00540T1=O1
00550PRINT "T2 (# OF GROUPS FOR VARIANCE) = ";
00560GOSUB 9000
00570 IF O1<5 THEN 580
00575 IF O1<=100 THEN 595
00580PRINT "REENTER. INPUT MUST BE BETWEEN 5 AND 100."
00590GOTO 550
00595T2=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=0 THEN 310
00622IF O1=1 THEN 625
00623PRINT "REENTER. INPUT MUST BE 0 OR 1."
00624GOTO 620
00625SCRATCH#2
00626 PRINT # 2,A8,B8,T1,T2
00640PRINT L$
00650PRINT "INPUT THE ACTUAL NUMBER OF GROUPS (M) IN YOUR ANALYSIS."
00670PRINT "M = ";
00680GOSUB 9000
00690IF O1<>INT(O1) THEN 700
00695 IF O1>=2 THEN 720
00700PRINT "REENTER. INPUT MUST BE AN INTEGER, 2 OR GREATER."
00710GOTO 670
00720M0=O1
00730PRINT
00740PRINT "NOW INPUT NUMBER OF OBSERVATIONS (N) PER GROUP."
00760PRINT "N = ";
00770GOSUB 9000
00780IF O1<>INT(O1) THEN 790
00785 IF O1>=6 THEN 810
00790PRINT "REENTER. INPUT MUST BE AN INTEGER, 6 OR GREATER."
00800GOTO 760
00810N3=O1
00820PRINT
00830PRINT "NOW ENTER YOUR SAMPLE DATA AS A FREQUENCY DISTRIBUTION OF THE"
00840:NUMBER OF GROUPS (OUT OF A TOTAL OF###) HAVING X SUCCESSES
00850:FOR EACH VALUE OF X FROM 0 TO### WHICH HAS BEEN OBSERVED
00860PRINT USING 840,M0
00870PRINT USING 850,N3
00880PRINT "FOR AT LEAST ONE GROUP. (I.E. YOU MAY SKIP ALL X VALUES WITH"
00890PRINT "ZERO ASSOCIATED FREQUENCIES.)"
00900PRINT "NOTE: MODULE WILL ACCEPT UP TO 25 DISTINCT FREQUENCIES"
00905PRINT
00910PRINT "INPUT NUMBER OF SUCCESSES, FOLLOWED BY FREQUENCY (NUMBER OF"
00920PRINT "GROUPS). THUS, IF 3 GROUPS HAD 0 SUCCESSES, YOU WOULD TYPE"
00930PRINT " 0,3 ."
00940PRINT
00950S=0
00960N1=0
00970N1=N1+1
00980GOSUB 9050
00990 IF O1<>INT(O1)THEN1000
00995 IF O1<0 THEN 1000
00996 IF O1<=N3 THEN 1030
01000PRINT "REENTER BOTH VALUES. THE FIRST INPUT MUST BE AN INTEGER BETWEEN"
01010PRINT "0 AND ";N3
01020GOTO 980
01030IF O2<>INT(O2) THEN 1040
01033IF O2<0 THEN 1040
01036 IF O2<=M0 THEN 1070
01040PRINT "REENTER BOTH VALUES. THE SECOND INPUT MUST BE AN INTEGER BETWEEN"
01050PRINT "0 AND ";M0
01060GOTO 980
01070X(N1)=O1
01080F(N1)=O2
01090S=S+O2
01100IF S<M0-.1 THEN 970
01110IF S>M0+.1 THEN 1230
01120PRINT
01130PRINT "THE SUM OF YOUR FREQUENCIES NOW EQUALS THE TOTAL NUMBER"
01140:OF GROUPS (###). IF YOU ARE READY TO GO ON TYPE '1'.
01150PRINT"OF GROUPS (";M0;"). IF YOU ARE READY TO CONTINUE TYPE '1'."
01160PRINT "IF YOU WISH TO REENTER YOUR DATA, TYPE '0'."
01170PRINT
01180GOSUB 9000
01190IF O1=1 THEN 1310
01200IF O1=0 THEN 640
01210PRINT "REENTER. VALUE MUST BE EITHER 0 OR 1.";
01220GOTO 1180
01230PRINT
01240PRINT "THE SUM OF YOUR FREQUENCIES NOW EXCEEDS THE TOTAL"
01250:NUMBER OF GROUPS (###) YOU SPECIFIED FOR THIS ANALYSIS.
01260PRINT USING 1250,M0
01270PRINT "WHEN YOU ARE READY TO REENTER YOUR DATA, PLEASE TYPE '1'.";
01280GOSUB 9000
01290GOTO 640
01300REM ***************************
01310V0=SQR(A8/(A8+B8))
01320GOSUB 3930
01330T3=V4
01340REM T3 IS THETA
01350N0=T2-1
01360REM N0 IS NU
01370G1=.25*T1/(T1+1)*(N0-2)/(A8+B8+1)
01380REM G1 IS LAMBDA
01390REM T1 IS ETA
01400REM ***************************
01410P4=1/(4*N3+2)
01420REM P4 IS V, THE VARIANCE OF THE TRANSFORMED PROPORTIONS.
01430N4=1/(N3+1)
01460G0=0
01470G2=0
01480FOR O=1 TO N1
01490N5=X(O)*N4
01500V0=SQR(N5)
01510GOSUB 3930
01520G(O)=V4
01530V0=SQR(N5+N4)
01540GOSUB 3930
01550 G(O)=(G(O)+V4)/2
01560G0=G0+G(O)*F(O)
01570NEXT O
01580G0=G0/M0
01590FOR O=1 TO N1
01600G2=G2+(G(O)-G0)*(G(O)-G0)*F(O)
01610NEXT O
01620G3=G2/M0
01630M2=N0+M0
01640M3=T1+M0
01650M4=(G0-T3)*(G0-T3)
01660M5=M4*T1*M0/M3
01670M6=(N0+M0-2)/2
01680M7=2*P4
01690M8=M4*T1*M0
01695PRINT
01696PRINT "BEFORE GOING ON, SOME LENGTHY COMPUTATIONS ARE REQUIRED."
01697PRINT "PLEASE BE PATIENT."
01700GOSUB 4160
01710Q1=-(N0+2)/2
01720Q2=-.5
01730Q3=1
01740F8=-1
01750X7=.00001
01760GOSUB 3820
01770F7=Q9
01780FOR X7=.009 TO .999 STEP .01
01790GOSUB 3820
01800IF F7>Q9 THEN 1830
01810F7=Q9
01820NEXT X7
01830F6=SGN(4-F7)
01840F7=ABS(4-F7)
01850F7=F6*INT(F7)
01860F8=1
01870Q1=Q1-1
01880FOR O=1 TO 6
01890IF O>3 THEN 1920
01900Q1=Q1+1
01910GOTO 1980
01920IF O=5 THEN 1970
01930Q2=Q2-1
01940Q3=Q3*M3
01950IF O=6 THEN 1900
01960GOTO 1980
01970Q1=Q1-1
01980M9=0
01990X0=.00001
02000FOR X8=.25 TO .75 STEP .25
02010GOSUB 9600
02020M9=M9+Y1
02030X0=X8
02040NEXT X8
02050X8=.99999
02060GOSUB 9600
02070J(O)=M9+Y1
02080IF O=1 THEN 2100
02090J(O)=J(O)/J(1)
02100NEXT O
02110REM J(2) IS EXPECTED VALUE OF RO
02120REM J(5) IS EXPECTED VALUE OF TAU
02130J(3)=J(3)-J(2)*J(2)
02140REM J(3) IS VARIANCE OF RO
02150J(4)=J(4)-J(2)*J(5)
02160REM J(4) IS COVARIANCE BETWEEN RO AND TAU
02170J(6)=J(6)-J(5)*J(5)
02180REM J(6) IS VARIANCE OF TAU
02190H3=J(2)
02200H4=(J(5)-J(2))*G0+(1-J(5))*T3
02210H5=J(2)*P4+(J(5)-J(2))*P4/M0+(G0-T3)*(G0-T3)*J(6)
02220H6=J(3)
02230H7=2*(G0-T3)*J(4)
02240FOR O=1 TO N1
02250H(O)=H3*G(O)+H4
02260H8=G(O)-G0
02270S(O)=H5+H7*H8+H6*H8*H8
02280S(O)=SQR(S(O))
02290NEXT O
02300GOSUB 2320
02310GOTO 2420
02320PRINT L$
02330PRINT " JOINT AND MARGINAL POINT ESTIMATES"
02360PRINT
02370: 'C X FREQUENCY X/N PI(JOINT)* PI(MARG.)**
02380PRINT USING 2370," "
02400PRINT" ----------------------------------------------------------"
02410RETURN
02420C2=1/(2*N3)
02430A9=0
02440FOR O=1 TO N1
02450P2=X(O)/N3
02460G5=T3+R0*(G(O)-G0)+T0*(G0-T3)
02470P0=SIN(G5)*SIN(G5)
02475REM COULD ALSO USE (1+C2)*SIN(G5)**2-C2/2
02480P1=SIN(H(O))*SIN(H(O))
02485REM DITTO
02490: ### ### ##.### ##.### ##.###
02500PRINT USING 2490,X(O),F(O),P2,P0,P1
02510A9=A9+1
02520IF A9<10 THEN 2590
02532GOSUB 2535
02534GOTO 2560
02535PRINT" -------------------------------------------------------"
02537PRINT " * JOINT ESTIMATES ARE BASED ON THE JOINT MODE OF THE"
02538PRINT " TRANSFORMED PI VALUES. THEY ARE NOT IDENTICAL TO"
02540PRINT " THE JOINT MODE OF THE PI VALUES THEMSELVES."
02541PRINT
02542PRINT " ** MARGINAL ESTIMATES ARE BASED ON THE MEANS OF THE"
02544PRINT " TRANSFORMED PI VALUES. THEY ARE NOT IDENTICAL TO"
02546PRINT " THE MEANS OF THE PI VALUES THEMSELVES."
02547PRINT
02550PRINT " WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
02552GOSUB 9000
02554RETURN
02560IF A9=N1 THEN 2630
02570GOSUB 2320
02580A9=0
02590NEXT O
02605GOSUB 2535
02630PRINT L$
02640V(1)=-1.2816
02650V(2)=-.6745
02660V(3)=0
02670V(4)=.6745
02680V(5)=1.2816
02690A9=0
02700GOSUB 2720
02710GOTO 2800
02720PRINT L$
02730PRINT " APPROXIMATE PERCENTILES FOR THE PI VALUES"
02740PRINT " (BASED ON A NORMAL APPROXIMATION TO THE MARGINAL POSTERIOR"
02750PRINT " DISTRIBUTIONS OF THE TRANSFORMED PI VALUES)"
02760PRINT
02770PRINT " X FREQ. X/N 10TH 25TH 50TH 75TH 90TH"
02780PRINT "-------------------------------------------------------------"
02790RETURN
02800FOR O=1 TO N1
02810:"### ### ##.###"
02815 SCRATCH#9
02820PRINT#9USING 2810,X(O),F(O),X(O)/N3
02822RESTORE#9
02824INPUT#9,D$
02825PRINTD$;
02830FOR K5=1 TO 5
02840P7=V(K5)*S(O)+H(O)
02841IF P7>0 THEN 2844
02842P7=0
02843GOTO 2900
02844IF P7<1.57079 THEN 2850
02845P7=1
02846GOTO 2900
02850P7=SIN(P7)
02860P7=P7*P7
02870:" ##.###"
02872SCRATCH#9
02900PRINT#9USING 2870,P7
02905RESTORE#9
02910INPUT#9,D$
02912PRINTD$;
02915 NEXT K5
02920PRINT
02930IF A9<14 THEN 3000
02940IF A9+1=N1 THEN 3020
02950PRINT "-------------------------------------------------------------"
02960PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
02970GOSUB 9000
02980A9=0
02990GOSUB 2720
03000A9=A9+1
03010NEXT O
03020PRINT "-------------------------------------------------------------"
03030PRINT "IF YOU WANT THE PROBABILITIES THAT PI EXCEEDS CERTAIN VALUES"
03040PRINT "WHICH YOU SPECIFY (UP TO 5 AT A TIME), TYPE THE NUMBER OF"
03045PRINT "VALUES YOU WISH TO SPECIFY, ELSE '0' "
03050GOSUB 9000
03060IF O1=0 THEN 3810
03070 IF O1<>INT(O1) THEN 3073
03071IF O1<1 THEN 3073
03072 IF O1<=5 THEN 3080
03073PRINT "REENTER. INPUT MUST BE AN INTEGER BETWEEN 0 AND 5."
03076GOTO 3050
03080PRINT L$
03150GOTO 3180
03160PRINT "REENTER. VALUE MUST BE BETWEEN 0 AND 1."
03170GOTO 3230
03180K6=O1
03190PRINT
03200FOR K5=1 TO K6
03210:"VALUE ## "
03215SCRATCH#9
03220PRINT#9USING 3210,K5
03222RESTORE#9
03223INPUT#9,D$
03225PRINTD$;
03230GOSUB 9000
03240 IF O1<0 THEN 3160
03241IF O1>1 THEN 3160
03260V(K5)=O1
03270NEXT K5
03280A9=0
03290GOSUB 3310
03300GOTO 3450
03310PRINT L$
03320PRINT " APPROXIMATE PROBABILITY THAT PI EXCEEDS PI(0)"
03330PRINT " (BASED ON A NORMAL APPROXIMATION TO THE MARGINAL POSTERIOR"
03340PRINT " DISTRIBUTIONS OF THE TRANSFORMED PI VALUES)"
03350PRINT
03358PRINT " VALUES OF PI(0)"
03359PRINT " ---------------------------------------------"
03360 PRINT" X FREQ. X/N";
03380FOR K5=1 TO K6
03390:" ##.###"
03395SCRATCH#9
03400PRINT#9USING 3390,V(K5)
03405 RESTORE#9
03406INPUT#9,D$
03407PRINTD$;
03410NEXT K5
03420PRINT
03425PRINT"-----------------------------------------------------------"
03440RETURN
03450FOR K7=1 TO N1
03460:"### ### ##.###"
03465SCRATCH#9
03470PRINT#9USING 3460,X(K7),F(K7),X(K7)/N3
03472RESTORE#9
03473INPUT#9,D$
03474PRINTD$;
03480FOR K5=1 TO K6
03490V0=SQR(V(K5))
03500GOSUB 3930
03510G6=V4
03520Z0=(G6-H(K7))/S(K7)
03530W0=Z0
03540GOSUB 4050
03550U0=1-W4
03560:" ##.##"
03565SCRATCH#9
03590PRINT #9USING 3560,U0
03592RESTORE#9
03594INPUT#9,D$
03596PRINTD$;
03600NEXT K5
03610PRINT
03620A9=A9+1
03630IF A9<10 THEN 3700
03640IF A9=N1 THEN 3710
03650PRINT"---------------------------------------------------------------"
03660PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
03670GOSUB 9000
03680A9=0
03690GOSUB 3310
03700NEXT K7
03710PRINT"------------------------------------------------------------- "
03720PRINT "HOW MANY MORE VALUES DO YOU WISH TO SPECIFY (0-5)?"
03725PRINT "(IF YOU ARE THROUGH WITH THIS ANALYSIS, TYPE '0'.)"
03730GOTO 3050
03810CHAIN "RSTRT"
03820REM OBTAINS FUNCTION TO BE INTEGRATED INPUT X7 OUTPUT Y1
03830Q4=T1+M0*X7
03840Q5=1-X7
03850Q9=Q1*LOG(X7)+Q2*LOG(Q4)+M6*LOG(Q5)-Q5/M7*(G1/X7+G2+M8/Q4)
03860IF F8<0 THEN 3890
03870IF Q9+F7<-80 THEN 3900
03880Y1=EXP(Q9+F7)*Q3
03890RETURN
03900Y1=0
03910RETURN
03920REM FINDS ARC-SINE OF V0 AND RETURNS IT AS V4
03930V1=ATN(1)*4
03940IF ABS(V0)<1.E-10 THEN 3970
03950V2=ATN(SQR(1-V0*V0)/V0)+.5*V1*(1-SGN(V0))
03960GOTO 3980
03970V2=V1/2-V0
03980IF ABS(V0)<ABS(V1/2-V2) THEN 4010
03990V3=SGN(V0)*ABS(V0)
04000GOTO 4020
04010V3=SGN(V0)*ABS(V1/2-V2)
04020V4=V3
04030RETURN
04040REM FINDS CUMULATIVE NORMAL PROBABILITIES. INPUT W0, OUTPUT W4.
04050W1=ABS(W0)
04060IF W1 >= 6 THEN 4140
04070W2=.398942*EXP(-W1*W1/2)
04080W3=1/(1+.231642*W1)
04090W4=((((1.33027*W3-1.82126)*W3+1.78148)*W3-.356564)*W3+.319381)
04100W4=1-W4*W3*W2
04110IF W0>0 THEN 4130
04120W4=1-W4
04130IF ABS(W1)<6 THEN 4150
04140W4=.5+SGN(W0)*.5
04150RETURN
04160REM FINDS RO TO GIVE JOINT POSTERIOR MODAL VALUES FOR GAMMAS
04170R0=-1
04180R4=1
04190T4=1
04200FOR C4=1 TO 100
04210F4=(R4*R4*G2+G1+T4*T4*M5)/M2
04220R5=F4/(F4+P4)
04230T5=M3*R5/(T1+M0*R5)
04240IF ABS(R4-R5)<.00001 THEN 4310
04250R4=R5
04260T4=T5
04270NEXT C4
04280R0=R5
04290T0=T5
04300GOTO 4380
04310IF R0 <> -1 THEN 4370
04320R0=R5
04330T0=T5
04340R4=0
04350T4=0
04360GOTO 4200
04370IF ABS(R0-R5)<.001 THEN 4420
04380PRINT
04390PRINT "THERE ARE CONVERGENCE PROBLEMS WHICH WILL MAKE THE JOINT"
04400PRINT "ESTIMATES OF THE PI VALUES SUSPECT. THIS WILL NOT AFFECT"
04410PRINT "OTHER RESULTS AND THE ANALYSIS WILL CONTINUE."
04420RETURN
09000REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED.
09005INPUT O1
09010IF O1=-9999 THEN 9020
09015RETURN
09020CHAIN "RSTRT"
09025REM*************END ROUTINE
09050REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED. 2 INPUTS
09055INPUT O1,O2
09060IF O1=-9999 THEN 9075
09065IF O2=-9999 THEN 9075
09070RETURN
09075CHAIN "RSTRT"
09080REM*************END ROUTINE
09600REM ROMBERG INTEGRATION ROUTINE
09605REM REQUIRES FUNCTION F(X) DEFINED AT 9800
09610REM MUST THEN GOSUB TO 9800
09615REM MAY CHANGE TO DEF IF YOU CAN FIT FUNCTION ON ONE LINE
09620REM X0=LOWER BOUND X8=UPPER BOUND
09625 DIM A(100)
09630E0=.002
09635REM EPSILON VALUE IN E0
09640X7=X0
09645GOSUB 9895
09650A(1)=Y1
09655X7=X8
09660GOSUB 9895
09665A(1)=(Y1+A(1))*.5
09670H0=X8-X0
09675IF ABS(H0)<.000001 THEN 9880
09680H1=H0
09685E1=E0/ABS(H0)
09690D3=0
09695P0=1
09700J0=1
09705FOR I9=2 TO 100
09710Y0=A(1)
09715D2=D3
09720H2=H1
09725H1=.5*H1
09730P0=.5*P0
09735X7=X0+H1
09740S0=0
09745FOR J1=1 TO J0
09750GOSUB 9895
09755S0=S0+Y1
09760X7=X7+H2
09765NEXT J1
09770A(I9)=.5*A(I9-1)+P0*S0
09775Q0=1
09780J9=I9-1
09785FOR J1=1 TO J9
09790I8=I9-J1
09795Q0=4*Q0
09800A(I8)=A(I8+1)+(A(I8+1)-A(I8))/(Q0-1)
09805NEXT J1
09810D3=ABS(Y0-A(1))
09815IF (I9-5)<0 THEN 9830
09820IF (D3-E1) <= 0 THEN 9880
09825IF (D3-D2) >= 0 THEN 9845
09830J0=J0+J0
09835NEXT I9
09840GOTO 9880
09845PRINT "CONVERGENCE PROBLEMS, PROBABLY DUE TO ROUNDING ERROR"
09850PRINT "CRITERION = ";E0,"MINIMUM ATTAINED =";D2*H0
09855PRINT "ANALYSIS WILL CONTINUE. IF, HOWEVER, DISCREPANCY BETWEEN"
09860PRINT "THE ABOVE TWO VALUES IS GREAT, SUBSEQUENT RESULTS MAY BE IN"
09865PRINT "ERROR."
09870Y1=Y0*H0
09875GOTO 9885
09880Y1=H0*A(1)
09885REM AFTER I9 ITERATIONS, INTEGRAL = Y1
09890RETURN
09895GOSUB 3820
09900RETURN
09999 END