C*************************************************************** C*** * C*** MAIN PROGRAM OVERALS * C*** CANONICAL ANALYSIS OF NS SETS * C*** INPUT AND OUTPUT OF THE PARAMETERS AND DYNAMICAL * C*** ARRAY ALLOCATION * C*** * C*** INPUT FROM UNIT IK, AUXILLIARY DATA SET ON UNIT IS * C*** AND OUTPUT TO UNIT IO * C*** IK,IS AND IO ARE FIXED TO RESPECTIVILY 5,12 AND 6 * C*** * C*** SUBROUTINES CALLED: DYNAM * C*** * C*************************************************************** C C OVERALS80 was prepared from OVERALS V1.0 by Stef van Buuren C at the Dept. of Psychology, UCLA, Los Angeles. C April 1988 C C*************************************************************** PROGRAM OVERALS EXTERNAL INDEXO CHARACTER*80 ITIT, IFMT, IF1, IF2, IF3, IF4, IF5 DIMENSION ITEM(16) IK=5 IO=6 IS=10 OPEN(UNIT = 10, STATUS = 'SCRATCH') WRITE(IO,9109) 9109 FORMAT('we will first build the parameter file interactively',//) WRITE(IO,9110) 9110 FORMAT('a title for the job') READ (IK,2) ITIT WRITE(IS,2) ITIT WRITE(IO,9111) 9111 FORMAT('the number of observations') READ (IK,*) N WRITE(IO,9112) 9112 FORMAT('the number of sets') READ (IK,*) NS WRITE(IO,9113) 9113 FORMAT('the number of passive variables') READ(IK,*) MM WRITE(IS,3)N,NS,MM WRITE(IO,9114) 9114 FORMAT('the number of dimensions') READ(IK,*) NP WRITE(IO,9115) 9115 FORMAT('the maximum number of iterations') READ(IK,*) IR WRITE(IO,9116) 9116 FORMAT('the tolerance for stress differences',/ + ,'default is 5 (i.e. 10**-5) and',/ 1 ,'maximum precision is 6 (i.e. 10**-6)') READ(IK,*) IEP WRITE(IO,9117) 9117 FORMAT('the initial configuration',/ + ,'zero is random, one is rational') READ(IK,*) INI WRITE(IS,3)NP,IR,IEP,INI WRITE(IO,9118) 9118 FORMAT('the number of active variables per set. input is',/ + ,'list directed but use exactly 16 numbers per record') NT=NS M=0 DO 200 K=1,NS,16 READ(IK,*)ITEM WRITE(IS,3)ITEM MIN=MIN0(NT,16) DO 150 J=1,MIN 150 M=M+ITEM(J) NT=NT-16 200 CONTINUE WRITE(IO,9119) 9119 FORMAT('print options: data',/ + ,'zero is ten rows, one is whole matrix') READ(IK,*)I1 WRITE(IO,9200) 9200 FORMAT('print options: object scores',/ + ,'zero is no, one is yes') READ(IK,*)I2 WRITE(IO,9201) 9201 FORMAT('print options: weights and loadings',/ + ,'zero is no, one is yes') READ(IK,*)I3 WRITE(IO,9202) 9202 FORMAT('print options: category quantifications',/ + ,'zero is no, one is yes') READ(IK,*)I4 WRITE(IO,9203) 9203 FORMAT('print options: loss partitioning',/ + ,'zero is no, one is yes') READ(IK,*)I5 WRITE(IS,3)I1,I2,I3,I4,I5 i1=0 WRITE(IO,9204) 9204 FORMAT('plot options: object scores',/ + ,'zero is no, one is yes') READ(IK,*)I2 WRITE(IO,9205) 9205 FORMAT('plot options: component loadings',/ + ,'zero is no, one is yes') READ(IK,*)I3 WRITE(IO,9206) 9206 FORMAT('plot options: quantifications',/ + ,'zero is no, one is yes') READ(IK,*)I4 WRITE(IS,3)I1,I2,I3,I4 WRITE(IO,9207) 9207 FORMAT('i/o options: name input file',/ + 'this must be an existing ascii file') READ(IK,2)IF1 OPEN(UNIT = 11, FILE = IF1, STATUS = 'OLD') WRITE(IO,9208) 9208 FORMAT('i/o options: name score output file') READ(IK,2)IF2 OPEN(UNIT = 12, FILE = IF2, STATUS = 'NEW') WRITE(IO,9209) 9209 FORMAT('i/o options: name loading output file') READ(IK,2)IF3 OPEN(UNIT = 13, FILE = IF3, STATUS = 'NEW') WRITE(IO,9210) 9210 FORMAT('i/o options: name quantification output file') READ(IK,2)IF4 OPEN(UNIT = 14, FILE = IF4, STATUS = 'NEW') WRITE(IO,9211) 9211 FORMAT('i/o options: name loss output file') READ(IK,2)IF5 OPEN(UNIT = 15, FILE = IF5, STATUS = 'NEW') WRITE(IS,3)11,12,13,14,15 WRITE(IO,9120) 9120 FORMAT('the number of active variables per set. input is',/ + ,'list directed but use exactly 16 numbers per record') IC=0 MC=0 MTOT=M+MM MMA=M MMB=MTOT DO 300 K=1,MTOT,16 READ(IK,*)ITEM WRITE(IS,3)ITEM MIN=MIN0(MMA,16) MIN2=MIN0(MMB,16) DO 250 J=1,MIN 250 MC=MC+ITEM(J) DO 255 J=1,MIN2 255 IC=MAX0(IC,ITEM(J)) MMA=MMA-16 MMB=MMB-16 300 CONTINUE WRITE(IO,9121) 9121 FORMAT('measurement level per active variable. input is',/ + ,'list directed but use exactly 16 numbers per record') DO 400 J=1,M,16 READ(IK,*)ITEM 400 WRITE(IS,3)ITEM WRITE(IO,9122) 9122 FORMAT('the number of active variables per set. input is',/ + ,'list directed but use exactly 16 numbers per record') DO 500 J=1,M,16 READ(IK,*)ITEM 500 WRITE(IS,3)ITEM WRITE(IO,9123) 9123 FORMAT('input format for variables') READ(IK,2)IFMT WRITE(IS,2)IFMT 2 FORMAT(A80) 3 FORMAT(16I5) REWIND IS IQ=N*MTOT+3*NP*MC+5*N*NP+2*MC+4*M+2*NS+6*IC+2*NP*NP+2*NP*M+NP 1+2*N+M*NP*NP+MTOT CALL DYNAM(INDEXO,A,IQ,N,M,MTOT,NP,MC,IC,NS,IS,IO) STOP END !!S DYNAM SUBROUTINE DYNAM(INDEXO,B,IQ,N,M,MTOT,NP,MC,IC,NS,IV,IO) C************************************************* C*** * C*** ALLOCATION OF FIXED STORAGE INSTEAD OF * C*** THE DYNAMIC ALLOCATION BY THE ASSEMBLER * C*** ROUTINE DYNAM * C*** * C************************************************* C*** INTEGER IQ,N,M,MMM,MM,NP,MC,IC,NS,IV,IO REAL A,B C-----ALLOCATION OF STORAGE IF DYNAM IS NOT LINKED DIMENSION A(50000) C*** ATTENTION, IQ HAS TO BE COMPARED WITH THE SIZE OF ARRAY A IF(IQ.LE.50000)GOTO 100 WRITE(IO,10) 10 FORMAT('Not enough storage for this problem') STOP 100 CALL INDEXO(A,IQ,N,M,MTOT,NP,MC,IC,NS,IV,IO) RETURN END !!S INDEXO SUBROUTINE INDEXO(A,IQ,N,M,MTOT,NP,MC,IC,NS,IS,IO) C************************************************* C*** * C*** COMPUTATION OF INDEX NUMBERS FOR * C*** MEMORY ALLOCATION * C*** * C*** SUBROUTINES CALLED: OVRALS * C*** * C************************************************* REAL A DIMENSION A(IQ) I1=1 I2=I1+N*MTOT I3=I2+MC I4=I3+NP*MC I5=I4+NP*MC I6=I5+NP*M I7=I6+MC I8=I7+N*NP I9=I8+N*NP I10=I9+N*NP I11=I10+N*NP I12=I11+M I13=I12+NP I14=I13+NP*NP I15=I14+NP*NP*M I16=I15+MTOT I17=I16+M I18=I17+NS I19=I18+NS I20=I19+M I21=I20+IC I22=I21+IC I23=I22+NP*NP I24=I23+M I25=I24+IC*2 I26=I25+IC*2 I27=I26+MC*NP I28=I27+NP*M I29=I28+N I30=I29+N I31=I30+N*NP C CALL OVRALS(A(I1),A(I2),A(I3),A(I4),A(I5),A(I6),A(I7),A(I8), 1 A(I9),A(I9),A(I10),A(I11),A(I12),A(I13),A(I14),A(I15), 2 A(I16),A(I17),A(I18),A(I19),A(I20),A(I20),A(I21), 3 A(I22),A(I23),A(I24),A(I25),A(I26),A(I27),A(I28), 4 A(I29),A(I30),N,M,MTOT,NP,MC,IC,NS,IS,IO,MC*NP,IC*2) RETURN END !!S MAINSUB CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C IN INDICATOR MATRIX IN REDUCED FORM (N,M) C C ID VECTOR OF MARGINAL FREQUENCIES (MC) C C EE MATRIX OF CATEGORY SCORES AFTER LAST ITERATION (NP,MC) C C V SUM OF INDIVIDU SCORES PER GROUP (N,NP) C C W SUM OF INDIVIDU SCORES OVER ALL GROUPS (N,NP) C C X NEW INDIVIDU SCORES C C H (N,NP) AND H7(IC,NP) WORKARRAYS C C E THE MULTIPLE QUANTIFICATIONS (NP,MC) C C A MATRIX OF WEIGHTS (M,NP) C C Z VECTOR OF CATEGORY SCORES FOR SINGLE VARIABLES (MC) C C NC NUMBER OF CATEGORIES PER VARIABLE (M) C C KC CUMULATIVE NUMBER OF CATEGORIES PER VARIABLE (M) C C NV NUMBER OF VARIABLES PER SET (NS) C C KV CUMULATIVE NUMBER OF VARIABLES PER SET (NS) C C IT TYPE NUMBER PER VARIABLE (M) C C HZ(IC), HS(NP,NP), HD(IC), HH(IC*2), HR(NP,NP); ALL WORK ARRAY C C H2(NP), H3(IC*2), H4(IC*2), H5(MCNP), HP(NP); ALL WORKARRAY C C DI TOTAL DISPERSION (NP,NP,M) C C CO COMPONENT LOADINGS (NP,M) C C MK 1=NON-MISSING DATA FOR OBJECT(I), 0=MISSING (N) C C MS MEAN OF THE SUM OVR MK(I) (N) C C N NUMBER OF INDIVIDUALS C C M NUMBER OF ACTIVE VARIABLES C C MM NUMBER OF PASSIVE VARIABLES C C MTOT TOTAL NUMBER OF VARIABLES C C NP NUMBER OF DIMENSIONS C C MC TOTAL NUMBER OF CATEGORIES C C IC MAXIMUM NUMBER OF CATEGORIES C C NS NUMBER OF SETS C C IR MAXIMUM NUMBER OF ITERATIONS C C EP CRITERIUM FOR STRESS DIFFERENCE C C IV UNIT NUMBER FOR INPUT PARAMETERS, INTERNALLY FIXED ON 5 C C IO UNIT NUMBER FOR OUTPUT ON PRINTER,INTERNALLY FIXED ON 6 C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C IMPORTANT|||| WORKARRAY H7(IC,NP) OVERLAPS W(N,NP),H(N,NP) C C WORKARRAY HH(IC*2) OVERLAPS HZ(IC),HD(IC) C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE OVRALS(IN,ID,E,EE,A,Z,X,V,H7,W,H,H2,HP,HR,DI,NC,KC,NV, 1 KV,IT,HH,HZ,HD,HS,IW,H3,H4,H5,CO,MK,MS,H8, 2 N,M,MTOT,NP,MC,IC,NS,IS,IO,MCNP,IC2) C********************************************************************** C*** * C*** OVERALS CANONICAL ANALYSIS OF TWO OR MORE SETS * C*** * C*** SUBROUTINES CALLED: INPUDA, MARGIN, RANDMA, DEVIA2, ININUM, * C*** PROIN2, PLUSM2, MINUS2, DISPER, INDPRO, DIAPRO, DOUBLE,STRESS * C*** PRINAX, OUTPAR, PASSIV, OUTUNI * C*** * C********************************************************************** REAL E,V,W,X,H,H2,H7,A,Z,HZ,HD,H3,H4,EE,MS,H12 DOUBLE PRECISION S1,S2,H11,H22,STR INTEGER IN,ID,NC,KC,NV,KV,IT,N,M,NP,MC,IC,NS,IS,IO,IR,MCNP,IC2,INI LOGICAL LAST DIMENSION IN(N,MTOT),ID(MC),E(MCNP),EE(MCNP),V(N,NP),W(N,NP), 1 X(N,NP),H7(IC,NP),H(N,NP),H2(M),A(NP,M),Z(MC),NC(MTOT),KC(M), 2 NV(NS),KV(NS),IT(M),HH(IC2),HZ(IC),HD(IC),HS(NP,NP),IW(M), 3 H3(IC2),H4(IC2),H5(MCNP),CO(NP,M),HP(NP),HR(NP,NP), 4 DI(NP,NP,M),MK(N),MS(N),H8(N,NP) COMMON IX(5),IY(5),IZ(5) C - INPUT PARAMETERS AND DATA MATRIx CALL INPUDA(IN,NV,KV,NC,KC,IT,IW,N,M,MTOT,NS,IR,IS,IO,EP,INI,MM) C - IA=0 AT LEAST ONE VARIABLE SINGLE, IB=1 ALL VARS MULTIPLE OR NUMERICAL IA=1 IB=1 DO 120 J=1,M IF(IT(J).NE.3)IA=0 120 IF(IT(J).LE.1)IB=0 EEP=AMAX1(EP,1.E-4) IG=1 IF(INI.EQ.0.OR.IB.EQ.1)IG=0 C - COMPUTE THE MARGINAL FREQUENCIES CALL MARGIN(IN,ID,MK,MS,H2,NV,KV,NC,KC,N,M,MC,NS,IC,IO,F,IOPT) C - RANDOM START VALUES FOR MATRIX W CALL RANDMA(W,N*NP) DO 45 I=1,MCNP 45 E(I)=0.0 C - S1 STRESS PREVIOUS UPDATE, S2 STRESS CURRENT UPDATE WRITE(IO,50) 50 FORMAT(///3H1 ,16HIteration Stress,7X,3HFit/) CALL ININUM(IN,ID,E,A,Z,X,V,W,H,HP, + HR,HS,NC,KC,NV,KV,IT,MK,MS,HZ,HD,N,M,NP,MC,IC,NS, + IS,IO,MCNP,IC2,EEP,S2,IG,F,IOPT) LAST=.FALSE. IF(IR.EQ.1)GOTO 850 IF(IR.EQ.2)LAST=.TRUE. II=0 100 S1=S2 S2=0.0 C - II ITERATION NUMBER II=II+1 C - X = THE DEVIATIONS FROM THE COLUMN MEANS OF W CALL DEVIA2(W,X,MS,N,NP,F,IOPT) C - PROCRUSTERS ORTOGONALIZATION OF MATRIX X CALL PROCRU(MS,X,X,H,HP,HR,HS,N,NP,IO,IOPT) C - STANDARDISATION ON N*I DO 190 L=1,NP DO 190 I=1,N 190 X(I,L)=X(I,L)*SQRT(FLOAT(N)) C - START WITH ZERO VALUES FOR W FOR EVERY ITERATION DO 195 J=1,NP DO 195 I=1,N 195 W(I,J)=0.0 C+++++++++++++++++++ LOOP OVER SETS +++++++++++++++++++++++++++++++ DO 700 K=1,NS DO 200 I=1,N 200 MK(I)=1 DO 205 J=1,NP DO 205 I=1,N 205 V(I,J)=0.0 J1=KV(K)-NV(K)+1 J2=KV(K) DO 220 I=1,N DO 210 J=J1,J2 IF(IN(I,J).LE.NC(J)) GOTO 210 MK(I)=0 GOTO 220 210 CONTINUE 220 CONTINUE C - LOOP OVER VARIABLES PER SET DO 300 J=J1,J2 C - V(K) = V(K) + IN(J) * E(J) SUM OVER J OF SET K CALL PROIN2(IN(1,J),E((KC(J)-NC(J))*NP+1),H,N,NC(J),NP,IOPT) CALL PLUSM2(V,H,V,MK,N,NP,IOPT) 300 CONTINUE J1=KV(K)-NV(K)+1 J2=KV(K) C******************** LOOP OVER VARS PER SET ************************* DO 500 J=J1,J2 I1=KC(J)-NC(J)+1 I2=(KC(J)-NC(J))*NP+1 C - V(K) = V(K) - IN(J) * E(J) CALL PROIN2 (IN(1,J),E(I2),H,N,NC(J),NP,IOPT) CALL MINUS2(V,H,V,MK,N,NP,IOPT) C - E(J) = (1/D) * IN(TRANSPOSED) * (X-V(K)) CALL MINUS2(X,V,H,MK,N,NP,IOPT) C - COMPUTE TOTAL DISPERSION MATRIX IF (LAST) CALL DISPER(H,H8,DI(1,1,J),HS,N,NP,J) CALL INDPRO(IN(1,J),H,E(I2),N,NC(J),NP) CALL DIAPRO(ID(I1),E(I2),E(I2),NC(J),NP) C - THE MULTIPLE CATEGORY QUANTIFICATIONS ARE STORED IN E AND EE IF (LAST) CALL DOUBLE(E(I2),EE(I2),NC(J),NP) IF(IT(J).EQ.3) GOTO 450 C - UPDATE E,A AND Z FOR SINGLE VARIABLES CALL SINGLE(E(I2),A(1,J),Z(I1),ID(I1),HZ,HD,NC(J),NP,N,IT(J)) C - V(K) = V(K) + IN(J) * E(J) 450 CALL PROIN2(IN(1,J),E(I2),H,N,NC(J),NP,IOPT) CALL PLUSM2(V,H,V,MK,N,NP,IOPT) 500 CONTINUE C*********************** END LOOP OVER VARS *************************** C - W = W + V(K) DO 460 L=1,NP DO 460 I=1,N 460 W(I,L)=W(I,L)+V(I,L) C - STRESS: S2=S2(-TR(2*X*V(K))+TR(V(K)*V(K))+N*NP)/NP*NS*N C SUM OVER K CALL STRESS(X,V,N,NP,NS,STR) S2=S2+STR 700 CONTINUE C+++++++++++++++++++++++++ END LOOP OVER SETS ++++++++++++++++++++++++ C - OUTPUT OF STRESS AND FIT H11=S2*FLOAT(NS) H22=S1*FLOAT(NS) H12=FLOAT(NP)*(1.-S2) WRITE(IO,750)II,H11,H12 750 FORMAT(3H ,I5,2F12.7) C - TEST ON MAXIMUM NUMBER OF ITERATIONS IF(LAST) GOTO 850 IF(II.GE.(IR-1)) GOTO 800 C - TEST ON STRESS DIFFERENCE IF(DABS(H22-H11).GT.EP) GOTO 100 C - A LAST TIME FOR COMPUTING AXES ROTATION 800 LAST=.TRUE. GOTO 100 C - STRESS PER DIMENSION 850 IF(IG.EQ.1)WRITE(IO,51) 51 FORMAT(///3H ,40HThe first stress is the stress of a solu, + 'tion with all single variables'/' treated as single', + ' numerical (with stress difference .0001 and maximum'/ + 29H number of iterations is 50.)) C - ROTATION OF OBJECT SCORES, MULTIPLE CATEGORY QUANTIFICATIONS C - AND ,FOR SINGLE VARIABLES, COMPONENT LOADINGS CALL PRINAX(X,E,EE,A,IT,HP,HR,HS,H,KC,NC,N,NP,MCNP,M,NS,IO) C - OUTPUT OF THE ESTIMATED PARAMETERS CALL OUTPAR(IN,X,E,EE,A,Z,NC,KC,IT,ID,H,H2,HH,HZ,HD,H3,H4,H5, + H7,HS,CO,MK,HP,HR,DI,M,N,NP,MCNP,MC,IO,IC,IC2,IW,IA, + IOPT) C - OUTPUT FOR PASSIVE VARIABLES IF(MM.NE.0) + CALL PASSIV(IN,X,H7,H,HP,HS,HZ,NP,NC,IC,N,IO,MM,MTOT,IC*NP) CALL OUTUNI(IO) RETURN END !!S COMPO SUBROUTINE COMPO(E,X,ID,H,HP,HS,N,NP,NC,J,IO) C********************************************************* C*** COMPONENT LOADING-MATRICES FOR MULTIPLE VARIABLES * C********************************************************* COMMON IX(5),IY(5),IZ(5) REAL E,X,H,HP,HS INTEGER ID,N,NP,J DIMENSION E(NC,NP),X(N,NP),ID(NC),H(N,NP),HP(NP),HS(NP,NP) C - COMPUTE OMEGA (=DISCRIMINATION MEASURES) DO 7 L=1,NP HP(L)=0.0 DO 5 K=1,NC 5 HP(L)=HP(L)+E(K,L)*ID(K)*E(K,L) 7 HP(L)=SQRT(HP(L)/FLOAT(N)) DO 30 L2=1,NP DO 30 L1=1,NP HS(L1,L2)=0.0 DO 30 I=1,N 30 HS(L1,L2)=HS(L1,L2)+X(I,L1)*H(I,L2) DO 40 L1=1,NP DO 40 L2=1,NP 40 HS(L1,L2)=HS(L1,L2)/(FLOAT(N)*HP(L2)) IF(IX(3).EQ.0)GOTO 50 DO 45 L2=1,NP 45 WRITE(IO,906)L2,(HS(L1,L2),L1=1,NP) 50 IF(IZ(3).EQ.0)GOTO 60 IU=IZ(3) DO 47 L2=1,NP 47 WRITE(IU,907)J,L2,(HS(L1,L2),L1=1,NP) 906 FORMAT(3H ,I5,10F7.2/(8X,10F7.2)) 907 FORMAT(2I5,12F8.3,(12X,12F8.3)) 60 RETURN END !!S DEVIA2 SUBROUTINE DEVIA2(W,X,MS,N,NP,F,IOPT) C************************************************* C*** * C*** X = DEVIATIONS FROM THE COLUMN MEANS OF W * C*** * C************************************************* REAL W,X,S,F,MS DIMENSION W(N,NP),X(N,NP),MS(N) GOTO (5,40), IOPT 5 DO 30 L=1,NP S=0.0 DO 10 I=1,N 10 S=S+W(I,L) S=S/F DO 20 I=1,N 20 X(I,L)=W(I,L)-(S*MS(I)) 30 CONTINUE RETURN 40 DO 70 L=1,NP S=0.0 DO 50 I=1,N 50 S=S+W(I,L) S=S/FLOAT(N) DO 60 I=1,N 60 X(I,L)=W(I,L)-S 70 CONTINUE RETURN END !!S DIAPRO SUBROUTINE DIAPRO(ID,A,B,NC,NP) C************************************************* C*** * C*** B = (ID**-1) * A * C*** ID DIAGONAL MATRIX (NC,NC) GIVEN AS * C*** VECTOR (NC) * C*** * C************************************************* REAL A,B INTEGER ID DIMENSION A(NC,NP),B(NC,NP),ID(NC) DO 10 I=1,NC DO 10 J=1,NP IF(ID(I).EQ.0) GOTO 10 B(I,J)=A(I,J)/FLOAT(ID(I)) 10 CONTINUE RETURN END !!S DISPER SUBROUTINE DISPER(H,H8,DI,HS,N,NP,J) C************************************************* C*** DISPERSION = (X-V(K))'(X-V(K))/N * C*** AFTER ROTATION WITH THE EIGENVECTORS HS * C************************************************* REAL H,DI INTEGER N,NP DIMENSION H(N,NP),DI(NP,NP),H8(N,NP),HS(NP,NP) DO 10 L=1,NP DO 10 I=1,N H8(I,L)=0.0 DO 10 K=1,NP 10 H8(I,L)=H8(I,L)+H(I,K)*HS(K,L) DO 30 L2=1,NP DO 30 L1=1,NP DI(L1,L2)=0.0 DO 20 I=1,N 20 DI(L1,L2)=DI(L1,L2)+H8(I,L1)*H8(I,L2) 30 DI(L1,L2)=DI(L1,L2)/FLOAT(N) RETURN END !!S DOUBLE SUBROUTINE DOUBLE(E,EE,NC,NP) C************************************************* C*** EE=E * C************************************************* REAL E,EE INTEGER NC,NP DIMENSION E(NC,NP),EE(NC,NP) DO 10 L=1,NP DO 10 I=1,NC 10 EE(I,L)=E(I,L) RETURN END !!S IMTQL2 SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C ***************************************************************** C * * C * I M T Q L 2 * C * * C * PURPOSE: DETERMINES THE EIGENVALUES AND EIGENVECTORS OF A * C * SYMMETRIC TRIDIAGONAL MATRIX USING THE IMPLICIT * C * QL METHOD. * C * THE EIGENVECTORS ARE STORED IN Z AND * C * THE EIGENVALUES ARE STORED IN D * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK-ORIGINAL VERSION NOT IN DOUBLE PRECISION* C * * C ***************************************************************** DIMENSION D(N),E(N),Z(NM,N) REAL D,E,Z,B,C,F,G,P,R,S,MACHEP REAL SQRT,ABS,SIGN C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPCIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** C MACHEP = 2.0**(-20) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0 C DO 240 L = 1, N J = 0 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 IF (ABS(E(M)) .LE. MACHEP* (ABS(D(M)) + ABS(D(M+1)))) X GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ.30) GO TO 1000 J = J + 1 C ********** FORM SHIFT ********** G = (D(L+1) - P) / (2.0 * E(L)) R = SQRT(G*G+1.0) G = D(M) - P + E(L) / (G + SIGN(R,G)) S = 1.0 C = 1.0 P = 0.0 NML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 200 II = 1, NML I = M - II F = S * E(I) B = C * E(I) IF (ABS(F) .LT. ABS(G)) GO TO 150 C = G / F R = SQRT(C*C+1.0) E(I+1) = F * R S = 1.0 / R C = C * S GO TO 160 150 S = F / G R = SQRT(S*S+1.0) E(I+1) = G* R C = 1.0 / R S = S * C 160 G = D(I+1) - P R = (D(I) - G) * S + 2.0 * C * B P = S * R D(I+1) = G + P G = C * R - B C ********** FORM VECTOR ********** DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0 GO TO 105 240 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P ) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 1000 IERR = L 1001 RETURN END !!S INDPRO SUBROUTINE INDPRO(IN,H,E,N,NC,NP) C************************************************* C*** * C*** E = (IN TRANSPOSED) * H * C*** IN INDICATOR MATRIX (N,NC) GIVEN * C*** IN REDUCED FORM, IE VECTOR (N) * C*** WITH CATEGORY SCORES * C*** * C************************************************* INTEGER IN REAL H,E DIMENSION IN(N),H(N,NP),E(NC,NP) DO 10 K=1,NP DO 10 J=1,NC E(J,K)=0.0 DO 10 I=1,N 10 IF(IN(I).EQ.J) E(J,K)=E(J,K)+H(I,K) RETURN END !!S ININUM SUBROUTINE ININUM(IN,ID,E,A,Z,X,V,W,H,HP,HR,HS,NC,KC,NV,KV,IT, + MK,MS,HZ,HD,N,M,NP,MC,IC,NS,IS,IO,MCNP,IC2, + EEP,S2,IG,F,IOPT) C ******************************************************************* C * COMPUTING INITIAL CONFIGURATION WHERE ALL SINGLE VARIABLES ARE * C * TREATED AS S I N G L E N U M E R I C A L AND ALL MULTIPLE * C * VARIABLES M U L T I P L E N O M I N A L * C ******************************************************************* LOGICAL ITER DOUBLE PRECISION S1,S2,STR DIMENSION IN(N,M),ID(MC),E(MCNP),V(N,NP),W(N,NP),X(N,NP), 1 H(N,NP),A(NP,M),Z(MC),NC(M),KC(M),NV(NS), 2 KV(NS),IT(M),HZ(IC),HD(IC),HS(NP,NP), 3 HP(NP),HR(NP,NP),MK(N),MS(N) S2=10000 II=0 ITER=.FALSE. IF(IG.EQ.0)ITER=.TRUE. 100 S1=S2 S2=0.0 C - II ITERATION NUMBER II=II+1 C - X = THE DEVIATIONS FROM THE COLUMN MEANS OF W CALL DEVIA2(W,X,MS,N,NP,F,IOPT) C - PROCRUSTERS ORTOGONALIZATION OF MATRIX X CALL PROCRU(MS,X,X,H,HP,HR,HS,N,NP,IO,IOPT) C - STANDARDISATION ON N*I DO 190 L=1,NP DO 190 I=1,N 190 X(I,L)=X(I,L)*SQRT(FLOAT(N)) C - START WITH ZERO VALUES FOR W FOR EVERY ITERATION DO 195 J=1,NP DO 195 I=1,N 195 W(I,J)=0.0 C+++++++++++++++++++ LOOP OVER SETS +++++++++++++++++++++++++++++++ DO 700 K=1,NS DO 200 I=1,N 200 MK(I)=1 DO 205 J=1,NP DO 205 I=1,N 205 V(I,J)=0.0 J1=KV(K)-NV(K)+1 J2=KV(K) DO 220 I=1,N DO 210 J=J1,J2 IF(IN(I,J).LE.NC(J)) GOTO 210 MK(I)=0 GOTO 220 210 CONTINUE 220 CONTINUE C - LOOP OVER VARIABLES PER SET DO 300 J=J1,J2 C - V(K) = V(K) + IN(J) * E(J) SUM OVER J OF SET K CALL PROIN2(IN(1,J),E((KC(J)-NC(J))*NP+1),H,N,NC(J),NP,IOPT) CALL PLUSM2(V,H,V,MK,N,NP,IOPT) 300 CONTINUE J1=KV(K)-NV(K)+1 J2=KV(K) C******************** LOOP OVER VARS PER SET ************************* DO 500 J=J1,J2 I1=KC(J)-NC(J)+1 I2=(KC(J)-NC(J))*NP+1 C - V(K) = V(K) - IN(J) * E(J) CALL PROIN2 (IN(1,J),E(I2),H,N,NC(J),NP,IOPT) CALL MINUS2(V,H,V,MK,N,NP,IOPT) C - E(J) = (1/D) * IN(TRANSPOSED) * (X-V(K)) CALL MINUS2(X,V,H,MK,N,NP,IOPT) CALL INDPRO(IN(1,J),H,E(I2),N,NC(J),NP) CALL DIAPRO(ID(I1),E(I2),E(I2),NC(J),NP) IF(IT(J).EQ.3) GOTO 450 C - UPDATE E,A AND Z FOR SINGLE VARIABLES IF(II.EQ.1) CALL SINGIN(E(I2),A(1,J),Z(I1),ID(I1),NC(J),NP,N,IO) IF(II.GT.1) + CALL SINGLE(E(I2),A(1,J),Z(I1),ID(I1),HZ,HD,NC(J),NP,N,2) C - V(K) = V(K) + IN(J) * E(J) 450 CALL PROIN2(IN(1,J),E(I2),H,N,NC(J),NP,IOPT) CALL PLUSM2(V,H,V,MK,N,NP,IOPT) 500 CONTINUE C*********************** END LOOP OVER VARS *************************** C - W = W + V(K) DO 460 L=1,NP DO 460 I=1,N 460 W(I,L)=W(I,L)+V(I,L) C - STRESS: S2=S2(-TR(2*X*V(K))+TR(V(K)*V(K))+N*NP)/NP*NS*N C SUM OVER K CALL STRESS(X,V,N,NP,NS,STR) S2=S2+STR 700 CONTINUE C+++++++++++++++++++++++++ END LOOP OVER SETS ++++++++++++++++++++++++ C - OUTPUT OF STRESS AND FIT C - TEST ON MAXIMUM NUMBER OF ITERATIONS IF(ITER) GOTO 850 IF(II.GE.49) GOTO 800 C - TEST ON STRESS DIFFERENCE IF(DABS(S1-S2).GT.EEP) GOTO 100 800 ITER=.TRUE. GOTO 100 850 RETURN END !!S INPUDA SUBROUTINE INPUDA(IN,NV,KV,NC,KC,IT,IW,N,M,MTOT,NS,IR,IS,IO,EP, + INI,MM) C********************************************************************* C*** * C*** INPUT OF PARAMETERS AND DATA MATRIX * C*** * C*** SUBROUTINES CALLED: PRINDA * C*** * C********************************************************************* REAL EP CHARACTER*80 FT,TEXT INTEGER IN,NV,KV,NC,KC,IT,N,M,NS,IR,IS,IO,IX,IY,IZ,IW,IP, + IEP,INI,MM,MTOT DIMENSION IN(N,MTOT),NV(NS),KV(NS),NC(MTOT),KC(M),IT(M), + IW(M) COMMON IX(5),IY(5),IZ(5) C C - INPUT PARAMETERS FROM UNIT IS READ(IS,3)TEXT READ(IS,1)N,NS,MM READ(IS,1)NP,IR,IEP,INI READ(IS,1)NV READ(IS,1)IX READ(IS,1)IY READ(IS,1)IZ READ(IS,1)NC READ(IS,1)IT READ(IS,1)IW READ(IS,3)FT C - CHECK ON MEASUREMENT LEVEL DO 10 J=1,M 10 IF(IT(J).GT.3) GOTO 50 WRITE(IO,16) WRITE(IO,17) C - COMPUTATION CUMULATIVE VARIABLES KV(1)=NV(1) DO 60 K=2,NS 60 KV(K)=KV(K-1)+NV(K) KC(1)=NC(1) DO 70 J=2,M 70 KC(J)=KC(J-1)+NC(J) C C - COMPUTATION OF MAXIMUM RANK WITHIN EACH SET OF DATA MATRIX IRANK=0 DO 30 J=1,M IF(IT(J).EQ.3) GOTO 20 IRANK=IRANK+1 GOTO 30 20 NCJ=NC(J) IRANK=IRANK+NCJ-1 30 CONTINUE NN=N-1 IRAK=MIN0(IRANK,NN) IF(NP.GT.IRAK)GOTO 75 C C - INPUT DATA MATRIX FROM UNIT IP IP=IZ(1) DO 80 I=1,N 80 READ(IP,FT)(IN(I,J),J=1,MTOT) C C - REWRITE AND DEFAULT MINIMUM STRESS DIFFERENCE IIP=IEP IF(IEP.GT.6.OR.IEP.LE.0)IEP=5 EP=1.0000/10**IEP C - OUTPUT PARAMETERS TO PRINTER WRITE(IO,11)TEXT WRITE(IO,504)N,NS,MM WRITE(IO,505)NP,IR,EP,INI WRITE(IO,506)NV WRITE(IO,507)FT WRITE(IO,508)IX(1),IZ(1),(IX(I),IY(I),IZ(I),I=2,4),IX(5),IZ(5) WRITE(IO,509)(I,NC(I),IT(I),IW(I),I=1,M) IF(MM.EQ.0) GOTO 105 M1=M+1 WRITE(IO,510) DO 90 I=M1,MTOT 90 WRITE(IO,511)I,NC(I) C C - PRINT DATA MATRIX 105 CALL PRINDA(IN,N,MTOT,IX(1),IO) C - ECHO OF PARAMETER CARDS WRITE(IO,12)N,NS,MM,NP,IR,IIP,INI WRITE(IO,13)NV WRITE(IO,13)IX WRITE(IO,15)(IY(J),J=2,5) WRITE(IO,13)IZ WRITE(IO,13)NC WRITE(IO,13)IT WRITE(IO,13)IW C DO 100 J=1,MTOT DO 100 I=1,N 100 IF(IN(I,J).LE.0.OR.IN(I,J).GT.NC(J)) IN(I,J)=NC(J)+1 C 1 FORMAT(16I5) 3 FORMAT(A80) 504 FORMAT( /37H Number of objects ,I5 + //37H Number of sets ,I5 + //37H Number of passive variables ,I5) 505 FORMAT( /37H Number of dimensions ,I5 + //37H Maximum number of iterations ,I5 + //36H Minimum stress difference ,F12.6 + //37H Initial configuration ,I5) 508 FORMAT(1H1 + ///1H ,38X,15HPrint Plot Unit + //37H Data matrix ,I5,4X,1H-,I5, + //37H Object scores ,3I5 + //37H Weights & Component loadings ,3I5 + //37H Multiple & Single & Projected ,3I5 + /37H Categories & Centroids , + //37H Partitioning loss ,I5,4X,1H-,I5) 506 FORMAT( /37H Number of active variables for ,8I5, + /37H each set ,8I5/(37X,8I5)) 507 FORMAT(//30H Input format ,/ + 3(1H ,A80/)) 509 FORMAT(//////// + ' Information for active variables:'// + ' Cat = Number of categories'/ + ' Type = Variable type 0 = Single nominal'/ + ' 1 = Single ordinal'/ + ' 2 = Single numerical'/ + ' 3 = Multiple nominal'/ + ' Plot = Plot specifications (O = no, 1 = yes)', + ///,11X,11H Variable,22H Cat Type Plot/, + (/,10X,3H ,I5,6X,3I6)) 510 FORMAT(// + ' Number of categories for passive variables'/) 511 FORMAT(/,10X,3H ,I5,6X,I6) 11 FORMAT(////////////////// + ' Title'/1H ,A80) 12 FORMAT(/// + ' Echo of parameter cards (without text and format)', + //1H ,3I5//1H ,4I5) 13 FORMAT(/1H ,16I5) 15 FORMAT(/1H ,5X,4I5) 16 FORMAT (1H1,//////, +10X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++/ +10X,58H + +/ +10X,58H + April 1988 O V E R A L S 8 0 Version 1.0 +/ +10X,58H + +/ +10X,58H + OVERALS80 is an adapted version of OVERALS V1.0 in */ +10X,58H + which all output is rescaled to 80 columns. Use the */ +10X,58H + OVERALS program to produce smooth printer output. */ +10X,58H + +/ +10X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++) 17 FORMAT(////////// + ' O V E R A L S 8 0 Canonical analysis of two or ', + 'more Renee Verdegaal '/ + ' sets for discrete data with m', + 'ixed Eeke van der Burg '/ + ' Version 1.0 measurement levels ', + ' Jan de Leeuw '/ + 1H ,59X,'Dept. of Data Theory'/ + ' April 1988 ', + ' University of Leiden'/ + 1H ,59X,'The Netherlands'////) RETURN 50 WRITE(IO,19)J 19 FORMAT(///46H Measurement level not correct for variable,I3) STOP 75 WRITE(IO,29)IRAK 29 FORMAT(///' The number of dimensions is greater than the'/ + ' maximum rank of the data matrix.'/ + ' The maximum rank is:',I5/ + ' Adjust the number of dimensions and resubmit', + ' the job.'//) STOP END !!S LIDICA SUBROUTINE LIDICA(MARFRE,YBLANK,EPLUS,DPLUS,KJ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C*** C C*** THIS ROUTINE COMPUTES A WEIGTHED LINEAR REGRESSION OF THE C C*** MODEL VECTOR YBLANK ON CATEGORICAL DISCRETE DATA C C*** C C*** NO SUBROUTINES CALLED J.VAN RIJCKEVORSEL MAY 78 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION MARFRE(KJ),YBLANK(KJ),EPLUS(KJ),DPLUS(KJ) REAL CONST,DELTA,GAMMA,HULP1,HULP2,ALPHA,BETA,DPLUS,EPLUS,YBLANK INTEGER MARFRE C*** CONST=0. DELTA=0. GAMMA=0. ALPHA=0. BETA =0. C*** DO 1 K=1,KJ CONST=CONST+MARFRE(K) DELTA=DELTA+MARFRE(K)*K GAMMA=GAMMA+MARFRE(K)*YBLANK(K) 1 CONTINUE C*** IF (CONST.GT.1.E-20)CONST=1./CONST DELTA=DELTA*CONST GAMMA=GAMMA*CONST HULP1=0. HULP2=0. C*** DO 2 K=1,KJ EPLUS(K)=YBLANK(K)-GAMMA DPLUS(K)=K-DELTA HULP1=HULP1+EPLUS(K)*MARFRE(K)*DPLUS(K) HULP2=HULP2+DPLUS(K)*MARFRE(K)*DPLUS(K) 2 CONTINUE C*** IF (HULP2.GT.1.E-20) ALPHA=HULP1/HULP2 IF(ALPHA.LT.0.) ALPHA=-ALPHA BETA=GAMMA-ALPHA*DELTA C*** DO 3 K=1,KJ IF(MARFRE(K).EQ.0) GOTO 3 YBLANK(K)=ALPHA*K+BETA 3 CONTINUE RETURN END !!S MARGIN SUBROUTINE MARGIN(IN,ID,MK,MS,IMIS,NV,KV,NC,KC,N,M,MC,NS,IC,IO, + F,IOPT) C**************************************************************** C*** * C*** COMPUTE MARGINAL FREQUENCIES (ID) OF MATRIX IN * C*** ID=IN(TRANSPOSED) * MK * IN * C*** COMPUTE M*=(SUM OF MK)/NUMBER OF SETS * C*** COMPUTE NUMBER OF MISSING PER VARIABLE (IMIS) * C*** * C**************************************************************** INTEGER IN,ID,NV,KV,NC,KC,N,M,MC,NS,IC,IO,IM,IMIS REAL F,MS DIMENSION IN(N,M),ID(MC),NV(NS),KV(NS),NC(M),KC(M),IMIS(M), + MK(N),MS(N) DATA AA/4H----/ DATA BB/4H I/ DO 10 L=1,MC 10 ID(L)=0 DO 15 I=1,N 15 MS(I)=0.0 C - COMPUTE MS(=M*) AND MK FOR EACH OBJECT DO 22 K=1,NS KK=K DO 17 I=1,N 17 MK(I)=1 J1=KV(K)-NV(K)+1 J2=KV(K) DO 21 I=1,N DO 20 J=J1,J2 IF(IN(I,J).LE.NC(J)) GOTO 20 MK(I)=0 GOTO 21 20 CONTINUE 21 MS(I)=MS(I)+MK(I) DO 19 J=J1,J2 IMIS(J)=0 DO 19 I=1,N IF(IN(I,J).LE.NC(J)) GOTO 7 IMIS(J)=IMIS(J)+1 GOTO 19 7 IR=KC(J)-NC(J)+IN(I,J) ID(IR)=ID(IR)+(1*MK(I)) 19 IF(MK(I).EQ.0) IN(I,J)=NC(J)+1 22 CONTINUE F=0 C - IOPT=1 : THERE ARE MISSING DATA, IOPT=2 :NO MISSING DATA IOPT=1 IG=0 DO 28 I=1,N IF(MS(I).EQ.0)WRITE(IO,11)I IF(MS(I).EQ.0)IG=1 MS(I)=MS(I)/NS 28 F=F+MS(I) IF(IG.EQ.1) GOTO 40 IF(F.EQ.N) IOPT=2 WRITE(IO,33) WRITE(IO,24)(J,J=1,IC) MIN=MIN0(IC,10) WRITE(IO,26)(AA,J=1,MIN) DO 30 J=1,M J1=KC(J)-NC(J)+1 J2=KC(J) ISOMID=0 DO 35 K=J1,J2 35 ISOMID=ISOMID+ID(K) NN=N-ISOMID WRITE(IO,29)BB,BB 30 WRITE(IO,27)J,BB,IMIS(J),NN,BB,(ID(L),L=J1,J2) 29 FORMAT(1H ,7X,A4,13X,A4) 33 FORMAT(/////// + ' Marginal frequencies left for analysis'/ + ' --------------------------------------'/ + ' Mis var = Number of missing per variable'/ + ' Mis set = Number of missing within the set where the ', + 'variable belongs to'// + /17X,4H *,/19X,4H *,/17X,26H * Categories, + /9X,18H Mis Mis *) 24 FORMAT(10X,19H var set *,10I5,(/29X,10I5)) 26 FORMAT(30H Variables ------------------,10(1X,A4)) 27 FORMAT(3H ,I5,A4,I5,3X,I5,A4,10I5/(29X,10I5)) 11 FORMAT(/14H Object number,5X,I5) 41 FORMAT(//' This object(s) must be removed from the analysis'/ + ' while this/these observation(s) has/have missing'/ + ' values in all sets.'/////) RETURN 40 WRITE(IO,41) STOP END !!S MINUS2 SUBROUTINE MINUS2(A,B,C,MK,N,NP,IOPT) C************************************************* C*** * C*** C = A - B * C*** * C************************************************* REAL A,B,C DIMENSION A(N,NP),B(N,NP),C(N,NP),MK(N) GOTO (10,30), IOPT 10 DO 20 L=1,NP DO 20 I=1,N 20 C(I,L)=(A(I,L)-B(I,L))*MK(I) RETURN 30 DO 40 L=1,NP DO 40 I=1,N 40 C(I,L)=A(I,L)-B(I,L) RETURN END !!S MTR C**** ****************************************************************** C**** * * C**** * SUBROUTINE MTR * C**** * * C**** * PURPOSE: * C**** * PERFORM A WEIGHTED MONOTONIC REGRESSION, * C**** * IN ASCENDING ORDER. * C**** * * C**** * USAGE: * C**** * CALL MTR (Z,W,N) * C**** * * C**** * SUBPROGRAMS CALLED: * C**** * NONE. * C**** * * C**** * * C**** * PARAMETERS: * C**** * * C**** * NAME TYPE DESCRIPTION * C**** * * C**** * Z DEFAULT INPUT ARRAY ON WHICH A MONOTONIC**** * C**** * REGRESSION IS TO BE PERFORMED; AND * C**** * OUTPUT ARRAY CONTAINING THE MONOTONIC * C**** * LEAST-SQUARES TRANSFORMATION OF Z. * C**** * W DEFAULT INPUT ARRAY OF WEIGHTS OF CORRESPONDING * C**** * ELEMENTS OF Z. * C**** * N INTEGER NUMBER OF ELEMENTS IN Z ; W. * C**** * * C**** * * C**** * LOCAL VARIABLES: * C**** * * C**** * NAME TYPE DESCRIPTION * C**** * * C**** * SDS DEFAULT IN LOOP 3: NUMBER EVENTUALLY CONTAINING * C**** * THE WEIGHTED SUM OF AN * C**** * ACTIVE (ORDER-VIOLATING) * C**** * BLOCK. * C**** * SW DEFAULT IN LOOP 3: NUMBER EVENTUALLY CONTAINING * C**** * THE SUM OF THE WEIGHTS OF * C**** * AN ORDER-VIOLATING BLOCK. * C**** * M INTEGER EQUALS N, BUT HALF ITS SIZE. * C**** * I INTEGER INDEX OF MAIN LOOP. * C**** * IT INTEGER INDEX OF ORDER-VIOLATING BLOCK * C**** * TRANSFORMATION LOOP. * C**** * * C**** * * C**** * REMARK: * C**** * TO AVOID EXECUTION-TIME INEFFICIENCIES OR IMPRECISE * C**** * RESULTS SDS, SW, Z AND W SHOULD BE OF THE SAME * C**** * TYPE, DOUBLE OR EXTENDED PRECISION. * C**** * * C**** * (ERNST VAN WANING). * C**** * * C**** ****************************************************************** SUBROUTINE MTR (Z,W,N) C************************************************* C*** * C*** WEIGHTED MONOTONE REGRESSION * C*** * C************************************************* INTEGER I,M,J,IT,IT1 REAL Z,W,SDS,SW DIMENSION Z(N),W(N) M=N DO 1 I=2,M IF (Z(I) .GE. Z(I-1)) GOTO 1 SDS=Z(I)*W(I) SW=W(I) IT=I 3 IT=IT-1 SDS=SDS+Z(IT)*W(IT) SW=SW+W(IT) Z(IT)=SDS/SW IF (IT .EQ. 1) GOTO 2 IF (Z(IT-1) .GT. Z(IT)) GOTO 3 2 IT1=IT+1 DO 4 J=IT1,I 4 Z(J)=Z(IT) 1 CONTINUE RETURN END !!S ORDINA SUBROUTINE ORDINA(ID,Z,HD,HZ,NC) C*********************************************************** C*** * C*** MONOTONE REGRESSION FOR ORDINAL VARIABLES * C*** IF THE QUANTIFICATIONS OF THE FIRST AND THE LAST * C*** ARE IDENTICAL,THE SIGN OF THE NON SCALED DATA VECTOR * C*** IS CHANGED ONCE,AND THE REGRESSION IS DONE AGAIN. * C*** * C*** SUBROUTINES CALLED: MTR * C*** * C*********************************************************** DIMENSION ID(NC),HD(NC),HZ(NC),Z(NC) INTEGER ID REAL Z,HZ,HD LOGICAL LOG,TEST LOG=.FALSE. TEST=.FALSE. C IT=0 4 I=0 IT=IT+1 DO 1 K=1,NC IF (ID(K).EQ.0) GOTO 1 IF(LOG)Z(K)=(-1.0)*Z(K) I=I+1 HD(I)=FLOAT(ID(K)) HZ(I)=Z(K) 1 CONTINUE IF(I.LE.1) GOTO 6 C CALL MTR(HZ,HD,I) C LOG=.FALSE. IF (TEST) GOTO 5 IF (HZ(1).EQ.HZ(I)) LOG=.TRUE. IF (LOG) TEST=.TRUE. IF (LOG) GOTO 4 5 CONTINUE I=0 DO 2 K=1,NC IF(ID(K).EQ.0) GOTO 2 I=I+1 Z(K)=HZ(I) 2 CONTINUE 6 RETURN END !!S OUTEZA SUBROUTINE OUTEZA(IN,X,E,EE,Z,A,ID,CO,HH,HZ,HD,H3,H4,H5,H7,MK, 1 NC,NP,J,IT,IW,N,IO,NC2,IOPT) C********************************************************************* C * C OUTPUT OF MULTIPLE CATEGORY QUANTIFICATIONS EE AND CATEGORY * C CENTROIDS; FOR SINGLE VARIABLES ALSO OUTPUT OF THE SINGLE * C CATEGORY QUANTIFICATIONS Z, RANK-ONE QUANTIFICATIONS E AND * C AVERAGE RANK-ONE QUANTIFICATIONS Z*CO * C J=VARIABLE NUMBER,NC=NO OF CATEGORIES,IT=VARIABLE TYPE * C * C SUBROUTINES CALLED INDPRO,DIAPRO,PLSCA2 * C * C********************************************************************* DIMENSION E(NC,NP),Z(NC),A(NP),ID(NC),EE(NC,NP),HZ(NC),HD(NC), + HH(NC2),H3(NC2),H4(NC2),H5(NC,NP),IN(N),X(N,NP),CO(NP), + H7(NC,NP),MK(N) COMMON IX(5),IY(5),IZ(5) C - COMPUTATION OF THE CATEGORY CENTROIDS CALL INDPRO(IN,X,H5,N,NC,NP) CALL DIAPRO(ID,H5,H5,NC,NP) IF(IT.EQ.3)GOTO 20 DO 10 L=1,NP DO 10 I=1,NC 10 H7(I,L)=Z(I)*CO(L) 20 IF(IX(4).EQ.0) GOTO 80 C - OUTPUT CATEGORY QUANTIFICATIONS IF(IT.EQ.0)WRITE(IO,600)J IF(IT.EQ.1)WRITE(IO,601)J IF(IT.EQ.2)WRITE(IO,602)J IF(IT.EQ.3)WRITE(IO,603)J WRITE(IO,500) C - EE: THE MULTIPLE CATEGORY QUANTIFICATIONS FOR SINGLE AND MULTIPLE C VARS. DO 30 I=1,NC 30 WRITE(IO,501)I,ID(I),(EE(I,L),L=1,NP) C - THE CATEGORY CENTROIDS WRITE(IO,513) DO 40 I=1,NC 40 WRITE(IO,501)I,ID(I),(H5(I,L),L=1,NP) IF(IT.EQ.3) GOTO 70 WRITE(IO,510)(I,ID(I),Z(I),I=1,NC) C - THE RANK-ONE QUANTIFICATIONS WRITE(IO,505) DO 50 I=1,NC 50 WRITE(IO,501)I,ID(I),(E(I,L),L=1,NP) C - THE AVERAGE RANK-ONE QUANTIFICATIONS WRITE(IO,506) DO 60 I=1,NC 60 WRITE(IO,501)I,ID(I),(H7(I,L),L=1,NP) 70 WRITE(IO,504) 80 IF(IZ(4).EQ.0) GOTO 130 IU=IZ(4) DO 90 I=1,NC 90 WRITE(IU,502)J,I,(EE(I,L),L=1,NP) DO 100 I=1,NC 100 WRITE(IU,502)J,I,(H5(I,L),L=1,NP) IF(IT.EQ.3) GOTO 130 WRITE(IU,512)(J,I,Z(I),I=1,NC) DO 110 I=1,NC 110 WRITE(IU,502)J,I,(E(I,L),L=1,NP) DO 120 I=1,NC 120 WRITE(IU,502)J,I,(H7(I,L),L=1,NP) 130 IF(IY(4).EQ.0) GOTO 150 C - SEVERAL PLOTS FOR CATEGORY QUANTIFICATIONS CALL PLSCA2(Z,E,EE,HH,H3,H4,H5,H7,IO,N,NP,NC,IT,IW,J,NC2) C**** PREPARE ARRAY H5 FOR PLOT OF CENTROIDS FOR MULTIPLE VARIABLES C**** AND AVERAGE RANK-ONE FOR SINGLE VARIABLES FOR ALL VARIABLES IN C**** ONE PLOT 150 IF(IT.EQ.3) RETURN DO 160 L=1,NP DO 160 I=1,NC 160 H5(I,L)=H7(I,L) C - FORMATS 500 FORMAT(//' Cat Freq Multiple category quantifications'/) 501 FORMAT(1H ,2I5,10F7.2/(11X,10F7.2)) 502 FORMAT(2I5,12F8.3/,(10X,12F8.3)) 504 FORMAT(//3H ,77(1H-)) 510 FORMAT(//' Cat Freq Single category quantifications', + //(1H ,2I5,F7.2)) 512 FORMAT(2I5,F8.3) 513 FORMAT(//' Cat Freq Category centroids'/) 505 FORMAT(//' Cat Freq Rank one quantifications'/) 506 FORMAT(//' Cat Freq Average rank one quantifications'/) 600 FORMAT(//14H Variable no,I3,5X,20HType: single nominal) 601 FORMAT(//14H Variable no,I3,5X,20HType: single ordinal) 602 FORMAT(//14H Variable no,I3,5X,22HType: single numerical) 603 FORMAT(//14H Variable no,I3,5X,22HType: multiple nominal) RETURN END !!S OUTLOS SUBROUTINE OUTLOS(IN,X,E,EE,A,ID,HS,HR,DI,NC,MCNP,NP, + N,M,J,IT,IO,IOPT) C********************************************************************* C * C DISCRIMINATION MATRICES; * C TOTAL LOSS; MULTIPLE FIT; MULTIPLE LOSS AND SINGLE LOSS. * C * C SUBROUTINES CALLED:PROIN2 * C * C********************************************************************* DIMENSION IN(N),X(N,NP),E(NC,NP),EE(NC,NP),A(NP),ID(NC), + HS(NP,NP),HR(NP,NP),DI(NP,NP) COMMON IX(5),IY(5),IZ(5) C - COMPUTE DISCRIMINATION MATRIX DO 20 L2=1,NP DO 20 L1=1,NP HR(L1,L2)=0.0 DO 10 K=1,NC 10 HR(L1,L2)=HR(L1,L2)+E(K,L1)*ID(K)*E(K,L2) 20 HR(L1,L2)=HR(L1,L2)/FLOAT(N) C - COMPUTE TOTAL LOSS DO 30 L2=1,NP DO 30 L1=1,NP 30 HS(L1,L2)=DI(L1,L2)-HR(L1,L2) IF(IX(5).EQ.0)GOTO 70 IF(IT.EQ.0)WRITE(IO,600)J IF(IT.EQ.1)WRITE(IO,601)J IF(IT.EQ.2)WRITE(IO,602)J IF(IT.EQ.3)WRITE(IO,603)J 600 FORMAT(//14H Variable no,I3,5X,20HType: single nominal) 601 FORMAT(//14H Variable no,I3,5X,20HType: single ordinal) 602 FORMAT(//14H Variable no,I3,5X,22HType: single numerical) 603 FORMAT(//14H Variable no,I3,5X,22HType: multiple nominal) C - PRINT DISCRIMINATION MATRIX (=MULTIPLE FIT FOR MULTIPLE C VARIABLES AND SINGLE FIT FOR SINGLE VARIABLES) WRITE(IO,810) DO 40 L1=1,NP 40 WRITE(IO,906)L1,(HR(L1,L2),L2=1,NP) C - PRINT TOTAL DISPERSION AS COMPUTED IN SUBROUTINE OVRALS (DI) WRITE(IO,815) DO 50 L1=1,NP 50 WRITE(IO,906)L1,(DI(L1,L2),L2=1,NP) C - PRINT TOTAL LOSS WRITE(IO,811) DO 60 L1=1,NP 60 WRITE(IO,906)L1,(HS(L1,L2),L2=1,NP) 70 IF(IZ(5).EQ.0)GOTO 110 IU=IZ(5) DO 80 L1=1,NP 80 WRITE(IU,910)J,L1,(HR(L1,L2),L2=1,NP) DO 90 L1=1,NP 90 WRITE(IU,910)J,L1,(DI(L1,L2),L2=1,NP) DO 100 L1=1,NP 100 WRITE(IU,910)J,L1,(HS(L1,L2),L2=1,NP) C - COMPUTE MULTIPLE FIT FOR SINGLE VARS 110 IF(IT.EQ.3) GOTO 230 DO 130 L2=1,NP DO 130 L1=1,NP HS(L1,L2)=0.0 DO 120 K=1,NC 120 HS(L1,L2)=HS(L1,L2)+EE(K,L1)*ID(K)*EE(K,L2) 130 HS(L1,L2)=HS(L1,L2)/FLOAT(N) IF(IX(5).EQ.0)GOTO 165 WRITE(IO,830) DO 160 L1=1,NP 160 WRITE(IO,906)L1,(HS(L1,L2),L2=1,NP) 165 IF(IZ(5).EQ.0)GOTO 169 IU=IZ(5) DO 167 L1=1,NP 167 WRITE(IU,910)J,L1,(HS(L1,L2),L2=1,NP) C - COMPUTE MULTIPLE LOSS FOR SINGLE VARIABLES 169 DO 140 L2=1,NP DO 140 L1=1,NP 140 DI(L1,L2)=DI(L1,L2)-HS(L1,L2) C - COMPUTE SINGLE LOSS FOR SINGLE VARIABLES DO 150 L2=1,NP DO 150 L1=1,NP 150 HS(L1,L2)=HS(L1,L2)-HR(L1,L2) IF(IX(5).EQ.0)GOTO 190 WRITE(IO,812) DO 170 L1=1,NP 170 WRITE(IO,906)L1,(DI(L1,L2),L2=1,NP) WRITE(IO,814) DO 180 L1=1,NP 180 WRITE(IO,906)L1,(HS(L1,L2),L2=1,NP) 190 IF(IZ(5).EQ.0) GOTO 230 DO 210 L1=1,NP 210 WRITE(IU,910)J,L1,(DI(L1,L2),L2=1,NP) DO 220 L1=1,NP 220 WRITE(IU,910)J,L1,(HS(L1,L2),L2=1,NP) 230 IF(IX(5).EQ.1)WRITE(IO,500) 810 FORMAT(//24H Discrimination matrix/) 500 FORMAT(//3H ,77(1H-)) 811 FORMAT(//13H Total loss/) 812 FORMAT(//16H Multiple loss/) 814 FORMAT(//14H Single loss/) 815 FORMAT(//19H Total dispersion/) 830 FORMAT(//15H Multiple fit/) 906 FORMAT(3H ,I5,10F7.2/(1H ,7X,10F7.2)) 910 FORMAT(2I5,12F8.3/(10X,12F8.3)) RETURN END !!S OUTPAR SUBROUTINE OUTPAR(IN,X,E,EE,A,Z,NC,KC,IT,ID,H,H2,HH,HZ,HD,H3, 1 H4,H5,H7,HS,CO,MK,HP,HR,DI,M,N,NP,MCNP,MC,IO, 2 IC,IC2,IW,IA,IOPT) C******************************************************************* C * C OUTPUT OF OBJECT SCORES X, MULTIPLE CATEGORY QUANTIFICATIONS * C EE, SINGLE CATEGORY QUANTIFICATIONS Z, RANK-ONE QUANTIFICATI- * C ONS E=Z*A , WEIGHTS A, AVERAGE RANK-ONE QUANTIFICATIONS Z*B * C AND COMPONENT LOADINGS B FOR EACH VARIABLE AND DIMENSION. * C * C SUBROUTINES CALLED: OUTLOS,PROIN2,PLTCOM,COMPO,OUTEZA, * C REORDE, PLACA2, PLTOBV * C * C******************************************************************* DIMENSION X(N,NP),A(NP,M),E(MCNP),EE(MCNP),NC(M),KC(M),IT(M), 1 IW(M),Z(MC),ID(MC),H(N,NP),HS(NP,NP),CO(NP,M),IN(N,M), 2 H2(M),HH(IC2),HZ(IC),HD(IC),H3(IC2),H4(IC2),H7(IC,NP), 3 HP(NP),H5(MCNP),HR(NP,NP),DI(NP,NP,M),MK(N) COMMON IX(5),IY(5),IZ(5) DATA AA/4H - / C - DISCRIMINATION-MATRIX; TOTAL LOSS; C - MULTIPLE FIT; MULTIPLE LOSS; SINGLE LOSS. IF(IX(5).EQ.0.AND.IZ(5).EQ.0)GOTO 51 IF(IX(5).EQ.1)WRITE(IO,814) DO 222 J=1,M I1=KC(J)-NC(J)+1 I2=(KC(J)-NC(J))*NP+1 222 CALL OUTLOS(IN(1,J),X,E(I2),EE(I2),A(1,J),ID(I1),HS, 1 HR,DI(1,1,J),NC(J),MCNP,NP,N,M,J,IT(J),IO,IOPT) C 51 IF(IA.EQ.1) GOTO 65 C - COMPUTE COMPONENT LOADINGS FOR SINGLE VARIABLES DO 29 J=1,M IF(IT(J).EQ.3) GOTO 29 I2=(KC(J)-NC(J))*NP+1 CALL PROIN2(IN(1,J),E(I2),H,N,NC(J),NP,IOPT) DO 20 L1=1,NP 20 CO(L1,J)=0.0 DO 25 L=1,NP DO 25 I=1,N 25 CO(L,J)=CO(L,J)+X(I,L)*H(I,L) DO 27 L=1,NP 27 CO(L,J)=CO(L,J)/(FLOAT(N)*A(L,J)) 29 CONTINUE C - WRITE WEIGHTS AND COMPONENT LOADINGS TO UNIT IF(IZ(3).EQ.0) GOTO 223 IU=IZ(3) DO 180 J=1,M IF(IT(J).EQ.3)WRITE(IU,909)J,(AA,L=1,NP) 180 IF(IT(J).NE.3)WRITE(IU,910)J,(A(L,J),L=1,NP) DO 226 J=1,M IF(IT(J).EQ.3)WRITE(IU,909)J,(AA,L=1,NP) 226 IF(IT(J).NE.3)WRITE(IU,910)J,(CO(L,J),L=1,NP) C - PLOT 223 IF(NP.EQ.1)GOTO 225 IF (IY(3).EQ.1) CALL PLTCOM(CO,H5,H2,M,NP,IO,IT) C - PRINT WEIGHTS AND COMPONENT LOADINGS 225 IF(IX(3).EQ.0) GOTO 65 WRITE(IO,891) DO 140 J=1,M IF(IT(J).EQ.3)WRITE(IO,907)J,(AA,L=1,NP) 140 IF(IT(J).NE.3)WRITE(IO,906)J,(A(L,J),L=1,NP) WRITE(IO,810) DO 224 J=1,M IF(IT(J).EQ.3)WRITE(IO,907)J,(AA,L=1,NP) 224 IF(IT(J).NE.3)WRITE(IO,906)J,(CO(L,J),L=1,NP) C C - COMPUTE COMPONENT LOADINGS FOR MULTIPLE VARIABLES 65 IF(IX(3).EQ.0.AND.IZ(3).EQ.0)GOTO 165 DO 129 J=1,M IF(IT(J).NE.3) GOTO 129 IF(IX(3).NE.0)WRITE(IO,888)J I1=KC(J)-NC(J)+1 I2=(KC(J)-NC(J))*NP+1 CALL PROIN2(IN(1,J),E(I2),H,N,NC(J),NP,IOPT) CALL COMPO(E(I2),X,ID(I1),H,HP,HS,N,NP,NC(J),J,IO) 129 CONTINUE C - MULTIPLE AND SINGLE CATEGORIES 165 IF((IX(4).EQ.0).AND.(IY(4).EQ.0).AND.(IZ(4).EQ.0)) GOTO 40 C WRITE(IO,908) IF(IX(4).EQ.1)WRITE(IO,815) DO 220 J=1,M I1=KC(J)-NC(J)+1 I2=(KC(J)-NC(J))*NP+1 C - COMPUTATION AND OUTPUT CATEGORY QUANTIFICATIONS 220 CALL OUTEZA(IN(1,J),X,E(I2),EE(I2),Z(I1),A(1,J),ID(I1),CO(1,J), 1 HH,HZ,HD,H3,H4,H5(I2),H7,MK,NC(J),NP,J,IT(J), 2 IW(J),N,IO,NC(J)*2,IOPT) C-PLOT OF ALL MULTIPLE QUANTIFICATIONS & RANK-ONE QUANTIFICATIONS IF(NP.EQ.1)GOTO 40 DO 300 J=1,M I2=(KC(J)-NC(J))*NP+1 300 CALL REORDE(H5(I2),EE,MCNP,MC,NC(J),KC(J),NP,IO) CALL PLACA2(EE,H5,E,MC,NP,IO) C - OBJECT SCORES C - PRINT 40 IF(IX(2).EQ.0) GOTO 70 WRITE(IO,890) DO 50 I=1,N 50 WRITE(IO,906)I,(X(I,L),L=1,NP) C - UNIT 70 IF(IZ(2).EQ.0) GOTO 90 IU=IZ(2) DO 80 I=1,N 80 WRITE(IU,910)I,(X(I,L),L=1,NP) C - PLOT 90 IF(IY(2).EQ.0) RETURN CALL PLTOBV(X,IN(1,1),H,N,NP,IO) C C 810 FORMAT (///42H Component loadings for single variables/) 888 FORMAT(///37H Component loadings for variable no,I3/) 890 FORMAT (///16H Object scores/) 891 FORMAT (1H1///10H Weights/) 906 FORMAT(3H ,I5,10F7.2/(1H ,7X,10F7.2)) 907 FORMAT(3H ,I5,10(3X,A4)/(1H ,7X,10(3X,A4))) 908 FORMAT(1H1) 909 FORMAT(I5,5X,12(4X,A4)/(10X,12(4X,A4))) 910 FORMAT(I5,5X,12F8.3/,(10X,12F8.3)) 814 FORMAT (///44H1 *****************************************, +25H*************************, + /40H * P A R T I T I O N I N G L O S S,27X, +2H *,/44H *****************************************, +25H*************************, + /44H * For single variables : Discrimination , +25Hmatrix = Single fit *, + /44H * Total loss = Total dispersion - Single, +4H Fit,19X,2H */34H * Total loss = Multiple loss +, +12H Single loss,21X,2H *,/27H * Multiple loss = Total, +26H dispersion - Multiple fit,14X,2H *, +/45H * Single loss = Multiple fit - Single fit,22X,2H *, +/54H * ------------------------------------------------, +15H------------ *,/31H * For multiple variables : , +38HDiscrimination matrix = Multiple fit *, + /44H * Total loss = Total dispersion - Multip, +6Hle fit,17X,2H *,/33H ******************************, +36H************************************,/) 815 FORMAT (/// + '1 **********************************************************', + '*************'/ + ' * O P T I M A L C A T E G O R Y Q U A N T I F I C A', + ' T I O N S *'/ + ' **********************************************************', + '*************'/) RETURN END !!S OUTUNI SUBROUTINE OUTUNI(IO) C******************************************* C*** * C*** COMMENT ON OUTPUT TO OTHER UNITS * C*** THAN THE PRINTER * C*** * C******************************************* C*** COMMON IX(5),IY(5),IZ(5) 10 FORMAT(50H For each object the object number and the object / + 52H scores for each dimension have been written to unit,I3,/ + 38H with format (I5,12F8.3/,(5X,12F8.3)).,/////) 20 FORMAT(49H The weights and the component loadings have been/ + 16H written to unit,I3,26H with format (10X,12F8.3)./ + /////) 30 FORMAT(/////49H For each variable the variable number, the categ, + 'ory number'/' and the multiple category quantifications'/ + 26H have been written to unit,I3, + 40H with format (2I5,12F8.3/,(10X,12F8.3))., + //49H For each variable the variable number, the categ, + 'ory number'/' and the category centroids have been'/ + 16H written to unit,I3,24H with format (2I5,F8.3)., + //56H For each single variable the variable number, the categ, + 'ory number'/' and the single category quantifications'/ + 26H have been written to unit,I3,24H with format (2I5,F8.3)., + //56H For each single variable the variable number, the categ, + 'ory number'/' and the rank-one quantifications'/ + 26H have been written to unit,I3,24H with format (2I5,F8.3)., + //56H For each single variable the variable number, the categ, + 'ory number'/' and the average rank-one quantifications'/ + 26H have been written to unit,I3, + 40H with format (2I5,12F8.3/,(10X,12F8.3)).,/////) 40 FORMAT(/////48H For each single variable the discrimination mat, + 'rix'/' total dispersion, total loss, multiple fit, '/ + ' multiple loss and single loss'/ + ' have been written to unit,I3,'/ + ' with format (I5,12F8.3/,(5X,12F8.3)).'/ + //54H For multiple variables the discrimination matrix, tot, + 'al dispersion'/' and total loss have been written to', + ' unit,I3,'/ + 38H with format (I5,12F8.3/,(5X,12F8.3)).,/////) C*** IF(IZ(2).GT.1)WRITE(IO,10)IZ(2) IF(IZ(3).GT.1)WRITE(IO,20)IZ(3) IF(IZ(4).GT.1)WRITE(IO,30)IZ(4),IZ(4),IZ(4),IZ(4) IF(IZ(5).GT.1)WRITE(IO,40)IZ(5),IZ(5) RETURN END !!S PASSIV SUBROUTINE PASSIV(IN,X,HULP,H,HP,HS,ID,NP,NC,IC,N,IO,MM,MTOT, + ICNP) C********************************************************************* C*** * C*** OUTPUT FOR PASSIVE VARIABLES (COMPONENT-LOADINGS & CENTROIDS) * C*** * C*** SUBROUTINES CALLED: INDPRO,DIAPRO,PROIN2, COMPO * C*** * C********************************************************************* REAL X,H,HP,HS,HULP INTEGER IN,ID,NP,NC,N,M,IO,IX,IY,IZ,IP,MM,MTOT,M1 DIMENSION IN(N,MTOT),NC(MTOT),ID(IC),HULP(ICNP), + X(N,NP),H(N,NP),HP(NP),HS(NP,NP) COMMON IX(5),IY(5),IZ(5) IL=1 IM=1 IF(IX(4).EQ.0.AND.IZ(4).EQ.0)IL=0 IF(IX(3).EQ.0.AND.IZ(3).EQ.0)IM=0 WRITE(IO,503) IF(IX(3).NE.0.OR.IX(4).NE.0)WRITE (IO,504) C COMPUTATION OF MARGINAL FREQUENCIES FOR PASSIVE VARIABLES M1=(MTOT-MM)+1 DO 50 J=M1,MTOT NCP=NC(J) DO 5 K=1,NCP 5 ID(K)=0 IOPT=2 DO 35 I=1,N IF(IN(I,J).LE.NC(J)) GOTO 7 IOPT=1 GOTO 35 7 IR=IN(I,J) ID(IR)=ID(IR)+1 35 CONTINUE IF(IL.EQ.0)GOTO 45 C - COMPUTATION OF THE CATEGORY CENTROIDS FOR PASSIVE VARS. CALL INDPRO(IN(1,J),X,HULP,N,NC(J),NP) CALL DIAPRO(ID,HULP,HULP,NC(J),NP) LL=0 NCNP=NCP*NP IF(IX(4).EQ.0) GOTO 41 WRITE(IO,514)J DO 40 I=1,NCP 40 WRITE(IO,501)I,ID(I),(HULP(LL+L),L=I,NCNP,NCP) 41 IF(IZ(4).EQ.0)GOTO 45 IU=IZ(4) DO 42 I=1,NCP 42 WRITE(IU,502)J,I,(HULP(LL+L),L=I,NCNP,NCP) C C - COMPUTE COMPONENT LOADINGS FOR PASSIVE VARIABLES 45 IF(IM.EQ.0)GOTO 50 IF(IX(3).NE.0)WRITE(IO,516)J CALL PROIN2(IN(1,J),HULP,H,N,NC(J),NP,IOPT) CALL COMPO(HULP,X,ID,H,HP,HS,N,NP,NC(J),J,IO) 50 CONTINUE 501 FORMAT(1H ,2I7,8F7.2/(1H ,14X,8F7.2)) 502 FORMAT(2I5,12F8.3,(10X,12F8.3)) 503 FORMAT(1H1) 504 FORMAT (///44H1 *****************************************, +25H*************************, + /40H * OUTPUT FOR PASSIVE VARIABLES ,27X, +2H *,/44H *****************************************, +25H*************************///) 516 FORMAT(//37H Component loadings for variable no,I3/) 514 FORMAT(//'0 Category Freq Category centroids: variable no',I3/) RETURN END !!S PLACA2 SUBROUTINE PLACA2(A5,HA,HB,MC,NP,IO) C*************************************************************** C*** * C*** PLOT ALL CENTROIDS FOR MULTIPLE VARIABLE AND AVERAGE * C*** RANK-ONE QUANTIFICATIONS FOR SINGLE VARIABLES * C*** HA=H5(USED AS WORKARRAY) HB=E(USED AS WORKARRAY) * C*** A5 CONTAINS REORDED CENTROIDS AND AVERAGE RANK-ONE QUAN. * C*** SUBROUTINES CALLED: PLOT * C*** * C*************************************************************** C*** INTEGER MC,NP,IO REAL HA,HB,A5 DIMENSION A5(MC,NP),HA(MC,2),HB(MC,1) MP=NP-1 DO 110 L=1,MP DO 103 I=1,MC 103 HA(I,1)=A5(I,L) LP=L+1 DO 110 K=LP,NP WRITE(IO,5)L,K DO 105 I=1,MC 105 HA(I,2)=A5(I,K) CALL PLOT(HA(1,1),HA(1,2),HB(1,1),MC,MC,1,IO) 110 CONTINUE 5 FORMAT(1H1,8X,'Centroids for multiple, average rank one for', + ' single vars H:',I2,' V:',I2) RETURN END !!S PLSCA2 SUBROUTINE PLSCA2(Z,E,EE,H1,H2,H3,H5,H7,IO,N,NP,NC,IT,IW,J,NC2) C*********************************************************** C*** * C*** SUBROUTINES CALLED: PLOT * C*** * C*********************************************************** C*** DIMENSION E(NC,NP),Z(NC),H1(NC2),H2(NC2),H3(NC2), 1 H7(NC,NP),EE(NC,NP),H5(NC,NP) REAL H1,H2,H3,H5,E,EE,Z INTEGER IW,IO,IT,J,L,K C*** IF(IW.EQ.0) RETURN IF(IT.EQ.3) GOTO 500 C PLOT OF THE CATEGORY NUMBERS AND THE SINGLE C CATEGORY QUANTIFICATIONS. DO 20 K=1,NC 20 H1(K)=K WRITE(IO,1)J C IF(IT.EQ.0)WRITE(IO,2) C IF(IT.EQ.1)WRITE(IO,3) C IF(IT.EQ.2)WRITE(IO,4) C IF(IT.EQ.3)WRITE(IO,5) CALL PLOT(H1,Z,H3,NC,NC,1,IO) IF(NP.EQ.1)RETURN C PLOT OF THE CENTROIDS AND AVERAGE RANK-ONE QUANTIFICATIONS MM=NP-1 DO 210 L=1,MM DO 203 I=1,NC H1(I)=H5(I,L) 203 H1(NC+I)=H7(I,L) LL=L+1 DO 210 K=LL,NP WRITE(IO,16)L,K,J C WRITE(IO,7)L,K,J C WRITE(IO,19) DO 205 I=1,NC H2(I)=H5(I,K) 205 H2(NC+I)=H7(I,K) CALL PLOT(H1,H2,H3,NC2,NC,1,IO) 210 CONTINUE C PLOT OF THE MULTIPLE AND THE RANK-ONE QUANTIFICATIONS MP=NP-1 DO 110 L=1,MP DO 103 I=1,NC H1(I)=EE(I,L) 103 H1(NC+I)=E(I,L) LP=L+1 DO 110 K=LP,NP C WRITE(IO,6) WRITE(IO,7)L,K,J C WRITE(IO,9) DO 105 I=1,NC H2(I)=EE(I,K) 105 H2(NC+I)=E(I,K) CALL PLOT(H1,H2,H3,NC2,NC,1,IO) 110 CONTINUE C PLOT OF THE MULTIPLE QUANTIFICATIONS AND THE CENTROIDS 500 IF(NP.EQ.1)RETURN MP=NP-1 DO 125 L=1,MP DO 115 I=1,NC H1(I)=EE(I,L) 115 H1(NC+I)=H5(I,L) LP=L+1 DO 125 K=LP,NP WRITE(IO,10)L,K,J C WRITE(IO,7)L,K,J C WRITE(IO,11) DO 120 I=1,NC H2(I)=EE(I,K) 120 H2(NC+I)=H5(I,K) CALL PLOT(H1,H2,H3,NC2,NC,1,IO) 125 CONTINUE 1 FORMAT(1H1,14X,'Category number (N) vs. quantification (C)', + ' Var',I4) 7 FORMAT(1H1,12X,'Multiple (N) & rank one quantifications (C)', + ' H:',I2,' V:',I2,' Var',I4) C 7 FORMAT(1H+,41X,23H (HORIZONTAL VARIATE NO,I2,24H AND VERTICAL VARI C *ATE NO,I2,16H) OF VARIABLE NO,I5) 10 FORMAT(1H1,11X,'Multiple quantifications (N) & centroids (C)', + ' H:',I2,' V:',I2,' Var',I4) 16 FORMAT(1H1,8X,'Category centroids (N) & average rank one ', + 'qnt. (C) H:',I2,' V:',I2,' Var',I4) RETURN END !!S PLOT SUBROUTINE PLOT(X,Y,IND,NITEM,NRC,IDENT,IWRITE) C ****************************************************************** C * * C * P L O T * C * * C * PURPOSE: PRINTPLOT OF Y VERSUS X * C * IF IDENT.GT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * INTEGER NUMBERS 1-9,0,1-9... , THE REST OF THE POINTS* C * BY CHARACTERS A-I,J,A-I... * C * IF IDENT.EQ.0, ALL POINTS ARE PRINTED AS STARS * C * IF IDENT.LT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * STARS, THE OTHERS BY A 'D' * C * * C * SUBROUTINES CALLED: NONE * C * * C * AUTHOR: WILLEM HEISER RELEASED: MARCH 1978 * C * ADAPTED : NOVEMBER 1982 * C * * C ****************************************************************** C C TO ADJUST THE PLOT CHANGE THE CONSTANTS UNDER WHICH C ## IS PLACED ACCORDING TO THE COMMENTS C DIMENSION LINE(64),XPR(10),LABEL(21),X(NITEM),Y(NITEM),IND(NITEM) C ##=NN+1 (NN=INT(NC/7)) C ##=NUMBER OF COLUMNS PER LINE=NC INTEGER BLANK,STAR DATA BLANK,STAR,MORE/1H ,1H*,1HM/,NL/20/,NC/64/,NN/9/ C ##=NC C ##=NUMBER OF LINES PER PAGE DATA LABEL/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, + 1HJ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1H0/ C C For some mysterious reason, LABEL(1) does not always initialize C properly, causing sloppy plots. An explicit assign is made. SvB. C LABEL(1)=LABEL(21) C 1 FORMAT(9X, 3H.--,9 (7H+------),5H+---.) C ##=NN 2 FORMAT(1X,F7.2,1X,1H|,2X,64A1,3X,1H|) C ##=NC 3 FORMAT(1X,F7.2,1X,1H|,69X,1H|) C ##=NC+5 4 FORMAT(8X, 10F7.2) C ##=NN+1 C DO 10 I=1,NITEM 10 IND(I)=I C C THE INDICES ARE ORDERED SUCH THAT THEY WOULD ORDER Y IN DESCENDING C ORDER C M=NITEM 20 M=M/2 IF (M.LT.1) GO TO 40 K=NITEM-M J=1 41 I=J 49 L=I+M IF (Y(IND(I)).GE.Y(IND(L))) GO TO 60 II=IND(I) IND(I)=IND(L) IND(L)=II I=I-M IF (I.GE.1) GO TO 49 60 J=J+1 IF (J-K) 41,41,20 40 CONTINUE C C THE PLOT IS ADJUSTED TO THE RANGE OF THE 'LONGEST' AXIS C XMAX=X(1) XMIN=X(1) DO 70 I=2,NITEM IF (X(I).GT.XMAX) XMAX=X(I) IF (X(I).LT.XMIN) XMIN=X(I) 70 CONTINUE SPANX=XMAX-XMIN SPANY=Y(IND(1))-Y(IND(NITEM)) IF (SPANY.LE.SPANX) GO TO 75 XMIN=XMIN-(SPANY-SPANX)/2.0 SPANX=SPANY 75 STEPX=SPANX/(7.0*NN-1.0) STEPY=(STEPX*NC)/(NL-0.5) HSTEP=STEPY/2.0 TOP=(Y(IND(NITEM))+Y(IND(1))+SPANX)/2.0 C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE(IWRITE,1) L=1 I=1 YPR=TOP+STEPY 80 YPR=YPR-STEPY IF (I.GT.NITEM) GO TO 110 IF ((YPR-Y(IND(I))).GT.HSTEP) GO TO 110 90 DO 95 J=1,NC 95 LINE(J)=BLANK 100 JP=(X(IND(I))-XMIN)/STEPX+1.0 IF (LINE(JP).EQ.BLANK) GO TO 101 LINE(JP)=MORE GO TO 103 101 IF (IDENT.EQ.0) GO TO 102 IF (IDENT.GT.0) GO TO 104 IF (IND(I).LE.NRC) LINE(JP)=STAR IF (IND(I).GT.NRC) LINE(JP)=LABEL(15) GO TO 103 104 IF (IND(I).LE.NRC) LINE(JP)=LABEL(MOD(IND(I),10)+1) IF (IND(I).GT.NRC) X LINE(JP)=LABEL(MOD((IND(I)-NRC),10)+11) GO TO 103 102 LINE(JP)=STAR 103 I=I+1 IF (I.GT.NITEM) GO TO 105 IF ((YPR-Y(IND(I))).LT.HSTEP) GO TO 100 105 WRITE(IWRITE,2) YPR,LINE GO TO 115 110 WRITE(IWRITE,3) YPR 115 L=L+1 IF (L.LE.NL) GO TO 80 WRITE(IWRITE,1) C C VALUES OF X-AXIS ARE PRINTED C XPR(1)=XMIN DO 130 J=1,10 130 XPR(J+1)=XPR(J)+STEPX*7.0 WRITE(IWRITE,4) XPR RETURN END !!S PLTCOM SUBROUTINE PLTCOM(B,H5,H2,M,NP,IO,IT) C************************************************* C*** * C*** PLOT COMPONENT LOADINGS * C*** H2 & H5 WORK ARRAY * C*** MULTIPLE VARS HAVE VALUE 0.0 * C*** * C*** SUBROUTINES CALLED: PLOT * C*** * C************************************************* C*** INTEGER M,NP,IO REAL H2,H5,B DIMENSION H5(M,2),H2(M),B(NP,M),IT(M) MP=NP-1 DO 110 L=1,MP DO 103 J=1,M IF(IT(J).EQ.3)H5(J,1)=0.0 103 IF(IT(J).NE.3)H5(J,1)=B(L,J) LP=L+1 DO 110 K=LP,NP C WRITE(IO,3) WRITE(IO,5)L,K DO 105 J=1,M IF(IT(J).EQ.3)H5(J,2)=0.0 105 IF(IT(J).NE.3)H5(J,2)=B(K,J) CALL PLOT(H5(1,1),H5(1,2),H2(1),M,M,1,IO) 110 CONTINUE 3 FORMAT(23H1 Component loadings) 5 FORMAT(1H1,14X,'Component loadings -- H: Variate no',I2, + ' V: Variate no',I2) RETURN END !!S PLTOBV SUBROUTINE PLTOBV(X,HA,HB,N,NP,IO) C************************************************* C*** * C*** PLOT OBJECT-SCORES * C*** HB WORK ARRAY H(N,NP) * C*** X OBJECT-SCORES * C*** HA WORK ARRAY IN(N,M) * C*** * C*** SUBROUTINES CALLED: PLOT * C*** * C************************************************* C*** INTEGER N,NP,IO COMMON IX(5),IY(5),IZ(5) REAL HA,X DIMENSION X(N,NP),HA(N,2),HB(N,1) IF(IY(2).EQ.0)RETURN IF(NP.EQ.1)RETURN MP=NP-1 DO 110 L=1,MP DO 103 I=1,N 103 HA(I,1)=X(I,L) LP=L+1 DO 110 K=LP,NP WRITE(IO,5)L,K DO 105 I=1,N 105 HA(I,2)=X(I,K) CALL PLOT(HA(1,1),HA(1,2),HB(1,1),N,N,0,IO) 110 CONTINUE 5 FORMAT(1H1,15X,'Object scores -- H: Variate no',I2, + ' V: Variate no',I2) RETURN END !!S PLUSM2 SUBROUTINE PLUSM2(A,B,C,MK,N,NP,IOPT) C************************************************* C*** * C*** C = A + B * C*** * C************************************************* REAL A,B,C DIMENSION A(N,NP),B(N,NP),C(N,NP),MK(N) GOTO (10,30), IOPT 10 DO 20 L=1,NP DO 20 I=1,N 20 C(I,L)=(A(I,L)+B(I,L))*MK(I) RETURN 30 DO 40 L=1,NP DO 40 I=1,N 40 C(I,L)=A(I,L)+B(I,L) RETURN END !!S PRINAX SUBROUTINE PRINAX (X,E,EE,A,IT,HP,HR,HS,H,KC,NC,N,NP,MCNP, 1 M,NS,IO) C************************************************************ C * C ROTATION OF THE OBJECT SCORES X AND THE CATEGORY * C QUANTIFICATIONS E TO THE PRINCIPAL AXES SOLUTION AND * C ROTATION OF THE WEIGHTS A FOR SINGLE VARS * C * C SUBROUTINES CALLED: PRODAB * C * C************************************************************ DIMENSION X(N,NP),E(MCNP),A(NP,M),IT(M),HS(NP,NP), + HR(NP),HP(NP),H(N,NP),KC(M),NC(M),EE(MCNP) C - THE MATRIX HP CONTAINS EIGENVALUES C - WRITE EIGENVALUES OF THE DIMENSIONS WRITE(IO,701) DO 105 L=1,NP HP(L)=1/HP(L) HP(L)=HP(L)/(SQRT(FLOAT(N))*FLOAT(NS)) 105 WRITE(IO,702)L,HP(L) 701 FORMAT (///' Eigenvalues'//) 702 FORMAT(I5,3X,F8.3) C - THE MATRIX HS CONTAINS THE EIGENVECTORS (K) OF C (1/SQRT(MS)* X )'(1/SQRT(MS) * X) C - ROTATE OBJECT SCORES X WITH THE EIGENVECTORS K (HS) CALL PRODAB(X,HS,H,N,NP) C - ROTATE THE MULTIPLE CATEGORY QUANTIFICATIONS DO 200 J=1,M IF(N.LT.NC(J)) GOTO 300 CALL PRODAB(EE((KC(J)-NC(J))*NP+1),HS,H,NC(J),NP) 200 CALL PRODAB(E((KC(J)-NC(J))*NP+1),HS,H,NC(J),NP) C - ROTATE THE WEIGHTS A DO 250 J=1,M IF(IT(J).EQ.3) GOTO 250 CALL PRODAB(A(1,J),HS,HR,1,NP) 250 CONTINUE RETURN 300 WRITE(IO,501) STOP 501 FORMAT( 1H1, + ' The number of categories of one variable is greater than', + ' than the number of objects.'// + ' Subroutine PRINAX will give wrong results.'//) END !!S PRINDA SUBROUTINE PRINDA(IN,N,M,I1,IO) C******************************************* C*** * C*** PRINT OUTPUT DATA MATRIX * C*** * C******************************************* C DIMENSION IN(N,M) DATA AA/4H----/ LN=N IF(I1.EQ.0)LN=MIN0(N,10) WRITE(IO,502)LN J1=1 J2=MIN0(M,14) DO 100 K=1,M,14 WRITE(IO,500)(J,J=J1,J2) WRITE(IO,503)(AA,J=J1,J2) DO 75 I=1,LN 75 WRITE(IO,501)I,(IN(I,J),J=J1,J2) J1=J1+14 J2=J2+14 J2=MIN0(M,J2) 100 CONTINUE 500 FORMAT(/3H ,5X,14I5,/) 501 FORMAT(3H ,I5,14I5) 502 FORMAT(1H1//18H Data matrix -,I5,18H rows are printed) 503 FORMAT(9X,14(1X,A4)) RETURN END !!S PROCRU SUBROUTINE PROCRU(MS,X,XX,H,HP,HR,HS,N,NP,IO,IOPT) C ********************************************************** C * * C * OBJECT SCORES CORRECTED FOR MISSING ENTRIES * C * PROCRUSTERS ORTHOGONALISATION OF OBJECT SCORES * C * SUBROUTINES CALLED : SUMPRO, TRED2, IMTQL2, PRODAB * C * * C ********************************************************** DIMENSION MS(N),X(N,NP),XX(N,NP),H(N,NP),HP(NP),HR(NP,NP), + HS(NP,NP) REAL MS,H,X,XX,HP,HR,HS DO 5 L2=1,NP DO 5 L1=1,NP 5 HS(L1,L2)=0.0 GOTO (10,50), IOPT 10 DO 20 I=1,N H(I,1)=0.0 20 IF(MS(I).NE.0)H(I,1)=1/(SQRT(MS(I))) DO 30 L=1,NP DO 30 I=1,N 30 XX(I,L)=H(I,1)*X(I,L) CALL SUMPRO(XX,HS,N,NP) CALL TRED2(NP,NP,HP,HR,HS) CALL IMTQL2(NP,NP,HP,HR,HS,IR) IF(IR.GT.0)WRITE(IO,11)IR C - COMPUTE HR = EIG * LABDA(** -1) * EIG(TRANSPOSED) DO 33 K=1,NP 33 HP(K)=1.0/SQRT(HP(K)) DO 35 L2=1,NP DO 35 L1=1,NP HR(L1,L2)=0.0 DO 35 L3=1,NP 35 HR(L1,L2)=HR(L1,L2)+HS(L1,L3)*HP(L3)*HS(L2,L3) DO 40 L=1,NP DO 40 I=1,N 40 X(I,L)=H(I,1)*XX(I,L) CALL PRODAB(X,HR,H,N,NP) RETURN 50 CALL SUMPRO(X,HS,N,NP) CALL TRED2(NP,NP,HP,HR,HS) CALL IMTQL2(NP,NP,HP,HR,HS,IR) IF(IR.GT.0)WRITE(IO,11)IR C - COMPUTE HR = EIG * LABDA(** -1) * EIG(TRANSPOSED) DO 53 K=1,NP 53 HP(K)=1.0/SQRT(HP(K)) DO 55 L2=1,NP DO 55 L1=1,NP HR(L1,L2)=0.0 DO 55 L3=1,NP 55 HR(L1,L2)=HR(L1,L2)+HS(L1,L3)*HP(L3)*HS(L2,L3) CALL PRODAB(X,HR,H,N,NP) 11 FORMAT(39H1Eigenvalue computation has failed,IR =,I3) RETURN END !!S PROIN2 SUBROUTINE PROIN2(IN,A,B,N,NC,NP,IOPT) C************************************************* C*** * C*** B = IN * A * C*** IN INDICATOR MATRIX (N,NC) GIVEN IN * C*** REDUCED FORM, I.E. AS VECTOR (N) OF * C*** CATEGORY SCORES * C*** * C************************************************* REAL A,B INTEGER IN DIMENSION IN(N),A(NC,NP),B(N,NP) GOTO (10,30), IOPT 10 DO 20 L=1,NP DO 20 I=1,N B(I,L)=0.0 20 IF(IN(I).LE.NC)B(I,L)=A(IN(I),L) RETURN 30 DO 40 L=1,NP DO 40 I=1,N 40 B(I,L)=A(IN(I),L) RETURN END !!S PRODAB SUBROUTINE PRODAB(A,B,H,N1,N2) C******************************************* C * C A = A * B * C H AUXILLIAIRY ARRAY * C * C******************************************* DIMENSION A(N1,N2),B(N2,N2),H(N1,N2) DO 10 L=1,N2 DO 10 I=1,N1 H(I,L)=0.0 DO 10 K=1,N2 10 H(I,L)=H(I,L)+A(I,K)*B(K,L) DO 20 L=1,N2 DO 20 I=1,N1 20 A(I,L)=H(I,L) RETURN END !!S RANDMA SUBROUTINE RANDMA(A,N) C********************************************************* C*** * C*** MACHINE INDEPENDENT RANDOM GENERATOR * C*** N RANDOM NUMBERS ARE PLACED IN ARRAY A * C*** WARNING:THE NECESSARY AVAILABLE INTEGER ARITHMETIC * C*** IS 125*2796203=348525375 (<2**29). * C*** * C********************************************************* DIMENSION A(N) REAL A,Z IY=315747 M=2796203 Z=2.0D0/M DO 10 I=1,N IY=IY*125 IY=IY-(IY/M)*M 10 A(I)=-1.0D0+IY*Z RETURN END !!S REORDE SUBROUTINE REORDE(H5,EE,MCNP,MC,NC,KCJ,NP,IO) C*************************************************************** C*** * C*** REORDER MULTIPLE & SINGLE CATECORIES * C*** EE CONTAINS REORDED CENTROIDS AND AVERAGE RANK-ONE QUANT.* C*** * C*************************************************************** C*** INTEGER MC,NP,IO,NC,KCJ REAL H5,EE DIMENSION H5(NC,NP),EE(MCNP) DO 110 L=1,NP DO 110 I=1,NC 110 EE((MC*(L-1))+(KCJ-NC)+I)=H5(I,L) RETURN END !!S SINGIN SUBROUTINE SINGIN(E,A,Z,ID,NC,NP,N,IO) C*********************************************************** C*** * C*** UPDATE OF E AND A FOR A SINGLE VARIABLE * C*** STANDARDIZATION OF Z (AS A COLUMN OF MATRIX IN) * C*** * C*** SUBROUTINES CALLED: WEDEVI,WESTAN * C*** * C*********************************************************** DIMENSION E(NC,NP),A(NP),Z(NC),ID(NC) C - STANDARDIZATION OF Z (AS A COLUMN OF IN) ISOMID=0 DO 11 I= 1,NC 11 ISOMID=ISOMID+ID(I) C - CHECK IF STANDARDDEVIATION NOT EQUALS ZERO DO 31 K=1,NC 31 IF(ISOMID.EQ.ID(K)) GOTO 43 CALL WEDEVI(Z,ID,NC,ISOMID) CALL WESTAN(Z,ID,NC,N) C - A= E(TRANSPOSED) * ID * Z/FLOAT(N) DO 25 L=1,NP A(L)=0.0 DO 20 I=1,NC 20 A(L)=A(L)+E(I,L)*Z(I)*ID(I) 25 A(L)=A(L)/FLOAT(N) C - E= Z * A(TRANSPOSED) DO 30 L=1,NP DO 30 I=1,NC 30 E(I,L)=Z(I)*A(L) RETURN 43 WRITE(IO,45) 45 FORMAT(//43H The analysis stops because each variable/ + 42H must have at least two filled categories./////) STOP END !!S SINGLE SUBROUTINE SINGLE(E,A,Z,ID,HZ,HD,NC,NP,N,IT) C*********************************************************** C*** * C*** UPDATE OF E,A AND Z FOR A SINGLE VARIABLE * C*** * C*** SUBROUTINES CALLED: WESTAN,ORDINA, LIDICA * C*** * C*********************************************************** DIMENSION E(NC,NP),A(NP),Z(NC),ID(NC),HZ(NC),HD(NC) C - Z= E * A 5 DO 10 I= 1,NC Z(I)=0.0 DO 10 L=1,NP 10 Z(I)=Z(I)+E(I,L)*A(L) C - REGRESSION FOR ORDINAL AND NUMERICAL VARIABLES IF(IT.EQ.1) CALL ORDINA(ID,Z,HD,HZ,NC) IF(IT.EQ.2) CALL LIDICA(ID,Z,HD,HZ,NC) C - WEIGHTED STANDARDIZATION OF Z 15 CALL WESTAN(Z,ID,NC,N) C - A= E(TRANSPOSED) * ID * Z /FLOAT(N) DO 25 L=1,NP A(L)=0.0 DO 20 I=1,NC 20 A(L)=A(L)+E(I,L)*Z(I)*ID(I) 25 A(L)=A(L)/FLOAT(N) C - E= Z * A(TRANSPOSED) DO 30 L=1,NP DO 30 I=1,NC 30 E(I,L)=Z(I)*A(L) RETURN END !!S STRESS SUBROUTINE STRESS(X,V,N,NP,NS,STR) C************************************************* C*** * C*** COMPUTATION OF STRESS PER VARIABLE * C*** * C************************************************* DOUBLE PRECISION T1,T2,STR REAL X,V INTEGER N,NP,NS DIMENSION X(N,NP),V(N,NP) T1=0.0 T2=0.0 DO 10 L=1,NP DO 10 I=1,N T1=T1+X(I,L)*V(I,L) 10 T2=T2+V(I,L)*V(I,L) STR=(-2*T1+T2+FLOAT(N*NP))/FLOAT(NS*NP*N) RETURN END !!S SUMPRO SUBROUTINE SUMPRO(A,C,N1,N2) C******************************************* C * C C = C + A(TRANSPOSED) * A * C * C******************************************* DIMENSION A(N1,N2),C(N2,N2) DO 10 L=1,N2 DO 10 K=1,N2 DO 10 I=1,N1 10 C(K,L)=C(K,L)+A(I,K)*A(I,L) RETURN END !!S TRED2 SUBROUTINE TRED2(NM,N,D,E,Z) C ****************************************************************** C * * C * T R E D 2 * C * * C * PURPOSE: REDUCES A REAL SYMMETRIC MATRIX TO A SYMMETRIC TRI- * C * DIAGONAL MATRIX USING ORTHOGONAL SIMILARITY TRANS- * C * FORMATIONS; THE ACCUMULATED ORTHOGONAL TRANSFORMA- * C * TIONS ARE RETAINED IN Z. * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK-ORIGINAL VERSION NOT IN DOUBLE PRECISION * C * * C ****************************************************************** REAL D,E,Z,H,SCALE,F,G,HH REAL ABS,SQRT,SIGN DIMENSION Z(NM,N),D(N),E(N) C IF (N .EQ. 1) GO TO 320 C ********** FOR I=N STEP -1 UNTIL 2 DO -- ********** DO 300 II =2, N I = N + 2 - II L = I - 1 H = 0.0 SCALE = 0.0 IF (L .LT.2) GO TO 130 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** DO 120 K =1, L 120 SCALE = SCALE + ABS(Z(I,K)) C IF (SCALE .NE. 0.0) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L Z(I,K) = Z(I,K) / SCALE H = H + Z(I,K) * Z(I,K) 150 CONTINUE C F = Z(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0 C DO 240 J = 1, L Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0 C ********** FORM ELEMENT OF A*U ********** DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C ********** FORM ELEMENT OF P ********** 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C ********** FORM REDUCED A ********** DO 260 J =1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) -G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0 E(1) = 0.0 C ********** ACCUMULATION OF TRANSFORMATION MATRICES ********** DO 500 I = 1, N L = I - 1 IF (D(I) .EQ. 0.0) GO TO 380 C C DO 360 J = 1, L G = 0.0 DO 340 K =1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0 IF (L .LT. 1) GO TO 500 C DO 400 J =1, L Z(I,J) = 0.0 Z(J,I) = 0.0 400 CONTINUE C 500 CONTINUE C RETURN END !!S WEDEVI SUBROUTINE WEDEVI(Z,ID,NC,ISOMID) C************************************************* C*** * C*** THE VECTOR Z IS REPLACED BY THE DEVIATIONS * C*** FROM THE WEIGHTED MEAN; ISOMID IS THE TO- * C*** TAL NUMBER OF NON-MISSING INDIVIDUALS * C*** VECTOR ID CONTAINS THE WEIGHTS * C*** * C************************************************* REAL Z,S INTEGER ID,IS,NC,ISOMID DIMENSION Z(NC),ID(NC) IS=0 DO 10 L=1,NC 10 IS=IS+L*ID(L) S=FLOAT(IS)/FLOAT(ISOMID) DO 20 L=1,NC Z(L)=0.0 20 IF(ID(L).NE.0) Z(L)=FLOAT(L)-S RETURN END !!S WESTAN SUBROUTINE WESTAN(Z,ID,NC,N) C************************************************* C*** * C*** WEIGHTED STANDARDIZATON OF VECTOR Z ON N * C*** * C************************************************* REAL Z,S INTEGER ID DIMENSION Z(NC),ID(NC) S=0.0 DO 10 I=1,NC 10 S=S+Z(I)*Z(I)*ID(I) S=SQRT(FLOAT(N)/S) DO 20 I=1,NC 20 Z(I)=Z(I)*S RETURN END