home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / statcorr / misr.for < prev    next >
Text File  |  1985-12-26  |  7KB  |  231 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE MISR
  5. C
  6. C     PURPOSE
  7. C        COMPUTE MEANS, STANDARD DEVIATIONS, SKEWNESS AND KURTOSIS,
  8. C        CORRELATION COEFFICIENTS, REGRESSION COEFFICIENTS, AND
  9. C        STANDARD ERRORS OF REGRESSION COEFFICIENTS WHEN THERE ARE
  10. C        MISSING DATA POINTS.  THE USER IDENTIFIES THE MISSING DATA
  11. C        BY MEANS OF A NUMERIC CODE.  THOSE VALUES HAVING THIS CODE
  12. C        ARE SKIPPED IN COMPUTING THE STATISTICS.  IN THE CASE OF THE
  13. C        CORRELATION COEFFICIENTS, ANY PAIR OF VALUES ARE SKIPPED IF
  14. C        EITHER ONE OF THEM ARE MISSING.
  15. C
  16. C     USAGE
  17. C        CALL MISR (NO,M,X,CODE,XBAR,STD,SKEW,CURT,R,N,A,B,S,IER)
  18. C
  19. C     DESCRIPTION OF PARAMETERS
  20. C        NO     - NUMBER OF OBSERVATIONS
  21. C        M     - NUMBER OF VARIABLES
  22. C        X     - INPUT DATA MATRIX OF SIZE NO X M.
  23. C        CODE - INPUT VECTOR OF LENGTH M, WHICH CONTAINS A NUMERIC
  24. C           MISSING DATA CODE FOR EACH VARIABLE. ANY OBSERVATION
  25. C           FOR A GIVEN VARIABLE HAVING A VALUE EQUAL TO THE CODE
  26. C           WILL BE DROPPED FOR THE COMPUTATIONS.
  27. C        XBAR - OUTPUT VECTOR OF LENGTH M CONTAINING MEANS
  28. C        STD  - OUTPUT VECTOR OF LENGTH M CONTAINING STANDARD DEVI-
  29. C           ATIONS
  30. C        SKEW - OUTPUT VECTOR OF LENGTH M CONTAINING SKEWNESS
  31. C        CURT - OUTPUT VECTOR OF LENGTH M CONTAINING KURTOSIS
  32. C        R     - OUTPUT MATRIX OF PRODUCT-MOMENT CORRELATION
  33. C           COEFFICIENTS.  THIS WILL BE THE UPPER TRIANGULAR
  34. C           MATRIX ONLY, SINCE THE M X M MATRIX OF COEFFICIENTS
  35. C           IS SYMMETRIC. (STORAGE MODE 1)
  36. C        N     - OUTPUT MATRIX OF NUMBER OF PAIRS OF OBSERVATIONS USED
  37. C           IN COMPUTING THE CORRELATION COEFFICIENTS.  ONLY THE
  38. C           UPPER TRIANGULAR PORTION OF THE MATRIX IS GIVEN.
  39. C           (STORAGE MODE 1)
  40. C        A     - OUTPUT MATRIX (M BY M)  CONTAINING INTERCEPTS OF
  41. C           REGRESSION LINES (A) OF THE FORM Y=A+BX.  THE FIRST
  42. C           SUBSCRIPT OF THIS MATRIX REFERS TO THE INDEPENDENT
  43. C           VARIABLE AND THE SECOND TO THE DEPENDENT VARIABLE.
  44. C           FOR EXAMPLE, A(1,3) CONTAINS THE INTERCEPT OF THE
  45. C           REGRESSION LINE FOR TWO VARIABLES WHERE VARIABLE 1
  46. C           IS INDEPENDENT AND VARIABLE 3 IS DEPENDENT.    NOTE
  47. C           THAT MATRIX A IS STORED IN A VECTOR FORM.
  48. C        B     - OUTPUT MATRIX (M BY M)  CONTAINING REGRESSION
  49. C           COEFFICIENTS (B) CORRESPONDING TO THE VALUES OF
  50. C           INTERCEPTS CONTAINED IN THE OUTPUT MATRIX A.
  51. C        S     - OUTPUT MATRIX (M BY M)  CONTAINING STANDARD ERRORS
  52. C           OF REGRESSION COEFFICIENTS CORRESPONDING TO THE
  53. C           COEFFICIENTS CONTAINED IN THE OUTPUT MATRIX B.
  54. C        IER  - 0, NO ERROR.
  55. C           1, IF NUMBER OF NON-MISSING DATA ELEMENTS FOR J-TH
  56. C              VARIABLE IS TWO OR LESS.    IN THIS CASE, STD(J),
  57. C              SKEW(J), AND CURT(J) ARE SET TO 10**75.  ALL
  58. C              VALUES OF R, A, B, AND S RELATED TO THIS VARIABLE
  59. C              ARE ALSO SET TO 10**75.
  60. C           2, IF VARIANCE OF J-TH VARIABLE IS LESS THAN
  61. C              10**(-20).  IN THIS CASE, STD(J), SKEW(J), AND
  62. C              CURT(J) ARE SET TO 10**75.  ALL VALUES OF R, A,
  63. C              B, AND S RELATED TO THIS VARIABLE ARE ALSO SET TO
  64. C              10**75.
  65. C
  66. C     REMARKS
  67. C        THIS SUBROUTINE CANNOT DISTINGUISH A BLANK AND A ZERO.
  68. C        THEREFORE, IF A BLANK IS SPECIFIED AS A MISSING DATA CODE IN
  69. C        INPUT CARDS, IT WILL BE TREATED AS 0 (ZERO).
  70. C
  71. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  72. C        NONE
  73. C
  74. C     METHOD
  75. C        LEAST SQUARES REGRESSION LINES AND PRODUCT-MOMENT CORRE-
  76. C        LATION COEFFICIENTS ARE COMPUTED.
  77. C
  78. C     ..................................................................
  79. C
  80.       SUBROUTINE MISR (NO,M,X,CODE,XBAR,STD,SKEW,CURT,R,N,A,B,S,IER)
  81. C
  82.       DIMENSION X(1),CODE(1),XBAR(1),STD(1),SKEW(1),CURT(1),R(1),N(1)
  83.       DIMENSION A(1),B(1),S(1)
  84. C
  85. C     COMPUTE MEANS
  86. C
  87.       IER=0
  88.       L=0
  89.       DO 20 J=1,M
  90.       FN=0.0
  91.       XBAR(J)=0.0
  92.       DO 15 I=1,NO
  93.       L=L+1
  94.       IF(X(L)-CODE(J)) 12, 15, 12
  95.    12 FN=FN+1.0
  96.       XBAR(J)=XBAR(J)+X(L)
  97.    15 CONTINUE
  98.       IF(FN) 16, 16, 17
  99.    16 XBAR(J)=0.0
  100.       GO TO 20
  101.    17 XBAR(J)=XBAR(J)/FN
  102.    20 CONTINUE
  103. C
  104. C     SET-UP WORK AREAS AND TEST WHETHER DATA IS MISSING
  105. C
  106.       L=0
  107.       DO 55 J=1,M
  108.       LJJ=NO*(J-1)
  109.       SKEW(J)=0.0
  110.       CURT(J)=0.0
  111.       KI=M*(J-1)
  112.       KJ=J-M
  113.       DO 54 I=1,J
  114.       KI=KI+1
  115.       KJ=KJ+M
  116.       SUMX=0.0
  117.       SUMY=0.0
  118.       TI=0.0
  119.       TJ=0.0
  120.       TII=0.0
  121.       TJJ=0.0
  122.       TIJ=0.0
  123.       NIJ=0
  124.       LI=NO*(I-1)
  125.       LJ=LJJ
  126.       L=L+1
  127.       DO 38 K=1,NO
  128.       LI=LI+1
  129.       LJ=LJ+1
  130.       IF(X(LI)-CODE(I)) 30, 38, 30
  131.    30 IF(X(LJ)-CODE(J)) 35, 38, 35
  132. C
  133. C     BOTH DATA ARE PRESENT
  134. C
  135.    35 XX=X(LI)-XBAR(I)
  136.       YY=X(LJ)-XBAR(J)
  137.       TI=TI+XX
  138.       TII=TII+XX**2
  139.       TJ=TJ+YY
  140.       TJJ=TJJ+YY**2
  141.       TIJ=TIJ+XX*YY
  142.       NIJ=NIJ+1
  143.       SUMX=SUMX+X(LI)
  144.       SUMY=SUMY+X(LJ)
  145.       IF(I-J) 38, 37, 37
  146.    37 SKEW(J)=SKEW(J)+YY**3
  147.       CURT(J)=CURT(J)+YY**4
  148.    38 CONTINUE
  149. C
  150. C     COMPUTE SUM OF CROSS-PRODUCTS OF DEVIATIONS
  151. C
  152.       IF(NIJ) 40, 40, 39
  153.    39 FN=NIJ
  154.       R(L)=TIJ-TI*TJ/FN
  155.       N(L)=NIJ
  156.       TII=TII-TI*TI/FN
  157.       TJJ=TJJ-TJ*TJ/FN
  158. C
  159. C     COMPUTE STANDARD DEVIATION, SKEWNESS, AND KURTOSIS
  160. C
  161.    40 IF(I-J) 47, 41, 47
  162.    41 IF(NIJ-2) 42,42,43
  163.    42 IER=1
  164.       R(L)=1.0E38
  165.       A(KI)=1.0E38
  166.       B(KI)=1.0E38
  167.       S(KI)=1.0E38
  168.       GO TO 45
  169. C
  170.    43 STD(J)=R(L)
  171.       R(L)=1.0
  172.       A(KI)=0.0
  173.       B(KI)=1.0
  174.       S(KI)=0.0
  175. C
  176.       IF(STD(J)-(1.0E-20)) 44,44,46
  177.    44 IER=2
  178.    45 STD(J)=1.0E38
  179.       SKEW(J)=1.0E38
  180.       CURT(J)=1.0E38
  181.       GO TO 55
  182. C
  183.    46 WORK=STD(J)/FN
  184.       SKEW(J)=(SKEW(J)/FN)/(WORK*SQRT(WORK))
  185.       CURT(J)=((CURT(J)/FN)/WORK**2)-3.0
  186.       STD(J)=SQRT(STD(J)/(FN-1.0))
  187.       GO TO 55
  188. C
  189. C     COMPUTE REGRESSION COEFFICIENTS
  190. C
  191.    47 IF(NIJ-2) 48,48,50
  192.    48 IER=1
  193.    49 R(L)=1.0E38
  194.       A(KI)=1.0E38
  195.       B(KI)=1.0E38
  196.       S(KI)=1.0E38
  197.       A(KJ)=1.0E38
  198.       B(KJ)=1.0E38
  199.       S(KJ)=1.0E38
  200.       GO TO 54
  201. C
  202.    50 IF(TII-(1.0E-20)) 52,52,51
  203.    51 IF(TJJ-(1.0E-20)) 52,52,53
  204.    52 IER=2
  205.       GO TO 49
  206. C
  207.    53 SUMX=SUMX/FN
  208.       SUMY=SUMY/FN
  209.       B(KI)=R(L)/TII
  210.       A(KI)=SUMY-B(KI)*SUMX
  211.       B(KJ)=R(L)/TJJ
  212.       A(KJ)=SUMX-B(KJ)*SUMY
  213. C
  214. C        COMPUTE CORRELATION COEFFICIENTS
  215. C
  216.       R(L)=R(L)/(SQRT(TII)*SQRT(TJJ))
  217. C
  218. C        COMPUTE STANDARD ERRORS OF REGRESSION COEFFICIENTS
  219. C
  220.       RR=R(L)**2
  221.       SUMX=(TJJ-TJJ*RR)/(FN-2)
  222.       S(KI)=SQRT(SUMX/TII)
  223.       SUMY=(TII-TII*RR)/(FN-2)
  224.       S(KJ)=SQRT(SUMY/TJJ)
  225. C
  226.    54 CONTINUE
  227.    55 CONTINUE
  228. C
  229.       RETURN
  230.       END
  231. A(KJ)=