rdirichlet: Generate Random Samples from a Dirichlet Distribution

View source: R/rdirichlet.R

rdirichletR Documentation

Generate Random Samples from a Dirichlet Distribution

Description

Generates random samples from a Dirichlet distribution using the Gamma representation: if Y_i \sim \mathrm{Gamma}(\alpha_i, 1) independently for i = 1, \ldots, K, then (Y_1 / S, \ldots, Y_K / S) \sim \mathrm{Dirichlet}(\alpha_1, \ldots, \alpha_K), where S = \sum_{i=1}^{K} Y_i.

Usage

rdirichlet(n, alpha)

Arguments

n

A positive integer specifying the number of random vectors to generate.

alpha

A numeric vector of length K \ge 2 containing positive concentration parameters of the Dirichlet distribution. All elements must be strictly positive.

Details

The Dirichlet distribution is a multivariate generalisation of the Beta distribution and is commonly used as a conjugate prior for multinomial proportions in Bayesian statistics.

The probability density function is:

f(x_1, \ldots, x_K) = \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}

where x_i > 0 and \sum_{i=1}^{K} x_i = 1.

Key properties:

  • Each marginal follows a Beta distribution: X_i \sim \mathrm{Beta}\!\left(\alpha_i,\, \sum_{l \neq i} \alpha_l\right).

  • E[X_i] = \alpha_i / \sum_{l=1}^{K} \alpha_l.

  • Components are negatively correlated unless one component dominates.

Implementation steps:

  1. Generate independent Y_i \sim \mathrm{Gamma}(\alpha_i, 1) for each i = 1, \ldots, K.

  2. Normalise: X_i = Y_i / \sum_{l=1}^{K} Y_l.

Value

A numeric matrix of dimensions n x K where each row is one random draw from the Dirichlet distribution, with all elements in [0, 1] and each row summing to 1. When n = 1, a numeric vector of length K is returned.

Examples

# Example 1: Generate 5 samples from Dirichlet(1, 1, 1) - uniform on simplex
samples <- rdirichlet(5, c(1, 1, 1))
print(samples)
rowSums(samples)  # Each row should sum to 1

# Example 2: Generate samples with unequal concentrations
samples <- rdirichlet(1000, c(2, 5, 3))
colMeans(samples)  # Expected values: approximately c(0.2, 0.5, 0.3)

# Example 3: Sparse Dirichlet (small alpha values)
samples <- rdirichlet(100, c(0.1, 0.1, 0.1, 0.1))
head(samples)  # Most weight concentrated on one component

# Example 4: Concentrated Dirichlet (large alpha values)
samples <- rdirichlet(100, c(100, 100, 100))
colMeans(samples)  # Concentrated around c(1/3, 1/3, 1/3)

# Example 5: Bayesian update with Jeffreys prior for 4 categories
prior_alpha     <- c(0.5, 0.5, 0.5, 0.5)
observed_counts <- c(10, 5, 8, 7)
posterior_samples <- rdirichlet(1000, prior_alpha + observed_counts)
colMeans(posterior_samples)  # Posterior mean


BayesianQDM documentation built on April 22, 2026, 1:09 a.m.