DCDFLIB Library of Fortran Routines for Cumulative Distribution Functions, Inverses, and Other Parameters (February, 1994) Full Documentation of Each Routine Compiled and Written by: Barry W. Brown James Lovato Kathy Russell Department of Biomathematics, Box 237 The University of Texas, M.D. Anderson Cancer Center 1515 Holcombe Boulevard Houston, TX 77030 This work was supported by grant CA-16672 from the National Cancer Institute. Summary Documentation of Each Routine WHICH and STATUS are integer; all other arguments are Double Precision. -------------------- DISTRIBUTION WHICH PARAMETERS INPUT RANGE SEARCH RANGE REQUIREMENTS -------------------- Beta SUBROUTINE CDFBET( WHICH, P, Q, X, Y, A, B, STATUS, BOUND ) 1 P and Q [0,1];[0,1] ----------- SUM to 1.0 2 X and Y [0,1];[0,1] [0,1],[0,1] SUM to 1.0 3 A (0,infinity) [1E-300,1E300] 4 B (0,infinity) [1E-300,1E300] -------------------- Binomial SUBROUTINE CDFBIN ( WHICH, P, Q, S, XN, PR, OMPR, STATUS, BOUND ) 1 P and Q [0,1];[0,1] ----------- SUM to 1.0 2 S [0,XN] [0,XN] 3 XN (0,infinity) [1E-300,1E300] 4 PR and OMPR [0,1];[0,1] [0,1];[0,1] SUM to 1.0 -------------------- Chi-square SUBROUTINE CDFCHI( WHICH, P, Q, X, DF, STATUS, BOUND ) 1 P and Q [0,1],(0,1] ----------- SUM to 1.0 2 X [0,infinity] [0,1E300] 3 DF (0,infinity) [1E-300,1E300] -------------------- Noncentral Chi-square SUBROUTINE CDFCHN( WHICH, P, Q, X, DF, PNONC, STATUS, BOUND ) 1 P and Q [0,1-1E-16],none ----------- 2 X [0,infinity] [0,1E300] 3 DF (0,infinity) [1E-300,1E300] 4 PNONC [0,infinity) [0,1E4] NOTE: We do not yet have a method to calculation the Noncentral Chi-Square distribution acurately near 1; therefore, Q is not used by CDFCHN. There are no input requirements of Q, and when WHICH is 1, Q is returned as 1-P. -------------------- F SUBROUTINE CDFF( WHICH, P, Q, F, DFN, DFD, STATUS, BOUND ) 1 P and Q [0,1],(0,1] ----------- SUM to 1.0 2 F [0,infinity) [0,1E300] 3 DFN (0,infinity) [1E-300,1E300] 4 DFD (0,infinity) [1E-300,1E300] -------------------- Noncentral F SUBROUTINE CDFFNC( WHICH, P, Q, F, DFN, DFD, PNONC, STATUS, BOUND ) 1 P and Q [0,1-1E-16],none ----------- 2 F [0,infinity) [0,1E300] 3 DFN (0,infinity) [1E-300,1E300] 4 DFD (0,infinity) [1E-300,1E300] 5 PNONC [0,infinity) [0,1E4] NOTE: We do not yet have a method to calculation the Noncentral F distribution acurately near 1; therefore, Q is not used by CDFF. There are no input requirements of Q, and when WHICH is 1, Q is returned as 1-P. -------------------- Gamma SUBROUTINE CDFGAM( WHICH, P, Q, X, SHAPE, SCALE, STATUS, BOUND ) 1 P and Q [0,1],(0,1] ----------- SUM to 1.0 2 X [0,infinity) [0,1E300] 3 SHAPE (0,infinity) [1E-300,1E300] 4 SCALE (0,infinity) [1E-300,1E300] -------------------- Negative Binomial SUBROUTINE CDFNBN ( WHICH, P, Q, S, XN, PR, OMPR, STATUS, BOUND ) 1 P and Q [0,1];(0,1] ----------- SUM to 1.0 2 S [0,infinity) [0,1E300] 3 XN [0,infinity) [0,1E300] 4 PR and OMPR [0,1];[0,1] [0,1];[0,1] SUM to 1.0 -------------------- Normal SUBROUTINE CDFNOR( WHICH, P, Q, X, MEAN, SD, STATUS, BOUND ) 1 P and Q (0,1],(0,1] ----------- SUM to 1.0 2 X (-inf.,inf.) ----------- 3 MEAN (-inf.,inf.) ----------- 4 SD (0,infinity) ----------- -------------------- Poisson SUBROUTINE CDFPOI( WHICH, P, Q, S, XLAM, STATUS, BOUND ) 1 P and Q [0,1],(0,1] ----------- SUM to 1.0 2 S [0,infinity) [0,1E300] 3 XLAM [0,infinity) [0,1E300] -------------------- Student's t SUBROUTINE CDFT( WHICH, P, Q, T, DF, STATUS, BOUND ) 1 P and Q (0,1],(0,1] ----------- SUM to 1.0 2 T (-inf.,inf.) [-1E300,1E300] 3 DF (0,infinity) [1E-300,1E10] SUMMARY OF DCDFLIB This library contains routines to compute cumulative distribution functions, inverses, and parameters of the distribution for the following set of statistical distributions: (1) Beta (2) Binomial (3) Chi-square (4) Noncentral Chi-square (5) F (6) Noncentral F (7) Gamma (8) Negative Binomial (9) Normal (10) Poisson (11) Student's t Given values of all but one parameter of a distribution, the other is computed. These calculations are done with FORTRAN Double Precision variables. -------------------- WARNINGS -------------------- The F and Noncentral F distribution are not necessarily monotone in either degree of freedom argument. Consequently, there may be two degree of freedom arguments that satisfy the specified condition. An arbitrary one of these will be found by the cdf routines. The amount of computation required for the noncentral chisquare and noncentral F distribution is proportional to the value of the noncentrality parameter. Very large values of this parameter can require immense numbers of computation. Consequently, when the noncentrality parameter is to be calculated, the upper limit searched is 10,000. -------------------- END WARNINGS -------------------- DOCUMENTATION This file contains an overview of the library and is the primary documentation. Other documentation is in directory 'doc' on the distribution as character (ASCII) files. A summary of all of the available routines is contained in dcdflib.chs (chs is an abbreviation of 'cheat sheet'). The 'chs' file will probably be the primary reference. The file, dcdflib.fdoc, contains the header comments for each routine intended for direct use. INSTALLATION The Fortran source routines are contained in directory src. A few routines use machine dependent constants. Lists of such constants for different machines are found in ipmpar.f. Uncomment the ones appropriate to your machine. The distributed version uses the IEEE arithmetic that is used by the IBM PC, Macintosh, and most Unix workstations. Ignore compilation warnings that lines of code are not reachable. We write in a Fortran structured preprocessor (FLECS) that is similar in spirit to Ratfor. Sometimes our coding practices in FLECS lead to unreachable lines. Also, FLECS inserts lines of code: STOP "CODE FLOWING INTO FLECS PROCEEDURES". All such lines should be unreachable. SOURCES The following routines, written by others, are incorporated into DCDFLIB. Beta Distribution DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant Digit Computation of the Incomplete Beta Function Ratios. ACM Trans. Math. Softw. 18 (1993), 360-373. Gamma Distribution and It's Inverse DiDinato, A. R. and Morris, A. H. Computation of the Incomplete Gamma Function Ratios and their Inverse. ACM Trans. Math. Softw. 12 (1986), 377-393. Normal Distribution Kennedy and Gentle, Statistical Computing, Marcel Dekker, NY, 1980. The rational function approximations from pages 90-95 are used during the calculation of the inverse normal. Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN Package of Special Function Routines and Test Drivers", acm Transactions on Mathematical Software. 19, 22-32. A slightly modified version of Cody's function anorm is used for the cumultive normal. Zero Finder J. C. P. Bus and T. J. Dekker. Two Efficient Algorithms with Guaranteed Convergence for Finding a Zero of a Function. ACM Trans. Math. Softw. 4 (1975), 330. We transliterated Algoritm R of this paper from Algol to Fortran. General Reference Abramowitz, M. and Stegun, I. A. Handbook of Mathematical Functions With Formulas, Graphs, and Mathematical Tables. (1964) National Bureau of Standards. This book has been reprinted by Dover and others. LEGALITIES Code that appeared in an ACM publication is subject to their algorithms policy: Submittal of an algorithm for publication in one of the ACM Transactions implies that unrestricted use of the algorithm within a computer is permissible. General permission to copy and distribute the algorithm without fee is granted provided that the copies are not made or distributed for direct commercial advantage. The ACM copyright notice and the title of the publication and its date appear, and notice is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission. Krogh, F. Algorithms Policy. ACM Tran. Math. Softw. 13(1987), 183-186. We place the DCDFLIB code that we have written in the public domain. NO WARRANTY WE PROVIDE ABSOLUTELY NO WARRANTY OF ANY KIND EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THIS PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. IN NO EVENT SHALL THE UNIVERSITY OF TEXAS OR ANY OF ITS COMPONENT INSTITUTIONS INCLUDING M. D. ANDERSON HOSPITAL BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY LOST PROFITS, LOST MONIES, OR OTHER SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA OR ITS ANALYSIS BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY THIRD PARTIES) THE PROGRAM. (Above NO WARRANTY modified from the GNU NO WARRANTY statement.) HOW TO USE THE ROUTINES The calling sequence for each routine is of the form: SUBROUTINE CDF( WHICH, P, Q, X, , STATUS, BOUND ). WHICH and STATUS are INTEGER, all other arguments are DOUBLE PRECISION. is a one to three character name identifying the distribution. WHICH is an input integer value that identifies what parameter value is to be calculated from the values of the other parameters. P is always the cdf evaluated at X, Q is always the compliment of the cdf evaluated at X, i.e. 1-P, and X is always the value at which the cdf is evaluated. The auxiliary parameters, , of the distribution differ by distribution. If WHICH is 1, P and Q are to be calculated, i.e., the cdf; if WHICH is 2, X is to be calculated, i.e., the inverse cdf. The value of one auxiliary parameter in can also be the value calculated. STATUS returns 0 if the calculation completes correctly. --------------------WARNING-------------------- If STATUS is not 0, no meaningful answer is returned. -------------------- END WARNING -------------------- STATUS returns -I if the I'th input parameter was not in the legal range (see below). Parameters are counted with WHICH being the first in these return values. A STATUS value of 1 indicates that the desired answer was apparently lower than the lower bound on the search interval. A return code of 2 indicates that the answer was apparently higher than the upper bound on the search interval. A return code of 3 indicates that P and Q did not sum to 1. Other positive codes are routine specific. BOUND is not set if STATUS is returned as 0. If STATUS is -I then BOUND is the bound illegally exceeded by input parameter I, where WHICH is counted as 1, P as 2, Q as 3, X as 4, etc. If STATUS is returned as 1 or 2 then BOUND is returned as the lower or upper bound on the search interval respectively. BOUNDS Below are the rules that we used in determining bounds on quantities to be calculated. Those who don't care can find a summary of the bounds in dcdflib.chs. Input bounds are checked for legality of input. The search range is the range of values searched for an answer. Input Bounds Bounds on input parameters are checked by the CDF* routines. These bounds were set according to the following rules. P: If the domain of the cdf (X) extends to -infinity then P must be greater than 0 otherwise P must be greater than or equal to 0. P must always be less than or equal to 1. Q: If the domain of the cdf (X) extends to +infinity then Q must be greater than 0 otherwise Q must be greater than or equal to 0. Q must always be less than or equal to 1. Further, P and Q must sum to 1. The smaller of the two P and Q will be used in calculations to increase accuracy X: If the domain is infinite in either the positive or negative direction, no check is performed in that direction. If the left end of the domain is 0, then X is checked to assure non-negativity. DF, SD, etc.: Some auxiliary parameters must be positive. The lowest input values accepted for these parameters is 1E-300. Search Bounds These are the ranges searched for an answer. If the domain of the parameter in the cdf is closed at some finite value, e.g., 0, then this value is the same endpoint of the search range. If the domain is open at some finite endpoint (which only occurs for 0 -- some parameters must be strictly positive) then the endpoint is 1E-300. If the domain is infinite in either direction then +/- 1E300 is used as the endpoint of the search range. HOW THE ROUTINES WORK The cumulative distribution functions are computed directly. The normal, gamma, and beta functions use the code from the references cited. Other cdfs are calculated by relating them to one of these distributions. For example, the binomial and negative binomial cdfs can be converted to a beta cdf. This is how fractional observations are handled. The formula from Abramowitz and Stegun for converting the cdfs is cited in the fdoc file. (We think the formula for the negative binomial in A&S is wrong, but there is a correct one which we used.) The inverse normal and gamma are also taken from the references. For all other parameters, a search is made for the value that provides the desired P. Initial values are chosen crudely for the search (e.g., 5). If the domain of the cdf for the parameter being calculated is infinite, a step doubling strategy is used to bound the desired value then the zero finder is employed to refine the answer. The zero finder attempts to obtain the answer accurately to about eight decimal places. C********************************************************************** C C SUBROUTINE CDFBET( WHICH, P, Q, X, Y, A, B, STATUS, BOUND ) C Cumulative Distribution Function C BETa Distribution C C C Function C C C Calculates any one parameter of the beta distribution given C values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next four argument C values is to be calculated from the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from X,Y,A and B C iwhich = 2 : Calculate X and Y from P,Q,A and B C iwhich = 3 : Calculate A from P,Q,X,Y and B C iwhich = 4 : Calculate B from P,Q,X,Y and A C C INTEGER WHICH C C P <--> The integral from 0 to X of the chi-square C distribution. C Input range: [0, 1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: [0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C X <--> Upper limit of integration of beta density. C Input range: [0,1]. C Search range: [0,1] C DOUBLE PRECISION X C C Y <--> 1-X. C Input range: [0,1]. C Search range: [0,1] C X + Y = 1.0. C DOUBLE PRECISION Y C C A <--> The first parameter of the beta density. C Input range: (0, +infinity). C Search range: [1D-300,1D300] C DOUBLE PRECISION A C C B <--> The second parameter of the beta density. C Input range: (0, +infinity). C Search range: [1D-300,1D300] C DOUBLE PRECISION B C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C 4 if X + Y .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Cumulative distribution function (P) is calculated directly by C code associated with the following reference. C C DiDinato, A. R. and Morris, A. H. Algorithm 708: Significant C Digit Computation of the Incomplete Beta Function Ratios. ACM C Trans. Math. Softw. 18 (1993), 360-373. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C Note C C C The beta density is proportional to C t^(A-1) * (1-t)^(B-1) C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFBIN ( WHICH, P, Q, S, XN, PR, OMPR, STATUS, BOUND ) C Cumulative Distribution Function C BINomial distribution C C C Function C C C Calculates any one parameter of the binomial C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next four argument C values is to be calculated from the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR C iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR C iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR C iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN C INTEGER WHICH C C P <--> The cumulation from 0 to S of the binomial distribution. C (Probablility of S or fewer successes in XN trials each C with probability of success PR.) C Input range: [0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: [0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C S <--> The number of successes observed. C Input range: [0, XN] C Search range: [0, XN] C DOUBLE PRECISION S C C XN <--> The number of binomial trials. C Input range: (0, +infinity). C Search range: [1E-300, 1E300] C DOUBLE PRECISION XN C C PR <--> The probability of success in each binomial trial. C Input range: [0,1]. C Search range: [0,1] C DOUBLE PRECISION PR C C OMPR <--> 1-PR C Input range: [0,1]. C Search range: [0,1] C PR + OMPR = 1.0 C DOUBLE PRECISION OMPR C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C 4 if PR + OMPR .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.5.24 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce the binomial C distribution to the cumulative incomplete beta distribution. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFCHI( WHICH, P, Q, X, DF, STATUS, BOUND ) C Cumulative Distribution Function C CHI-Square distribution C C C Function C C C Calculates any one parameter of the chi-square C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next three argument C values is to be calculated from the others. C Legal range: 1..3 C iwhich = 1 : Calculate P and Q from X and DF C iwhich = 2 : Calculate X from P,Q and DF C iwhich = 3 : Calculate DF from P,Q and X C INTEGER WHICH C C P <--> The integral from 0 to X of the chi-square C distribution. C Input range: [0, 1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C X <--> Upper limit of integration of the non-central C chi-square distribution. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION X C C DF <--> Degrees of freedom of the C chi-square distribution. C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DF C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C 10 indicates error returned from cumgam. See C references in cdfgam C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.4.19 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce the chisqure C distribution to the incomplete distribution. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFCHN( WHICH, P, Q, X, DF, PNONC, STATUS, BOUND ) C Cumulative Distribution Function C Non-central Chi-Square C C C Function C C C Calculates any one parameter of the non-central chi-square C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next three argument C values is to be calculated from the others. C Input range: 1..4 C iwhich = 1 : Calculate P and Q from X and DF C iwhich = 2 : Calculate X from P,DF and PNONC C iwhich = 3 : Calculate DF from P,X and PNONC C iwhich = 3 : Calculate PNONC from P,X and DF C INTEGER WHICH C C P <--> The integral from 0 to X of the non-central chi-square C distribution. C Input range: [0, 1-1E-16). C DOUBLE PRECISION P C C Q <--> 1-P. C Q is not used by this subroutine and is only included C for similarity with other cdf* routines. C DOUBLE PRECISION Q C C X <--> Upper limit of integration of the non-central C chi-square distribution. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION X C C DF <--> Degrees of freedom of the non-central C chi-square distribution. C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DF C C PNONC <--> Non-centrality parameter of the non-central C chi-square distribution. C Input range: [0, +infinity). C Search range: [0,1E4] C DOUBLE PRECISION PNONC C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.4.25 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to compute the cumulative C distribution function. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C WARNING C C The computation time required for this routine is proportional C to the noncentrality parameter (PNONC). Very large values of C this parameter can consume immense computer resources. This is C why the search range is bounded by 10,000. C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFF( WHICH, P, Q, F, DFN, DFD, STATUS, BOUND ) C Cumulative Distribution Function C F distribution C C C Function C C C Calculates any one parameter of the F distribution C given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next four argument C values is to be calculated from the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from F,DFN and DFD C iwhich = 2 : Calculate F from P,Q,DFN and DFD C iwhich = 3 : Calculate DFN from P,Q,F and DFD C iwhich = 4 : Calculate DFD from P,Q,F and DFN C INTEGER WHICH C C P <--> The integral from 0 to F of the f-density. C Input range: [0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C F <--> Upper limit of integration of the f-density. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION F C C DFN < --> Degrees of freedom of the numerator sum of squares. C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DFN C C DFD < --> Degrees of freedom of the denominator sum of squares. C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DFD C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.6.2 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce the computation C of the cumulative distribution function for the F variate to C that of an incomplete beta. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C WARNING C C The value of the cumulative F distribution is not necessarily C monotone in either degrees of freedom. There thus may be two C values that provide a given CDF value. This routine assumes C monotonicity and will find an arbitrary one of the two values. C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFFNC( WHICH, P, Q, F, DFN, DFD, PNONC, STATUS, BOUND ) C Cumulative Distribution Function C Non-central F distribution C C C Function C C C Calculates any one parameter of the Non-central F C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next five argument C values is to be calculated from the others. C Legal range: 1..5 C iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC C iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC C iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC C iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC C iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD C INTEGER WHICH C C P <--> The integral from 0 to F of the non-central f-density. C Input range: [0,1-1E-16). C DOUBLE PRECISION P C C Q <--> 1-P. C Q is not used by this subroutine and is only included C for similarity with other cdf* routines. C DOUBLE PRECISION Q C C F <--> Upper limit of integration of the non-central f-density. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION F C C DFN < --> Degrees of freedom of the numerator sum of squares. C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DFN C C DFD < --> Degrees of freedom of the denominator sum of squares. C Must be in range: (0, +infinity). C Input range: (0, +infinity). C Search range: [ 1E-300, 1E300] C DOUBLE PRECISION DFD C C PNONC <-> The non-centrality parameter C Input range: [0,infinity) C Search range: [0,1E4] C DOUBLE PRECISION PHONC C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.6.20 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to compute the cumulative C distribution function. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C WARNING C C The computation time required for this routine is proportional C to the noncentrality parameter (PNONC). Very large values of C this parameter can consume immense computer resources. This is C why the search range is bounded by 10,000. C C WARNING C C The value of the cumulative noncentral F distribution is not C necessarily monotone in either degrees of freedom. There thus C may be two values that provide a given CDF value. This routine C assumes monotonicity and will find an arbitrary one of the two C values. C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFGAM( WHICH, P, Q, X, SHAPE, SCALE, STATUS, BOUND ) C Cumulative Distribution Function C GAMma Distribution C C C Function C C C Calculates any one parameter of the gamma C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next four argument C values is to be calculated from the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from X,SHAPE and SCALE C iwhich = 2 : Calculate X from P,Q,SHAPE and SCALE C iwhich = 3 : Calculate SHAPE from P,Q,X and SCALE C iwhich = 4 : Calculate SCALE from P,Q,X and SHAPE C INTEGER WHICH C C P <--> The integral from 0 to X of the gamma density. C Input range: [0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C C X <--> The upper limit of integration of the gamma density. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION X C C SHAPE <--> The shape parameter of the gamma density. C Input range: (0, +infinity). C Search range: [1E-300,1E300] C DOUBLE PRECISION SHAPE C C C SCALE <--> The scale parameter of the gamma density. C Input range: (0, +infinity). C Search range: (1E-300,1E300] C DOUBLE PRECISION SCALE C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C 10 if the gamma or inverse gamma routine cannot C compute the answer. Usually happens only for C X and SHAPE very large (gt 1E10 or more) C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Cumulative distribution function (P) is calculated directly by C the code associated with: C C DiDinato, A. R. and Morris, A. H. Computation of the incomplete C gamma function ratios and their inverse. ACM Trans. Math. C Softw. 12 (1986), 377-393. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C Note C C C C The gamma density is proportional to C T**(SHAPE - 1) * EXP(- SCALE * T) C C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFNBN ( WHICH, P, S, XN, PR, STATUS, BOUND ) C Cumulative Distribution Function C Negative BiNomial distribution C C C Function C C C Calculates any one parameter of the negative binomial C distribution given values for the others. C C The cumulative negative binomial distribution returns the C probability that there will be F or fewer failures before the C XNth success in binomial trials each of which has probability of C success PR. C C The individual term of the negative binomial is the probability of C S failures before XN successes and is C Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S C C C Arguments C C C WHICH --> Integer indicating which of the next four argument C values is to be calculated from the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR C iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR C iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR C iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN C INTEGER WHICH C C P <--> The cumulation from 0 to S of the negative C binomial distribution. C Input range: [0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C S <--> The upper limit of cumulation of the binomial distribution. C There are F or fewer failures before the XNth success. C Input range: [0, +infinity). C Search range: [0, 1E300] C DOUBLE PRECISION S C C XN <--> The number of successes. C Input range: [0, +infinity). C Search range: [0, 1E300] C DOUBLE PRECISION XN C C PR <--> The probability of success in each binomial trial. C Input range: [0,1]. C Search range: [0,1]. C DOUBLE PRECISION PR C C OMPR <--> 1-PR C Input range: [0,1]. C Search range: [0,1] C PR + OMPR = 1.0 C DOUBLE PRECISION OMPR C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C 4 if PR + OMPR .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.5.26 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce calculation of C the cumulative distribution function to that of an incomplete C beta. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFNOR( WHICH, P, Q, X, MEAN, SD, STATUS, BOUND ) C Cumulative Distribution Function C NORmal distribution C C C Function C C C Calculates any one parameter of the normal C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which of the next parameter C values is to be calculated using values of the others. C Legal range: 1..4 C iwhich = 1 : Calculate P and Q from X,MEAN and SD C iwhich = 2 : Calculate X from P,Q,MEAN and SD C iwhich = 3 : Calculate MEAN from P,Q,X and SD C iwhich = 4 : Calculate SD from P,Q,X and MEAN C INTEGER WHICH C C P <--> The integral from -infinity to X of the normal density. C Input range: (0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C X < --> Upper limit of integration of the normal-density. C Input range: ( -infinity, +infinity) C DOUBLE PRECISION X C C MEAN <--> The mean of the normal density. C Input range: (-infinity, +infinity) C DOUBLE PRECISION MEAN C C SD <--> Standard Deviation of the normal density. C Input range: (0, +infinity). C DOUBLE PRECISION SD C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C C C A slightly modified version of ANORM from C C Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN C Package of Special Function Routines and Test Drivers" C acm Transactions on Mathematical Software. 19, 22-32. C C is used to calulate the cumulative standard normal distribution. C C The rational functions from pages 90-95 of Kennedy and Gentle, C Statistical Computing, Marcel Dekker, NY, 1980 are used as C starting values to Newton's Iterations which compute the inverse C standard normal. Therefore no searches are necessary for any C parameter. C C For X < -15, the asymptotic expansion for the normal is used as C the starting value in finding the inverse standard normal. C This is formula 26.2.12 of Abramowitz and Stegun. C C C Note C C C The normal density is proportional to C exp( - 0.5 * (( X - MEAN)/SD)**2) C C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFPOI( WHICH, P, Q, S, XLAM, STATUS, BOUND ) C Cumulative Distribution Function C POIsson distribution C C C Function C C C Calculates any one parameter of the Poisson C distribution given values for the others. C C C Arguments C C C WHICH --> Integer indicating which argument C value is to be calculated from the others. C Legal range: 1..3 C iwhich = 1 : Calculate P and Q from S and XLAM C iwhich = 2 : Calculate A from P,Q and XLAM C iwhich = 3 : Calculate XLAM from P,Q and S C INTEGER WHICH C C P <--> The cumulation from 0 to S of the poisson density. C Input range: [0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C S <--> Upper limit of cumulation of the Poisson. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION S C C XLAM <--> Mean of the Poisson distribution. C Input range: [0, +infinity). C Search range: [0,1E300] C DOUBLE PRECISION XLAM C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.4.21 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce the computation C of the cumulative distribution function to that of computing a C chi-square, hence an incomplete gamma function. C C Cumulative distribution function (P) is calculated directly. C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C C********************************************************************** C********************************************************************** C C SUBROUTINE CDFT( WHICH, P, Q, T, DF, STATUS, BOUND ) C Cumulative Distribution Function C T distribution C C C Function C C C Calculates any one parameter of the t distribution given C values for the others. C C C Arguments C C C WHICH --> Integer indicating which argument C values is to be calculated from the others. C Legal range: 1..3 C iwhich = 1 : Calculate P and Q from T and DF C iwhich = 2 : Calculate T from P,Q and DF C iwhich = 3 : Calculate DF from P,Q and T C INTEGER WHICH C C P <--> The integral from -infinity to t of the t-density. C Input range: (0,1]. C DOUBLE PRECISION P C C Q <--> 1-P. C Input range: (0, 1]. C P + Q = 1.0. C DOUBLE PRECISION Q C C T <--> Upper limit of integration of the t-density. C Input range: ( -infinity, +infinity). C Search range: [ -1E300, 1E300 ] C DOUBLE PRECISION T C C DF <--> Degrees of freedom of the t-distribution. C Input range: (0 , +infinity). C Search range: [1e-300, 1E10] C DOUBLE PRECISION DF C C STATUS <-- 0 if calculation completed correctly C -I if input parameter number I is out of range C 1 if answer appears to be lower than lowest C search bound C 2 if answer appears to be higher than greatest C search bound C 3 if P + Q .ne. 1 C INTEGER STATUS C C BOUND <-- Undefined if STATUS is 0 C C Bound exceeded by parameter number I if STATUS C is negative. C C Lower search bound if STATUS is 1. C C Upper search bound if STATUS is 2. C C C Method C C C Formula 26.5.27 of Abramowitz and Stegun, Handbook of C Mathematical Functions (1966) is used to reduce the computation C of the cumulative distribution function to that of an incomplete C beta. C C Computation of other parameters involve a seach for a value that C produces the desired value of P. The search relies on the C monotinicity of P with the other parameter. C C**********************************************************************