Statistics for Data Science – Probability Distributions Explained 2026

Probability distributions are the mathematical foundation of statistics and machine learning. Every model you build makes assumptions about the distributions of its inputs and errors. Understanding distributions means understanding when model assumptions hold, how to simulate data, how to choose priors in Bayesian analysis, and how to interpret model outputs correctly. This guide covers the distributions every data scientist must know, with Python code and real-world intuition.

Normal (Gaussian) Distribution

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

fig, axes = plt.subplots(2, 3, figsize=(15, 9))

# Normal distribution
mu, sigma = 0, 1
x  = np.linspace(-4, 4, 300)
rv = stats.norm(loc=mu, scale=sigma)

axes[0, 0].plot(x, rv.pdf(x), 'b-', lw=2)
axes[0, 0].fill_between(x, rv.pdf(x), where=(x > -1) & (x < 1),
                          alpha=0.3, label='68% (±1σ)')
axes[0, 0].fill_between(x, rv.pdf(x), where=(x > -2) & (x < 2),
                          alpha=0.2, label='95% (±2σ)')
axes[0, 0].set_title('Normal Distribution N(0,1)'); axes[0, 0].legend()

# Key properties
print('Normal N(0,1):')
print(f'  PDF at x=0:    {rv.pdf(0):.4f}')
print(f'  CDF at x=1.96: {rv.cdf(1.96):.4f}')  # 97.5th percentile
print(f'  PPF at 0.975:  {rv.ppf(0.975):.4f}')  # 1.96 (z-score)
print(f'  95% CI:        {rv.interval(0.95)}')

# Sampling
samples = rv.rvs(size=10000, random_state=42)
print(f'  Sample mean:   {samples.mean():.4f}')
print(f'  Sample std:    {samples.std():.4f}')

# Applications
# Heights: N(170, 7) for adult males in India
heights = stats.norm(loc=170, scale=7)
print(f'
P(height > 185cm): {1 - heights.cdf(185):.4f}')
print(f'90th percentile:    {heights.ppf(0.90):.1f} cm')

Binomial Distribution

# Binomial: number of successes in n independent trials, each with probability p
n, p = 20, 0.3   # 20 trials, 30% success probability
rv_b = stats.binom(n=n, p=p)

k = np.arange(0, n+1)
axes[0, 1].bar(k, rv_b.pmf(k), color='steelblue', edgecolor='white')
axes[0, 1].axvline(rv_b.mean(), color='red', linestyle='--',
                    label=f'Mean={rv_b.mean():.1f}')
axes[0, 1].set_title(f'Binomial B(n={n}, p={p})')
axes[0, 1].legend()

print(f'
Binomial B(20, 0.3):')
print(f'  P(X=6):       {rv_b.pmf(6):.4f}')
print(f'  P(X≤6):       {rv_b.cdf(6):.4f}')
print(f'  P(X>10):      {1 - rv_b.cdf(10):.4f}')
print(f'  Mean:         {rv_b.mean():.1f} (= np = {n*p})')
print(f'  Std:          {rv_b.std():.4f} (= sqrt(np(1-p)) = {np.sqrt(n*p*(1-p)):.4f})')

# Application: A/B test — 300 visitors, 30% conversion rate
# What's the probability of getting >= 100 conversions?
print(f'
P(>=100 conversions in 300 visits | p=0.3): '
      f'{1 - stats.binom(300, 0.3).cdf(99):.4f}')

Poisson Distribution

# Poisson: number of events in a fixed interval when events occur at a constant rate
# Parameter λ = expected number of events

lambdas = [1, 4, 10]
k_range = np.arange(0, 25)

for lam in lambdas:
    rv_p = stats.poisson(mu=lam)
    axes[0, 2].plot(k_range, rv_p.pmf(k_range), 'o-',
                     markersize=4, label=f'λ={lam}')
axes[0, 2].set_title('Poisson Distribution')
axes[0, 2].legend()

rv_p4 = stats.poisson(mu=4)
print(f'
Poisson(λ=4): server errors per hour')
print(f'  P(X=0):  {rv_p4.pmf(0):.4f}  (no errors)')
print(f'  P(X=4):  {rv_p4.pmf(4):.4f}  (expected number)')
print(f'  P(X≥8):  {1 - rv_p4.cdf(7):.4f}  (unusually high load)')

# Application: call centre receives 10 calls/minute
rv_calls = stats.poisson(mu=10)
print(f'
Call centre (10 calls/min):')
print(f'  P(>15 calls) = {1 - rv_calls.cdf(15):.4f}')

Exponential Distribution

# Exponential: time between Poisson events; memoryless
# If events occur at rate λ, wait time ~ Exp(λ), mean = 1/λ
lam_exp = 0.5   # 0.5 events per unit time → mean wait = 2 units

rv_e = stats.expon(scale=1/lam_exp)
x_e  = np.linspace(0, 12, 200)
axes[1, 0].plot(x_e, rv_e.pdf(x_e), 'g-', lw=2)
axes[1, 0].fill_between(x_e, rv_e.pdf(x_e), alpha=0.3)
axes[1, 0].set_title(f'Exponential(λ={lam_exp})')

print(f'
Exponential(λ=0.5): time between events')
print(f'  Mean wait:  {rv_e.mean():.1f} units')
print(f'  P(wait<1):  {rv_e.cdf(1):.4f}')
print(f'  P(wait>4):  {1 - rv_e.cdf(4):.4f}')
print(f'  Median:     {rv_e.median():.4f}')

Beta Distribution

# Beta: models probabilities (values between 0 and 1)
# Alpha = prior successes + 1, Beta = prior failures + 1

x_b = np.linspace(0, 1, 200)
configs = [(1,1,'Uniform prior'), (2,5,'Prior: 1 success, 4 failures'),
           (10,30,'After 9 successes, 29 failures'), (50,150,'After 49s, 149f')]

for a, b, label in configs:
    rv_beta = stats.beta(a, b)
    axes[1, 1].plot(x_b, rv_beta.pdf(x_b), lw=2, label=label)
axes[1, 1].set_title('Beta Distributions (Bayesian A/B Testing)')
axes[1, 1].legend(fontsize=7)
axes[1, 1].set_xlabel('Conversion rate p')

# Bayesian A/B test
control_a, control_b       = 50, 150   # 50 conversions from 200 visitors
treatment_a, treatment_b   = 65, 135   # 65 conversions from 200 visitors

samples_ctrl = stats.beta(control_a, control_b).rvs(100_000, random_state=42)
samples_trt  = stats.beta(treatment_a, treatment_b).rvs(100_000, random_state=42)

p_treatment_better = (samples_trt > samples_ctrl).mean()
expected_lift       = (samples_trt / samples_ctrl - 1).mean()
print(f'
Bayesian A/B result:')
print(f'  Control rate:    {control_a / (control_a + control_b):.1%}')
print(f'  Treatment rate:  {treatment_a / (treatment_a + treatment_b):.1%}')
print(f'  P(treatment > control): {p_treatment_better:.3f}')
print(f'  Expected lift:   {expected_lift:.1%}')

Choosing the Right Distribution

# Fitting a distribution to data
data = np.random.gamma(2, 3, 1000)   # pretend we don't know the distribution

# Test multiple distributions
candidates = ['norm', 'gamma', 'lognorm', 'expon', 'weibull_min']
results    = []

for dist_name in candidates:
    dist  = getattr(stats, dist_name)
    params = dist.fit(data)
    ks_stat, p_value = stats.kstest(data, dist_name, args=params)
    results.append({'distribution': dist_name, 'KS stat': ks_stat,
                    'p-value': p_value, 'params': params})

results_df = pd.DataFrame(results).sort_values('KS stat')
print('
Distribution fit results (lower KS stat = better fit):')
print(results_df[['distribution', 'KS stat', 'p-value']].to_string(index=False))

Conclusion

Know which distribution matches your data’s generating process, not just its shape. Count data with a fixed maximum → Binomial. Count of rare events → Poisson. Time until next event → Exponential. Continuous symmetric data → Normal (after CLT kicks in). Proportions and probabilities → Beta. Positive skewed data → Log-normal or Gamma. Knowing the right distribution guides your choice of likelihood in Bayesian models, informs which statistical tests apply, and helps you simulate realistic data for testing your pipelines.