Trailing-Edge
-
PDP-10 Archives
-
decuslib10-06
-
43,50422/cmodg.bas
There are 2 other files named cmodg.bas in the archive. Click here to see a list.
00020REM*************************************************************
00030REM CMODG CMODG CMODG CMODG CMODG CMODG
00040REM***********************************************************
00050REM
00060REM PRIOR FOR MULTINOMIAL DIRCHLET
00070REM
00080REM*******************************************************
00090FILES RFILE1,RFILE2,RFILE3
00130RESTORE#1
00131 INPUT# 1,I1,I2,I3
00140 SCRATCH#1
00141 PRINT # 1,16,I2,I3
00150DIM Q(10)
00160DIM R(10)
00170FOR K5=1 TO 10
00180R(K5)=0
00190NEXT K5
00200V9=0
00210K3=0
00220GOSUB 5160
00230X=0
00240V8=0
00250V7=0
00260S8=0
00270S7=0
00280K6=0
00290DIM P(10)
00300I2=1
00310PRINT L$
00320PRINT " PRIOR DISTRIBUTION - MULTINOMIAL MODEL"
00330PRINT
00340PRINT "THIS MODULE WILL ASSIST YOU IN FITTING A PRIOR DISTRIBUTION TO"
00350PRINT "YOUR BELIEFS ABOUT THE PROPORTIONS IN A MULTINOMIAL ANALYSIS."
00360PRINT "THE MODULE ATTEMPTS TO FIT A DIRICHLET DISTRIBUTION TO YOUR"
00370PRINT "BELIEFS ABOUT THE PROPORTIONS CONSIDERED JOINTLY, AND BETA"
00380PRINT "DISTRIBUTIONS TO YOUR BELIEFS ABOUT THE PROPORTIONS CONSIDERED"
00390PRINT "SEPARATELY. THERE CAN BE AS MANY AS 10 CATEGORIES."
00400PRINT
00410PRINT "HOW MANY CATEGORIES ARE THERE IN YOUR ANALYSIS";
00420GOSUB 9000
00430S0=O1
00440PRINT
00450IF ABS(S0-INT(S0))<.0001 THEN 480
00460PRINT "REENTER. NUMBER OF CATEGORIES MUST BE AN INTEGER."
00470GOTO 420
00480IF S0<3 THEN 500
00490IF S0<11 THEN 630
00500PRINT "THE NUMBER OF CATEGORIES MUST BE GREATER THAN 2 BUT NOT GREATER"
00510PRINT "THAN 10. IF THERE ARE ONLY TWO CATEGORIES THEN YOU SHOULD SELECT"
00520PRINT "A BINARY MODEL FOR YOUR ANALYSIS."
00530PRINT
00540PRINT "IF YOU WANT TO RESPECIFY THE NUMBER OF CATEGORIES TYPE '1'."
00550PRINT "IF YOU WANT TO EXIT THE MODULE TYPE '0'."
00560GOSUB 9000
00570IF O1=1 THEN 410
00580IF O1=0 THEN 600
00590GOTO 410
00600CHAIN "RSTRT"
00610PRINT "REENTER. INPUT MUST 0 OR 1."
00620GOTO 560
00630GOSUB 700
00640GOTO 990
00650REM ****************************************************************
00660REM
00670REM ROUTINE FOR INPUTTING JOINT MODE ESTIMATE
00680REM
00690REM ****************************************************************
00700PRINT "INPUT YOUR ESTIMATES OF THE PROPORTION OF THE POPULATION IN"
00710PRINT "EACH CATEGORY. THESE ESTIMATES SHOULD SUM TO 1.0"
00720S1=0
00730FOR I1=1 TO S0
00740PRINT
00750REM
00760 PRINT "CATEGORY ";I1;
00770GOSUB 9000
00780P(I1)=O1
00790IF P(I1) <= 0 THEN 820
00800IF P(I1) >= 1 THEN 820
00810GOTO 840
00820PRINT "REENTER. ESTIMATE MUST BE BETWEEN 0 AND 1."
00830GOTO 750
00840S1=S1+P(I1)
00850IF I1=1 THEN 880
00860IF I1=S0 THEN 880
00870 PRINT " ACCUMULATED PROPORTION=";S1;
00880NEXT I1
00890: ACCUMULATED PROPORTION = ##.##
00900IF ABS(S1-1)>.0001 THEN 920
00910RETURN
00920PRINT
00930PRINT "PROPORTIONS DO NOT SUM TO 1.0. RESPECIFY."
00940PRINT
00950GOTO 720
00960REM ************* END OF ROUTINE - INPUT JOINT MODE ***************
00970REM
00980REM
00990PRINT L$
01000PRINT " SPECIFY THE 25TH AND 75TH PERCENTILES OF YOUR PRIOR"
01010PRINT "DISTRIBUTION ON THE PROPORTION IN EACH OF THE CATEGORIES."
01020PRINT "WE ARE ASSUMING YOUR ESTIMATE WAS A MEASURE OF THE CENTRAL"
01030PRINT "TENDENCY OF THE PRIOR DISTRIBUTION."
01040PRINT
01050PRINT "CATEGORY ESTIMATE"
01060FOR K5=1 TO S0
01070PRINT USING 1080,K5,P(K5)
01080: ## ##.##
01090PRINT " 25TH";
01100GOSUB 9000
01110Q1=O1
01120IF O1 <= 0 THEN 1140
01130IF O1<P(K5) THEN 1160
01140PRINT "REENTER. 25TH MUST BE LESS THAN ESTIMATE AND MORE THAN 0."
01150GOTO 1100
01160PRINT " 75TH";
01170GOSUB 9000
01180Q3=O1
01190IF O1 >= 1 THEN 1210
01200IF O1>P(K5) THEN 1230
01210PRINT "REENTER. 75TH MUST BE MORE THAN ESTIMATE AND LESS THAN 1."
01220GOTO 1170
01230Q2=P(K5)
01240GOSUB 1360
01250K3=K3+W8
01260NEXT K5
01270W8=K3/S0
01280IF W8<8000 THEN 1300
01290W8=8000
01300GOTO 1550
01310REM *****************************************************************
01320REM
01330REM ESTIMATE A AND B PARAMETERS FROM PERCENTILES
01340REM
01350REM *******************************************************************
01360D1=SQR(Q2*(1-Q1))-SQR(Q1*(1-Q2))
01370D1=D1*D1
01380D3=SQR(Q2*(1-Q3))-SQR(Q3*(1-Q2))
01390D3=D3*D3
01400C=.056*(1/D1+1/D3)
01410A=C*Q2+1/3
01420IF A>1.15 THEN 1510
01430A=2
01440W8=A/P(K5)
01450W7=(A-1)/P(K5)+S0
01460IF W8>W7 THEN 1480
01470W8=W7
01480IF W8-A>1.15 THEN 1540
01490A=A+1
01500GOTO 1440
01510B=C*(1-Q2)+1/3
01520IF B <= 1.15 THEN 1430
01530W8=A+B
01540RETURN
01550PRINT L$
01560PRINT " HERE ARE THE PERCENTILES OF THE MARGINAL BETA DISTRIBUTIONS"
01570PRINT "FITTED TO YOUR SPECIFICATIONS."
01580PRINT
01590PRINT USING 1600,W8
01600: HYPOTHETICAL SAMPLE SIZE (A) =####.##
01610IF S0>8 THEN 1630
01620PRINT
01630PRINT"CATEGORY JOINT ESTIMATE 25TH 50TH 75TH"
01640REM
01650REM ******************************************************************
01660FOR K5=1 TO S0
01670R(K5)=.5*(P(K5)*W8+(W8-S0)*P(K5)+1)
01680A=R(K5)
01690B=W8-A
01700GOSUB 5220
01710P1=25
01720GOSUB 3030
01730Q(1)=J2
01740P1=50
01750GOSUB 3050
01760Q(2)=J2
01770P1=75
01780GOSUB 3050
01790IF J2<.99 THEN 1810
01800J2=.99
01810Q(3)=J2
01820PRINT USING 1830,K5,P(K5),Q(1),Q(2),Q(3)
01830: ## ##.## ##.## ##.## ##.##
01840NEXT K5
01850PRINT
01860IF V8=1 THEN 1940
01870PRINT "IF YOU DO NOT FEEL THAT THE HYPOTHETICAL SAMPLE SIZE (A) REFLECTS"
01880PRINT "YOUR PRIOR INFORMATION ABOUT THE PROPORTIONS YOU CAN SPECIFY A"
01890PRINT "DIFFERENT (A). THIS WILL NOT AFFECT THE JOINT ESTIMATE BUT WILL"
01900PRINT "CHANGE THE MARGINAL PERCENTILES. A LARGER (A) WILL RESULT IN "
01910PRINT "SMALLER INTERPERCENTILE DIFFERENCES AND A SMALLER (A) IN LARGER"
01920PRINT "ONES."
01930PRINT
01940PRINT "IF YOU WANT TO CHANGE (A) TYPE NEW VALUE ELSE '0'.";
01950GOSUB 9000
01960IF O1>9999 THEN 2140
01970IF O1=0 THEN 2160
01980V8=1
01990K7=1
02000FOR K5=2 TO S0
02010IF P(K7)<P(K5) THEN 2030
02020K7=K5
02030NEXT K5
02040K7=INT(.5/P(K7)+S0+1)
02050IF K7<W8 THEN 2070
02060K7=W8
02070IF O1<K7 THEN 2100
02080W8=O1
02090GOTO 1550
02100PRINT
02110PRINT USING 2120,K7
02120:REENTER. MINIMUM (A)=###.### MAX (A)=9999
02130GOTO 1950
02140PRINT "REENTER. MAXIMUM (A) = 9999"
02150GOTO 1950
02160PRINT
02170PRINT "IF YOU WANT TO SPECIFY A DIFFERENT JOINT ESTIMATE TYPE '1'."
02180PRINT "IF YOU DO NOT TYPE '0'."
02190GOSUB 9000
02200IF O1=0 THEN 2360
02210IF O1=1 THEN 2240
02220PRINT "REEENTER. INPUT MUST BE 0 OR 1."
02230GOTO 2190
02240PRINT L$
02250PRINT
02260PRINT "PRESENT JOINT ESTIMATE."
02270FOR K5=1 TO S0
02280PRINTP(K5);
02300NEXT K5
02310PRINT " "
02320PRINT "INPUT YOUR NEW JOINT ESTIMATE."
02330S8=1
02340GOSUB 720
02350GOTO 990
02360PRINT L$
02370PRINT " SUMMARY OF PRIOR DISTRIBUTIONS"
02380PRINT
02390PRINT USING 2410,W8
02400PRINT
02410: HYPOTHETICAL SAMPLE SIZE (A)=####.##
02420PRINT" JOINT MARGINAL PERCENTILES"
02430REM
02440: ## #####.## ##.## ##.## ##.## ##.## ##.## ##.## ##.##
02450 PRINT "CATEGORY P MEAN MODE";
02452 PRINT " 10TH 25TH 50TH 75TH 90TH"
02470FOR K5=1 TO S0
02480A=R(K5)
02490B=W8-A
02500GOSUB 5220
02510M4=(A-1)/(W8-S0)
02520M5=A/W8
02530P1=5
02540GOSUB 3030
02550Q(1)=J2
02560P1=10
02570GOSUB 3050
02580Q(2)=J2
02590P1=25
02600GOSUB 3050
02610Q(3)=J2
02620P1=50
02630GOSUB 3050
02640Q(4)=J2
02650P1=75
02660GOSUB 3050
02670IF J2<.99 THEN 2690
02680J2=.99
02690Q(5)=J2
02700P1=90
02710GOSUB 3050
02720IF J2<.99 THEN 2740
02730J2=.99
02740Q(6)=J2
02750GOTO 2810
02760P1=95
02770GOSUB 3050
02780IF J2<.99 THEN 2800
02790J2=.99
02800Q(7)=J2
02810PRINT USING 2440,K5,A,M5,M4,Q(2),Q(3),Q(4),Q(5),Q(6)
02820NEXT K5
02830SCRATCH#2
02831FOR I=1 TO 10
02833 PRINT #2,R(I)
02834NEXTI
02840SCRATCH#3
02841 PRINT # 3,S0
02850PRINT "* P DENOTES THE PARAMETERS OF THE DIRICHLET PRIOR"
02860PRINT
02870PRINT "IF YOU ARE NOT GOING DIRECTLY TO THE POSTERIOR ANALYSIS"
02880PRINT "YOU SHOULD COPY DOWN THE PARAMETERS OF THE JOINT DISTRIBUTION."
02890PRINT "IF YOU WANT TO GO DIRECTLY TO THE POSTERIOR ANALYSIS TYPE '1'."
02900PRINT "IF YOU DO NOT TYPE '0'."
02910GOSUB 9000
02920PRINT
02930IF O1=1 THEN 2990
02940IF O1=0 THEN 2970
02950PRINT "REENTER. INPUT MUST BE 0 OR 1."
02960GOTO 2910
02970REM
02980CHAIN "RSTRT"
02990CHAIN "CMODH"
03000PRINT "YOU MUST SELECT A VALUE BETWEEN .01 AND .99."
03010GOSUB 5160
03020GOSUB 7000
03030E1=0
03040A7=0
03050IF P1>50 THEN 3140
03060E2=A/(A+B)
03070IF E2>(A-1)/(A+B-2) THEN 3090
03080E2=(A-1)/(A+B-2)
03090J1=0
03100J2=E2
03110GOSUB 5000
03120A8=P
03130GOTO 3160
03140A8=.999999
03150E2=1
03160A6=P1/100
03170S9=.85
03180IF S9<(A6-A7)/(A8-A7) THEN 3200
03190S9=(A6-A7)/(A8-A7)
03200IF S9>.15 THEN 3220
03210S9=.15
03220J2=E1+S9*(E2-E1)
03230J1=0
03240GOSUB 5000
03250IF ABS(E1-E2)<.0001 THEN 3340
03260IF ABS(P-A6)<.001 THEN 3340
03270IF P>A6 THEN 3310
03280E1=J2
03290A7=P
03300GOTO 3170
03310E2=J2
03320A8=P
03330GOTO 3170
03340RETURN
03350REM*************************************************************
03360REM APPENDED GOSUBS FOLLOW
03370REM
05000REM ***************************************************
05005REM BETA CDF ROUTINE
05010REM INPUT A B J2
05015REM OUTPUT P
05020REM GOSUB'S TO BE CALLED PRIOR 5160 AND 5220
05025DIM W(16),O(16)
05030GOSUB 5340
05035IF A+B>85 THEN 5280
05040P=0
05045C6=0
05050IF A>1 THEN 5080
05055C6=A
05060C7=B
05065A=C7
05070B=C6
05075J2=1-J2
05080D0=(J2-J1)*.5
05085D1=(J1+J2)*.5
05090FOR I1=1 TO 16
05095D9=D0*O(I1)+D1
05100IF D9=0 THEN 5115
05105IF D9=1 THEN 5115
05107D9=LOG(D9)*(A-1)+LOG(1-D9)*(B-1)
05108IF D9<-80 THEN 5115
05110P=P+W(I1)*EXP(D9)
05115NEXT I1
05120P=P*F0
05125P=P*D0
05130IF C6=0 THEN 5155
05135A=C6
05140B=C7
05145P=1-P
05150J2=1-J2
05155RETURN
05160FOR I1=1 TO 16
05165READ W(I1),O(I1)
05170NEXT I1
05175DATA 2.71525E-02,-.989401
05180DATA 6.22535E-02,-.944575,9.51585E-02,-.865631
05185DATA .124629,-.755404,.149596,-.617876
05190DATA .169156,-.458017,.182603,-.281604,.189451,-9.50125E-02
05195DATA .189451,9.50125E-02,.182603,.281604,.169156,.458017
05200DATA .149596,.617876,.124629,.755404
05205DATA 9.51585E-02,.865631,6.22535E-02,.944575,2.71525E-02
05210DATA .989401
05215RETURN
05220G9=A+B
05225GOSUB 5850
05230F0=G0
05235G9=A
05240GOSUB 5850
05245F0=F0-G0
05250G9=B
05255GOSUB 5850
05260F0=F0-G0
05265IF A+B>85 THEN 5275
05270F0=EXP(F0)
05275RETURN
05280W1=(B*J2)^(1/3)
05285W2=(A*(1-J2))^(1/3)
05290GOSUB 5325
05295I1=P
05300W1=(B*J1)^(1/3)
05305W2=(A*(1-J1))^(1/3)
05310GOSUB 5325
05315P=I1-P
05320RETURN
05325Y3=3*(W1*(1-1/9/B)-W2*(1-1/9/A))/SQR(W1*W1/B+W2*W2/A)
05330GOSUB 8000
05335RETURN
05340REM 2/16/76 CHANGED TO ALL LOG
05345D2=F0+(A-1)*LOG(J2)+(B-1)*LOG(1-J2)
05350IF D2<-80 THEN 5370
05355IF D2>85 THEN 5380
05360D2=EXP(D2)
05365RETURN
05370D2=1.E-37
05375RETURN
05380D2=1.E+37
05385RETURN
05390REM
05395REM END OF BETA CDF ROUTINE
05400REM *******************************************************
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 ****************************************************
07000REM**************************************************************
07002REM BETA HDR ROUTINE
07004REM INPUT A B J5
07006REM OUTPUT J1 J2
07008U0=0
07010IF A<B THEN 7020
07012U0=1
07014J7=A
07016A=B
07018B=J7
07020J8=2
07022GOSUB 5220
07024J=(A-1)/(B-1)
07026J7=1
07028GOSUB 7096
07030IF J8<-99 THEN 7040
07032P9=0
07034GOSUB 5000
07036P0=P
07038IF ABS(P0-J5)>.001 THEN 7042
07040GOTO 7108
07042IF P0>J5 THEN 7052
07044J7=J8
07046J8=J8*2
07048P9=P0
07050GOTO 7028
07052J9=J7
07054J0=J8
07056J8=(J5-P9)/(P0-P9)
07058IF J8>.15 THEN 7064
07060J8=.15
07062GOTO 7068
07064IF J8<.85 THEN 7068
07066J8=.85
07068J8=J8*(J0-J9)+J9
07070GOSUB 7096
07072IF J8<-99 THEN 7080
07074GOSUB 5000
07076J3=P
07077IF ABS(J1-J2)<.0001 THEN 7108
07078IF ABS(J3-J5)>.001 THEN 7082
07080GOTO 7108
07082IF J3>J5 THEN 7090
07084J9=J8
07086P9=J3
07088GOTO 7056
07090J0=J8
07092P0=J3
07094GOTO 7056
07096J1=(J8^J-1)/(J8^(J+1)-1)
07098J2=J8*J1
07100IF J2<1 THEN 7106
07102J8=J8*.95
07104GOTO 7100
07106RETURN
07108IF U0=0 THEN 7122
07110J7=A
07112A=B
07114B=J7
07116J7=J2
07118J2=1-J1
07120J1=1-J7
07122RETURN
07124REM END OF BETA HDR ROUTINE
07126REM*************************************************
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