tulpa_ep: Expectation-Propagation fit for a GLM

In tulpa: Templated Unified Library for Posterior Approximation in Bayesian Hierarchical Models

Expectation-Propagation fit for a GLM

Description

Fits a generalized linear model with a Gaussian coefficient prior by Expectation Propagation: the posterior is approximated by a Gaussian whose per-observation sites match the moments of the tilted distribution (via Gauss-Hermite quadrature). EP is exact for a Gaussian likelihood and typically more accurate than Laplace on skewed likelihoods, since it matches marginal moments rather than the mode curvature.

Usage

tulpa_ep(
  formula,
  data,
  family = "binomial",
  phi = 1,
  phi2 = NULL,
  n_trials = NULL,
  beta_prior = .tulpa_default_beta_prior("ep"),
  control = list()
)

Arguments

formula	Model formula.
data	A data frame.
family	Character family name (see family_names()).
phi	Dispersion passed to the family, held fixed. One convention at every door: for gaussian / lognormal this is the residual VARIANCE (the SD is sqrt(phi)), for neg_binomial_2 the size, gamma the shape, beta the precision, t the scale; binomial and poisson ignore it. The compiled kernels parameterize the two variance families by the residual SD and are handed sqrt(phi) at the boundary. Defaulted, it conditions at 1 and says so with a warning, since for a family that reads a dispersion that is a modelling choice rather than a neutral value. fit$phi_estimated records whether the value on the fit was estimated or conditioned on.
phi2	Optional second dispersion (Student-t degrees of freedom for family = "t"; default 4 when NULL).
n_trials	Binomial denominators (length nrow(data)), or NULL (= 1).
beta_prior	Fixed-effect prior as list(mean, sd): a mean-zero (mean = 0) Gaussian on every coefficient with SD sd (default the engine default, prior_normal(0, 2.5)). EP's site parameterisation assumes a mean-zero coefficient prior, so a non-zero mean errors – use a sampler (mode = "mala") for a shifted prior.
control	List: max_sweeps (default 50), tol (default 1e-6), damping (default 0.8), n_quad (Gauss-Hermite nodes, default 20), n_draws (default 2000), seed. An offset(...) term in formula is honoured: it shifts each observation's linear predictor before the likelihood is evaluated, and is not itself part of the fitted posterior.

Value

A tulpa_fit with means (posterior mean) and cov (posterior covariance) of the EP Gaussian, which is the posterior coef(), vcov(), summary.tulpa_fit() and confint() report; draws sampled from that Gaussian; log_marginal (the EP approximation); converged (the sweep tolerance was met AND every site's adaptive quadrature mode search converged AND no site's tilted variance floored – a site that is not well resolved by the quadrature turns this FALSE even when the outer sweep loop itself settled); n_site_not_converged, n_site_floored.

References

Minka (2001). Expectation Propagation for approximate Bayesian inference. UAI. Rasmussen & Williams (2006). Gaussian Processes for Machine Learning, Algorithm 3.5.

See Also

tulpa() (Laplace / sampler tiers), pathfinder() (VI).

Examples

set.seed(1)
d <- data.frame(x = rnorm(200))
d$y <- rbinom(200, 1, plogis(-0.3 + 0.8 * d$x))
fit <- tulpa_ep(y ~ x, data = d, family = "binomial")
coef(fit)

tulpa documentation built on Oct. 19, 2026, 5:09 p.m.