randomdistmod module

Description

Random number generators and distributions module

This module contains subroutines and functions to compute the gamma and normal distribution, as well as a random number generator. Most of the subroutines have been translated from other languages to FORTRAN. The main references are

Quick access

Types:randomstate
Variables:zero, half, rng1, rng2, one, qsiz, coffs, cmul, vsmall
Routines:normal_cdf_inv(), gamma_cdf_inv(), qchisq_appr(), normal_01_cdf_inv(), ran_normal(), ran_seed(), gamma_cdf(), gamma_pdf(), normal_01_cdf(), r8poly_value_horner(), gamma_inc(), ranur(), refill(), ranu(), ran_gp(), gamma_log(), r8_gamma(), ran_gamma(), ran_gamma_gp()

Needed modules

Types

  • type randomdistmod/randomstate
    Type fields:
    • % q (qsiz) [integer]
    • % have [logical]
    • % carry [integer]
    • % xs [integer]
    • % indx [integer]
    • % xcng [integer]

Variables

  • randomdistmod/cmul [integer,parameter=69609]
  • randomdistmod/one [real,parameter=1.0]
  • randomdistmod/rng1 [real,parameter=1./(2.*real(huge(i4)))]
  • randomdistmod/zero [real,parameter=0.0]
  • randomdistmod/half [real,parameter=0.5]
  • randomdistmod/coffs [integer,parameter=123]
  • randomdistmod/rng2 [real,parameter=1./real(huge(i4))]
  • randomdistmod/qsiz [integer,parameter=10]

    41265_i4

  • randomdistmod/vsmall [real,parameter=tiny(1.)]

Subroutines and functions

function randomdistmod/ranu(state)

Generates a uniformly distributed random 4 byte integer with the range (-huge(i4),+huge(i4)) based on the 32-bit super KISS random number generator by George Marsaglia, published online and translated to Fortran 90 by user “mecej4” and Marsaglia, http://forums.silverfrost.com/viewtopic.php?t=1480 Further modifications to pass the complete state of the generator as an argument by J.O. Kaplan, 2011

Parameters:state [randomstate,inout,target] :: state of the uniform random number generator
Return:ranu [integer]
Called from:ranur(), ran_normal()
Call to:refill()
function randomdistmod/ranur(state)

generate a random number in the range (0,1)

Parameters:state [randomstate,inout,target] :: state of the uniform random number generator
Return:ranur [real]
Called from:ran_gp(), ran_gamma(), weathergen()
Call to:ranu()
function randomdistmod/refill(state)

reset a random state

Parameters:state [randomstate,inout,target] :: state of the uniform random number generator
Return:s [integer]
Called from:ranu()
subroutine randomdistmod/ran_seed(sval, state)

Set the seed of the random state

Parameters:
  • sval [integer,in] :: state of the uniform random number generator
  • state [randomstate,inout,target] :: the random state
Called from:

gwgen
subroutine randomdistmod/ran_normal(state, nval)

Sampler for the normal distribution centered at 0 with std. dev. of unity, based on Marsaglia polar method

Parameters:
  • state [randomstate,inout] :: state of the uniform random number generator
  • nval [real,out] :: output: The random number from the normal distribution
Called from:

ran_gamma(), weathergen()
Call to:

ranu()
function randomdistmod/ran_gamma(state, first, shape_bn, scale)

Select a random number from a Gamma distribution

adapted from the cpp adaptation of the Marsaglia & Tsang random gamma algorithm in: http://www.johndcook.com/SimpleRNG.cpp

Uses the algorithm in Marsaglia, G. and Tsang, W.W. (2000), A simple method for generating gamma variables, Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.

Parameters:
  • state [randomstate,inout] :: state of the uniform random number generator
  • first [logical,in] :: flag if this is the first call to the distribution with this shape
  • shape_bn [real,in]
  • scale [real,in] :: scale parameter of the Gamma distribution (theta = 1/beta)
Return:

ret
Called from:

ran_gamma_gp(), ran_gamma()
Call to:

ran_normal(), ranur(), ran_gamma()
function randomdistmod/ran_gamma_gp(state, first, shape_bn, scale, thresh, shape_gp, scale_gp)

Select a random number from a hybrid Gamma-GP distribution

Parameters:
  • state [randomstate,inout] :: state of the uniform random number generator
  • first [logical,in] :: flag if this is the first call to the distribution with this shape
  • shape_bn [real,in]
  • scale [real,in] :: scale parameter of the Gamma distribution (theta = 1/beta)
  • thresh [real,in] :: the threshold above which to choose the GP distribution
  • shape_gp [real,in] :: shape parameter of the GP distribution
  • scale_gp [real,in] :: scale parameter of the GP distribution
Return:

ret
Called from:

weathergen()
Call to:

ran_gamma(), ran_gp()
function randomdistmod/ran_gp(state, shape_bn, scale, loc)

Select a random number from a generalized pareto (GP) distribution

Parameters:
  • state [randomstate,inout] :: state of the uniform random number generator
  • shape_bn [real,in]
  • scale [real,in] :: scale parameter of the GP distribution (theta = 1/beta)
  • loc [real,in] :: the location of the GP distribution
Return:

ran_gp [real]
Called from:

ran_gamma_gp()
Call to:

ranur()
function randomdistmod/qchisq_appr(p, nu, g, tol)

chi-square approximation for the gamma_cdf_inv() function

Parameters:
  • p [real,in] :: the quantile
  • nu [real,in] :: twice the gamma shape
  • g [real,in] :: the logarithm of the gamma function at the gamma shape
  • tol [real,in] :: the tolerance for the approximation
Return:

qchisq_appr [real]
Called from:

gamma_cdf_inv()
Call to:

gamma_log(), normal_cdf_inv()
function randomdistmod/gamma_cdf_inv(p, alpha, scale)

Compute the quantile function of the gamma distribution.

This function is based on the Applied Statistics Algorithm AS 91 (“ppchi2”) and via pgamma(.) AS 239.

References
Best, D. J. and D. E. Roberts (1975). Percentage Points of the Chi-Squared Distribution. Applied Statistics 24, page 385.

Note

Compared to the original R function, we do not use the final newton step which might lead to values going to infinity for quantiles close to 1

Parameters:
  • p [real,in] :: the quantile between 0 and 1
  • alpha [real,in] :: the shape of the gamma distribution
  • scale [real,in] :: the scale of the gamma distribution
Return:

gamma_cdf_inv [real]
Called from:

weathergen()
Call to:

gamma_log(), qchisq_appr(), gamma_cdf()
subroutine randomdistmod/gamma_cdf(x, a, b, c, cdf)

Evaluate the Gamma CDF.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
02 January 2000 Extracted: June, 2016
Author:
John Burkardt Extracted by Philipp Sommer
Parameters:
  • x [real,in] :: the input value for which to compute the CDF
  • a [real,in] :: the location (< x) of the gamma distribution (usually 0)
  • b [real,in] :: the shape (> 0.0) of the distribution
  • c [real,in] :: the scale (>0.0) of the distribution
  • cdf [real,out] :: the returned value of the CDF
Called from:

gamma_cdf_inv(), weathergen()
Call to:

gamma_inc()
function randomdistmod/gamma_inc(p, x)

Compute the incomplete Gamma function.

Formulas:

\[\Gamma_{inc}(P, 0) = 0\]
\[\Gamma_{inc}(P, \infty) = 1.\]
\[\Gamma_{inc}(P,X) = \int_0^x{T^{P-1} \exp{(-T)} \mathrm{d}t} / \Gamma(P)\]
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 01 May 2001
  • Extracted: June, 2016
Author:
  • Original FORTRAN77 version by B L Shea.
  • FORTRAN90 version by John Burkardt
  • Extracted by Philipp Sommer
Reference:
BL Shea, Chi-squared and Incomplete Gamma Integral, Algorithm AS239, Applied Statistics, Volume 37, Number 3, 1988, pages 466-473.
Parameters:
  • p [real,in] :: the exponent parameter (0.0 < P)
  • x [real,in] :: the integral limit parameter. If X is less than or equal to 0, GAMMA_INC is returned as 0.
Return:

gamma_inc [real]
Called from:

gamma_cdf()
Call to:

normal_01_cdf(), gamma_log()
function randomdistmod/gamma_log(x)

Calculate the natural logarithm of GAMMA ( X ).

Computation is based on an algorithm outlined in references 1 and 2. The program uses rational functions that theoretically approximate \(\log(\Gamma(X))\) to at least 18 significant decimal digits. The approximation for 12 < X is from Hart et al, while approximations for X < 12.0D+00 are similar to those in Cody and Hillstrom, but are unpublished.

The accuracy achieved depends on the arithmetic system, the compiler, intrinsic functions, and proper selection of the machine dependent constants.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 16 June 1999
  • Extracted June, 2016
Author:
  • Original FORTRAN77 version by William Cody, Laura Stoltz.
  • FORTRAN90 version by John Burkardt.
  • Extracted by Philipp Sommer
Reference:
  • William Cody, Kenneth Hillstrom, Chebyshev Approximations for the Natural Logarithm of the Gamma Function, Mathematics of Computation, Volume 21, 1967, pages 198-203.
  • Kenneth Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May 1969.
  • John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, John Rice, Henry Thacher, Christoph Witzgall, Computer Approximations, Wiley, 1968.
Parameters:x [real,in] :: the argument of the Gamma function (> 0.0)
Return:gamma_log [real]
Called from:gamma_cdf_inv(), qchisq_appr(), gamma_inc()
subroutine randomdistmod/gamma_pdf(x, a, b, c, pdf)

Evaluate the Gamma PDF.

\[PDF(a,b,c;x) = \exp({-(x-a)/b}) \cdot ((x-a)/b)^{c-1} / (b \cdot \Gamma(c))\]
  • GAMMA_PDF(A,B,C;X), where C is an integer, is the Erlang PDF.
  • GAMMA_PDF(A,B,1;X) is the Exponential PDF.
  • GAMMA_PDF(0,2,C/2;X) is the Chi Squared PDF with C degrees of freedom.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 02 January 2000
  • Extracted: June, 2016
Author:
  • John Burkardt
  • Extracted by Philipp Sommer
Parameters:
  • x [real,in] :: the argument of the PDF. A <= X
  • a [real,in] :: the location of the peak; A is often chosen to be 0.0.
  • b [real,in] :: the “scale” parameter; 0.0 < B, and is often 1.0.
  • c [real,in] :: the “shape” parameter; 0.0 < C, and is often 1.0.
  • pdf [real,out] :: the returned value of the PDF.
Called from:

weathergen()
Call to:

r8_gamma()
function randomdistmod/r8_gamma(x)

Evaluate Gamma(X) for a real argument.

This routine calculates the gamma function for a real argument X.

Computation is based on an algorithm outlined in reference 1. The program uses rational functions that approximate the gamma function to at least 20 significant decimal digits. Coefficients for the approximation over the interval (1,2) are unpublished. Those for the approximation for 12 <= X are from reference 2.

Modified:
  • 11 February 2008
  • Extracted: June, 2016
Author:
  • Original FORTRAN77 version by William Cody, Laura Stoltz.
  • FORTRAN90 version by John Burkardt.
  • Extracted by Philipp Sommer
Reference:
  • William Cody, An Overview of Software Development for Special Functions, in Numerical Analysis Dundee, 1975, edited by GA Watson, Lecture Notes in Mathematics 506, Springer, 1976.
  • John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, John Rice, Henry Thatcher, Christoph Witzgall, Computer Approximations, Wiley, 1968, LC: QA297.C64.
Parameters:x [real,in] :: the argument of the function.
Return:r8_gamma [real]
Called from:gamma_pdf()
subroutine randomdistmod/normal_cdf_inv(cdf, a, b, x)

Invert the Normal CDF.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 23 February 1999
  • Extracted: November, 2016
Author:
  • John Burkardt
  • Extracted by Philipp Sommer
Parameters:
  • cdf [real,in] :: the value of the CDF. 0.0 <= CDF <= 1.0.
  • a [real,in] :: the mean of the pdf
  • b [real,in] :: the standard deviation of the pdf
  • x [real,out] :: the corresponding argument
Called from:

qchisq_appr()
Call to:

normal_01_cdf_inv()
subroutine randomdistmod/normal_01_cdf_inv(p, x)

Invert the standard normal CDF.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 05 June 2007
  • Extracted: November, 2016
Author:
  • Original FORTRAN77 version by Michael Wichura.
  • FORTRAN90 version by John Burkardt.
  • Extracted by Philipp Sommer
Reference:
Michael Wichura, Algorithm AS241: The Percentage Points of the Normal Distribution, Applied Statistics, Volume 37, Number 3, pages 477-484, 1988.

Note

The result is accurate to about 1 part in 10^16.

Parameters:
  • p [real,in] :: the value of the cumulative probability densitity function. 0 < P < 1.
  • x [real,out] :: the normal deviate value with the property that the probability of a
Called from:

normal_cdf_inv()
Call to:

r8poly_value_horner()
function randomdistmod/r8poly_value_horner(m, c, x)

Evaluate a polynomial using Horner’s method.

The polynomial

\[p(x) = c_0 + c_1 * x + c_2 * x^2 + ... + c_m * x^m\]

is to be evaluated at the value X.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 02 January 2014
  • Extracted: November, 2016
Author:
  • John Burkardt
  • Extracted by Philipp Sommer
Parameters:
  • m [integer,in,] :: the polynomial coefficients. C(I) is the coefficient of \(X^I\)
  • c (m + 1) [real,in]
  • x [real,in] :: the polynomial value
Return:

r8poly_value_horner [real]
Called from:

normal_01_cdf_inv()
subroutine randomdistmod/normal_01_cdf(x, cdf)

evaluate the Normal 01 CDF.

Licensing:
This code is distributed under the GNU LGPL license.
Modified:
  • 10 February 1999
  • Extracted: June, 2016
Author:
  • John Burkardt
  • Extracted by Philipp Sommer
Reference:
AG Adams, Algorithm 39, Areas Under the Normal Curve, Computer Journal, Volume 12, pages 197-198, 1969.
Parameters:
  • x [real,in] :: the argument of the CDF.
  • cdf [real,out] :: the value of the CDF.
Called from:

gamma_inc()