Gamma Distribution in Python: Complete Guide

Key Insights

The gamma distribution’s two parameters (shape and scale) give it remarkable flexibility for modeling positive continuous data like wait times, insurance claims, and equipment lifetimes
SciPy’s scipy.stats.gamma uses a shape/scale parameterization, while NumPy’s numpy.random.gamma uses shape/scale—but watch out for rate vs. scale confusion in different libraries
Parameter estimation via gamma.fit() makes it straightforward to fit real-world data, but always validate your fit visually and statistically before drawing conclusions

Introduction to the Gamma Distribution

The gamma distribution is one of the most versatile continuous probability distributions in statistics. It models positive real numbers and appears constantly in applied work: customer wait times, insurance claim amounts, rainfall accumulation, and equipment failure analysis.

What makes the gamma distribution powerful is its flexibility. With two parameters—shape (α or k) and scale (θ or β)—it can represent everything from exponential decay to bell-shaped curves shifted away from zero. The exponential distribution is actually a special case where shape equals 1.

If you’ve worked with Poisson processes, you’ve implicitly worked with gamma distributions. The waiting time until the k-th event in a Poisson process follows a gamma distribution. This connection makes it fundamental to queuing theory, reliability engineering, and survival analysis.

Mathematical Foundation

The gamma distribution’s probability density function (PDF) is:

$$f(x; \alpha, \theta) = \frac{x^{\alpha-1} e^{-x/\theta}}{\theta^\alpha \Gamma(\alpha)}$$

Where:

α (alpha) is the shape parameter
θ (theta) is the scale parameter
Γ(α) is the gamma function

The key statistics are straightforward:

Mean: μ = αθ
Variance: σ² = αθ²
Mode: (α - 1)θ for α ≥ 1

The shape parameter controls the distribution’s form. When α < 1, the PDF is J-shaped with a pole at zero. When α = 1, you get the exponential distribution. When α > 1, the distribution becomes bell-shaped and increasingly symmetric as α grows.

The scale parameter stretches or compresses the distribution along the x-axis without changing its fundamental shape.

Let’s visualize how these parameters affect the distribution:

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

x = np.linspace(0, 20, 500)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Varying shape parameter
shapes = [0.5, 1, 2, 5, 9]
for alpha in shapes:
    y = stats.gamma.pdf(x, a=alpha, scale=2)
    axes[0].plot(x, y, label=f'α={alpha}, θ=2')

axes[0].set_xlabel('x')
axes[0].set_ylabel('Probability Density')
axes[0].set_title('Effect of Shape Parameter (α)')
axes[0].legend()
axes[0].set_ylim(0, 0.5)

# Varying scale parameter
scales = [0.5, 1, 2, 4]
for theta in scales:
    y = stats.gamma.pdf(x, a=3, scale=theta)
    axes[1].plot(x, y, label=f'α=3, θ={theta}')

axes[1].set_xlabel('x')
axes[1].set_ylabel('Probability Density')
axes[1].set_title('Effect of Scale Parameter (θ)')
axes[1].legend()

plt.tight_layout()
plt.show()

This visualization reveals the intuition: shape controls whether the distribution is skewed or symmetric, while scale controls the spread.

Implementing Gamma Distribution with SciPy

SciPy’s scipy.stats.gamma provides a complete implementation. The key parameter mapping: SciPy uses a for shape and scale for the scale parameter.

from scipy import stats

# Create a gamma distribution with shape=3, scale=2
gamma_dist = stats.gamma(a=3, scale=2)

# Calculate PDF at specific points
x_values = np.array([1, 2, 5, 10])
pdf_values = gamma_dist.pdf(x_values)
print("PDF values:", pdf_values)

# Calculate CDF - probability that X <= x
cdf_values = gamma_dist.cdf(x_values)
print("CDF values:", cdf_values)

# Inverse CDF (percent point function) - find x for given probability
probabilities = np.array([0.25, 0.5, 0.75, 0.95])
quantiles = gamma_dist.ppf(probabilities)
print("Quantiles:", quantiles)

# Distribution statistics
print(f"Mean: {gamma_dist.mean()}")
print(f"Variance: {gamma_dist.var()}")
print(f"Standard deviation: {gamma_dist.std()}")

The ppf (percent point function) is particularly useful. It answers questions like “what value does 95% of the distribution fall below?” This is essential for setting thresholds and confidence intervals.

You can also calculate survival function (1 - CDF) directly:

# Probability that X > 10
survival_prob = gamma_dist.sf(10)
print(f"P(X > 10) = {survival_prob:.4f}")

# Log-probability (useful for numerical stability)
log_pdf = gamma_dist.logpdf(x_values)
print("Log-PDF values:", log_pdf)

Generating Random Samples

Both SciPy and NumPy can generate gamma-distributed random numbers. They’re functionally equivalent, but NumPy is faster for large samples.

import numpy as np
from scipy import stats

# Set seed for reproducibility
np.random.seed(42)

shape, scale = 3, 2
n_samples = 10000

# SciPy method
samples_scipy = stats.gamma.rvs(a=shape, scale=scale, size=n_samples)

# NumPy method
samples_numpy = np.random.gamma(shape=shape, scale=scale, size=n_samples)

# Visualize both
fig, axes = plt.subplots(1, 2, figsize=(12, 4))

for ax, samples, title in zip(axes, 
                               [samples_scipy, samples_numpy],
                               ['SciPy gamma.rvs()', 'NumPy random.gamma()']):
    ax.hist(samples, bins=50, density=True, alpha=0.7, label='Samples')
    
    # Overlay theoretical PDF
    x = np.linspace(0, 20, 200)
    ax.plot(x, stats.gamma.pdf(x, a=shape, scale=scale), 
            'r-', linewidth=2, label='Theoretical PDF')
    
    ax.set_xlabel('x')
    ax.set_ylabel('Density')
    ax.set_title(title)
    ax.legend()

plt.tight_layout()
plt.show()

# Compare sample statistics to theoretical values
print(f"Theoretical mean: {shape * scale}")
print(f"Sample mean (SciPy): {samples_scipy.mean():.3f}")
print(f"Sample mean (NumPy): {samples_numpy.mean():.3f}")

A warning about parameterization: some textbooks and software use a rate parameter (β = 1/θ) instead of scale. Always check documentation. SciPy and NumPy both use scale, but if you’re translating formulas from academic papers, you may need to invert the parameter.

Parameter Estimation from Data

When you have real data that you suspect follows a gamma distribution, gamma.fit() estimates the parameters using maximum likelihood estimation (MLE).

from scipy import stats
import numpy as np

# Generate synthetic "real" data
np.random.seed(123)
true_shape, true_scale = 4.5, 3.2
observed_data = np.random.gamma(true_shape, true_scale, size=500)

# Fit gamma distribution to data
# Returns: shape, loc, scale (loc is shift parameter, usually 0)
fitted_shape, fitted_loc, fitted_scale = stats.gamma.fit(observed_data)

print(f"True parameters: shape={true_shape}, scale={true_scale}")
print(f"Fitted parameters: shape={fitted_shape:.3f}, scale={fitted_scale:.3f}")
print(f"Location parameter: {fitted_loc:.3f}")

# If you know the distribution starts at 0, fix loc=0
fitted_shape_fixed, _, fitted_scale_fixed = stats.gamma.fit(observed_data, floc=0)
print(f"\nWith fixed loc=0: shape={fitted_shape_fixed:.3f}, scale={fitted_scale_fixed:.3f}")

# Visualize the fit
plt.figure(figsize=(10, 5))
plt.hist(observed_data, bins=40, density=True, alpha=0.7, label='Observed data')

x = np.linspace(0, observed_data.max(), 200)
plt.plot(x, stats.gamma.pdf(x, fitted_shape_fixed, scale=fitted_scale_fixed),
         'r-', linewidth=2, label=f'Fitted gamma (α={fitted_shape_fixed:.2f}, θ={fitted_scale_fixed:.2f})')

plt.xlabel('Value')
plt.ylabel('Density')
plt.title('Gamma Distribution Fit to Data')
plt.legend()
plt.show()

# Goodness-of-fit test (Kolmogorov-Smirnov)
ks_stat, p_value = stats.kstest(observed_data, 'gamma', 
                                 args=(fitted_shape_fixed, 0, fitted_scale_fixed))
print(f"\nKS test: statistic={ks_stat:.4f}, p-value={p_value:.4f}")

The floc=0 argument is important. By default, fit() estimates a location parameter that shifts the distribution. For most gamma applications, you want the distribution to start at zero.

Practical Application: Equipment Failure Analysis

Let’s work through a complete example: analyzing equipment failure times to predict maintenance needs.

import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt

# Simulate equipment failure data (time in days until failure)
np.random.seed(42)
failure_times = np.concatenate([
    np.random.gamma(shape=3.5, scale=45, size=150),  # Normal wear
    np.random.gamma(shape=2, scale=80, size=50)       # Some outliers
])

# Basic exploratory statistics
print("Failure Time Statistics")
print("-" * 30)
print(f"Sample size: {len(failure_times)}")
print(f"Mean: {failure_times.mean():.1f} days")
print(f"Median: {np.median(failure_times):.1f} days")
print(f"Std Dev: {failure_times.std():.1f} days")
print(f"Min: {failure_times.min():.1f} days")
print(f"Max: {failure_times.max():.1f} days")

# Fit gamma distribution
shape, loc, scale = stats.gamma.fit(failure_times, floc=0)
fitted_dist = stats.gamma(a=shape, scale=scale)

print(f"\nFitted Gamma Parameters:")
print(f"Shape (α): {shape:.3f}")
print(f"Scale (θ): {scale:.3f}")
print(f"Implied mean: {shape * scale:.1f} days")

# Create comprehensive visualization
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 1. Histogram with fitted PDF
ax1 = axes[0, 0]
ax1.hist(failure_times, bins=30, density=True, alpha=0.7, label='Observed')
x = np.linspace(0, failure_times.max(), 200)
ax1.plot(x, fitted_dist.pdf(x), 'r-', linewidth=2, label='Fitted Gamma')
ax1.set_xlabel('Days to Failure')
ax1.set_ylabel('Density')
ax1.set_title('Failure Time Distribution')
ax1.legend()

# 2. CDF comparison (empirical vs fitted)
ax2 = axes[0, 1]
sorted_data = np.sort(failure_times)
empirical_cdf = np.arange(1, len(sorted_data) + 1) / len(sorted_data)
ax2.step(sorted_data, empirical_cdf, label='Empirical CDF', where='post')
ax2.plot(x, fitted_dist.cdf(x), 'r-', linewidth=2, label='Fitted CDF')
ax2.set_xlabel('Days to Failure')
ax2.set_ylabel('Cumulative Probability')
ax2.set_title('CDF Comparison')
ax2.legend()

# 3. Q-Q plot for goodness of fit
ax3 = axes[1, 0]
theoretical_quantiles = fitted_dist.ppf(empirical_cdf)
ax3.scatter(theoretical_quantiles, sorted_data, alpha=0.5)
max_val = max(theoretical_quantiles.max(), sorted_data.max())
ax3.plot([0, max_val], [0, max_val], 'r--', label='Perfect fit')
ax3.set_xlabel('Theoretical Quantiles')
ax3.set_ylabel('Sample Quantiles')
ax3.set_title('Q-Q Plot')
ax3.legend()

# 4. Survival function (reliability)
ax4 = axes[1, 1]
ax4.plot(x, fitted_dist.sf(x), 'b-', linewidth=2)
ax4.axhline(y=0.5, color='gray', linestyle='--', alpha=0.5)
ax4.axhline(y=0.1, color='gray', linestyle='--', alpha=0.5)
median_life = fitted_dist.ppf(0.5)
ax4.axvline(x=median_life, color='green', linestyle='--', label=f'Median: {median_life:.0f} days')
ax4.set_xlabel('Days')
ax4.set_ylabel('Survival Probability')
ax4.set_title('Reliability Function')
ax4.legend()

plt.tight_layout()
plt.show()

# Practical maintenance calculations
print("\n" + "=" * 50)
print("MAINTENANCE PLANNING RECOMMENDATIONS")
print("=" * 50)

# Key percentiles for maintenance scheduling
percentiles = [0.05, 0.10, 0.25, 0.50]
for p in percentiles:
    days = fitted_dist.ppf(p)
    print(f"{int(p*100)}% of equipment fails by day {days:.0f}")

# Probability of failure within warranty period
warranty_days = 90
failure_prob = fitted_dist.cdf(warranty_days)
print(f"\nProbability of failure within {warranty_days}-day warranty: {failure_prob:.1%}")

# Recommended preventive maintenance interval (before 10% failure rate)
pm_interval = fitted_dist.ppf(0.10)
print(f"Recommended PM interval (90% reliability): {pm_interval:.0f} days")

This example demonstrates the full workflow: exploratory analysis, parameter estimation, visual validation, and practical decision-making based on the fitted model.

Summary and Further Resources

The gamma distribution is your go-to choice for modeling positive continuous data with right skew. Its flexibility through shape and scale parameters makes it applicable across domains.

When to use gamma vs. alternatives:

Exponential: When shape ≈ 1 (constant hazard rate)
Weibull: When you need more flexible hazard functions
Log-normal: When log-transformed data appears normal
Gamma: When you’re modeling sums of exponential processes or need a flexible positive distribution

Key takeaways:

Always fix floc=0 when fitting unless you have reason to shift the distribution
Validate fits visually with Q-Q plots, not just statistical tests
Use the survival function for reliability and time-to-event problems
Remember the parameterization differences between libraries

For deeper exploration, consult the SciPy documentation and consider studying the gamma distribution’s role in Bayesian statistics, where it serves as a conjugate prior for several likelihood functions.