subroutine anglit_cdf ( x, cdf ) ! !******************************************************************************* ! !! ANGLIT_CDF evaluates the Anglit CDF. ! ! ! Modified: ! ! 29 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( x < - 0.25E+00 * pi ) then cdf = 0.0E+00 else if ( x < 0.25E+00 * pi ) then cdf = 0.5E+00 - 0.5E+00 * cos ( 2.0E+00 * x + pi / 2.0E+00 ) else cdf = 1.0E+00 end if return end subroutine anglit_cdf_inv ( cdf, x ) ! !******************************************************************************* ! !! ANGLIT_CDF_INV inverts the Anglit CDF. ! ! ! Modified: ! ! 29 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Output, real X, the corresponding argument. ! implicit none ! real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'ANGLIT_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if x = 0.5E+00 * ( acos ( 1.0E+00 - 2.0E+00 * cdf ) - pi / 2.0E+00 ) return end subroutine anglit_mean ( mean ) ! !******************************************************************************* ! !! ANGLIT_MEAN returns the mean of the Anglit PDF. ! ! ! Modified: ! ! 28 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real mean ! mean = 0.0E+00 return end subroutine anglit_pdf ( x, pdf ) ! !******************************************************************************* ! !! ANGLIT_PDF evaluates the Anglit PDF. ! ! ! Formula: ! ! PDF(X) = SIN ( 2 * X + PI / 2 ) for -PI/4 <= X <= PI/4 ! ! Modified: ! ! 28 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( x <= - 0.25E+00 * pi .or. x >= 0.25E+00 * pi ) then pdf = 0.0E+00 else pdf = sin ( 2.0E+00 * x + 0.25E+00 * pi ) end if return end subroutine anglit_sample ( x ) ! !******************************************************************************* ! !! ANGLIT_SAMPLE samples the Anglit PDF. ! ! ! Modified: ! ! 28 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real X, a sample of the PDF. ! implicit none ! real cdf real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call anglit_cdf_inv ( cdf, x ) return end subroutine anglit_variance ( variance ) ! !******************************************************************************* ! !! ANGLIT_VARIANCE returns the variance of the Anglit PDF. ! ! ! Discussion: ! ! Variance = ! Integral ( -PI/4 <= X <= PI/4 ) X**2 * SIN ( 2 * X + PI / 2 ) ! ! Antiderivative = ! 0.5E+00 * X * SIN ( 2 * X + PI / 2 ) ! + ( 0.25 - 0.5E+00 * X**2 ) * COS ( 2 * X + PI / 2 ) ! ! Modified: ! ! 29 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real variance ! variance = 0.0625E+00 * pi**2 - 0.5E+00 return end subroutine arcsin_cdf ( x, cdf ) ! !******************************************************************************* ! !! ARCSIN_CDF evaluates the Arcsin CDF. ! ! ! Modified: ! ! 24 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( x < 0.0E+00 ) then cdf = 0.0E+00 else if ( x < 1.0E+00 ) then cdf = 2.0E+00 * asin ( x ) / pi else cdf = 1.0E+00 end if return end subroutine arcsin_cdf_inv ( cdf, x ) ! !******************************************************************************* ! !! ARCSIN_CDF_INV inverts the Arcsin CDF. ! ! ! Modified: ! ! 24 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Output, real X, the corresponding argument. ! implicit none ! real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'ARCSIN_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if x = sin ( pi * cdf / 2.0E+00 ) return end subroutine arcsin_mean ( mean ) ! !******************************************************************************* ! !! ARCSIN_MEAN returns the mean of the Arcsin PDF. ! ! ! Modified: ! ! 17 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real mean ! mean = 0.5E+00 return end subroutine arcsin_pdf ( x, pdf ) ! !******************************************************************************* ! !! ARCSIN_PDF evaluates the Arcsin PDF. ! ! ! Discussion: ! ! The LOGISTIC EQUATION has the form: ! ! X(N+1) = 4.0E+00 * LAMBDA * ( 1.0E+00 - X(N) ). ! ! where 0 < LAMBDA <= 1. This nonlinear difference equation maps ! the unit interval into itself, and is a simple example of a system ! exhibiting chaotic behavior. Ulam and von Neumann studied the ! logistic equation with LAMBDA = 1, and showed that iterates of the ! function generated a sequence of pseudorandom numbers with ! the Arcsin probability density function. ! ! The derived sequence ! ! Y(N) = ( 2 / PI ) * Arcsin ( SQRT ( X(N) ) ) ! ! is a pseudorandom sequence with the uniform probability density ! function on [0,1]. For certain starting values, such as X(0) = 0, 0.75, ! or 1.0E+00, the sequence degenerates into a constant sequence, and for ! values very near these, the sequence takes a while before becoming ! chaotic. ! ! Formula: ! ! PDF(X) = 1 / ( PI * Sqrt ( X * ( 1 - X ) ) ) ! ! Reference: ! ! Daniel Zwillinger and Stephen Kokoska, ! CRC Standard Probability and Statistics Tables and Formulae, ! Chapman and Hall/CRC, 2000, pages 114-115. ! ! Modified: ! ! 24 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! 0.0E+00 < X < 1.0. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! if ( x <= 0.0E+00 .or. x >= 1.0E+00 ) then pdf = 0.0E+00 else pdf = 1.0E+00 / ( pi * sqrt ( x * ( 1.0E+00 - x ) ) ) end if return end subroutine arcsin_sample ( x ) ! !******************************************************************************* ! !! ARCSIN_SAMPLE samples the Arcsin PDF. ! ! ! Modified: ! ! 17 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real X, a sample of the PDF. ! implicit none ! real cdf real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call arcsin_cdf_inv ( cdf, x ) return end subroutine arcsin_variance ( variance ) ! !******************************************************************************* ! !! ARCSIN_VARIANCE returns the variance of the Arcsin PDF. ! ! ! Modified: ! ! 17 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real variance ! variance = 0.125E+00 return end subroutine benford_pdf ( x, pdf ) ! !******************************************************************************* ! !! BENFORD_PDF returns the Benford probability of one or more significant digits. ! ! ! Comments: ! ! Benford's law is an empirical formula explaining the observed ! distribution of initial digits in lists culled from newspapers, ! tax forms, stock market prices, and so on. It predicts the observed ! high frequency of the initial digit 1, for instance. ! ! Note that the probabilities of digits 1 through 9 are guaranteed ! to add up to 1, since ! LOG10 ( 2/1 ) + LOG10 ( 3/2) + LOG10 ( 4/3 ) + ... + LOG10 ( 10/9 ) ! = LOG10 ( 2/1 * 3/2 * 4/3 * ... * 10/9 ) = LOG10 ( 10 ) = 1. ! ! Formula: ! ! PDF(X) = LOG10 ( ( X + 1 ) / X ). ! ! Reference: ! ! F Benford, ! The Law of Anomalous Numbers, ! Proceedings of the American Philosophical Society, ! Volume 78, pages 551-572, 1938. ! ! T P Hill, ! The First Digit Phenomenon, ! American Scientist, ! Volume 86, July/August 1998, pages 358 - 363. ! ! R Raimi, ! The Peculiar Distribution of First Digits, ! Scientific American, ! December 1969, pages 109-119. ! ! Modified: ! ! 13 August 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the string of significant digits to be checked. ! If X is 1, then we are asking for the Benford probability that ! a value will have first digit 1. If X is 123, we are asking for ! the probability that the first three digits will be 123, and so on. ! ! Output, real PDF, the Benford probability that an item taken ! from a real world distribution will have the initial digits X. ! implicit none ! real pdf integer x ! if ( x <= 0 ) then pdf = 0.0E+00 else pdf = log10 ( real ( x + 1 ) / real ( x ) ) end if return end subroutine bernoulli_cdf ( x, a, cdf ) ! !******************************************************************************* ! !! BERNOULLI_CDF evaluates the Bernoulli CDF. ! ! ! Modified: ! ! 24 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the number of successes on a single trial. ! X = 0 or 1. ! ! Input, real A, the probability of success on one trial. ! 0.0E+00 <= A <= 1.0. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real cdf integer x ! if ( x < 0 ) then cdf = 0.0E+00 else if ( x == 0 ) then cdf = 1.0E+00 - a else cdf = 1.0E+00 end if return end subroutine bernoulli_cdf_inv ( cdf, a, x ) ! !******************************************************************************* ! !! BERNOULLI_CDF_INV inverts the Bernoulli CDF. ! ! ! Modified: ! ! 24 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, real A, the parameter of the PDF. ! 0.0E+00 <= A <= 1.0. ! ! Output, integer X, the corresponding argument. ! implicit none ! real a real cdf integer x ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BERNOULLI_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if if ( cdf <= 1.0E+00 - a ) then x = 0 else x = 1 end if return end subroutine bernoulli_check ( a ) ! !******************************************************************************* ! !! BERNOULLI_CHECK checks the parameter of the Bernoulli CDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the PDF. ! 0.0E+00 <= A <= 1.0. ! implicit none ! real a ! if ( a < 0.0E+00 .or. a > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BERNOULLI_CHECK - Fatal error!' write ( *, * ) ' A < 0 or 1 < A.' stop end if return end subroutine bernoulli_mean ( a, mean ) ! !******************************************************************************* ! !! BERNOULLI_MEAN returns the mean of the Bernoulli PDF. ! ! ! Modified: ! ! 16 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the probability of success. ! 0.0E+00 <= A <= 1.0. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real mean ! mean = a return end subroutine bernoulli_pdf ( x, a, pdf ) ! !******************************************************************************* ! !! BERNOULLI_PDF evaluates the Bernoulli PDF. ! ! ! Formula: ! ! PDF(X)(A) = A**X * ( 1.0E+00 - A )**( X - 1 ) ! ! X = 0 or 1. ! ! Discussion: ! ! The Bernoulli PDF describes the simple case in which a single trial ! is carried out, with two possible outcomes, called "success" and ! "failure"; the probability of success is A. ! ! Modified: ! ! 24 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the number of successes on a single trial. ! X = 0 or 1. ! ! Input, real A, the probability of success on one trial. ! 0.0E+00 <= A <= 1.0. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real pdf integer x ! if ( x < 0 ) then pdf = 0.0E+00 else if ( x == 0 ) then pdf = 1.0E+00 - a else if ( x == 1 ) then pdf = a else pdf = 0.0E+00 end if return end subroutine bernoulli_sample ( a, x ) ! !******************************************************************************* ! !! BERNOULLI_SAMPLE samples the Bernoulli PDF. ! ! ! Modified: ! ! 23 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the probability of success on one trial. ! 0.0E+00 <= A <= 1.0. ! ! Output, integer X, a sample of the PDF. ! implicit none ! real a real cdf integer x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call bernoulli_cdf_inv ( cdf, a, x ) return end subroutine bernoulli_variance ( a, variance ) ! !******************************************************************************* ! !! BERNOULLI_VARIANCE returns the variance of the Bernoulli PDF. ! ! ! Modified: ! ! 16 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the probability of success on one trial. ! 0.0E+00 <= A <= 1.0. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real variance ! variance = a * ( 1.0E+00 - a ) return end function beta ( a, b ) ! !******************************************************************************* ! !! BETA returns the value of the Beta function. ! ! ! Formula: ! ! BETA(A,B) = ( GAMMA ( A ) * GAMMA ( B ) ) / GAMMA ( A + B ) ! = Integral ( 0 <= T <= 1 ) T**(A-1) (1-T)**(B-1) dT. ! ! Modified: ! ! 10 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the function. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real BETA, the value of the function. ! implicit none ! real a real b real beta real gamma_log ! if ( a <= 0.0E+00 .or. b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA - Fatal error!' write ( *, * ) ' Both A and B must be greater than 0.' stop end if beta = exp ( gamma_log ( a ) + gamma_log ( b ) - gamma_log ( a + b ) ) return end subroutine beta_binomial_cdf ( x, a, b, c, cdf ) ! !******************************************************************************* ! !! BETA_BINOMIAL_CDF evaluates the Beta Binomial CDF. ! ! ! Discussion: ! ! A simple-minded summing approach is used. ! ! Modified: ! ! 07 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the argument of the CDF. ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real beta integer c real cdf real pdf integer x integer y ! if ( x < 0 ) then cdf = 0.0E+00 else if ( x < c ) then cdf = 0.0E+00 do y = 0, x pdf = beta ( a + real ( y ), b + real ( c - y ) ) / ( real ( c + 1 ) & * beta ( real ( y + 1), real ( c - y + 1 ) ) * beta ( a, b ) ) cdf = cdf + pdf end do else if ( x >= c ) then cdf = 1.0E+00 end if return end subroutine beta_binomial_cdf_inv ( cdf, a, b, c, x ) ! !******************************************************************************* ! !! BETA_BINOMIAL_CDF_INV inverts the Beta Binomial CDF. ! ! ! Discussion: ! ! A simple-minded discrete approach is used. ! ! Modified: ! ! 07 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! ! Output, integer X, the smallest X whose cumulative density function ! is greater than or equal to CDF. ! implicit none ! real a real b real beta integer c real cdf real cum real pdf integer x integer y ! if ( cdf <= 0.0E+00 ) then x = 0 else cum = 0.0E+00 do y = 0, c pdf = beta ( a + real ( y ), b + real ( c - y ) ) / ( real ( c + 1 ) & * beta ( real ( y + 1), real ( c - y + 1 ) ) * beta ( a, b ) ) cum = cum + pdf if ( cdf <= cum ) then x = y return end if end do x = c end if return end subroutine beta_binomial_check ( a, b, c ) ! !******************************************************************************* ! !! BETA_BINOMIAL_CHECK checks the parameters of the Beta Binomial PDF. ! ! ! Modified: ! ! 07 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! implicit none ! real a real b integer c ! if ( a <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!' write ( *, * ) ' A <= 0.' stop end if if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!' write ( *, * ) ' B <= 0.' stop end if if ( c < 0 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_BINOMIAL_CHECK - Fatal error!' write ( *, * ) ' C < 0.' stop end if return end subroutine beta_binomial_mean ( a, b, c, mean ) ! !******************************************************************************* ! !! BETA_BINOMIAL_MEAN returns the mean of the Beta Binomial PDF. ! ! ! Modified: ! ! 26 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= N. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b integer c real mean ! mean = real ( c ) * a / ( a + b ) return end subroutine beta_binomial_pdf ( x, a, b, c, pdf ) ! !******************************************************************************* ! !! BETA_BINOMIAL_PDF evaluates the Beta Binomial PDF. ! ! ! Formula: ! ! PDF(X)(A,B,C) = Beta(A+X,B+C-X) ! / ( (C+1) * Beta(X+1,C-X+1) * Beta(A,B) ) for 0 <= X <= C. ! ! This PDF can be reformulated as: ! ! The beta binomial probability density function for X successes ! out of N trials is ! ! PDF2(X)( N, MU, THETA ) = ! C(N,X) * Product ( 0 <= R <= X - 1 ) ( MU + R * THETA ) ! * Product ( 0 <= R <= N - X - 1 ) ( 1 - MU + R * THETA ) ! / Product ( 0 <= R <= N - 1 ) ( 1 + R * THETA ) ! ! where ! ! C(N,X) is the combinatorial coefficient; ! MU is the expectation of the underlying Beta distribution; ! THETA is a shape parameter. ! ! A THETA value of 0 ( or A+B --> Infinity ) results in the binomial ! distribution: ! ! PDF2(X) ( N, MU, 0 ) = C(N,X) * MU**X * ( 1 - MU )**(N-X) ! ! Given A, B, C for PDF, then the equivalent PDF2 has: ! ! N = C ! MU = A / ( A + B ) ! THETA = 1 / ( A + B ) ! ! Given N, MU, THETA for PDF2, the equivalent PDF has: ! ! A = MU / THETA ! B = ( 1 - MU ) / THETA ! C = N ! ! Discussion: ! ! BETA_BINOMIAL_PDF(X)(1,1,C) = UNIFORM_DISCRETE_PDF(X)(0,C-1) ! ! Modified: ! ! 18 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the argument of the PDF. ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real beta integer c real pdf integer x ! if ( x < 0 ) then pdf = 0.0E+00 else if ( x <= c ) then pdf = beta ( a + real ( x ), b + real ( c - x ) ) / ( real ( c + 1 ) & * beta ( real ( x + 1 ), real ( c - x + 1 ) ) * beta ( a, b ) ) else if ( x > c ) then pdf = 0.0E+00 end if return end subroutine beta_binomial_sample ( a, b, c, x ) ! !******************************************************************************* ! !! BETA_BINOMIAL_SAMPLE samples the Beta Binomial CDF. ! ! ! Modified: ! ! 07 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! ! Output, integer X, a sample of the PDF. ! implicit none ! real a real b integer c real cdf integer x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call beta_binomial_cdf_inv ( cdf, a, b, c, x ) return end subroutine beta_binomial_variance ( a, b, c, variance ) ! !******************************************************************************* ! !! BETA_BINOMIAL_VARIANCE returns the variance of the Beta Binomial PDF. ! ! ! Modified: ! ! 26 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, integer C, a parameter of the PDF. ! 0 <= C. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b integer c real variance ! variance = ( real ( c ) * a * b ) * ( a + b + real ( c ) ) & / ( ( a + b )**2 * ( a + b + 1.0E+00 ) ) return end subroutine beta_cdf ( x, a, b, cdf ) ! !******************************************************************************* ! !! BETA_CDF evaluates the Beta CDF. ! ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real beta_inc real cdf real x ! if ( x <= 0.0E+00 ) then cdf = 0.0E+00 else if ( x <= 1.0E+00 ) then cdf = beta_inc ( a, b, x ) else cdf = 1.0E+00 end if return end subroutine beta_cdf_inv ( cdf, a, b, x ) ! !******************************************************************************* ! !! BETA_CDF_INV inverts the Beta CDF. ! ! ! Reference: ! ! Algorithm 724, ! ACM Transactions on Mathematical Software, ! Volume 19, Number 4, December 1993, pages 481-483. ! ! Modified: ! ! 21 April 2001 ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real X, the argument of the CDF. ! implicit none ! integer, parameter :: maxk = 20 ! real a real b real bcoeff real cdf real cdf_x real d(2:maxk,0:maxk-2) real, parameter :: error = 0.0001E+00 real, parameter :: errapp = 0.01E+00 integer i integer j integer k integer loopct real pdf_x real q real s1 real s2 real sum2 real t real tail real x real xold ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if ! ! Estimate the solution. ! x = a / ( a + b ) xold = 0.0E+00 loopct = 2 do while ( abs ( ( x - xold ) / x ) >= errapp .and. loopct /= 0 ) xold = x loopct = loopct - 1 ! ! CDF_X = PROB { BETA(A,B) <= X }. ! Q = ( CDF - CDF_X ) / PDF_X. ! call beta_cdf ( x, a, b, cdf_x ) call beta_pdf ( x, a, b, pdf_x ) q = ( cdf - cdf_x ) / pdf_x ! ! D(N,K) = C(N,K) * Q**(N+K-1) / (N-1)! ! t = 1.0E+00 - x s1 = q * ( b - 1.0E+00 ) / t s2 = q * ( 1.0E+00 - a ) / x d(2,0) = s1 + s2 tail = d(2,0) * q / 2.0E+00 x = x + q + tail k = 3 do while ( abs ( tail / x ) > error .and. k <= maxk ) ! ! Find D(2,K-2). ! s1 = q * ( real ( k ) - 2.0E+00 ) * s1 / t s2 = q * ( 2.0E+00 - real ( k ) ) * s2 / x d(2,k-2) = s1 + s2 ! ! Find D(3,K-3), D(4,K-4), D(5,K-5), ... , D(K-1,1). ! do i = 3, k-1 sum2 = d(2,0) * d(i-1,k-i) bcoeff = 1.0E+00 do j = 1, k-i bcoeff = ( bcoeff * real ( k - i - j + 1 ) ) / real ( j ) sum2 = sum2 + bcoeff * d(2,j) * d(i-1,k-i-j) end do d(i,k-i) = sum2 + d(i-1,k-i+1) / real ( i - 1 ) end do ! ! Compute D(K,0) and use it to expand the series. ! d(k,0) = d(2,0) * d(k-1,0) + d(k-1,1) / real ( k - 1 ) tail = d(k,0) * q / real ( k ) x = x + tail ! ! Check for divergence. ! if ( x <= 0.0E+00 .or. x >= 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_CDF_INV - Fatal error!' write ( *, * ) ' The series has diverged.' write ( *, * ) ' X = ', x x = - 1.0E+00 return end if k = k + 1 end do end do return end subroutine beta_check ( a, b ) ! !******************************************************************************* ! !! BETA_CHECK checks the parameters of the Beta PDF. ! ! ! Modified: ! ! 08 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! implicit none ! real a real b ! if ( a <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_CHECK - Fatal error!' write ( *, * ) ' A <= 0.' stop end if if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_CHECK - Fatal error!' write ( *, * ) ' B <= 0.' stop end if return end function beta_inc ( a, b, x ) ! !******************************************************************************* ! !! BETA_INC returns the value of the incomplete Beta function. ! ! ! Formula: ! ! BETA_INC(A,B,X) ! ! = Integral ( 0 <= T <= X ) T**(A-1) (1-T)**(B-1) dT ! / Integral ( 0 <= T <= 1 ) T**(A-1) (1-T)**(B-1) dT ! ! = Integral ( 0 <= T <= X ) T**(A-1) (1-T)**(B-1) dT ! / BETA(A,B) ! ! Reference: ! ! Algorithm AS63, ! Applied Statistics, ! 1973, volume 22, number 3. ! ! Modified: ! ! 06 May 2001 ! ! Parameters: ! ! Input, real A, B, the parameters of the function. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Input, real X, the argument of the function. ! Normally, 0.0E+00 <= X <= 1.0. ! ! Output, real BETA_INC, the value of the function. ! implicit none ! real a real b real beta real beta_inc real cx integer i logical indx integer ns real pp real psq real qq real rx real temp real term real, parameter :: tol = 1.0E-07 real x real xx ! if ( a <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_INC - Fatal error!' write ( *, * ) ' A <= 0.' stop end if if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_INC - Fatal error!' write ( *, * ) ' B <= 0.' stop end if if ( x <= 0.0E+00 ) then beta_inc = 0.0E+00 return else if ( x >= 1.0E+00 ) then beta_inc = 1.0E+00 return end if ! ! Change tail if necessary and determine S. ! psq = a + b if ( a < ( a + b ) * x ) then xx = 1.0E+00 - x cx = x pp = b qq = a indx = .true. else xx = x cx = 1.0E+00 - x pp = a qq = b indx = .false. end if term = 1.0E+00 i = 1 beta_inc = 1.0E+00 ns = int ( qq + cx * ( a + b ) ) ! ! Use Soper's reduction formulas. ! rx = xx / cx temp = qq - real ( i ) if ( ns == 0 ) then rx = xx end if do term = term * temp * rx / ( pp + real ( i ) ) beta_inc = beta_inc + term temp = abs ( term ) if ( temp <= tol .and. temp <= tol * beta_inc ) then exit end if i = i + 1 ns = ns - 1 if ( ns >= 0 ) then temp = qq - real ( i ) if ( ns == 0 ) then rx = xx end if else temp = psq psq = psq + 1.0E+00 end if end do ! ! Finish calculation. ! beta_inc = beta_inc * exp ( pp * log ( xx ) + ( qq - 1.0E+00 ) * log ( cx ) ) & / ( beta ( a, b ) * pp ) if ( indx ) then beta_inc = 1.0E+00 - beta_inc end if return end subroutine beta_inc_values ( n, a, b, x, fx ) ! !******************************************************************************* ! !! BETA_INC_VALUES returns some values of the incomplete Beta function. ! ! ! Discussion: ! ! The incomplete Beta function may be written ! ! BETA_INC(A,B,X) = Integral (0 to X) T**(A-1) * (1-T)**(B-1) dT ! / BETA(A,B) ! ! Thus, ! ! BETA_INC(A,B,0.0) = 0.0 ! BETA_INC(A,B,1.0) = 1.0 ! ! Modified: ! ! 09 May 2001 ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer N. ! On input, if N is 0, the first test data is returned, and N is set ! to the index of the test data. On each subsequent call, N is ! incremented and that test data is returned. When there is no more ! test data, N is set to 0. ! ! Output, real A, B, X, the arguments of the function. ! ! Output, real FX, the value of the function. ! implicit none ! integer, parameter :: nmax = 20 ! real a real, save, dimension ( nmax ) :: avec = (/ & 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, & 1.0E+00, 1.0E+00, 1.0E+00, 5.0E+00, & 10.0E+00, 10.0E+00, 10.0E+00, 10.0E+00, & 20.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, & 20.0E+00, 30.0E+00, 30.0E+00, 40.0E+00 /) real b real, save, dimension ( nmax ) :: bvec = (/ & 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, & 0.5E+00, 0.5E+00, 1.0E+00, 5.0E+00, & 0.5E+00, 5.0E+00, 5.0E+00, 10.0E+00, & 5.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, & 20.0E+00, 10.0E+00, 10.0E+00, 20.0E+00 /) real fx real, save, dimension ( nmax ) :: fxvec = (/ & 0.0637686E+00, 0.2048328E+00, 1.0000000E+00, 0.0050126E+00, & 0.0513167E+00, 1.0000000E+00, 0.5000000E+00, 0.5000000E+00, & 0.1516409E+00, 0.0897827E+00, 1.0000000E+00, 0.5000000E+00, & 0.4598773E+00, 0.2146816E+00, 0.9507365E+00, 0.5000000E+00, & 0.8979414E+00, 0.2241297E+00, 0.7586405E+00, 0.7001783E+00 /) integer n real x real, save, dimension ( nmax ) :: xvec = (/ & 0.01E+00, 0.10E+00, 1.00E+00, 0.01E+00, & 0.10E+00, 1.00E+00, 0.50E+00, 0.50E+00, & 0.90E+00, 0.50E+00, 1.00E+00, 0.50E+00, & 0.80E+00, 0.60E+00, 0.80E+00, 0.50E+00, & 0.60E+00, 0.70E+00, 0.80E+00, 0.70E+00 /) ! if ( n < 0 ) then n = 0 end if n = n + 1 if ( n > nmax ) then n = 0 a = 0.0E+00 b = 0.0E+00 x = 0.0E+00 fx = 0.0E+00 return end if a = avec(n) b = bvec(n) x = xvec(n) fx = fxvec(n) return end subroutine beta_mean ( a, b, mean ) ! !******************************************************************************* ! !! BETA_MEAN returns the mean of the Beta PDF. ! ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b real mean ! mean = a / ( a + b ) return end subroutine beta_pascal_cdf ( x, a, b, c, cdf ) ! !******************************************************************************* ! !! BETA_PASCAL_CDF evaluates the Beta Pascal CDF. ! ! ! Discussion: ! ! A simple-minded summing approach is used. ! ! Modified: ! ! 09 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the argument of the CDF. ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, real CDF, the value of the CDF. ! implicit none ! integer a real b real c real cdf real pdf integer x integer y ! cdf = 0.0E+00 do y = a, x call beta_pascal_pdf ( y, a, b, c, pdf ) cdf = cdf + pdf end do return end subroutine beta_pascal_cdf_inv ( cdf, a, b, c, x ) ! !******************************************************************************* ! !! BETA_PASCAL_CDF_INV inverts the Beta Pascal CDF. ! ! ! Discussion: ! ! A simple-minded discrete approach is used. ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, integer X, the smallest X whose cumulative density function ! is greater than or equal to CDF. ! implicit none ! integer a real b real c real cdf real cum real pdf integer x integer, parameter :: xmax = 1000 ! ! if ( cdf <= 0.0E+00 ) then x = a else if ( cdf < 1.0E+00 ) then x = a call beta_pascal_pdf ( x, a, b, c, pdf ) cum = pdf do while ( cdf > cum .and. x < xmax ) x = x + 1 call beta_pascal_pdf ( x, a, b, c, pdf ) cum = cum + pdf end do else write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_CDF_INV - Fatal error!' write ( *, * ) ' Input CDF >= 1.' stop end if return end subroutine beta_pascal_check ( a, b, c ) ! !******************************************************************************* ! !! BETA_PASCAL_CHECK checks the parameters of the Beta Pascal PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! implicit none ! integer a real b real c ! if ( a <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_CHECK - Fatal error!' write ( *, * ) ' A <= 0.' stop end if if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_CHECK - Fatal error!' write ( *, * ) ' B <= 0.' stop end if if ( c <= b ) then write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_CHECK - Fatal error!' write ( *, * ) ' C <= B.' stop end if return end subroutine beta_pascal_mean ( a, b, c, mean ) ! !******************************************************************************* ! !! BETA_PASCAL_MEAN returns the mean of the Beta Pascal PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! integer a real b real c real mean ! if ( b <= 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_MEAN - Fatal error!' write ( *, * ) ' The mean is undefined for B <= 1.' stop end if mean = ( real ( a ) * ( c - 1.0E+00 ) ) / ( b - 1.0E+00 ) return end subroutine beta_pascal_pdf ( x, a, b, c, pdf ) ! !******************************************************************************* ! !! BETA_PASCAL_PDF evaluates the Beta Pascal PDF. ! ! ! Discussion: ! ! Something is wrong. If A = C, the results seem OK, but as ! A increases, the PDF gets too small. ! (5,3,4), the CDF seems to get stuck at 0.1428 = 1/7 ! (6,3,4), the CDF gets stuck at 1/56. ! ! Formula: ! ! PDF(X)(A,B,C) = Gamma(X) * Gamma(C) * Gamma(B+C) * Gamma(X-A+C-B) ! / ( Gamma(A) * Gamma(X-A+1) * Gamma(B) * Gamma(C-B) * Gamma(X+C) ) ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the argument of the PDF. ! A <= X. ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! integer a real b real c real gamma_log real pdf real pdf_log real ra real rx integer x ! if ( x < a ) then pdf = 0.0E+00 else ra = real ( a ) rx = real ( x ) pdf_log = & + gamma_log ( rx ) & + gamma_log ( c ) & + gamma_log ( b + c ) & + gamma_log ( rx - ra + c - b ) & - gamma_log ( ra ) & - gamma_log ( rx + 1.0E+00 - ra ) & - gamma_log ( b ) & - gamma_log ( c - b ) & - gamma_log ( rx + c ) pdf = exp ( pdf_log ) end if return end subroutine beta_pascal_sample ( a, b, c, x ) ! !******************************************************************************* ! !! BETA_PASCAL_SAMPLE samples the Beta Pascal CDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, integer X, a sample of the PDF. ! implicit none ! integer a real b real c real cdf integer x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call beta_pascal_cdf_inv ( cdf, a, b, c, x ) return end subroutine beta_pascal_variance ( a, b, c, variance ) ! !******************************************************************************* ! !! BETA_PASCAL_VARIANCE returns the variance of the Beta Pascal PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, real B, C, parameters of the PDF. ! 0 < A, ! 0.0E+00 < B < C. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! integer a real b real bot real c real top real variance ! if ( b <= 2.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BETA_PASCAL_VARIANCE - Fatal error!' write ( *, * ) ' The mean is undefined for B <= 2.' stop end if top = real ( a ) * ( real ( a ) + b - 1.0E+00 ) * ( c - 1.0E+00 ) * ( c - b ) bot = ( b - 1.0E+00 )**2 * ( b - 2.0E+00 ) variance = top / bot return end subroutine beta_pdf ( x, a, b, pdf ) ! !******************************************************************************* ! !! BETA_PDF evaluates the Beta PDF. ! ! ! Formula: ! ! PDF(X)(A,B) = X**(A-1) * (1-X)**(B-1) / BETA(A,B). ! ! Discussion: ! ! A = B = 1 yields the Uniform distribution on [0,1]. ! A = B = 1/2 yields the Arcsin distribution. ! B = 1 yields the power function distribution. ! A = B -> Infinity tends to the Normal distribution. ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! 0.0E+00 <= X <= 1.0. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real beta real pdf real x ! if ( x < 0.0E+00 .or. x > 1.0E+00 ) then pdf = 0.0E+00 else pdf = x**( a - 1.0E+00 ) * ( 1.0E+00 - x )**( b - 1.0E+00 ) / beta ( a, b ) end if return end subroutine beta_sample ( a, b, x ) ! !******************************************************************************* ! !! BETA_SAMPLE samples the Beta PDF. ! ! ! Reference: ! ! Algorithm BN, ! William Kennedy and James Gentle, ! Statistical Computing, ! Dekker, 1980. ! ! Modified: ! ! 05 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real mu real stdev real test real u real x real y ! mu = ( a - 1.0E+00 ) / ( a + b - 2.0E+00 ) stdev = 0.5E+00 / sqrt ( a + b - 2.0E+00 ) do call normal_01_sample ( y ) x = mu + stdev * y if ( x < 0.0E+00 .or. x > 1.0E+00 ) then cycle end if call r_random ( 0.0E+00, 1.0E+00, u ) test = ( a - 1.0E+00 ) * log ( x / ( a - 1.0E+00 ) ) & + ( b - 1.0E+00 ) * log ( ( 1.0E+00 - x ) / ( b - 1.0E+00 ) ) & + ( a + b - 2.0E+00 ) * log ( a + b - 2.0E+00 ) + 0.5E+00 * y**2 if ( log ( u ) <= test ) then exit end if end do return end subroutine beta_variance ( a, b, variance ) ! !******************************************************************************* ! !! BETA_VARIANCE returns the variance of the Beta PDF. ! ! ! Modified: ! ! 27 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < A, ! 0.0E+00 < B. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b real variance ! variance = ( a * b ) / ( ( a + b )**2 * ( 1.0E+00 + a + b ) ) return end subroutine binomial_cdf ( x, a, b, cdf ) ! !******************************************************************************* ! !! BINOMIAL_CDF evaluates the Binomial CDF. ! ! ! Definition: ! ! CDF(X)(A,B) is the probability of at most X successes in A trials, ! given that the probability of success on a single trial is B. ! ! Discussion: ! ! A sequence of trials with fixed probability of success on ! any trial is known as a sequence of Bernoulli trials. ! ! Modified: ! ! 28 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the desired number of successes. ! 0 <= X <= A. ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, real CDF, the value of the CDF. ! implicit none ! integer a real b integer cnk real cdf integer j real pr integer x ! if ( x < 0 ) then cdf = 0.0E+00 else if ( x >= a ) then cdf = 1.0E+00 else if ( b == 0.0E+00 ) then cdf = 1.0E+00 else if ( b == 1.0E+00 ) then cdf = 0.0E+00 else cdf = 0.0E+00 do j = 0, x call binomial_coef ( a, j, cnk ) pr = real ( cnk ) * b**j * ( 1.0E+00 - b )**( a - j ) cdf = cdf + pr end do end if return end subroutine binomial_cdf_values ( n, a, b, x, fx ) ! !******************************************************************************* ! !! BINOMIAL_CDF_VALUES returns some values of the binomial CDF. ! ! ! Discussion: ! ! CDF(X)(A,B) is the probability of at most X successes in A trials, ! given that the probability of success on a single trial is B. ! ! Modified: ! ! 27 May 2001 ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Daniel Zwillinger, ! CRC Standard Mathematical Tables and Formulae, ! 30th Edition, CRC Press, 1996, pages 651-652. ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer N. ! On input, if N is 0, the first test data is returned, and N is set ! to the index of the test data. On each subsequent call, N is ! incremented and that test data is returned. When there is no more ! test data, N is set to 0. ! ! Output, integer A, real B, integer X, the arguments of the function. ! ! Output, real FX, the value of the function. ! implicit none ! integer, parameter :: nmax = 17 ! integer a integer, save, dimension ( nmax ) :: avec = (/ & 2, 2, 2, 2, & 2, 4, 4, 4, & 4, 10, 10, 10, & 10, 10, 10, 10, & 10 /) real b real, save, dimension ( nmax ) :: bvec = (/ & 0.05E+00, 0.05E+00, 0.05E+00, 0.50E+00, & 0.50E+00, 0.25E+00, 0.25E+00, 0.25E+00, & 0.25E+00, 0.05E+00, 0.10E+00, 0.15E+00, & 0.20E+00, 0.25E+00, 0.30E+00, 0.40E+00, & 0.50E+00 /) real fx real, save, dimension ( nmax ) :: fxvec = (/ & 0.9025E+00, 0.9975E+00, 1.0000E+00, 0.2500E+00, & 0.7500E+00, 0.3164E+00, 0.7383E+00, 0.9492E+00, & 0.9961E+00, 0.9999E+00, 0.9984E+00, 0.9901E+00, & 0.9672E+00, 0.9219E+00, 0.8497E+00, 0.6331E+00, & 0.3770E+00 /) integer n integer x integer, save, dimension ( nmax ) :: xvec = (/ & 0, 1, 2, 0, & 1, 0, 1, 2, & 3, 4, 4, 4, & 4, 4, 4, 4, & 4 /) ! if ( n < 0 ) then n = 0 end if n = n + 1 if ( n > nmax ) then n = 0 a = 0 b = 0.0E+00 x = 0 fx = 0.0E+00 return end if a = avec(n) b = bvec(n) x = xvec(n) fx = fxvec(n) return end subroutine binomial_cdf_inv ( cdf, a, b, x ) ! !******************************************************************************* ! !! BINOMIAL_CDF_INV inverts the Binomial CDF. ! ! ! Modified: ! ! 06 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, integer X, the corresponding argument. ! implicit none ! integer a real b real cdf real cdf2 real pdf integer x integer x2 ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BINOMIAL_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if cdf2 = 0.0E+00 do x2 = 0, a call binomial_pdf ( x2, a, b, pdf ) cdf2 = cdf2 + pdf if ( cdf <= cdf2 ) then x = x2 return end if end do return end subroutine binomial_check ( a, b ) ! !******************************************************************************* ! !! BINOMIAL_CHECK checks the parameter of the Binomial PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! implicit none ! integer a real b ! if ( a < 1 ) then write ( *, * ) ' ' write ( *, * ) 'BINOMIAL_CHECK - Fatal error!' write ( *, * ) ' A < 1.' stop end if if ( b < 0.0E+00 .or. b > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BINOMIAL_CHECK - Fatal error!' write ( *, * ) ' B < 0 or 1 < B.' stop end if return end subroutine binomial_coef ( n, k, cnk ) ! !******************************************************************************* ! !! BINOMIAL_COEF computes the Binomial coefficient C(N,K). ! ! ! Discussion: ! ! The value is calculated in such a way as to avoid overflow and ! roundoff. The calculation is done in integer arithmetic. ! ! Formula: ! ! CNK = C(N,K) = N! / ( K! * (N-K)! ) ! ! Reference: ! ! M L Wolfson and H V Wright, ! Combinatorial of M Things Taken N at a Time, ! ACM Algorithm 160, ! Communications of the ACM, ! April, 1963. ! ! Modified: ! ! 17 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, K, are the values of N and K. ! ! Output, integer CNK, the number of combinations of N ! things taken K at a time. ! implicit none ! integer cnk integer i integer k integer mn integer mx integer n ! mn = min ( k, n-k ) if ( mn < 0 ) then cnk = 0 else if ( mn == 0 ) then cnk = 1 else mx = max ( k, n-k ) cnk = mx + 1 do i = 2, mn cnk = ( cnk * ( mx + i ) ) / i end do end if return end subroutine binomial_coef_log ( n, k, cnk_log ) ! !******************************************************************************* ! !! BINOMIAL_COEF_LOG computes the logarithm of the Binomial coefficient. ! ! ! Formula: ! ! CNK_LOG = LOG ( C(N,K) ) = LOG ( N! / ( K! * (N-K)! ) ). ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, K, are the values of N and K. ! ! Output, real CNK_LOG, the logarithm of C(N,K). ! implicit none ! real cnk_log real factorial_log integer k integer n ! cnk_log = factorial_log ( n ) - factorial_log ( k ) - factorial_log ( n - k ) return end subroutine binomial_mean ( a, b, mean ) ! !******************************************************************************* ! !! BINOMIAL_MEAN returns the mean of the Binomial PDF. ! ! ! Modified: ! ! 28 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, real MEAN, the expected value of the number of ! successes in A trials. ! implicit none ! integer a real b real mean ! mean = real ( a ) * b return end subroutine binomial_pdf ( x, a, b, pdf ) ! !******************************************************************************* ! !! BINOMIAL_PDF evaluates the Binomial PDF. ! ! ! Definition: ! ! PDF(X)(A,B) is the probability of exactly X successes in A trials, ! given that the probability of success on a single trial is B. ! ! Formula: ! ! PDF(X)(A,B) = C(N,X) * B**X * ( 1.0E+00 - B )**( A - X ) ! ! Discussion: ! ! Binomial_PDF(X)(1,B) = Bernoulli_PDF(X)(B). ! ! Modified: ! ! 28 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the desired number of successes. ! 0 <= X <= A. ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, real PDF, the value of the PDF. ! implicit none ! integer a real b integer cnk real pdf integer x ! if ( a < 1 ) then pdf = 0.0E+00 else if ( x < 0 .or. x > a ) then pdf = 0.0E+00 else if ( b == 0.0E+00 ) then if ( x == 0 ) then pdf = 1.0E+00 else pdf = 0.0E+00 end if else if ( b == 1.0E+00 ) then if ( x == a ) then pdf = 1.0E+00 else pdf = 0.0E+00 end if else call binomial_coef ( a, x, cnk ) pdf = real ( cnk ) * b**x * ( 1.0E+00 - b )**( a - x ) end if return end subroutine binomial_sample ( a, b, x ) ! !******************************************************************************* ! !! BINOMIAL_SAMPLE samples the Binomial PDF. ! ! ! Reference: ! ! Algorithm BU, ! William Kennedy and James Gentle, ! Statistical Computing, ! Dekker, 1980. ! ! Modified: ! ! 02 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, integer X, a sample of the PDF. ! implicit none ! integer a real b integer i real u integer x ! x = 0 do i = 1, a call r_random ( 0.0E+00, 1.0E+00, u ) if ( u <= b ) then x = x + 1 end if end do return end subroutine binomial_variance ( a, b, variance ) ! !******************************************************************************* ! !! BINOMIAL_VARIANCE returns the variance of the Binomial PDF. ! ! ! Modified: ! ! 28 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of trials. ! 1 <= A. ! ! Input, real B, the probability of success on one trial. ! 0.0E+00 <= B <= 1.0. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! integer a real b real variance ! variance = real ( a ) * b * ( 1.0E+00 - b ) return end subroutine bradford_cdf ( x, a, b, c, cdf ) ! !******************************************************************************* ! !! BRADFORD_CDF evaluates the Bradford CDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real c real cdf real x ! if ( x <= a ) then cdf = 0.0E+00 else if ( x <= b ) then cdf = log ( 1.0E+00 + c * ( x - a ) / ( b - a ) ) / log ( c + 1.0E+00 ) else if ( x > b ) then cdf = 1.0E+00 end if return end subroutine bradford_cdf_inv ( cdf, a, b, c, x ) ! !******************************************************************************* ! !! BRADFORD_CDF_INV inverts the Bradford CDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real X, the corresponding argument of the CDF. ! implicit none ! real a real b real c real cdf real x ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BRADFORD_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if if ( cdf <= 0.0E+00 ) then x = a else if ( cdf < 1.0E+00 ) then x = a + ( b - a ) * ( ( c + 1.0E+00 )**cdf - 1.0E+00 ) / c else if ( cdf >= 1.0E+00 ) then x = b end if return end subroutine bradford_check ( a, b, c ) ! !******************************************************************************* ! !! BRADFORD_CHECK checks the parameters of the Bradford PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! implicit none ! real a real b real c ! if ( a >= b ) then write ( *, * ) ' ' write ( *, * ) 'BRADFORD_CHECK - Fatal error!' write ( *, * ) ' A >= B.' stop end if if ( c <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BRADFORD_CHECK - Fatal error!' write ( *, * ) ' C <= 0.' stop end if return end subroutine bradford_mean ( a, b, c, mean ) ! !******************************************************************************* ! !! BRADFORD_MEAN returns the mean of the Bradford PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b real c real mean ! mean = ( c * ( b - a ) + log ( c + 1.0E+00 ) * ( a * ( c + 1.0E+00 ) - b ) ) & / ( c * log ( c + 1.0E+00 ) ) return end subroutine bradford_pdf ( x, a, b, c, pdf ) ! !******************************************************************************* ! !! BRADFORD_PDF evaluates the Bradford PDF. ! ! ! Formula: ! ! PDF(X)(A,B,C) = ! C / ( ( C * ( X - A ) + B - A ) * log ( C + 1 ) ) ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! A <= X ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real c real pdf real x ! if ( x <= a ) then pdf = 0.0E+00 else if ( x <= b ) then pdf = c / ( ( c * ( x - a ) + b - a ) * log ( c + 1.0E+00 ) ) else if ( x > b ) then pdf = 0.0E+00 end if return end subroutine bradford_sample ( a, b, c, x ) ! !******************************************************************************* ! !! BRADFORD_SAMPLE samples the Bradford PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real c real cdf real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) x = a + ( b - a ) * ( ( c + 1.0E+00 )**cdf - 1.0E+00 ) / c return end subroutine bradford_variance ( a, b, c, variance ) ! !******************************************************************************* ! !! BRADFORD_VARIANCE returns the variance of the Bradford PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! A < B, ! 0.0E+00 < C. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b real c real variance ! variance = ( b - a )**2 * & ( c * ( log ( c + 1.0E+00 ) - 2.0E+00 ) + 2.0E+00 * log ( c + 1.0E+00 ) ) & / ( 2.0E+00 * c * ( log ( c + 1.0E+00 ) )**2 ) return end subroutine burr_cdf ( x, a, b, c, d, cdf ) ! !******************************************************************************* ! !! BURR_CDF evaluates the Burr CDF. ! ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real c real cdf real d real x ! if ( x <= a ) then cdf = 0.0E+00 else cdf = 1.0E+00 / ( 1.0E+00 + ( b / ( x - a ) )**c )**d end if return end subroutine burr_cdf_inv ( cdf, a, b, c, d, x ) ! !******************************************************************************* ! !! BURR_CDF_INV inverts the Burr CDF. ! ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real X, the corresponding argument. ! implicit none ! real a real b real c real cdf real d real x ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BURR_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if x = a + b / ( ( 1.0E+00 / cdf )**(1.0E+00 / d ) - 1.0E+00 )**( 1.0E+00 / c ) return end subroutine burr_check ( a, b, c, d ) ! !******************************************************************************* ! !! BURR_CHECK checks the parameters of the Burr CDF. ! ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! implicit none ! real a real b real c real d ! if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BURR_CHECK - Fatal error!' write ( *, * ) ' B <= 0.' stop end if if ( c <= 0 ) then write ( *, * ) ' ' write ( *, * ) 'BURR_CHECK - Fatal error!' write ( *, * ) ' C <= 0.' stop end if return end subroutine burr_mean ( a, b, c, d, mean ) ! !******************************************************************************* ! !! BURR_MEAN returns the mean of the Burr PDF. ! ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b real c real d real gamma real mean ! mean = a + b * gamma ( 1.0E+00 - 1.0E+00 / c ) & * gamma ( d + 1.0E+00 / c ) / gamma ( d ) return end subroutine burr_pdf ( x, a, b, c, d, pdf ) ! !******************************************************************************* ! !! BURR_PDF evaluates the Burr PDF. ! ! ! Formula: ! ! PDF(X)(A,B,C,D) = ( C * D / B ) * ( ( X - A ) / B )**( - C - 1 ) ! * ( 1 + ( ( X - A ) / B )**( - C ) )**( - D - 1 ) ! ! Reference: ! ! M E Johnson, ! Multivariate Statistical Simulation, ! Wiley, New York, 1987. ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! A <= X ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real c real d real pdf real x real y ! if ( x <= a ) then pdf = 0.0E+00 else y = ( x - a ) / b pdf = ( c * d / b ) * y**( - c - 1.0E+00 ) & * ( 1.0E+00 + y**( - c ) )**( - d - 1.0E+00 ) end if return end subroutine burr_sample ( a, b, c, d, x ) ! !******************************************************************************* ! !! BURR_SAMPLE samples the Burr PDF. ! ! ! Modified: ! ! 01 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real c real cdf real d real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call burr_cdf_inv ( cdf, a, b, c, d, x ) return end subroutine burr_variance ( a, b, c, d, variance ) ! !******************************************************************************* ! !! BURR_VARIANCE returns the variance of the Burr PDF. ! ! ! Modified: ! ! 02 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, D, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b real c real d real gamma real k real variance ! if ( c <= 2.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'BURR_VARIANCE - Warning!' write ( *, * ) ' Variance undefined for C <= 2.' variance = huge ( variance ) else k = gamma ( d ) * gamma ( 1.0E+00 - 2.0E+00 / c ) & * gamma ( d + 2.0E+00 / c ) & - ( gamma ( 1.0E+00 - 1.0E+00 / c ) * gamma ( d + 1.0E+00 / c ) )**2 variance = k * b**2 / ( gamma ( d ) )**2 end if return end subroutine c_normal_01_sample ( x ) ! !******************************************************************************* ! !! C_NORMAL_01_SAMPLE samples the complex Normal 01 PDF. ! ! ! Modified: ! ! 30 May 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, complex X, a sample of the PDF. ! implicit none ! real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real v1 real v2 complex x real x_c real x_r ! call random_number ( harvest = v1 ) call random_number ( harvest = v2 ) x_r = sqrt ( - 2.0E+00 * log ( v1 ) ) * cos ( 2.0E+00 * pi * v2 ) x_c = sqrt ( - 2.0E+00 * log ( v1 ) ) * sin ( 2.0E+00 * pi * v2 ) x = cmplx ( x_r, x_c ) return end subroutine cauchy_cdf ( x, a, b, cdf ) ! !******************************************************************************* ! !! CAUCHY_CDF evaluates the Cauchy CDF. ! ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x real y ! y = ( x - a ) / b cdf = 0.5E+00 + atan ( y ) / PI return end subroutine cauchy_cdf_inv ( cdf, a, b, x ) ! !******************************************************************************* ! !! CAUCHY_CDF_INV inverts the Cauchy CDF. ! ! ! Modified: ! ! 12 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real X, the corresponding argument. ! implicit none ! real a real b real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! x = a + b * tan ( pi * ( cdf - 0.5E+00 ) ) return end subroutine cauchy_check ( a, b ) ! !******************************************************************************* ! !! CAUCHY_CHECK checks the parameters of the Cauchy CDF. ! ! ! Modified: ! ! 09 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! implicit none ! real a real b ! if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CAUCHY_CHECK - Fatal error!' write ( *, * ) ' B <= 0.' stop end if return end subroutine cauchy_mean ( a, b, mean ) ! !******************************************************************************* ! !! CAUCHY_MEAN returns the mean of the Cauchy PDF. ! ! ! Modified: ! ! 12 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b real mean ! mean = a return end subroutine cauchy_pdf ( x, a, b, pdf ) ! !******************************************************************************* ! !! CAUCHY_PDF evaluates the Cauchy PDF. ! ! ! Formula: ! ! PDF(X)(A,B) = 1 / ( PI * B * ( 1 + ( ( X - A ) / B )**2 ) ) ! ! Discussion: ! ! The Cauchy PDF is also known as the Breit-Wigner PDF. It ! has some unusual properties. In particular, the integrals for the ! expected value and higher order moments are "singular", in the ! sense that the limiting values do not exist. A result can be ! obtained if the upper and lower limits of integration are set ! equal to +T and -T, and the limit as T=>INFINITY is taken, but ! this is a very weak and unreliable sort of limit. ! ! Modified: ! ! 09 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x real y ! y = ( x - a ) / b pdf = 1.0E+00 / ( pi * b * ( 1.0E+00 + y**2 ) ) return end subroutine cauchy_sample ( a, b, x ) ! !******************************************************************************* ! !! CAUCHY_SAMPLE samples the Cauchy PDF. ! ! ! Modified: ! ! 11 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real cdf real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call cauchy_cdf_inv ( cdf, a, b, x ) return end subroutine cauchy_variance ( a, b, variance ) ! !******************************************************************************* ! !! CAUCHY_VARIANCE returns the variance of the Cauchy PDF. ! ! ! Discussion: ! ! The variance of the Cauchy PDF is not well defined. This routine ! is made available for completeness only, and simply returns ! a "very large" number. ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real VARIANCE, the mean of the PDF. ! implicit none ! real a real b real variance ! variance = huge ( variance ) return end subroutine chi_cdf ( x, a, b, c, cdf ) ! !******************************************************************************* ! !! CHI_CDF evaluates the Chi CDF. ! ! ! Modified: ! ! 05 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real c real cdf real gamma_inc real p2 real x real x2 real y ! if ( x <= a ) then cdf = 0.0E+00 else y = ( x - a ) / b x2 = 0.5E+00 * y**2 p2 = 0.5E+00 * c cdf = gamma_inc ( p2, x2 ) end if return end subroutine chi_cdf_inv ( cdf, a, b, c, x ) ! !******************************************************************************* ! !! CHI_CDF_INV inverts the Chi CDF. ! ! ! Discussion: ! ! A simple bisection method is used. ! ! Modified: ! ! 05 January 2000 ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real X, the corresponding argument of the CDF. ! implicit none ! real a real b real c real cdf real cdf1 real cdf2 real cdf3 integer it integer, parameter :: it_max = 10 real, parameter :: tol = 0.0001E+00 real x real x1 real x2 real x3 ! if ( cdf <= 0.0E+00 ) then x = a return else if ( cdf >= 1.0E+00 ) then x = huge ( x ) return end if x1 = a cdf1 = 0.0E+00 x2 = a + 1.0E+00 do call chi_cdf ( x2, a, b, c, cdf2 ) if ( cdf2 > cdf ) then exit end if x2 = a + 2.0E+00 * ( x2 - a ) end do ! ! Now use bisection. ! it = 0 do it = it + 1 x3 = 0.5E+00 * ( x1 + x2 ) call chi_cdf ( x3, a, b, c, cdf3 ) if ( abs ( cdf3 - cdf ) < tol ) then x = x3 return end if if ( it > it_max ) then write ( *, * ) ' ' write ( *, * ) 'CHI_CDF_INV - Fatal error!' write ( *, * ) ' Iteration limit exceeded.' stop end if if ( sign ( 1.0E+00, cdf3 - cdf ) == sign ( 1.0E+00, cdf1 - cdf ) ) then x1 = x3 cdf1 = cdf3 else x2 = x3 cdf2 = cdf3 end if end do return end subroutine chi_check ( a, b, c ) ! !******************************************************************************* ! !! CHI_CHECK checks the parameters of the Chi CDF. ! ! ! Modified: ! ! 31 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! implicit none ! real a real b real c ! if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CHI_CHECK - Fatal error!' write ( *, * ) ' B <= 0.0.' stop end if if ( c <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CHI_CHECK - Fatal error!' write ( *, * ) ' C <= 0.0.' stop end if return end subroutine chi_mean ( a, b, c, mean ) ! !******************************************************************************* ! !! CHI_MEAN returns the mean of the Chi PDF. ! ! ! Modified: ! ! 31 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real MEAN, the mean value. ! implicit none ! real a real b real c real gamma real mean ! mean = a + sqrt ( 2.0E+00 ) * b * gamma ( 0.5E+00 * ( c + 1.0E+00 ) ) & / gamma ( 0.5E+00 * c ) return end subroutine chi_pdf ( x, a, b, c, pdf ) ! !******************************************************************************* ! !! CHI_PDF evaluates the Chi PDF. ! ! ! Formula: ! ! PDF(X)(A,B,C) = EXP ( - 0.5E+00 * ( ( X - A ) / B )**2 ) ! * ( ( X - A ) / B )**( C - 1 ) / ! ( 2**( 0.5E+00 * C - 1 ) * B * GAMMA ( 0.5E+00 * C ) ) ! ! Discussion: ! ! CHI(A,B,1) is the Half Normal PDF; ! CHI(0,B,2) is the Rayleigh PDF; ! CHI(0,B,3) is the Maxwell PDF. ! ! Modified: ! ! 31 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! A <= X ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real c real gamma real pdf real x real y ! if ( x <= a ) then pdf = 0.0E+00 else y = ( x - a ) / b pdf = exp ( - 0.5E+00 * y**2 ) * y**( c - 1.0E+00 ) / & ( 2.0E+00**( 0.5E+00 * c - 1.0E+00 ) * b * gamma ( 0.5E+00 * c ) ) end if return end subroutine chi_sample ( a, b, c, x ) ! !******************************************************************************* ! !! CHI_SAMPLE samples the Chi PDF. ! ! ! Modified: ! ! 02 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real c real x ! call chisquare_central_sample ( c, x ) x = a + b * sqrt ( x ) return end subroutine chi_variance ( a, b, c, variance ) ! !******************************************************************************* ! !! CHI_VARIANCE returns the variance of the Chi PDF. ! ! ! Modified: ! ! 31 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the PDF. ! 0 < B, ! 0 < C. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b real c real gamma real variance ! variance = b**2 * ( c - 2.0E+00 * & ( gamma ( 0.5E+00 * ( c + 1.0E+00 ) ) / gamma ( 0.5E+00 * c ) )**2 ) return end subroutine chisquare_central_cdf ( x, a, cdf ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_CDF evaluates the chi-squared CDF. ! ! ! Modified: ! ! 30 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the value of the random deviate. ! ! Input, real A, the parameter of the distribution, usually ! the number of degrees of freedom. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real a2 real b2 real c2 real cdf real x real x2 ! x2 = 0.5E+00 * x a2 = 0.0E+00 b2 = 1.0E+00 c2 = 0.5E+00 * a call gamma_cdf ( x2, a2, b2, c2, cdf ) cdf = 1.0E+00 - cdf return end subroutine chisquare_central_cdf_inv ( cdf, a, x ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_CDF_INV inverts the chi-squared PDF. ! ! ! Reference: ! ! Best and Roberts, ! The Percentage Points of the Chi-Squared Distribution, ! Algorithm AS 91, ! Applied Statistics, ! Volume 24, Number ?, pages 385-390, 1975. ! ! Modified: ! ! 30 November 1999 ! ! Parameters: ! ! Input, real CDF, a value of the chi-squared cumulative probability ! density function. ! 0.000002 <= CDF <= 0.999998. ! ! Input, real A, the parameter of the chi-squared probability density ! function. A > 0. ! ! Output, real X, the value of the chi-squared random deviate ! with the property that the probability that a chi-squared random ! deviate with parameter A is less than or equal to PPCHI2 is P. ! implicit none ! real a real a2 real, parameter :: aa = 0.6931471806E+00 real b real c real, parameter :: c1 = 0.01E+00 real, parameter :: c2 = 0.222222E+00 real, parameter :: c3 = 0.32E+00 real, parameter :: c4 = 0.4E+00 real, parameter :: c5 = 1.24E+00 real, parameter :: c6 = 2.2E+00 real, parameter :: c7 = 4.67E+00 real, parameter :: c8 = 6.66E+00 real, parameter :: c9 = 6.73E+00 real, parameter :: c10 = 13.32E+00 real, parameter :: c11 = 60.0E+00 real, parameter :: c12 = 70.0E+00 real, parameter :: c13 = 84.0E+00 real, parameter :: c14 = 105.0E+00 real, parameter :: c15 = 120.0E+00 real, parameter :: c16 = 127.0E+00 real, parameter :: c17 = 140.0E+00 real, parameter :: c18 = 175.0E+00 real, parameter :: c19 = 210.0E+00 real, parameter :: c20 = 252.0E+00 real, parameter :: c21 = 264.0E+00 real, parameter :: c22 = 294.0E+00 real, parameter :: c23 = 346.0E+00 real, parameter :: c24 = 420.0E+00 real, parameter :: c25 = 462.0E+00 real, parameter :: c26 = 606.0E+00 real, parameter :: c27 = 672.0E+00 real, parameter :: c28 = 707.0E+00 real, parameter :: c29 = 735.0E+00 real, parameter :: c30 = 889.0E+00 real, parameter :: c31 = 932.0E+00 real, parameter :: c32 = 966.0E+00 real, parameter :: c33 = 1141.0E+00 real, parameter :: c34 = 1182.0E+00 real, parameter :: c35 = 1278.0E+00 real, parameter :: c36 = 1740.0E+00 real, parameter :: c37 = 2520.0E+00 real, parameter :: c38 = 5040.0E+00 real cdf real, parameter :: cdf_max = 0.999998E+00 real, parameter :: cdf_min = 0.000002E+00 real ch real, parameter :: e = 0.0000005E+00 real g real gamma_inc real gamma_log integer i integer, parameter :: it_max = 20 real p1 real p2 real q real s1 real s2 real s3 real s4 real s5 real s6 real t real x real x2 real xx ! if ( cdf < cdf_min ) then x = - 1.0E+00 write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_CENTRAL_CDF_INV - Fatal error!' write ( *, * ) ' CDF < CDF_MIN.' stop end if if ( cdf > cdf_max ) then x = - 1.0E+00 write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_CENTRAL_CDF_INV - Fatal error!' write ( *, * ) ' CDF > CDF_MAX.' stop end if xx = 0.5E+00 * a c = xx - 1.0E+00 ! ! Compute Log ( Gamma ( A/2 ) ). ! g = gamma_log ( a / 2.0E+00 ) ! ! Starting approximation for small chi-squared. ! if ( a < - c5 * log ( cdf ) ) then ch = ( cdf * xx * exp ( g + xx * aa ) )**( 1.0E+00 / xx ) if ( ch < e ) then x = ch return end if ! ! Starting approximation for A less than or equal to 0.32. ! else if ( a <= c3 ) then ch = c4 a2 = log ( 1.0E+00 - cdf ) do q = ch p1 = 1.0E+00 + ch * ( c7 + ch ) p2 = ch * ( c9 + ch * ( c8 + ch ) ) t = - 0.5E+00 + ( c7 + 2.0E+00 * ch ) / p1 & - ( c9 + ch * ( c10 + 3.0E+00 * ch ) ) / p2 ch = ch - ( 1.0E+00 - exp ( a2 + g + 0.5E+00 * ch + c * aa ) & * p2 / p1 ) / t if ( abs ( q / ch - 1.0E+00 ) <= c1 ) then exit end if end do ! ! Call to algorithm AS 111. ! Note that P has been tested above. ! AS 241 could be used as an alternative. ! else call normal_01_cdf_inv ( cdf, x2 ) ! ! Starting approximation using Wilson and Hilferty estimate. ! p1 = c2 / a ch = a * ( x2 * sqrt ( p1 ) + 1.0E+00 - p1 )**3 ! ! Starting approximation for P tending to 1. ! if ( ch > c6 * a + 6.0E+00 ) then ch = - 2.0E+00 * ( log ( 1.0E+00 - cdf ) - c * log ( 0.5E+00 * ch ) + g ) end if end if ! ! Call to algorithm AS 239 and calculation of seven term Taylor series. ! do i = 1, it_max q = ch p1 = 0.5E+00 * ch ! ! I'm guessing these GAMMA_INC arguments are in the correct order...? ! p2 = cdf - gamma_inc ( p1, xx ) t = p2 * exp ( xx * aa + g + p1 - c * log ( ch ) ) b = t / ch a2 = 0.5E+00 * t - b * c s1 = ( c19 + a2 & * ( c17 + a2 & * ( c14 + a2 & * ( c13 + a2 & * ( c12 + a2 & * c11 ) ) ) ) ) / c24 s2 = ( c24 + a2 & * ( c29 + a2 & * ( c32 + a2 & * ( c33 + a2 & * c35 ) ) ) ) / c37 s3 = ( c19 + a2 & * ( c25 + a2 & * ( c28 + a2 & * c31 ) ) ) / c37 s4 = ( c20 + a2 & * ( c27 + a2 & * c34 ) + c & * ( c22 + a2 & * ( c30 + a2 & * c36 ) ) ) / c38 s5 = ( c13 + c21 * a2 + c * ( c18 + c26 * a2 ) ) / c37 s6 = ( c15 + c * ( c23 + c16 * c ) ) / c38 ch = ch + t * ( 1.0E+00 + 0.5E+00 * t * s1 - b * c & * ( s1 - b & * ( s2 - b & * ( s3 - b & * ( s4 - b & * ( s5 - b & * s6 ) ) ) ) ) ) if ( abs ( q / ch - 1.0E+00 ) > e ) then x = ch return end if end do x = ch write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_CENTRAL_CDF_INV - Warning!' write ( *, * ) ' Convergence not reached.' return end subroutine chisquare_central_cdf_values ( n, a, x, fx ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_CDF_VALUES returns values of the Chi-Square central CDF. ! ! ! Modified: ! ! 29 May 2001 ! ! Reference: ! ! Milton Abramowitz and Irene Stegun, ! Handbook of Mathematical Functions, ! US Department of Commerce, 1964. ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, integer N. ! On input, if N is 0, the first test data is returned, and N is set ! to the index of the test data. On each subsequent call, N is ! incremented and that test data is returned. When there is no more ! test data, N is set to 0. ! ! Output, integer A, real X, the arguments of the function. ! ! Output, real FX, the value of the function. ! implicit none ! integer, parameter :: nmax = 18 ! integer a integer, save, dimension ( nmax ) :: avec = (/ & 1, 2, 1, 2, & 1, 2, 3, 4, & 1, 2, 3, 4, & 5, 3, 3, 3, & 3, 3 /) real fx real, save, dimension ( nmax ) :: fxvec = (/ & 0.92034E+00, 0.99501E+00, 0.88754E+00, 0.99005E+00, & 0.52709E+00, 0.81873E+00, 0.94024E+00, 0.98248E+00, & 0.31731E+00, 0.60653E+00, 0.80125E+00, 0.90980E+00, & 0.96257E+00, 0.57241E+00, 0.39163E+00, 0.26146E+00, & 0.17180E+00, 0.11161E+00 /) integer n real x real, save, dimension ( nmax ) :: xvec = (/ & 0.01E+00, 0.01E+00, 0.02E+00, 0.02E+00, & 0.40E+00, 0.40E+00, 0.40E+00, 0.40E+00, & 1.00E+00, 1.00E+00, 1.00E+00, 1.00E+00, & 1.00E+00, 2.00E+00, 3.00E+00, 4.00E+00, & 5.00E+00, 6.00E+00 /) ! if ( n < 0 ) then n = 0 end if n = n + 1 if ( n > nmax ) then n = 0 a = 0 x = 0.0E+00 fx = 0.0E+00 return end if a = avec(n) x = xvec(n) fx = fxvec(n) return end subroutine chisquare_central_check ( a ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_CHECK checks the parameter of the central Chi-Squared PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the distribution. ! 1 <= A. ! implicit none ! real a ! if ( a < 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_CENTRAL_CHECK - Fatal error!' write ( *, * ) ' A < 1.0.' stop end if return end subroutine chisquare_central_mean ( a, mean ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_MEAN returns the mean of the central Chi-Squared PDF. ! ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the distribution. ! 1 <= A. ! ! Output, real MEAN, the mean value. ! implicit none ! real a real mean ! mean = a return end subroutine chisquare_central_pdf ( x, a, pdf ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_PDF evaluates the central Chi-Squared PDF. ! ! ! Formula: ! ! PDF(X)(A) = ! EXP ( - X / 2 ) * X**((A-2)/2) / ( 2**(A/2) * GAMMA ( A/2 ) ) ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! 0.0E+00 <= X ! ! Input, real A, the parameter of the PDF. ! 1 <= A. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real gamma real pdf real x ! if ( x < 0.0E+00 ) then pdf = 0.0E+00 else b = a / 2.0E+00 pdf = exp ( - 0.5E+00 * x ) * x**( b - 1.0E+00 ) & / ( 2.0E+00**b * gamma ( b ) ) end if return end subroutine chisquare_central_sample ( a, x ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_SAMPLE samples the central Chisquare PDF. ! ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the PDF. ! 1 <= A. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real a2 real b2 real c2 integer i integer, parameter :: it_max = 100 integer n real x real x2 ! n = int ( a ) if ( real ( n ) == a .and. n <= it_max ) then x = 0.0E+00 do i = 1, n call normal_01_sample ( x2 ) x = x + x2**2 end do else a2 = 0.0E+00 b2 = 1.0E+00 c2 = a / 2.0E+00 call gamma_sample ( a2, b2, c2, x ) x = 2.0E+00 * x end if return end subroutine chisquare_central_variance ( a, variance ) ! !******************************************************************************* ! !! CHISQUARE_CENTRAL_VARIANCE returns the variance of the central Chisquare PDF. ! ! ! Modified: ! ! 10 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the distribution. ! 1 <= A. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real variance ! variance = 2.0E+00 * a return end subroutine chisquare_noncentral_check ( a, b ) ! !******************************************************************************* ! !! CHISQUARE_NONCENTRAL_CHECK checks the parameters of the noncentral Chi-Squared PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the parameter of the PDF. ! 1.0E+00 <= A. ! ! Input, real B, the noncentrality parameter of the PDF. ! 0.0E+00 <= B. ! implicit none ! real a real b ! if ( a < 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_NONCENTRAL_CHECK - Fatal error!' write ( *, * ) ' A < 1.' stop end if if ( b < 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CHISQUARE_NONCENTRAL_CHECK - Fatal error!' write ( *, * ) ' B < 0.' stop end if return end subroutine chisquare_noncentral_mean ( a, b, mean ) ! !******************************************************************************* ! !! CHISQUARE_NONCENTRAL_MEAN returns the mean of the noncentral Chi-Squared PDF. ! ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the parameter of the PDF. ! 1.0E+00 <= A. ! ! Input, real B, the noncentrality parameter of the PDF. ! 0.0E+00 <= B. ! ! Output, real MEAN, the mean value. ! implicit none ! real a real b real mean ! mean = a + b return end subroutine chisquare_noncentral_sample ( a, b, x ) ! !******************************************************************************* ! !! CHISQUARE_NONCENTRAL_SAMPLE samples the noncentral Chisquare PDF. ! ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the parameter of the PDF. ! 1.0E+00 <= A. ! ! Input, real B, the noncentrality parameter of the PDF. ! 0.0E+00 <= B. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real a1 real a2 real b real b2 real x real x1 real x2 ! a1 = a - 1.0E+00 call chisquare_central_sample ( a1, x1 ) a2 = sqrt ( b ) b2 = 1.0E+00 call normal_sample ( a2, b2, x2 ) x = x1 + x2**2 return end subroutine chisquare_noncentral_variance ( a, b, variance ) ! !******************************************************************************* ! !! CHISQUARE_NONCENTRAL_VARIANCE returns the variance of the noncentral Chi-Squared PDF. ! ! ! Modified: ! ! 30 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, the parameter of the PDF. ! 1 <= A. ! ! Input, real B, the noncentrality parameter of the PDF. ! 0.0E+00 <= B. ! ! Output, real VARIANCE, the variance value. ! implicit none ! real a real b real variance ! variance = 2.0E+00 * ( a + 2.0E+00 * b ) return end subroutine circle_sample ( a, b, c, x1, x2 ) ! !******************************************************************************* ! !! CIRCLE_SAMPLE samples points from a circle. ! ! ! Modified: ! ! 03 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, C, the parameters of the circle. ! The circle is centered at (A,B) and has radius C. ! 0.0E+00 < C. ! ! Output, real X1, X2, a sampled point of the circle. ! implicit none ! real a real angle real b real c real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real radius_frac real x1 real x2 ! call r_random ( 0.0E+00, 1.0E+00, radius_frac ) radius_frac = sqrt ( radius_frac ) call r_random ( 0.0E+00, 2.0E+00 * pi, angle ) x1 = a + c * radius_frac * cos ( angle ) x2 = b + c * radius_frac * sin ( angle ) return end subroutine circular_normal_01_mean ( mean ) ! !******************************************************************************* ! !! CIRCULAR_NORMAL_01_MEAN returns the mean of the Circular Normal 01 PDF. ! ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real MEAN(2), the mean of the PDF. ! implicit none ! real mean(2) ! mean(1) = 0.0E+00 mean(2) = 0.0E+00 return end subroutine circular_normal_01_pdf ( x, pdf ) ! !******************************************************************************* ! !! CIRCULAR_NORMAL_01_PDF evaluates the Circular Normal 01 PDF. ! ! ! Formula: ! ! PDF(X) = EXP ( - 0.5E+00 * ( X(1)**2 + X(2)**2 ) ) / ( 2 * PI ) ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X(2), the argument of the PDF. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x(2) ! pdf = exp ( - 0.5E+00 * ( x(1)**2 + x(2)**2 ) ) / ( 2.0E+00 * pi ) return end subroutine circular_normal_01_sample ( x ) ! !******************************************************************************* ! !! CIRCULAR_NORMAL_01_SAMPLE samples the Circular Normal 01 PDF. ! ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real X(2), a sample of the PDF. ! implicit none ! real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real v1 real v2 real x(2) ! call r_random ( 0.0E+00, 1.0E+00, v1 ) call r_random ( 0.0E+00, 1.0E+00, v2 ) x(1) = sqrt ( - 2.0E+00 * log ( v1 ) ) * cos ( 2.0E+00 * pi * v2 ) x(2) = sqrt ( - 2.0E+00 * log ( v1 ) ) * sin ( 2.0E+00 * pi * v2 ) return end subroutine circular_normal_01_variance ( variance ) ! !******************************************************************************* ! !! CIRCULAR_NORMAL_01_VARIANCE returns the variance of the Circular Normal 01 PDF. ! ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real VARIANCE(2), the variance of the PDF. ! implicit none ! real variance(2) ! variance(1) = 1.0E+00 variance(2) = 1.0E+00 return end subroutine cosine_cdf ( x, a, b, cdf ) ! !******************************************************************************* ! !! COSINE_CDF evaluates the Cosine CDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Input, real A, B, the parameter of the PDF. ! 0.0E+00 < B. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x real y ! if ( x <= a - pi * b ) then cdf = 0.0E+00 else if ( x <= a + pi * b ) then y = ( x - a ) / b cdf = 0.5E+00 + ( y + sin ( y ) ) / ( 2.0E+00 * pi ) else if ( x > a + pi * b ) then cdf = 1.0E+00 end if return end subroutine cosine_cdf_inv ( cdf, a, b, x ) ! !******************************************************************************* ! !! COSINE_CDF_INV inverts the Cosine CDF. ! ! ! Discussion: ! ! A simple bisection method is used on the interval ! [ A - PI * B, A + PI * B ]. ! ! Modified: ! ! 30 December 1999 ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real X, the corresponding argument of the CDF. ! implicit none ! real a real b real cdf real cdf1 real cdf2 real cdf3 integer it integer, parameter :: it_max = 100 real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real, parameter :: tol = 0.0001E+00 real x real x1 real x2 real x3 ! if ( cdf <= 0.0E+00 ) then x = a - pi * b return else if ( cdf >= 1.0E+00 ) then x = a + pi * b return end if ! x1 = a - pi * b cdf1 = 0.0E+00 x2 = a + pi * b cdf2 = 1.0E+00 ! ! Now use bisection. ! it = 0 do it = 1, it_max x3 = 0.5E+00 * ( x1 + x2 ) call cosine_cdf ( x3, a, b, cdf3 ) if ( abs ( cdf3 - cdf ) < tol ) then x = x3 return end if if ( sign ( 1.0E+00, cdf3 - cdf ) == sign ( 1.0E+00, cdf1 - cdf ) ) then x1 = x3 cdf1 = cdf3 else x2 = x3 cdf2 = cdf3 end if end do write ( *, * ) ' ' write ( *, * ) 'COSINE_CDF_INV - Fatal error!' write ( *, * ) ' Iteration limit exceeded.' stop end subroutine cosine_check ( a, b ) ! !******************************************************************************* ! !! COSINE_CHECK checks the parameters of the Cosine CDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameter of the PDF. ! 0.0E+00 < B. ! implicit none ! real a real b ! if ( b <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'COSINE_CHECK - Fatal error!' write ( *, * ) ' B <= 0.0' stop end if return end subroutine cosine_mean ( a, b, mean ) ! !******************************************************************************* ! !! COSINE_MEAN returns the mean of the Cosine PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! real a real b real mean ! mean = a return end subroutine cosine_pdf ( x, a, b, pdf ) ! !******************************************************************************* ! !! COSINE_PDF evaluates the Cosine PDF. ! ! ! Formula: ! ! PDF(X)(A,B) = ( 1 / ( 2 * PI * B ) ) * COS ( ( X - A ) / B ) ! for A - PI * B <= X <= A + PI * B ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x real y ! if ( x < a - pi * b ) then pdf = 0.0E+00 else if ( x <= a + pi * b ) then y = ( x - a ) / b pdf = 1.0E+00 / ( 2.0E+00 * pi * b ) * cos ( y ) else if ( x > a + pi * b ) then pdf = 0.0E+00 end if return end subroutine cosine_sample ( a, b, x ) ! !******************************************************************************* ! !! COSINE_SAMPLE samples the Cosine PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real b real cdf real x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call cosine_cdf_inv ( cdf, a, b, x ) return end subroutine cosine_variance ( a, b, variance ) ! !******************************************************************************* ! !! COSINE_VARIANCE returns the variance of the Cosine PDF. ! ! ! Modified: ! ! 30 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! 0.0E+00 < B. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! real a real b real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real variance ! variance = ( pi**2 / 3.0E+00 - 2.0E+00 ) * b**2 return end subroutine coupon_mean ( j, n, mean ) ! !******************************************************************************* ! !! COUPON_MEAN returns the variance of the Coupon PDF. ! ! ! Modified: ! ! 14 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer J, the number of distinct coupons to be collected. ! J must be between 1 and N. ! ! Input, integer N, the number of distinct coupons. ! ! Output, real MEAN, the mean number of coupons that must be collected ! in order to just get J distinct kinds. ! implicit none ! integer i integer j integer n real mean ! if ( j > n ) then write ( *, * ) ' ' write ( *, * ) 'COUPON_MEAN - Fatal error!' write ( *, * ) ' Number of distinct coupons desired must be no more' write ( *, * ) ' than the total number of distinct coupons.' stop end if mean = 0.0E+00 do i = 1, j mean = mean + 1.0E+00 / real ( n - i + 1 ) end do mean = mean * real ( n ) return end subroutine coupon_simulate ( n_type, coupon, n_coupon ) ! !******************************************************************************* ! !! COUPON_SIMULATE simulates the coupon collector's problem. ! ! ! Discussion: ! ! The coupon collector needs to collect one of each of N_TYPE ! coupons. The collector may draw one coupon on each trial, ! and takes as many trials as necessary to complete the task. ! On each trial, the probability of picking any particular type ! of coupon is always 1 / N_TYPE. ! ! The most interesting question is, what is the expected number ! of drawings necessary to complete the collection? ! how does this number vary as N_TYPE increases? What is the ! distribution of the numbers of each type of coupon in a typical ! collection when it is just completed? ! ! As N increases, the number of coupons necessary to be ! collected in order to get a complete set in any simulation ! strongly tends to the value N_TYPE * LOG ( N_TYPE ). ! ! If N_TYPE is 1, the simulation ends with a single drawing. ! ! If N_TYPE is 2, then we may call the coupon taken on the first drawing ! a "Head", say, and the process then is similar to the question of the ! length, plus one, of a run of Heads or Tails in coin flipping. ! ! Modified: ! ! 31 March 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N_TYPE, the number of types of coupons. ! ! Output, integer COUPON(N_TYPE), the number of coupons of each type ! that were collected during the simulation. ! ! Output, integer N_COUPON, the total number of coupons collected. ! implicit none ! integer n_type ! integer coupon(n_type) integer i integer, parameter :: max_coupon = 2000 integer n_coupon integer straight ! coupon(1:n_type) = 0 straight = 0 n_coupon = 0 ! ! Draw another coupon. ! do while ( n_coupon < max_coupon ) call i_random ( 1, n_type, i ) ! ! Increment the number of I coupons. ! coupon(i) = coupon(i) + 1 n_coupon = n_coupon + 1 ! ! If I is the next one we needed, increase STRAIGHT by 1. ! if ( i == straight + 1 ) then do straight = straight + 1 ! ! If STRAIGHT = N_TYPE, we have all of them. ! if ( straight >= n_type ) then return end if ! ! If the next coupon has not been collected, our straight is over. ! if ( coupon(straight+1) <= 0 ) then exit end if end do end if end do write ( *, * ) ' ' write ( *, * ) 'COUPON_SIMULATE - Fatal error!' write ( *, * ) ' Maximum number of coupons drawn without success.' stop end subroutine coupon_variance ( j, n, variance ) ! !******************************************************************************* ! !! COUPON_VARIANCE returns the variance of the Coupon PDF. ! ! ! Modified: ! ! 14 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer J, the number of distinct coupons to be collected. ! J must be between 1 and N. ! ! Input, integer N, the number of distinct coupons. ! ! Output, real VARIANCE, the variance of the number of ! coupons that must be collected in order to just get J distinct kinds. ! implicit none ! integer i integer j integer n real variance ! if ( j > n ) then write ( *, * ) ' ' write ( *, * ) 'COUPON_VARIANCE - Fatal error!' write ( *, * ) ' Number of distinct coupons desired must be no more' write ( *, * ) ' than the total number of distinct coupons.' stop end if variance = 0.0E+00 do i = 1, j variance = variance + real ( i - 1 ) / real ( n - i + 1 )**2 end do variance = variance * real ( n ) return end function csc ( theta ) ! !******************************************************************************* ! !! CSC returns the cosecant of X. ! ! ! Definition: ! ! CSC ( THETA ) = 1.0E+00 / SIN ( THETA ) ! ! Discussion: ! ! CSC is not a built-in function in FORTRAN, and occasionally it ! is handier, or more concise, to be able to refer to it directly ! rather than through its definition in terms of the sine function. ! ! Modified: ! ! 01 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real THETA, the angle, in radians, whose cosecant is desired. ! It must be the case that SIN ( THETA ) is not zero. ! ! Output, real CSC, the cosecant of THETA. ! implicit none ! real csc real theta ! csc = sin ( theta ) if ( csc == 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'CSC - Fatal error!' write ( *, * ) ' CSC undefined for THETA = ', theta stop end if csc = 1.0E+00 / csc return end subroutine deranged_cdf ( x, a, cdf ) ! !******************************************************************************* ! !! DERANGED_CDF evaluates the Deranged CDF. ! ! ! Modified: ! ! 06 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the maximum number of items in their correct places. ! 0 <= X <= A. ! ! Input, integer A, the number of items. ! 1 <= A. ! ! Output, real CDF, the value of the CDF. ! implicit none ! integer a real cdf integer cnk integer deranged_enum integer dnmk integer i_factorial integer nfact integer sum2 integer x integer x2 ! if ( x < 0 .or. x > a ) then cdf = 0.0E+00 else sum2 = 0 do x2 = 0, x call binomial_coef ( a, x2, cnk ) dnmk = deranged_enum ( a-x2 ) sum2 = sum2 + cnk * dnmk end do nfact = i_factorial ( a ) cdf = real ( sum2 ) / real ( nfact ) end if return end subroutine deranged_cdf_inv ( cdf, a, x ) ! !******************************************************************************* ! !! DERANGED_CDF_INV inverts the Deranged CDF. ! ! ! Modified: ! ! 06 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! 0.0E+00 <= CDF <= 1.0. ! ! Input, integer A, the number of items. ! 1 <= A. ! ! Output, integer X, the corresponding argument. ! implicit none ! integer a real cdf real cdf2 real pdf integer x integer x2 ! if ( cdf < 0.0E+00 .or. cdf > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DERANGED_CDF_INV - Fatal error!' write ( *, * ) ' CDF < 0 or 1 < CDF.' stop end if cdf2 = 0.0E+00 do x2 = 0, a call deranged_pdf ( x2, a, pdf ) cdf2 = cdf2 + pdf if ( cdf <= cdf2 ) then x = x2 return end if end do x = a return end subroutine deranged_check ( a ) ! !******************************************************************************* ! !! DERANGED_CHECK checks the parameter of the Deranged PDF. ! ! ! Modified: ! ! 18 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the total number of items. ! 1 <= A. ! implicit none ! integer a ! if ( a < 1 ) then write ( *, * ) ' ' write ( *, * ) 'DERANGED_CHECK - Fatal error!' write ( *, * ) ' A < 1.' stop end if return end function deranged_enum ( n ) ! !******************************************************************************* ! !! DERANGED_ENUM returns the number of derangements of N objects. ! ! ! Definition: ! ! A derangement of N objects is a permutation with no fixed ! points. If we symbolize the permutation operation by "P", ! then for a derangment, P(I) is never equal to I. ! ! Recursion: ! ! D(0) = 1 ! D(1) = 0 ! D(2) = 1 ! D(N) = (N-1) * ( D(N-1) + D(N-2) ) ! ! or ! ! D(0) = 1 ! D(1) = 0 ! D(N) = N * D(N-1) + (-1)**N ! ! Formula: ! ! D(N) = N! * ( 1 - 1/1! + 1/2! - 1/3! ... 1/N! ) ! ! Based on the inclusion/exclusion law. ! ! Comments: ! ! D(N) is the number of ways of placing N non-attacking rooks on ! an N by N chessboard with one diagonal deleted. ! ! Limit ( N -> Infinity ) D(N)/N! = 1 / e. ! ! The number of permutations with exactly K items in the right ! place is COMB(N,K) * D(N-K). ! ! First values: ! ! N D(N) ! 0 1 ! 1 0 ! 2 1 ! 3 2 ! 4 9 ! 5 44 ! 6 265 ! 7 1854 ! 8 14833 ! 9 133496 ! 10 1334961 ! ! Modified: ! ! 01 February 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of objects to be permuted. ! ! Output, integer DERANGED_ENUM, the number of derangements of N objects. ! implicit none ! integer deranged_enum integer dn integer dnm1 integer dnm2 integer i integer n ! if ( n < 0 ) then dn = 0 else if ( n == 0 ) then dn = 1 else if ( n == 1 ) then dn = 0 else if ( n == 2 ) then dn = 1 else dnm1 = 0 dn = 1 do i = 3, n dnm2 = dnm1 dnm1 = dn dn = ( i - 1 ) * ( dnm1 + dnm2 ) end do end if deranged_enum = dn return end subroutine deranged_mean ( a, mean ) ! !******************************************************************************* ! !! DERANGED_MEAN returns the mean of the Deranged CDF. ! ! ! Discussion: ! ! The mean is computed by straightforward summation. ! ! Modified: ! ! 05 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of items. ! 1 <= A. ! ! Output, real MEAN, the mean of the PDF. ! implicit none ! integer a real mean real pdf integer x ! mean = 0.0E+00 do x = 1, a call deranged_pdf ( x, a, pdf ) mean = mean + pdf * x end do return end subroutine deranged_pdf ( x, a, pdf ) ! !******************************************************************************* ! !! DERANGED_PDF evaluates the Deranged PDF. ! ! ! Definition: ! ! PDF(X)(A) is the probability that exactly X items will occur in ! their proper place after a random permutation of A items. ! ! Modified: ! ! 06 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X, the number of items in their correct places. ! 0 <= X <= A. ! ! Input, integer A, the total number of items. ! 1 <= A. ! ! Output, real PDF, the value of the PDF. ! implicit none ! integer a integer cnk integer deranged_enum integer dnmk integer i_factorial integer nfact real pdf integer x ! if ( x < 0 .or. x > a ) then pdf = 0.0E+00 else call binomial_coef ( a, x, cnk ) dnmk = deranged_enum ( a-x ) nfact = i_factorial ( a ) pdf = real ( cnk * dnmk ) / real ( nfact ) end if return end subroutine deranged_sample ( a, x ) ! !******************************************************************************* ! !! DERANGED_SAMPLE samples the Deranged PDF. ! ! ! Modified: ! ! 06 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of items. ! 1 <= A. ! ! Output, integer X, a sample of the PDF. ! implicit none ! integer a real cdf integer x ! call r_random ( 0.0E+00, 1.0E+00, cdf ) call deranged_cdf_inv ( cdf, a, x ) return end subroutine deranged_variance ( a, variance ) ! !******************************************************************************* ! !! DERANGED_VARIANCE returns the variance of the Deranged CDF. ! ! ! Discussion: ! ! The variance is computed by straightforward summation. ! ! Modified: ! ! 06 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer A, the number of items. ! 1 <= A. ! ! Output, real VARIANCE, the variance of the PDF. ! implicit none ! integer a real mean real pdf integer x real variance ! call deranged_mean ( a, mean ) variance = 0.0E+00 do x = 1, a call deranged_pdf ( x, a, pdf ) variance = variance + pdf * ( x - mean )**2 end do return end function digamma ( x ) ! !******************************************************************************* ! !! DIGAMMA calculates the digamma or Psi function = d ( LOG ( GAMMA ( X ) ) ) / dX ! ! ! Reference: ! ! J Bernardo, ! Psi ( Digamma ) Function, ! Algorithm AS 103, ! Applied Statistics, ! Volume 25, Number 3, pages 315-317, 1976. ! ! Modified: ! ! 03 January 2000 ! ! Parameters: ! ! Input, real X, the argument of the digamma function. ! 0 < X. ! ! Output, real DIGAMMA, the value of the digamma function at X. ! implicit none ! real, parameter :: c = 8.5E+00 real, parameter :: d1 = -0.5772156649E+00 real digamma real r real, parameter :: s = 0.00001E+00 real, parameter :: s3 = 0.08333333333E+00 real, parameter :: s4 = 0.0083333333333E+00 real, parameter :: s5 = 0.003968253968E+00 real x real y ! ! The argument must be positive. ! if ( x <= 0.0E+00 ) then digamma = 0.0E+00 write ( *, * ) ' ' write ( *, * ) 'DIGAMMA - Fatal error!' write ( *, * ) ' X <= 0.' stop ! ! Use approximation if argument <= S. ! else if ( x <= s ) then digamma = d1 - 1.0E+00 / x ! ! Reduce the argument to DIGAMMA(X + N) where (X + N) >= C. ! else digamma = 0.0E+00 y = x do while ( y < c ) digamma = digamma - 1.0E+00 / y y = y + 1.0E+00 end do ! ! Use Stirling's (actually de Moivre's) expansion if argument > C. ! r = 1.0E+00 / y**2 digamma = digamma + log ( y ) - 0.5E+00 / y & - r * ( s3 - r * ( s4 - r * s5 ) ) end if return end subroutine dipole_cdf ( x, a, b, cdf ) ! !******************************************************************************* ! !! DIPOLE_CDF evaluates the Dipole CDF. ! ! ! Modified: ! ! 28 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the CDF. ! ! Input, real A, B, the parameters of the PDF. ! A is arbitrary, but represents an angle, so only 0 <= A <= 2 * PI ! is interesting, and -1.0E+00 <= B <= 1.0. ! ! Output, real CDF, the value of the CDF. ! implicit none ! real a real b real cdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! cdf = 0.5E+00 + ( 1.0E+00 / pi ) * atan ( x ) & + b**2 * ( x * cos ( 2.0E+00 * a ) & - sin ( 2.0E+00 * a ) ) / ( pi * ( 1.0E+00 + x**2 ) ) return end subroutine dipole_cdf_inv ( cdf, a, b, x ) ! !******************************************************************************* ! !! DIPOLE_CDF_INV inverts the Dipole CDF. ! ! ! Discussion: ! ! A simple bisection method is used. ! ! Modified: ! ! 04 January 2000 ! ! Parameters: ! ! Input, real CDF, the value of the CDF. ! ! Input, real A, B, the parameters of the PDF. ! -1.0E+00 <= B <= 1.0. ! ! Output, real X, the corresponding argument of the CDF. ! implicit none ! real a real b real cdf real cdf1 real cdf2 real cdf3 integer it integer, parameter :: it_max = 100 real, parameter :: tol = 0.0001E+00 real x real x1 real x2 real x3 ! ! Take care of horrible input. ! if ( cdf <= 0.0E+00 ) then x = - huge ( x ) return else if ( cdf >= 1.0E+00 ) then x = huge ( x ) return end if ! ! Seek X1 < X < X2. ! x1 = - 1.0E+00 do call dipole_cdf ( x1, a, b, cdf1 ) if ( cdf1 <= cdf ) then exit end if x1 = 2.0E+00 * x1 end do x2 = 1.0E+00 do call dipole_cdf ( x2, a, b, cdf2 ) if ( cdf2 >= cdf ) then exit end if x2 = 2.0E+00 * x2 end do ! ! Now use bisection. ! it = 0 do it = it + 1 x3 = 0.5E+00 * ( x1 + x2 ) call dipole_cdf ( x3, a, b, cdf3 ) if ( abs ( cdf3 - cdf ) < tol ) then x = x3 exit end if if ( it > it_max ) then write ( *, * ) ' ' write ( *, * ) 'DIPOLE_CDF_INV - Fatal error!' write ( *, * ) ' Iteration limit exceeded.' stop end if if ( sign ( 1.0E+00, cdf3 - cdf ) == sign ( 1.0E+00, cdf1 - cdf ) ) then x1 = x3 cdf1 = cdf3 else x2 = x3 cdf2 = cdf3 end if end do return end subroutine dipole_check ( a, b ) ! !******************************************************************************* ! !! DIPOLE_CHECK checks the parameters of the Dipole CDF. ! ! ! Modified: ! ! 28 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! A is arbitrary, but represents an angle, so only 0 <= A <= 2 * PI ! is interesting, and -1.0E+00 <= B <= 1.0. ! implicit none ! real a real b ! if ( b < -1.0E+00 .or. b > 1.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DIPOLE_CHECK - Fatal error!' write ( *, * ) ' -1.0E+00 <= B <= 1.0E+00 is required.' write ( *, * ) ' The input B = ', b stop end if return end subroutine dipole_pdf ( x, a, b, pdf ) ! !******************************************************************************* ! !! DIPOLE_PDF evaluates the Dipole PDF. ! ! ! Formula: ! ! PDF(X)(A,B) = ! 1 / ( PI * ( 1 + X**2 ) ) ! + B**2 * ( ( 1 - X**2 ) * cos ( 2 * A ) + 2.0E+00 * X * sin ( 2 * A ) ) ! / ( PI * ( 1 + X )**2 ) ! ! Discussion: ! ! Densities of this kind commonly occur in the analysis of resonant ! scattering of elementary particles. ! ! DIPOLE_PDF(X)(A,0) = CAUCHY_PDF(X)(A) ! A = 0, B = 1 yields the single channel dipole distribution. ! ! Reference: ! ! Robert Knop, ! Algorithm 441, ! ACM Transactions on Mathematical Software. ! ! Modified: ! ! 28 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X, the argument of the PDF. ! ! Input, real A, B, the parameters of the PDF. ! A is arbitrary, but represents an angle, so only 0 <= A <= 2 * PI ! is interesting, ! and -1.0E+00 <= B <= 1.0. ! ! Output, real PDF, the value of the PDF. ! implicit none ! real a real b real pdf real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real x ! pdf = 1.0E+00 / ( pi * ( 1.0E+00 + x**2 ) ) & + b**2 * ( ( 1.0E+00 - x**2 ) * cos ( 2.0E+00 * a ) & + 2.0E+00 * x * sin ( 2.0E+00 * x ) ) / ( pi * ( 1.0E+00 + x**2 )**2 ) return end subroutine dipole_sample ( a, b, x ) ! !******************************************************************************* ! !! DIPOLE_SAMPLE samples the Dipole PDF. ! ! ! Reference: ! ! Robert Knop, ! Algorithm 441, ! ACM Transactions on Mathematical Software. ! ! Modified: ! ! 04 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A, B, the parameters of the PDF. ! A is arbitrary, but represents an angle, so only 0 <= A <= 2 * PI ! is interesting, ! and -1.0E+00 <= B <= 1.0E+00. ! ! Output, real X, a sample of the PDF. ! implicit none ! real a real a2 real b real b2 real c2 real x real x1 real x2 ! ! Find (X1,X2) at random in a circle. ! a2 = b * sin ( a ) b2 = b * cos ( a ) c2 = 1.0E+00 call circle_sample ( a2, b2, c2, x1, x2 ) ! ! The dipole variate is the ratio X1 / X2. ! x = x1 / x2 return end subroutine dirichlet_check ( n, a ) ! !******************************************************************************* ! !! DIRICHLET_CHECK checks the parameters of the Dirichlet PDF. ! ! ! Modified: ! ! 08 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components. ! ! Input, real A(N), the probabilities for each component. ! Each A(I) should be nonnegative, and at least one should be positive. ! implicit none ! integer n ! real a(n) integer i logical positive ! positive = .false. do i = 1, n if ( a(i) < 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_CHECK - Fatal error!' write ( *, * ) ' A(I) < 0.' write ( *, * ) ' For I = ', i write ( *, * ) ' A(I) = ', a(i) stop else if ( a(i) > 0.0E+00 ) then positive = .true. end if end do if ( .not. positive ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_CHECK - Fatal error!' write ( *, * ) ' All parameters are zero!' stop end if return end subroutine dirichlet_mean ( n, a, mean ) ! !******************************************************************************* ! !! DIRICHLET_MEAN returns the means of the Dirichlet PDF. ! ! ! Modified: ! ! 23 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components. ! ! Input, real A(N), the probabilities for each component. ! Each A(I) should be nonnegative, and at least one should be positive. ! ! Output, real MEAN(N), the means of the PDF. ! implicit none ! integer n ! real a(n) real mean(n) ! mean(1:n) = a(1:n) call rvec_unit_sum ( n, mean ) return end subroutine dirichlet_mix_check ( comp_num, elem_max, elem_num, a, comp_weight ) ! !******************************************************************************* ! !! DIRICHLET_MIX_CHECK checks the parameters of a Dirichlet mixture PDF. ! ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer COMP_NUM, the number of components in the Dirichlet ! mixture density, that is, the number of distinct Dirichlet PDF's ! that are mixed together. ! ! Input, integer ELEM_MAX, the leading dimension of A, which must ! be at least ELEM_NUM. ! ! Input, integer ELEM_NUM, the number of elements of an observation. ! ! Input, real A(ELEM_MAX,COMP_NUM), the probabilities for element ELEM_NUM ! in component COMP_NUM. ! Each A(I,J) should be greater than or equal to 0.0. ! ! Input, integer COMP_WEIGHT(COMP_NUM), the mixture weights of the densities. ! These do not need to be normalized. The weight of a given component is ! the relative probability that that component will be used to generate ! the sample. ! implicit none ! integer comp_num integer elem_max integer elem_num ! real a(elem_max,comp_num) integer comp_i real comp_weight(comp_num) integer elem_i logical positive ! do comp_i = 1, comp_num do elem_i = 1, elem_num if ( a(elem_i,comp_i) < 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_MIX_CHECK - Fatal error!' write ( *, * ) ' A(ELEM,COMP) < 0.' write ( *, * ) ' COMP = ', comp_i write ( *, * ) ' ELEM = ', elem_i write ( *, * ) ' A(COMP,ELEM) = ', a(elem_i,comp_i) stop end if end do end do positive = .false. do comp_i = 1, comp_num if ( comp_weight(comp_i) < 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_MIX_CHECK - Fatal error!' write ( *, * ) ' COMP_WEIGHT(COMP) < 0.' write ( *, * ) ' COMP = ', comp_i write ( *, * ) ' COMP_WEIGHT(COMP) = ', comp_weight(comp_i) stop else if ( comp_weight(comp_i) > 0.0E+00 ) then positive = .true. end if end do if ( .not. positive ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_MIX_CHECK - Fatal error!' write ( *, * ) ' All component weights are zero.' stop end if return end subroutine dirichlet_mix_mean ( comp_num, elem_max, elem_num, a, comp_weight, & mean ) ! !******************************************************************************* ! !! DIRICHLET_MIX_MEAN returns the means of a Dirichlet mixture PDF. ! ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer COMP_NUM, the number of components in the Dirichlet ! mixture density, that is, the number of distinct Dirichlet PDF's ! that are mixed together. ! ! Input, integer ELEM_MAX, the leading dimension of A, which must ! be at least ELEM_NUM. ! ! Input, integer ELEM_NUM, the number of elements of an observation. ! ! Input, real A(ELEM_MAX,COMP_NUM), the probabilities for element ELEM_NUM ! in component COMP_NUM. ! Each A(I,J) should be greater than or equal to 0.0. ! ! Input, integer COMP_WEIGHT(COMP_NUM), the mixture weights of the densities. ! These do not need to be normalized. The weight of a given component is ! the relative probability that that component will be used to generate ! the sample. ! ! Output, real MEAN(ELEM_NUM), the means for each element. ! implicit none ! integer comp_num integer elem_max integer elem_num ! real a(elem_max,comp_num) integer comp_i real comp_mean(elem_num) real comp_weight(comp_num) real comp_weight_sum integer elem_i real mean(elem_num) ! comp_weight_sum = sum ( comp_weight ) mean(1:elem_num) = 0.0E+00 do comp_i = 1, comp_num call dirichlet_mean ( elem_num, a(1,comp_i), comp_mean ) mean(1:elem_num) = mean(1:elem_num) & + comp_weight(comp_i) * comp_mean(1:elem_num) end do mean(1:elem_num) = mean(1:elem_num) / comp_weight_sum return end subroutine dirichlet_mix_pdf ( x, comp_num, elem_max, elem_num, a, & comp_weight, pdf ) ! !******************************************************************************* ! !! DIRICHLET_MIX_PDF evaluates a Dirichlet mixture PDF. ! ! ! Discussion: ! ! The PDF is a weighted sum of Dirichlet PDF's. Each PDF is a ! "component", with an associated weight. ! ! Modified: ! ! 13 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X(ELEM_NUM), the argument of the PDF. ! ! Input, integer COMP_NUM, the number of components in the Dirichlet ! mixture density, that is, the number of distinct Dirichlet PDF's ! that are mixed together. ! ! Input, integer ELEM_MAX, the leading dimension of A, which must ! be at least ELEM_NUM. ! ! Input, integer ELEM_NUM, the number of elements of an observation. ! ! Input, real A(ELEM_MAX,COMP_NUM), the probabilities for element ELEM_NUM ! in component COMP_NUM. ! Each A(I,J) should be greater than or equal to 0.0. ! ! Input, integer COMP_WEIGHT(COMP_NUM), the mixture weights of the densities. ! These do not need to be normalized. The weight of a given component is ! the relative probability that that component will be used to generate ! the sample. ! ! Output, real PDF, the value of the PDF. ! implicit none ! integer comp_num integer elem_max integer elem_num ! real a(elem_max,comp_num) integer comp_i real comp_pdf real comp_weight(comp_num) real comp_weight_sum real pdf real x(elem_num) ! comp_weight_sum = sum ( comp_weight ) pdf = 0.0E+00 do comp_i = 1, comp_num call dirichlet_pdf ( x, elem_num, a(1,comp_i), comp_pdf ) pdf = pdf + comp_weight(comp_i) * comp_pdf / comp_weight_sum end do return end subroutine dirichlet_mix_sample ( comp_num, elem_max, elem_num, a, & comp_weight, comp, x ) ! !******************************************************************************* ! !! DIRICHLET_MIX_SAMPLE samples a Dirichlet mixture PDF. ! ! ! Modified: ! ! 04 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer COMP_NUM, the number of components in the Dirichlet ! mixture density, that is, the number of distinct Dirichlet PDF's ! that are mixed together. ! ! Input, integer ELEM_MAX, the leading dimension of A, which must ! be at least ELEM_NUM. ! ! Input, integer ELEM_NUM, the number of elements of an observation. ! ! Input, real A(ELEM_MAX,COMP_NUM), the probabilities for element ELEM_NUM ! in component COMP_NUM. ! Each A(I,J) should be greater than or equal to 0.0. ! ! Input, integer COMP_WEIGHT(COMP_NUM), the mixture weights of the densities. ! These do not need to be normalized. The weight of a given component is ! the relative probability that that component will be used to generate ! the sample. ! ! Output, integer COMP, the index of the component of the Dirichlet ! mixture that was chosen to generate the sample. ! ! Output, real X(ELEM_NUM), a sample of the PDF. ! implicit none ! integer comp_num integer elem_max integer elem_num ! real a(elem_max,comp_num) integer comp real comp_weight(comp_num) real x(elem_num) ! ! Choose a particular density component COMP. ! call discrete_sample ( comp_num, comp_weight, comp ) ! ! Sample the density number COMP. ! call dirichlet_sample ( elem_num, a(1,comp), x ) return end subroutine dirichlet_moment2 ( n, a, m2 ) ! !******************************************************************************* ! !! DIRICHLET_MOMENT2 returns the second moments of the Dirichlet PDF. ! ! ! Modified: ! ! 23 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components. ! ! Input, real A(N), the probabilities for each component. ! Each A(I) should be nonnegative, and at least one should be positive. ! ! Output, real M2(N,N), the second moments of the PDF. ! implicit none ! integer n ! real a(n) real a_sum real m2(n,n) integer i integer j ! a_sum = sum ( a(1:n) ) do i = 1, n do j = 1, n if ( i == j ) then m2(i,j) = a(i) * ( a(i) + 1.0E+00 ) / ( a_sum * ( a_sum + 1.0E+00 ) ) else m2(i,j) = a(i) * a(j) / ( a_sum * ( a_sum + 1.0E+00 ) ) end if end do end do return end subroutine dirichlet_multinomial_pdf ( x, a, b, c, pdf ) ! !******************************************************************************* ! !! DIRICHLET_MULTINOMIAL_PDF evaluates a Dirichlet Multinomial PDF. ! ! ! Formula: ! ! PDF(X)(A,B,C) = Comb(A,B,X) * ( Gamma(C_Sum) / Gamma(C_Sum+A) ) ! Product ( 1 <= I <= B ) Gamma(C(I)+X(I)) / Gamma(C(I)) ! ! where: ! ! Comb(A,B,X) is the multinomial coefficient C( A; X(1), X(2), ..., X(B) ), ! C_Sum = Sum ( 1 <= I <= B ) C(I) ! ! Reference: ! ! Kenneth Lange, ! Mathematical and Statistical Methods for Genetic Analysis, ! Springer, 1997, page 45. ! ! Modified: ! ! 17 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer X(B); X(I) counts the number of occurrences of ! outcome I, out of the total of A trials. ! ! Input, integer A, the total number of trials. ! ! Input, integer B, the number of different possible outcomes on ! one trial. ! ! Input, integer C(B); C(I) is the Dirichlet parameter associated ! with outcome I. ! ! Output, real PDF, the value of the Dirichlet multinomial PDF. ! implicit none ! integer b ! integer a real c(b) real c_sum real gamma_log integer i real pdf real pdf_log integer x(b) ! c_sum = sum ( c(1:b) ) pdf_log = - gamma_log ( c_sum + real ( a ) ) + gamma_log ( c_sum ) & + gamma_log ( real ( a + 1 ) ) do i = 1, b pdf_log = pdf_log + gamma_log ( c(i) + real ( x(i) ) ) & - gamma_log ( c(i) ) - gamma_log ( real ( x(i) + 1 ) ) end do pdf = exp ( pdf_log ) return end subroutine dirichlet_pdf ( x, n, a, pdf ) ! !******************************************************************************* ! !! DIRICHLET_PDF evaluates the Dirichlet PDF. ! ! ! Definition: ! ! PDF(X)(N,A) = Product ( 1 <= I <= N ) X(I)**( A(I) - 1 ) ! * Gamma ( A_SUM ) / A_PROD ! ! where ! ! 0 <= A(I) for all I; ! 0 <= X(I) for all I; ! Sum ( 1 <= I <= N ) X(I) = 1; ! A_SUM = Sum ( 1 <= I <= N ) A(I). ! A_PROD = Product ( 1 <= I <= N ) Gamma ( A(I) ) ! ! Modified: ! ! 06 December 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real X(N), the argument of the PDF. Each X(I) should ! be greater than 0.0E+00, and the X(I)'s must add up to 1.0. ! ! Input, integer N, the number of components. ! ! Input, real A(N), the probabilities for each component. ! Each A(I) should be nonnegative, and at least one should be ! positive. ! ! Output, real PDF, the value of the PDF. ! implicit none ! integer n ! real a(n) real a_prod real a_sum real gamma integer i real pdf real, parameter :: tol = 0.0001E+00 real x(n) real x_sum ! do i = 1, n if ( x(i) <= 0.0E+00 ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_PDF - Fatal error!' write ( *, * ) ' X(I) <= 0.' end if end do x_sum = sum ( x(1:n) ) if ( abs ( x_sum - 1.0E+00 ) > tol ) then write ( *, * ) ' ' write ( *, * ) 'DIRICHLET_PDF - Fatal error!' write ( *, * ) ' SUM X(I) =/= 1.' end if a_sum = sum ( a(1:n) ) a_prod = 1.0E+00 do i = 1, n a_prod = a_prod * gamma ( a(i) ) end do pdf = gamma ( a_sum ) / a_prod do i = 1, n pdf = pdf * x(i)**( a(i) - 1.0E+00 ) end do return end subroutine dirichlet_sample ( n, a, x ) ! !******************************************************************************* ! !! DIRICHLET_SAMPLE samples the Dirichlet PDF. ! ! ! Reference: ! ! Jerry Banks, editor, ! Handbook of Simulation, ! Engineering and Management Press Books, 1998, page 169. ! ! Modified: ! ! 23 November 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components. ! ! Input, real A(N), the probabilities for each component. !