home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol269 / prbsub.for < prev    next >
Encoding:
Text File  |  1986-05-22  |  6.4 KB  |  225 lines
  1. C [PRBSUB.FOR]
  2. C    ** Function Subprograms on Satistical Distributions **
  3. C
  4. C    Written by Yshio MONMA (JUG-CP/M N.43)
  5. C
  6.       FUNCTION     PF(Q,N1,N2)
  7. C
  8. C     ** Percentage Point of F-Distribution **
  9. C
  10. C     * REFERENCE, T.Haga & S.Hashimoto: ó─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐ú, P.94
  11. C
  12. C     * Parameters:
  13. C     Q           Upper Probability
  14. C     N1,N2       Degree of Freedom
  15. C
  16. C                                       * (6.37)
  17.        F(X,D1,D2,DIM2) = (X/D1)*(1.0+(X-DIM2)/(2.0*D2)
  18.      1   +(4.0*X*X-11.0*DIM2*X+DIM2*(7.0*D1-10.0))/(24.0*D2*D2)
  19.      2   +(2.0*X*X*X-10.0*DIM2*X*X+DIM2*(17.0*D1-26.0)*X
  20.      3   -DIM2*DIM2*(9.0*D1-6.0))/(48.0*D2*D2*D2))
  21. C
  22.       DF1 = N1
  23.       DF2 = N2
  24.       IF (N2.EQ.1)                      GO TO 10
  25.       IF (N2.EQ.2)                      GO TO 30
  26.       IF (N1.EQ.1)                      GO TO 20
  27.       IF (N1.GT.N2)                     GO TO 30
  28.                                         GO TO 40
  29. C                                       * 1/T**2 (6.36)
  30.    10 W = 1.0
  31.       IF (MOD(N1,2).EQ.1) W = 2.0/3.14159
  32.       IF (N1.LE.2)                      GO TO 14
  33.       II = (N1-1)/2
  34.       G = (DF1-1.0)/2.0
  35.       DO 12 I=1,II
  36.          W = W*G/(G-0.5)
  37.          G = G-1.0
  38.    12 CONTINUE
  39.    14 PF = (W/Q)**2/DF1
  40.       RETURN
  41. C                                       * T**2 (6.31)
  42.    20 PF = PT(Q/2.0,N2)**2
  43.       RETURN
  44. C                                       * 1/Chi-2 (6.33) & (6.37)
  45.    30 PF = 1.0/F(PCHI2(1.0-Q,N2),DF2,DF1,DF2-2.0)
  46.       RETURN
  47. C                                       * Chi-2 (6.37)
  48.    40 PF = F(PCHI2(Q,N1),DF1,DF2,DF1-2.0)
  49.       IF (N2.GT.10) RETURN
  50. C
  51. C    Mean of two approximation      * (6.38)
  52.       A1 = 2.0/(9.0*DF1)
  53.       AA1 = 1.0-A1
  54.       A2 = 2.0/(9.0*DF2)
  55.       AA2 = 1.0-A2
  56.       U = PNOR(Q)
  57.       PFF = ((AA1*AA2+U*SQRT(AA1**2*A2+AA2**2*A1-A1*A2*U**2))
  58.      1        /(AA2**2-A2*U**2))**3.0
  59.       PF = (PF+PFF)/2.0
  60. C
  61.       RETURN
  62.       END
  63.       FUNCTION     PT(Q,NDF)
  64. C
  65. C     ** Percentage Point of t-Distribution **
  66. C
  67. C     * REFERENCE, T.Haga & S.Hashimoto ─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐, P.91
  68. C
  69. C     * Parameters:
  70. C     Q       Upper Probability of t-Distribution
  71. C     NDF     Degree of Freedom
  72. C
  73. C
  74.       DF = NDF
  75.       U = PNOR(Q)
  76.       U2 = U*U
  77.       PT = U*(1.0+(U2+1.0)/(4.0*DF)
  78.      1    +((5.0*U2+16.0)*U2+3.0)/(96.0*DF*DF)
  79.      2    +(((3.0*U2+19.0)*U2+17.0)*U2-15.0)/(384.0*DF*DF*DF)
  80.      3    +((((79.0*U2+776.0)*U2+1482.0)*U2-1920.0)*U2-945.0)
  81.      4         /(92160.0*DF*DF*DF*DF)
  82.      5    +(((((27.0*U2+339.0)*U2+930.0)*U2-1782.0)*U2-765.0)*U2
  83.      6          +17955.0)/(368640.0*DF*DF*DF*DF*DF))
  84. C
  85.       RETURN
  86.       END
  87.       FUNCTION     PCHI2(Q,NDF)
  88. C
  89. C     ** PERCENTAGE POINT OF CHI-SQUARE DISTRIBUTION **
  90. C
  91. C     * REFERENCE, T.HAGA & S.HASHIMOTO: ─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐, P.88
  92. C
  93. C     * PARAMETERS:
  94. C     Q     Upper probability
  95. C     NDF   Degree of freedom
  96. C
  97. C
  98.       IF (NDF.EQ.1)                     GO TO 10
  99.       IF (NDF.EQ.2)                     GO TO 20
  100.                                         GO TO 30
  101. C                                       * N = 1 (6.23)
  102.    10 PCHI2 = PNOR(Q/2.0)**2
  103.       RETURN
  104. C                                       * N = 2 (6.23)
  105.    20 PCHI2 = -2.0*ALOG(Q)
  106.       RETURN
  107. C                                       * N > 2 (6.24)
  108.    30 U = PNOR(Q)
  109.       W = 2.0/(9.0*NDF)
  110.       PCHI2 = NDF*(1.0-W+U*SQRT(W))**3
  111.       IF (PCHI2.LT.0.0) PCHI2 = 0.0
  112.       IF (NDF.GT.10)                    RETURN
  113.       G = 1.0
  114.       IF (MOD(NDF,2).EQ.1) G = 1.772454
  115.       N2 = (NDF+1)/2-1
  116.       DO 40 I=1,N2
  117.          G = G*(NDF/2.0-I)
  118.    40 CONTINUE
  119.       F = (PCHI2/2.0)**(NDF/2.0-1.0)*EXP(-PCHI2/2.0)/2.0
  120.       PCHI2 = PCHI2+(QCHI2(PCHI2,NDF)-Q)/(F/G)
  121. C
  122.       RETURN
  123.       END
  124.     FUNCTION     QCHI2(CHI2,NDF)
  125. C
  126. C     ** Upper Probability of Chi-Square Distribution **
  127. C
  128. C     * REFERENCE, T.Haga & S.Hashimoto: ó─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐ú, P.87
  129. C
  130. C     * Parameters:
  131. C     CHI2    Chi-Square
  132. C     NDF     Degree of freedom
  133. C
  134. C
  135.       IF (CHI2.LE.0.0)                  GO TO 90
  136.       IF (NDF.GT.10)                    GO TO 50
  137.       QCHI2 = 0.0
  138.       IF (CHI2.GT.400.0)                RETURN
  139.       A = EXP(-CHI2/2.0)
  140.       IF (MOD(NDF,2).EQ.0)              GO TO 10
  141. C                                       * NDF is odd (6.20)
  142.       CHI = SQRT(CHI2)
  143.       QCHI2 = 2.0*QNOR(CHI)
  144.       A = 0.7978845*A/CHI
  145.       I1 = 1
  146.                                         GO TO 20
  147. C                                       * NDF is even (6.20)
  148.    10 QCHI2 = A
  149.       I1 = 2
  150.    20 IF (NDF.LE.2)                     RETURN
  151.       I2 = NDF-2
  152.       DO 30 I=I1,I2,2
  153.          A = A*CHI2/I
  154.          QCHI2 = QCHI2+A
  155.    30 CONTINUE
  156.       RETURN
  157. C                                       * N > 10 (6.22)
  158.    50 W = 2.0/(9.0*NDF)
  159.       U = ((CHI2/NDF)**(1.0/3.0)-1.0+W)/SQRT(W)
  160.       CHI2 = QNOR(U)
  161.       RETURN
  162. C
  163.    90 QCHI2 = 1.0
  164. C
  165.       RETURN
  166.       END
  167.       FUNCTION     QNOR(U)
  168. C
  169. C     ** Upper Probability of Normal Distribution **
  170. C
  171. C     * REFERENCE, T.Haga & S.Hashimoto: ó─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐ú, P.81
  172. C
  173. C     * Parameters:
  174. C       U = (X-MU)/SIGMA for the normal distribution function;
  175. C
  176. C    F(X) = (1/(SQRT(2*PI)*SIGMA))*EXP(-(X-MU)**2/(2*SIGMA**2))
  177. C
  178.       DATA         N/12/
  179. C
  180.       FJ = N
  181.       QNOR = 0.0
  182.       X = ABS(U)
  183.       X2 = X*X
  184.       IF (X2.GT.300.0)                  GO TO 40
  185.       F = EXP(-0.5*X2)*0.398942
  186.       IF (X.LT.2.4)                     GO TO 20
  187.       DO 10 I=1,N
  188.          QNOR = FJ/(X+QNOR)
  189.          FJ = FJ-1.0
  190.    10 CONTINUE
  191.       QNOR = F/(X+QNOR)
  192.                                         GO TO 40
  193. C
  194.    20 S = -1.0
  195.       DO 30 I=1,N
  196.          QNOR = FJ*X2/(2.0*FJ+1.0+S*QNOR)
  197.          S = -S
  198.          FJ = FJ-1.0
  199.    30 CONTINUE
  200.       QNOR = 0.5-F*X/(1.0-QNOR)
  201. C
  202.    40 IF (U.LT.0.0) QNOR = 1.0-QNOR
  203. C
  204.       RETURN
  205.       END
  206.       FUNCTION     PNOR(Q)
  207. C
  208. C     ** Percentage Point of Normal Distribution **
  209. C
  210. C     * REFERENCE, T.Haga & S.Hashimoto: ó─│╣▓╢▓╛╖ PROGRAM ╔ ╖┐ú, P.83
  211. C
  212. C     * PARAMETER:
  213. C     Q     Upper probability of normal distribution
  214. C
  215.       QQ = Q
  216.       IF (Q.GT.0.5) QQ = 1.0-Q
  217.       Z = SQRT(-2.0*ALOG(QQ))
  218.       PNOR = Z-(0.27061*Z+2.30753)/((0.04481*Z+0.99229)*Z+1.0)
  219.       D = QNOR(PNOR)-QQ
  220.       PNOR = PNOR+D/(0.398942*EXP(-0.5*PNOR*PNOR))
  221.       IF (Q.GT.0.5) PNOR = -PNOR
  222. C
  223.       RETURN
  224.       END
  225.