C Chuong trinh xac dinh ham phan lop C Day la chuong trinh chinh, can phai C co 3 chuong trinh con la C DMATX, DISCR, MINV C Program phanlop PARAMETER (KMAX=5,MMAX=20,NMAX=200) DIMENSION N(KMAX),CMEAN(MMAX),XBAR(KMAX*MMAX), & C((MMAX+1)*KMAX),D(MMAX*MMAX) DIMENSION P(NMAX*KMAX),LG(NMAX*KMAX),X(NMAX*KMAX*MMAX) ! dimension N(5),CMEAN(10),XBAR(50),C(55),D(100) ! DIMENSION P(250),LG(250),X(2500) character *12 fsl,fra C Kich thuoc duoc mo ta nhu sau C Cua N : >= so nhom K C Cua CMEAN : >= so bien M C Cua XBAR : >= tich so M x K C Cua C : >= tich (M+1) x K C Cua D : >= tich M x M C cua P,LG : >= kich thuoc chung cua tap mau K nhom, T C (T=N(1)+N(2)+...+N(k)) C Cua X : >= tich T x M C================================================================== C...CAU TRUC FILE SO LIEU NHU SAU: C...DONG 1: GOM K+1 SO NGUYEN DUONG, XEP THEO THU TU LA C K (SO NHOM), M (SO BIEN), N1,.., NK (DUNG LUONG C MAU CUA TUNG NHOM) C...DONG 2 TRO DI GOM N1+...+NK DONG, MOI DONG GOM M SO LA C SO LIEU QUAN TRAC CUA M BIEN. CAC DONG DUOC XEP THEO THU C TU TU TREN XUONG DUOI TUONG UNG TU LOP 1 (N1) DEN LOP K (NK) C write(*,'(a\)')' Cho ten File SL : ' read(*,'(a)')fsl write(*,'(a\)')' Cho ten File Ra : ' read(*,'(a)') fra open(1,File=fsl,status='old') open(3,File=fra,status='New') do 400 i=1,25 400 write(*,*) 100 read(1,*) K,M,(N(i),i=1,K) write(3,2) K,M 2 format(' Chuong trinh Ham phan lop '/'So nhom',I5/'So bien',I5) C do 110 i=1,k C 110 write(3,'(2i6)') I,N(i) L=0 do 130 i=1,k n1=n(i) do 120 j=1,n1 read(1,*)(cmean(ij),ij=1,M) L=L+1 n2=L-n1 do 120 ij=1,M n2=n2+n1 120 x(n2)=cmean(ij) 130 L=n2 call DMATX(k,m,n,x,xbar,d,cmean) L=0 do 150 i=1,k do 140 j=1,m L=L+1 140 cmean(j)=xbar(L) C 150 write(3,'('' Nhom'',i4,'' Tbinh ''/8f15.5)') I,(cmean(j),j=1,M) C write(3,'('' Ma tran hiep phuong sai '')') 150 continue do 170 i=1,m L=i-m do 160 j=1,m L=L+m 160 cmean(j)=D(L) C 170 write(3,8)i,(CMEAN(j),j=1,M) 170 continue 8 format(' Hang ',i3/(8f15.5)) Call MINV(D,M,DET,CMEAN,C) Call DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG) C write(3,9)(CMEAN(I),I=1,M) 9 format(' Trung binh chung '/(8f15.5)) C write(3,10) V 10 format(' Hang so D-Mahalanov ',f15.5) n1=1 n2=M+1 do 180 i=1,k write(3,11) I,(C(J),J=N1,N2) 11 FORMAT(' Ham phan lop ',i3/6x,'Hang so Cac he so'/9f12.3) n1=n1+(M+1) 180 n2=n2+(M+1) C write(3,12) 12 format('Evoluation of classification function' \' for each observation') n1=1 n2=N(1) do 210 i=1,k write(3,13) I 13 format(' Nhom ',i3) C 13 format(' Nhom ',i3/14x,' Xac suat ung voi ham So thu tu'/ C \' Quan trac phan lop lon nhat ham lon nhat') L=0 do 190 j=n1,n2 L=L+1 190 write(3,14) L,LG(j) 14 format(1x,i4,2x,i6) if (i-k) 200,300,300 200 n1=n1+N(i) n2=n2+N(i+1) 210 continue 300 continue end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C --------- CAC CHUONG TRINH CON ------------- C SUBROUTINE DMATX(K,M,N,X,XBAR,D,CMEAN) dimension N(1),X(1),XBAR(1),D(1),CMEAN(1) MM=M*M do 100 i=1,MM 100 D(i)=0 N4=0 L=0 LM=0 do 160 NG=1,K N1=N(NG) FN=N1 do 130 j=1,M LM=LM+1 XBAR(LM)=0 do 120 i=1,N1 L=L+1 120 XBAR(LM)=XBAR(LM)+X(L) 130 XBAR(LM)=XBAR(LM)/FN LMEAN=LM-M do 150 i=1,N1 LL=N4+i-N1 do 140 j=1,M LL=LL+N1 N2=LMEAN+j 140 CMEAN(j)=X(LL)-XBAR(N2) LL=0 do 150 j=1,M do 150 jj=1,M LL=LL+1 150 D(LL)=D(LL)+CMEAN(J)*CMEAN(JJ) 160 N4=N4+N1*M LL=-K DO 170 I=1,K 170 LL=LL+N(I) FN=LL do 180 I=1,MM 180 D(I)=D(I)/FN RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG) dimension N(1),X(1),XBAR(1),D(1),CMEAN(1),C(1),P(1),LG(1) N1=N(1) do 100 i=2,k 100 N1=N1+N(I) FN1=N1 do 110 i=1,k 110 P(i)=N(I) do 130 i=1,M CMEAN(I)=0 N1=I-M do 120 j=1,k N1=N1+M 120 CMEAN(I)=CMEAN(I)+P(J)*XBAR(N1) 130 CMEAN(I)=CMEAN(I)/FN1 L=0 do 140 i=1,k do 140 j=1,M L=L+1 140 C(L)=XBAR(L)-CMEAN(J) V=0 L=0 do 160 j=1,M do 160 i=1,M N1=I-M N2=J-M Sum=0 do 150 ij=1,k N1=N1+M N2=N2+M 150 Sum=Sum+P(IJ)*C(N1)*C(N2) L=L+1 160 V=V+D(L)*Sum N2=0 do 190 KA=1,k do 170 i=1,M N2=N2+1 170 P(I)=XBAR(N2) IQ=(M+1)*(KA-1)+1 Sum=0 do 180 j=1,M N1=J-M do 180 L=1,M N1=N1+M 180 Sum=Sum+D(N1)*P(J)*P(L) C(IQ)=-(Sum/2) do 190 i=1,M N1=I-M IQ=IQ+1 C(IQ)=0 do 190 j=1,M N1=N1+M 190 C(IQ)=C(IQ)+D(N1)*P(J) LBASE=0 N1=0 do 270 KG=1,K NN=N(KG) do 260 i=1,NN L=I-NN+LBASE do 200 j=1,M L=L+NN 200 D(J)=X(L) N2=0 do 220 KA=1,K N2=N2+1 Sum=C(N2) do 210 j=1,M N2=N2+1 210 Sum=Sum+C(N2)*D(J) 220 XBAR(KA)=Sum L=1 Sum=XBAR(1) do 240 j=2,k IF(Sum-XBAR(j)) 230,240,240 230 L=J Sum=XBAR(J) 240 CONTINUE PL=0 do 250 j=1,k 250 PL=PL+EXP(XBAR(J)-Sum) N1=N1+1 LG(N1)=L 260 P(N1)=1/PL 270 LBASE=LBASE+NN*M RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE MINV(A,N,D,L,M) DIMENSION A(1),L(1),M(1) D=1. NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE J=L(K) IF(J-K)35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI)=HOLD 35 I=M(K) IF(I-K)45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) 40 A(JI)=HOLD 45 IF(BIGA)48,46,48 46 D=0. RETURN 48 DO 55 I=1,N IF(I-K)50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K)60,65,60 60 IF(J-K)62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K)70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE D=D*BIGA A(KK)=1./BIGA 80 CONTINUE K=N 100 K=(K-1) IF(K)150,150,105 105 I=L(K) IF(I-K)120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) 110 A(JI)=HOLD 120 J=M(K) IF(J-K)100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) 130 A(JI)=HOLD GOTO 100 150 RETURN END