C polyfac.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), " Polynomial factor analysis" in Proceedings of C the 77th Annual Convention of the American Psychological Association, C 4,103-104. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ *PFAC POLYNOMIAL FACTOR ANALYSIS * * INPUT * * 1) N,NF,M,MAXIT,NOPTY,NV,KPUNX,ID,KAL,NRAND,CRIT,CALF(9I4,I6,2F4.0) * N = NUMBER OF OBSERVATIONS. * NF = NUMBER OF FACTORS. * M = TOTAL NUMBER OF POLYNOMIAL COMPONENTS. * MAXIT = MAXIMUM NUMBER OF ITERATIONS. * NOPTY = CONTROL PARAMETER DEFINING INPUT MATRIX. * IF=0, SQUARE SYMMETRIC MATRIX C IS INPUT. * =1, NV X N MATRIX Y IS INPUT. * C=Y(TRANSPOSED)*Y, IS COMPUTED. * =2, Y IS INPUT, ROW MEANS ARE SUBTRACTED, * C=Y(TRANSPOSED)*Y. * =3, Y IS INPUT, COLUMN MEANS ARE SUBTRACTED AND * COLUMN VARIANCES STANDARDIZED, C=Y(TRANSPOSED)*Y. * =4, Y IS INPUT, COLUMN MEANS ARE SUBTRACTED, * C=Y(TRANSPOSED)*Y. * =5, Y IS INPUT, ROW MEANS ARE SUBTRACTED AND ROW * VARIANCES STANDARDIZED, C=Y(TRANSPOSED)*Y. * NV = NUMBER OF VARIABLES (IF NOPTY IS GREATER THAN 0, * OTHERWISE NOT RELEVANT). CANNOT BE GREATER THAN N. * KPUNX = PUNCH OPTION FOR MATRIX X. * IF=0, MATRIX X IS PUNCHED. * =1, MATRIX X IS NOT PUNCHED. * KAL = OPTION FOR ALPHA * IF=0, ALPHA=CALF*(4.0**(COS(IT)**3))/SQRT(GRL(IT)). * IF=1, ALPHA=(A-B)/C. * ID = IDENTIFICATION OPTION. * IF=0, NUMBERS ARE PLOTTED. * IF=1, READ IN IDENTIFICATION. * NRAND = RANDOM NUMBER OPTION. * IF=0, READ IN X MATRIX OF COORDINATES. * OTHERWISE = SIX DIGIT ODD INTEGER, USED TO GENERATE * RANDOM INITIAL X MATRIX. * CRIT = CRITERION FOR STOPPING. * IF=0, PROGRAM ASSIGNS .01 TO CRIT. * CALF = FLOATING POINT CONSTANT, USED TO CALCULATE ALPHA. * PRESENT ONLY IF KAL = 0 OR BLANK. * * 2) TITLE (12A6) * * 3) FORMAT FOR MATRIX C. (12A6) * * 4) MATRIX C * N X N MATRIX OF CROSS PRODUCTS. * Y(TRANSPOSED)*Y (INPUT ONLY IF NOPTY=0, COMPUTED OTHERWISE). * * 5) FORMAT FOR MATRIX Y. (12A6) * * 6) MATRIX Y. * NV X N MATRIX OF "SCORES" (RELEVANT ONLY IF NOPTY IS GREATER * THAN 0). * * 7) FORMAT FOR MATRIX X. (12A6) * * 8) MATRIX X. * NF X N MATRIX OF "FACTOR LOADINGS". * (INPUT ONLY IF NRAND = 0, OR BLANK, OTHERWISE RANDOM STARTING * CONFIGURATION IS GENERATED.) * * 9) FORMAT FOR MATRIX NE. (12A6) * * 10) MATRIX NE * M X NF MATRIX OF EXPONENTS FOR M POLYNOMIAL COMPONENTS WITH * RESPECT TO NF FACTORS. * * 11) FORMAT FOR IDENTIFICATION. * * 12) IDENTIFICATION ARRAY. * (INPUT ONLY IF ID = 1, OTHERWISE, NUMBERS ARE PLOTTED). * * REPEAT NUMBERS 1 THRU 12 FOR EACH ADDITIONAL SET OF DATA. * * FINAL CARD MUST BE BLANK. * * DIMENSION TITLE(13),FMTC(12),FMTX(12),FMTNE(12) DIMENSION VAF(101),GRL(101),ALF(101),COS(101),NBACK(101),X(5,100), 1DX(5,100),Z(20,100),DZ(20,100),NE(20,5),C(100,100),R(20,20),U(20,1 200),Y(50,100) DATA TITLE (1)/6H75H / CALL DIMENS (R,U) 111 READ 501,N,NF,M,MAXIT,NOPTY,NV,KPUNX,ID,KAL,NRAND,CRIT,CALF 501 FORMAT (9I4,I6,2F4.0) IF(N.EQ.0) GO TO 999 IF(NRAND.NE.0) PRINT 600, NRAND 600 FORMAT (9H1NRAND = ,I6) IF(CRIT.EQ.0.) CRIT=.01 READ 504,(TITLE(K),K=2,13) 504 FORMAT (12A6) READ 504,(FMTC(I),I=1,12) IF(NOPTY.GT.0) GO TO 503 * * READ MATRIX C. (IF NOPTY=0). * DO 510 I=1,N 510 READ FMTC,(C(I,J),J=1,N) GO TO 502 * * READ MATRIX Y. (IF NOPTY IS GREATER THAN 0). * 503 DO 520 I=1,NV 520 READ FMTC,(C(I,J),J=1,N) PRINT 705 705 FORMAT (17H0MATRIX Y (INPUT)) DO 710 I=1,NV 710 PRINT 521,I,(C(I,J),J=1,N) 502 IF(NRAND.EQ.0) GO TO 505 * * GENERATE RANDOM INITIAL X MATRIX. * CALL RANUM (X,N,NF,NRAND) DO 95 J=1,N SSQ=0.0 DO 90 I=1,NF 90 SSQ = SSQ+X(I,J)**2 95 X(I,J) = ABS(X(I,J)/SQRT(SSQ)) XN = N DO 80 I=1,NF SUM = 0. DO 85 J=1,N 85 SUM = SUM + X(I,J) XMEAN = SUM/XN DO 80 J=1,N 80 X(I,J) = X(I,J) - XMEAN GO TO 506 * * READ MATRIX X. * 505 READ 504,(FMTX(I),I=1,12) DO 530 I=1,NF 530 READ FMTX,(X(I,J),J=1,N) * * READ MATRIX NE. * 506 READ 504,(FMTNE(K),K=1,12) DO 540 I=1,M 540 READ FMTNE,(NE(I,J),J=1,NF) PRINT 543 543 FORMAT (19H0MATRIX NE (M X NF)) DO 545 I=1,M 545 PRINT 544,I,(NE(I,J),J=1,NF) 544 FORMAT (1H0,I2,3X,5I10) * * NORMALIZE MATRIX Y. * IF(NOPTY.EQ.0) GO TO 1100 IF(NOPTY.EQ.1) GO TO 1000 IF(NOPTY.EQ.3) GO TO 900 DO 820 J=1,NV YMEAN=0.0 DO 810 K=1,N 810 YMEAN=YMEAN+C(J,K) YMEAN=YMEAN/FLOAT(N) DO 820 K=1,N 820 C(J,K)=C(J,K)-YMEAN IF(NOPTY.EQ.2) GO TO 1000 DO 830 J=1,NV CNORM=0.0 DO 825 K=1,N 825 CNORM=CNORM+C(J,K)**2 CNORM=1./SQRT(CNORM) DO 830 K=1,N 830 C(J,K)=CNORM*C(J,K) GO TO 1000 900 DO 920 K=1,N YMEAN=0.0 DO 910 J=1,NV 910 YMEAN=YMEAN+C(J,K) YMEAN=YMEAN/FLOAT(NV) DO 920 J=1,NV 920 C(J,K)=C(J,K)-YMEAN IF(NOPTY.EQ.4) GO TO 1000 DO 940 K=1,N CNORM=0.0 DO 930 J=1,NV 930 CNORM=CNORM+C(J,K)**2 CNORM=1./SQRT(CNORM) DO 940 J=1,NV 940 C(J,K)=CNORM*C(J,K) 1000 PRINT 951 951 FORMAT (1H0,19HNORMALIZED Y MATRIX) DO 952I=1,NV 952 PRINT 521,I,(C(I,J),J=1,N) IF(NV.LE.50) NNV=NV IF(NV.GT.50) NNV=50 DO 950 K=1,N DO 950 J=1,NNV 950 Y(J,K)=C(J,K) * * COMPUTE MATRIX C. * DO 1002 K=1,N DO 1001 KK=K,N VAF(KK)=0.0 DO 1001 J=1,NV 1001 VAF(KK)=VAF(KK)+C(J,K)*C(J,KK) DO 1002 KK=K,N 1002 C(KK,K)=VAF(KK) DO 1003 K=1,N DO 1003 KK=K,N 1003 C(K,KK)=C(KK,K) 1100 CDSUM=0.0 DO 1101 K=1,N 1101 CDSUM=CDSUM+C(K,K) DO 1111 K=1,N DO 1111 KK=1,N 1111 C(K,KK)=C(K,KK)/CDSUM PRINT 715 715 FORMAT (9H0MATRIX C) DO 720 I=1,N 720 PRINT 521,I,(C(I,J),J=1,N) ENF=NF ALPHA=0.0 MAXIT1=MAXIT+1 DO 500 IT=1,MAXIT1 TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9999,TIME A ***** 9999 FORMAT(14H0A TIME = ,F8.5) VAF(IT)=0.0 GRL(IT)=0.0 ALF(IT)=0.0 COS(IT)=0.0 NBACK(IT)=0 ITM1=IT-1 * * TAKE STEP OF SIZE ALPHA, AND NORMALIZE X MATRIX. * 1 DO 6 J=1,NF CNORM=0.0 DO 5 K=1,N X(J,K) = X(J,K)+ALPHA*DX(J,K) 5 CNORM = CNORM+X(J,K)**2 CNORM = 1./SQRT(CNORM*ENF) DO 6 K=1,N 6 X(J,K) = CNORM*X(J,K) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9998,TIME B ***** 9998 FORMAT(14H0B TIME = ,F8.5) PRINT 430 430 FORMAT (9H0X MATRIX) DO 435 J=1,NF 435 PRINT 521,J,(X(J,K),K=1,N) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9997,TIME C ***** 9997 FORMAT(14H0C TIME = ,F8.5) * * DEFINE Z * DO 10 L=1,M DO 10 K=1,N Z(L,K) = 1.0 DO 10 J=1,NF IF(NE(L,J).EQ.0) GO TO 10 Z(L,K) = Z(L,K)*X(J,K)**NE(L,J) 10 CONTINUE TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9996,TIME D ***** 9996 FORMAT(14H0D TIME = ,F8.5) * * COMPUTE Z*Z(TRANSPOSED)=R * DO 20 L=1,M DO 20 LL=1,L R(L,LL)=0.0 DO 15 K=1,N 15 R(L,LL) = R(L,LL)+Z(L,K)*Z(LL,K) 20 R(LL,L)=R(L,LL) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9995,TIME E ***** 9995 FORMAT(14H0E TIME = ,F8.5) * * INVERT R=Z*Z(TRANSPOSED) * MM=1 CALL MATINV (R,M,MM,U) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9994,TIME F ***** 9994 FORMAT(14H0F TIME = ,F8.5) IF(MM.GT.1) GO TO 888 * * COMPUTE U=(Z*Z(TRANSPOSED)) INVERSED*Z * DO 30 L=1,M DO 30 K=1,N U(L,K)=0.0 DO 30 LL=1,M 30 U(L,K) = U(L,K) + R(L,LL)*Z(LL,K) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9993,TIME G ***** 9993 FORMAT(14H0G TIME = ,F8.5) * * COMPUTE DZ = 2.0((UC*U(TRANSPOSED))Z-UC) * DO 40 L=1,M DO 40 J=1,N DZ(L,J)=0.0 DO 40 I=1,N 40 DZ(L,J) = DZ(L,J)+U(L,I)*C(I,J) DO 50 L=1,M DO 50 LL=L,M R(L,LL) = 0.0 DO 45 J=1,N 45 R(L,LL) = R(L,LL)+DZ(L,J)*U(LL,J) 50 R(LL,L) = R(L,LL) DO 100 K=1,N DO 100 L=1,M A=0.0 DO 60 I=1,M 60 A = A+R(L,I)*Z(I,K) 100 DZ(L,K) = 2.0*(A-DZ(L,K)) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9992,TIME H ***** 9992 FORMAT(14H0H TIME = ,F8.5) * * COMPUTE TRACE U*C*Z(TRANSPOSED)=VAF(IT) * VAF(IT)=0.0 DO 35 L=1,M DO 35 K=1,N DO 35 KK=1,N 35 VAF(IT) = VAF(IT)+U(L,K)*C(K,KK)*Z(L,KK) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9991,TIME I ***** 9991 FORMAT(14H0I TIME = ,F8.5) PRINT 470,IT,VAF(IT) 470 FORMAT (5H0VAF(,I3,4H) = ,F10.5) * 38 ALF(IT) = ALPHA PRINT 475,ITM1,NBACK(ITM1) 475 FORMAT (7H0NBACK(,I2,4H) = ,I3) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9990,TIME J ***** 9990 FORMAT(14H0J TIME = ,F8.5) * * COMPUTE DX * 39 CRPR=0.0 GR=0.0 DO 200 J=1,NF DO 200 K=1,N DXC=0.0 DO 170 L=1,M 170 DXC = DXC+DZ(L,K)*FLOAT(NE(L,J))*Z(L,K)/X(J,K) CRPR = CRPR+DXC*DX(J,K) GR = GR+DXC**2 200 DX(J,K)=DXC TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9989,TIME K ***** 9989 FORMAT(14H0K TIME = ,F8.5) GRL(IT) = SQRT(GR) PRINT 485,IT,GRL(IT) 485 FORMAT (5H0GRL(,I3,4H) = ,F10.5) IF(IT.EQ.1) GO TO 299 COS(IT)=CRPR/(GRL(ITM1)*GRL(IT)) 299 IF(GRL(IT).LE.CRIT) GO TO 777 PRINT 490 490 FORMAT (10H0DX MATRIX) DO 495 J=1,NF 495 PRINT 521,J,(DX(J,K),K=1,N) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9988,TIME L ***** 9988 FORMAT(14H0L TIME = ,F8.5) * * BEGIN COMPUTATION OF "STEP-SIZE" ALPHA. * * COMPUTE DELTA Z (CALLED DZ). * IF (KAL.EQ.0) GO TO 265 IF(IT.EQ.1) GO TO 263 ALPHA = ALPHA*ABS(GRL(ITM1)/(GRL(ITM1)-GRL(IT)*COS(IT))) GO TO 500 263 ALPHA = CALF GO TO 500 265 DO 300 L=1,M DO 300 K=1,N DDZ=0.0 DO 270 J=1,NF 270 DDZ = DDZ+FLOAT(NE(L,J))*DX(J,K)/X(J,K) 300 DZ(L,K) = Z(L,K)*DDZ TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9987,TIME M ***** 9987 FORMAT(14H0M TIME = ,F8.5) * * COMPUTE * A = 1/2 TRACE(Z*DZ(TRANSPOSED)+DZ*Z(TRANSPOSED)) * *(U*C*U(TRANSPOSED)), * B = TRACE DZ*C*U(TRANSPOSED), AND * CC = TRACE (DZ*DZ(TRANSPOSED))*(U*C*U(TRANSPOSED)). * ALPHA = (A-B)/C * A=0.0 B=0.0 CC=0.0 DO 400 L=1,M DO 400 LL=1,L C1=0.0 C2=0.0 C3=0.0 C4=0.0 DO 350 K=1,N C1 = C1+Z(L,K)*DZ(LL,K)+DZ(L,K)*Z(LL,K) C2 = C2+DZ(L,K)*DZ(LL,K) DO 350 KK=1,N C3 = C3+U(L,K)*C(K,KK)*U(LL,KK) 350 C4 = C4+DZ(L,K)*C(K,KK)*U(LL,KK) IF(L.EQ.LL) GO TO 380 A=A+C1*C3 CC=CC+2.0*C2*C3 GO TO 400 380 A=A+C1*C3/2.0 B=B+C4 CC=CC+C2*C3 400 CONTINUE TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9986,TIME N ***** 9986 FORMAT(14H0N TIME = ,F8.5) ALPHA=(B-A)/CC 500 CONTINUE 5001 IF(NOPTY.EQ.0) GO TO 606 * * COMPUTE BHAT TRANSPOSED = U * Y (TRANSPOSED). * DO 601 J=1,M DO 601 K=1,NNV DZ(J,K)=0.0 DO 601 L=1,N 601 DZ(J,K) = DZ(J,K)+U(J,L)*Y(K,L) PRINT 603 603 FORMAT (1H0,58HMATRIX BHAT TRANSPOSED (MATRIX OF REGRESSION COEFFI 1CIENTS)) DO 610 I=1,M 610 PRINT 521,I,(DZ(I,J),J=1,NNV) * * COMPUTE YHAT = BHAT * Z. * DO 602 J=1,NNV DO 602 K=1,N Y(J,K)=0.0 DO 602 L=1,M 602 Y(J,K) = Y(J,K)+DZ(L,J)*Z(L,K) PRINT 604 604 FORMAT (1H0,11HMATRIX YHAT) DO 605 I=1,NNV 605 PRINT 521,I,(Y(I,J),J=1,N) 606 TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9985,TIME O ***** 9985 FORMAT(14H0O TIME = ,F8.5) * * PLOT "FINAL CONFIGURATION." * CALL PLOT (ID,N,NF,TITLE,X) * * OUTPUT AND FINIS. * 777 PRINT 519,(TITLE(K),K=2,13) 519 FORMAT (1H1,12A6) PRINT 511,ITM1 511 FORMAT (14H0 ITERATION = ,I3) PRINT 512 512 FORMAT (11H0 ITERATION,7X,3HVAF,7X,3HGRL,7X,3HCOS,7X,3HALF) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9984,TIME P ***** 9984 FORMAT(14H0P TIME = ,F8.5) DO 550 K=1,IT KK=K-1 550 PRINT 513,KK,VAF(K),GRL(K),COS(K),ALF(K) 513 FORMAT (1H ,4X,I3,3X,4F10.5) TIME = ITIMEZ(0) TIME = TIME/10. 0000 PRINT 9983,TIME Q ***** 9983 FORMAT(14H0Q TIME = ,F8.5) PRINT 514 514 FORMAT (56H0 MATRIX X (NUMBER OF FACTORS BY NUMBER OF OBSERVATIONS 1)) DO 560 I=1,NF 560 PRINT 521,I,(X(I,J),J=1,N) PRINT 516 516 FORMAT (57H0 MATRIX DX (NUMBER OF FACTORS BY NUMBER OF OBSERVATION 1S)) DO 570 I=1,NF 570 PRINT 521,I,(DX(I,J),J=1,N) * 521 FORMAT (1H ,I3,10E12.4/(1H ,3X,10E12.4)) * PRINT 517 517 FORMAT (59H0 MATRIX Z (NUMBER OF COMPONENTS BY NUMBER OF OBSERVATI 1ONS)) DO 580 I=1,M 580 PRINT 521,I,(Z(I,J),J=1,N) IF(KPUNX.EQ.1) GO TO 999 PUNCH 504,(TITLE(K),K=2,13) PUNCH 522 522 FORMAT (11H(2X,5E13.5)) DO 590 I=1,NF 590 PUNCH 518,I,(X(I,J),J=1,N) 518 FORMAT (I2,5E13.5/(2X,5E13.5)) GO TO 111 888 PRINT 889 889 FORMAT (29H0 R REMAINS UNINVERTED. STOP.) 999 STOP END