C VERSION PURGE5 -- THIS IS A DOUBLE PRECISION VERSION OF C =================================================== C CD PURGING WAS ADDED BY YU XIE AT THE C CETER FOR DEOGRAPHY AND ECOLOGY, UNIVERSITY OF C WISCONSIN-MADISON, 1988 C =================================================== C PURGE4 -- ALL ELSE REMAINS THE SAME C C NEW PURGE PROGRAM OF JUNE 1986 C METHOD DESCRIBED IN CLOGG (1978), DEMOGRAPHY; IN CLOGG AND C SHOCKEY (1985),DEMOGRAPHY; IN CLOGG, SHOCKEY, AND ELIASON (1985) C POP. ASSOC. PAPER. C C MODIFICATIONS OF JUNE 1986 FOR FULL JACKKNIFING ON CONTRASTS C C FURTHER MODIFIED BY L. SANTI AT UW/CDE IN DECEMBER, 1987. C THESE MODIFICATIONS EXPANDED THE PROGRAM'S LIMITATIONS FROM C 20 COMPOSITION CATEGORIES, 5 GROUPS, AND 5 LEVELS OF THE C DEPENDENT (RATE) VARIABLE TO 100 C, 100 G, AND 10 D. C THESE PROGRAM MODIFICATIONS NECESSITATED A CHANGE IN THE FORMAT C OF CONTROL CARD 2; INSTEAD OF READING THE 8 PARAMETERS IN (8I2) C FORMAT, THE PROGRAM WILL BE LOOKING FOR (2I3,6I2). C C CONTROL CARDS: C CARD 1: TITLE. C CARD 2: 8 PARAMETERS IN I2 FORMAT: C COL 2: NO. OF COMPOSITION CATEGORIES, C COL 4: NO. OF GROUPS, C COL 6: NO. OF LEVELS OF DEPENDENT VARIABLE C COL 8: SUPPRESS FREQUENCIES (1=YES), C COL 10: SUPPRESS TAUS FOR CG INTERACTION (1=YES), C COL 12: USE FITS FROM MODEL OF NO 3-FACTOR (1=YES), C COL 14: JACKKNIFE ESTIMATES OF RATES AND ST.DEV (1=YES), C COL 16: NUMBER OF CONTRAST VECTORS TO BE INPUTTED C CARD 3: INPUT FORMAT (E.G., (10F4.0)). FREQUENCIES C VARY FIRST BY LEVEL OF DEP, THEN BY LEVEL OF COMP, C THEN BY GROUP. C CARD 4: DATA IN ABOVE FORMAT, AS MANY "CARDS" AS NECESSARY C CARD 5: CONTRAST COEFFICIENTS IN (10F8.6) FORMAT-- C INPUT WITH GROUPS = CARDS, COLS = CON COEFF. NO C CARDS HERE IF NO CONTRASTS REQUESTED. C CARD 6: TWO PARAMETERS IN I2 FORMAT. FIRST TELLS TYPE C OF RATE (0=CRUDE;1=PARTIAL CG; 2=PARTIAL CG AND C CGD; 3=MARGINAL CG, 4=DIRECT STANDARDIZATION, 5= C MARGINAL CG AND 3-FACTOR CGD), 6=PARTIAL CD, C 7=PARTIAL CD AND CGD. SECOND PARM TELLS C STANDARD GROUP TO BE USED (0=NO STANDARD). REPEAT C CARD 6'S AS MUCH AS DESIRED. C PRESENT RE-CALCULATION OF RATES FOR EACH CARD 5 DONE SO THAT C JACKKNIFING CAN BE DONE EXPEDIENTLY. DOUBLE PRECISION FMT(20),TITLE(20),F(100,100,10), *R(100,10),XG(100),FP(100,100,10), *AF(100,100,10),RS(100,10),RSS(100,10),RSL(100,10), *RSSL(100,10),CON(100,10),RL(100,10), *CONV(100,10),CONVL(100,10),CONVS(100,10), *CONVSS(100,10),CONVLS(100,10), *CNVLSS(100,10),XN C MOST ARRAYS ABOVE ARE FOR JACKKNIFE WORK. SOME SIMPLIFICATION IS C POSSIBLE. USE OF ARRAYS AS FOLLOWS: F, FREQUENCIES; R, RATES; XG, C GROUP TOTALS (FOR RESCALING); FP, PURGED F'S; AF, WORKING FREQS; C CON, CONTRAST COEFF;CONV, VALUE OF CONTRASTS; ARRAYS ENDING IN S C FOR JACKKNIFE ON RATE ESTIMATE, ARRAYS ENDING IN SS FOR SUM OF SQS C USED IN CALCULATING STANDARD DEVIATIONS DATA XG/100*0.0/,CON/1000*0.0/,CONV/1000*0.0/,CONVL/1000*0.0/ C IOFF=1 WRITE(6,10) 10 FORMAT(//,' PROGRAM FOR PURGED AND STANDARDIZED RATES') READ(5,15)TITLE 15 FORMAT(20A4) WRITE(6,16)TITLE 16 FORMAT(//,1X,20A4) C READ DIMENSIONS OF TABLE AND SET OPTIONS FOR RUN READ(5,20)IC,IG,ID,NOF,NOP,IFIT,JACK,NCON 20 FORMAT(2I3,6I2) WRITE(6,25)IC,IG,ID,NOF,NOP,IFIT,JACK,NCON 25 FORMAT(/,' INPUT DATA HAS: ', */,' COMPOSITION CATEGORIES............',I4, */,' GROUPS............................',I4, */,' LEVELS OF DEPENDENT VARIABLE......',I4, */,' OPTIONS SPECIFIED FOR RUN: ', */,' SUPPRESS FREQUENCIES (1=Y)........',I4, */,' SUPPRESS TAUS AND LAMBDAS (1=Y)...',I4, */,' USE FITS FROM NO 3-FACTOR MODEL...',I4, */,' JACKKNIFE CALCULATION OF VARIANCE.',I4, */,' NUMBER OF CONTRASTS = ............',I4) READ(5,15)FMT READ(5,FMT)(((F(I,J,K),K=1,ID),I=1,IC),J=1,IG) C DATA ARE IN ARRAY F C CORRECT FOR ZERO ENTRIES--NOTE THAT THIS WILL APPLY TO ALL TYPES DO 40 I=1,IC DO 40 J=1,IG DO 40 K=1,ID AF(I,J,K)=F(I,J,K) IF(F(I,J,K).GT.0.0)GO TO 40 IF(IFIT.EQ.0)AF(I,J,K)=.5 WRITE(6,41)I,J,K 41 FORMAT(/,' ZERO FREQUENCY FOR COMP =',I4,/,20X,' GROUP =',I4,/, *20X,' DEP LEVEL = ',I4) 40 CONTINUE 220 CONTINUE IF(NCON.EQ.0)GO TO 55 WRITE(6,8) 8 FORMAT(/,' GROUP CONTRASTS',//) DO 4 J=1,IG READ(5,*)(CON(I,J),I=1,NCON) WRITE(6,7)J,(CON(I,J),I=1,NCON) 7 FORMAT(//,' GROUP ',I4,10F8.4) 4 CONTINUE 55 CONTINUE C GET GROUP TOTALS FOR LATER USE DO 5 I=1,IC DO 5 K=1,ID DO 5 J=1,IG XG(J)=XG(J) + AF(I,J,K) 5 CONTINUE KOUNT=0 IF(IFIT.NE.0)CALL FIT(AF,IC,IG,ID,IOFF) 200 CONTINUE READ(5,30,END=999)ITYPE,ISTD 30 FORMAT(5I2) IF(ITYPE.EQ.0)WRITE(6,36) 36 FORMAT('1UNADJUSTED OR CRUDE RATES') IF(ITYPE.EQ.1)WRITE(6,37)ISTD 37 FORMAT('1PARTIAL CG PURGED RATES',/,' STANDARD GROUP = ',I4) IF(ITYPE.EQ.2)WRITE(6,38)ISTD 38 FORMAT('1PARTIAL CG AND CGD PURGED RATES',/,' STANDARD GROUP=', *I4) IF(ITYPE.EQ.3)WRITE(6,31)ISTD 31 FORMAT('1MARGINAL CG PURGED RATES',/,' STANDARD GROUP = ',I4) IF(ITYPE.EQ.4)WRITE(6,39)ISTD 39 FORMAT('1DIRECT STANDARDIZATION',/,' STANDARD GROUP = ',I4) IF(ITYPE.EQ.5)WRITE(6,32)ISTD 32 FORMAT('1RATES PURGED OF MARGINAL CG AND 3-FACTOR CGD', */,' STANDARD GROUP = ',I4) IF(ITYPE.EQ.6)WRITE(6,321)ISTD 321 FORMAT('1PARTIAL CD PURGED RATES',/,' STANDARD GROUP = ',I4) IF(ITYPE.EQ.7)WRITE(6,322)ISTD 322 FORMAT('1PARTIAL CD-CGD PURGED RATES',/,' STANDARD GROUP = ',I4) CALL PURGE(AF,FP,IC,IG,ID,ITYPE,ISTD,NOP) CALL RESCAL(FP,XG,IC,IG,ID) IF(NOF.EQ.0)CALL PRINT(FP,XG,IC,IG,ID) CALL RATE(FP,R,RL,CONV,CONVL,CON,XG,IC,IG,ID,NCON) WRITE(6,35) 35 FORMAT(//,' SUMMARY GROUP RATES') CALL PRATE(R,IG,ID) WRITE(6,305) 305 FORMAT(//,' SUMMARY GROUP LOG-RATES') CALL PRATE(RL,IG,ID) WRITE(6,306) 306 FORMAT(//,' CONTRAST VALUES: ARITHMETIC, THEN LOGS:GROUP=CON') IF(NCON.NE.0)CALL PRATE(CONV,NCON,ID) IF(NCON.NE.0)CALL PRATE(CONVL,NCON,ID) KOUNT=KOUNT+1 IF(JACK.EQ.0)GO TO 201 C BEGIN JACKKNIFE CALCULATIONS AS A LOOP THROUGH ABOVE PROCEDURE WRITE(6,101) 101 FORMAT(/,' JACKKNIFE ESTIMATES FOLLOW-REGARDLESS OF TYPE OF RATE') C SET NOP=NOPJ=1 TO SUPPRESS PARAMETER PRINT THROUGH STAGES OF JACK NOPJ=1 IOFF=0 C INITIALIZE RATE AND SS VECTORS DO 98 J=1,IG DO 98 K=1,ID RS(J,K)=0.0 RSL(J,K)=0.0 RSSL(J,K)=0.0 CONVS(J,K)=0.0 CONVSS(J,K)=0.0 CONVLS(J,K)=0.0 CNVLSS(J,K)=0.0 RSS(J,K)=0.0 98 CONTINUE XN=0.0 DO 100 I=1,IC DO 100 J=1,IG DO 100 K=1,ID XN=XN+F(I,J,K) IF(IFIT.EQ.0)GO TO 120 DO 119 II=1,IC DO 119 JJ=1,IG DO 119 KK=1,ID AF(II,JJ,KK)=F(II,JJ,KK) 119 CONTINUE IF(AF(I,J,K).LT.1.0)GO TO 100 GO TO 121 120 IF(F(I,J,K).LE.1.0)GO TO 100 121 AF(I,J,K)=F(I,J,K) - 1.0 C CORRECT GROUP TOTAL, AND CORRECT THIS AFTER LOOP XG(J)=XG(J) - 1.0 IF(IFIT.NE.0)CALL FIT(AF,IC,IG,ID,IOFF) CALL PURGE(AF,FP,IC,IG,ID,ITYPE,ISTD,NOPJ) CALL RESCAL(FP,XG,IC,IG,ID) CALL RATE(FP,R,RL,CONV,CONVL,CON,XG,IC,IG,ID,NCON) DO 105 JJ=1,IG DO 105 KK=1,ID RS(JJ,KK)=RS(JJ,KK) + F(I,J,K)*R(JJ,KK) RSS(JJ,KK)=RSS(JJ,KK) + F(I,J,K)*(R(JJ,KK)**2) RSL(JJ,KK)=RSL(JJ,KK) + F(I,J,K)*RL(JJ,KK) RSSL(JJ,KK)=RSSL(JJ,KK) + F(I,J,K)*(RL(JJ,KK)**2) IF(NCON.EQ.0)GO TO 105 CONVS(JJ,KK)=CONVS(JJ,KK) + F(I,J,K)*CONV(JJ,KK) CONVSS(JJ,KK)=CONVSS(JJ,KK) + F(I,J,K)*(CONV(JJ,KK)**2) CONVLS(JJ,KK)=CONVLS(JJ,KK) + F(I,J,K)*CONVL(JJ,KK) CNVLSS(JJ,KK)=CNVLSS(JJ,KK) + F(I,J,K)*(CONVL(JJ,KK)**2) 105 CONTINUE XG(J)=XG(J)+1.0 AF(I,J,K)=AF(I,J,K) + 1.0 100 CONTINUE DO 110 J=1,IG DO 110 K=1,ID RS(J,K)=RS(J,K)/XN RSS(J,K)=((XN-1.0)*(RSS(J,K) - XN*(RS(J,K)**2.0)))/XN RSS(J,K)=DSQRT(RSS(J,K)) RSL(J,K)=RSL(J,K)/XN RSSL(J,K)=((XN-1.0)*(RSSL(J,K) - XN*(RSL(J,K)**2)))/XN RSSL(J,K)=DSQRT(RSSL(J,K)) 110 CONTINUE IF (NCON.EQ.0)GO TO 113 DO 112 J=1,NCON DO 112 K=1,ID CONVS(J,K)=CONVS(J,K)/XN CONVSS(J,K)=((XN-1.0)*(CONVSS(J,K)-XN*(CONVS(J,K)**2)))/XN CONVSS(J,K)=DSQRT(CONVSS(J,K)) CONVLS(J,K)=CONVLS(J,K)/XN CNVLSS(J,K)=((XN-1.0)*(CNVLSS(J,K)-XN*(CONVLS(J,K)**2)))/XN CNVLSS(J,K)=DSQRT(CNVLSS(J,K)) 112 CONTINUE 113 CONTINUE WRITE(6,115) 115 FORMAT(//,' JACKKNIFED ESTIMATES OF RATES') CALL PRATE(RS,IG,ID) WRITE(6,116) 116 FORMAT(//,' JACKKNIFED ESTIMATES OF STANDARD DEVIATIONS') CALL PRATE(RSS,IG,ID) WRITE(6,117) 117 FORMAT(//,' JACKKNIFE ESTIMATES OF LOG-RATES') CALL PRATE(RSL,IG,ID) WRITE(6,116) CALL PRATE(RSSL,IG,ID) IF(NCON.EQ.0)GO TO 201 WRITE(6,123) 123 FORMAT(//,' JACKKNIFE ESTIMATES OF CONTRASTS') CALL PRATE(CONVS,NCON,ID) WRITE(6,116) CALL PRATE(CONVSS,NCON,ID) WRITE(6,124) 124 FORMAT(//,' JACKKNIFE ESTIMATES OF CONTRASTS OF LOG-RATES') CALL PRATE(CONVLS,NCON,ID) WRITE(6,116) CALL PRATE(CNVLSS,NCON,ID) 201 GO TO 200 999 STOP END SUBROUTINE PRINT(F,XG,IC,IG,ID) DOUBLE PRECISION F(100,100,10),XG(100),XD(10),T C PRINTS ARRAY IN F WRITE(6,10) 10 FORMAT(//,' CELL FREQUENCIES') T=0.0 DO 21 J=1,IG T=T+XG(J) WRITE(6,25)J 25 FORMAT(//,' GROUP = ',I4) 31 FORMAT(/,' I= ',I4,4X,7F10.2) DO 20 I=1,IC WRITE(6,31)I,(F(I,J,K),K=1,ID) 20 CONTINUE WRITE(6,32)XG(J) 21 CONTINUE 32 FORMAT(/,' GROUP TOTAL IS = ',10X,F10.2) RETURN END SUBROUTINE RATE(F,R,RL,CONV,CONVL,CON,XG,IC,IG,ID,NCON) DOUBLE PRECISION F(100,100,10),XG(100),R(100,10), +RL(100,10),CONV(100,10), +CON(100,10),CONVL(100,10),T,TT C GET GROUP RATES. FIRST GET GROUP TOTALS. DO 10 J=1,IG XG(J)=0.0 10 CONTINUE DO 20 J=1,IG DO 20 I=1,IC DO 20 K=1,ID XG(J)=XG(J) + F(I,J,K) 20 CONTINUE DO 30 J=1,IG DO 30 K=1,ID T = 0.0 DO 31 I=1,IC T = T + F(I,J,K) 31 CONTINUE R(J,K)=T/XG(J) RL(J,K)=DLOG(R(J,K)) 30 CONTINUE IF(NCON.EQ.0)GO TO 50 C FILL IN CONTRAST VALUES, ABSOLUTE IN CONV, LOGS IN CONVL DO 40 K=1,ID DO 40 I=1,NCON T=0.0 TT=0.0 DO 43 J=1,IG T=T + CON(I,J)*R(J,K) TT=TT + CON(I,J)*RL(J,K) 43 CONTINUE CONV(I,K)=T CONVL(I,K)=TT 40 CONTINUE 50 RETURN END SUBROUTINE PRATE(R,IG,ID) DOUBLE PRECISION R(100,10) DO 10 J=1,IG WRITE(6,15)J,(R(J,K),K=1,ID) 10 CONTINUE 15 FORMAT(//,' GROUP = ',I4,/,4X, 8F10.6) RETURN END SUBROUTINE PURGE(F,FP,IC,IG,ID,ITYPE,ISTD,NOP) DOUBLE PRECISION F(100,100,10),AF(100,100,10), *FP(100,100,10),VI(100),VJ(100), *VIJ(100,100),VIK(100,10),VJK(100,10),VK(10), *T,V,X,Y,Z,XZ,YZ,XY,XYZ IF(ITYPE.NE.0) GO TO 4 DO 3 I=1,IC DO 3 J=1,IG DO 3 K=1,ID FP(I,J,K)=F(I,J,K) 3 CONTINUE GO TO 500 4 X=(IC) Y=(IG) Z=(ID) XZ=X*Z YZ=Y*Z XY=X*Y XYZ=XY*Z DO 5 I=1,IC VI(I)=0.0 DO 5 J=1,IG VIJ(I,J)=0.0 5 CONTINUE DO 6 J=1,IG VJ(J)=0.0 6 CONTINUE V=0.0 C DOES THE PURGING. ITYPE TELLS WHAT TO DO. IF(ITYPE.NE.1)GO TO 50 C CALCULATE NU'S DO 10 I=1,IC DO 10 J=1,IG DO 10 K=1,ID AF(I,J,K)=DLOG(F(I,J,K)) V=V + AF(I,J,K) 10 CONTINUE DO 12 I=1,IC DO 12 J=1,IG DO 12 K=1,ID VIJ(I,J)=VIJ(I,J) + AF(I,J,K) 12 CONTINUE DO 13 I=1,IC DO 13 J=1,IG VI(I)=VI(I) + VIJ(I,J) 13 CONTINUE DO 14 J=1,IG DO 14 I=1,IC VJ(J)=VJ(J) + VIJ(I,J) 14 CONTINUE DO 20 I=1,IC DO 20 J=1,IG IF(ISTD.EQ.0)FP(I,J,1)=(XY*VIJ(I,J) - X*VI(I) - Y*VJ(J) + V)/XYZ IF(ISTD.GT.0)FP(I,J,1)=(X*(VIJ(I,J)-VIJ(I,ISTD)) - (VJ(J)- *VJ(ISTD)))/XZ FP(I,J,2)=DEXP(FP(I,J,1)) 20 CONTINUE IF(NOP.NE.0)GO TO 40 WRITE (6,35) 35 FORMAT(//,' PARTIAL COMPOSITION-GROUP TAUS') DO 38 I=1,IC WRITE(6,39)I,(FP(I,J,2),J=1,IG) 38 CONTINUE 39 FORMAT(/,' COMP LEVEL =',2X,I3,8F10.3) 40 CONTINUE C NOW PURGE OF PARTIAL CG DO 45 I=1,IC DO 45 J=1,IG T=FP(I,J,2) DO 45 K=1,ID FP(I,J,K)=F(I,J,K)/T 45 CONTINUE GO TO 500 C INSERT ADJUSTMENT FOR PARTIAL CG AND 3-FACTOR CGD HERE 50 IF(ITYPE.NE.2)GO TO 80 C ADJUSTMENT FOR PARTIAL CG AND CGD DO 55 K=1,ID VK(K)=0.0 DO 56 I=1,IC 56 VIK(I,K)=0.0 DO 55 J=1,IG 55 VJK(J,K)=0.0 DO 58 I=1,IC DO 58 J=1,IG DO 58 K=1,ID 58 AF(I,J,K)=DLOG(F(I,J,K)) DO 60 I=1,IC DO 60 K=1,ID DO 60 J=1,IG 60 VIK(I,K)=VIK(I,K) +AF(I,J,K) DO 62 K=1,ID DO 62 I=1,IC 62 VK(K)=VK(K) + VIK(I,K) DO 64 J=1,IG DO 64 K=1,ID DO 64 I=1,IC 64 VJK(J,K)=VJK(J,K) + AF(I,J,K) DO 70 I=1,IC DO 70 J=1,IG DO 70 K=1,ID IF(ISTD.EQ.0)FP(I,J,K)=(XZ*VIK(I,K) + YZ*VJK(J,K) - Z*VK(K))/XYZ IF(ISTD.NE.0)FP(I,J,K)=AF(I,ISTD,K) -(VJK(ISTD,K) - VJK(J,K))/X 70 FP(I,J,K)=DEXP(FP(I,J,K)) GO TO 500 80 IF(ITYPE.NE.3)GO TO 100 C MARGINAL CG ADJUSTMENT 81 CONTINUE DO 85 I=1,IC DO 85 J=1,IG DO 85 K=1,ID FP(I,J,K)=0.0 85 AF(I,J,K)=0.0 DO 88 I=1,IC DO 88 J=1,IG DO 89 K=1,ID 89 FP(I,J,1)=FP(I,J,1) + F(I,J,K) AF(I,J,1)=DLOG(FP(I,J,1)) 88 V = V + AF(I,J,1) DO 90 I=1,IC DO 90 J=1,IG 90 VI(I)=VI(I) + AF(I,J,1) DO 91 J=1,IG DO 91 I=1,IC 91 VJ(J)=VJ(J) + AF(I,J,1) DO 93 I=1,IC DO 93 J=1,IG IF(ISTD.EQ.0)VIJ(I,J)=AF(I,J,1) - VI(I)/Y - VJ(J)/X + V/XY IF(ISTD.NE.0)VIJ(I,J)=(AF(I,J,1)-AF(I,ISTD,1)) -(VJ(J)-VJ(ISTD))/X 93 VIJ(I,J)=DEXP(VIJ(I,J)) IF(NOP.NE.0)GO TO 97 WRITE(6,83) 83 FORMAT(//,' MARGINAL CG TAUS') DO 96 I=1,IC WRITE(6,39)I,(VIJ(I,J),J=1,IG) 96 CONTINUE 97 CONTINUE DO 98 I=1,IC DO 98 J=1,IG DO 98 K=1,ID 98 FP(I,J,K)=F(I,J,K)/VIJ(I,J) IF(ITYPE.EQ.5)GO TO 151 GO TO 500 100 IF(ITYPE.NE.4)GO TO 150 C CALCULATE FREQUENCIES UNDER DIRECT STANDARDIZATION NEXT DO 105 I=1,IC DO 105 J=1,IG 105 VIJ(I,J)=0.0 DO 106 I=1,IC DO 106 J=1,IG DO 106 K=1,ID 106 VIJ(I,J)=VIJ(I,J) + F(I,J,K) C NOW GET STD DISTRIBUTION DO 108 I=1,IC 108 VI(I)=0.0 IF(ISTD.NE.0)GO TO 115 C CALCULATE AVERAGE COMPOSITION C FIRST GET GROUP TOTALS DO 109 J=1,IG DO 109 I=1,IC 109 VJ(J)=0.0 DO 103 J=1,IG DO 103 I=1,IC 103 VJ(J)=VJ(J)+VIJ(I,J) DO 110 I=1,IC T=0.0 DO 111 J=1,IG 111 T=T + VIJ(I,J)/VJ(J) 110 VI(I)=T*VJ(1)/Y GO TO 125 115 CONTINUE DO 120 I=1,IC 120 VI(I)=VIJ(I,ISTD) 125 CONTINUE DO 130 I=1,IC DO 130 J=1,IG DO 130 K=1,ID 130 FP(I,J,K)=VI(I)*F(I,J,K)/VIJ(I,J) GO TO 500 150 IF(ITYPE.NE.5)GO TO 502 GO TO 81 151 CONTINUE C ADJUSTMENT FOR MARGINAL CG AND 3-FACTOR CGD, IN 2 STEPS C FP HAS FREQUENCIES PURGED FOR MARGINAL CG V=0.0 DO 152 I=1,IC DO 152 J=1,IG DO 152 K=1,ID AF(I,J,K)=DLOG(F(I,J,K)) 152 V=V + AF(I,J,K) DO 153 I=1,IC VI(I)=0.0 DO 153 J=1,IG 153 VIJ(I,J)=0.0 DO 154 K=1,ID VK(K)=0.0 DO 154 J=1,IG 154 VJK(J,K)=0.0 DO 155 J=1,IG 155 VJ(J)=0.0 DO 156 I=1,IC DO 156 K=1,ID 156 VIK(I,K)=0.0 DO 160 I=1,IC DO 160 J=1,IG DO 160 K=1,ID 160 VIJ(I,J)=VIJ(I,J) + AF(I,J,K) DO 161 I=1,IC DO 161 J=1,IG 161 VI(I)=VI(I) + VIJ(I,J) DO 162 J=1,IG DO 162 I=1,IC 162 VJ(J)=VJ(J) + VIJ(I,J) DO 163 I=1,IC DO 163 K=1,ID DO 163 J=1,IG 163 VIK(I,K)=VIK(I,K) + AF(I,J,K) DO 164 K=1,ID DO 164 I=1,IC 164 VK(K)=VK(K) + VIK(I,K) DO 165 J=1,IG DO 165 K=1,ID DO 165 I=1,IC 165 VJK(J,K)=VJK(J,K) + AF(I,J,K) V=V/XYZ DO 170 I=1,IC DO 170 J=1,IG DO 170 K=1,ID IF(ISTD.EQ.0)T=AF(I,J,K)-VIJ(I,J)/Z-VIK(I,K)/Y-VJK(J,K)/X+ *VI(I)/YZ+VJ(J)/XZ+VK(K)/XY + V IF(ISTD.NE.0)T=(AF(I,J,K)-AF(I,ISTD,K)) -(VIJ(I,J)-VIJ(I,ISTD)) *-(VJK(J,K)-VJK(ISTD,K)) + (VJ(J)-VJ(ISTD)) C T HAS LAMBDA CGD-IJK T=DEXP(T) 170 FP(I,J,K)=FP(I,J,K)/T GO TO 500 C CD PARTIAL PURGING 502 IF(ITYPE.NE.6)GO TO 501 DO 102 K=1,ID VK(K)=0.0 DO 102 I=1,IC VIK(I,K)=0.0 102 CONTINUE C CALCULATE NU'S DO 101 I=1,IC DO 101 J=1,IG DO 101 K=1,ID AF(I,J,K)=DLOG(F(I,J,K)) V=V + AF(I,J,K) 101 CONTINUE DO 121 I=1,IC DO 121 K=1,ID DO 121 J=1,IG VIK(I,K)=VIK(I,K) + AF(I,J,K) 121 CONTINUE DO 131 I=1,IC DO 131 K=1,ID VI(I)=VI(I) + VIK(I,K) 131 CONTINUE DO 141 K=1,ID DO 141 I=1,IC VK(K)=VK(K) + VIK(I,K) 141 CONTINUE DO 201 I=1,IC DO 201 K=1,ID IF(ISTD.EQ.0)FP(I,1,K)=(XZ*VIK(I,K) - X*VI(I) - Z*VK(K) + V)/XYZ IF(ISTD.GT.0) GO TO 5001 FP(I,2,K)=DEXP(FP(I,1,K)) 201 CONTINUE IF(NOP.NE.0)GO TO 401 WRITE (6,351) 351 FORMAT(//,' PARTIAL COMPOSITION-DEPENDENT TAUS') DO 381 I=1,IC WRITE(6,391)I,(FP(I,2,K),K=1,ID) 381 CONTINUE 391 FORMAT(/,' COMP LEVEL =',2X,I3,8F10.3) 401 CONTINUE C NOW PURGE OF PARTIAL CD DO 451 I=1,IC DO 451 K=1,ID T=FP(I,2,K) DO 451 J=1,IG FP(I,J,K)=F(I,J,K)/T 451 CONTINUE GO TO 500 C INSERT ADJUSTMENT FOR PARTIAL CD AND 3-FACTOR CGD HERE 501 IF(ITYPE.NE.7)GO TO 500 C ADJUSTMENT FOR PARTIAL CD AND CGD DO 561 J=1,IG VJ(J)=0.0 DO 561 I=1,IC VIJ(I,J)=0.0 561 CONTINUE DO 551 K=1,ID VK(K)=0.0 DO 551 J=1,IG VJK(J,K)=0.0 551 CONTINUE DO 581 I=1,IC DO 581 J=1,IG DO 581 K=1,ID 581 AF(I,J,K)=DLOG(F(I,J,K)) DO 601 I=1,IC DO 601 J=1,IG DO 601 K=1,ID 601 VIJ(I,J)=VIJ(I,J) +AF(I,J,K) DO 621 J=1,IG DO 621 I=1,IC 621 VJ(J)=VJ(J) + VIJ(I,J) DO 641 J=1,IG DO 641 K=1,ID DO 641 I=1,IC 641 VJK(J,K)=VJK(J,K) + AF(I,J,K) DO 701 I=1,IC DO 701 J=1,IG DO 701 K=1,ID IF(ISTD.EQ.0)FP(I,J,K)=(XY*VIJ(I,J) + YZ*VJK(J,K) - Y*VJ(J))/XYZ IF(ISTD.NE.0) GOTO 5001 701 FP(I,J,K)=DEXP(FP(I,J,K)) GO TO 500 5001 WRITE (6,5002) 5002 FORMAT(//,'ERROR, CD PURGING DOES NOT ALLOW STD. GROUP') GO TO 500 500 RETURN END SUBROUTINE FIT(F,IC,IG,ID,IOFF) C FITS MODEL OF NO 3-FACTOR INTERACTION DOUBLE PRECISION F(100,100,10),FF(100,100,10),CG(100,100), * CGI(100,100),CD(100,10), * CDI(100,10),GD(100,10),GDI(100,10),X,XD,XL,XN,XP,Y DATA FF/100000*1.D0/ TOL=.001 IT=40 CALL XMCG(F,CG,IC,IG,ID) CALL XMCD(F,CD,IC,IG,ID) CALL XMGD(F,GD,IC,IG,ID) DO 50 M= 1,IT CALL XMCG(FF,CGI,IC,IG,ID) DO 10 I=1,IC DO 10 J=1,IG DO 10 K=1,ID 10 FF(I,J,K)=FF(I,J,K)*CG(I,J)/CGI(I,J) CALL XMCD(FF,CDI,IC,IG,ID) DO 20 I=1,IC DO 20 J=1,IG DO 20 K=1,ID 20 FF(I,J,K)=FF(I,J,K)*CD(I,K)/CDI(I,K) CALL XMGD(FF,GDI,IC,IG,ID) DO 30 I=1,IC DO 30 J=1,IG DO 30 K=1,ID 30 FF(I,J,K)=FF(I,J,K)*GD(J,K)/GDI(J,K) C CHECK CONVERGENCE CALL XMCG(FF,CGI,IC,IG,ID) X=0.0 DO 40 I=1,IC DO 40 J=1,IG Y=(CG(I,J) - CGI(I,J)) 40 X= X + DABS(Y) X= X/DFLOAT((IC*IG)) IF(X.LE.TOL)GO TO 51 50 CONTINUE 51 CONTINUE C CALCULATION OF FIT STATISTICS. PUT FITS IN F TO RETURN TO MAIN IF(IOFF.EQ.0)GO TO 60 XD=0.0 XP=0.0 XL=0.0 XN=0.0 DO 61 I=1,IC DO 61 J=1,IG DO 61 K=1,ID XN=XN+F(I,J,K) XD= XD + DABS(F(I,J,K) - FF(I,J,K)) XP= XP + ((F(I,J,K) - FF(I,J,K))**2)/FF(I,J,K) IF(F(I,J,K).LT..1)GO TO 61 XL=XL + F(I,J,K)*DLOG(F(I,J,K)/FF(I,J,K)) 61 CONTINUE IDF=(IC-1)*(IG-1)*(ID-1) XD=XD/(2.0*XN) XL=2.0*XL WRITE(6,70)XP,XL,XD,IDF,M,X 70 FORMAT(//,' FIT STATISTICS FOR MODEL OF NO 3-FACTOR INTERACTION 1',/,' PEARSON =',F10.2,/,' LIKELIHOOD-RATIO=',F10.2,/, 2' INDEX OF DISSIMILARITY=',F10.5,/,' DEGREES OF FREEDOM=',I10,/, 3' ITERATIONS =',I10,/,' DEVIATION = ',F10.7) 60 CONTINUE DO 62 I=1,IC DO 62 J=1,IG DO 62 K=1,ID 62 F(I,J,K)=FF(I,J,K) 99 RETURN END SUBROUTINE XMCG(F,CG,IC,IG,ID) C CALCULATES A 2-WAY MARGINAL FROM F AND RETURNS IN CG DOUBLE PRECISION F(100,100,10),CG(100,100),X DO 10 I=1,IC DO 10 J=1,IG X=0.0 DO 11 K=1,ID 11 X=X + F(I,J,K) CG(I,J)=X 10 CONTINUE RETURN END SUBROUTINE XMCD(F,CD,IC,IG,ID) DOUBLE PRECISION F(100,100,10),CD(100,10),X DO 10 I=1,IC DO 10 K=1,ID X=0.0 DO 11 J=1,IG 11 X=X + F(I,J,K) CD(I,K)=X 10 CONTINUE RETURN END SUBROUTINE XMGD(F,GD,IC,IG,ID) DOUBLE PRECISION F(100,100,10),GD(100,10),X DO 10 J=1,IG DO 10 K=1,ID X=0.0 DO 11 I=1,IC 11 X=X + F(I,J,K) GD(J,K)=X 10 CONTINUE RETURN END SUBROUTINE RESCAL (F,XG,IC,IG,ID) 00011450 C RESCALES PURGED F'S SO THAT GROUP TOTALS ARE PRESERVED 00011500 DOUBLE PRECISION F(100,100,10),XG(100),SUM 00011550 DO 10 J=1,IG 00011600 SUM=0.0 DO 11 I=1,IC 00011700 DO 11 K=1,ID 00011750 11 SUM= SUM + F(I,J,K) 00011800 DO 12 I=1,IC 00011850 DO 12 K=1,ID 00011900 12 F(I,J,K)=F(I,J,K)*(XG(J)/SUM) 00011950 10 CONTINUE 00012000 RETURN 00012050 END 00012100