PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian inference for logistic models using Polya-Gamma latent variables

Nicholas G. Polson, James G. Scott, Jesse Windle

arXiv:1205.0310v3stat.MEstat.COstat.ML

TL;DR

Bayesian logistic regression is difficult because its likelihood is analytically inconvenient, and existing logit augmentation approaches can be approximate or involve multiple latent-variable layers. The paper introduces an exact Polya-Gamma augmentation for binomial log-odds models, and reports that it is most efficient in virtually all cases beyond simple logit models with abundant data and no hierarchical structure.

  • Problem

    Bayesian logistic regression has an analytically inconvenient likelihood, while prior logit augmentation approaches could be approximate or substantially more complex than the probit model’s latent-variable method.

    Source: §1 · pp. 1–2

  • Method

    The paper introduces an exact data-augmentation algorithm using Polya-Gamma variables to represent binomial log-odds likelihoods as Gaussian mixtures.

    Source: §1 · p. 2

  • Results

    The method is most efficient in virtually all cases except simple logit models with abundant data and no hierarchy, where it is a close second to carefully tuned independence Metropolis-Hastings.

    Source: §1 · p. 3

  • Takeaways & Limitations

    The approach provides simple Gibbs-sampling methods for Bayesian logistic models that exploit standard normal linear-model theory and are simpler than previous schemes.

    Source: §6 · p. 21

    With relatively high counts, summing one Polya-Gamma draw per count at every MCMC iteration becomes a burden; the alternative sampler is beyond this paper’s scope.

    Source: §4.4 · p. 17, §5 · p. 20

Abstract

from arXiv · show

We propose a new data-augmentation strategy for fully Bayesian inference in models with binomial likelihoods. The approach appeals to a new class of Polya-Gamma distributions, which are constructed in detail. A variety of examples are presented to show the versatility of the method, including logistic regression, negative binomial regression, nonlinear mixed-effects models, and spatial models for count data. In each case, our data-augmentation strategy leads to simple, effective methods for posterior inference that: (1) circumvent the need for analytic approximations, numerical integration, or Metropolis-Hastings; and (2) outperform other known data-augmentation strategies, both in ease of use and in computational efficiency. All methods, including an efficient sampler for the Polya-Gamma distribution, are implemented in the R package BayesLogit. In the technical supplement appended to the end of the paper, we provide further details regarding the generation of Polya-Gamma random variables; the empirical benchmarks reported in the main manuscript; and the extension of the basic data-augmentation framework to contingency tables and multinomial outcomes.

1 Introduction

Bayesian logistic regression has an analytically inconvenient likelihood, and earlier logit augmentation methods could be approximate or involve multiple latent-variable layers. The paper introduces an exact, simple Polya-Gamma augmentation and reports that it is most efficient in virtually all cases beyond simple, abundant-data models without hierarchy.

1 Introduction

  • The method represents binomial likelihoods parametrized by log odds as Gaussian mixtures, enabling Gibbs sampling with Gaussian parameter draws and a single layer of Polya-Gamma latent variables.It applies to Bayesian logistic regression and negative-binomial models for overdispersed counts.Source: §1 · p. 2
  • An exact accept/reject sampler for Polya-Gamma draws avoids truncating the infinite-sum representation and requires only exponential and inverse-Gaussian draws for PG(1, c).For PG(1, c), the acceptance probability is uniformly bounded below at 0.99919, with no tuning required.Source: §1 · pp. 2–3
  • The Polya-Gamma method is a close second to carefully tuned independence Metropolis-Hastings for simple logit models with abundant data and no hierarchy, and is most efficient in virtually all other cases.The authors note one exception: negative-binomial regression with many counts per observation and no hierarchical prior.Source: §1 · p. 3
  • For negative-binomial regression with many counts per observation and no hierarchical prior, summing independent PG(1, c) draws can reduce the method’s effective sampling rate.The authors report that effective sample size remains best in this exception; a faster approximation is available at the cost of a more complex latent-variable structure.Source: §1 · p. 3
  • The authors report speed and simplicity advantages across structured Bayesian models for binary and count data, with larger efficiency advantages as models become more complex and require expensive operations.They identify Gaussian-process spatial models as a setting with especially large differences.Source: §1 · p. 3

2 The P´olya-Gamma distribution

The paper constructs the Polya-Gamma family from infinite convolutions of gamma distributions. It develops the zero-tilt case, extends it by exponential tilting, and records properties useful for simulation and inference.

2.1 The case PG(b, 0)

  • For PG(b, 0), the paper defines a family through an infinite convolution of gamma distributions, with an equivalent infinite-sum representation using independent gamma variables.The PG(b, 0) family is related to the J∗(b) distributions surveyed by Biane et al. (2001), with PG(b, 0) distributed as J∗(b)/4.Source: §2.1 · p. 4

2.2 The general PG(b, c) class

  • The general PG(b, c) family arises by exponentially tilting the PG(b, 0) density and can also be written as an infinite convolution of gamma distributions.The normalizing constant follows from the PG(b, 0) Laplace transform.Source: §2.2 · p. 5

2.3 Further properties

  • The PG(b, c) density can be expressed as an alternating-sign sum of inverse-Gaussian densities, a characterization used in the Polya-Gamma simulation method.The exponential tilt is normalized using the PG(b, 0) Laplace transform.Source: §2.3 · pp. 5–6
  • All finite moments of Polya-Gamma variables are available in closed form, and the family is closed under convolution when independent variables share the same scale parameter.The closed-form expectation supports EM algorithms; convolution adds the shape parameters.Source: §2.3 · p. 6

3 The data-augmentation strategy

The Pólya–Gamma augmentation represents binomial log-odds likelihoods as Gaussian mixtures and yields Gibbs updates using Gaussian parameter draws and Pólya–Gamma latent draws. The section distinguishes this exact, single-layer construction from earlier schemes and illustrates its use in a mixed model.

3.1 Main result

  • The integral identity makes binomial log-odds likelihoods Gaussian mixtures, with Pólya–Gamma conditional latent variables and conditionally Gaussian updates for regression parameters.This gives a two-step Gibbs sampler using Gaussian draws for the main parameters and Pólya–Gamma draws for one latent layer.Source: §3.1 · pp. 7–8

3.2 Existing data-augmentation schemes

  • Unlike earlier random-utility approaches, the method uses an infinite gamma convolution, gives exact posterior draws, and requires only one layer of latent variables.The compared schemes use additional latent-variable layers, and one discrete normal-mixture approach is an approximation.Source: §3.2 · pp. 8–10
  • For negative-binomial regression, the Pólya–Gamma method outperforms the discrete-mixture approach in some cases but not others.The paper relates this variation to the inner workings of its sampler.Source: §3.2 · p. 10

3.3 Mixed model example

  • Shrinkage toward a global mean matters because the unbalanced design would otherwise give several district random intercepts unrealistic values.The figure displays posterior distributions for district intercepts.Source: §3.3 · p. 11, Figure 2 · p. 11
  • In the Bangladesh contraception example, a mixed model accommodates district effects despite sparse districts, and 2,000 draws take 27.3 seconds versus 26.1 seconds for logistic regression.The authors present the analysis as a proof of concept, not a definitive model of the data.Source: §3.3 · pp. 10–11

4 Simulating P´olya-Gamma random variables

The paper develops an exact accept–reject sampler for PG(1,z) using an alternating-series representation of a tilted Jacobi distribution. It analyzes acceptance and series-checking costs, then extends sampling to integer shape parameters by adding PG(1,z) draws.

4.1 The PG(1,z) sampler

  • The exact PG(1,z) sampler draws from a tilted Jacobi distribution using an alternating-series accept–reject procedure, then rescales the draw by one quarter.Its proposal can be sampled from a mixture of inverse-Gaussian and exponential distributions.Source: §4.1 · pp. 12–13

4.2 Analysis of acceptance rate

  • The optimized truncation point is independent of z, and the sampler rejects no more than 9 proposals per 10,000 on average across z.The reported lower bound on acceptance probability is 0.99919, attained near z ≃ 1.378.Source: §4.2 · p. 14

4.3 Analysis of tail probabilities

  • At the worst-case envelope z ≃ 1.378, the algorithm averages 1.0016 proposals, while the probability of needing more partial sums decays rapidly.Together, low rejection and rapidly declining partial-sum tail probabilities support the sampler’s efficiency.Source: §4.3 · pp. 15–16

4.4 The general PG(b, z) case

  • For integer b, PG(b,z) draws are formed by summing b independent PG(1,z) draws, making runtime roughly linear in total counts.The authors use this approach in their benchmarks; summing many draws can become a bottleneck for large shape parameters.Source: §4.4 · pp. 16–17

5 Experiments

Across binary-logit datasets, the Pólya–Gamma method outperformed other data-augmentation schemes, while comparisons with Metropolis samplers depended on model complexity and proposal choice. Results also favor the method in mixed and spatial models, with negative-binomial performance varying with count size and sampling burden.

Benchmarking approach

  • The benchmarks compare effective sample size and effective sampling rate, defined as effective sample size per second of runtime, using R timings and median effective sample size across model parameters.The reported timings include compiled C calls for Pólya–Gamma draws, like R routines for common distributions.Source: §5 · p. 17

Binary logit models

  • Across 6 real and 2 simulated binary-logit datasets, Pólya–Gamma was always more efficient than the next-best data-augmentation scheme, typically by a factor of 200%–500%.The comparison covers both effective sample size and effective sampling rate.Source: Table 1 · p. 18, §5 · p. 18
  • Pólya–Gamma always had higher effective sample size than the two default Metropolis samplers, but the best Metropolis algorithm had a higher effective sampling rate on 5 of 8 datasets.On Australian credit data, the Gaussian proposal yielded 122 samples per second versus 2.6 for a similar multivariate t6 proposal; performance depended on proposal choice.Source: §5 · p. 18

Binary logistic mixed models

  • For mixed models, Pólya–Gamma was slightly better than a Laplace-based independence Metropolis sampler on balanced synthetic data and much better on two highly unbalanced real datasets.The authors note that better-tuned Metropolis samplers are possible; Pólya–Gamma requires no tuning and had optimal or near-optimal performance across a range of cases.Source: §5 · p. 19, Table 2 · p. 19

Negative-binomial models

  • Pólya–Gamma consistently gave the best effective sample sizes in negative-binomial regression, but its effective sampling rate suffered with large counts or non-integral over-dispersion.For integer b, the sampler sums b PG(1, ψ) draws; in negative-binomial models this requires at least Σ_i y_i draws per iteration, which becomes burdensome when counts are high.Source: §5 · pp. 19–20
  • With 400 observations and 3 regressors, Pólya–Gamma had a better effective sampling rate than FS09 for Sim1 (average count 8), while FS09 was faster for Sim2 (average count 24).Both comparisons restrict the sampler to integer over-dispersion parameters.Source: Table 3 · p. 19, §5 · p. 20

Spatial models

  • In two spatial Gaussian-process examples, Pólya–Gamma was more efficient by factors of 3 and 5, respectively, while both methods required costly Cholesky decompositions at each iteration.Its higher effective sample size wastes fewer of those expensive draws for the main parameter.Source: Table 3 · p. 19, §5 · pp. 20–21

6 Discussion

The discussion presents Pólya–Gamma augmentation as a simple Gibbs-sampling route for Bayesian logistic models, supported by an efficient sampler and a general posterior-tilting characterization. Its benchmarks are serial, and large integral shape parameters remain a potential target for GPU acceleration.

Discussion

  • The Pólya–Gamma augmentation yields simple Gibbs samplers for logistic models that exploit standard normal linear-model theory, alongside an accept/reject sampler with strong efficiency guarantees.The paper identifies the scheme as simpler than previous approaches.Source: §6 · p. 21
  • The evidence suggests Pólya–Gamma is the best current method for complex Bayesian hierarchical models with binomial likelihoods and may support exact Bayesian treatment of mixture-of-logits classifiers.Applications to those machine-learning problems are described as an active area of research.Source: §6 · p. 21
  • Exponential tilting links the PG(b, 0) prior to the conditional PG(b, c) distribution, which the authors characterize using its moment-generating function rather than Bayes’ rule for densities.They relate this reasoning to the Lévy-penalty framework and Brownian-motion sojourn times.Source: §6 · p. 21
  • The benchmarks use serial computation; GPU sampling is less straightforward but could address sluggish draws for PG(n, z) with large integral n.That sampling burden was the largest source of inefficiency in the negative-binomial results.Source: §6 · p. 21

S1 Details of P´olya-Gamma sampling algorithm

The supplement describes an exact accept/reject sampler for PG(1, z), using alternating-series decisions and truncated inverse-Gaussian proposals. For other shape parameters, it gives a sum construction for integers and a finite-sum approximation for fractional parts.

PG(1, z) sampler

  • The PG(1, z) sampler uses an alternating-series accept/reject procedure, with truncated inverse-Gaussian sampling strategies chosen according to whether 1/z exceeds the truncation point t.For 1/z > t, it uses an inverse-chi-square-based proposal and gives the acceptance probability as exp(−xz²/2).Source: S1 Details of P´olya-Gamma sampling… · pp. 25–26

Other shape parameters

S2 Benchmarks: overview

The benchmarks compare the Polya-Gamma method with alternatives for binary logistic and negative-binomial regression, using effective sampling rate as the primary measure. The authors report favorable performance against several data-augmentation methods, while cautioning that the empirical comparisons depend on implementation and hardware.

S2 Benchmarks: overview

  • The Polya-Gamma method performs better than the methods of O’Brien and Dunson (2004), Gramacy and Polson (2012), and Frühwirth-Schnatter and Frühwirth (2010).The comparisons cover binary logistic regression and negative-binomial regression for count data.Source: S2 Benchmarks: overview · p. 27
  • The empirical comparisons warrant caution because ESR can vary with implementation, programming language, and hardware.Experiments used R 2.15.1 on Ubuntu with an 8GB machine and an Intel Core i5 quad-core processor; the authors’ C code did not use parallelism.Source: S2 Benchmarks: overview · pp. 27–28
  • Effective sampling rate (ESR), the effective sample size per second, is the primary comparison metric.The reported rates exclude initialization, preprocessing, and burn-in time.Source: S2 Benchmarks: overview · p. 28
  • Each dataset was evaluated in 10 MCMC runs of 12,000 samples, with the first 2,000 discarded as burn-in.The retained samples formed 10 batches of 10,000 samples.Source: S2 Benchmarks: overview · p. 28

S3 Benchmarks: binary logistic regression

The supplement benchmarks binary logistic regression across real and synthetic datasets and describes the methods compared. Results are reported in Tables 4–11, averaged over 10 runs.

S3.1 Data Sets

  • The datasets include prostate nodal status, diabetes, heart disease, credit-risk data, and synthetic designs with orthogonal or collinear predictors.The dataset descriptions also note data cleaning and predictor encoding choices, including removing zero-coded Pima measurements and using two German Credit predictor encodings.Source: S3.1 Data Sets · pp. 28–29

S3.2 Methods

  • The comparisons include Polya-Gamma Gibbs sampling, independence Metropolis samplers, and several latent-variable augmentation methods.The listed alternatives include approximate discrete-mixture methods and the methods of O’Brien and Dunson (2004) and Gramacy and Polson (2012).Source: S3.2 Methods · pp. 29–30

S3.3 Results

S4 Benchmarks: logit mixed models

The supplement evaluates Polya-Gamma sampling for binary logistic mixed models with group-dependent random intercepts. The authors describe it as a sensible default for these models and report that independence Metropolis can perform poorly in the considered mixed-model settings.

S4 Benchmarks: logit mixed models

  • The Polya-Gamma method is presented as a sensible default choice for nonlinear mixed-effect models.The benchmark set includes synthetic, presidential-polling, and Xerop data, with random intercepts and varying fixed-effect dimensions.Source: S4 Benchmarks: logit mixed models · p. 33
  • Independence Metropolis, which usually works well for binary logistic regression, does not work well for the mixed models considered.For the polls dataset, at least two heuristics indicate that the Laplace approximation will be a poor proposal.Source: S4 Benchmarks: logit mixed models · p. 33
  • Table 12 reports effective sample sizes for (α, β, m), excluding φ; two competing methods are thinned every 10 samples and still have lower ESS than Polya-Gamma.The table covers three benchmarks with random effects limited to group-dependent intercepts.Source: Table 12 · p. 34
  • The posterior mode and mean differ most in φ, where the posterior is especially flat and mass concentrates at relatively small values.The reported minimum eigenvalue of the Hessian at the mode points predominantly in the φ direction; the authors say the behavior is not due to initialization or stopping rules.Source: S4 Benchmarks: logit mixed models · p. 34

S5 Benchmarks: negative-binomial models

The negative-binomial benchmarks show that Polya-Gamma sampling can have higher effective sample size but lower effective sampling rate when counts are larger. In spatial models, where sampling the main parameter is expensive, the method is competitive or better than the comparison approach.

S5 Benchmarks: negative-binomial models

  • The synthetic negative-binomial models use N = 400 data points and estimate d with a random-walk Metropolis-Hastings step over integers.The parameter is not efficiently handled at non-integer values by the Polya-Gamma method or the Fussl (2012) R package.Source: S5 Benchmarks: negative-binomial models · pp. 34, 37
  • Figure 4 shows that the posterior mode and mean differ most in φ, while marginal posterior mass is concentrated at relatively small φ values.The log posterior is very flat when moving along φ.Source: Figure 4 · p. 35
  • Polya-Gamma sampling has superior effective sample size in both simulated negative-binomial cases, but a lower effective sampling rate in the second case.The rate disadvantage is attributed to drawing PG(n, z) by summing n PG(1, z) variables; an efficient sampler for arbitrary n remains an open challenge.Source: Table 13 · p. 36, S5 Benchmarks: negative-binomial models · p. 37
  • In spatial models, Polya-Gamma sampling is competitive or sometimes better when expensive Gaussian-process draws reduce the relative time spent sampling auxiliary variables.The examples use 256 and 1,000 spatial points; the larger problem shows a greater performance improvement over Frühwirth-Schnatter et al. (2009).Source: S5 Benchmarks: negative-binomial models · p. 37

S6 Extensions

The supplement extends Polya-Gamma augmentation to multi-center tables and multinomial logistic regression. In glass-identification data, it produces reasonable posterior distributions where an independent Metropolis-Hastings proposal is hampered by ill-conditioning and perfect separation.

S6.1 2 × 2 × N tables

  • For the eight-center topical-cream study, the Gibbs sampler estimates center-specific log-odds ratios and displays posterior means and 95% credible intervals alongside maximum-likelihood estimates.The maximum-likelihood estimate for the control group is −∞ in centers 5 and 6 because neither had observed successes.Source: S6.1 2 × 2 × N tables · p. 38, Figure 5 · p. 39

S6.2 Higher-order tables

  • For higher-order tables, Polya-Gamma augmentation yields normal-kernel conditional updates and a Gibbs sampler that cycles through centers and response-specific parameter vectors.The extension considers multi-center studies with multiple treatment arms and possible outcomes.Source: S6.2 Higher-order tables · pp. 39, 41

S6.3 Multinomial logistic regression

  • For multinomial logistic regression, the method introduces Polya-Gamma variables into category-wise conditional likelihoods and iterates Gibbs updates across categories.The paper notes that this avoids analytic posterior approximations, an advantage when observations are not much more numerous than parameters.Source: S6.3 Multinomial logistic regression · pp. 41–42
  • In glass-identification data, the Polya-Gamma method produces reasonable posterior distributions despite perfect separation in category 6, which impairs the independent Metropolis-Hastings approach.The joint proposal is also challenged by a poorly conditioned Hessian under a vague prior.Source: S6.3 Multinomial logistic regression · p. 42

Related papers

Loading 1205.0310v3…