stdlib_stats_distribution_gamma

rvs_gamma - gamma distribution random variates

Status

Experimental

Description

The gamma distribution is a continuous probability distribution widely used in modeling waiting times, queue lengths, and other non-negative quantities. It is defined by a shape parameter (\(k\) or \(\alpha\)) that controls the distribution’s form and a rate parameter (\(\lambda\)) that controls its scale.

Without arguments, the function returns a random sample from the standard gamma distribution \(\Gamma(k=1, \lambda=1)\).

With a single argument for shape parameter, the function returns a random sample from the gamma distribution \(\Gamma(k=\text{shape}, \lambda=1)\).

With two arguments (shape, rate), the function returns a random sample from the gamma distribution \(\Gamma(k=\text{shape}, \lambda=\text{rate})\).

The optional loc parameter specifies the location (shift) of the distribution.

With three or more arguments including array_size, the function returns a rank-1 array of gamma distributed random variates.

For complex shape and rate parameters, the real and imaginary parts are sampled independently of each other.

Note

The algorithm used for generating gamma random variates is fundamentally limited to double precision.[1]

Syntax

result = rvs_gamma([[shape] , rate , loc , array_size])

Class

Impure elemental function

Arguments

shape: optional argument has intent(in) and is a scalar of type real or complex. If shape is real, its value must be positive. If shape is complex, both the real and imaginary components must be positive. This is the shape (or “form”) parameter \(k\) of the distribution. Default value is 1.0.

rate: optional argument has intent(in) and is a scalar of type real or complex. If rate is real, its value must be positive. If rate is complex, both the real and imaginary components must be positive. This is the rate parameter (\lambda) (inverse of scale). Default value is 1.0.

loc: optional argument has intent(in) and is a scalar of type real or complex. Specifies the location (shift) of the distribution with default value 0.0. The distribution support is x > loc.

array_size: optional argument has intent(in) and is a scalar of type integer with default kind.

Return value

The result is a scalar or rank-1 array with a size of array_size, and the same type as shape. If shape or rate is non-positive, the result is NaN.

Example

program example_gamma_rvs
  use stdlib_stats, only: random_seed
  use stdlib_stats, only: rgamma => rvs_gamma

  implicit none
  real :: shape_arr(2, 3, 4)
  complex :: cshape, crate
  integer :: seed_put, seed_get

  seed_put = 1234567
  call random_seed(seed_put, seed_get)

  ! single standard gamma random variate with shape=2.0
  print *, rgamma(2.0)
  ! 2.50538206

  ! gamma random variate with shape=3.0, rate=2.0 (loc=0.0 by default)
  print *, rgamma(3.0, 2.0)
  ! 1.30591583

  ! gamma random variate with shape=3.0, rate=2.0, loc=1.5
  print *, rgamma(3.0, 2.0, 1.5)
  ! 2.80591583

  ! a rank-3 array of 60 standard gamma random variates with shape=0.5
  shape_arr(:, :, :) = 0.5
  print *, rgamma(shape_arr)
  ! 1.03841162  1.33044529  0.912742674  0.131288037  0.638593793
  ! 1.03565669E-02  0.624804378  1.12179172  4.91380468E-02  6.69969944E-03
  ! 6.67014271E-02  0.132111162  0.101102419  0.648416579  1.14922595
  ! 2.29003578E-02  1.85964716E-04  1.21213868E-02  1.69112933
  ! 7.30440915E-02  0.395139128  0.182758048  0.427981257  0.985665262

  ! an array of 10 random variates with shape=0.5, rate=1.0 (loc=0.0 by default)
  print *, rgamma(0.5, 1.0, 10)
  ! 1.39297554E-04  0.296419382  0.352113068  2.80515051  3.65264394E-04
  ! 0.197743446  5.54569438E-02  9.30598825E-02  1.02596343  1.85311246

  ! an array of 10 random variates with shape=0.5, rate=1.0, loc=2.0
  print *, rgamma(0.5, 1.0, 10, 2.0)
  ! 2.00013971  2.29641938  2.35211301  4.80515051  2.00036526
  ! 2.19774345  2.05545694  2.09305988  3.02596343  3.85311246

  cshape = (3.0, 4.0)
  crate = (2.0, 0.7)
  ! single complex gamma random variate with real part of shape=3.0, rate=2.0;
  ! imaginary part of shape=4.0, rate=0.7 (loc=(0,0) by default)
  print *, rgamma(cshape, crate)
  ! (0.826188326,3.54749799)

end program example_gamma_rvs

pdf_gamma - gamma distribution probability density function

Status

Experimental

Description

The probability density function (pdf) of the single real variable gamma distribution is:

\[ f(x)= \frac{\lambda^{k}}{\Gamma (k)}x^{k-1}e^{-\lambda x} ,\quad x>0,\ k>0,\ \lambda>0 \]

where \(k\) is the shape parameter, \(\lambda\) is the rate parameter, and \(\Gamma(k)\) is the gamma function.

An optional loc parameter specifies the location (shift) of the distribution.

For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginary \(y\) parts, the joint probability density function is the product of the corresponding real and imaginary marginal pdfs:[2]

\[f(x+\mathit{i}y)=f(x)f(y)\]

Syntax

result = pdf_gamma(x, shape, rate , loc)

Class

Impure elemental function

Arguments

x: has intent(in) and is a scalar of type real or complex. The point at which to evaluate the pdf.

shape: has intent(in) and is a scalar of type real or complex. The shape parameter \(k\). If shape is real, its value must be positive. If shape is complex, both the real and imaginary components must be positive.

rate: has intent(in) and is a scalar of type real or complex. The rate parameter \(\lambda\). If rate is real, its value must be positive. If rate is complex, both the real and imaginary components must be positive.

loc: optional argument has intent(in) and is a scalar of type real or complex. The location (shift) parameter with default value 0.0.

All arguments must have the same type.

Return value

The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If shape or rate is non-positive, the result is NaN.

Example

program example_gamma_pdf
  use stdlib_stats, only: random_seed
  use stdlib_stats, only: rgamma => rvs_gamma, &
                                               gamma_pdf => pdf_gamma

  implicit none
  real :: x(2, 3, 4), shape(2, 3, 4), rate(2, 3, 4), loc
  complex :: cshape, crate
  integer :: seed_put, seed_get

  seed_put = 1234567
  call random_seed(seed_put, seed_get)

  ! probability density at x=1.0 with shape=1.0, rate=1.0 (loc=0.0 by default)
  print *, gamma_pdf(1.0, 1.0, 1.0)
  ! 0.367879450

  ! probability density at x=2.0 with shape=2.0, rate=2.0 (loc=0.0 by default)
  print *, gamma_pdf(2.0, 2.0, 2.0)
  ! 0.270670563

  ! probability density at x=1.0 with shape=2.0, rate=1.0, loc=0.5
  print *, gamma_pdf(1.0, 2.0, 1.0, 0.5)
  ! 0.303265381

  ! probability density at x=0.5 with shape=2.0, rate=1.0, loc=1.0 (x < loc)
  print *, gamma_pdf(0.5, 2.0, 1.0, 1.0)
  ! 0.00000000

  ! probability density at x=1.0 with invalid shape (shape <= 0)
  print *, gamma_pdf(1.0, -1.0, 1.0)
  ! NaN

  ! gamma random variates array with shape=2.0, rate=1.0
  shape(:, :, :) = 2.0
  rate(:, :, :) = 1.0
  x = reshape(rgamma(2.0, 1.0, 24), [2, 3, 4])

  ! a rank 3 gamma probability density array (loc=0.0 by default)
  print *, gamma_pdf(x, shape, rate)
  ! 0.204550430  0.320178866  0.274986655  0.348611295  0.101865448
  ! 0.102199331  0.358981341  0.223676488  0.254329354  0.356714427
  ! 0.267390072  0.305148095  0.367848188  7.26194456E-02  1.49471285E-02
  ! 0.246272027  0.360770017  0.339665830  0.101558588  0.358678699
  ! 0.224196941  0.359253854  7.56355673E-02  0.251869917

  ! a rank 3 gamma probability density array with loc=1.5
  loc = 1.5
  print *, gamma_pdf(x, shape, rate, loc)
  ! (different values due to shift)

  cshape = (1.0, 1.5)
  crate = (1.0, 2.0)
  ! complex gamma probability density function at (1.5, 1.0) with real part
  ! of shape=1.0, rate=1.0 and imaginary part of shape=1.5, rate=2.0
  ! (loc=(0,0) by default)
  print *, gamma_pdf((1.5, 1.0), cshape, crate)
  ! 9.63761061E-02

end program example_gamma_pdf

cdf_gamma - gamma distribution cumulative distribution function

Status

Experimental

Description

Cumulative distribution function (cdf) of the single real variable gamma distribution:

\[ F(x)= \frac{\gamma (k, \lambda x)}{\Gamma (k)},\quad x>0,\ k>0,\ \lambda>0 \]

where \(\gamma(k, z)\) is the lower incomplete gamma function, and \(\Gamma(k)\) is the gamma function. This is often referred to as the regularized gamma P function.

An optional loc parameter specifies the location (shift) of the distribution.

For a complex variable \(z=(x + y i)\) with independent real \(x\) and imaginary \(y\) parts, the joint cumulative distribution function is the product of the corresponding real and imaginary marginal cdfs:[2]

\[F(x+\mathit{i}y)=F(x)F(y)\]

Syntax

result = cdf_gamma(x, shape, rate , loc)

Class

Impure elemental function

Arguments

x: has intent(in) and is a scalar of type real or complex. The point at which to evaluate the cdf.

shape: has intent(in) and is a scalar of type real or complex. The shape parameter (k). If shape is real, its value must be positive. If shape is complex, both the real and imaginary components must be positive.

rate: has intent(in) and is a scalar of type real or complex. The rate parameter (\lambda). If rate is real, its value must be positive. If rate is complex, both the real and imaginary components must be positive.

loc: optional argument has intent(in) and is a scalar of type real or complex. The location (shift) parameter with default value 0.0.

All arguments must have the same type.

Return value

The result is a scalar or an array, with a shape conformable to the arguments, and the same type as the input arguments. If shape or rate is non-positive, the result is NaN.

Example

program example_gamma_cdf
  use stdlib_stats, only: random_seed
  use stdlib_stats, only: rgamma => rvs_gamma, &
                                               gamma_cdf => cdf_gamma

  implicit none
  real :: x(2, 3, 4), shape(2, 3, 4), rate(2, 3, 4), loc
  complex :: cshape, crate
  integer :: seed_put, seed_get

  seed_put = 1234567
  call random_seed(seed_put, seed_get)

  ! a standard gamma cumulative at x=1.0 with shape=0.5, rate=1.0 (loc=0.0 by default)
  print *, gamma_cdf(1.0, 0.5, 1.0)
  ! 0.842700839

  ! a cumulative at x=2.0 with shape=1.5, rate=2.0 (loc=0.0 by default)
  print *, gamma_cdf(2.0, 1.5, 2.0)
  ! 0.953988254

  ! cumulative at x=1.5 with shape=2.0, rate=1.0, loc=0.5
  print *, gamma_cdf(1.5, 2.0, 1.0, 0.5)
  ! 0.518520014

  ! cumulative at x=0.5 with shape=2.0, rate=1.0, loc=1.0 (x < loc)
  print *, gamma_cdf(0.5, 2.0, 1.0, 1.0)
  ! 0.00000000

  ! cumulative at x=1.0 with invalid rate (rate <= 0)
  print *, gamma_cdf(1.0, 2.0, -1.0)
  ! NaN

  ! gamma random variates array with shape=1.0, rate=1.0
  shape(:, :, :) = 1.0
  rate(:, :, :) = 1.0
  x = reshape(rgamma(1.0, 1.0, 24), [2, 3, 4])

  ! a rank 3 standard gamma cumulative array (loc=0.0 by default)
  print *, gamma_cdf(x, shape, rate)
  ! 0.710880339  0.472411335  0.578345954  0.383050948  0.870905757
  ! 0.870430350  0.170215249  0.677347481  0.620089889  0.161825046
  ! 4.17549349E-02  0.510665894  0.252201647  0.911497891  0.984424412
  ! 0.635621786  0.177783430  0.414842933  0.871342421  0.338317066
  ! 2.06879266E-02  0.335232288  0.907408893  0.624871135

  ! a rank 3 standard gamma cumulative array with loc=0.5
  loc = 0.5
  print *, gamma_cdf(x, shape, rate, loc)
  ! (different values due to shift)

  cshape = (0.7, 2.1)
  crate = (0.5, 1.0)
  ! complex gamma cumulative distribution at (0.5, 0.5) with real part of
  ! shape=0.7, rate=0.5 and imaginary part of shape=2.1, rate=1.0
  ! (loc=(0,0) by default)
  print *, gamma_cdf((0.5, 0.5), cshape, crate)
  ! 2.87349485E-02

end program example_gamma_cdf

---

[1]

Marsaglia, George, and Wai Wan Tsang. “The ziggurat method for generating random variables.” Journal of statistical software 5 (2000): 1-7.

[2] (1,2)

Miller, Scott, and Donald Childers. Probability and random processes: With applications to signal processing and communications. Academic Press, 2012 (p. 197).