Probabilistic Regression Model

26 min read
Updated
Suggest editHistoryTalk
RawGraph

Last edited

Fact-checked

In review queue

Sources

20 citations

Revision

v6 · 5,205 words

Fact-checks are independent of edits: a reviewer re-verifies the article against its sources and stamps the date. How we verify

A probabilistic regression model (also called distributional regression) is a regression model that outputs a full probability distribution over possible target values rather than a single point estimate. By modeling the conditional distribution p(yx)p(y \mid x) of a continuous target variable y given input features x, these models capture the uncertainty inherent in predictions. A point prediction states what value a model expects; only a predicted distribution also states how confident the model is in that value. This makes them especially valuable in applications where understanding the range and likelihood of outcomes is as important as the prediction itself, such as in medical diagnosis, weather forecasting, financial risk assessment, and engineering safety analysis.

Unlike standard regression methods (such as ordinary least squares) that return only a best-guess value, probabilistic regression models provide prediction intervals, variance estimates, or entire density functions. This additional information allows practitioners to make risk-aware decisions and to quantify how confident a model is in each of its predictions.

Explain like I'm 5 (ELI5)

Imagine you want to guess how tall a sunflower will grow. A regular guess gives you one number, like "6 feet." But a probabilistic guess is more like saying, "It will probably be between 5 and 7 feet tall, and most likely around 6 feet." The guess includes a range, and it tells you how sure you are. If you have seen lots of sunflowers before, your range is small because you know what to expect. If you have never grown one, your range is bigger because you are less certain. That is what a probabilistic regression model does: it gives not just a single answer but a whole picture of what could happen and how likely each outcome is.

How does a probabilistic regression model differ from standard regression?

Standard (point) regression, such as ordinary least squares, returns one number for each input: a single best estimate of the target. A probabilistic regression model returns a distribution, from which a point estimate (the mean or median), a variance, a prediction interval, or the entire density can all be read off. The two approaches also differ in what they optimize and in what a practitioner can do with the output.

AspectStandard (point) regressionProbabilistic regression
Output for each inputA single value (point estimate)A full probability distribution over the target
Typical training objectiveSquared or absolute errorNegative log-likelihood, pinball loss, or a proper scoring rule
Uncertainty informationNoneVariance, prediction intervals, quantiles, or the full density
Downstream useOne fixed decisionRisk-aware decisions, interval-based planning, active learning
Example question answered"What is the predicted value?""How likely is the value to exceed a threshold?"

Because a probabilistic model reports both a central estimate and its uncertainty, it can flag when a prediction is unreliable, for example in regions of the input space with little training data. This property is the main reason probabilistic regression is preferred in safety-critical settings.

Mathematical foundation

In a standard regression setting, we model the expected value of y given x:

y^=f(x)\hat{y} = f(x)

A probabilistic regression model instead estimates the full conditional distribution:

p(yx,θ)p(y \mid x, \theta)

where θ\theta represents the model parameters. For many probabilistic regression methods, the output distribution is assumed to be Gaussian:

p(yx)=N(y;μ(x),σ2(x))p(y \mid x) = \mathcal{N}(y; \mu(x), \sigma^2(x))

Here, both the mean μ(x)\mu(x) and the variance σ2(x)\sigma^2(x) are functions of the input. The model is typically trained by maximizing the log-likelihood:

L(θ)=i=1Nlogp(yixi,θ)\mathcal{L}(\theta) = \sum_{i=1}^{N} \log p(y_i \mid x_i, \theta)

For a Gaussian output, the negative log-likelihood reduces to:

L(θ)=12i=1N[logσ2(xi)+(yiμ(xi))2σ2(xi)]+const-\mathcal{L}(\theta) = \frac{1}{2} \sum_{i=1}^{N} \left[ \log \sigma^2(x_i) + \frac{(y_i - \mu(x_i))^2}{\sigma^2(x_i)} \right] + \text{const}

This formulation naturally penalizes both inaccurate mean predictions and poorly calibrated variance estimates. Maximizing the log-likelihood is equivalent to optimizing the logarithmic score, one of the proper scoring rules discussed below.

What types of uncertainty do these models capture?

Probabilistic regression models distinguish between two fundamental sources of uncertainty, a taxonomy formalized for machine learning by Hullermeier and Waegeman (2021) [13]:

TypeAlso known asSourceReducible?Example
Aleatoric uncertaintyData uncertainty, statistical uncertaintyInherent noise in the data-generating processNoMeasurement sensor noise, natural variability in patient outcomes
Epistemic uncertaintyModel uncertainty, knowledge uncertaintyLack of knowledge due to limited training dataYes (with more data)Predictions in regions of feature space with few training examples

Aleatoric uncertainty is irreducible and stems from the randomness inherent in the system being modeled. For example, even with perfect knowledge of all relevant factors, there is still random variation in daily stock prices. Aleatoric uncertainty can be further divided into homoscedastic uncertainty (constant noise across the input space) and heteroscedastic uncertainty (noise that varies with input).

Epistemic uncertainty arises from incomplete knowledge and can be reduced by collecting more training data or by improving the model. Bayesian methods naturally capture epistemic uncertainty through posterior distributions over model parameters. In regions where few training examples exist, the posterior is broad, signaling high epistemic uncertainty.

The terminology reflects these two sources: the word aleatoric derives from the Latin alea, meaning dice or chance, while epistemic derives from the Greek episteme, meaning knowledge. Distinguishing the two matters in practice because only epistemic uncertainty can be driven down with more data, which is why it guides strategies such as active learning and experimental design [13].

A well-designed probabilistic regression model should ideally capture both types. Total predictive uncertainty is often expressed as the sum of aleatoric and epistemic components:

Var[yx]=E[σ2(x)]aleatoric+Var[μ(x)]epistemic\text{Var}[y \mid x] = \underbrace{\mathbb{E}[\sigma^2(x)]}_{\text{aleatoric}} + \underbrace{\text{Var}[\mu(x)]}_{\text{epistemic}}

Classical probabilistic regression methods

Bayesian linear regression

Bayesian linear regression extends classical linear regression by treating the regression coefficients as random variables with prior distributions rather than fixed unknown parameters. Given a linear model y=Xβ+ϵy = X\beta + \epsilon with Gaussian noise ϵN(0,σ2I)\epsilon \sim \mathcal{N}(0, \sigma^2 I), a Gaussian prior is placed on the weights:

βN(μ0,Λ01)\beta \sim \mathcal{N}(\mu_0, \Lambda_0^{-1})

Using Bayes' theorem, the posterior distribution over the weights given observed data (X, y) is:

p(βX,y)=N(μn,Λn1)p(\beta \mid X, y) = \mathcal{N}(\mu_n, \Lambda_n^{-1})

where:

  • Λn=XX/σ2+Λ0\Lambda_n = X^\top X / \sigma^2 + \Lambda_0 (posterior precision)
  • μn=Λn1(Xy/σ2+Λ0μ0)\mu_n = \Lambda_n^{-1}(X^\top y / \sigma^2 + \Lambda_0 \mu_0) (posterior mean)

The posterior mean is a weighted combination of the prior mean and the maximum likelihood estimate, with the weighting determined by the relative precision of each. The posterior predictive distribution for a new input xx^* is:

p(yx,X,y)=N(xTμn,  σ2+xTΛn1x)p(y^* \mid x^*, X, y) = \mathcal{N}(x^{*T} \mu_n, \; \sigma^2 + x^{*T} \Lambda_n^{-1} x^*)

This prediction captures two sources of variance: the noise variance σ2\sigma^2 (aleatoric) and the parameter uncertainty xΛn1xx^{*\top} \Lambda_n^{-1} x^* (epistemic). As more data is observed, Λn\Lambda_n grows and the epistemic component shrinks.

Bayesian linear regression also enables principled model comparison through the marginal likelihood (model evidence), which can be computed in closed form for conjugate priors.

Gaussian process regression

Gaussian process (GP) regression, also known as kriging, is a nonparametric Bayesian approach that places a prior distribution directly over functions rather than over a finite set of parameters. The method originated in geostatistics: the South African mining engineer Danie G. Krige introduced it empirically in 1951 to improve ore-reserve estimation on the Witwatersrand [19], and the French mathematician Georges Matheron gave it a formal theoretical basis in 1963 through his theory of regionalized variables [20]. The standard machine learning reference is the textbook by Rasmussen and Williams (2006) [4]. A Gaussian process is fully specified by a mean function m(x)m(x) and a covariance (kernel) function k(x,x)k(x, x'):

f(x)GP(m(x),k(x,x))f(x) \sim \mathcal{GP}(m(x), k(x, x'))

Given n training observations (X, y) and assuming Gaussian noise with variance σ2\sigma^2, the predictive distribution at a new test point xx^* is Gaussian:

p(fx,X,y)=N(fˉ,cov(f))p(f^* \mid x^*, X, y) = \mathcal{N}(\bar{f}^*, \text{cov}(f^*))

with:

  • Posterior mean: fˉ=k(x,X)[k(X,X)+σ2I]1y\bar{f}^* = k(x^*, X)[k(X, X) + \sigma^2 I]^{-1} y
  • Posterior variance: cov(f)=k(x,x)k(x,X)[k(X,X)+σ2I]1k(X,x)\text{cov}(f^*) = k(x^*, x^*) - k(x^*, X)[k(X, X) + \sigma^2 I]^{-1} k(X, x^*)

The posterior mean acts as a weighted combination of the observed outputs, while the posterior variance depends only on the input locations and decreases in regions near the training data.

Common kernel functions

The choice of kernel encodes assumptions about the underlying function, such as smoothness, periodicity, and length scale.

KernelFormulaProperties
Squared Exponential (RBF)k(x,x)=σf2exp(xx2/22)k(x, x') = \sigma_f^2 \exp(-\lVert x - x' \rVert^2 / 2\ell^2)Infinitely differentiable; very smooth functions
Matern (ν=3/2\nu = 3/2)k(x,x)=σf2(1+3d/)exp(3d/)k(x, x') = \sigma_f^2 (1 + \sqrt{3}\, d/\ell) \exp(-\sqrt{3}\, d/\ell)Once differentiable; realistic for many physical processes
Matern (ν=5/2\nu = 5/2)k(x,x)=σf2(1+5d/+5d2/32)exp(5d/)k(x, x') = \sigma_f^2 (1 + \sqrt{5}\, d/\ell + 5d^2/3\ell^2) \exp(-\sqrt{5}\, d/\ell)Twice differentiable; good balance of smoothness
Periodick(x,x)=σf2exp(2sin2(πxx/p)/2)k(x, x') = \sigma_f^2 \exp(-2 \sin^2(\pi \lVert x - x' \rVert / p) / \ell^2)Captures periodic patterns with period p
Lineark(x,x)=σb2+σv2(xc)(xc)k(x, x') = \sigma_b^2 + \sigma_v^2 (x - c)(x' - c)Equivalent to Bayesian linear regression
Rational Quadratick(x,x)=σf2(1+xx2/2α2)αk(x, x') = \sigma_f^2 (1 + \lVert x - x' \rVert^2 / 2\alpha\ell^2)^{-\alpha}Mixture of RBF kernels with different length scales

Kernel hyperparameters (such as the length scale \ell and signal variance σf2\sigma_f^2) are typically optimized by maximizing the log marginal likelihood. Because this objective can be non-convex, multiple restarts of the optimizer are often used to avoid local optima.

A major limitation of standard GP regression is its computational cost: training requires O(n3)O(n^3) time and O(n2)O(n^2) memory due to the inversion of the n-by-n covariance matrix. This has motivated the development of sparse GP approximations, including inducing point methods and the Vecchia approximation, which reduce complexity to O(nm2)O(nm^2) where m is the number of inducing points (mnm \ll n).

Generalized linear models

Generalized linear models (GLMs) extend linear regression by allowing the response variable to follow any distribution from the exponential family (such as Gaussian, Poisson, Gamma, or Inverse-Gaussian). A GLM has three components:

  1. A probability distribution for y from the exponential family
  2. A linear predictor η=Xβ\eta = X\beta
  3. A link function gg that relates the mean of y to the linear predictor: g(μ)=ηg(\mu) = \eta

For probabilistic regression, GLMs with appropriate link functions and distributions can model non-Gaussian response variables while still providing distributional predictions. For instance, a Gamma GLM with a log link is commonly used for modeling strictly positive, right-skewed continuous outcomes.

Quantile regression

Quantile regression, introduced by Koenker and Bassett in 1978 [1], estimates conditional quantiles of the response variable rather than the conditional mean. For a given quantile level τ\tau (where 0<τ<10 < \tau < 1), the τ\tau-th conditional quantile function qτ(x)q_\tau(x) is estimated by minimizing the pinball loss (also called the check loss):

Lτ(y,q^)={τ(yq^)if yq^(1τ)(q^y)if y<q^L_\tau(y, \hat{q}) = \begin{cases} \tau (y - \hat{q}) & \text{if } y \geq \hat{q} \\ (1 - \tau)(\hat{q} - y) & \text{if } y < \hat{q} \end{cases}

By fitting separate quantile functions for multiple values of τ\tau (for example, τ=0.05,0.5,0.95\tau = 0.05, 0.5, 0.95), quantile regression constructs prediction intervals without assuming any parametric distribution for the errors. This makes it particularly useful for heteroscedastic data, where the spread of y varies with x. Unlike methods that assume Gaussian errors, quantile regression is robust to outliers and non-normal error distributions.

The special case of τ=0.5\tau = 0.5 yields the conditional median, which minimizes the sum of absolute deviations rather than squared errors.

Deep learning approaches

Heteroscedastic neural networks

Nix and Weigend (1994) proposed one of the earliest neural network approaches to probabilistic regression [2]. The network outputs two values for each input: a predicted mean μ(x)\mu(x) and a predicted variance σ2(x)\sigma^2(x). The model is trained by minimizing the Gaussian negative log-likelihood:

L=12Ni=1N[logσ2(xi)+(yiμ(xi))2σ2(xi)]\mathcal{L} = \frac{1}{2N} \sum_{i=1}^{N} \left[ \log \sigma^2(x_i) + \frac{(y_i - \mu(x_i))^2}{\sigma^2(x_i)} \right]

This allows the model to learn input-dependent (heteroscedastic) noise, producing wider prediction intervals where the data is noisier. The variance output is typically parameterized through a softplus or exponential activation to ensure positivity.

One known issue with this approach is that the variance predictions can become overconfident, especially in regions with limited data. Regularization techniques and careful initialization help mitigate this problem.

Mixture density networks

Mixture density networks (MDNs), introduced by Christopher Bishop in 1994 [3], combine a neural network with a Gaussian mixture model. The network outputs the parameters of a mixture of K Gaussians:

p(yx)=k=1Kπk(x)N(y;μk(x),σk2(x))p(y \mid x) = \sum_{k=1}^{K} \pi_k(x) \, \mathcal{N}(y; \mu_k(x), \sigma_k^2(x))

where πk(x)\pi_k(x) are the mixing coefficients (summing to 1), and μk(x)\mu_k(x) and σk2(x)\sigma_k^2(x) are the mean and variance of the k-th component. MDNs can represent multimodal predictive distributions, which is especially important for inverse problems and situations where multiple valid outputs exist for a single input. For example, in robotics, a single target position might be reachable through multiple joint configurations.

The network is trained end-to-end using maximum likelihood, and the mixing coefficients are typically produced via a softmax layer.

Monte Carlo dropout

Gal and Ghahramani (2016) showed that dropout, a regularization technique commonly used during training, can be reinterpreted as approximate Bayesian inference [7]. Their theoretical framework casts "dropout training in deep neural networks (NNs) as approximate Bayesian inference in deep Gaussian processes" [7]. By keeping dropout active during prediction and performing T stochastic forward passes, one obtains a set of predictions {y^1,,y^T}\{\hat{y}_1, \ldots, \hat{y}_T\} for each input. The predictive mean and variance are estimated as:

μ^(x)=1Tt=1Ty^t(x),σ^2(x)=1Tt=1T(y^t(x)μ^(x))2\hat{\mu}(x) = \frac{1}{T} \sum_{t=1}^{T} \hat{y}_t(x), \qquad \hat{\sigma}^2(x) = \frac{1}{T} \sum_{t=1}^{T} (\hat{y}_t(x) - \hat{\mu}(x))^2

This method is attractive because it requires no changes to the network architecture and can be applied to any model that already uses dropout. However, the quality of the uncertainty estimates depends on the dropout rate and may not fully capture the true posterior, particularly during extrapolation.

Deep ensembles

Lakshminarayanan, Pritzel, and Blundell (2017) proposed deep ensembles as a simple and scalable alternative to Bayesian neural networks [8]. The authors described their approach as "simple to implement, readily parallelizable, requires very little hyperparameter tuning, and yields high quality predictive uncertainty estimates" [8]. An ensemble of M independently trained networks (each with different random initializations) produces predictions that are combined as a mixture of Gaussians:

p(yx)=1Mm=1MN(y;μm(x),σm2(x))p(y \mid x) = \frac{1}{M} \sum_{m=1}^{M} \mathcal{N}(y; \mu_m(x), \sigma_m^2(x))

The predictive mean and total variance are:

μ(x)=1Mm=1Mμm(x)\mu^*(x) = \frac{1}{M} \sum_{m=1}^{M} \mu_m(x)

σ2(x)=1Mm=1M[σm2(x)+μm2(x)]μ2(x)\sigma^{*2}(x) = \frac{1}{M} \sum_{m=1}^{M} [\sigma_m^2(x) + \mu_m^2(x)] - \mu^{*2}(x)

In practice, a deep ensemble uses only a small number of independently trained networks, commonly around five, which already delivers strong calibration. Deep ensembles have been shown to produce well-calibrated uncertainty estimates and have become a widely used baseline for predictive uncertainty in deep learning [8]. They tend to outperform Monte Carlo dropout on out-of-distribution detection and overall calibration; the original study specifically evaluated predictive uncertainty on test examples from both known and unknown distributions and showed higher uncertainty on out-of-distribution inputs [8]. The main drawback is the increased computational cost of training and storing multiple models.

Bayesian neural networks

Bayesian neural networks (BNNs) place probability distributions over the network weights rather than learning fixed point estimates. The posterior over weights p(wD)p(w \mid D) is combined with the likelihood to produce predictive distributions:

p(yx,D)=p(yx,w)p(wD)dwp(y \mid x, D) = \int p(y \mid x, w) \, p(w \mid D) \, dw

Since this integral is intractable for all but the simplest networks, approximate inference methods are used:

  • Variational inference: Graves (2011) [5] and Blundell et al. (2015, "Bayes by Backprop") [6] approximate the posterior with a parameterized distribution q(w;ϕ)q(w; \phi) by minimizing the KL divergence between q and the true posterior. The variational parameters ϕ\phi are learned alongside the model through backpropagation.
  • Markov chain Monte Carlo (MCMC): Methods such as stochastic gradient Langevin dynamics (SGLD) and Hamiltonian Monte Carlo (HMC) sample from the posterior. These methods are asymptotically exact but computationally expensive for large networks.

BNNs naturally decompose uncertainty into aleatoric and epistemic components and provide principled uncertainty estimates, but they remain more difficult to train and scale than ensembles or dropout-based methods.

Deep evidential regression

Amini et al. (2020) proposed deep evidential regression, which places an evidential prior (Normal Inverse-Gamma distribution) over the parameters of the Gaussian likelihood [12]. The network outputs four parameters (γ,ν,α,β)(\gamma, \nu, \alpha, \beta) that parameterize this higher-order distribution, allowing simultaneous estimation of aleatoric and epistemic uncertainty in a single forward pass without sampling. The authors emphasize that the method "does not rely on sampling during inference or on out-of-distribution (OOD) examples for training, thus enabling efficient and scalable uncertainty learning" [12].

The predictive distribution has mean γ\gamma and variance β(1+1/ν)/(α1)\beta(1 + 1/\nu) / (\alpha - 1). Epistemic uncertainty is captured through ν\nu (inversely), while aleatoric uncertainty is captured through β/α\beta/\alpha. This approach is computationally efficient since it avoids the need for multiple forward passes or ensemble training.

Gradient boosting approaches

NGBoost

NGBoost (Natural Gradient Boosting), developed by Duan et al. (2020) at Stanford, extends gradient boosting to probabilistic prediction [11]. Traditional gradient boosting optimizes a scalar loss function, while NGBoost treats the parameters of a conditional probability distribution as the targets of a multiparameter boosting algorithm.

The key insight is the use of the natural gradient (rather than the ordinary gradient) to update the distribution parameters at each boosting iteration. As the authors put it, "the Natural Gradient is required to correct the training dynamics of our multiparameter boosting approach" [11]. The natural gradient accounts for the geometry of the parameter space of probability distributions, leading to more stable and efficient optimization. NGBoost is flexible in that it can be used with any base learner, any parametric distribution family, and any proper scoring rule [11].

CatBoost and LightGBM

Other popular gradient boosting frameworks, such as CatBoost and LightGBM, also support quantile regression objectives, enabling probabilistic predictions through estimated conditional quantiles. These approaches do not model a full parametric distribution but provide prediction intervals by fitting separate models for different quantile levels.

Distribution-free methods

Conformal prediction

Conformal prediction, pioneered by Vladimir Vovk and collaborators and given a book-length treatment by Vovk, Gammerman, and Shafer in 2005 [16], is a framework for constructing prediction intervals with finite-sample coverage guarantees, regardless of the underlying data distribution or model architecture. Given a desired coverage level 1α1 - \alpha (for example, 90%), conformal prediction produces intervals C(x)C(x) such that:

P(yC(x))1αP(y \in C(x)) \geq 1 - \alpha

The method works by computing nonconformity scores on a held-out calibration set and using the empirical quantile of these scores to set the interval width. The only assumption required is exchangeability of the data (a weaker condition than the i.i.d. assumption). As Angelopoulos and Bates (2023) summarize, conformal sets "are valid in a distribution-free sense: they possess explicit, non-asymptotic guarantees even without distributional assumptions or model assumptions" [14].

Conformalized quantile regression (CQR), proposed by Romano, Patterson, and Candes (2019), combines quantile regression with conformal calibration [10]. The base quantile regression model produces initial interval estimates, and conformal prediction adjusts these intervals to guarantee marginal coverage. CQR intervals can adapt to heteroscedastic noise because the underlying quantile regression model produces intervals of varying width.

Conformal prediction is model-agnostic, computationally efficient (no retraining needed), and provides valid coverage guarantees even when the model is misspecified. However, it guarantees only marginal coverage (averaged over the test distribution) rather than conditional coverage for each individual input.

How are probabilistic regression models evaluated?

Evaluating probabilistic regression models requires metrics that assess the quality of the entire predicted distribution, not just the point prediction. A proper scoring rule is a function S(F,y)S(F, y) that assigns a score to a predictive distribution F when the true outcome is y, and for which the expected score is optimized when the forecaster reports the true distribution [9]. A scoring rule is strictly proper if the optimum is unique.

Scoring ruleFormula (simplified)Properties
Log-likelihood (logarithmic score)S(F,y)=logf(y)S(F, y) = \log f(y)Strictly proper; local (depends only on density at y); sensitive to tail behavior
Continuous Ranked Probability Score (CRPS)CRPS(F,y)=(F(z)1{yz})2dz\mathrm{CRPS}(F, y) = \int (F(z) - \mathbf{1}\{y \le z\})^2 \, dzStrictly proper; generalizes MAE to distributional forecasts; robust to outliers
Interval ScoreISα=(ul)+(2/α)(ly)1{y<l}+(2/α)(yu)1{y>u}\mathrm{IS}_\alpha = (u - l) + (2/\alpha)(l - y)\mathbf{1}\{y < l\} + (2/\alpha)(y - u)\mathbf{1}\{y > u\}Proper for prediction intervals at level (1α)(1 - \alpha); penalizes width and miscoverage
Energy ScoreES(F,y)=EYy(1/2)EYY\mathrm{ES}(F, y) = \mathbb{E}\lVert Y - y \rVert - (1/2)\mathbb{E}\lVert Y - Y' \rVertStrictly proper; multivariate generalization of CRPS

Calibration

A probabilistic regression model is calibrated if the predicted probability levels match the observed frequencies. For example, if a model predicts that 90% of outcomes will fall within a particular interval, then approximately 90% of observed outcomes should indeed fall within that interval. Calibration is typically assessed using:

  • Probability Integral Transform (PIT) histograms: If the model is well-calibrated, the PIT values (F(y)F(y) evaluated at observed y) should be uniformly distributed.
  • Calibration plots: These compare the nominal coverage level with the observed coverage across multiple probability thresholds.
  • Reliability diagrams: Similar to calibration plots but binned by predicted confidence level.

A model can have good calibration but poor sharpness (overly wide intervals), or vice versa. Gneiting, Balabdaoui, and Raftery (2007) framed the goal of probabilistic forecasting as "maximizing the sharpness of the predictive distributions subject to calibration" [15]. In other words, the ideal model is both sharp (narrow intervals) and calibrated (correct coverage), with calibration as a hard constraint and sharpness as the quantity to optimize within it.

How do the methods compare?

MethodDistributional assumptionUncertainty type capturedScalabilityMultimodal supportKey advantage
Bayesian linear regressionGaussian prior and likelihoodBothO(d3)O(d^3) in dimensionsNoClosed-form posterior
Gaussian process regressionPrior over functionsBothO(n3)O(n^3) trainingNoNonparametric flexibility
Quantile regressionDistribution-freeAleatoricO(n)O(n) per quantilePartial (via crossing quantiles)No distributional assumptions
Heteroscedastic NNGaussian outputAleatoric onlyGPU-scalableNoSimple to implement
Mixture density networkGaussian mixture outputAleatoricGPU-scalableYesModels multimodal targets
Monte Carlo dropoutApproximate BayesianBoth (approximate)Same as base modelNoNo architecture changes
Deep ensembleMixture of predictionsBothM times base costPartialStrong calibration
BNN (variational)Full weight posteriorBothGPU-scalable (approximate)Depends on likelihoodPrincipled uncertainty
Deep evidential regressionNormal Inverse-Gamma priorBothSingle forward passNoFast inference
NGBoostParametric (user-chosen)BothO(n)O(n) per iterationDepends on familyWorks with any base learner
Conformal predictionDistribution-freeCoverage guaranteeO(nlogn)O(n \log n) calibrationN/A (interval only)Finite-sample guarantee

What is probabilistic regression used for?

Probabilistic regression models are used across many domains where quantifying predictive uncertainty is important.

Weather and climate forecasting

Weather prediction is inherently uncertain due to the chaotic nature of atmospheric dynamics. Probabilistic weather forecasts provide ensemble predictions or distributional forecasts (for example, "there is a 70% chance that tomorrow's temperature will be between 18 and 24 degrees Celsius"). Modern numerical weather prediction systems use ensemble methods, post-processed with techniques such as quantile regression forests or Bayesian model averaging, to produce calibrated probabilistic forecasts. These forecasts support decision-making in agriculture, energy, aviation, and disaster preparedness.

Finance and risk management

In financial applications, probabilistic regression is used for value-at-risk (VaR) estimation, portfolio optimization, and option pricing. Quantile regression is especially common for estimating tail risks, where the focus is on the extreme percentiles of return distributions. Heteroscedastic models like GARCH (Generalized Autoregressive Conditional Heteroskedasticity), introduced by Bollerslev (1986) [17] as a generalization of Engle's (1982) ARCH model [18], capture the time-varying volatility of financial returns, providing distributional forecasts that inform risk management strategies.

Healthcare and medical diagnosis

Probabilistic models in healthcare support clinical decision-making by estimating uncertainty in patient outcomes. Growth charts, used in pediatrics, are a classic example of quantile regression applied to medical data. Survival analysis models, such as Cox proportional hazards with Bayesian extensions, provide probabilistic predictions of patient prognosis. Gaussian processes have been applied to model patient trajectories in intensive care, where principled uncertainty estimates help clinicians identify deteriorating patients.

Autonomous systems and robotics

Self-driving vehicles and robotic systems use probabilistic regression to predict the future trajectories of other agents (pedestrians, vehicles). Mixture density networks are well-suited to this task because trajectory prediction is inherently multimodal: a pedestrian at a crosswalk might continue walking, stop, or turn. Uncertainty-aware predictions enable safer planning and control decisions.

Engineering and materials science

Bayesian optimization, which relies on Gaussian process regression as a surrogate model, is widely used for experimental design and hyperparameter tuning. The GP posterior provides both a predicted value and an uncertainty estimate at each candidate point, and the acquisition function (such as expected improvement) uses both quantities to balance exploration and exploitation. In materials science, probabilistic models guide the search for materials with desired properties while minimizing costly experiments.

Energy systems

Probabilistic forecasting of renewable energy generation (solar irradiance, wind power) is critical for grid management and energy trading. Quantile regression and ensemble methods provide prediction intervals that help grid operators maintain supply-demand balance and plan reserve capacity under uncertainty.

What software and tools are available?

Several widely used libraries provide implementations of probabilistic regression methods.

LibraryLanguageMethods supported
scikit-learnPythonGaussian process regression, quantile regression (GradientBoostingRegressor)
GPyTorchPythonScalable Gaussian process regression with GPU support
PyTorchPythonHeteroscedastic NNs, MDNs, BNNs, deep ensembles (via custom code)
TensorFlow ProbabilityPythonBNNs, variational inference, mixture density layers
Stan / PyStanPython / RFull Bayesian regression with MCMC and variational inference
PyMCPythonBayesian regression, Gaussian processes, hierarchical models
NGBoostPythonNatural gradient boosting for probabilistic prediction
MAPIEPythonConformal prediction for regression and classification
brmsRBayesian regression with Gaussian process support
statsmodelsPythonQuantile regression, GLMs

History and key publications

The development of probabilistic regression has drawn from statistics, machine learning, and geostatistics over several decades.

YearContributionAuthors
1763Bayes' theorem formulated (published posthumously)Thomas Bayes
1951Kriging introduced for ore-reserve estimationDanie G. Krige
1963Kriging given a theoretical basis (theory of regionalized variables)Georges Matheron
1978Quantile regression introducedRoger Koenker, Gilbert Bassett
1982ARCH model for time-varying varianceRobert F. Engle
1986GARCH model introducedTim Bollerslev
1994Heteroscedastic neural network regression proposedDavid Nix, Andreas Weigend
1994Mixture density networks introducedChristopher Bishop
2005Conformal prediction (book-length treatment)Vladimir Vovk, Alexander Gammerman, Glenn Shafer
2006"Gaussian Processes for Machine Learning" textbook publishedCarl Edward Rasmussen, Christopher Williams
2007Sharpness-subject-to-calibration paradigm; strictly proper scoring rulesTilmann Gneiting, Fadoua Balabdaoui, Adrian Raftery
2011Practical variational inference for neural networksAlex Graves
2015Bayes by Backprop (weight uncertainty in neural networks)Charles Blundell et al.
2016Monte Carlo dropout as Bayesian approximationYarin Gal, Zoubin Ghahramani
2017Deep ensembles for uncertainty estimationBalaji Lakshminarayanan, Alexander Pritzel, Charles Blundell
2019Conformalized quantile regressionYaniv Romano, Evan Patterson, Emmanuel Candes
2020NGBoost: natural gradient boosting for probabilistic predictionTony Duan, Anand Avati, et al.
2020Deep evidential regressionAlexander Amini et al.

See also

References

  1. Koenker, R., & Bassett, G. (1978). "Regression Quantiles." *Econometrica*, 46(1), 33-50.
  2. Nix, D. A., & Weigend, A. S. (1994). "Estimating the mean and variance of the target probability distribution." *Proceedings of the IEEE International Conference on Neural Networks (ICNN'94)*, 1, 55-60.
  3. Bishop, C. M. (1994). "Mixture Density Networks." Technical Report NCRG/94/004, Aston University.
  4. Rasmussen, C. E., & Williams, C. K. I. (2006). *Gaussian Processes for Machine Learning*. MIT Press.
  5. Graves, A. (2011). "Practical Variational Inference for Neural Networks." *Advances in Neural Information Processing Systems*, 24.
  6. Blundell, C., Cornebise, J., Kavukcuoglu, K., & Wierstra, D. (2015). "Weight Uncertainty in Neural Networks." *Proceedings of the 32nd International Conference on Machine Learning (ICML)*.
  7. Gal, Y., & Ghahramani, Z. (2016). "Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning." *Proceedings of the 33rd International Conference on Machine Learning (ICML)*.
  8. Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017). "Simple and Scalable Predictive Uncertainty Estimation using Deep Ensembles." *Advances in Neural Information Processing Systems*, 30.
  9. Gneiting, T., & Raftery, A. E. (2007). "Strictly Proper Scoring Rules, Prediction, and Estimation." *Journal of the American Statistical Association*, 102(477), 359-378.
  10. Romano, Y., Patterson, E., & Candes, E. (2019). "Conformalized Quantile Regression." *Advances in Neural Information Processing Systems*, 32.
  11. Duan, T., Avati, A., Ding, D. Y., Thai, K. K., Basu, S., Ng, A., & Schuler, A. (2020). "NGBoost: Natural Gradient Boosting for Probabilistic Prediction." *Proceedings of the 37th International Conference on Machine Learning (ICML)*.
  12. Amini, A., Schwarting, W., Soleimany, A., & Rus, D. (2020). "Deep Evidential Regression." *Advances in Neural Information Processing Systems*, 33.
  13. Hullermeier, E., & Waegeman, W. (2021). "Aleatoric and Epistemic Uncertainty in Machine Learning: An Introduction to Concepts and Methods." *Machine Learning*, 110, 457-506.
  14. Angelopoulos, A. N., & Bates, S. (2023). "Conformal Prediction: A Gentle Introduction." *Foundations and Trends in Machine Learning*, 16(4), 494-591.
  15. Gneiting, T., Balabdaoui, F., & Raftery, A. E. (2007). "Probabilistic Forecasts, Calibration and Sharpness." *Journal of the Royal Statistical Society: Series B (Statistical Methodology)*, 69(2), 243-268.
  16. Vovk, V., Gammerman, A., & Shafer, G. (2005). *Algorithmic Learning in a Random World*. Springer.
  17. Bollerslev, T. (1986). "Generalized Autoregressive Conditional Heteroskedasticity." *Journal of Econometrics*, 31(3), 307-327.
  18. Engle, R. F. (1982). "Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation." *Econometrica*, 50(4), 987-1008.
  19. Krige, D. G. (1951). "A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand." *Journal of the Chemical, Metallurgical and Mining Society of South Africa*, 52(6), 119-139.
  20. Matheron, G. (1963). "Principles of Geostatistics." *Economic Geology*, 58(8), 1246-1266.

Improve this article

Add missing citations, update stale details, or suggest a clearer explanation. Every suggestion is reviewed for sourcing before it goes live.

5 revisions by 1 contributors · full history

Suggest edit