Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0113/cmodh.bas
There are 2 other files named cmodh.bas in the archive. Click here to see a list.
00015REM ************************************************************
00020REM
00030REM POSTERIOR MULTINOMIAL
00040REM
00050REM*********************************************************
00060X=0
00070REM--NECESSARY FOR NORMAL APPROX
00080FILES RFILE1,RFILE2,RFILE3
00120RESTORE#1
00121 INPUT# 1,I1,I2,I3
00130 SCRATCH#1
00131 PRINT # 1,17,I2,I3
00140DIM Q(10)
00150DIM R(10)
00155RESTORE#2
00156FORI=1 TO 10
00157INPUT#2,R(I)
00160NEXTI
00170RESTORE#3
00171 INPUT# 3,S0
00180GOSUB 5160
00190GOTO 600
00200PRINT" JOINT MARGINAL PERCENTILES"
00220: ## ###.## ##.## ##.## ##.## ##.## ##.## ##.## ##.##
00230PRINT"CATEGORY P* MEAN MODE ";
00240 PRINT "10TH 25TH 50TH 75TH 90TH"
00250FOR K5=1 TO S0
00260A=R(K5)
00270B=W8-A
00280GOSUB 5220
00290M4=(A-1)/(W8-S0)
00300M5=A/W8
00310P1=5
00320GOSUB 1360
00330Q(1)=J2
00340P1=10
00350GOSUB 1380
00360Q(2)=J2
00370P1=25
00380GOSUB 1380
00390Q(3)=J2
00400P1=50
00410GOSUB 1380
00420Q(4)=J2
00430P1=75
00440GOSUB 1380
00450Q(5)=J2
00460P1=90
00470GOSUB 1380
00472IF J2<.99 THEN 480
00474J2=.99
00480Q(6)=J2
00490GOTO 530
00500P1=95
00510GOSUB 1380
00512IF J2<.99 THEN 520
00514J2=.99
00520Q(7)=J2
00530PRINT USING 220,K5,A,M5,M4,Q(2),Q(3),Q(4),Q(5),Q(6)
00540NEXT K5
00550PRINT
00552PRINT "* PARAMETERS OF THE POSTERIOR DIRICHLET JOINT DISTRIBUTION"
00554PRINT
00560PRINT "THIS COMPLETES THE MULTINOMIAL ANALYSIS."
00570PRINT "WHEN YOU ARE READY TO CONTINUE TYPE '1'.";
00580GOSUB 9000
00590CHAIN "RSTRT"
00600PRINT L$
00610PRINT " POSTERIOR ANALYSIS-MULTINOMIAL MODEL"
00620PRINT
00630PRINT "THIS MODULE COMPUTES THE JOINT AND MARGINAL POSTERIOR"
00640PRINT "DISTRIBUTIONS FOR THE PROPORTIONS."
00650PRINT
00660IF S0 <> 0 THEN 990
00670PRINT "INPUT NUMBER OF CATEGORIES IN YOUR ANAYLSIS.";
00680GOSUB 9000
00690S0=O1
00700IF ABS(S0-INT(S0))<.0001 THEN 730
00710PRINT "REENTER. NUMBER OF CATEGORIES MUST BE INTEGER."
00720GOTO 680
00730IF S0>2 THEN 800
00740PRINT "IF THERE ARE ONLY TWO CATEGORIES IN YOUR ANALYSIS YOU"
00750PRINT "SHOULD USE A BINARY MODEL."
00760PRINT "IF YOU HAVE MORE THAN TWO THEN INPUT AGAIN."
00770PRINT "IF NOT TYPE '0'."
00780GOSUB 9000
00790IF O1=0 THEN 590
00795GOTO 690
00800IF S0 <= 10 THEN 840
00810PRINT "THIS MODULE CAN ONLY HANDLE 10 OR FEWER. "
00820PRINT "IF YOU WANT TO RESPECIFY TYPE THE NEW NUMBER."
00830PRINT "IF YOU DO NOT TYPE '0'."
00834IF O1=0 THEN 590
00836GOTO 690
00840PRINT
00850PRINT "INPUT THE PARAMETERS OF YOUR PRIOR DIRICHLET DISTRIBUTION"
00860PRINT "BY CATEGORY."
00870PRINT
00880FOR K5=1 TO S0
00890:CATEGORY ##
00900PRINT"CATEGORY ";K5;
00910GOSUB 9000
00920PRINT
00930IF O1>1.1 THEN 960
00940PRINT "REENTER. EACH PARAMETER MUST BE GREATER THAN 1.1."
00950GOTO 910
00960R(K5)=O1
00970NEXT K5
00980PRINT
00990PRINT
01000PRINT "ENTER THE NUMBER OF OBSERVATIONS IN EACH CATEGORY."
01010W8=0
01020PRINT
01030FOR K5=1 TO S0
01040REM
01050PRINT"CATEGORY ";K5;
01060GOSUB 9000
01063IF O1 >= 0 THEN 1070
01065PRINT "REENTER. CAN NOT BE LESS THAN 0."
01067GOTO 1040
01070PRINT
01080N(K5)=O1
01090NEXT K5
01100PRINT
01110PRINT "HERE IS WHAT YOU ENTERED."
01120FOR K5=1 TO S0
01125:CATEGORY ## #####
01130PRINT USING 1125,K5,N(K5)
01140NEXT K5
01150PRINT "IF YOU WANT TO REENTER THE SAMPLE DATA TYPE '1' ELSE '0'.";
01170GOSUB 9000
01175PRINT L$
01180IF O1=0 THEN 1210
01190IF O1=1 THEN 990
01210PRINT " SUMMARY OF POSTERIOR DISTRIBUTIONS"
01220PRINT
01230FOR K5=1 TO S0
01240R(K5)=R(K5)+N(K5)
01250W8=W8+R(K5)
01260NEXT K5
01270GOTO 200
01360E1=0
01370A7=0
01380IF P1>50 THEN 1450
01390E2=A/(A+B)
01392IF E2>(A-1)/(A+B-2) THEN 1400
01395E2=(A-1)/(A+B-2)
01400J1=0
01410J2=E2
01420GOSUB 5000
01430A8=P
01440GOTO 1470
01450A8=.999999
01460E2=1
01470A6=P1/100
01480S9=(A6-A7)/(A8-A7)
01485IF S9<.85 THEN 1490
01488S9=.85
01489GOTO 1500
01490IF S9>.15 THEN 1500
01495S9=.15
01500J2=E1+S9*(E2-E1)
01510J1=0
01520GOSUB 5000
01530IF ABS(P-A6)<.001 THEN 1610
01540IF P>A6 THEN 1580
01550E1=J2
01560A7=P
01570GOTO 1480
01580E2=J2
01590A8=P
01600GOTO 1480
01610RETURN
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