Power-Programmierung/CD1.mdf/fortran/library/ssp/statdscr/discr.for

Text File | 1985-11-29 | 6KB | 191 lines

C
C     ..................................................................
C
C        SUBROUTINE DISCR
C
C        PURPOSE
C           COMPUTE A SET OF LINEAR FUNCTIONS WHICH SERVE AS INDICES
C           FOR CLASSIFYING AN INDIVIDUAL INTO ONE OF SEVERAL GROUPS.
C           NORMALLY THIS SUBROUTINE IS USED IN THE PERFORMANCE OF
C           DISCRIMINANT ANALYSIS.
C
C        USAGE
C           CALL DISCR (K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)
C
C        DESCRIPTION OF PARAMETERS
C           K     - NUMBER OF GROUPS. K MUST BE GREATER THAN ONE.
C           M     - NUMBER OF VARIABLES
C           N     - INPUT VECTOR OF LENGTH K CONTAINING SAMPLE SIZES OF
C                   GROUPS.
C           X     - INPUT VECTOR CONTAINING DATA IN THE MANNER EQUIVA-
C                   LENT TO A 3-DIMENSIONAL FORTRAN ARRAY, X(1,1,1),
C                   X(2,1,1), X(3,1,1), ETC.  THE FIRST SUBSCRIPT IS
C                   CASE NUMBER, THE SECOND SUBSCRIPT IS VARIABLE NUMBER
C                   AND THE THIRD SUBSCRIPT IS GROUP NUMBER.  THE
C                   LENGTH OF VECTOR X IS EQUAL TO THE TOTAL NUMBER OF
C                   DATA POINTS, T*M, WHERE T = N(1)+N(2)+...+N(K).
C           XBAR  - INPUT MATRIX (M X K) CONTAINING MEANS OF M VARIABLES
C                   IN K GROUPS
C           D     - INPUT MATRIX (M X M) CONTAINING THE INVERSE OF
C                   POOLED DISPERSION MATRIX.
C           CMEAN - OUTPUT VECTOR OF LENGTH M CONTAINING COMMON MEANS.
C           V     - OUTPUT VARIABLE CONTAINING GENERALIZED MAHALANOBIS
C                   D-SQUARE.
C           C     - OUTPUT MATRIX (M+1 X K) CONTAINING THE COEFFICIENTS
C                   OF DISCRIMINANT FUNCTIONS.  THE FIRST POSITION OF
C                   EACH COLUMN (FUNCTION) CONTAINS THE VALUE OF THE
C                   CONSTANT FOR THAT FUNCTION.
C           P     - OUTPUT VECTOR CONTAINING THE PROBABILITY ASSOCIATED
C                   WITH THE LARGEST DISCRIMINANT FUNCTIONS OF ALL CASES
C                   IN ALL GROUPS.  CALCULATED RESULTS ARE STORED IN THE
C                   MANNER EQUIVALENT TO A 2-DIMENSIONAL AREA (THE
C                   FIRST SUBSCRIPT IS CASE NUMBER, AND THE SECOND
C                   SUBSCRIPT IS GROUP NUMBER).  VECTOR P HAS LENGTH
C                   EQUAL TO THE TOTAL NUMBER OF CASES, T (T = N(1)+N(2)
C                   +...+N(K)).
C           LG    - OUTPUT VECTOR CONTAINING THE SUBSCRIPTS OF THE
C                   LARGEST DISCRIMINANT FUNCTIONS STORED IN VECTOR P.
C                   THE LENGTH OF VECTOR LG IS THE SAME AS THE LENGTH
C                   OF VECTOR P.
C
C        REMARKS
C           THE NUMBER OF VARIABLES MUST BE GREATER THAN OR EQUAL TO
C           THE NUMBER OF GROUPS.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           REFER TO 'BMD COMPUTER PROGRAMS MANUAL', EDITED BY W. J.
C           DIXON, UCLA, 1964, AND T. W. ANDERSON, 'INTRODUCTION TO
C           MULTIVARIATE STATISTICAL ANALYSIS', JOHN WILEY AND SONS,
C           1958, SECTION 6.6-6.8.
C
C     ..................................................................
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)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION XBAR,D,CMEAN,V,C,SUM,P,PL
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  EXP IN STATEMENT
C        250 MUST BE CHANGED TO DEXP.
C
C        ...............................................................
C
C     CALCULATE COMMON MEANS
C
      N1=N(1)
      DO 100 I=2,K
  100 N1=N1+N(I)
      FNT=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)/FNT
C
C     CALCULATE GENERALIZED MAHALANOBIS D SQUARE
C
      L=0
      DO 140 I=1,K
      DO 140 J=1,M
      L=L+1
  140 C(L)=XBAR(L)-CMEAN(J)
      V=0.0
      L=0
      DO 160 J=1,M
      DO 160 I=1,M
      N1=I-M
      N2=J-M
      SUM=0.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
C
C     CALCULATE THE COEFFICIENTS OF DISCRIMINANT FUNCTIONS
C
      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.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.0)
      DO 190 I=1,M
      N1=I-M
      IQ=IQ+1
      C(IQ)=0.0
      DO 190 J=1,M
      N1=N1+M
  190 C(IQ)=C(IQ)+D(N1)*P(J)
C
C     FOR EACH CASE IN EACH GROUP, CALCULATE..
C
C        DISCRIMINANT FUNCTIONS
C
      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
C
C        THE LARGEST DISCRIMINANT FUNCTION
C
      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
C
C        PROBABILITY ASSOCIATED WITH THE LARGEST DISCRIMINANT FUNCTION
C
      PL=0.0
      DO 250 J=1,K
  250 PL=PL+ EXP(XBAR(J)-SUM)
      N1=N1+1
      LG(N1)=L
  260 P(N1)=1.0/PL
  270 LBASE=LBASE+NN*M
C
      RETURN
      END