!------------------------------------------------------------------------------------------------------------ SUBROUTINE HQ (NMax, MMax, N,M,NS,PCT,X,IDX, & XMEAN,DX, RR, YORG,YHQ, NSTEPS,LS,KSTEP,BS,ANSS) ! For outputs ! Ten, XMEAN,DX, RR, YORG,YHQ, NSTEPS,LS,KSTEP,BS,ANSS) ! For outputs ! Chuong trinh xay dung phuong trinh Hoi quy tuyen tinh giua bien Phu thuoc Y va cac bien Doc lap X1, X2,... ! bang phuong phap Hoi quy tung buoc (Stepwise Regression) ! Input: ! NMax, MMax: Tuong ung la Dung luong mau va So bien cuc dai (bao gom ca bien phu thuoc) ! N, M: Tuong ung la Dung luong mau va So bien thuc te (bao gom ca mot bien phu thuoc) ! NS: =1 (Neu NS<=0 se khong thuc hien chuong trinh) ! PCT: Nguong xac dinh do sut giam sai so khi dua them 1 bien moi ! X(N*M): Manng 1 chieu chua so lieu dau vao (luu theo cot) ! IDX(M): Mang chua cac gia tri xac dinh thuoc tinh cac bien: 3-bien phu thuoc, hoac 2,1,0-cac bien doc lap ! Output: ! XMEAN(M), DX(M): Tuong ung la vector Trung binh va do lech chuan cac bien ! RR(M,M): Ma tran tuong quan giua cac bien ! YORG(N),YHQ(N): Vector chua cac gia tri bien phu thuoc (YORG) va uoc luong cua no (YHQ) ! NSTEPS,LS,KSTEP,BS,ANSS: Cac mang chua ket qua tinh cac he so HQ va cac thong tin khac ! ----------------- Cac bien Output quan trong: -------------------- ! K=NSTEPS(NN,4) : So bien doc lap duoc chon. Mang NSTEPS duoc gan ban dau =-999 ! Neu K>1 : NSTEPS(i,4), i=1,K la thu tu cac buoc ung voi i bien duoc chon ! NN : So phan tu /=-999 cua mang 1 chieu NSTEPS(J,4). NN duoc xac dinh boi cac dong lenh sau ! DO J=1,MMAX ! IF (NSTEPS(J,4).EQ.-999) THEN ! NN=J-1; EXIT ! ENDIF ! ENDDO ! NN = 0 : Khong co bien nao duoc chon !!!! ! LS(I,J): Thu tu vi tri bien J duoc chon o buoc I (I=1,NN, K=NSTEPS(I,4), j=1,K) ! BS(I,J): He so HQ cua bien J o buoc I (I=1,NN, K=NSTEPS(I,4), j=1,K) ! ANSS(I,7): FISHER-F o buoc I ! ANSS(I,8): Standard Error o buoc I ! ANSS(I,6): MulCoef R o buoc I ! ANSS(I,9): In-cep A0 o buoc I ! ! WRITE(3,"(/' No ',' FISHER-F',' Error S',' MulCoef R',' In-cep A0',30A10)") (TEN(LS(K,J)),J=1,K) ! DO I=NN,NN ! 1,NN ! K=NSTEPS(I,4) ! HSO=' ----' ! DO J=1,K ! WRITE(HSO(LS(I,J)),'(F10.4)') BS(I,J) ! Convert He so BS ve dang Character ! ENDDO ! WRITE(3,"(1X,I3,F10.2,3F10.4,30A10)") K,ANSS(I,7),ANSS(I,8),ANSS(I,6),ANSS(I,9), (HSO(LS(I,J)),J=1,K) ! ENDDO !-------------------------------------------------------------------------------------------------------------- Integer :: NMax, MMax REAL :: X(N*M),XBAR(M),STD(M),D(M), B(M), RX(M*M),R(M*M),ANS(11),T(M) INTEGER :: IDX(M), L(M), NSTEP(5) Real :: PCT ! Common OUT: INTEGER :: NSTEPS(MMAX,5), LS(MMAX,MMAX), KSTEP REAL :: BS(MMAX,MMAX), ANSS(MMAX,11) ! Common TQHQ: REAL :: XMEAN(MMAX),DX(MMAX), RR(MMAX,MMAX), YORG(NMAX), YHQ(NMAX) ! Common VARN: ! Character (Len=8) :: Ten(MMAX) ! Variable's name KSTEP=0 NSTEPS=-999 print*,"Calling HQ" IO=1 CALL CORRE(NMax, MMax, N,M,IO,X,XBAR,STD,RX,R,B,D,T) DO 105 I=1,M XMEAN(I)=XBAR(I) DX(I)=STD(I) 105 CONTINUE JK = 0 DO J = 1,M DO K = 1,J JK = JK + 1 RR(J,K) = R(JK) RR(K,J) = RR(J,K) ENDDO ENDDO IF(NS)135,135,140 135 CONTINUE GoTO 200 140 CONTINUE NSEL=1 GOTO 150 145 CALL MSTR(R,RX,M,1,0) 150 CONTINUE N35=0 DO 155 K=1,M IF (IDX(K))152,153,153 152 CONTINUE GOTO 185 153 IF(IDX(K)-3)155,154,152 154 N35=N35+1 155 CONTINUE IF (N35-1)156,157,156 156 CONTINUE GOTO 185 157 CALL STEP (NMax, MMax, M,N,RX,XBAR,IDX,PCT,NSTEP,ANS,L,B,STD,TT,D,IER, & NSTEPS,LS,KSTEP,BS,ANSS) print*,"Calling STEP done" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 160 CONTINUE !! 160 WRITE(3,*)' Bang phan du ' !! WRITE(3,*)' Mau Gia tri Y Uoc luong cua Y Phan du' MM=NSTEP(1) I1 = 0 DO 180 I=1,N DO 181 J=1,M JJ=N*(j-1)+I 181 D(J)=X(JJ) YEST=ANS(9) K=NSTEP(4) DO 170 J=1,K KK=L(J) 170 YEST=YEST+B(J)*D(KK) RESI=D(MM)-YEST I1 = I1 + 1 YORG(I1) = D(MM) YHQ (I1) = YEST !! 180 WRITE(3,15)I,D(MM),YEST,RESI 180 CONTINUE 15 FORMAT(1X,I8,F11.3,F17.3,F9.3) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 185 IF(NSEL-NS)190,200,200 190 NSEL=NSEL+1 GOTO 145 200 CONTINUE return END SUBROUTINE HQ !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CORRE(NMax, MMax, N,M,I0,X,XBAR,STD,RX,R,B,D,T) ! INCLUDE 'HQTB.INC' ! DUNG LUONG MAU VA SO BIEN CUC DAI Integer :: NMax, MMax REAL X(N*M),XBAR(M),STD(M),RX(M*M), R(M*M),B(M),D(M),T(M) !print*,"Calling CORRE" DO 100 J=1,M B(J)=0. 100 T(J)=0. K=(M*M+M)/2 DO 102 I=1,K 102 R(I)=0. FN=N L=0 IF(I0)105,127,105 105 DO 108 J=1,M DO 107 I=1,N L=L+1 107 T(J)=T(J)+X(L) XBAR(J)=T(J) 108 T(J)=T(J)/FN DO 115 I=1,N JK=0 L=I-N DO 110 J=1,M L=L+N D(J)=X(L)-T(J) 110 B(J)=B(J)+D(J) DO 115 J=1,M DO 115 K=1,J JK=JK+1 115 R(JK)=R(JK)+D(J)*D(K) GOTO 205 127 IF(N-M)130,130,135 130 KK=N GOTO 137 135 KK=M 137 DO 140 I=1,KK CALL DATA(M,D) DO 140 J=1,M T(J)=T(J)+D(J) L=L+1 140 RX(L)=D(J) FRR=KK DO 150 J=1,M XBAR(J)=T(J) 150 T(J)=T(J)/FRR L=0 DO 180 I=1,KK JK=0 DO 170 J=1,M L=L+1 170 D(J)=RX(L)-T(J) DO 180 J=1,M B(J)=B(J)+D(J) DO 180 K=1,J JK=JK+1 180 R(JK)=R(JK)+D(J)*D(K) IF(N-KK) 205,205,185 185 KK=N-KK DO 200 I=1,KK JK=0 CALL DATA(M,D) DO 190 J=1,M XBAR(J)=XBAR(J)+D(J) D(J)=D(J)-T(J) 190 B(J)=B(J)+D(J) DO 200 J=1,M DO 200 K=1,M JK=JK+1 200 R(JK)=R(JK)+D(J)*D(K) 205 JK=0 DO 210 J=1,M XBAR(J)=XBAR(J)/FN DO 210 K=1,J JK=JK+1 210 R(JK)=R(JK)-B(J)*B(K)/FN JK=0 DO 220 J=1,M JK=JK+J 220 STD(J)=SQRT(ABS(R(JK))) DO 230 J=1,M DO 230 K=J,M JK=J+(K*K-K)/2 L=M*(J-1)+K RX(L)=R(JK) L=M*(K-1)+J RX(L)=R(JK) IF(STD(J)-STD(K)) 225,222,225 222 R(JK)=1. GOTO 230 225 STD2=STD(J)*STD(K) R(JK)=R(JK)/STD2 230 CONTINUE FN=SQRT(FN-1.) DO 240 J=1,M 240 STD(J)=STD(J)/FN L=-M DO 250 I=1,M L=L+M+1 250 B(I)=RX(L) RETURN END SUBROUTINE CORRE !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE DATA(M,D) RETURN END SUBROUTINE DATA !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE LOC(I,J,IR,N,M,MS) IX=I JX=J IF (MS-1)10,20,30 10 IRX=N*(JX-1)+IX GOTO 36 20 IF(IX-JX)22,24,24 22 IRX=IX+(JX*JX-JX)/2 GOTO 36 24 IRX=JX+(IX*IX-IX)/2 GOTO 36 30 IRX=0 IF(IX-JX)36,32,36 32 IRX=IX 36 IR=IRX RETURN END SUBROUTINE LOC !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE MSTR(A,R,N,MSA,MSR) REAL A(1),R(1) DO 20 I=1,N DO 20 J=1,N IF(MSR)5,10,5 5 IF(I-J)10,10,20 10 CALL LOC(I,J,IR,N,N,MSR) IF (IR)20,20,15 15 R(IR)=0 CALL LOC(I,J,IA,N,N,MSA) IF(IA)20,20,18 18 R(IR)=A(IA) 20 CONTINUE RETURN END SUBROUTINE MSTR !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE STEP (NMax, MMax, M,N,D,XBAR,IDX,PCT,NSTEP,ANS,L,B,S,T,LL,IER, & NSTEPS,LS,KSTEP,BS,ANSS) ! Ten,NSTEPS,LS,KSTEP,BS,ANSS) ! INCLUDE 'HQTB.INC' ! DUNG LUONG MAU VA SO BIEN CUC DAI Integer :: NMax, MMax REAL D(M*M),XBAR(M), ANS(11), B(M),S(M), T(M), LL(M) INTEGER IDX(M), L(M), NSTEP(5) ! Common OUT: INTEGER :: NSTEPS(MMAX,5), LS(MMAX,MMAX), KSTEP REAL :: BS(MMAX,MMAX), ANSS(MMAX,11) ! Common VARN: ! Character (Len=8) :: Ten(MMAX) ! Variable's name !print*,"Calling STEP" IER=0 ONM=N-1 NFO=0 NSTEP(3)=0 ANS(3)=0 ANS(4)=0 NSTOP=0 DO 30 I=1,M LL(I)=1 IF(IDX(I))30,30,10 10 IF(IDX(I)-2)15,20,25 15 NFO=NFO+1 IDX(NFO)=I GOTO 30 20 NSTEP(3)=NSTEP(3)+1 LL(I)=-1 GOTO 30 25 MY=I NSTEP(1)=MY LY=M*(MY-1) LYP=LY+MY ANS(5)=D(LYP) 30 CONTINUE NSTEP(2)=NFO MX=M-NSTEP(3)-1 DO 140 NL=1,MX RD=0 IF(NL-NFO) 35,35,55 35 DO 50 I=1,NFO K=IDX(I) IF(LL(K))50,50,40 40 LYP=LY+K IP=M*(K-1)+K RE=D(LYP)**2/D(IP) IF(RD-RE) 45,50,50 45 RD=RE NEW=K 50 CONTINUE GOTO 75 55 DO 70 I=1,M IF(I-MY) 60,70,60 60 IF(LL(I))70,70,62 62 LYP=LY+I IP=M*(I-1)+I RE=D(LYP)**2/D(IP) IF(RD-RE) 64,70,70 64 RD=RE NEW=I 70 CONTINUE 75 IF(RD) 77,77,76 76 IF(ANS(5)-(ANS(3)+RD))77,77,78 77 IER=1 GOTO 150 78 RE=RD/ANS(5) IF(RE-PCT) 150,80,80 80 LL(NEW)=0 L(NL)=NEW ANS(1)=RD ANS(2)=RE ANS(3)=ANS(3)+RD ANS(4)=ANS(4)+RE NSTEP(4)=NL NSTEP(5)=NEW 85 ANS(6)=SQRT(ANS(4)) RD=NL RE=ONM-RD RE=(ANS(5)-ANS(3))/RE ANS(7)=(ANS(3)/RD)/RE 90 CONTINUE IF (RE.GE.0) THEN ANS(8)=SQRT(RE) ELSE ANS(8)=-999 ENDIF ! 90 ANS(8)=SQRT(RE) IP=M*(NEW-1)+NEW RD=D(IP) LYP=NEW-M DO 100 J=1,M LYP=LYP+M IF(LL(J))100,94,97 94 IF(J-NEW) 96,98,96 96 IJ=M*(J-1)+J D(IJ)=D(IJ)+D(LYP)*D(LYP)/RD 97 D(LYP)=D(LYP)/RD GOTO 100 98 D(IP)=1/RD 100 CONTINUE LYP=LY+NEW B(NL)=D(LYP) IF(NL-1)112,112,105 105 ID=NL-1 DO 110 J=1,ID IJ=NL-J KK=L(IJ) LYP=LY+KK B(IJ)=D(LYP) DO 110 K=1,J IK=NL-K+1 MK=L(IK) LYP=M*(MK-1)+KK 110 B(IJ)=B(IJ)-D(LYP)*B(IK) 112 ANS(9)=XBAR(MY) DO 115 I=1,NL KK=L(I) ANS(9)=ANS(9)-B(I)*XBAR(KK) IJ=M*(KK-1)+KK 114 S(I)=ANS(8)*SQRT(D(IJ)) 115 T(I)=B(I)/S(I) IP=M*(NEW-1) DO 130 I=1,M IJ=I-M IK=NEW-M IP=IP+1 IF(LL(I)) 130,130,120 120 DO 126 J=1,M IJ=IJ+M IK=IK+M IF(LL(J)) 126,122,122 122 IF(J-NEW)124,126,124 124 D(IJ)=D(IJ)-D(IP)*D(IK) 126 CONTINUE D(IP)=D(IP)/(-RD) 130 CONTINUE RD=N-NSTEP(4) RD=ONM/RD 132 ANS(10)=SQRT(1-(1-ANS(6)**2)*RD) 134 IF(RD.GE.0) GOTO 905 ANS(11)=999. GOTO 906 905 ANS(11)=ANS(8)*SQRT(RD) 906 CONTINUE CALL INHQ (NMax, MMax, NSTEP,ANS,L,B,S,T,NSTOP, & NSTEPS, LS, KSTEP, BS, ANSS) IF (NSTEP(4).GE.9) NSTOP = 1 ! LAY TOI DA 9 BUOC !print*,"After Calling STEP" IF (NSTOP) 140,140,150 140 CONTINUE 150 RETURN END SUBROUTINE STEP !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE INHQ (NMax, MMax, NSTEP,ANS,L,B,S,TT,NSTOP, & NSTEPS, LS, KSTEP, BS, ANSS) ! For outputs ! Ten, NSTEPS, LS, KSTEP, BS, ANSS) ! For outputs ! INCLUDE 'HQTB.INC' ! DUNG LUONG MAU VA SO BIEN CUC DAI Integer :: NMax, MMax REAL :: ANS(11), B(MMAX),S(MMAX),TT(MMAX) ! Character *8 Ten(MMAX) ! Common /VARN/ Ten INTEGER :: NSTEP(5), L(MMAX) ! INTEGER NSTEPS(MMAX,5), LS(MMAX,MMAX),KSTEP ! REAL BS(MMAX,MMAX), ANSS(MMAX,11) ! Common OUT: INTEGER :: NSTEPS(MMAX,5), LS(MMAX,MMAX), KSTEP, i, k REAL :: BS(MMAX,MMAX), ANSS(MMAX,11) ! Common VARN: ! Character (Len=8) :: Ten(MMAX) ! Variable's name ! COMMON /OUT/ NSTEPS,ANSS,LS,BS,KSTEP !print*,"Calling INHQ: KSTEP", KSTEP KSTEP=KSTEP+1 DO I=1,5 NSTEPS(KSTEP,I)=NSTEP(I) !print*,"1: ",i ENDDO DO I=1,11 ANSS(KSTEP,I)=ANS(I) !print*,"2: ",i ENDDO DO I=1,MMAX LS(KSTEP,I)=L(I) BS(KSTEP,I)=B(I) ENDDO K=NSTEP(4) !print*,"K=",K RETURN END SUBROUTINE INHQ !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC