Cauchy Distribution in R: Complete Guide

Key Insights

The Cauchy distribution has no defined mean or variance, causing sample means to fail as estimators and breaking the Central Limit Theorem—use median and interquartile range instead.
R provides four built-in functions (dcauchy, pcauchy, qcauchy, rcauchy) that handle all standard distribution operations with location and scale parameters.
Maximum likelihood estimation works for Cauchy parameters, but requires numerical optimization since closed-form solutions don’t exist.

Introduction to the Cauchy Distribution

The Cauchy distribution is the troublemaker of probability theory. It looks innocent enough—a bell-shaped curve similar to the normal distribution—but it breaks nearly every statistical rule you’ve learned.

Mathematically, the Cauchy distribution has a probability density function:

$$f(x) = \frac{1}{\pi\gamma\left[1 + \left(\frac{x-x_0}{\gamma}\right)^2\right]}$$

Where $x_0$ is the location parameter (similar to mean) and $\gamma$ is the scale parameter (similar to standard deviation). The critical difference from normal distributions: the Cauchy has no defined mean or variance. The integrals simply don’t converge.

This “pathological” behavior makes Cauchy distributions appear in unexpected places: resonance phenomena in physics, the ratio of two independent standard normal variables, and as a robust prior in Bayesian statistics. Understanding how to work with Cauchy distributions in R prepares you for edge cases that break conventional statistical assumptions.

Let’s see why the Cauchy is problematic:

library(ggplot2)

# Compare Cauchy and Normal PDFs
x <- seq(-10, 10, length.out = 1000)
df <- data.frame(
  x = rep(x, 2),
  density = c(dnorm(x), dcauchy(x)),
  distribution = rep(c("Normal", "Cauchy"), each = length(x))
)

ggplot(df, aes(x = x, y = density, color = distribution)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Cauchy vs Normal: Heavy Tails Matter",
       x = "x", y = "Density") +
  theme_minimal() +
  scale_color_manual(values = c("Cauchy" = "#E74C3C", "Normal" = "#3498DB"))

The Cauchy’s heavier tails mean extreme values occur far more frequently than you’d expect from a normal distribution.

R’s Built-in Cauchy Functions

R provides four functions following the standard distribution naming convention:

Function	Purpose	Returns
`dcauchy()`	Probability density	PDF value at x
`pcauchy()`	Cumulative distribution	P(X ≤ x)
`qcauchy()`	Quantile function	x for given probability
`rcauchy()`	Random generation	Random samples

All functions accept location (default 0) and scale (default 1) parameters:

# Density at x = 0 for standard Cauchy
dcauchy(0)  # 0.3183099

# Density at x = 2 with location = 1, scale = 2
dcauchy(2, location = 1, scale = 2)  # 0.1273240

# Cumulative probability P(X <= 0)
pcauchy(0)  # 0.5

# Probability between -1 and 1 for standard Cauchy
pcauchy(1) - pcauchy(-1)  # 0.5 (only 50%!)

# 95th percentile
qcauchy(0.95)  # 6.313752

# Compare to normal: 95th percentile
qnorm(0.95)  # 1.644854

# Generate 5 random Cauchy values
set.seed(42)
rcauchy(5, location = 0, scale = 1)

Notice that the 95th percentile of a standard Cauchy (6.31) is nearly four times larger than the normal (1.64). This illustrates the heavy-tail problem.

Generating and Visualizing Cauchy Random Variables

When you generate Cauchy samples, prepare for outliers:

set.seed(123)
n <- 1000
cauchy_samples <- rcauchy(n)
normal_samples <- rnorm(n)

# Summary statistics reveal the problem
cat("Cauchy range:", range(cauchy_samples), "\n")
cat("Normal range:", range(normal_samples), "\n")

# Count extreme values (|x| > 10)
cat("Cauchy extremes:", sum(abs(cauchy_samples) > 10), "\n")
cat("Normal extremes:", sum(abs(normal_samples) > 10), "\n")

Visualizing requires careful handling of outliers:

library(ggplot2)

set.seed(456)
samples <- data.frame(
  value = c(rcauchy(5000), rnorm(5000)),
  distribution = rep(c("Cauchy", "Normal"), each = 5000)
)

# Truncate for visualization (Cauchy outliers would destroy the plot)
samples$value_truncated <- pmax(pmin(samples$value, 10), -10)

ggplot(samples, aes(x = value_truncated, fill = distribution)) +
  geom_histogram(bins = 100, alpha = 0.6, position = "identity") +
  facet_wrap(~distribution, scales = "free_y") +
  labs(title = "Cauchy vs Normal Samples (truncated at ±10)",
       x = "Value", y = "Count") +
  theme_minimal() +
  xlim(-10, 10)

A density comparison works better:

ggplot(samples, aes(x = value_truncated, color = distribution)) +
  geom_density(linewidth = 1) +
  labs(title = "Empirical Density Comparison",
       x = "Value", y = "Density") +
  theme_minimal() +
  xlim(-10, 10)

Working with Cauchy CDFs and Quantiles

The CDF and quantile functions behave normally—it’s the moments that cause problems:

# CDF plot
x <- seq(-10, 10, length.out = 500)
cdf_df <- data.frame(
  x = x,
  cauchy_cdf = pcauchy(x),
  normal_cdf = pnorm(x)
)

ggplot(cdf_df) +
  geom_line(aes(x = x, y = cauchy_cdf, color = "Cauchy"), linewidth = 1) +
  geom_line(aes(x = x, y = normal_cdf, color = "Normal"), linewidth = 1) +
  labs(title = "CDF Comparison", x = "x", y = "F(x)") +
  theme_minimal()

Q-Q plots verify whether data follows a Cauchy distribution:

set.seed(789)
test_samples <- rcauchy(500, location = 2, scale = 3)

# Q-Q plot against theoretical Cauchy
theoretical_quantiles <- qcauchy(ppoints(500), location = 2, scale = 3)
sample_quantiles <- sort(test_samples)

qq_df <- data.frame(
  theoretical = theoretical_quantiles,
  sample = sample_quantiles
)

# Truncate for visualization
qq_df <- qq_df[abs(qq_df$theoretical) < 50 & abs(qq_df$sample) < 50, ]

ggplot(qq_df, aes(x = theoretical, y = sample)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red", linewidth = 1) +
  labs(title = "Q-Q Plot: Sample vs Theoretical Cauchy",
       x = "Theoretical Quantiles", y = "Sample Quantiles") +
  theme_minimal()

For formal testing, use the Kolmogorov-Smirnov test:

set.seed(101)
samples <- rcauchy(200, location = 1, scale = 2)

# Test against hypothesized Cauchy(1, 2)
ks.test(samples, "pcauchy", location = 1, scale = 2)

Parameter Estimation Challenges

Here’s where the Cauchy distribution earns its reputation. The sample mean doesn’t converge:

set.seed(2024)
sample_sizes <- c(10, 100, 1000, 10000, 100000)
means <- sapply(sample_sizes, function(n) mean(rcauchy(n)))

cat("Sample means by size:\n")
print(data.frame(n = sample_sizes, sample_mean = round(means, 3)))

Run this multiple times—the sample mean jumps wildly regardless of sample size. The Law of Large Numbers fails.

Use median and IQR instead:

set.seed(303)
samples <- rcauchy(10000, location = 5, scale = 2)

# Robust estimators
estimated_location <- median(samples)
estimated_scale <- IQR(samples) / 2  # IQR = 2*scale for Cauchy

cat("True location: 5, Estimated:", round(estimated_location, 3), "\n")
cat("True scale: 2, Estimated:", round(estimated_scale, 3), "\n")

For maximum likelihood estimation, use MASS::fitdistr():

library(MASS)

set.seed(404)
samples <- rcauchy(1000, location = 3, scale = 1.5)

# MLE estimation
fit <- fitdistr(samples, "cauchy")
print(fit)

Or implement custom MLE with optim():

# Negative log-likelihood for Cauchy
neg_log_lik <- function(params, data) {
  location <- params[1]
  scale <- params[2]
  if (scale <= 0) return(Inf)
  -sum(dcauchy(data, location, scale, log = TRUE))
}

set.seed(505)
samples <- rcauchy(500, location = -2, scale = 4)

result <- optim(
  par = c(median(samples), IQR(samples)/2),
  fn = neg_log_lik,
  data = samples,
  method = "L-BFGS-B",
  lower = c(-Inf, 0.001)
)

cat("MLE estimates - Location:", round(result$par[1], 3), 
    "Scale:", round(result$par[2], 3), "\n")

Practical Applications and Simulations

Central Limit Theorem Failure:

The CLT states that sample means converge to normal. Not with Cauchy:

set.seed(606)
n_simulations <- 1000
sample_size <- 100

# Sample means from Cauchy samples
cauchy_means <- replicate(n_simulations, mean(rcauchy(sample_size)))
normal_means <- replicate(n_simulations, mean(rnorm(sample_size)))

# The distribution of Cauchy means is still Cauchy, not normal!
par(mfrow = c(1, 2))
hist(normal_means, breaks = 50, main = "Normal: Means → Normal", xlim = c(-1, 1))
hist(cauchy_means[abs(cauchy_means) < 10], breaks = 50, 
     main = "Cauchy: Means → Cauchy (truncated)")

Ratio of Normal Variables:

The ratio of two independent standard normal variables follows a Cauchy distribution:

set.seed(707)
n <- 50000

z1 <- rnorm(n)
z2 <- rnorm(n)
ratio <- z1 / z2

# Compare empirical to theoretical
x <- seq(-10, 10, length.out = 1000)

ggplot() +
  geom_histogram(aes(x = ratio[abs(ratio) < 10], y = after_stat(density)), 
                 bins = 100, fill = "steelblue", alpha = 0.7) +
  geom_line(aes(x = x, y = dcauchy(x)), color = "red", linewidth = 1) +
  labs(title = "Ratio of Normals = Cauchy Distribution",
       x = "Z1/Z2", y = "Density") +
  theme_minimal()

Bayesian Prior Example:

Cauchy priors are popular in Bayesian regression for their robustness:

# Simulate prior draws for a regression coefficient
set.seed(808)
prior_draws <- rcauchy(10000, location = 0, scale = 2.5)

# Cauchy(0, 2.5) is a common weakly informative prior
quantile(prior_draws, c(0.025, 0.5, 0.975))

Summary and Best Practices

Working with Cauchy distributions in R requires abandoning some statistical reflexes:

Never use sample mean or variance for estimation. Use median and IQR.
MLE works but requires numerical optimization—MASS::fitdistr() handles this cleanly.
Visualizations need truncation. Outliers will destroy your plots otherwise.
The CLT doesn’t apply. Averages of Cauchy variables remain Cauchy.
Use Cauchy when appropriate: ratio statistics, robust priors, modeling phenomena with extreme outliers.

Avoid Cauchy distributions when you need convergent moments or when standard statistical procedures must apply. Use them when heavy tails genuinely represent your data-generating process or when you want priors that accommodate extreme parameter values without overwhelming your likelihood.