Probability Distributions for Data Science – Complete Guide with Python

Probability distributions are the mathematical language of uncertainty — every statistical inference, every machine learning model’s loss function, every A/B test, and every generative model rests on probability distributions. A data scientist who truly understands distributions can choose the right model for any data type, derive likelihood functions from first principles, diagnose when model assumptions are violated, and communicate uncertainty correctly. This guide covers the essential distributions every data scientist must know, their relationships, their applications in ML, and the Python patterns for working with them.

Probability distributions underpin the statistical tests in our Hypothesis Testing guide and A/B Testing guide, the Bayesian priors in our Bayesian Statistics guide, and the likelihood functions behind the models in our Machine Learning Interview Q&A. The Statistics Interview Q&A tests knowledge of these distributions directly — expected values, variance, moment-generating functions, and relationships between distributions appear frequently in technical screens.

Discrete Distributions

Bernoulli(p): Models a single binary trial — success (1) with probability p, failure (0) with probability 1-p. Mean = p, Variance = p(1-p). The fundamental building block: binary classification models predict Bernoulli probabilities; logistic regression’s log-likelihood is a sum of Bernoulli log-likelihoods.

Binomial(n, p): The number of successes in n independent Bernoulli(p) trials. P(X=k) = C(n,k) * p^k * (1-p)^(n-k). Mean = np, Variance = np(1-p). Used for: A/B test sample size calculations, quality control (defects in n items), click modelling (k clicks out of n impressions). As n grows large and p is small, Binomial approximates Poisson.

Poisson(lambda): Models the number of events in a fixed interval when events occur at a constant average rate lambda and independently. P(X=k) = e^{-lambda} * lambda^k / k!. Mean = Variance = lambda (equal mean and variance is the diagnostic fingerprint of Poisson). Used for: event count data (website visits per hour, customer arrivals per minute, errors per day). Over-dispersion (variance > mean) suggests Negative Binomial instead.

Geometric(p): Number of trials until the first success. Mean = 1/p. Memoryless property: given failure so far, the distribution of remaining trials is the same as the original. Used for: customer churn modelling (trials until a customer churns), session length modelling.

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

# --- Discrete distributions ---
n, p = 100, 0.15
binom = stats.binom(n=n, p=p)
print('Binomial mean:', binom.mean(), 'std:', binom.std())
print('P(X >= 20):', 1 - binom.cdf(19))
print('90th percentile:', binom.ppf(0.90))

lam = 4.2  # average arrivals per hour
poisson = stats.poisson(mu=lam)
print('P(X = 0):', poisson.pmf(0))       # probability of no arrivals
print('P(X >= 8):', 1 - poisson.cdf(7))  # probability of 8+ arrivals

# Fit Poisson to data
counts = np.array([3, 5, 2, 7, 4, 4, 6, 3, 5, 2])
lam_hat = counts.mean()   # MLE for Poisson is just the sample mean
print('Fitted lambda:', lam_hat)

# Test if data is Poisson (chi-squared goodness of fit)
observed = np.bincount(counts, minlength=12)
expected = stats.poisson.pmf(np.arange(12), lam_hat) * len(counts)
chi2, p_val = stats.chisquare(observed[observed > 0], expected[observed > 0])
print('Goodness of fit p-value:', round(p_val, 4))

Continuous Distributions

Normal (Gaussian) N(mu, sigma^2): The most important distribution in statistics. PDF: f(x) = (1/sqrt(2*pi*sigma^2)) * exp(-(x-mu)^2 / (2*sigma^2)). Properties: symmetric, characterised entirely by mean and variance, sum of n independent normals is normal (reproductive property), Central Limit Theorem says sample means converge to normal regardless of population distribution. Used for: linear regression residuals, feature distributions after standardisation, many natural phenomena. Z-score = (x – mu) / sigma converts to standard normal N(0,1).

Exponential(lambda): Models time between Poisson events. PDF: f(x) = lambda * e^{-lambda*x}. Mean = 1/lambda, Variance = 1/lambda^2. Memoryless: P(X > s+t | X > s) = P(X > t). Used for: survival analysis (time until failure), inter-arrival times, service times in queuing theory.

Log-Normal: If X ~ N(mu, sigma^2) then Y = e^X is log-normal. Right-skewed, positive-valued. Used for: income distributions, stock prices, website session durations, file sizes — quantities that are products of many small independent factors. Log-transform a log-normal variable to get a normal, making standard regression valid.

Beta(alpha, beta): Defined on [0,1], parameterised by shape parameters alpha and beta. Mean = alpha/(alpha+beta). Used for: modelling probabilities (prior for conversion rates in Bayesian A/B testing — see our Bayesian Statistics guide), representing proportions. Beta(1,1) = Uniform[0,1]. As alpha and beta grow large, Beta approaches Normal.

Fitting Distributions and Diagnostic Plots

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

# --- Fit and compare distributions to data ---
data = np.random.lognormal(mean=2.0, sigma=0.8, size=1000)

candidate_distributions = ['norm', 'lognorm', 'expon', 'gamma', 'beta']
results = []

for dist_name in candidate_distributions:
    dist = getattr(stats, dist_name)
    params = dist.fit(data)
    # Kolmogorov-Smirnov test: small p-value rejects the distribution
    ks_stat, ks_p = stats.kstest(data, dist_name, args=params)
    results.append((dist_name, ks_stat, ks_p, params))
    print(dist_name + ': KS=' + str(round(ks_stat, 4)) + ' p=' + str(round(ks_p, 4)))

# Best fit: distribution with lowest KS statistic (highest p-value)
best = min(results, key=lambda x: x[1])
print('Best fit:', best[0])

# --- Q-Q plot: visual normality check ---
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Raw data Q-Q plot
stats.probplot(data, dist='norm', plot=axes[0])
axes[0].set_title('Q-Q Plot: Raw Data vs Normal')

# Log-transformed
stats.probplot(np.log(data), dist='norm', plot=axes[1])
axes[1].set_title('Q-Q Plot: log(Data) vs Normal')

plt.tight_layout(); plt.show()

# --- Central Limit Theorem demonstration ---
population = stats.expon(scale=2.0)   # exponential population
sample_means = [population.rvs(size=30).mean() for _ in range(5000)]
# sample means are approximately N(2.0, 2.0/sqrt(30)) regardless of population shape
stats.probplot(sample_means, dist='norm', plot=plt); plt.show()

For the hypothesis tests that use these distributions as test statistics (t-distribution for t-tests, chi-squared for goodness-of-fit and independence tests, F-distribution for ANOVA), our Hypothesis Testing guide covers the mechanics. The Bayesian use of distributions as priors and posteriors — Beta-Binomial, Gaussian-Gaussian, Dirichlet-Multinomial conjugate pairs — is covered in our Bayesian Statistics guide. Interview questions on distributions that appear in technical screens at top companies are in our Statistics Interview Q&A.