Probability for Machine Learning: From Foundations to Production

Probability isn't just a math class you survived—it's the operating system of machine learning. Every prediction, every confidence score, every uncertainty estimate traces back to probability theory. In 2026, as models move into high-stakes domains like healthcare and autonomous driving, understanding probability isn't optional; it's a safety requirement.

Most ML tutorials gloss over the foundations, jumping straight to scikit-learn. That works until your model fails in production because of data drift or miscalibrated probabilities. This article bridges that gap: from the axioms of Kolmogorov to debugging a production classifier that's suddenly overconfident.

We'll cover discrete and continuous distributions, Bayes' theorem, maximum likelihood estimation, and the central limit theorem—all with code examples and production war stories. You'll learn not just the 'what' but the 'why' and the 'how to fix when it breaks.'

By the end, you'll think probabilistically about your models: not just 'what's the prediction?' but 'how certain are we?' That's the difference between a junior and a senior ML engineer.

The Axioms of Probability: Sample Space, Events, and Measures

Probability theory starts with a sample space Ω, the set of all possible outcomes of an experiment. For a fair six-sided die, Ω = {1,2,3,4,5,6}. An event is any subset of Ω — for example, "rolling an even number" corresponds to {2,4,6}. The probability measure P assigns a number between 0 and 1 to each event, with P(Ω) = 1. Kolmogorov's axioms require countable additivity: for disjoint events A₁, A₂, ..., P(∪Aᵢ) = ΣP(Aᵢ). This is the bedrock upon which all of machine learning's uncertainty quantification is built.

In practice, you rarely work with raw sample spaces. Instead, you define a sigma-algebra (a collection of events closed under complement and countable union) and a measure on it. For discrete spaces, the power set suffices. For continuous spaces like ℝ, you use the Borel sigma-algebra generated by open intervals. The probability measure then integrates a density function over those intervals. This measure-theoretic formalism lets you treat discrete and continuous distributions under one roof — critical when you mix Bernoulli and Gaussian components in a latent variable model.

A common pitfall: confusing probability of an event with probability density. P(X = x) = 0 for continuous random variables, but the density f_X(x) can be > 0. The probability that X falls in [a, b] is ∫_a^b f_X(x) dx. In ML, you evaluate likelihoods — which are densities for continuous outputs — and treat them as probabilities in loss functions. This works because maximizing the density at observed points is equivalent to minimizing cross-entropy, but you must remember that density values can exceed 1 and are not probabilities.

The law of total probability and Bayes' rule both follow directly from the axioms. The chain rule for joint distributions — P(A ∩ B) = P(A|B)P(B) — is the foundation of probabilistic graphical models. Every time you factor a joint distribution in a Bayesian network, you're applying this axiomatically derived identity. Ignore the axioms at your peril: they enforce consistency, preventing nonsense like probabilities summing to 1.1.

import itertools

def empirical_probability(outcomes, event_condition):
    """Estimate P(event) from a list of outcomes."""
    total = len(outcomes)
    favorable = sum(1 for o in outcomes if event_condition(o))
    return favorable / total if total > 0 else 0.0

# Simulate 10000 die rolls
import random
random.seed(42)
rolls = [random.randint(1, 6) for _ in range(10000)]

# Event: even number
even_prob = empirical_probability(rolls, lambda x: x % 2 == 0)
print(f"P(even) ≈ {even_prob:.4f} (theoretical: 0.5)")

# Event: number > 4
high_prob = empirical_probability(rolls, lambda x: x > 4)
print(f"P(>4) ≈ {high_prob:.4f} (theoretical: 1/3 ≈ 0.3333)")

# Check additivity: P(even ∪ >4) = P(even) + P(>4) - P(even ∩ >4)
even_and_high = empirical_probability(rolls, lambda x: x % 2 == 0 and x > 4)
union_prob = even_prob + high_prob - even_and_high
print(f"P(even ∪ >4) ≈ {union_prob:.4f}")

Random Variables: Discrete vs. Continuous

A random variable X is a function from the sample space Ω to ℝ. It's not random in the sense of unpredictable; it's a deterministic mapping that assigns a number to each outcome. For a coin flip, you might define X(heads)=1, X(tails)=0. The randomness lives in the underlying probability measure on Ω, not in X itself. This distinction matters when you're building generative models: you choose the mapping, but the data-generating process determines the distribution.

Discrete random variables take on a countable set of values. The probability mass function (PMF) p(x) = P(X=x) sums to 1 over all x. Examples: Bernoulli (0 or 1), Binomial (number of successes in n trials), Poisson (count of events in a fixed interval). In ML, discrete variables appear as class labels, word counts, or latent cluster assignments. When you compute cross-entropy loss for classification, you're evaluating the negative log-likelihood under a categorical PMF.

Continuous random variables take on uncountably many values. The probability density function (PDF) f(x) is not a probability — it's a density. P(a ≤ X ≤ b) = ∫_a^b f(x) dx. The PDF can be > 1 (e.g., a narrow Gaussian with σ=0.1 has f(μ) ≈ 4). Common continuous distributions: Gaussian (real-valued data), Exponential (waiting times), Beta (probabilities). In regression, you typically assume the target is Gaussian conditioned on features: y|x ~ N(μ(x), σ²). The loss is negative log-likelihood, which for a Gaussian reduces to mean squared error up to a constant.

Mixed random variables combine discrete and continuous components. Example: a zero-inflated model for insurance claims — P(X=0) = π, and for X>0, X follows a continuous distribution. In practice, you handle these by splitting the likelihood into a point mass at zero and a density for positive values. Pyro and TensorFlow Probability support such distributions natively. Ignoring the mixed nature leads to biased estimates and poor calibration.

import numpy as np
from scipy import stats

# Discrete: Bernoulli PMF
p = 0.7
bernoulli = stats.bernoulli(p)
print(f"P(X=1) = {bernoulli.pmf(1):.4f}")
print(f"P(X=0) = {bernoulli.pmf(0):.4f}")

# Continuous: Gaussian PDF (not a probability!)
mu, sigma = 0.0, 0.5
gaussian = stats.norm(mu, sigma)
x_vals = np.array([-0.5, 0.0, 0.5])
print(f"PDF at x=0: {gaussian.pdf(0):.4f} (can be >1)")
print(f"P(-0.1 < X < 0.1) = {gaussian.cdf(0.1) - gaussian.cdf(-0.1):.4f}")

# Mixed: zero-inflated Gaussian (manual)
def zero_inflated_log_prob(x, pi, mu, sigma):
    if x == 0:
        return np.log(pi + (1-pi) * stats.norm.pdf(0, mu, sigma))
    else:
        return np.log(1-pi) + stats.norm.logpdf(x, mu, sigma)

print(f"Log-prob of 0: {zero_inflated_log_prob(0.0, 0.3, 0.0, 1.0):.4f}")
print(f"Log-prob of 1.5: {zero_inflated_log_prob(1.5, 0.3, 0.0, 1.0):.4f}")

Probability Distributions: Bernoulli, Gaussian, Poisson, and Beyond

Probability distributions are parametric families that describe how probability mass or density is spread over the support of a random variable. The Bernoulli distribution models a single binary trial: X ~ Bern(p) with P(X=1)=p, P(X=0)=1-p. Its mean is p, variance p(1-p). In ML, the Bernoulli is the building block for binary classification — the model outputs logits, and the loss is binary cross-entropy: -[y log(p̂) + (1-y) log(1-p̂)]. This is exactly the negative log-likelihood of a Bernoulli.

The Gaussian (normal) distribution N(μ, σ²) is the workhorse of continuous modeling. Its PDF is (1/√(2πσ²)) exp(-(x-μ)²/(2σ²)). The Gaussian is the maximum entropy distribution given a fixed mean and variance — meaning it makes the fewest assumptions beyond those constraints. This is why linear regression with squared error implicitly assumes Gaussian noise. The central limit theorem explains why sums of independent random variables converge to a Gaussian, justifying its ubiquity.

The Poisson distribution models count data: X ~ Pois(λ) with PMF e^{-λ} λ^x / x!. It's used for rare events: number of website visits per minute, defects per unit, etc. Its mean equals its variance (λ). In practice, real count data often has variance > mean (overdispersion), requiring a Negative Binomial distribution instead. Never blindly use Poisson without checking the mean-variance relationship.

Beyond these, you'll encounter: Beta (prior for Bernoulli probabilities), Gamma (waiting times, conjugate prior for Poisson rate), Dirichlet (prior for categorical probabilities), and Student-t (robust to outliers). The exponential family unifies many of these: distributions of the form p(x|θ) = h(x) exp(η(θ)·T(x) - A(θ)). This structure enables efficient learning via generalized linear models and variational inference. Understanding which distribution matches your data's generative process is the single most impactful modeling decision you'll make.

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

# Bernoulli
p = 0.3
bern = stats.bernoulli(p)
samples = bern.rvs(10000)
print(f"Bernoulli mean: {samples.mean():.4f} (theoretical: {p})")

# Gaussian
mu, sigma = 5.0, 2.0
norm = stats.norm(mu, sigma)
samples = norm.rvs(10000)
print(f"Gaussian mean: {samples.mean():.4f}, std: {samples.std():.4f}")

# Poisson
lam = 3.0
pois = stats.poisson(lam)
samples = pois.rvs(10000)
print(f"Poisson mean: {samples.mean():.4f}, var: {samples.var():.4f} (theoretical: {lam})")

# Negative Binomial for overdispersed counts
r, p_nb = 5, 0.5
nbinom = stats.nbinom(r, p_nb)
samples = nbinom.rvs(10000)
print(f"NegBinomial mean: {samples.mean():.4f}, var: {samples.var():.4f} (theoretical mean: {r*(1-p_nb)/p_nb})")

Bayes' Theorem: The Engine of Inference

Bayes' theorem is a direct consequence of the definition of conditional probability: P(A|B) = P(B|A)P(A) / P(B). In ML, we write it as P(θ|D) = P(D|θ)P(θ) / P(D), where θ are parameters and D is data. The posterior P(θ|D) combines prior beliefs P(θ) with the likelihood P(D|θ). The marginal likelihood P(D) = ∫ P(D|θ)P(θ) dθ normalizes the posterior. This framework is the foundation of Bayesian inference, probabilistic programming, and many regularization techniques.

Consider a simple example: you have a coin that may be biased. Prior: P(bias = 0.5) = 0.8, P(bias = 0.7) = 0.2. You flip it 10 times and get 8 heads. The likelihood under bias=0.5 is C(10,8)(0.5)^10 ≈ 0.044; under bias=0.7 it's C(10,8)(0.7)^8(0.3)^2 ≈ 0.233. The posterior: P(bias=0.7|data) = (0.233 0.2) / (0.0440.8 + 0.233*0.2) ≈ 0.57. The data shifted your belief from 20% to 57% for the biased coin. This is Bayesian updating in action.

In practice, the marginal likelihood is often intractable. Markov Chain Monte Carlo (MCMC) and variational inference approximate the posterior without computing P(D). Modern probabilistic programming languages (Pyro, Stan, TensorFlow Probability) handle this automatically. For large-scale ML, Bayesian neural networks place priors on weights and use variational dropout or Laplace approximation to estimate uncertainty.

Bayes' theorem also underpins the concept of conjugate priors: if the posterior has the same form as the prior, inference simplifies dramatically. For example, Beta is conjugate to Bernoulli: Beta(α, β) prior + Binomial data → Beta(α + successes, β + failures) posterior. This is why Beta-Bernoulli models are so common in A/B testing and bandit algorithms. Understanding conjugacy lets you derive closed-form updates without numerical integration.

import numpy as np
from scipy import stats

def bayesian_coin_update(prior_alpha, prior_beta, heads, total_flips):
    """Beta-Bernoulli conjugate update."""
    posterior_alpha = prior_alpha + heads
    posterior_beta = prior_beta + (total_flips - heads)
    return posterior_alpha, posterior_beta

# Prior: Beta(2,2) — weak belief in fair coin
prior_a, prior_b = 2, 2
print(f"Prior mean: {prior_a/(prior_a+prior_b):.4f}")

# Observe 8 heads in 10 flips
heads, n = 8, 10
post_a, post_b = bayesian_coin_update(prior_a, prior_b, heads, n)
print(f"Posterior mean: {post_a/(post_a+post_b):.4f}")
print(f"Posterior 95% CI: {stats.beta.interval(0.95, post_a, post_b)}")

# Predictive probability of heads on next flip
predictive_mean = post_a / (post_a + post_b)
print(f"Predictive P(heads) = {predictive_mean:.4f}")

# Compare with MLE
mle = heads / n
print(f"MLE: {mle:.4f}")

Maximum Likelihood and Maximum A Posteriori Estimation

Maximum Likelihood Estimation (MLE) and Maximum A Posteriori (MAP) estimation are the two workhorses for parameter estimation in probabilistic models. MLE finds the parameter values that maximize the likelihood of the observed data: θ̂_MLE = argmax_θ P(X|θ). For a dataset of i.i.d. samples, this becomes argmax_θ ∏_i P(x_i|θ), often computed as the sum of log-likelihoods for numerical stability. In practice, MLE is straightforward: for a Bernoulli distribution, the MLE of the success probability is simply the empirical frequency of successes. For a Gaussian, the MLE of the mean is the sample mean, and the variance is the biased sample variance (dividing by n, not n-1). MLE is consistent and asymptotically efficient under regularity conditions, but it can overfit with small data because it treats the parameters as fixed unknowns.

MAP estimation injects a prior distribution over the parameters, turning the optimization into θ̂_MAP = argmax_θ P(θ|X) = argmax_θ P(X|θ)P(θ). The prior acts as a regularizer. For example, a Beta prior on a Bernoulli probability yields a posterior that is a Beta distribution with parameters updated by the counts of successes and failures. The MAP estimate is then (α + successes - 1) / (α + β + n - 2) for a Beta(α,β) prior, which shrinks the estimate toward the prior mean when data is scarce. As sample size grows, the prior's influence diminishes and MAP converges to MLE. In deep learning, weight decay (L2 regularization) is equivalent to a Gaussian prior on weights under MAP estimation.

Choosing between MLE and MAP depends on your data size and prior knowledge. With millions of samples, MLE is fine and computationally cheaper. With hundreds of samples or high-dimensional parameters, MAP prevents degenerate solutions. In practice, many production systems use MAP with weak priors (e.g., Laplace smoothing in Naive Bayes) as a safety net. Never assume MLE will generalize—always validate on held-out data. The bias-variance tradeoff is real: MLE has lower bias but higher variance; MAP introduces bias but reduces variance, often improving out-of-sample performance.

import numpy as np
from scipy.stats import beta

# Simulate biased coin flips: true p = 0.7
np.random.seed(42)
n = 20
p_true = 0.7
samples = np.random.binomial(1, p_true, size=n)
successes = samples.sum()
failures = n - successes

# MLE
p_mle = successes / n
print(f"MLE estimate: {p_mle:.3f}")

# MAP with Beta(2,2) prior (weakly informative, centered at 0.5)
alpha_prior, beta_prior = 2, 2
alpha_post = alpha_prior + successes
beta_post = beta_prior + failures
p_map = (alpha_post - 1) / (alpha_post + beta_post - 2)
print(f"MAP estimate (Beta(2,2) prior): {p_map:.3f}")

# Compare with posterior mean (Bayesian estimate)
p_bayes = alpha_post / (alpha_post + beta_post)
print(f"Posterior mean: {p_bayes:.3f}")
print(f"True p: {p_true}")

The Central Limit Theorem and Why It Matters in ML

The Central Limit Theorem (CLT) states that the sum (or average) of a large number of independent, identically distributed random variables, each with finite mean μ and variance σ², is approximately normally distributed, regardless of the original distribution's shape. Formally, let X₁, X₂, ..., Xₙ be i.i.d. with mean μ and variance σ². Then (1/√n) Σ (Xᵢ - μ) converges in distribution to N(0, σ²). In practice, for n ≥ 30, the sample mean X̄ is approximately N(μ, σ²/n). This is not a theorem about convergence of the data itself—it's about the sampling distribution of the mean.

In ML, the CLT justifies the use of Gaussian assumptions in many statistical methods. For example, confidence intervals for model accuracy rely on the CLT: if you have n test samples, the accuracy is approximately normal with variance p(1-p)/n, where p is the true accuracy. This lets you compute error bars. Gradient-based optimization also benefits: the stochastic gradient is an average over a mini-batch, and the CLT ensures that the gradient noise is approximately Gaussian, which is why momentum and Adam work well. The CLT also underpins the asymptotic normality of MLE estimates, enabling hypothesis testing and model comparison via likelihood ratio tests.

However, the CLT has assumptions that are often violated in practice. Heavy-tailed distributions (e.g., power laws in click-through rates) converge slowly—you may need thousands of samples, not 30. Non-i.i.d. data (e.g., time series with autocorrelation) invalidates the independence assumption. In such cases, the sample mean may not be normal, and confidence intervals will be misleading. Always check the empirical distribution of your estimator via bootstrapping rather than blindly applying CLT-based formulas.

In production, the CLT is used to set monitoring thresholds. For example, if you track the average prediction confidence over a sliding window, the CLT tells you that the distribution of that average is approximately normal. You can then set an alert when the average deviates by more than 3 standard deviations from the historical mean. This is a simple but effective drift detection method. Just remember: the CLT is an asymptotic result—finite-sample behavior can be poor, especially with skewed or multimodal distributions.

import numpy as np
import matplotlib.pyplot as plt

# Simulate CLT: sample means from exponential distribution (skewed)
np.random.seed(42)
lam = 0.5  # rate parameter, mean = 1/lam = 2
n_samples = 1000
sample_size = 50

# Generate many sample means
sample_means = []
for _ in range(n_samples):
    data = np.random.exponential(scale=1/lam, size=sample_size)
    sample_means.append(np.mean(data))

sample_means = np.array(sample_means)
print(f"Mean of sample means: {np.mean(sample_means):.3f} (expected: {1/lam})")
print(f"Std of sample means: {np.std(sample_means):.3f} (expected: {1/(lam*np.sqrt(sample_size)):.3f})")

# Check normality via histogram (code omitted for brevity, but would show bell shape)

Model Calibration: From Softmax to True Probabilities

Model calibration measures how well a model's predicted probabilities match the true empirical frequencies. A perfectly calibrated model means that among all predictions with confidence p, approximately p fraction are correct. For example, if a model predicts 0.8 for 100 samples, roughly 80 should be true positives. In practice, neural networks with softmax outputs are often miscalibrated—they tend to be overconfident, especially with modern architectures like deep ResNets. This is a known issue: temperature scaling, a simple post-processing method, can dramatically improve calibration.

Temperature scaling divides the logits by a scalar T > 0 before applying softmax: softmax(z/T). T is learned by minimizing the negative log-likelihood on a held-out validation set. T > 1 softens the probabilities (reducing overconfidence), T < 1 sharpens them. Temperature scaling does not change the model's accuracy—it only adjusts confidence. Other calibration methods include Platt scaling (logistic regression on logits) and isotonic regression (non-parametric). For multiclass problems, temperature scaling is the most common due to its simplicity and effectiveness.

Evaluation of calibration is done via reliability diagrams (plots of predicted confidence vs. empirical accuracy) and metrics like Expected Calibration Error (ECE). ECE bins predictions by confidence (e.g., 10 bins) and computes the average absolute difference between accuracy and confidence within each bin. A well-calibrated model has ECE near 0. In production, poor calibration can lead to bad decisions: a fraud detection system that outputs 0.99 confidence but only 80% accuracy will cause unnecessary false positives. Always calibrate if your downstream system uses probability thresholds.

Calibration is not just for classification. Regression models can be calibrated by checking that the predicted quantiles match empirical quantiles (e.g., in quantile regression). For probabilistic forecasts (e.g., weather prediction), the Probability Integral Transform (PIT) histogram is the standard diagnostic. In ML pipelines, calibrate after training but before deployment, and re-calibrate if data distribution shifts. Remember: calibration does not fix discrimination—a model can be perfectly calibrated but still biased. Use separate fairness metrics.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
from sklearn.calibration import calibration_curve

# Generate synthetic data
X, y = make_classification(n_samples=1000, n_features=10, random_state=42)
X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.3, random_state=42)

# Train a model (use a simple one for demo)
model = LogisticRegression()
model.fit(X_train, y_train)

# Get raw probabilities
probs = model.predict_proba(X_val)[:, 1]

# Compute calibration curve
prob_true, prob_pred = calibration_curve(y_val, probs, n_bins=10)
print("Bin centers (predicted):", prob_pred)
print("Empirical accuracy:", prob_true)

# Compute ECE
ece = np.mean(np.abs(prob_true - prob_pred))
print(f"Expected Calibration Error: {ece:.4f}")

# Temperature scaling (simple grid search)
from scipy.optimize import minimize

def nll_temperature(T, logits, y_true):
    scaled = logits / T
    probs = 1 / (1 + np.exp(-scaled))
    eps = 1e-12
    return -np.mean(y_true * np.log(probs + eps) + (1 - y_true) * np.log(1 - probs + eps))

logits = model.decision_function(X_val)
res = minimize(nll_temperature, x0=1.0, args=(logits, y_val), method='L-BFGS-B')
T_opt = res.x[0]
print(f"Optimal temperature: {T_opt:.3f}")

# Recompute calibration after scaling
probs_scaled = 1 / (1 + np.exp(-logits / T_opt))
prob_true_s, prob_pred_s = calibration_curve(y_val, probs_scaled, n_bins=10)
ece_scaled = np.mean(np.abs(prob_true_s - prob_pred_s))
print(f"ECE after temperature scaling: {ece_scaled:.4f}")

Production Pitfalls: Data Drift, Rare Events, and Debugging

Data drift occurs when the statistical properties of the input data change between training and production. This can be covariate drift (P(X) changes), label drift (P(Y) changes), or concept drift (P(Y|X) changes). In practice, covariate drift is the most common: user behavior shifts seasonally, new features are introduced, or data pipelines break silently. For example, a model trained on click data from 2023 may fail in 2024 if user demographics shift. Detection methods include monitoring feature distributions via KS tests, population stability index (PSI), or simply tracking prediction confidence averages over time.

Rare events (e.g., fraud, equipment failure) pose a different challenge. With class imbalance, models often predict the majority class, leading to low recall. Standard remedies like oversampling, SMOTE, or cost-sensitive learning help, but they introduce bias. In production, rare events are often non-stationary—fraud patterns evolve. A model trained on last year's fraud data will miss new tactics. The key is to treat rare event detection as an anomaly detection problem: use unsupervised methods (e.g., isolation forests) or train a separate model on the minority class. Always monitor precision-recall curves, not accuracy, because accuracy is misleading with imbalance.

Debugging probabilistic models in production requires a systematic approach. Start with input validation: check for NaN, infinite values, or out-of-range features. Then examine prediction distributions: are they too sharp (overconfident) or too flat (underconfident)? Use calibration curves and ECE. If a model's average confidence drops suddenly, it could indicate drift or a data pipeline bug. Log all predictions with their features and ground truth (when available) to enable root cause analysis. Tools like SHAP or LIME can help explain individual predictions, but they are slow for real-time debugging—use them offline on sampled data.

A common production pitfall is the 'silent failure' where the model still outputs numbers but they are meaningless. For example, a missing feature encoded as 0 might cause the model to always predict a certain class. Implement automated monitoring: track prediction entropy, confidence, and feature statistics. Set alerts for when these metrics deviate by more than 3 standard deviations from the training baseline. Also, have a fallback model (e.g., a simpler logistic regression) that can be deployed quickly if the primary model degrades. Remember: in production, the model is only as good as the data pipeline feeding it.

import numpy as np
from scipy.stats import ks_2samp

# Simulate training and production feature distributions
np.random.seed(42)
train_feature = np.random.normal(loc=0.0, scale=1.0, size=1000)
prod_feature = np.random.normal(loc=0.5, scale=1.2, size=1000)  # drift!

# Kolmogorov-Smirnov test
stat, p_value = ks_2samp(train_feature, prod_feature)
print(f"KS statistic: {stat:.3f}, p-value: {p_value:.5f}")
if p_value < 0.05:
    print("ALERT: Significant drift detected (p < 0.05)")
else:
    print("No significant drift")

# Population Stability Index (PSI)
def psi(expected, actual, bins=10):
    # Simple PSI calculation
    expected_counts, _ = np.histogram(expected, bins=bins, range=(-3,3))
    actual_counts, _ = np.histogram(actual, bins=bins, range=(-3,3))
    expected_pct = expected_counts / expected_counts.sum() + 1e-6
    actual_pct = actual_counts / actual_counts.sum() + 1e-6
    return np.sum((actual_pct - expected_pct) * np.log(actual_pct / expected_pct))

psi_val = psi(train_feature, prod_feature)
print(f"PSI: {psi_val:.3f} (threshold: >0.1 indicates drift)")