Tweedie’s Formula: A Demo (spike-and-slab priors)

1. Tweedie’s formula

Setup

Suppose we have a latent “true” parameter \(\mu\) drawn from some prior \(g\), and we observe \[ X = \mu + Z, \qquad Z \sim N(0,1), \qquad Z \perp \mu. \]

Given \(X = x\), we want the posterior mean \(E[\mu \mid X = x]\). Bayes’ rule gives the posterior density \[ p(\mu \mid x) \;=\; \frac{g(\mu)\,\varphi(x-\mu)}{f(x)}, \] where \(\varphi\) is the standard normal density and \(f(x) = \int g(\mu)\,\varphi(x-\mu)\,d\mu\) is the marginal density of \(X\).

Derivation

The key observation is \[ \frac{\partial}{\partial x}\,\varphi(x-\mu) \;=\; -(x-\mu)\,\varphi(x-\mu) \;=\; (\mu - x)\,\varphi(x-\mu). \]

Multiply by \(g(\mu)\) and integrate in \(\mu\): \[ \int (\mu - x)\, g(\mu)\,\varphi(x-\mu)\,d\mu \;=\; \frac{d}{dx}\int g(\mu)\,\varphi(x-\mu)\,d\mu \;=\; f'(x). \]

The left-hand side is \(f(x)\cdot\bigl(E[\mu \mid x] - x\bigr)\). Dividing by \(f(x)\), \[ \boxed{\,E[\mu \mid X = x] \;=\; x \;+\; \frac{f'(x)}{f(x)} \;=\; x \;+\; \frac{d}{dx}\log f(x).\,} \]

This is Tweedie’s formula (Robbins 1956; Efron 2011). The remarkable feature: the posterior mean depends on the prior \(g\) only through the marginal \(f\). If you can estimate \(\log f\) from data and differentiate it, you can do Bayes without ever writing down \(g\).

Spike-and-slab priors

Today we’ll use a spike-and-slab prior: \[ g \;=\; \pi_0\,\delta_0 \;+\; (1 - \pi_0)\,h, \qquad \pi_0 = 0.9, \] where \(h\) is the “slab” — Exp(1), Logistic, or \(t_2\) below. With 90% mass on \(\mu = 0\), most observations are pure noise; the interesting Bayesian question is which observations look real.

The marginal \(f\) inherits the same mixture structure (since convolution is linear): \[ f(x) \;=\; \pi_0\,\varphi(x) \;+\; (1-\pi_0)\,(h * \varphi)(x). \]

Tweedie’s formula applies unchanged — we just plug in this \(f\).

2. Numerical implementation

tweedie_spike_slab <- function(slab_pdf, x_grid, pi0 = 0.9,
                               mu_range = c(-40, 40), n_mu = 4001) {
  mu_grid <- seq(mu_range[1], mu_range[2], length.out = n_mu)
  dmu     <- mu_grid[2] - mu_grid[1]
  h_vals  <- slab_pdf(mu_grid)

  # slab contribution: (h * φ)(x) = ∫ h(μ) φ(x - μ) dμ
  slab_marg <- vapply(x_grid,
                      function(x) sum(h_vals * dnorm(x - mu_grid)) * dmu,
                      numeric(1))

  # spike contribution is exact: π₀ φ(x)
  f_x <- pi0 * dnorm(x_grid) + (1 - pi0) * slab_marg

  # centered finite difference of log f
  log_f  <- log(f_x)
  n      <- length(x_grid)
  dlogf  <- rep(NA_real_, n)
  dlogf[2:(n-1)] <- (log_f[3:n] - log_f[1:(n-2)]) /
                   (x_grid[3:n] - x_grid[1:(n-2)])

  data.frame(x = x_grid, f = f_x, Emu = x_grid + dlogf)
}

Sampler for the spike-and-slab: a Bernoulli(\(1-\pi_0\)) “is it real?” then draw from \(h\) if so.

spike_slab_sample <- function(n, slab_sampler, pi0 = 0.9) {
  is_signal <- rbinom(n, 1, prob = 1 - pi0) == 1
  mu        <- numeric(n)
  mu[is_signal] <- slab_sampler(sum(is_signal))
  mu
}

End-to-end demo helper.

demo_prior <- function(name, slab_sampler, slab_pdf,
                       x_lim, mu_range, pi0 = 0.9, n = 2000) {
  mu <- spike_slab_sample(n, slab_sampler, pi0 = pi0)
  x  <- mu + rnorm(n)

  x_grid <- seq(x_lim[1], x_lim[2], length.out = 801)
  tw     <- tweedie_spike_slab(slab_pdf, x_grid,
                               pi0 = pi0, mu_range = mu_range)

  p_hist <- ggplot(data.frame(x = x), aes(x)) +
    geom_histogram(aes(y = after_stat(density)),
                   bins = 60, fill = "grey70", color = "white") +
    geom_line(data = tw, aes(x, f), color = "steelblue", linewidth = 1) +
    labs(title = sprintf("Marginal of X  (prior: 0.9·δ₀ + 0.1·%s)", name),
         x = "x", y = "density") +
    theme_minimal()

  p_scatter <- ggplot(data.frame(mu = mu, x = x), aes(x, mu)) +
    geom_point(alpha = 0.25, size = 0.9) +
    geom_abline(slope = 1, intercept = 0,
                linetype = "dashed", color = "grey50") +
    geom_hline(yintercept = 0, linetype = "dotted", color = "grey60") +
    geom_line(data = tw, aes(x, Emu),
              color = "firebrick", linewidth = 1) +
    coord_cartesian(xlim = x_lim) +
    labs(title = sprintf("Joint (X, μ) with E[μ | X]  —  slab = %s", name),
         x = "x", y = "μ",
         caption = "red = E[μ | X] via Tweedie;  dashed = identity") +
    theme_minimal()

  grid.arrange(p_hist, p_scatter, ncol = 2, widths = c(1, 1.2))
}

3. Three slabs, each with a 90% spike at 0

(i) Slab = \(\mathrm{Exp}(1)\)

A one-sided slab. For \(x\) near 0 or negative, the spike at 0 dominates the posterior and \(E[\mu\mid x] \approx 0\). Only sufficiently large positive \(x\) “escape” the spike.

demo_prior(
  name         = "Exp(1)",
  slab_sampler = function(n) rexp(n, rate = 1),
  slab_pdf     = function(mu) dexp(mu, rate = 1) * (mu >= 0),
  x_lim        = c(-4, 8),
  mu_range     = c(0, 30)
)

(ii) Slab = Logistic(0, 1)

Symmetric, slightly heavier than Gaussian. The shrinkage curve is now roughly flat near 0 (everything in the bulk gets pulled to the spike) and then bends up away from 0 — a smooth, soft-thresholding shape.

demo_prior(
  name         = "Logistic(0, 1)",
  slab_sampler = function(n) rlogis(n),
  slab_pdf     = function(mu) dlogis(mu),
  x_lim        = c(-10, 10),
  mu_range     = c(-40, 40)
)

(iii) Slab = \(t_2\)

Heavy-tailed slab. Strong shrinkage to 0 in the bulk, but the heavy tail means large \(|x|\) get believed almost verbatim — Tweedie’s curve approaches the identity line in the tails.

demo_prior(
  name         = "t_2",
  slab_sampler = function(n) rt(n, df = 2),
  slab_pdf     = function(mu) dt(mu, df = 2),
  x_lim        = c(-10, 10),
  mu_range     = c(-60, 60)
)

4. Why two different tail behaviors?

Across the three plots above, two distinct asymptotic shapes appear:

• Exp(1) and Logistic (exponential tails): the curve \(\hat{\mu}(x) = E[\mu\mid x]\) is parallel to the identity in the tail, \(\hat{\mu}(x) \sim x - b\) for some constant \(b > 0\) — a soft-thresholding rule.
• \(t_2\) (polynomial tail): the curve asymptotes to the identity, \(\hat{\mu}(x) \sim x\).

Here is why, working from Tweedie’s formula \(\hat{\mu}(x) = x + (\log f)'(x)\). The question is the behavior of \((\log f)'(x)\) for large \(x\). (For our spike-and-slab, the spike contributes \(\pi_0 \varphi(x)\) to \(f\), which decays super-exponentially and is negligible in the tail — the asymptotics are governed entirely by the slab.)

Tail of the marginal

Write \(f(x) \propto (h * \varphi)(x) = \int h(\mu)\,\varphi(x-\mu)\,d\mu\) in the tail, and substitute \(u = x - \mu\): \[ (h * \varphi)(x) \;=\; \int h(x-u)\,\varphi(u)\,du. \] The behavior of this integral depends on how fast \(h\) decays relative to the Gaussian \(\varphi\).

Case 1: exponential tail, \(h(\mu) \sim C e^{-\lambda \mu}\)

Then \(h(x - u) = C e^{-\lambda x} e^{\lambda u}\), and \[ (h * \varphi)(x) \;\sim\; C e^{-\lambda x} \int e^{\lambda u}\,\varphi(u)\,du \;=\; C e^{-\lambda x}\cdot e^{\lambda^2/2}. \] The slab is light enough that the Gaussian kernel does meaningful smoothing, but the result is still exponential. Taking logs and differentiating, \[ (\log f)'(x) \;\longrightarrow\; -\lambda, \qquad\text{so}\qquad \hat{\mu}(x) \;\sim\; x - \lambda. \]

• Exp(1) has \(\lambda = 1\), giving \(\hat{\mu}(x) \sim x - 1\).
• Logistic has \(h(\mu) \sim e^{-|\mu|}\) in both tails, also \(\lambda = 1\), giving \(\hat{\mu}(x) \sim x - 1\) for large \(x\) and \(\hat{\mu}(x) \sim x + 1\) for \(x \to -\infty\) — soft thresholding.

Case 2: polynomial tail, \(h(\mu) \sim C |\mu|^{-(\alpha+1)}\)

For \(t_\nu\), \(\alpha = \nu\); for \(t_2\), \(\alpha = 2\). Here the slab decays so much slower than the Gaussian that, in the tail, the integral is dominated by \(\mu \approx x\), where \(h\) is essentially constant on the \(O(1)\) scale of the Gaussian. To leading order, \[ (h * \varphi)(x) \;\sim\; h(x), \] i.e. the convolution just reproduces the slab. Then \(\log f(x) \sim -(\alpha + 1)\log|x|\), so \[ (\log f)'(x) \;\sim\; -\frac{\alpha + 1}{x} \;\longrightarrow\; 0, \qquad\text{so}\qquad \hat{\mu}(x) \;\sim\; x \;-\; \frac{\alpha+1}{x} \;\sim\; x. \]

The correction vanishes: heavy tails ⇒ believe the observation almost verbatim.

What the posterior looks like

A cleaner way to see this dichotomy: look directly at the posterior. The posterior does not concentrate to a point at large \(|x|\) — its variance stays of order 1, because that is the noise variance and the prior cannot override it. What changes is the posterior mean.

Exponential tail: posterior is exactly Gaussian

For \(h(\mu) \propto e^{-\lambda\mu}\), complete the square: \[ \log p(\mu \mid x) \;=\; -\tfrac12(x-\mu)^2 - \lambda\mu + \text{const} \;=\; -\tfrac12\bigl(\mu - (x - \lambda)\bigr)^2 + \text{const}. \] So \[ (\mu \mid X = x) \;\sim\; N\bigl(x - \lambda,\;1\bigr). \] The exponential prior simply shifts the Gaussian likelihood; the posterior mean is \(x - \lambda\) and the posterior variance is still 1. (For one-sided priors like Exp\((1)\) this is exact only for \(x\) large enough that the truncation at 0 is negligible; for two-sided exponential-tailed priors like Logistic, it’s the leading-order asymptotic.)

Polynomial tail: posterior is approximately N\((x,1)\)

For \(h(\mu) \sim C|\mu|^{-(\alpha+1)}\), write the posterior in terms of the deviation \(z = \mu - x\): \[ p(z \mid x) \;\propto\; h(x + z)\,\varphi(z). \] For large \(x\), \(h\) varies slowly on the \(O(1)\) scale where \(\varphi(z)\) has mass. Expanding, \[ \frac{h(x+z)}{h(x)} \;\approx\; 1 + (\log h)'(x)\,z + O(z^2/x^2) \;=\; 1 - \frac{\alpha+1}{x}\,z + \cdots \] So the posterior is a small tilt of \(N(0,1)\) — to leading order, \[ (\mu \mid X = x) \;\sim\; N(x,\,1), \] with the mean shifted by \((\log h)'(x) = -(\alpha+1)/x\). The prior is so flat in the tail that, locally, it adds essentially no information beyond what the likelihood already provides.

The dichotomy

Both cases give a posterior with variance 1. The difference is the mean:

• Exponential slab: mean is \(x - \lambda\) — a constant shift, fixed by the rate of decay.
• Polynomial slab: mean is \(x - (\alpha+1)/x\) — a vanishing correction; the data passes through unshrunk in the tail.

The competition is between the Gaussian likelihood and the slab: if the slab decays at least as fast as exponential, it “tugs” the posterior by a constant amount; if it decays only polynomially, the likelihood dominates and the posterior centers on the observation.

Numerical check

We can read these constants off the plots above. For example, for the Logistic slab at \(x = 10\), Tweedie gives:

Tail gap x − Ê[μ|x]: converges to a constant for exponential slabs, to 0 for polynomial slabs.

x	x − μ̂ (Exp1)	x − μ̂ (Logistic)	x − μ̂ (t₂)
5	1.005	0.953	0.648
10	1.000	1.000	0.306
15	1.000	1.000	0.202

For Exp(1) and Logistic, the gap \(x - \hat{\mu}(x)\) stabilizes near \(\lambda = 1\). For \(t_2\), it shrinks toward zero like \((\alpha+1)/x = 3/x\).

5. Takeaways

• \(E[\mu \mid X = x] = x + (\log f)'(x)\). Tweedie’s formula is unchanged for spike-and-slab — the spike just contributes \(\pi_0\,\varphi(x)\) to \(f\) exactly.
• The 90% spike at 0 induces soft thresholding: the posterior mean is approximately 0 in the bulk and only “turns on” once \(|x|\) is large enough that the slab outweighs the spike at that location.
• The shape of the slab controls how the curve bends past the threshold: Exp(1) gives a one-sided ramp; Logistic gives a smooth S; \(t_2\) approaches the identity in the tails (heavy tail ⇒ believe extremes).

References. Robbins (1956), An empirical Bayes approach to statistics; Efron (2011), Tweedie’s formula and selection bias, JASA.