Bayesian Logistic Regression

Bayesian Logistic Regression

The Bayesian logistic regression model is
$y_i\mid \mathbf{\beta },{\mathbf{x}}_i\sim \text{Bernoulli}(\sigma ({\mathbf{x}}_i^{\top }\mathbf{\beta })),\mathbf{\beta }\sim N({\mathbf{\beta }}_0,{\mathbf{\Sigma }}_0),$
with $\sigma (u)=(1+e^{-u}{)}^{-1}$ the logistic sigmoid. The likelihood is
$p(\mathbf{y}\mid \mathbf{\beta })={\prod }_{i=1}^n\sigma ({\mathbf{x}}_i^{\top }\mathbf{\beta }{)}^{y_i}(1-\sigma ({\mathbf{x}}_i^{\top }\mathbf{\beta }){)}^{1-y_i}.$
The Bernoulli likelihood is not conjugate to any standard prior on $\mathbf{\beta }$, so the posterior $p(\mathbf{\beta }\mid \mathbf{y})$ has no closed form. Three standard approaches:

1. Laplace approximation — approximate the posterior by a Gaussian around the MAP:
   $p(\mathbf{\beta }\mid \mathbf{y})\approx N\left(\hat{\mathbf{\beta }},{\mathbf{\Sigma }}_{\text{Lap}}\right),\hat{\mathbf{\beta }}=\arg {\max }_{\mathbf{\beta }}\log p(\mathbf{\beta }\mid \mathbf{y}),{\mathbf{\Sigma }}_{\text{Lap}}=[-{\nabla }^2\log p(\hat{\mathbf{\beta }}\mid \mathbf{y}){]}^{-1}.$
   The Hessian is closed form:
   $-{\nabla }^2\log p(\mathbf{\beta }\mid \mathbf{y})={\mathbf{X}}^{\top }\mathbf{W}(\mathbf{\beta })\mathbf{X}+{\mathbf{\Sigma }}_0^{-1},$
   with $\mathbf{W}(\mathbf{\beta })=\mathrm{diag}(p_i(1-p_i))$ and $p_i=\sigma ({\mathbf{x}}_i^{\top }\mathbf{\beta })$.
2. Metropolis-Hastings with Laplace proposal — use the Laplace Gaussian as a proposal distribution for MH. See 3.11.1 Metropolis-Hastings.
3. Gibbs via Pólya-Gamma data augmentation — introduce latent ${\omega }_i\sim \text{PG}(1,0)$ variables so that the conditional posterior for $\mathbf{\beta }$ given $\mathbf{\omega }$ is Gaussian (Polson-Scott-Windle 2013). See 3.11.2 Gibbs Sampling.

MAP via Newton’s method (IRLS). The MAP optimization is convex in $\mathbf{\beta }$ (log-concave posterior) and solved by Newton’s method:
${\mathbf{\beta }}_{k+1}={\mathbf{\beta }}_k+[{\mathbf{X}}^{\top }{\mathbf{W}}_k\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{]}^{-1}[{\mathbf{X}}^{\top }(\mathbf{y}-{\mathbf{p}}_k)-{\mathbf{\Sigma }}_0^{-1}({\mathbf{\beta }}_k-{\mathbf{\beta }}_0)].$
Each Newton step is equivalent to a weighted least-squares fit against the working response $z_{k,i}={\mathbf{x}}_i^{\top }{\mathbf{\beta }}_k+(y_i-p_{k,i})/(p_{k,i}(1-p_{k,i}))$ with weights ${\mathbf{W}}_k$. This is iteratively reweighted least squares (IRLS), the standard frequentist logistic-regression solver, recovered here as a Bayesian MAP algorithm.

Posterior predictive for a new input ${\mathbf{x}}_{\star }$:
$P({\tilde{y}}_{\star }=1\mid \mathbf{y})=\int \sigma ({\mathbf{x}}_{\star }^{\top }\mathbf{\beta })p(\mathbf{\beta }\mid \mathbf{y})d\mathbf{\beta }.$
Under the Laplace approximation, a one-dimensional integral over $u={\mathbf{x}}_{\star }^{\top }\mathbf{\beta }\sim N({\hat{\mu }}_{\star },{\hat{\sigma }}_{\star }^2)$ remains, well-approximated by $\sigma ({\hat{\mu }}_{\star }/\sqrt{1+\pi {\hat{\sigma }}_{\star }^2/8})$ (MacKay’s probit approximation).

Prerequisites

• 3.1 Bayesian Framework — Bayesian mechanics.
• 3.8 Bayesian Regression — parent.
• 3.9 Laplace Approximation — the standard non-conjugate tool applied here.
• 2.4. Logistic Regression — frequentist counterpart for comparison.

Intuition

The sigmoid squashes a real-valued linear combination ${\mathbf{x}}^{\top }\mathbf{\beta }$ into a probability, and the Bernoulli likelihood pins each observed $y_i$ to that probability. Unlike the Gaussian likelihood, the sigmoid is not conjugate to Gaussian priors on $\mathbf{\beta }$ — the log-likelihood is ${\sum }_i[y_i{\mathbf{x}}_i^{\top }\mathbf{\beta }-\log (1+e^{{\mathbf{x}}_i^{\top }\mathbf{\beta }})]$, which is not quadratic in $\mathbf{\beta }$, so prior-likelihood multiplication does not yield another Gaussian. But the log-posterior is concave, so it has a unique maximum, and near that maximum it is well-approximated by a quadratic — which is exactly the setup that makes the Laplace approximation accurate.

IRLS emerges because each Newton step re-expresses the log-likelihood around the current estimate as a weighted least-squares problem. The weights $p_i(1-p_i)$ are smallest when $p_i$ is near 0 or 1 (confident predictions) and largest at $p_i=1/2$ (maximally uncertain), so ambiguous observations carry more Fisher information and pull the estimate harder. The Bayesian addition is the prior-precision term ${\mathbf{\Sigma }}_0^{-1}$: regularized logistic regression (${\ell }_2$-penalty on $\mathbf{\beta }$) is precisely the Bayesian MAP with a Gaussian prior, and the penalty strength is the prior precision.

The Laplace approximation is typically very accurate in this setting because the log-posterior is globally concave and the Hessian is well-defined and invertible wherever the design is not singular. The approximation fails when the posterior is multimodal (e.g., under separable data with strongly conflicting priors) or strongly skewed (rare in practice with $n>d$ and a proper Gaussian prior).

Properties

(Ridge-Regularized Logistic Regression as MAP) The ${\ell }_2$-penalized logistic regression is the Bayesian MAP under a Gaussian prior.

The objective $-{\sum }_i\log p(y_i\mid \mathbf{\beta })+\frac{\lambda }{2}\Vert \mathbf{\beta }{\Vert }^2$ equals $-\log p(\mathbf{y},\mathbf{\beta })$ up to constants when the prior is $\mathbf{\beta }\sim N(0,{\lambda }^{-1}\mathbf{I})$. The Bayesian interpretation turns the penalty strength into a prior precision and gives a principled way to choose it (Type-II MLE / evidence maximization via Laplace; see 3.7 Marginal Likelihood).

(IRLS as Newton’s Method on the MAP) Each Newton step is a weighted least-squares fit.

Starting from ${\mathbf{\beta }}_k$, define working response $z_{k,i}={\mathbf{x}}_i^{\top }{\mathbf{\beta }}_k+(y_i-p_{k,i})/(p_{k,i}(1-p_{k,i}))$ and weights $w_{k,i}=p_{k,i}(1-p_{k,i})$. Then
${\mathbf{\beta }}_{k+1}=({\mathbf{X}}^{\top }{\mathbf{W}}_k\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}({\mathbf{X}}^{\top }{\mathbf{W}}_k{\mathbf{z}}_k+{\mathbf{\Sigma }}_0^{-1}{\mathbf{\beta }}_0),$
which is the Bayesian weighted-ridge normal equation. Convergence is quadratic near the MAP; the algorithm is numerically stable as long as the Hessian is well-conditioned, which the prior term ${\mathbf{\Sigma }}_0^{-1}$ guarantees even when $\mathbf{X}$ is rank-deficient.

(Laplace Approximation is the Default Posterior) The Gaussian approximation around the MAP suffices for most inference tasks.

The posterior covariance is ${\mathbf{\Sigma }}_{\text{Lap}}=({\mathbf{X}}^{\top }\mathbf{W}(\hat{\mathbf{\beta }})\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}$. Credible intervals for individual coefficients come from the diagonal of ${\mathbf{\Sigma }}_{\text{Lap}}$, exactly mirroring the frequentist Wald interval with $({\mathbf{X}}^{\top }\mathbf{WX}{)}^{-1}$ as the inverse Fisher information. The Bayesian interpretation adds a coherent story for small samples, integrates naturally with hierarchical models, and provides proper posterior predictive distributions; the approximation’s accuracy is essentially identical to the Wald approximation in the frequentist literature.

(MacKay’s Probit Approximation for Posterior Predictive) The 1-D integral of a sigmoid against a Gaussian has a closed-form approximation.

$\int \sigma (u)N(u;\mu ,s^2)du\approx \sigma \left(\frac{\mu }{\sqrt{1+\pi s^2/8}}\right),$
accurate to about 0.2% for all $\mu ,s$. This reduces posterior prediction to a scalar sigmoid evaluation at a variance-corrected logit, avoiding any numerical integration. The “pull toward 0.5” as predictive variance grows is automatic: a Bayesian predictor is less confident than the MAP predictor precisely because it integrates over posterior uncertainty.

(Separable Data Gives Improper MLE, Proper Posterior) Classical logistic MLE is undefined when classes are linearly separable; the Bayesian posterior remains well-defined with any proper prior.

If there exists ${\mathbf{\beta }}^{\ast }$ with ${\mathbf{x}}_i^{\top }{\mathbf{\beta }}^{\ast }>0$ iff $y_i=1$, the log-likelihood is unbounded above (sending $\mathbf{\beta }\to c{\mathbf{\beta }}^{\ast }$ with $c\to \infty$ gives perfect classification with zero loss), so the MLE does not exist. The proper Gaussian prior bounds $\mathbf{\beta }$‘s magnitude and yields a unique MAP. This is one of the cleanest practical advantages of Bayesian logistic regression: regularization is built in and separable data is not a pathology.

(Pólya-Gamma Data Augmentation for Exact Gibbs) Introducing latent ${\omega }_i\sim \text{PG}(1,{\mathbf{x}}_i^{\top }\mathbf{\beta })$ variables restores conditional Gaussianity.

Polson-Scott-Windle (2013) showed that the Bernoulli likelihood factor $\sigma (u{)}^y(1-\sigma (u){)}^{1-y}=e^{\kappa u}/(1+e^u)$ (with $\kappa =y-1/2$) can be represented as $\int e^{\kappa u-\omega u^2/2}p(\omega )d\omega$ where $\omega \sim \text{PG}(1,0)$. Conditioning on $\omega$, the likelihood for $\mathbf{\beta }$ becomes Gaussian, enabling Gibbs updates with closed-form conditionals. This is the go-to MCMC method for Bayesian logistic regression when Laplace is not accurate enough.

(Posterior Variance Inflates Predictions Away From the Decision Boundary) Predictions are most uncertain near ${\mathbf{x}}^{\top }\mathbf{\beta }\approx 0$.

The posterior predictive probability curve, as a function of ${\mathbf{x}}_{\star }$, is flatter than the plug-in $\sigma ({\mathbf{x}}_{\star }^{\top }\hat{\mathbf{\beta }})$ curve — especially in regions where ${\hat{\sigma }}_{\star }^2$ is large (test inputs far from training). A Bayesian predictor interpolates more cautiously into covariate regions where the training data is sparse, giving predictions that are closer to 0.5 than the MLE plug-in.

Examples

Example 1: Single-Covariate Logistic Regression via Laplace

Data $\{(x_i,y_i){\}}_{i=1}^{100}$ with $x_i\sim N(0,1)$ and $y_i\sim \text{Bernoulli}(\sigma (1+0.5x_i))$. Prior $\mathbf{\beta }=({\beta }_0,{\beta }_1)\sim N(0,10\mathbf{I})$.

IRLS: starting from ${\mathbf{\beta }}_0=0$, iterate Newton steps until convergence. Typically reaches $\hat{\mathbf{\beta }}\approx (1.02,0.48)$ in 4–5 iterations.
Hessian at MAP: ${\mathbf{X}}^{\top }\text{diag}(p_i(1-p_i))\mathbf{X}+0.1\mathbf{I}$, invertible by regularization.
Laplace posterior: $N(\hat{\mathbf{\beta }},{\mathbf{\Sigma }}_{\text{Lap}})$; marginal standard deviations $\approx (0.20,0.21)$.
Credible interval for ${\beta }_1$: $0.48\pm 1.96\cdot 0.21=[0.07,0.89]$ — excludes zero, consistent with the true coefficient 0.5.
Predictive at $x_{\star }=1$: ${\hat{\mu }}_{\star }=1.02+0.48=1.50$, ${\hat{\sigma }}_{\star }^2={\mathbf{x}}_{\star }^{\top }{\mathbf{\Sigma }}_{\text{Lap}}{\mathbf{x}}_{\star }\approx 0.075$. MacKay approximation gives predictive probability $\sigma (1.50/\sqrt{1+\pi \cdot 0.075/8})=\sigma (1.48)=0.815$ — slightly pulled toward 0.5 compared with the plug-in $\sigma (1.50)=0.818$.

Example 2: Separable Data Saved by the Prior

Data $\{(x_i,y_i){\}}_{i=1}^{20}$ with $x_i$ a single covariate such that $y_i=1[x_i>0]$ perfectly. The MLE does not exist: the likelihood grows as ${\beta }_1\to \infty$ with ${\beta }_0=0$.

With prior $\mathbf{\beta }\sim N(0,4\mathbf{I})$: the log-posterior is ${\sum }_i\log \sigma (y_i^{\prime }({\beta }_0+{\beta }_1{x}_i))-\Vert \mathbf{\beta }{\Vert }^2/8$ with $y_i^{\prime }=2y_i-1$. The prior penalty prevents ${\beta }_1$ from running off to infinity, and the MAP exists (typically ${\beta }_1\approx 3$–$5$). The posterior is proper and concentrates around the MAP.
Predictions: at an input far from the training data (say $x_{\star }=5$), the Laplace posterior gives predictive probability well short of 1 (perhaps 0.95 instead of the MLE’s 1 - $ϵ$), reflecting genuine uncertainty about extrapolation.

Classical logistic regression software either fails or silently returns an enormous MLE; the Bayesian model handles it without comment.

Example 3: Evidence-Based Tuning of the Prior Scale

The prior precision $\tau$ in $\mathbf{\beta }\sim N(0,{\tau }^{-1}\mathbf{I})$ controls regularization strength. Evidence-based tuning picks $\tau$ to maximize the Laplace-approximated marginal likelihood:
$\hat{\tau }=\arg {\max }_{\tau }\log p(\mathbf{y}\mid \hat{\mathbf{\beta }}(\tau ))-\frac{\tau }{2}\Vert \hat{\mathbf{\beta }}(\tau ){\Vert }^2-\frac{1}{2}\log |{\mathbf{X}}^{\top }\mathbf{WX}+\tau \mathbf{I}|+\frac{d}{2}\log \tau +\text{const}.$

This is the logistic-regression analog of the evidence-framework procedure MacKay used for neural networks in the 1990s. For a typical dataset with $n=1000,d=20$, the evidence surface is unimodal in $\log \tau$ with maximum near $\tau \approx 0.1$–$1$; cross-validation agrees to within a factor of 2. Unlike CV, evidence tuning requires no held-out splits and integrates naturally with model comparison across different feature sets.

Derivation

Goal: derive the Laplace posterior covariance from the Hessian of the log-posterior.

The log-posterior is
$\log p(\mathbf{\beta }\mid \mathbf{y})={\sum }_{i=1}^n[y_i{\mathbf{x}}_i^{\top }\mathbf{\beta }-\log (1+e^{{\mathbf{x}}_i^{\top }\mathbf{\beta }})]-\frac{1}{2}(\mathbf{\beta }-{\mathbf{\beta }}_0{)}^{\top }{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0)+\text{const}.$
Differentiating twice with respect to $\mathbf{\beta }$: the first derivative is
$\nabla \log p(\mathbf{\beta }\mid \mathbf{y})={\mathbf{X}}^{\top }(\mathbf{y}-\mathbf{p}(\mathbf{\beta }))-{\mathbf{\Sigma }}_0^{-1}(\mathbf{\beta }-{\mathbf{\beta }}_0),$
with $\mathbf{p}(\mathbf{\beta })=(\sigma ({\mathbf{x}}_1^{\top }\mathbf{\beta }),\dots ,\sigma ({\mathbf{x}}_n^{\top }\mathbf{\beta }){)}^{\top }$. The second derivative is
${\nabla }^2\log p(\mathbf{\beta }\mid \mathbf{y})=-{\mathbf{X}}^{\top }\mathbf{W}(\mathbf{\beta })\mathbf{X}-{\mathbf{\Sigma }}_0^{-1},$
with $\mathbf{W}(\mathbf{\beta })=\mathrm{diag}(p_i(1-p_i))$ using ${\sigma }^{\prime }(u)=\sigma (u)(1-\sigma (u))$. The Hessian is negative-definite (log-posterior is concave), so the MAP $\hat{\mathbf{\beta }}$ is unique. The Laplace approximation is
$p(\mathbf{\beta }\mid \mathbf{y})\approx N(\hat{\mathbf{\beta }},{\mathbf{\Sigma }}_{\text{Lap}}),{\mathbf{\Sigma }}_{\text{Lap}}=({\mathbf{X}}^{\top }\mathbf{W}(\hat{\mathbf{\beta }})\mathbf{X}+{\mathbf{\Sigma }}_0^{-1}{)}^{-1}.$
The Newton iteration ${\mathbf{\beta }}_{k+1}={\mathbf{\beta }}_k-[{\nabla }^2\log p{]}^{-1}\nabla \log p$ applied to this log-posterior is IRLS with prior-regularized weighted normal equations.

MacKay’s predictive approximation: for $u\sim N(\mu ,s^2)$, approximate $\int \sigma (u)p(u)du\approx \sigma (\mu /\sqrt{1+\pi s^2/8})$. This arises from approximating $\sigma (u)\approx \Phi (\lambda u)$ for $\lambda =\sqrt{\pi /8}$, where $\Phi$ is the standard-normal CDF, and then using the exact identity $\int \Phi (\lambda u)N(u;\mu ,s^2)du=\Phi (\lambda \mu /\sqrt{1+{\lambda }^2{s}^2})$.

See Also

• 3.8 Bayesian Regression — parent index.
• 3.9 Laplace Approximation — underlying approximation technique.
• 3.8.4 Bayesian Poisson Regression — non-conjugate count-response counterpart.
• 3.11.1 Metropolis-Hastings — alternative inference via sampling, often with Laplace proposal.
• 3.11.2 Gibbs Sampling — Pólya-Gamma augmentation enables exact Gibbs.
• 3.7 Marginal Likelihood — evidence-based tuning of the prior precision.
• 2.4. Logistic Regression — frequentist counterpart.