44  Dirichlet Distribution

The Dirichlet distribution is the multivariate generalization of the Beta distribution. It is a probability distribution over the \((K-1)\)-simplex — the set of \(K\)-dimensional vectors whose components are non-negative and sum to one. The Dirichlet distribution arises naturally as the conjugate prior for the Multinomial distribution in Bayesian inference.

Formally, a random vector \(\mathbf{X} = (X_1, \ldots, X_K)\) is said to follow a Dirichlet distribution (i.e. \(\mathbf{X} \sim \text{Dir}(\boldsymbol{\alpha})\)) with concentration parameters \(\alpha_1, \ldots, \alpha_K > 0\) if the components satisfy \(X_i \geq 0\) and \(\sum_{i=1}^K X_i = 1\).

44.1 Probability Density Function

\[ f(\mathbf{x}; \boldsymbol{\alpha}) = \frac{\Gamma\!\left(\sum_{i=1}^K \alpha_i\right)}{\prod_{i=1}^K \Gamma(\alpha_i)} \prod_{i=1}^K x_i^{\alpha_i - 1}, \quad \mathbf{x} \in S_K \]

where \(S_K = \{\mathbf{x} \in \mathbb{R}^K : x_i \geq 0,\, \sum x_i = 1\}\) is the \(K\)-simplex.

The figure below shows the Dirichlet density on the 3-simplex (ternary plot) for different concentration parameter configurations.

Code 
dirichlet_density <- function(x, alpha) {
  x <- as.matrix(x)
  log_const <- lgamma(sum(alpha)) - sum(lgamma(alpha))
  exp(log_const + rowSums(sweep(log(x), 2, alpha - 1, `*`)))
}

draw_simplex <- function(alpha, title_text) {
  vals <- seq(0.02, 0.96, length.out = 50)
  pts <- do.call(rbind, lapply(vals, function(x1) {
    x2 <- vals[vals < 0.98 - x1]
    if (!length(x2)) return(NULL)
    x3 <- 1 - x1 - x2
    keep <- x3 > 0.02
    if (!any(keep)) return(NULL)
    cbind(x1 = x1, x2 = x2[keep], x3 = x3[keep])
  }))
  if (is.null(pts)) return()
  dens <- dirichlet_density(pts, alpha)
  tri_x <- c(0, 1, 0.5, 0)
  tri_y <- c(0, 0, sqrt(3)/2, 0)
  plot(tri_x, tri_y, type = "n", asp = 1, axes = FALSE, ann = FALSE,
       xlim = c(-0.05, 1.05), ylim = c(-0.05, 0.92), main = title_text)
  polygon(tri_x, tri_y, border = "#2c3e50", lwd = 2, col = "grey98")
  px <- pts[, 2] + 0.5 * pts[, 3]
  py <- (sqrt(3)/2) * pts[, 3]
  pal <- colorRampPalette(c("#f7fbff", "#74a9cf", "#0570b0", "#08306b"))(100)
  z <- log(pmax(dens, 1e-300))
  idx <- cut(z, breaks = 100, include.lowest = TRUE, labels = FALSE)
  points(px, py, pch = 15, cex = 0.7, col = pal[idx])
}

par(mfrow = c(2, 2), mar = c(1, 1, 3, 1))
draw_simplex(c(2, 2, 2), expression(alpha == "(2, 2, 2)"))
draw_simplex(c(5, 5, 5), expression(alpha == "(5, 5, 5)"))
draw_simplex(c(0.5, 0.5, 0.5), expression(alpha == "(0.5, 0.5, 0.5)"))
draw_simplex(c(10, 2, 1), expression(alpha == "(10, 2, 1)"))
par(mfrow = c(1, 1))
Figure 44.1: Dirichlet density heatmaps on the 3-simplex for different concentration parameters

44.2 Purpose

The Dirichlet distribution models uncertainty over compositions — vectors of proportions that must sum to one. Common applications include:

  • Bayesian analysis with categorical or multinomial data (as the conjugate prior)
  • Topic modeling (Latent Dirichlet Allocation), where topic proportions follow a Dirichlet
  • Compositional data analysis: market shares, diet proportions, soil composition, alloy mixtures
  • Election forecasting: modeling vote share distributions across parties
  • Bayesian A/B/C testing: posterior probabilities over conversion rates for multiple variants

44.3 Moment Generating Function

The Dirichlet distribution does not have a standard moment generating function, but its moments are well known and given below.

44.4 Expected Value

\[ \text{E}(X_i) = \frac{\alpha_i}{\alpha_0}, \quad \alpha_0 = \sum_{j=1}^K \alpha_j \]

44.5 Variance

\[ \text{V}(X_i) = \frac{\alpha_i(\alpha_0 - \alpha_i)}{\alpha_0^2(\alpha_0 + 1)} \]

44.6 Covariance

\[ \text{Cov}(X_i, X_j) = -\frac{\alpha_i \alpha_j}{\alpha_0^2(\alpha_0 + 1)}, \quad i \neq j \]

44.7 Mode

When all \(\alpha_i > 1\):

\[ \text{Mo}(X_i) = \frac{\alpha_i - 1}{\alpha_0 - K} \]

When any \(\alpha_i \leq 1\), the density concentrates at simplex boundaries (corners or edges) and the mode lies on the boundary of the simplex.

44.8 Parameter Estimation

The maximum likelihood estimators for \(\boldsymbol{\alpha}\) have no closed form and require iterative methods (e.g., Newton-Raphson on the log-likelihood). A method-of-moments estimator uses:

\[ \tilde{\alpha}_0 = \frac{\bar{x}_1(1 - \bar{x}_1)}{s_1^2} - 1, \qquad \tilde{\alpha}_i = \tilde{\alpha}_0 \bar{x}_i \]

set.seed(42)
alpha_true <- c(3, 5, 2)
n <- 200
x <- t(replicate(n, {
  g <- rgamma(3, shape = alpha_true)
  g / sum(g)
}))
xbar <- colMeans(x)
s2 <- apply(x, 2, var)
alpha0_hat <- xbar[1] * (1 - xbar[1]) / s2[1] - 1
alpha_hat <- alpha0_hat * xbar
cat("MoM alpha:", round(alpha_hat, 3), "\n")
cat("True alpha:", alpha_true, "\n")
MoM alpha: 2.947 4.682 2.001 
True alpha: 3 5 2

44.9 R Module

44.9.1 RFC

The Dirichlet Distribution module is available in RFC under the menu Distributions / Dirichlet Distribution.

44.9.3 R Code

alpha <- c(3, 5, 2)
alpha0 <- sum(alpha)

# Expected values
cat("E(X):", alpha / alpha0, "\n")

# Variances
cat("Var(X):", alpha * (alpha0 - alpha) / (alpha0^2 * (alpha0 + 1)), "\n")

# Mode (all alpha > 1)
cat("Mode:", (alpha - 1) / (alpha0 - length(alpha)), "\n")

# Generate random Dirichlet samples via Gamma
set.seed(42)
g <- matrix(rgamma(10 * length(alpha), shape = rep(alpha, each = 10)), ncol = length(alpha))
x <- g / rowSums(g)
round(head(x), 4)
E(X): 0.3 0.5 0.2 
Var(X): 0.01909091 0.02272727 0.01454545 
Mode: 0.2857143 0.5714286 0.1428571 
       [,1]   [,2]   [,3]
[1,] 0.3482 0.5234 0.1283
[2,] 0.1765 0.7130 0.1105
[3,] 0.2782 0.7090 0.0127
[4,] 0.1470 0.7749 0.0781
[5,] 0.2080 0.6041 0.1879
[6,] 0.2243 0.4201 0.3557

44.10 Example

A market analyst models the share of three brands using a Dirichlet prior with \(\boldsymbol{\alpha} = (3, 5, 2)\). After observing 30 purchases with counts \((12, 14, 4)\), the posterior is \(\text{Dir}(15, 19, 6)\).

prior <- c(3, 5, 2)
counts <- c(12, 14, 4)
posterior <- prior + counts

cat("Prior mean:", round(prior / sum(prior), 4), "\n")
cat("Posterior mean:", round(posterior / sum(posterior), 4), "\n")
cat("Observed proportions:", round(counts / sum(counts), 4), "\n")
Prior mean: 0.3 0.5 0.2 
Posterior mean: 0.375 0.475 0.15 
Observed proportions: 0.4 0.4667 0.1333
Interactive Shiny app (click to load).
Open in new tab

44.11 Property 1: Marginals are Beta

Each marginal \(X_i\) follows a Beta distribution:

\[ X_i \sim \text{Beta}(\alpha_i,\; \alpha_0 - \alpha_i) \]

See Chapter 30.

44.12 Property 2: Aggregation Property

Merging categories preserves the Dirichlet: if \(\mathbf{X} \sim \text{Dir}(\alpha_1, \ldots, \alpha_K)\) and we define \(Y = X_1 + X_2\), then \((Y, X_3, \ldots, X_K) \sim \text{Dir}(\alpha_1 + \alpha_2, \alpha_3, \ldots, \alpha_K)\).

44.13 Property 3: Conjugate Prior for Multinomial

If \(\boldsymbol{\theta} \sim \text{Dir}(\boldsymbol{\alpha})\) is the prior on category probabilities and \(\mathbf{n} = (n_1, \ldots, n_K)\) are observed multinomial counts, then the posterior is:

\[ \boldsymbol{\theta} \mid \mathbf{n} \sim \text{Dir}(\alpha_1 + n_1, \ldots, \alpha_K + n_K) \]

See Chapter 17 and Chapter 113.