C ccm.f C The authors of this software are Jih-Jie Chang and J D Carroll. C Copyright (c) 1993 by AT&T. C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire C notice is included in all copies of any software which is or includes C a copy or modification of this software and in all copies of the C supporting documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP- C LIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C This software comes from the FIRST MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C Carroll,J.D. (1969), "Categorical conjoint measurement", unpublished C manuscript, Bell Laboratories, Murray Hill, NJ. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTRAN $ FORTRAN C CATEGORICAL CONJOINT MEASUREMENT C JULY, 1969, CHANG FOR CARROLL DIMENSION A(2000),V(2000),P(2000),ANCL(2000) DIMENSION TIT(12),NWW(4) C NWAY - NUMBER OF WAYS (OR FACTORS). MAX.=4 C NC - NUMBER OF CATEGORIES. MAX.=20 C NL - NUMBER OF SOLUTIONS. MAX .LE. NC C NWW - NUMBER OF LEVELS WITHIN FACTORS. THE PRODUCT OF LEVELS C MUST BE .LE. 2000. C FLGEOF IS SYSTEM ROUTINE TO TEST EOF, IF SO SET NEND=1. NEND=0 CALL FLGEOF(5,NEND) 50 READ (5,100) NWAY,NC,NL IF(NEND.EQ.1) GOTO999 READ (5,101) TIT DO1I=1,4 1 NWW(I)=1 READ (5,100) (NWW(I),I=1,NWAY) N1=NWW(1) N2=NWW(2) N3=NWW(3) N4=NWW(4) CALL READA(A,N1,N2,N3,N4,NC,V) C FIND MAXIMUM OF NWW MNW=0 ND=0 DO2I=1,NWAY ND=ND&NWW(I) IF(NWW(I).LE.MNW) GOTO2 MNW=NWW(I) 2 CONTINUE PRINT102,TIT PRINT103,NWAY,NC PRINT104,(NWW(I),I=1,NWAY) CALL CCM(A,N1,N2,N3,N4,NWAY,NWW,NC,NL,V,MNW,ANCL,P,ND) GOTO50 100 FORMAT(20I4) 101 FORMAT(12A6) 102 FORMAT(33H1CATEGORICAL CONJOINT MEASUREMENT/1H012A6) 103 FORMAT(6H0NWAY=,I2/14H CATEGORIES - ,I3) 104 FORMAT(10H LEVELS - ,4I3) 999 STOP END $ FORTRAN SUBROUTINE READA(A,N1,N2,N3,N4,NC,V) INTEGER A DIMENSION A(N1,N2,N3,N4),SYM(20),FMT(12),V(N1,N2,N3,N4) READ (5,100) FMT 100 FORMAT(12A6) READ (5,FMT) (SYM(I),I=1,NC) DO5I4=1,N4 DO5I3=1,N3 DO5I2=1,N2 5 READ (5,FMT) (V(I1,I2,I3,I4),I1=1,N1) C CHANGE FROM SYMBOLS TO NUMERALS DO6I4=1,N4 DO6I3=1,N3 DO6I2=1,N2 DO6I1=1,N1 DO7M=1,NC IF(V(I1,I2,I3,I4).EQ.SYM(M)) GOTO8 7 CONTINUE PRINT101,I1,I2,I3,I4 101 FORMAT(6H0CELL(,4I3,35H) SYMBOL DOES NOT MATCH, CHECK DATA) STOP 8 A(I1,I2,I3,I4)=M 6 CONTINUE RETURN END $ FORTRAN SUBROUTINE CCM, R=P*PT, DOES ANOVA SUBROUTINE CCM(A,N1,N2,N3,N4,NWAY,NWW,NC,NL,V,MNW,ANCL,P,ND) DIMENSION A(N1,N2,N3,N4),V(N1,N2,N3,N4),NWW(1),CCT(20) DIMENSION ANCL(MNW,4,NC),P(ND,NC),R(20,20),B(20,20) DIMENSION AM(50,4,10),XNK(4),XNWW(4) DIMENSION SS(4,10),SSE(10),XMS(4,10),EMS(10),F(4,10),NDF(4) C NCT(I) - NUMBER OF OCCURENCES OF CATEGORY I. C ANCL(I,J,K) - NUMBER OF OCCURENCES OF CAT. K FOR FACTOR J LEV. I. INTEGER A NT=1 DO16I=1,NWAY NDF(I)=NWW(I)-1 XNWW(I)=NWW(I) 16 NT=NT*NWW(I) NDFE=NT-1-ND&NWAY EDF=NDFE XNT=NT DO17I=1,NWAY 17 XNK(I)=NT/NWW(I) PRINT113 DO30I4=1,N4 DO30I3=1,N3 DO30I1=1,N1 30 PRINT114,I1,(A(I1,I2,I3,I4),I2=1,N2) C COMPUTE NUMBER OF OCCURENCES CALL ZEROT(ANCL,MNW,4,NC) DO3I=1,NC 3 CCT(I)=0. DO5I4=1,N4 DO5I3=1,N3 DO5I2=1,N2 DO5I1=1,N1 M=A(I1,I2,I3,I4) CCT(M)=CCT(M)&1. GOTO(6,7,8,9),NWAY 9 ANCL(I4,4,M)=ANCL(I4,4,M)&1. 8 ANCL(I3,3,M)=ANCL(I3,3,M)&1. 7 ANCL(I2,2,M)=ANCL(I2,2,M)&1. 6 ANCL(I1,1,M)=ANCL(I1,1,M)&1. 5 CONTINUE PRINT101,(I,I=1,NC) PRINT102,( CCT(I),I=1,NC) PRINT103 DO1I=1,NWAY L=NWW(I) DO1J=1,L 1 PRINT104,I,J,(ANCL(J,I,M),M=1,NC) DO10M=1,NC L=0 DO10I=1,NWAY TEMP=SQRT(1./(XNK(I)*CCT(M))) TEM=CCT(M)/XNWW(I) K=NWW(I) DO10J=1,K L=L&1 10 P(L,M)=TEMP*(ANCL(J,I,M)-TEM) L=0 DO20I=1,NWAY 20 XNWW(I)=XNWW(I)-1. C COMPUTE R DO11I=1,NC DO11J=1,NC R(I,J)=0. DO11M=1,ND 11 R(I,J)=R(I,J)&P(M,I)*P(M,J) MM=1 CALL MTEVV(R,20,20,B,20,20,NC,MM) PRINT108,(R(I,I),I=1,NC) PRINT109 DO14I=1,NC CCT(I)=SQRT(CCT(I)) 14 PRINT107,I,(B(J,I),J=1,NC) DO42I=1,NC DO42J=1,NC 42 B(J,I)=B(J,I)/CCT(J) PRINT122 DO45I=1,NC 45 PRINT107,I,(B(J,I),J=1,NC) CALL ZEROT(AM,50,4,10) DO13L=1,NL DO12I4=1,N4 DO12I3=1,N3 DO12I2=1,N2 DO12I1=1,N1 M=A(I1,I2,I3,I4) 12 V(I1,I2,I3,I4)=B(M,L) C PUNCH REAL VALUED FACTORS PUNCH115,L PUNCH118 DO15I4=1,N4 DO15I3=1,N3 DO15I1=1,N1 15 PUNCH117,I1,(V(I1,I2,I3,I4),I2=1,N2) C COMPUTE MARGINAL MEANS DO21I4=1,N4 DO21I3=1,N3 DO21I2=1,N2 DO21I1=1,N1 Z=V(I1,I2,I3,I4) GOTO(22,23,24,25),NWAY 25 AM(I4,4,L)=AM(I4,4,L)&Z 24 AM(I3,3,L)=AM(I3,3,L)&Z 23 AM(I2,2,L)=AM(I2,2,L)&Z 22 AM(I1,1,L)=AM(I1,1,L)&Z 21 CONTINUE SSQ=0. DO18I=1,NWAY M=NWW(I) SS(I,L)=0. DO36J=1,M AM(J,I,L)=AM(J,I,L)/XNK(I) 36 SS(I,L)=SS(I,L)&AM(J,I,L)**2 SS(I,L)=SS(I,L)*XNK(I) XMS(I,L)=SS(I,L)/XNWW(I) 18 SSQ=SSQ&SS(I,L) SSE(L)=1.-SSQ EMS(L)=SSE(L)/EDF DO37I=1,NWAY 37 F(I,L)=XMS(I,L)/EMS(L) DO32I4=1,N4 DO32I3=1,N3 DO32I2=1,N2 DO32I1=1,N1 32 V(I1,I2,I3,I4)=AM(I1,1,L)&AM(I2,2,L)&AM(I3,3,L)&AM(I4,4,L) PRINT125,L DO33I4=1,N4 DO33I3=1,N3 DO33I1=1,N1 33 PRINT107,I1,(V(I1,I2,I3,I4),I2=1,N2) 13 CONTINUE PRINT110 DO19I=1,NWAY PRINT111,I PRINT112,(L,L=1,NL) M=NWW(I) DO19J=1,M 19 PRINT116,J,(AM(J,I,L),L=1,NL) PRINT119 DO38L=1,NL PRINT120,L PRINT121 DO40I=1,NWAY 40 PRINT123,I,SS(I,L),NDF(I),XMS(I,L),F(I,L) 38 PRINT124,SSE(L),NDFE,EMS(L) 101 FORMAT(38H0NUMBER OF OCCURENCES OF EACH CATEGORY/(10I6)/) 102 FORMAT(10F6.1) 103 FORMAT(43H0OCCURRENCES OF EACH CATEGORY ON EACH LEVEL/) 104 FORMAT(2I2,3X,10E12.4/(7X,10E12.4)) 107 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) 108 FORMAT(12H0EIGEN ROOTS/(2X,10E12.4)) 109 FORMAT(14H0EIGEN VECTORS) 110 FORMAT(15H1MARGINAL MEANS) 111 FORMAT(7H0FACTOR,I3) 112 FORMAT(2X,5HLEVEL,5X,8HSOLUTION,I2,I10,6I15) 113 FORMAT(5H0DATA) 114 FORMAT(I4,2X,40I3/(6X,40I3)) 115 FORMAT(36H REAL VALUED FACTORS - FROM SOLUTION,I2) 116 FORMAT(I5,1X,8E15.5) 117 FORMAT(I3,2X,5E13.5/(5X,5E13.5)) 118 FORMAT(11H(5X,5E13.5)) 119 FORMAT(21H1ANALYSIS OF VARIANCE) 120 FORMAT(9H0SOLUTIONI3) 121 FORMAT(4H0WAY,8X,3HSSQ,10X,2HDF,11X,2HMS,15X,1HF) 122 FORMAT(16H0CATEGORY VALUES) 123 FORMAT(I3,3X,E15.5,I6,2E18.5) 124 FORMAT(4H RES,2X,E15.5,I6,3X,E15.5) 125 FORMAT(26H0ADDITIVE TABLE - SOLUTION,I2) RETURN END $ FORTRAN SUBROUTINE ZEROS(A,N) DIMENSION A(1) K=N 10 DO5II=1,K 5 A(II)=0. GOTO50 ENTRY ZEROD(A,I,J) K=I*J GOTO10 ENTRY ZEROT(A,I,J,M) K=I*J*M GOTO10 50 RETURN END $ FORTRAN $ INCODE IBMF CMTINV MTINV SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MTPK1115 DIMENSION A(I1,I2), B(I3,I4) MTPK1116 C MTPK1117 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MTPK1118 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MTPK1119 C PRECISION ARITHMETIC MTPK1120 C MTPK1121 DOUBLE PRECISION B,BMAX,BMULT MTPK1122 EQUIVALENCE (BMAX,BMULT) MTPK1123 C MTPK1124 C SINGLE TO DOUBLE PRECISION TRANSFER MTPK1125 C MTPK1126 101 DO 103 I = 1,N MTPK1127 102 DO 103 J = 1,N MTPK1128 103 B(I,J) = A(I,J) MTPK1129 C MTPK1130 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. MTPK1133 C MTPK1134 200 DO 430 K = 1,N MTPK1131 A(K,1) = FLOAT(K) MTPK1135 A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) MTPK1137 210 DO 219 I = K,N MTPK1138 211 DO 219 J = K,N MTPK1139 IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) MTPK1141 214 A(K,1) = FLOAT(I) MTPK1142 215 A(K,2) = FLOAT(J) MTPK1143 219 CONTINUE MTPK1145 216 IF (BMAX - 1.D-38) 801,301,301 C MTPK1146 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MTPK1147 C MTPK1148 301 I = IFIX(A(K,1)) MTPK1149 302 J = IFIX(A(K,2)) MTPK1150 313 DO 314 M = 1,N MTPK1151 BMULT = B(I,M) MTPK1152 B(I,M) = B(K,M) MTPK1153 314 B(K,M) = BMULT MTPK1155 303 DO 304 M = 1,N MTPK1156 BMULT = B(M,J) MTPK1157 B(M,J) = B(M,K) MTPK1158 304 B(M,K) = BMULT MTPK1159 C MTPK1160 C ACTUAL INVERSION STEP. MTPK1161 C BEGIN ROW ITERATION MTPK1162 C MTPK1163 401 DO 425 I = 1,N MTPK1164 C IGNORE ROW K. MTPK1165 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. MTPK1167 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. MTPK1169 410 DO 420 J = 1,N MTPK1170 C IGNORE COLUMN K MTPK1171 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT MTPK1173 413 GO TO 420 MTPK1174 414 B(I,J) = B(I,J) - B(K,J) * BMULT MTPK1175 420 CONTINUE MTPK1176 425 CONTINUE MTPK1177 C MTPK1178 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MTPK1179 C MTPK1180 426 BMULT = 1.D0/ B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION MTPK1183 428 DO 429 J = 1,N MTPK1184 429 B(K,J) = B(K,J) * BMULT MTPK1185 430 CONTINUE MTPK1186 C MTPK1187 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MTPK1188 C MTPK1189 501 DO 512 IJK = 1,N MTPK1190 502 I = N-IJK&1 MTPK1191 503 L = IFIX(A(I,2)) MTPK1192 504 J = IFIX(A(I,1)) MTPK1193 505 DO 508 M = 1,N MTPK1194 506 BMULT = B(I,M) MTPK1195 507 B(I,M) = B(L,M) MTPK1196 508 B(L,M) = BMULT MTPK1197 509 DO 512 M = 1,N MTPK1198 510 BMULT = B(M,I) MTPK1199 511 B(M,I) = B(M,J) MTPK1200 512 B(M,J) = BMULT MTPK1201 C MTPK1202 C MOVE A INVERSE BACK INTO A MTPK1203 C MTPK1204 601 DO 603 I = 1,N MTPK1205 602 DO 603 J = 1,N MTPK1206 603 A(I,J) = B(I,J) 701 M = 1 MTPK1208 702 RETURN MTPK1209 C MTPK1210 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C MTPK1212 801 M = 2 MTPK1213 802 RETURN MTPK1214 999 END MTPK1215 $ FORTRAN $ INCODE IBMF *EIGENJ EIGENVALUE AND EIGENVECTOR SUBROUTINE EIGJ0020 * CD600D4.009 DATE 05/05/65 EIGJ0030 SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4) DIMENSION A(5050),V( I3, I4),E(NN) LOGICAL IV EIGJ0060 * EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF EIGJ0070 * A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD EIGJ0080 * ONLY THE UPPER TRIANGLE IS STORED EIGJ0090 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0100 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0110 N=NN EIGJ0120 IND=1 EIGJ0130 NM2=N-2 EIGJ0140 AMAX=0.0 EIGJ0150 NM=N-1 EIGJ0160 IF(.NOT.IV)GO TO 5 EIGJ0170 * SET UP VECTOR MATRIX EIGJ0180 1 DO 3 I=1,N EIGJ0190 DO 2 J=1,N EIGJ0200 2 V(I,J)=0.0 EIGJ0210 3 V(I,I)=1.0 EIGJ0220 5 K=2 EIGJ0230 DO 28 I=1,NM EIGJ0240 IP=I&1 EIGJ0250 DO 27 J=IP,N EIGJ0260 Y=A(K) EIGJ0270 X=ABS(Y) EIGJ0280 IF(X-AMAX)27,27,6 EIGJ0290 6 GO TO (7,8),IND EIGJ0300 7 AMAX=X EIGJ0310 GO TO 27 EIGJ0320 * TRANSFORMATION SETUP FOLLOWS EIGJ0330 8 II=K&I-J EIGJ0340 JJ=(J*(N&N-J&3))/2-N EIGJ0350 ITEST=1 EIGJ0360 X=A(II)-A(JJ) EIGJ0370 IF(X)10,11,10 EIGJ0380 10 X=Y/X EIGJ0390 Y=.5*X EIGJ0400 IF(ABS(Y)-.414213562)14,11,11 EIGJ0410 11 C=.707106781 EIGJ0420 IF(Y)12,12,13 EIGJ0430 12 S=-C EIGJ0440 GO TO 15 EIGJ0450 13 S=C EIGJ0460 GO TO 15 EIGJ0470 14 X=Y*Y EIGJ0480 C=1.&X EIGJ0490 S=2.*Y/C EIGJ0500 C=(1.-X)/C EIGJ0510 15 X=S*S EIGJ0520 Y=C*C EIGJ0530 XY=S*C EIGJ0540 AXY=2.*A(K)*XY EIGJ0550 Z=A(II)*Y&AXY&A(JJ)*X EIGJ0560 W=A(II)*X-AXY&A(JJ)*Y EIGJ0570 A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II)) EIGJ0580 A(II)=Z EIGJ0590 A(JJ)=W EIGJ0600 IF(NM2)24,24,16 EIGJ0610 16 IT=I EIGJ0620 JT=J EIGJ0630 * TRANSFORM A EIGJ0640 DO 23 M=1,NM2 EIGJ0650 NP=N-M EIGJ0660 IF(JT-K)20,17,18 EIGJ0670 17 IT=IT&1 EIGJ0680 JT=JT&NP EIGJ0690 ITEST=2 EIGJ0700 18 IF(IT-K)20,19,19 EIGJ0710 19 IT=IT&1 EIGJ0720 JT=JT&1 EIGJ0730 ITEST=3 EIGJ0740 20 X=A(IT)*C&A(JT)*S EIGJ0750 A(JT)=A(JT)*C-A(IT)*S EIGJ0760 A(IT)=X EIGJ0770 GO TO (21,22,23),ITEST EIGJ0780 21 IT=IT&NP EIGJ0790 JT=JT&NP EIGJ0800 GO TO 23 EIGJ0810 22 IT=IT&1 EIGJ0820 JT=JT&NP-1 EIGJ0830 23 CONTINUE EIGJ0840 24 IF(.NOT.IV)GO TO 27 EIGJ0850 * TRANSFORM V EIGJ0860 25 DO 26 M=1,N EIGJ0870 X=V(M,I)*C&V(M,J)*S EIGJ0880 V(M,J)=V(M,J)*C-V(M,I)*S EIGJ0890 26 V(M,I)=X EIGJ0900 27 K=K&1 EIGJ0910 28 K=K&1 EIGJ0920 * SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP EIGJ0930 GO TO (29,31),IND EIGJ0940 29 IF(AMAX)30,35,30 EIGJ0950 30 AT=AMAX EIGJ0960 IND=2 EIGJ0970 F=.25 EIGJ0980 AMAX=F*AT EIGJ0990 MAX=6 EIGJ1000 GO TO 5 EIGJ1010 31 MAX=MAX-1 EIGJ1020 IF(MAX)34,33,32 EIGJ1030 32 F=2.0*F*F EIGJ1040 AMAX=AT*F EIGJ1050 33 C=0.0 EIGJ1060 GO TO 5 EIGJ1070 34 IF(C)33,35,33 EIGJ1080 * COPY EIGENVALUES INTO E EIGJ1090 35 MM=1 EIGJ1100 NM=N&1 EIGJ1110 DO 36 M=1,N EIGJ1120 E(M)=A(MM) EIGJ1130 36 MM=MM&NM-M EIGJ1140 RETURN EIGJ1150 END $ FORTRAN $ INCODE IBMF SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M) DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100) K=1 DO1J=1,N DO1I=J,N A(K)=X(I,J) 1 K=K&1 CALL EIGENJ(A,D,N,M,E,I3,I4) CALL SORTFL(-1) IF(M)4,3,4 3 CALL SORT(E(1),E(2),N) DO 50 I=1,N 50 X(I,I)=E(I) RETURN 4 DO 5 I=1,N MS(I)=I 5 NS(I)=I CALL SORT(E(1),E(2),N,MS(1)) CALL SORTFX(1) CALL SORT(MS(1),MS(2),N,NS(1)) CALL ARREV(N,NS,D,I3,I4) DO 6 I=1,N 6 X(I,I)=E(I) RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE ARREV(N,NS,D,I3,I4) DIMENSION NS(N),D(I3,I4) DO 10 I=1,N 1 IF(NS(I)-I)2,10,2 2 CALL XCH(N,I,NS(I),D,NS,I3,I4) GOTO1 10 CONTINUE RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE XCH(N,L,M,D,NS,I3,I4) DIMENSION NS(N),D(I3,I4) I=L J=M DO 1 K=1,N T=D(K,I) D(K,I)=D(K,J) 1 D(K,J)=T NS(I)=NS(J) NS(J)=J RETURN END