Gamma distribution

The gamma distribution is a common choice for modeling continuous quantities that are non-negative: dollar amounts, waiting times, etc. It is governed by two positive parameters: the shape \(\alpha>0\) and the rate \(\beta>0\). Don’t worry about interpreting what those words mean. They’re just labels.

Basic properties

The gamma distribution gets its name because it is related to the so-called gamma function, which appears in the normalizing constant:

\[ \Gamma(\alpha) = \int_0^\infty x^{\alpha-1}e^{-x}\,\text{d}x. \]

We have a lil’ math explainer where you can read more about it.

Notation	\(X\sim\text{Gamma}(\alpha,\,\beta)\)
Range	\([0,\,\infty)\)
Parameter space	\(\alpha,\,\beta>0\)
PDF	\(f(x)=\begin{cases}\frac{\beta^\alpha}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x} &x\geq 0\\0 & x<0.\end{cases}\)
MGF	\(M(t)=\left(\frac{\beta}{\beta-t}\right)^\alpha\), \(t<\beta\)
Expectation	\(\alpha/\beta\)
Variance	\(\alpha/\beta^2\)

Special cases

Particular special cases of the gamma distribution get special nicknames:

• If \(X\sim\text{Gamma}(1,\,\beta)\), then \(X\) has the exponential distribution, denoted \(X\sim\text{Exp}(\beta)\) This special case is nice because we get closed-form expressions for the CDF and the quantile function:

\[ \begin{aligned} f(x)&=\beta e^{-\beta x}\\ F(x)&=1-e^{-\beta x}\\ F^{-1}(p)&=-\frac{\ln(1-p)}{\beta}. \end{aligned} \]

• If \(X\sim\text{Gamma}(\nu/2,\, 1/2)\), then \(X\) has the chi-squared distribution, denoted \(X\sim\chi^2_\nu\). In this context we call the parameter \(\nu>0\) the degrees of freedom of the chi-squared distribution.

• If \(X\sim\text{Gamma}(k,\,\beta)\) for some integer \(k\), then \(X\) has the Erlang distribution, denoted \(X\sim\text{Erlang}(k,\, \beta)\)

R commands

Here is the documentation for the suite of commands that let you work with the gamma distribution in R:

dgamma(x, shape, rate = 1) # PDF
pgamma(q, shape, rate = 1) # CDF: P(X <= q)
qgamma(p, shape, rate = 1) # quantile function (inverse CDF)
rgamma(n, shape, rate = 1) # random numbers

Play around!

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 700

library(shiny)

ui <- fluidPage(
  titlePanel("Gamma distribution CDF and PDF"),
  
  sidebarLayout(
    sidebarPanel(
      sliderInput("alpha", "Shape (α):",
                  min = 0.5, max = 5, value = 2, step = 0.1),
      sliderInput("beta", "Rate (β):",
                  min = 0.2, max = 2, value = 0.5, step = 0.05),
      hr(),
      verbatimTextOutput("moments")
    ),
    
    mainPanel(
      plotOutput("gammaPlot", height = "600px")
    )
  )
)

server <- function(input, output) {
  output$moments <- renderText({
    alpha <- input$alpha
    beta  <- input$beta
    mu <- alpha / beta
    sd <- sqrt(alpha) / beta
    paste0("E(X) = ", round(mu, 3),
           "\nsd(X) = ", round(sd, 3))
  })
  
  output$gammaPlot <- renderPlot({
    alpha <- input$alpha
    beta  <- input$beta
    
    # Fixed plotting range
    x_min <- -2
    x_max <- 25
    x <- seq(x_min, x_max, length.out = 2000)
    
    # PDF and CDF
    pdf_vals <- ifelse(x >= 0, dgamma(x, shape = alpha, rate = beta), 0)
    cdf_vals <- ifelse(x >= 0, pgamma(x, shape = alpha, rate = beta), 0)
    
    mu <- alpha / beta
    
    # Fixed y-limits
    pdf_ylim <- c(0, 0.8)
    cdf_ylim <- c(0, 1)
    
    par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
    
    # --- CDF ---
    plot(x, cdf_vals, type = "l", lwd = 2, col = "blue",
         xlim = c(x_min, x_max), ylim = cdf_ylim,
         main = "Cumulative Distribution Function (CDF)",
         xlab = "", ylab = "F(x)")
    abline(h = c(0, 1), col = "gray80", lty = 2)
    #abline(v = mu, col = "gray60", lty = 3)
    
    # --- PDF ---
    plot(x, pdf_vals, type = "l", lwd = 2, col = "darkred",
         xlim = c(x_min, x_max), ylim = pdf_ylim,
         main = "Probability Density Function (PDF)",
         xlab = "x", ylab = "f(x)")
    #abline(v = mu, col = "gray60", lty = 3)
    
    # Label E(X) with mtext along x-axis
    if (mu >= x_min && mu <= x_max) {
      mtext("E(X)",
            side = 1, line = 2.2, at = mu, cex = 1.0)
    }
  })
}

shinyApp(ui = ui, server = server)