29  Gamma Distribution

The Gamma distribution answers the question: how long until the \(k\)-th event? It generalizes the Exponential distribution — which covers the waiting time for a single event — to model the accumulated waiting time across multiple events, and provides a flexible family for any positive continuous measurement where right-skewness needs to be controlled.

Formally, the random variate \(X\) defined for the range \(X > 0\), is said to have a Gamma Distribution (i.e. \(X \sim \text{Gamma}(k, \lambda)\)) with shape parameter \(k > 0\) and rate parameter \(\lambda > 0\).

The Gamma distribution generalizes the Exponential distribution (shape \(k = 1\)) and the Erlang distribution (integer shape), and is related to the Chi-squared distribution. In R, the shape parameter is shape and the rate parameter is rate; the equivalent scale parameterization \(\theta = 1/\lambda\) is also accepted — pass scale \(= 1/\lambda\) as an argument instead of rate.

29.1 Probability Density Function

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

where \(\Gamma(k) = \int_0^\infty t^{k-1} e^{-t}\, dt\) is the Gamma function. For positive integers \(k\), \(\Gamma(k) = (k-1)!\).

The figure below shows examples of the Gamma Probability Density Function for different shape values with \(\lambda = 1\).

Code 
par(mfrow = c(2, 2))
x <- seq(0, 12, length = 500)

plot(x, dgamma(x, shape = 0.5, rate = 1), type = "l", lwd = 2, col = "blue",
     xlab = "x", ylab = "f(x)", main = expression(paste(k == 0.5, ",  ", lambda == 1)),
     ylim = c(0, 2))

plot(x, dgamma(x, shape = 1, rate = 1), type = "l", lwd = 2, col = "blue",
     xlab = "x", ylab = "f(x)", main = expression(paste(k == 1, ",  ", lambda == 1)))

plot(x, dgamma(x, shape = 2, rate = 1), type = "l", lwd = 2, col = "blue",
     xlab = "x", ylab = "f(x)", main = expression(paste(k == 2, ",  ", lambda == 1)))

plot(x, dgamma(x, shape = 5, rate = 1), type = "l", lwd = 2, col = "blue",
     xlab = "x", ylab = "f(x)", main = expression(paste(k == 5, ",  ", lambda == 1)))

par(mfrow = c(1, 1))
Figure 29.1: Gamma Probability Density Function for various shape values (rate = 1)

29.2 Purpose

The Gamma distribution models the total waiting time until the \(k\)-th event in a process where individual waiting times are independent and exponentially distributed. Its shape parameter \(k\) controls the degree of skewness, making it a flexible family for any positive continuous measurement. Common applications include:

  • Waiting time until the \(k\)-th customer arrival, failure, or detected signal
  • Aggregate insurance claims and total loss amounts
  • Rainfall amounts, flood volumes, and similar hydrology measurements
  • Processing and repair times in operations research and reliability
  • Prior distribution for rate and precision parameters in Bayesian models

Relation to the discrete setting. The Gamma\((k, \lambda)\) distribution is the continuous analog of the Negative Binomial distribution. The Negative Binomial counts the total number of discrete trials until the \(k\)-th success; the Gamma measures the continuous total waiting time until the \(k\)-th event in a Poisson process. When the shape \(k\) is a positive integer, this analogy is exact via the sum-of-Exponentials construction: just as a Negative Binomial count is the sum of \(k\) i.i.d. Geometric variates, a Gamma\((k,\lambda)\) variate is the sum of \(k\) i.i.d. Exp\((\lambda)\) variates.

29.3 Distribution Function

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

where \(\gamma(k, z) = \int_0^z t^{k-1} e^{-t}\, dt\) is the lower incomplete gamma function. The ratio \(\gamma(k, z)/\Gamma(k)\) is also called the regularized incomplete gamma function and is computed by pgamma() in R.

The figure below shows the Gamma Distribution Function for \(k = 3\) and \(\lambda = 1\).

Code 
x <- seq(0, 12, length = 500)
plot(x, pgamma(x, shape = 3, rate = 1), type = "l", lwd = 2, col = "blue",
     xlab = "x", ylab = "F(x)", main = "Gamma Distribution Function",
     sub = expression(paste(k == 3, ",  ", lambda == 1)))
Figure 29.2: Gamma Distribution Function (shape = 3, rate = 1)

29.4 Moment Generating Function

\[ M_X(t) = \left(\frac{\lambda}{\lambda - t}\right)^k, \quad t < \lambda \]

29.5 1st Uncentered Moment

\[ \mu_1' = \frac{k}{\lambda} \]

29.6 2nd Uncentered Moment

\[ \mu_2' = \frac{k(k+1)}{\lambda^2} \]

29.7 3rd Uncentered Moment

\[ \mu_3' = \frac{k(k+1)(k+2)}{\lambda^3} \]

29.8 4th Uncentered Moment

\[ \mu_4' = \frac{k(k+1)(k+2)(k+3)}{\lambda^4} \]

The general formula is \(\mu_n' = \lambda^{-n} \prod_{i=0}^{n-1}(k+i)\).

29.9 2nd Centered Moment

\[ \mu_2 = \frac{k}{\lambda^2} \]

29.10 3rd Centered Moment

\[ \mu_3 = \frac{2k}{\lambda^3} \]

29.11 4th Centered Moment

\[ \mu_4 = \frac{3k(k+2)}{\lambda^4} \]

29.12 Expected Value

\[ \text{E}(X) = \frac{k}{\lambda} \]

29.13 Variance

\[ \text{V}(X) = \frac{k}{\lambda^2} \]

29.14 Median

There is no closed-form expression for the median of the Gamma distribution in general. It is computed numerically via the qgamma function in R:

# Median for Gamma(k = 3, lambda = 1)
qgamma(0.5, shape = 3, rate = 1)
[1] 2.67406

29.15 Mode

\[ \text{Mo}(X) = \frac{k - 1}{\lambda} \quad \text{for } k \geq 1 \]

For \(0 < k < 1\), the density is strictly decreasing on \((0, \infty)\) and the mode is at the left boundary (\(x \to 0^+\)). For \(k = 1\), the mode coincides with the left boundary at 0 (Exponential distribution).

29.16 Coefficient of Skewness

\[ g_1 = \frac{2}{\sqrt{k}} \]

The Gamma distribution is always right-skewed. Skewness decreases as the shape parameter \(k\) increases, and the distribution approaches symmetry (and normality, by the CLT) as \(k \to \infty\).

29.17 Coefficient of Kurtosis

\[ g_2 = 3 + \frac{6}{k} \]

The excess kurtosis is \(6/k\), which is always positive, indicating heavier tails than the Normal distribution. As \(k \to \infty\), the kurtosis approaches 3 (Normal).

29.18 Parameter Estimation

The maximum likelihood estimators of \(k\) and \(\lambda\) do not have closed-form solutions. A method-of-moments starting point gives \(\tilde{k} = \bar{x}^2/s^2\) and \(\tilde{\lambda} = \bar{x}/s^2\), where \(s^2\) is the sample variance. The fitdistr function in R’s MASS package uses these as starting values for numerical optimization.

library(MASS)

# Method-of-moments starting values
set.seed(42)
x_data <- rgamma(200, shape = 3, rate = 0.5)

k_mom    <- mean(x_data)^2 / var(x_data)
lam_mom  <- mean(x_data) / var(x_data)
cat("Method-of-moments: shape =", round(k_mom, 3),
    "  rate =", round(lam_mom, 3), "\n")
Method-of-moments: shape = 2.807   rate = 0.491

29.19 R Module

29.19.1 RFC

The Gamma Distribution module is available in RFC under the menu “Distributions / Gamma Distribution”.

29.19.3 R Code

The following code demonstrates Gamma probability calculations:

# Probability density function: f(x)
dgamma(x = 6, shape = 3, rate = 0.5)

# Distribution function: P(X <= x)
pgamma(q = 6, shape = 3, rate = 0.5)

# Quantile function
qgamma(p = 0.5, shape = 3, rate = 0.5)

# Generate random Gamma numbers
set.seed(42)
rgamma(n = 10, shape = 3, rate = 0.5)
[1] 0.1120209
[1] 0.5768099
[1] 5.348121
 [1] 10.275115  3.373710  5.152571  2.117020  4.670036  3.136948 13.419833
 [8]  4.803647  2.155189  1.572504

To fit a Gamma distribution to observed data:

library(MASS)

# Example: waiting times (minutes) until 3rd customer arrival
set.seed(7)
waiting <- rgamma(100, shape = 3, rate = 0.5)

fit <- fitdistr(waiting, "gamma")
print(fit)
     shape         rate   
  4.12937589   0.63620274 
 (0.56198028) (0.09206918)

29.20 Example

A service counter has customers arriving according to a Poisson process at rate \(\lambda = 0.5\) arrivals per minute. The waiting time until the 3rd customer arrives follows \(X \sim \text{Gamma}(3, 0.5)\), with mean \(3/0.5 = 6\) minutes and standard deviation \(\sqrt{3}/0.5 \approx 3.46\) minutes.

k      <- 3
lambda <- 0.5

# P(X <= 8): probability the 3rd customer arrives within 8 minutes
cat("P(3rd arrival within 8 min):", round(pgamma(8, shape = k, rate = lambda), 4), "\n")

# P(X > 10): probability of waiting more than 10 minutes
cat("P(wait > 10 min):", round(1 - pgamma(10, shape = k, rate = lambda), 4), "\n")

# 95th percentile: 95% of the time the 3rd customer arrives before this time
cat("95th percentile (min):", round(qgamma(0.95, shape = k, rate = lambda), 2), "\n")
P(3rd arrival within 8 min): 0.7619 
P(wait > 10 min): 0.1247 
95th percentile (min): 12.59
Interactive Shiny app (click to load).
Open in new tab

29.21 Random Number Generator

A Gamma\((k, \lambda)\) variate with integer shape \(k\) can be generated as the sum of \(k\) independent Exp\((\lambda)\) variates (see Property 1 below). For non-integer \(k\), R uses Marsaglia and Tsang’s efficient acceptance-rejection algorithm.

set.seed(123)
n      <- 1000
k      <- 3
lambda <- 0.5

# Sum-of-exponentials method (integer k only)
x_sum <- rowSums(matrix(rexp(n * k, rate = lambda), nrow = n))

# Built-in function
x_rgamma <- rgamma(n, shape = k, rate = lambda)

cat("Sum-of-Exp: mean =", round(mean(x_sum), 4),
    "  var =", round(var(x_sum), 4), "\n")
cat("rgamma():   mean =", round(mean(x_rgamma), 4),
    "  var =", round(var(x_rgamma), 4), "\n")
cat("Theoretical: mean =", k/lambda,
    "  var =", k/lambda^2, "\n")
Sum-of-Exp: mean = 6.0951   var = 13.319 
rgamma():   mean = 5.9367   var = 10.9331 
Theoretical: mean = 6   var = 12
Code 
set.seed(123)
x <- rgamma(1000, shape = 3, rate = 0.5)
hist(x, breaks = 35, col = "steelblue", freq = FALSE,
     xlab = "x", main = "Gamma Random Numbers (n = 1000, shape = 3, rate = 0.5)")
curve(dgamma(x, shape = 3, rate = 0.5), add = TRUE, col = "red", lwd = 2)
legend("topright", legend = "Theoretical density", col = "red", lwd = 2)
Figure 29.3: Histogram of simulated Gamma random numbers (n = 1000, shape = 3, rate = 0.5)
Interactive Shiny app (click to load).
Open in new tab

29.22 Property 1: Sum of Independent Exponentials

The sum of \(k\) independent Exponential random variables, each with rate \(\lambda\), follows a Gamma distribution:

\[ \sum_{i=1}^{k} X_i \sim \text{Gamma}(k, \lambda) \quad \text{where } X_i \overset{\text{i.i.d.}}{\sim} \text{Exp}(\lambda) \]

This is the origin of the Gamma distribution’s interpretation as the waiting time until the \(k\)-th event in a Poisson process (see Chapter 18).

29.23 Property 2: Additivity

If \(X_1 \sim \text{Gamma}(k_1, \lambda)\) and \(X_2 \sim \text{Gamma}(k_2, \lambda)\) are independent, their sum is also Gamma distributed:

\[ X_1 + X_2 \sim \text{Gamma}(k_1 + k_2,\; \lambda) \]

This property requires both variables to share the same rate parameter \(\lambda\).

29.24 Property 3: Scaling

If \(X \sim \text{Gamma}(k, \lambda)\) and \(c > 0\), then:

\[ cX \sim \text{Gamma}\!\left(k,\; \frac{\lambda}{c}\right) \]

Scaling a Gamma variate changes only the rate (equivalently, the scale), leaving the shape parameter unchanged.