home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / statrand / cdtr.for < prev    next >
Text File  |  1985-12-31  |  6KB  |  251 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C     SUBROUTINE CDTR
  5. C
  6. C     PURPOSE
  7. C        COMPUTES P(X) = PROBABILITY THAT THE RANDOM VARIABLE U,
  8. C        DISTRIBUTED ACCORDING TO THE CHI-SQUARE DISTRIBUTION WITH G
  9. C        DEGREES OF FREEDOM, IS LESS THAN OR EQUAL TO X.  F(G,X), THE
  10. C        ORDINATE OF THE CHI-SQUARE DENSITY AT X, IS ALSO COMPUTED.
  11. C
  12. C     USAGE
  13. C        CALL CDTR(X,G,P,D,IER)
  14. C
  15. C     DESCRIPTION OF PARAMETERS
  16. C        X    - INPUT SCALAR FOR WHICH P(X) IS COMPUTED.
  17. C        G    - NUMBER OF DEGREES OF FREEDOM OF THE CHI-SQUARE
  18. C          DISTRIBUTION.  G IS A CONTINUOUS PARAMETER.
  19. C        P    - OUTPUT PROBABILITY.
  20. C        D    - OUTPUT DENSITY.
  21. C        IER - RESULTANT ERROR CODE WHERE
  22. C        IER= 0 --- NO ERROR
  23. C        IER=-1 --- AN INPUT PARAMETER IS INVALID.  X IS LESS
  24. C               THAN 0.0, OR G IS LESS THAN 0.5 OR GREATER
  25. C               THAN 2*10**(+5).  P AND D ARE SET TO -1.E75.
  26. C        IER=+1 --- INVALID OUTPUT.  P IS LESS THAN ZERO OR
  27. C               GREATER THAN ONE, OR SERIES FOR T1 (SEE
  28. C               MATHEMATICAL DESCRIPTION) HAS FAILED TO
  29. C               CONVERGE.  P IS SET TO 1.E75.
  30. C
  31. C     REMARKS
  32. C        SEE MATHEMATICAL DESCRIPTION.
  33. C
  34. C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  35. C        DLGAM
  36. C        NDTR
  37. C
  38. C     METHOD
  39. C        REFER TO R.E. BARGMANN AND S.P. GHOSH, STATISTICAL
  40. C        DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE,
  41. C        IBM RESEARCH REPORT RC-1094, 1963.
  42. C
  43. C     ..................................................................
  44. C
  45.       SUBROUTINE CDTR(X,G,P,D,IER)
  46.       DOUBLE PRECISION XX,DLXX,X2,DLX2,GG,G2,DLT3,THETA,THP1,
  47.      1GLG2,DD,T11,SER,CC,XI,FAC,TLOG,TERM,GTH,A2,A,B,C,DT2,DT3,THPI
  48. C
  49. C     TEST FOR VALID INPUT DATA
  50. C
  51.       IF(G-(.5-1.E-5)) 590,10,10
  52.    10 IF(G-2.E+5) 20,20,590
  53.    20 IF(X) 590,30,30
  54. C
  55. C     TEST FOR X NEAR 0.0
  56. C
  57.    30 IF(X-1.E-8) 40,40,80
  58.    40 P=0.0
  59.       IF(G-2.) 50,60,70
  60.    50 D=1.E38
  61.       GO TO 610
  62.    60 D=0.5
  63.       GO TO 610
  64.    70 D=0.0
  65.       GO TO 610
  66. C
  67. C     TEST FOR X GREATER THAN 1.E+6
  68. C
  69.    80 IF(X-1.E+6) 100,100,90
  70.    90 D=0.0
  71.       P=1.0
  72.       GO TO 610
  73. C
  74. C     SET PROGRAM PARAMETERS
  75. C
  76.   100 XX=DBLE(X)
  77.       DLXX=DLOG(XX)
  78.       X2=XX/2.D0
  79.       DLX2=DLOG(X2)
  80.       GG=DBLE(G)
  81.       G2=GG/2.D0
  82. C
  83. C     COMPUTE ORDINATE
  84. C
  85.       CALL DLGAM(G2,GLG2,IOK)
  86.       DD=(G2-1.D0)*DLXX-X2-G2*.6931471805599453 -GLG2
  87.       IF(DD-1.68D02) 110,110,120
  88.   110 IF(DD+1.68D02) 130,130,140
  89.   120 D=1.E38
  90.       GO TO 150
  91.   130 D=0.0
  92.       GO TO 150
  93.   140 DD=DEXP(DD)
  94.       D=SNGL(DD)
  95. C
  96. C     TEST FOR G GREATER THAN 1000.0
  97. C     TEST FOR X GREATER THAN 2000.0
  98. C
  99.   150 IF(G-1000.) 160,160,180
  100.   160 IF(X-2000.) 190,190,170
  101.   170 P=1.0
  102.       GO TO 610
  103.   180 A=DLOG(XX/GG)/3.D0
  104.       A=DEXP(A)
  105.       B=2.D0/(9.D0*GG)
  106.       C=(A-1.D0+B)/DSQRT(B)
  107.       SC=SNGL(C)
  108.       CALL NDTR(SC,P,DUMMY)
  109.       GO TO 490
  110. C
  111. C     COMPUTE THETA
  112. C
  113.   190 K= IDINT(G2)
  114.       THETA=G2-DFLOAT(K)
  115.       IF(THETA-1.D-8) 200,200,210
  116.   200 THETA=0.D0
  117.   210 THP1=THETA+1.D0
  118. C
  119. C     SELECT METHOD OF COMPUTING T1
  120. C
  121.       IF(THETA) 230,230,220
  122.   220 IF(XX-10.D0) 260,260,320
  123. C
  124. C     COMPUTE T1 FOR THETA EQUALS 0.0
  125. C
  126.   230 IF(X2-1.68D02) 250,240,240
  127.   240 T1=1.0
  128.       GO TO 400
  129.   250 T11=1.D0-DEXP(-X2)
  130.       T1=SNGL(T11)
  131.       GO TO 400
  132. C
  133. C     COMPUTE T1 FOR THETA GREATER THAN 0.0 AND
  134. C     X LESS THAN OR EQUAL TO 10.0
  135. C
  136.   260 SER=X2*(1.D0/THP1 -X2/(THP1+1.D0))
  137.       J=+1
  138.       CC=DFLOAT(J)
  139.       DO 270 IT1=3,30
  140.       XI=DFLOAT(IT1)
  141.       CALL DLGAM(XI,FAC,IOK)
  142.       TLOG= XI*DLX2-FAC-DLOG(XI+THETA)
  143.       TERM=DEXP(TLOG)
  144.       TERM=DSIGN(TERM,CC)
  145.       SER=SER+TERM
  146.       CC=-CC
  147.       IF(DABS(TERM)-1.D-9) 280,270,270
  148.   270 CONTINUE
  149.       GO TO 600
  150.   280 IF(SER) 600,600,290
  151.   290 CALL DLGAM(THP1,GTH,IOK)
  152.       TLOG=THETA*DLX2+DLOG(SER)-GTH
  153.       IF(TLOG+1.68D02) 300,300,310
  154.   300 T1=0.0
  155.       GO TO 400
  156.   310 T11=DEXP(TLOG)
  157.       T1=SNGL(T11)
  158.       GO TO 400
  159. C
  160. C     COMPUTE T1 FOR THETA GREATER THAN 0.0 AND
  161. C     X GREATER THAN 10.0 AND LESS THAN 2000.0
  162. C
  163.   320 A2=0.D0
  164.       DO 340 I=1,25
  165.       XI=DFLOAT(I)
  166.       CALL DLGAM(THP1,GTH,IOK)
  167.       T11=-(13.D0*XX)/XI +THP1*DLOG(13.D0*XX/XI) -GTH-DLOG(XI)
  168.       IF(T11+1.68D02) 340,340,330
  169.   330 T11=DEXP(T11)
  170.       A2=A2+T11
  171.   340 CONTINUE
  172.       A=1.01282051+THETA/156.D0-XX/312.D0
  173.       B=DABS(A)
  174.       C= -X2+THP1*DLX2+DLOG(B)-GTH-3.951243718581427
  175.       IF(C+1.68D02) 370,370,350
  176.   350 IF (A) 360,370,380
  177.   360 C=-DEXP(C)
  178.       GO TO 390
  179.   370 C=0.D0
  180.       GO TO 390
  181.   380 C=DEXP(C)
  182.   390 C=A2+C
  183.       T11=1.D0-C
  184.       T1=SNGL(T11)
  185. C
  186. C     SELECT PROPER EXPRESSION FOR P
  187. C
  188.   400 IF(G-2.) 420,410,410
  189.   410 IF(G-4.) 450,460,460
  190. C
  191. C     COMPUTE P FOR G GREATER THAN ZERO AND LESS THAN 2.0
  192. C
  193.   420 CALL DLGAM(THP1,GTH,IOK)
  194.       DT2=THETA*DLXX-X2-THP1*.6931471805599453 -GTH
  195.       IF(DT2+1.68D02) 430,430,440
  196.   430 P=T1
  197.       GO TO 490
  198.   440 DT2=DEXP(DT2)
  199.       T2=SNGL(DT2)
  200.       P=T1+T2+T2
  201.       GO TO 490
  202. C
  203. C     COMPUTE P FOR G GREATER THAN OR EQUAL TO 2.0
  204. C     AND LESS THAN 4.0
  205. C
  206.   450 P=T1
  207.       GO TO 490
  208. C
  209. C     COMPUTE P FOR G GREATER THAN OR EQUAL TO 4.0
  210. C     AND LESS THAN OR EQUAL TO 1000.0
  211. C
  212.   460 DT3=0.D0
  213.       DO 480 I3=2,K
  214.       THPI=DFLOAT(I3)+THETA
  215.       CALL DLGAM(THPI,GTH,IOK)
  216.       DLT3=THPI*DLX2-DLXX-X2-GTH
  217.       IF(DLT3+1.68D02) 480,480,470
  218.   470 DT3=DT3+DEXP(DLT3)
  219.   480 CONTINUE
  220.       T3=SNGL(DT3)
  221.       P=T1-T3-T3
  222. C
  223. C     SET ERROR INDICATOR
  224. C
  225.   490 IF(P) 500,520,520
  226.   500 IF(ABS(P)-1.E-7) 510,510,600
  227.   510 P=0.0
  228.       GO TO 610
  229.   520 IF(1.-P) 530,550,550
  230.   530 IF(ABS(1.-P)-1.E-7) 540,540,600
  231.   540 P=1.0
  232.       GO TO 610
  233.   550 IF(P-1.E-8) 560,560,570
  234.   560 P=0.0
  235.       GO TO 610
  236.   570 IF((1.0-P)-1.E-8) 580,580,610
  237.   580 P=1.0
  238.       GO TO 610
  239.   590 IER=-1
  240.       D=-1.E38
  241.       P=-1.E38
  242.       GO TO 620
  243.   600 IER=+1
  244.       P= 1.E38
  245.       GO TO 620
  246.   610 IER=0
  247.   620 RETURN
  248.       END
  249. 1.E-7) 510,510,600
  250.   510 P=0.0
  251.