home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol269 / histnorm.for < prev    next >
Encoding:
Text File  |  1986-05-22  |  8.8 KB  |  306 lines
  1. C [HISTNORM.FOR of JUGCPM Vol.12]
  2. C
  3. C     ** HISTOGRAM AND NORMAL DISTRIBUTION **   [╛▓╖╠▐▌╠▀ ╔ ╣▌├▓]
  4. C
  5. C        Yoshio MONMA                                     84-05-08   
  6. C
  7. C        This program has been written in FORTRAN-80 and the output
  8. C     printer has been assumed to be EPSON's RP-80.  The input data 
  9. C     must be in FORT0X.DAT (X=6-10).  So, you have to prepare the 
  10. C     data file as follows:
  11. C
  12. C             REN FORT0X.DAT=HISTNORM.DAT
  13. C     or
  14. C             PIP FORT0X.DAT=HISTNORM.DAT
  15. C
  16. C        Also note that the data file (FORT0X.DAT) must be placed on the
  17. C     same disk with the .COM file.
  18. C
  19. C     * Reference *
  20. C       T. Haga & S. Hashimoto: "Fundamentals of Programming for 
  21. C       Statistical Analysis" (─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐),
  22. C       JUSE (╞»╢╖▐┌▌), 1980.
  23. C
  24. C    Nomenclature is as follows:
  25. C
  26. C     INPT  [I*1]  File reference no. for input data (6<=INPT<=10)
  27. C        N  [I*2]  Data size
  28. C     KDTA  [I*2]  Code to specify output of data (KDTA=1: No)
  29. C    TITLE  [R*4]  Title (15A4)
  30. C    VFRMT  [R*4]  Variable FORMAT for data (15A4)
  31. C     X(N)  [R*4]  Array of data
  32. C       KK  [I*1]  No. of classes
  33. C     XMIN  [R*4]  Min. of X(N)
  34. C        D  [R*4]  Width of a class
  35. C
  36. C     * List of input data *
  37. C       INPT (I5)    Standard input file ref. no. of F80 
  38. C ---------------------------------------------------------
  39. C                    6<= INPT <= 10
  40. C       Following data are input from FORTXX.DAT (XX=INPT)
  41. C          N (I5)     Data size
  42. C       KDTA (I5)     Data output (KDTA=1: No)
  43. C      TITLE (15A4)   Title
  44. C      VFRMT (15A4)   Variable format for the array of data X(N)
  45. C       X(N) (VFRMT)  Array of the data
  46. C         KK (I5)     No. of classes
  47. C       XMIN (F10.5)  Lower bound of the data
  48. C          D (F10.5)  Width of the class
  49. C ---------------------------------------------------------
  50. C       INPT (I5)     [Next data set]
  51. C
  52.       PROGRAM      HSTNRM
  53. C
  54.       INTEGER*1    SI,INPT,KDTA,KK
  55.       REAL*4       X(500),TITLE(15),VFRMT(15)
  56. C
  57.       DATA         ESC,BIGM/Z'1B',Z'4D'/
  58. C
  59.       WRITE(1,1000)
  60.  1000   FORMAT(1H ,5X,'** Histogram and Normal Distribution ',
  61.      1         '(HISTNORM) **'/)
  62.       WRITE(2,2000) ESC,BIGM
  63.  2000   FORMAT(1H ,2A1)
  64. C
  65.    10 WRITE(1,1010)
  66.  1010   FORMAT(1H ,8X,'Specify input file [FORT0X.DAT], X in I5) =')
  67.       READ(1,1020) INPT
  68.  1020   FORMAT(16I5)
  69.       IF (INPT.LT.6.OR.INPT.GT.10)       GOTO 10
  70.    20 READ(INPT,1020) N,KDTA
  71.       IF (N.LE.0) STOP
  72. C
  73.       WRITE(2,1000)
  74.       READ(INPT,6010) TITLE
  75.  6010   FORMAT(1X,15A4)
  76.       READ(INPT,6010) VFRMT
  77.       WRITE(1,1030) N
  78.  1030   FORMAT(1H0,5X,'DATA SIZE =',I4,', Reading data from DK...')
  79.       READ(INPT,VFRMT) (X(I),I=1,N)
  80.       READ(INPT,6020) KK,XMIN,D
  81.  6020   FORMAT(I5,2F10.5)
  82.       WRITE(2,2010) TITLE
  83.  2010   FORMAT(1H0,5X,15A4)
  84.       WRITE(2,2020) N,KDTA
  85.  2020   FORMAT(1H0,10X,'N =',I4,5X,'KDTA =',I2)
  86.       IF (KDTA.NE.1) WRITE(2,2030) (X(I),I=1,N)
  87.  2030   FORMAT(1H ,1P5E14.5)
  88. C
  89.       WRITE(1,1040)
  90.  1040   FORMAT(1H0,5X,'* Call STAT0 *')
  91.       CALL STAT0(X,N,SX,AVGX,SXX,S,V,SIGX,SKEW,FKURT,GEARY)
  92.       WRITE(2,2040) AVGX,SIGX,SKEW,FKURT,GEARY
  93.  2040   FORMAT(1H0,5X,'Avg. =',1PE15.5,', SD =',E15.5/6X,
  94.      1         'Skewness =',0PF8.4,', Kurtosis =',F8.4,
  95.      2         ',  Geary =',F8.4/)
  96. C
  97.       WRITE(1,1050)
  98.  1050   FORMAT(1H0,5X,'* Call HISTGM *')
  99.       CALL HISTGM(X,N,KK,XMIN,D)
  100. C
  101.       WRITE(1,1060)
  102.  1060   FORMAT(1H0,5X,'* Sorting...(SORT1)')
  103.       CALL SORT1(X,N)
  104. C
  105.       WRITE(2,2050)
  106.  2050   FORMAT(1H1)
  107.       WRITE(2,2010) TITLE
  108.       WRITE(1,1070)
  109.  1070   FORMAT(1H0,5X,'* Call NRPLOT *')
  110.       CALL NRPLOT(X,N)
  111.       WRITE(2,2050)
  112. C
  113.                               GOTO 20
  114.       END
  115. C
  116.       SUBROUTINE STAT0 (X,N,SX,AVGX,SXX,S,V,SIGX,SKEW,FKURT,GEARY)
  117. C
  118. C     ** Fundamental Statistics **
  119. C
  120.       DIMENSION X(1)
  121. C
  122.       SX = 0.0
  123.       SXX = 0.0
  124.       DO 10 I=1,N
  125.         SX = SX+X(I)
  126.         SXX = SXX+X(I)*X(I)
  127.    10 CONTINUE
  128.       FN = N
  129.       AVGX = SX/FN
  130.       S = SXX-SX*AVGX
  131.       V = S/(FN-1.0)
  132.       SIGX = SQRT(V)
  133.       S3 = 0.0
  134.       S4 = 0.0
  135.       GEARY = 0.0
  136.       DO 20 I=1,N
  137.          D = X(I) - AVGX
  138.          S3 = S3 + D*D*D
  139.          S4 = S4 + D*D*D*D
  140.          GEARY = GEARY + ABS(D)
  141.    20 CONTINUE
  142.       SKEW = S3/(FN*SIGX**3)
  143.       FKURT = S4/(FN*SIGX**4)
  144.       GEARY = GEARY/(FN*SIGX)
  145.       RETURN
  146.       END
  147. C
  148.       SUBROUTINE   HISTGM(X,II,KK,XMIN,D)
  149. C
  150. C     * Histogram *
  151. C
  152. C     Arguments:           Input        |    Output
  153. C      DTAID     [R*8]  Data id. (A8)   |  <No change> 
  154. C      X(II)     [R*4]  Array of data   |  <No change>     
  155. C         II     [I*2]  Data size       |  <No change>
  156. C         KK     [I*1]  No. of classes  |  <No change>
  157. C       XMIN     [R*4]  Min. of data    |  <No change>   
  158. C          D     [R*4]  Width of class  |  <No change> 
  159. C   
  160.       INTEGER*1    K,KK,KKM,H(51),HSPACE,HVRTL,HSTAR,HPLUS,N(20)
  161.       REAL*4       X(1)
  162. C
  163.       DATA         HSPACE,HVRTL,HSTAR,HPLUS /' ','|',Z'EC','+'/
  164. C
  165.       WRITE(2,200) KK,XMIN,D
  166.   200   FORMAT(1H0,5X,'No. of classes (KK) =',I3,', XMIN =',F10.3,
  167.      1         ', Width (D) =',F7.3/)
  168.       KKM = KK - 2
  169.       DO 10 K=1,KK
  170.          N(K) = 0
  171.    10 CONTINUE
  172.       DO 15 I=1,II
  173.          K = (X(I)-XMIN)/D+2.0
  174.          IF (K.LT.1) K = 1
  175.          IF (K.GT.KK) K = KK
  176.          N(K) = N(K) + 1
  177.    15 CONTINUE
  178.       NMAX = 0
  179.       DO 20 K=1,KK
  180.          IF (N(K).GT.NMAX) NMAX = N(K)
  181.    20 CONTINUE
  182.       KEISU = (NMAX-1)/50+1
  183.       J5 = 5*KEISU
  184.       J50 = 10*J5
  185.       WRITE(2,210) (J,J=J5,J50,J5)
  186.   210   FORMAT(1H ,21X,1H|,10I5/4X,'From',6X,'To',6X,51('-'))
  187.       XH = XMIN
  188.       NT = 0
  189.       DO 70 K=1,KK
  190.          DO 30 J=1,51
  191.             H(J) = HSPACE
  192.    30    CONTINUE         
  193.          DO 35 J=1,51,10
  194.             H(J) = HVRTL
  195.    35    CONTINUE
  196.          NK = (N(K)+KIESU-1)/KEISU+2
  197.          IF (NK.EQ.1)         GOTO 45
  198.          DO 40 J=2,NK
  199.             H(J) = HSTAR
  200.    40    CONTINUE
  201.    45    NT = NT + N(K)
  202.          J = FLOAT(NT*50)/FLOAT(II)+1.999
  203.          H(J) = HPLUS
  204.          IF (K.EQ.1)          GOTO 50
  205.          IF (K.EQ.KK)         GOTO 55
  206.          WRITE(2,220) XL,XH,N(K),(H(I),I=1,51)
  207.   220      FORMAT(1H ,F7.3,' -',F7.3,I4,1X,51A1)
  208.                               GOTO 60
  209.    50    IF (N(1).NE.0) WRITE(2,230) XH,N(1),(H(I),I=1,51)
  210.   230                     FORMAT(1H ,7X,' -',F7.3,I4,1X,51A1)
  211.                               GOTO 60
  212.    55    IF (N(K).NE.0) WRITE(2,240) XL,N(K),(H(I),I=1,51)
  213.   240                     FORMAT(1H ,F7.3,' -',7X,I4,1X,51A1)
  214.    60    XL = XH
  215.          XH = XH + D
  216.    70 CONTINUE
  217.       WRITE(2,250) (J,J=10,100,10)
  218.   250   FORMAT(1H ,3X,'From',6X,'To',6X,51('-')/22X,1H|,10I5,1H%///)
  219.       RETURN
  220.       END
  221. C
  222.       SUBROUTINE SORT1 (X,N)
  223. C
  224. C     ** SORT IN ASCENDING ORDER **
  225. C
  226.       DIMENSION X(1)
  227. C
  228.       NM1 = N-1
  229.       DO 20 I=1,NM1
  230.         K = I
  231.         IP1 = I+1
  232.         DO 10 J=IP1,N
  233.           IF (X(J).LT.X(K)) K = J
  234.    10   CONTINUE
  235.         W = X(I)
  236.         X(I) = X(K)
  237.         X(K) = W
  238.    20 CONTINUE
  239.       RETURN
  240.       END
  241. C
  242.       SUBROUTINE   NRPLOT(X,II)
  243. C
  244. C     ** NORMAL PROBABILITY PLOT **
  245. C
  246. C      * Must be called after SORT1!
  247. C
  248.       INTEGER*1    H(37,73),HSPACE,HVRTL,HRIZN,HSTAR
  249.       DIMENSION    X(II),SCALE(7)
  250.       DATA         J1,K1,LL,HSPACE,HVRTL,HRIZN,HSTAR /6,12,7,
  251.      1             ' ','|','-',Z'EC'/
  252. C
  253.       WRITE(2,200)
  254.   200   FORMAT(1H0,5X,'* Normal Probability Plot *'/8X,'(Unit of ',
  255.      1         'ordinate = 1 sigma)'/)
  256.       FII = II
  257.       JJ = 6*J1 + 1
  258.       KK = 6*K1 + 1
  259.       FJ1 = J1
  260.       FK1 = K1
  261.       FJ2 = FLOAT(JJ)/2.0+1.0
  262.       FK2 = FLOAT(KK)/2.0+1.0
  263.       CALL STAT0(X,II,SX,AVGX,SXX,SSX,VX,SIGX,SKEW,FKURT,GEARY)
  264.       DO 30 J=1,JJ
  265.          DO 10 K=1,KK
  266.             H(J,K) = HSPACE
  267.    10    CONTINUE
  268.          DO 20 K=1,KK,K1
  269.             H(J,K) = HVRTL
  270.    20    CONTINUE
  271.    30 CONTINUE
  272.       DO 50 K=1,KK
  273.          DO 40 J=1,JJ,J1
  274.             H(J,K) = HRIZN
  275.    40    CONTINUE
  276.    50 CONTINUE
  277.       DO 60 I=1,II
  278.          Q = (FLOAT(I)-0.5)/FII
  279.          QQ = Q
  280.          IF (Q.GT.0.5) QQ = 1.0 - Q
  281.          Z = SQRT(ALOG(1.0/QQ**2))
  282.          P = Z - (0.27061*Z+2.30753)/((0.04481*Z+0.9929)*Z+1.0)
  283.          IF (Q.GT.0.5) P = -P
  284.          J = P*FJ1 + FJ2
  285.          IF (J.LT.1) J = 1
  286.          IF (J.GT.JJ) J = JJ
  287.          K = (X(I)-AVGX)/SIGX*FK1+FK2
  288.          IF (K.LT.1) K = 1
  289.          IF (K.GT.KK) K = KK
  290.          H(J,K) = HSTAR
  291.    60 CONTINUE
  292.       DO 70 L=1,LL
  293.          SCALE(L) = AVGX + SIGX*(L-4)
  294.    70 CONTINUE
  295.       WRITE(2,210) (SCALE(L),L=1,LL)
  296.   210   FORMAT(1H ,F7.2,6F12.2)
  297.       DO 80 J=1,JJ
  298.          WRITE(2,220) (H(J,K),K=1,KK)
  299.   220      FORMAT(1H ,4X,73A1)
  300.    80 CONTINUE
  301.       WRITE(2,210) (SCALE(L),L=1,LL)
  302.       WRITE(2,230)
  303.   230   FORMAT(1H0)
  304.       RETURN
  305.       END
  306.