Binomial Distribution in Python: Complete Guide

Key Insights

The binomial distribution models the number of successes in a fixed number of independent trials, making it essential for A/B testing, quality control, and any scenario with binary outcomes.
Python’s scipy.stats.binom provides a complete toolkit for binomial calculations, while numpy.random.binomial excels at large-scale simulations.
Understanding when to use binomial versus normal approximation or Poisson distribution prevents common statistical errors in production code.

Introduction to Binomial Distribution

The binomial distribution answers a simple question: if you flip a biased coin n times, how likely are you to get exactly k heads? This seemingly basic concept underlies critical business decisions—from determining whether your new checkout flow actually improves conversions to deciding if a manufacturing batch meets quality standards.

A binomial distribution has two parameters: n (number of trials) and p (probability of success on each trial). The distribution is discrete, meaning it only produces integer values from 0 to n. Each trial must be independent, and the probability p must remain constant across all trials.

Real-world applications include A/B testing (did 47 conversions out of 500 visitors happen by chance?), quality control (how many defective items should we expect in a batch of 1000?), and clinical trials (is this drug’s success rate significantly better than placebo?).

Mathematical Foundation

The probability mass function (PMF) gives the probability of observing exactly k successes in n trials:

$$P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}$$

The binomial coefficient $\binom{n}{k}$ counts the number of ways to choose k successes from n trials. The mean of a binomial distribution is $\mu = np$, and the variance is $\sigma^2 = np(1-p)$.

Here’s how to calculate binomial probability from scratch:

import math

def binomial_probability(n: int, k: int, p: float) -> float:
    """Calculate P(X = k) for binomial distribution."""
    coefficient = math.comb(n, k)
    return coefficient * (p ** k) * ((1 - p) ** (n - k))

# Example: probability of exactly 3 heads in 10 fair coin flips
prob = binomial_probability(n=10, k=3, p=0.5)
print(f"P(X = 3) = {prob:.4f}")  # Output: P(X = 3) = 0.1172

This manual approach helps build intuition, but you should use scipy for production code—it handles edge cases and numerical stability.

Using scipy.stats.binom

The scipy.stats.binom module is your primary tool for binomial calculations. It provides methods for every common operation:

from scipy.stats import binom

# Define distribution parameters
n, p = 20, 0.3

# Create a frozen distribution (optional but cleaner)
dist = binom(n=n, p=p)

# PMF: P(X = k) - probability of exactly k successes
print(f"P(X = 6) = {dist.pmf(6):.4f}")  # 0.1916

# CDF: P(X ≤ k) - probability of at most k successes
print(f"P(X ≤ 6) = {dist.cdf(6):.4f}")  # 0.6080

# SF (survival function): P(X > k) - probability of more than k successes
print(f"P(X > 6) = {dist.sf(6):.4f}")  # 0.3920

# PPF (percent point function): inverse of CDF
# Find k such that P(X ≤ k) ≥ 0.95
print(f"95th percentile: {dist.ppf(0.95)}")  # 10.0

# RVS: generate random samples
samples = dist.rvs(size=1000, random_state=42)
print(f"Sample mean: {samples.mean():.2f}, Expected: {n * p}")

A common pattern is calculating the probability of a range of values:

# P(3 ≤ X ≤ 8)
prob_range = dist.cdf(8) - dist.cdf(2)  # Note: cdf(2) not cdf(3)
print(f"P(3 ≤ X ≤ 8) = {prob_range:.4f}")  # 0.8484

# Alternative using pmf sum
prob_range_alt = sum(dist.pmf(k) for k in range(3, 9))
print(f"Verification: {prob_range_alt:.4f}")  # 0.8484

Visualization Techniques

Visualizing distributions builds intuition and communicates results effectively. Here’s how to create publication-quality plots:

import numpy as np
import matplotlib.pyplot as plt

def plot_binomial_comparison():
    """Compare binomial distributions with different parameters."""
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    
    # PMF comparison
    ax1 = axes[0]
    params = [(20, 0.3, 'n=20, p=0.3'), 
              (20, 0.5, 'n=20, p=0.5'),
              (20, 0.7, 'n=20, p=0.7')]
    
    for n, p, label in params:
        k = np.arange(0, n + 1)
        pmf = binom.pmf(k, n, p)
        ax1.bar(k, pmf, alpha=0.5, label=label, width=0.8)
    
    ax1.set_xlabel('Number of Successes (k)')
    ax1.set_ylabel('Probability')
    ax1.set_title('PMF: Effect of Probability p')
    ax1.legend()
    
    # CDF comparison
    ax2 = axes[1]
    for n, p, label in params:
        k = np.arange(0, n + 1)
        cdf = binom.cdf(k, n, p)
        ax2.step(k, cdf, where='mid', label=label, linewidth=2)
    
    ax2.set_xlabel('Number of Successes (k)')
    ax2.set_ylabel('Cumulative Probability')
    ax2.set_title('CDF: Cumulative Distribution')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('binomial_comparison.png', dpi=150)
    plt.show()

plot_binomial_comparison()

For presentations, consider using seaborn’s styling:

import seaborn as sns

sns.set_style("whitegrid")
n, p = 50, 0.4
k = np.arange(0, n + 1)
pmf = binom.pmf(k, n, p)

plt.figure(figsize=(10, 5))
plt.bar(k, pmf, color=sns.color_palette("viridis", 1)[0], alpha=0.8)
plt.axvline(x=n*p, color='red', linestyle='--', label=f'Mean = {n*p}')
plt.xlabel('Number of Successes')
plt.ylabel('Probability')
plt.title(f'Binomial Distribution (n={n}, p={p})')
plt.legend()
plt.show()

Practical Applications

A/B Test Significance

The most common application: determining if your test results are statistically significant.

def ab_test_significance(
    control_visitors: int,
    control_conversions: int,
    treatment_visitors: int,
    treatment_conversions: int,
    alpha: float = 0.05
) -> dict:
    """
    Determine if treatment conversion rate is significantly 
    different from control using binomial test.
    """
    control_rate = control_conversions / control_visitors
    treatment_rate = treatment_conversions / treatment_visitors
    
    # Under null hypothesis, treatment has same rate as control
    # Calculate probability of observing treatment_conversions or more extreme
    # Two-tailed test
    dist = binom(n=treatment_visitors, p=control_rate)
    
    if treatment_conversions > treatment_visitors * control_rate:
        # Treatment appears better - test upper tail
        p_value = 2 * dist.sf(treatment_conversions - 1)
    else:
        # Treatment appears worse - test lower tail
        p_value = 2 * dist.cdf(treatment_conversions)
    
    p_value = min(p_value, 1.0)  # Cap at 1.0
    
    return {
        'control_rate': control_rate,
        'treatment_rate': treatment_rate,
        'lift': (treatment_rate - control_rate) / control_rate,
        'p_value': p_value,
        'significant': p_value < alpha,
        'conclusion': 'Reject null' if p_value < alpha else 'Fail to reject null'
    }

# Example: New checkout flow test
results = ab_test_significance(
    control_visitors=1000,
    control_conversions=50,      # 5% conversion rate
    treatment_visitors=1000,
    treatment_conversions=68     # 6.8% conversion rate
)

print(f"Control rate: {results['control_rate']:.1%}")
print(f"Treatment rate: {results['treatment_rate']:.1%}")
print(f"Lift: {results['lift']:.1%}")
print(f"P-value: {results['p_value']:.4f}")
print(f"Conclusion: {results['conclusion']}")

Quality Control Analysis

Manufacturing often needs to determine acceptable defect thresholds:

def defect_analysis(
    batch_size: int,
    historical_defect_rate: float,
    observed_defects: int
) -> dict:
    """Analyze whether observed defects indicate a process problem."""
    dist = binom(n=batch_size, p=historical_defect_rate)
    
    # Expected defects
    expected = batch_size * historical_defect_rate
    
    # Probability of seeing this many or more defects
    p_value = dist.sf(observed_defects - 1)
    
    # Control limits (99% of normal variation)
    lower_limit = dist.ppf(0.005)
    upper_limit = dist.ppf(0.995)
    
    return {
        'expected_defects': expected,
        'observed_defects': observed_defects,
        'p_value': p_value,
        'control_limits': (lower_limit, upper_limit),
        'in_control': lower_limit <= observed_defects <= upper_limit
    }

result = defect_analysis(
    batch_size=500,
    historical_defect_rate=0.02,
    observed_defects=18
)
print(f"Expected: {result['expected_defects']:.1f}")
print(f"Observed: {result['observed_defects']}")
print(f"Control limits: {result['control_limits']}")
print(f"Process in control: {result['in_control']}")

NumPy for Simulation

When you need to simulate thousands of experiments, numpy.random.binomial is faster than scipy:

import numpy as np

def monte_carlo_verification(n: int, p: float, num_simulations: int = 10000):
    """Compare empirical simulation results with theoretical values."""
    # Generate samples
    rng = np.random.default_rng(42)
    samples = rng.binomial(n=n, p=p, size=num_simulations)
    
    # Theoretical values
    theoretical_mean = n * p
    theoretical_var = n * p * (1 - p)
    theoretical_std = np.sqrt(theoretical_var)
    
    # Empirical values
    empirical_mean = samples.mean()
    empirical_std = samples.std()
    
    # Compare distributions
    print(f"Theoretical mean: {theoretical_mean:.4f}")
    print(f"Empirical mean:   {empirical_mean:.4f}")
    print(f"Theoretical std:  {theoretical_std:.4f}")
    print(f"Empirical std:    {empirical_std:.4f}")
    
    # Verify PMF empirically
    print("\nPMF Verification (selected values):")
    for k in [int(theoretical_mean - theoretical_std), 
              int(theoretical_mean), 
              int(theoretical_mean + theoretical_std)]:
        theoretical_prob = binom.pmf(k, n, p)
        empirical_prob = np.mean(samples == k)
        print(f"  P(X={k}): theoretical={theoretical_prob:.4f}, "
              f"empirical={empirical_prob:.4f}")
    
    return samples

samples = monte_carlo_verification(n=50, p=0.3, num_simulations=100000)

This simulation approach is invaluable for validating analytical results and exploring scenarios where closed-form solutions are difficult.

Summary and Best Practices

When to use binomial distribution:

Fixed number of trials known in advance
Each trial is independent
Constant probability of success
Binary outcomes only

When to consider alternatives:

Normal approximation: When n is large (np > 5 and n(1-p) > 5), use normal distribution for faster computation
Poisson distribution: When n is large and p is small (rare events)
Beta-binomial: When probability varies between trials

Quick reference:

Task	Function
P(X = k)	`binom.pmf(k, n, p)`
P(X ≤ k)	`binom.cdf(k, n, p)`
P(X > k)	`binom.sf(k, n, p)`
Find k for given probability	`binom.ppf(q, n, p)`
Generate samples	`binom.rvs(n, p, size=N)`
Fast simulation	`np.random.binomial(n, p, size=N)`

Common pitfalls to avoid:

Forgetting that CDF is P(X ≤ k), not P(X < k)
Assuming independence when trials are correlated
Using binomial when sample size is a significant fraction of population (use hypergeometric instead)
Not checking if normal approximation conditions are met before using it

The binomial distribution is foundational. Master it, and you’ll have the tools to make data-driven decisions in product development, operations, and beyond.