Trailing-Edge
-
PDP-10 Archives
-
decus_20tap4_198111
-
decus/20-0113/cmod77.bas
There are 2 other files named cmod77.bas in the archive. Click here to see a list.
00040REM****************************************************************
00050REM CMOD77 CMOD77 CMOD777 CMOD77 CMOD77
00060REM***************************************************************
00070 DIM I(20),Y(30),U(30),S(30),Z(30,1)
00080 DIM J(12),V(202),G(202)
00090 FILES RFILE1,RFILE2,RFILE3
00130RESTORE#1
00131 INPUT# 1,I1,I2,I3
00140SCRATCH#1
00141 PRINT # 1,77,I2,I3
00150Z7=0
00160RESTORE#2
00161 INPUT# 2,V2,E0,L2,V1,L1,U0
00190X=0
00200PRINT L$
00210PRINT " POSTERIOR ANALYSIS - SIMULTANEOUS ESTIMATION OF MEANS"
00220PRINT
00230PRINT "THIS MODULE OBTAINS POSTERIOR DISTRIBUTIONS FOR THE MEANS OF"
00240PRINT "M (MAX=12) DIFFERENT GROUPS UNDER THE ASSUMPTION THAT YOUR "
00250PRINT "PRIOR BELIEFS ABOUT THE MEANS ARE EXCHANGEABLE."
00260PRINT
00270IF V2 <> 0 THEN 1040
00280PRINT "IN ORDER TO DO THE POSTERIOR ANALYSIS YOU MUST PROVIDE"
00290PRINT "CERTAIN INFORMATION ABOUT YOUR PRIOR BELIEFS."
00300PRINT
00310PRINT "IF YOU HAVE THIS INFORMATION TYPE '1', ELSE '0'.";
00320GOSUB 9000
00330IF O1=1 THEN 490
00340IF O1=0 THEN 370
00350PRINT "REENTER. INPUT MUST BE 1 OR 0."
00360GOTO 320
00370PRINT
00380PRINT "PRIOR DISTRIBUTIONS ARE FITTED TO YOUR BELIEFS IN MODULE 1"
00390PRINT "OF THIS MODEL."
00400PRINT
00410PRINT "IF YOU WANT TO FIT PRIORS TYPE '1', ELSE '0'.";
00420GOSUB 9000
00430IF O1=1 THEN 450
00440CHAIN "RSTRT"
00450CHAIN "CMOD75"
00460PRINT
00470PRINT "WHEN YOU ARE READY TO CONTINUE, TYPE '1'.";
00480GOSUB 9000
00490PRINT L$
00500PRINT "INPUT THE PARAMETERS OF YOUR PRIOR DISTRIBUTION ON THE GRAND"
00510Z7=0
00530X=0
00540PRINT "MEAN."
00550PRINT
00560PRINT "DEGREES OF FREEDOM=";
00570GOSUB 9000
00580IF O1 >= 4 THEN 610
00590PRINT "REENTER. MUST BE AT LEAST 4."
00600GOTO 570
00610V2=O1
00620PRINT
00630PRINT "MEAN=";
00640GOSUB 9000
00650U0=O1
00660U2=U0
00670PRINT
00680PRINT "ST. DEV.=";
00690GOSUB 9000
00700IF O1>0 THEN 730
00710PRINT "REENTER. STANDARD DEVAITION MUST BE GREATER THAN 0."
00720GOTO 690
00730E0=O1
00740PRINT
00750PRINT "INPUT THE STANDARD DEVIATION OF YOUR PRIOR ON THE MEAN OF"
00760PRINT "A RANDOMLY SELECTED GROUP.";
00770GOSUB 9000
00780IF O1>E0 THEN 810
00790PRINT "REENTER. MUST BE LARGER THAN ST. DEV. ON GRAND MEAN."
00800GOTO 770
00810REM
00820L2=SQR((V2-2)*(O1*O1-E0*E0))
00830PRINT
00840PRINT "INPUT PARAMETERS OF PRIOR ON RANDOMLY SAMPLED OBSERVATION"
00850PRINT "FROM A GROUP WITH KNOWN MEAN."
00860PRINT
00870PRINT "DEGREES OF FREEDOM=";
00880GOSUB 9000
00910IF O1 >= 4 THEN 940
00920PRINT "REENTER. MUST BE GREATER THAN 4."
00930GOTO 880
00940V1=O1
00950PRINT
00960PRINT "ST. DEV.=";
00970GOSUB 9000
00980IF O1>0 THEN 1010
00990PRINT "REENTER. MUST BE GREATER THAN 0."
01000GOTO 970
01010L1=SQR(O1*O1*(V1-2))
01020E0=L2*L2/((E0*E0)*(V2-2))
01030PRINT L$
01040PRINT "HOW MANY GROUPS ARE THERE IN YOUR DATA SET ";
01050GOSUB 9000
01060M0=O1
01070IF M0>12 THEN 1090
01080IF M0 >= 2 THEN 1110
01090PRINT "REENTER. MINIMUM IS 2. MAXIMUM IS 12."
01100GOTO 1050
01110M0=O1
01120PRINT
01130U2=U0
01150PRINT L$
01160PRINT "ENTER THE SAMPLE DATA:"
01170PRINT " 1. NUMBER OF OBSERVATIONS (N)"
01180PRINT " 2. SAMPLE MEAN (X.)"
01190PRINT " 3. STANDARD DEVIATION (DIVISOR N)"
01200PRINT
01210PRINT "GROUP N, X., ST.DEV."
01220PRINT "--------------------------------------------"
01230M7=0
01240S7=0
01250S2=0
01260N0=0
01270FOR O=1 TO M0
01280REM
01290PRINT" ";O;" ";
01300GOSUB 9100
01310IF O3>0 THEN 1340
01320PRINT "REENTER. STANDARD DEVIATION MUST BE GREATER THAN 0."
01330GOTO 1300
01340IF O1 >= 2 THEN 1370
01350PRINT "REENTER. MUST BE AT LEAST TWO OBSERVATIONS IN EACH GROUP."
01360GOTO 1300
01370S(O)=O1
01380Y(O)=O2
01390J(O)=O3
01400N0=N0+S(O)
01430S2=S2+S(O)*O3*O3
01440NEXT O
01450PRINT L$
01460PRINT "HERE ARE THE SAMPLE DATA YOU ENTERED."
01470PRINT
01480PRINT "GROUP N X. ST.DEV."
01490FOR I=1 TO M0
01500:#### ####### #######.## ########.##
01510PRINT USING 1500,I,S(I),Y(I),J(I)
01520NEXT I
01530PRINT
01540PRINT "IF THE DATA ARE CORRECT TYPE '1', ELSE '0'.";
01550GOSUB 9000
01560IF O1=1 THEN 1610
01570PRINT L$
01580GOTO 1210
01610G=N0+V1+V2
01630PRINT
01640PRINT "TOTAL SAMPLE SIZE =";N0
01650PRINT
01660PRINT "SOME LENGTHY COMPUTATIONS ARE REQUIRED. PLEASE BE PATIENT."
01662GOSUB 6060
01680X0=.00001
01690X8=.99999
01710F8=-1
01720V0=X0
01730GOSUB 3460
01740F7=I0
01750K=0
01760V7=(X8-X0)/50
01770FOR V0=X0 TO X8+.5*V7 STEP V7
01780GOSUB 3460
01790K=K+1
01800V(K)=I0
01810IF I0<F7 THEN 1840
01820F7=I0
01830V6=V0
01840NEXT V0
01850K=0
01860FOR J=1 TO 50
01870IF V(J)+10>F7 THEN 1890
01880NEXT J
01890IF J=1 THEN 1910
01900X0=X0+J*V7-V7
01910FOR J=50 TO 1 STEP -1
01920IF V(J)+10>F7 THEN 1940
01930NEXT J
01940IF J=50 THEN 1960
01950X8=J*V7+X8-49*V7
01960V7=(X8-X0)/200
01970T8=0
01980FOR V0=X0 TO X8+.5*V7 STEP V7
01990GOSUB 3460
02000K=K+1
02010V(K)=I0
02020G(K)=G4
02030T8=T8+EXP(I0-F7)
02040NEXT V0
02050F7=-F7
02060Y9=V7*T8
02090PRINT L$
02100PRINT "TYPE THE NUMBER OF THE OPTION YOU WANT."
02110PRINT " 1. MEANS OF THE POSTERIOR MARGINAL DISTRIBUTIONS ON THE"
02120PRINT " THE GROUP MEANS"
02130PRINT " 2. POSTERIOR MARGINAL PROBABILITIES FOR THE GROUP MEANS"
02140PRINT " 3. POSTERIOR PROBABILITIES FOR LINEAR COMBINATIONS"
02150PRINT " 4. EXIT"
02160GOSUB 9000
02170IF O1=4 THEN 2230
02180IF O1=1 THEN 2870
02190IF O1=3 THEN 2240
02200IF O1=2 THEN 3220
02210PRINT "REENTER. INPUT MUST BE NUMBER OF OPTION."
02220GOTO 2160
02230CHAIN "RSTRT"
02240PRINT L$
02250PRINT "POSTERIOR PROBABILITIES FOR LINEAR COMBINATIONS OF GROUP MEANS"
02260GOSUB 2280
02270GOTO 2410
02280PRINT
02290IF Z7=0 THEN 2340
02300PRINT "GROUP SAMPLE MEAN MEAN OF POSTERIOR MARGINAL"
02302FOR I=1 TO M0
02310:#### #######.## #######.##
02320PRINT USING 2310,I,Y(I),Z(I,1)
02330NEXT I
02332GOTO 2390
02340PRINT "GROUP SAMPLE MEAN"
02350FOR I=1 TO M0
02360:#### #######.##
02370PRINT USING 2360,I,Y(I)
02380NEXT I
02390PRINT
02400RETURN
02410PRINT "INPUT THE NUMBER OF GROUPS IN THE LINEAR COMBINATION (EXIT=0)";
02420GOSUB 9000
02430C7=O1
02440IF O1=0 THEN 2090
02450IF O1>M0 THEN 2470
02460IF O1 >= 2 THEN 2500
02470PRINT "REENTER. NUMBER OF GROUPS MUST BE AT LEAST 2 BUT NOT MORE"
02480PRINT "THE NUMBER OF GROUPS."
02490GOTO 2420
02500PRINT
02510PRINT "INPUT GROUP NUMBER AND COEFFICIENT FOR THAT GROUP."
02520PRINT
02530FOR K0=1 TO C7
02540PRINT "GROUP NUMBER,COEFFICIENT ";
02550GOSUB 9050
02560C(K0)=O2
02570I(K0)=O1
02580IF O1 >= 1 THEN 2610
02590PRINT "REENTER. YOU HAVE ";M0;" GROUPS."
02600GOTO 2540
02610IF O1>M0 THEN 2590
02620NEXT K0
02630T3=0
02640F8=3
02650PRINT L$
02660PRINT "HERE IS THE LINEAR COMBINATION YOU ARE EXAMINING."
02670PRINT
02680PRINT "GROUP COEFFICIENT MEAN"
02690FOR J=1 TO C7
02700:#### #######.## #######.##
02710PRINT USING 2700,J,C(J),Y(I(J))
02720NEXT J
02730PRINT
02740PRINT "MODULE WILL GIVE PROBABILITY LESS THAN SOME VALUE X."
02750PRINT
02760PRINT "INPUT X (EXIT=-7777)";
02770GOSUB 9000
02780T0=O1
02790IF T0=-7777 THEN 2850
02800T4=T0
02810GOSUB 2920
02820 : PROB( L.C. < ######.## ) = ##.##
02830PRINT USING 2820,T0,Y8
02840GOTO 2760
02850PRINT L$
02860GOTO 2260
02870PRINT L$
02880PRINT "HERE ARE THE MEANS OF THE POSTERIOR MARGINAL DISTRIBUTIONS ON THE"
02890PRINT "THE GROUP MEANS."
02900F8=0
02910GOTO 3060
02920REM CONTROL OF INTEGRATION
02930Y8=0
02940K=0
02950FOR V0=X0 TO X8+.5*V7 STEP V7
02951 IF F8<>3 THEN 2960
02952 FOR K2=1 TO C7
02953 K0=I(K2)
02954 U(K0)=(S(K0)*V0)/(1+(S(K0)-1)*V0)
02955 NEXT K2
02956 GOTO 2970
02960U(K2)=(S(K2)*V0)/(1+(S(K2)-1)*V0)
02970K=K+1
02980I0=V(K)
02990G4=G(K)
03000GOSUB 3660
03010Y8=Y8+Y1
03020NEXT V0
03030Y8=Y8*V7/Y9
03040RETURN
03050PRINT
03060PRINT
03070PRINT " GROUP N SAMPLE MEAN MEAN OF POSTERIOR"
03080PRINT "------------------------------------------------------------"
03090FOR I4=1 TO M0
03100K2=I4
03110IF Z7=1 THEN 3150
03120GOSUB 2920
03130:###### ##### ########.## #######.##
03140Z(I4,1)=Y8
03150PRINT USING 3130,K2,S(K2),Y(K2),Z(I4,1)
03160NEXT I4
03170Z7=1
03180PRINT "------------------------------------------------------------"
03190PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
03200GOSUB 9000
03210GOTO 2090
03220PRINT L$
03230PRINT " POSTERIOR PROBABILITIES FOR MEANS"
03235 GOSUB 2280
03240: PROB(MEAN <######.##)= ##.##
03250PRINT "POSTERIOR PROBABILITY THAT GROUP MEAN IS LESS THAN X."
03260PRINT
03270PRINT "INPUT THE GROUP NUMBER (EXIT=0).";
03280GOSUB 9000
03290IF O1 <= 0 THEN 2090
03300IF O1 >= 1 THEN 3330
03310PRINT "REENTER. MUST NUMBER OF GROUP."
03320GOTO 3280
03330IF O1>M0 THEN 3310
03340K2=O1
03350PRINT "INPUT X (EXIT=-7777)";
03360GOSUB 9000
03370IF O1=-7777 THEN 3260
03380T0=O1
03390T4=T0
03400F8=2
03410GOSUB 2920
03420 : PROB ( MEAN < #######.## ) = ##.##
03430PRINT USING 3420,T0,Y8
03440GOTO 3350
03450STOP
03460REM INTEGRAND SETRUP
03470MAT U=ZER
03480W0=0
03490X9=0
03500FOR I=1 TO M0
03510U(I)=(S(I)*V0)/(1+(S(I)-1)*V0)
03520W0=W0+U(I)
03530X9=X9+U(I)*Y(I)
03540NEXT I
03550X9=X9/W0
03560G4=L2*L2+W0*E0/(W0+E0)*(X9-U0)*(X9-U0)
03570P1=0
03580FOR I=1 TO M0
03590G4=G4+U(I)*(Y(I)-X9)*(Y(I)-X9)
03600P1=P1+LOG(U(I))
03610NEXT I
03620G4=G4*(1-V0)/V0+S2+L1*L1
03630I0=(M0+V2+2)*LOG(V0)-(M0+V2-2)*LOG(1-V0)+(N0+V1+V2)*LOG(G4)
03640I0=(P1-LOG(W0+E0)-I0)/2
03650IF F8=-1 THEN 4270
03660IF I0+F7<-80 THEN 4280
03670IF F8=0 THEN 4130
03680IF F8=1 THEN 4210
03690IF F8=2 THEN 4130
03700IF F8=3 THEN 4010
04010C9=0
04020F9=0
04030V9=0
04031 FOR K1=1 TO C7
04032 K0=I(K1)
04033 FOR K2=1 TO C7
04034 K3=I(K2)
04035 IF K1<>K2THEN 4037
04036 V9=V9+C(K1)*C(K2)*(1-U(K0))
04037 V9=V9+(1-U(K0))*(1-U(K3))*C(K1)*C(K2)/(W0+E0)
04038 NEXT K2
04039 NEXT K1
04040FOR K1=1 TO C7
04050K0=I(K1)
04060F9=F9+C(K1)*U(K0)*Y(K0)
04080C9=C9+C(K1)*(1-U(K0))
04090NEXT K1
04100 G4=G4*(V0/(1-V0))*V9
04110F9=F9+C9*(W0*X9+E0*U0)/(W0+E0)
04120GOTO 4180
04130X5=Y(K2)
04140W5=U(K2)
04150F9=W5*X5+(1-W5)*(X9*W0+E0*U0)/(W0+E0)
04160IF F8=0 THEN 4210
04170 G4=G4*V0/(1-V0)*(1-W5+(1-W5)*(1-W5)/(W0+E0))
04180J6=(T4-F9)*SQR(G/G4)
04182J2=ABS(J6)
04190GOSUB 6000
04192IF J6>0 THEN 4200
04193P=1-P
04200F9=P
04210Y1=EXP(I0+F7)
04220IF F8=3 THEN 4260
04230IF F8=2 THEN 4260
04240IF F8=0 THEN 4260
04250RETURN
04260Y1=Y1*F9
04270RETURN
04280Y1=0
04290RETURN
04300REM END OF FUNCTION CALCULATION
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/G5
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 ****************************************************
06000REM****************************************************************
06001REM STUDENT'S T CDF ROUTINE
06002REM INPUT G J2
06003REM OUTPUT P
06004REM PRIOR GOSUB 6060
06005P=0
06006IF J2<.0001 THEN 6034
06007IF G<6 THEN 6011
06008IF J2<6 THEN 6012
06009P=1
06010GOTO 6035
06011IF J2>12 THEN 6009
06012 DIM W(16),O(16)
06013Y3=J2
06014IF G<10 THEN 6019
06015Y3=(G-2/3+.1/G)*SQR(ABS(LOG(1+J2*J2/G))/(G-5/6))
06016GOSUB 8000
06017GOTO 6035
06018REM PEIZER PRATT APPROXIMATION
06019IF G=1 THEN 6098
06020J1=0
06021N=G
06022P=0
06023GOSUB 6084
06024D0=(J2-J1)*.5
06025D1=(J1+J2)*.5
06026FOR I1=1 TO 16
06027D9=D0*O(I1)+D1
06028IF D9=0 THEN 6031
06029IF D9=1 THEN 6031
06030P=P+W(I1)*(EXP(-(N+1)/2*LOG(1+D9*D9/N)))
06031NEXT I1
06032P=P*F0
06033P=P*D0
06034P=P+.5
06035RETURN
06060FOR I1=1 TO 16
06062READ W(I1),O(I1)
06064NEXT I1
06066DATA 2.71525E-02,-.989401
06068DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
06070DATA .124629,-.755404,.149596,-.617876
06072DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
06074DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
06076DATA .149596,.617876,.124629,.755404
06078DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
06080DATA .989401
06082RETURN
06084G9=(N+1)/2
06086GOSUB 5850
06088F0=G0
06090G9=N/2
06092GOSUB 5850
06094F0=EXP(F0-G0)/SQR(3.14159*N)
06096RETURN
06098REM FOLLOWING FOR NU=1
06100P=.5+1/3.14159*ATN(Y3)
06102RETURN
06104REM END OF STUDENT'S T CDF ROUTINE
06106REM*************************************************************
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
09050REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED. 2 INPUTS
09055INPUT O1,O2
09065IF O1=-9999 THEN 9080
09070IF O2=-9999 THEN 9080
09075RETURN
09080CHAIN "RSTRT"
09090REM*************END ROUTINE
09100REM--SUBROUTINE THAT DETERMINES IF RESTART HAS BEEN REQUESTED. 3 INPUTS
09105INPUT O1,O2,O3
09115IF O1=-9999 THEN 9135
09120IF O2=-9999 THEN 9135
09125IF O3=-9999 THEN 9135
09130RETURN
09135CHAIN "RSTRT"
09145REM*************END ROUTINE
09999END