Beta Distribution in R: Complete Guide

Key Insights

The beta distribution is the workhorse for modeling probabilities and proportions in R, with its two shape parameters (α and β) giving you remarkable flexibility to represent prior beliefs and observed data
R’s four core beta functions (dbeta, pbeta, qbeta, rbeta) follow a consistent pattern you’ll recognize from other distributions, making them intuitive once you understand the naming convention
Bayesian A/B testing with beta distributions provides a more interpretable framework than frequentist hypothesis testing, letting you directly answer questions like “What’s the probability that variant B is better than A?”

Introduction to the Beta Distribution

The beta distribution is a continuous probability distribution bounded between 0 and 1, making it ideal for modeling probabilities, proportions, and rates. If you’re working with conversion rates, click-through rates, or any metric that represents a proportion, the beta distribution should be your first choice.

Two shape parameters control the distribution: alpha (α) and beta (β). Think of α as the number of successes plus one, and β as the number of failures plus one. When α = β = 1, you get a uniform distribution. When α > β, the distribution skews right (higher probabilities more likely). When α < β, it skews left.

# Visualize different beta distributions
library(ggplot2)

x <- seq(0, 1, length.out = 1000)

params <- data.frame(
  alpha = c(1, 2, 5, 2, 8),
  beta = c(1, 2, 2, 5, 2)
)

plot_data <- do.call(rbind, lapply(1:nrow(params), function(i) {
  data.frame(
    x = x,
    density = dbeta(x, params$alpha[i], params$beta[i]),
    params = paste0("α=", params$alpha[i], ", β=", params$beta[i])
  )
}))

ggplot(plot_data, aes(x = x, y = density, color = params)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Beta Distributions with Different Parameters",
       x = "x", y = "Density", color = "Parameters") +
  theme_minimal()

Core R Functions for Beta Distribution

R provides four essential functions for working with the beta distribution. Master these and you can handle any beta-related task.

dbeta() - Probability Density Function

Use dbeta() when you need the height of the density curve at a specific point. This is useful for likelihood calculations and plotting.

# Density at x = 0.3 for Beta(5, 2)
dbeta(0.3, shape1 = 5, shape2 = 2)
# [1] 0.5765

# Vectorized: densities at multiple points
dbeta(c(0.2, 0.5, 0.8), shape1 = 5, shape2 = 2)
# [1] 0.2048 1.8750 2.4576

pbeta() - Cumulative Distribution Function

Use pbeta() to find the probability that a random variable falls below a threshold. This answers questions like “What’s the probability our conversion rate is less than 10%?”

# P(X < 0.3) for Beta(5, 2)
pbeta(0.3, shape1 = 5, shape2 = 2)
# [1] 0.0579

# P(X > 0.7) = 1 - P(X < 0.7)
1 - pbeta(0.7, shape1 = 5, shape2 = 2)
# [1] 0.4718

# Or use lower.tail = FALSE
pbeta(0.7, shape1 = 5, shape2 = 2, lower.tail = FALSE)
# [1] 0.4718

qbeta() - Quantile Function

The quantile function is the inverse of the CDF. Use it to find credible intervals and percentiles.

# Find the median (50th percentile)
qbeta(0.5, shape1 = 5, shape2 = 2)
# [1] 0.7356

# 95% credible interval
qbeta(c(0.025, 0.975), shape1 = 5, shape2 = 2)
# [1] 0.3397 0.9629

rbeta() - Random Number Generation

Generate random samples from a beta distribution for simulations and Monte Carlo methods.

set.seed(42)
samples <- rbeta(10000, shape1 = 5, shape2 = 2)

# Verify sample statistics match theoretical values
mean(samples)  # Theoretical: 5/(5+2) = 0.714
# [1] 0.7139

var(samples)   # Theoretical: (5*2)/((5+2)^2*(5+2+1)) = 0.0255
# [1] 0.0254

Visualizing Beta Distributions

Effective visualization helps communicate uncertainty. Here’s how to create publication-ready graphics with ggplot2.

library(ggplot2)
library(dplyr)

# Create PDF and CDF plots side by side
x <- seq(0, 1, length.out = 500)
alpha <- 10
beta <- 3

viz_data <- data.frame(
  x = rep(x, 2),
  y = c(dbeta(x, alpha, beta), pbeta(x, alpha, beta)),
  type = rep(c("PDF", "CDF"), each = length(x))
)

ggplot(viz_data, aes(x = x, y = y)) +
  geom_line(linewidth = 1, color = "#2563eb") +
  geom_area(alpha = 0.2, fill = "#2563eb") +
  facet_wrap(~type, scales = "free_y") +
  labs(title = paste0("Beta(", alpha, ", ", beta, ") Distribution"),
       x = "Probability", y = NULL) +
  theme_minimal() +
  theme(strip.text = element_text(size = 12, face = "bold"))

For comparing distributions, layer them with transparency:

# Compare prior vs posterior
prior_alpha <- 2; prior_beta <- 2
post_alpha <- 12; post_beta <- 5

compare_data <- data.frame(
  x = rep(x, 2),
  density = c(dbeta(x, prior_alpha, prior_beta),
              dbeta(x, post_alpha, post_beta)),
  Distribution = rep(c("Prior", "Posterior"), each = length(x))
)

ggplot(compare_data, aes(x = x, y = density, fill = Distribution)) +
  geom_area(alpha = 0.5, position = "identity") +
  scale_fill_manual(values = c("Prior" = "#94a3b8", "Posterior" = "#2563eb")) +
  labs(title = "Prior vs Posterior Distribution",
       x = "Conversion Rate", y = "Density") +
  theme_minimal()

Parameter Estimation and Fitting

When you have proportion data, you’ll often need to estimate beta parameters. Two approaches dominate: method of moments and maximum likelihood.

Method of Moments

This approach matches sample moments to theoretical moments. It’s fast and provides reasonable starting estimates.

# Method of moments estimation
estimate_beta_mom <- function(x) {
  m <- mean(x)
  v <- var(x)
  
  # Solve for alpha and beta
  alpha <- m * (m * (1 - m) / v - 1)
  beta <- (1 - m) * (m * (1 - m) / v - 1)
  
  list(alpha = alpha, beta = beta)
}

# Test with simulated data
set.seed(123)
true_alpha <- 8
true_beta <- 3
sample_data <- rbeta(500, true_alpha, true_beta)

mom_est <- estimate_beta_mom(sample_data)
cat("MoM Estimates: α =", round(mom_est$alpha, 2), 
    ", β =", round(mom_est$beta, 2), "\n")
# MoM Estimates: α = 7.89 , β = 2.94

Maximum Likelihood Estimation

MLE is more efficient but requires numerical optimization. Use the fitdistrplus package for a robust implementation.

library(fitdistrplus)

# Fit beta distribution via MLE
fit <- fitdist(sample_data, "beta", method = "mle")
summary(fit)

# Extract parameters
mle_alpha <- fit$estimate["shape1"]
mle_beta <- fit$estimate["shape2"]
cat("MLE Estimates: α =", round(mle_alpha, 2), 
    ", β =", round(mle_beta, 2), "\n")
# MLE Estimates: α = 7.91 , β = 2.95

# Goodness-of-fit
gofstat(fit)

Bayesian Applications with Beta Priors

The beta distribution shines in Bayesian analysis because it’s the conjugate prior for binomial data. This means the posterior is also a beta distribution, making updates analytically tractable.

The update rule is elegant: if your prior is Beta(α, β) and you observe k successes in n trials, your posterior is Beta(α + k, β + n - k).

# A/B Testing Example
# Prior belief: conversion rates are around 5-15%
prior_alpha <- 3
prior_beta <- 30  # Gives mean of ~9%

# Observed data
# Variant A: 45 conversions out of 500 visitors
# Variant B: 62 conversions out of 500 visitors
a_conversions <- 45; a_visitors <- 500
b_conversions <- 62; b_visitors <- 500

# Posterior distributions
post_a_alpha <- prior_alpha + a_conversions
post_a_beta <- prior_beta + a_visitors - a_conversions
post_b_alpha <- prior_alpha + b_conversions
post_b_beta <- prior_beta + b_visitors - b_conversions

# Probability that B > A (Monte Carlo estimation)
set.seed(42)
n_sim <- 100000
samples_a <- rbeta(n_sim, post_a_alpha, post_a_beta)
samples_b <- rbeta(n_sim, post_b_alpha, post_b_beta)

prob_b_better <- mean(samples_b > samples_a)
cat("P(B > A) =", round(prob_b_better, 3), "\n")
# P(B > A) = 0.977

# Expected lift
expected_lift <- mean((samples_b - samples_a) / samples_a)
cat("Expected lift:", round(expected_lift * 100, 1), "%\n")
# Expected lift: 36.2 %

Practical Use Case: Conversion Rate Analysis

Let’s walk through a complete analysis of marketing conversion data.

library(ggplot2)

# Simulated campaign data
set.seed(2024)
campaigns <- data.frame(
  campaign = c("Email", "Social", "Search"),
  conversions = c(127, 89, 203),
  visitors = c(1250, 980, 1800)
)

# Weakly informative prior
prior_a <- 1
prior_b <- 1

# Calculate posteriors and credible intervals
results <- lapply(1:nrow(campaigns), function(i) {
  post_alpha <- prior_a + campaigns$conversions[i]
  post_beta <- prior_b + campaigns$visitors[i] - campaigns$conversions[i]
  
  data.frame(
    campaign = campaigns$campaign[i],
    mean_rate = post_alpha / (post_alpha + post_beta),
    ci_lower = qbeta(0.025, post_alpha, post_beta),
    ci_upper = qbeta(0.975, post_alpha, post_beta),
    post_alpha = post_alpha,
    post_beta = post_beta
  )
})

results_df <- do.call(rbind, results)
print(results_df)

# Visualize posterior distributions
x <- seq(0, 0.2, length.out = 500)
plot_data <- do.call(rbind, lapply(1:nrow(results_df), function(i) {
  data.frame(
    x = x,
    density = dbeta(x, results_df$post_alpha[i], results_df$post_beta[i]),
    campaign = results_df$campaign[i]
  )
}))

ggplot(plot_data, aes(x = x, y = density, fill = campaign)) +
  geom_area(alpha = 0.6, position = "identity") +
  geom_vline(data = results_df, aes(xintercept = mean_rate, color = campaign),
             linetype = "dashed", linewidth = 1) +
  scale_x_continuous(labels = scales::percent) +
  labs(title = "Posterior Conversion Rate Distributions",
       x = "Conversion Rate", y = "Density") +
  theme_minimal()

Summary and Best Practices

Use the beta distribution when modeling proportions, probabilities, or rates bounded between 0 and 1. For count data or unbounded proportions, consider alternatives like the Dirichlet (multivariate extension) or beta regression for modeling proportion outcomes with covariates.

Common pitfalls to avoid:

Don’t confuse shape parameters with mean and variance—convert explicitly
Remember that Beta(1, 1) is uniform, not uninformative in all contexts
When fitting to data at boundaries (0 or 1), consider zero-one-inflated beta models

Useful packages:

betareg: Beta regression for proportion outcomes
bayestestR: Bayesian analysis utilities including credible intervals
fitdistrplus: Robust distribution fitting with diagnostics

The beta distribution’s conjugacy with binomial data makes it invaluable for Bayesian A/B testing, where you can directly compute the probability that one variant outperforms another—a far more actionable insight than a p-value.