Power-Programmierung/CD1.mdf/fortran/library/ssp/statmisc/tetra.for

Text File | 1985-12-31 | 5.0 KB | 171 lines

C
C     ..................................................................
C
C     SUBROUTINE TETRA
C
C     PURPOSE
C        COMPUTE A TETRACHORIC CORRELATION COEFFICIENT BETWEEN TWO
C        VARIABLES WHERE DATA IN BOTH VARIABLES HAVE BEEN REDUCED
C        ARTIFICIALLY TO TWO CATEGORIES.
C
C     USAGE
C        CALL TETRA (N,U,V,HU,HV,R,RS,IE)
C
C     DESCRIPTION OF PARAMETERS
C        N  - NUMBER OF OBSERVATIONS
C        U  - INPUT VECTOR OF LENGTH N CONTAINING THE FIRST VARIABLE
C         REDUCED TO TWO CATEGORIES
C        V  - INPUT VECTOR OF LENGTH N CONTAINING THE SECOND VARIABLE
C         REDUCED TO TWO CATEGORIES
C        HU - INPUT NUMERICAL CODE INDICATING THE HIGHER CATEGORY OF
C         THE FIRST VARIABLE.  IF ANY VALUE OF VARIABLE U IS
C         EQUAL TO OR GREATER THAN HU, IT WILL BE CLASSIFIED AS
C         THE HIGHER CATEGORY, OTHERWISE AS THE LOWER CATEGORY.
C        HV - SAME AS HU EXCEPT THAT HV IS FOR THE SECOND VARIABLE.
C        R  - TETRACHORIC CORRELATION COMPUTED
C        RS - STANDARD ERROR OF TETRACHORIC CORRELATION COMPUTED
C        IE - ERROR CODE
C         0 - NO ERROR
C         1 - UNABLE TO COMPUTE A TETRACHORIC CORRELATION DUE TO
C             THE FACT THAT AT LEAST ONE CELL SHOWS ZERO FRE-
C             QUENCY IN THE 2X2 CONTINGENCY TABLE CONSTRUCTED
C             FROM INPUT DATA.  IN THIS CASE, R AND RS ARE SET
C             TO 10**75.  (SEE GUILFORD, 1956)
C         2 - THE ROOT SOLVER GIVES MULTIPLE ROOTS, OR NO ROOTS,
C             R, IN THE INTERVAL (-1,1) INCLUSIVE. R AND RS ARE
C             SET TO 10**75.
C         3 - UNABLE TO COMPUTE A SATISFACTORY VALUE OF TETRA-
C             CHORIC CORRELATION USING NEWTON-RAPHSON METHOD OF
C             APPROXIMATION TO THE ROOT OF THE EQUATION.  R AND
C             RS ARE SET TO 10**75.  SEE SUBROUTINE POLRT ERROR
C             INDICATORS.
C         4 - HIGH ORDER COEFFICIENT OF THE POLYNOMIAL IS ZERO.
C             SEE SUBROUTINE POLRT ERROR INDICATORS.
C
C     REMARKS
C        VALUES OF VARIABLES U AND V MUST BE NUMERICAL, AND
C        ALPHABETIC AND SPECIAL CHARACTERS MUST NOT BE USED.
C        FOR A DEPENDABLE RESULT FOR TETRACHORIC CORRELATION,
C        IT IS RECOMMENDED THAT N BE AT LEAST 200 OR GREATER.
C
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C        NDTRI
C        POLRT--THIS POLYNOMIAL ROOT ROUTINE WAS SELECTED BECAUSE OF
C           ITS SMALL STORAGE REQUIREMENT.  OTHER SSP ROUTINES
C           WHICH COULD REPLACE POLRT ARE PRQD AND PRBM.  THEIR
C           USE WOULD REQUIRE MODIFICATION OF TETRA.
C
C     METHOD
C        REFER TO J. P. GUILFORD, 'FUNDAMENTAL STATISTICS IN PSYCHO-
C        LOGY AND EDUCATION', MCGRAW-HILL, NEW YORK, 1956, CHAPTER 13
C        AND W. P. ELDERTON, 'FREQUENCY CURVES AND CORRELATION' 4-TH
C        ED., CAMBRIDGE UNIVERSITY PRESS, 1953, CHAPTER 9.
C
C     ..................................................................
C
      SUBROUTINE TETRA (N,U,V,HU,HV,R,RS,IE)
C
      DIMENSION XCOF(8),COF(8),ROOTR(7),ROOTI(7)
      DIMENSION U(1),V(1)
      DOUBLE PRECISION X31,X32,X312,X322
C
C     CONSTRUCT A 2X2 CONTINGENCY TABLE
C
      A=0.0
      B=0.0
      C=0.0
      D=0.0
      DO 40 I=1,N
      IF(U(I)-HU) 10, 25, 25
   10 IF(V(I)-HV) 15, 20, 20
   15 D=D+1.0
      GO TO 40
   20 B=B+1.0
      GO TO 40
   25 IF(V(I)-HV) 30, 35, 35
   30 C=C+1.0
      GO TO 40
   35 A=A+1.0
   40 CONTINUE
C
C     TEST WHETHER ANY CELL IN THE CONTINGENCY TABLE IS ZERO.
C     IF SO, RETURN TO THE CALLING ROUTINE WITH R=0.0 AND IE=1.
C
      IE=0
      IF(A) 60, 60, 45
   45 IF(B) 60, 60, 50
   50 IF(C) 60, 60, 55
   55 IF(D) 60, 60, 70
   60 IE=1
      GO TO 86
C
C     COMPUTE P1, Q1, P2, AND Q2
C
   70 FN=N
      P1=(A+C)/FN
      Q1=(B+D)/FN
      P2=(A+B)/FN
      Q2=(C+D)/FN
C
C     FIND THE STANDARD NORMAL DEVIATES AT Q1 AND Q2, AND THE
C     ORDINATES AT THOSE POINTS
C
      CALL NDTRI (Q1,X1,Y1,ER)
      CALL NDTRI (Q2,X2,Y2,ER)
C
C     COMPUTE THE TETRACHORIC CORRELATION COEFFICIENT
C
      IF(X1) 76, 72, 76
   72 IF(X2) 76, 74, 76
   74 R=0.0
      GO TO 90
   76 XCOF(1)=-((A*D-B*C)/(Y1*Y2*FN*FN))
      XCOF(2)=1.0
      XCOF(3)=X1*X2/2.0
      XCOF(4)=(X1*X1-1.0)*(X2*X2-1.0)/6.0
      X31=DBLE(X1)
      X32=DBLE(X2)
      X312=X31**2
      X322=X32**2
      XCOF(5)=SNGL(X31*(X312-3.0D0)*X32*(X322-3.0D0)/24.0D0)
      XCOF(6)=SNGL((X312*(X312-6.0D0)+3.0D0)*(X322*(X322-6.0D0)+3.0D0)
     1          /120.0D0)
      XCOF(7)=SNGL(X31*(X312*(X312-10.0D0)+15.0D0)*X32*(X322*(X322-10.0
     1          D0)+15.0D0)/720.0D0)
      XCOF(8)=SNGL((((X312-15.0D0)*X312+45.0D0)*X312-15.0D0)*(((X322-
     1          15.0D0)*X322+45.0D0)*X322-15.0D0)/5040.0D0)
C
      CALL POLRT (XCOF,COF,7,ROOTR,ROOTI,IER)
C
      J=0
      IF(IER) 78, 78, 84
   78 DO 82 I=1,7
      IF(ABS(ROOTI(I))-.5*ABS(ROOTR(I))*1.0E-6)79,79,82
   79 R=ROOTR(I)
      IF(ABS(R)-1.0)81,81,80
   80 R=1.E38
      GO TO 82
   81 J=J+1
   82 CONTINUE
      IF(J-1)83,88,83
   83 IE=2
      GO TO 86
C
C     UNABLE TO COMPUTE R
C
   84 IE=IER
   86 R=1.0E38
      RS=R
      GO TO 100
   88 IF(R-1.0E38)90,83,83
C
C     STANDARD ERROR OF R=0.0
C
   90 RS= SQRT(P1*P2*Q1*Q2)/(Y1*Y2* SQRT(FN))
C
  100 RETURN
      END

COF,COF,7,ROOTR,ROOTI,IER)
C
      J=0