EconBase
← Back to paper

Polynomial Log-Marginals and Tweedie's Formula : When Is Bayes Possible?

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

33,601 characters · 9 sections · 86 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Polynomial Log-Marginals and Tweedie's Formula When Is Bayes Possible?

abstractMotivated by Tweedie's formula for the Compound Decision problem robbins1951asymptotically, we examine the theoretical foundations of empirical Bayes estimators that directly model the marginal density $m(y)$. Our main result shows that polynomial log-marginals of degree $k \ge 3 $ cannot arise from any valid prior distribution in exponential family models, while quadratic forms correspond exactly to Gaussian priors. This provides theoretical justification for why certain empirical Bayes decision rules, while practically useful, do not correspond to any formal Bayes procedures. We also strengthen the diagnostic by showing that a marginal is a Gaussian convolution only if it extends to a bounded solution of the heat equation in a neighborhood of the smoothing parameter, beyond the convexity of $c(y)=\tfrac12 y^2+\log m(y)$.

\noindentKeywords: Tweedie's formula; empirical Bayes; normal means; Gaussian convolution; heat equation.

Empirical Bayes and compound decision

Consider an experiment in which unknown parameters, say $\boldsymbol{\mu} = (\mu_1, \ldots, \mu_n)$, gives rise to observations $Y_i \stackrel{\mathrm{ind}}{\sim} f(y \mid \mu_i)$, $i = 1, 2, \ldots, n$, and the goal is to estimate $\boldsymbol{\mu}$ under a loss function $\ell(y, \mu)$, e.g., the squared-error loss $n^{-1}\left|\left|\hat\boldsymbol{\mu}-\boldsymbol{\mu}\right|\right|^2$. We assume that the $Y_i$'s are exchangeable, i.e., invariant under permutation of the indices, which implies conditional independence and identical nature by de Finetti's theorem. For Gaussian $f(\cdot)$, this problem is known as the normal means, or Gaussian sequence, or the Gaussian compound decision problem. In the compound decision problem robbins1951asymptotically, the estimates (or any decision $\delta_i$ about $\mu_i$) is allowed to depend on all observations ${\bf y} = (y_1, \ldots, y_n)$, under the compound risk ${\mathbb E}_{\boldsymbol{\mu}}\{\ell(\delta, \boldsymbol{\mu})\}$. The compound decision problem, together with Stein's shrinkage phenomenon stein1956inadmissibility, shows that by allowing individual decisions/estimates to depend on the entire sequence, substantial reduction of total loss can be achieved. It is easy to see that the means $\mu_i$ are typically less extreme than the observations $y_i$; as Efron efron2011tweedie emphasizes, such regression to the mean is essentially the selection-bias or “winner's curse" phenomenon. The empirical Bayes (EB) methodology, introduced by robbins1956empirical, provides a statistical procedure where the goal of approximating the ideal or oracle Bayes rule is nearly achieved without specifying a prior zhang2003compound, efron2011tweedie. {A key appeal here is observability using Tweedie's formula that depends only on the marginal $m(y)$, which is estimable from data without committing to a prior. Tweedie's formula is not only a convenient empirical Bayes identity but also sits at the core of Stein’s unbiased risk estimate (SURE), highlighting how shrinkage and risk estimation meet in the normal means model.} As noted by zhang2003compound, the compound decision problem and the empirical Bayes methodology have been called the `two breakthroughs' in post-war statistics by neyman1962two, as they have influenced a lot of modern statistical methods, especially in the context of high-dimensional data. {In practice, EB includes Type II MLE/MMLE fits and related procedures that target sparse-signal regimes common in modern applications.}

One of the main attractions of the empirical Bayes framework for the compound decision problem is that it can result in risk reduction without solving the deconvolution problem of estimating the common conditional distribution $F$. Rather, one could work with the observable marginal density $m(y)$, and use either parametric or nonparametric models to estimate it, and use the estimated $\hat{m}(y)$ as a plug-in in the empirical Bayes approach. {This observable-only route computes estimators directly from $m(y)$, avoiding explicit priors and high-dimensional integration.} For example, assuming $\mu_i$'s to be normally distributed with zero mean and unknown variance $\sigma_0^2$, leads to a Gaussian marginal for the $Y_i$'s. Using a method of moments estimator for the $\sigma_0^2$ leads to the famous, classical James--Stein estimator $\delta(y) = (1 - (n-2)/\sum_{i=1}^{n} y_i^2) y$ stein1956inadmissibility. Difficulty in choosing a parametric model for the marginal, leads the way to more flexible nonparametric models such as kiefer1956consistency, which is easily seen to be convex and can benefit hugely from recent advances in convex optimization tools, as shown by koenker2014convex. Yet another approach is proposed by Efron efron2008microarrays, efron2012large, efron2011tweedie, using `Lindsey's method', which consists of estimating the marginal distribution from a Poisson regression of an integer order $K$, based on Lindsey's Method that models the histogram bin counts of the observed ${\bf y}$-values as Poisson random variables. efron2011tweedie's formulation yields locally adaptive shrinkage based on the local shape of the marginal histogram, while making no parametric assumptions on the prior. This has inspired subsequent developments in the empirical Bayes/compound decision problem by nonparametric density estimates of the marginal simon2013estimating, wager2014geometric. As emphasized by efron2014two, empirical Bayes admits two complementary strategies: f-modeling, which models the marginal density $m(\cdot)$ on the observation ($y$) scale, and g-modeling, which models the prior $F$ on the parameter ($\mu$) scale. The former have been relied upon heavily in applications, especially after robbins1956empirical, while the latter have been predominant in the theoretical EB literature laird1978nonparametric, jiang2009general. {Beyond EB, the Kiefer--Wolfowitz NPMLE can be viewed directly as a compound-decision rule; “no oracle needed” results due to ritov2024no show its Tweedie plug-in is minimax without reducing to an oracle Bayes benchmark.} {In high-dimensional prediction problems, it is also appealing to bypass the challenges of prior specification on a large parameter vector and instead work with the observable one-dimensional marginal $m(y)$ via Tweedie's formula. Such a `predictive Bayes' viewpoint leverages $m(y)$ directly for calibrated shrinkage and prediction without specifying a high-dimensional prior padilla2018deconvolution. In a recent discussion, efron2020pea contrasts prediction with attribution, arguing that prediction is inherently easier and more forgiving than attribution, and is often `anti-parsimonious' and `a-probabilistic,' with theoretical efficiency replaced by empirical methods. He also shows that a two-step hierarchical approach—using Tweedie's formula for the local false discovery rate in wide data settings (large $p$, small $n$) can drastically reduce dimensionality while avoiding explicit prior modeling and remaining transparent about the observable target. This motivates our focus on which functional forms of $m(\cdot)$ can arise from a bona fide Bayes model, and when the observable-only route, though appealing, might not always satisfy essential Bayes constraints.}

Given the popularity and performance of EB methods, a natural question from a Bayesian perspective is: for what class of prior distributions could a particular empirical-decision rule arise? {Our focus is precisely on which functional forms of the observable $m$ can actually occur as genuine Bayes marginals via Tweedie's identity.} In this article, we provide partial answers. First, if $\log m(y)$ is a polynomial of degree $>2$, no prior can produce $m(y)$ via Tweedie's formula; the only polynomial cases compatible with Bayes are degrees $\le 2$, with $K=2$ corresponding to a Gaussian prior. Second, beyond the convexity constraint on $c(y)=\tfrac12 y^2+\log m(y)$ as required by koenker2014convex, we point out that a marginal is a Gaussian convolution if and only if it extends (in a neighborhood of the smoothing parameter) to a bounded solution of the heat equation $u_t=u_{xx}$ by the Hirschman-Widder Weierstrass representation theorem widder1955convolution.

Tweedie's formula

In the simplest normal means model, let $\mu_i \stackrel{\text{iid}}{\sim} F$ and $Y_i \mid \mu_i \stackrel{\text{ind}}{\sim} \mathcal{N}(\mu_i,1)$. A remarkable result is the Tweedie's formula\footnote{It seems the nomenclature for Tweedie's formula is an interesting case of Stigler's law of Eponymy. In his famous paper, robbins1956empirical credits Tweedie, and efron2011tweedie follows, but, in her blog, Jiaying Gu documents a chain from Laird--Louis to Tukey to dyson1926method, and comments that `Tweedie's formula perhaps should be called Eddington's formula'. On the other hand, west1982aspects provides the formulae for posterior moments for both location and scale and credits masreliez1975approximate.}:

equation[equation omitted — 152 chars of source]

Here, $s_F(y)$ is the `Bayesian correction term' and $m(y)=\int \varphi(y-\mu)\,dF(\mu)$ is the marginal density of $y$. An attractive property of the Tweedie--Eddington formula (ref) is that the Bayes decision $\delta(\cdot)$ is expressed as a function of the observations ${\bf y}$, not involving the unknown prior $F$ efron2011tweedie, ritov2024no. efron2011tweedie discusses the Tweedie's formula for the general exponential family. Start with a canonical parameter $\eta$ and cumulant generating function $\lambda(\cdot)$, and density at $y = 0$ being denoted by $m_0$, we have: \[ y \mid \eta \sim f_{\eta}(y) = {\rm e}^{\eta y - \lambda(\eta)} m_0(y), \qquad \eta \sim g(\cdot). \] The posterior distribution under this model would be: $\pi(\eta \mid y) = f_{\eta}(y) g(\eta)/ m(y)$, where $m(y)$ is the marginal density given by: $m(y) = \int f_{\eta}(y) g(\eta) d \eta$, where the integral is carried over an appropriate sample space of the exponential family. robbins1956empirical showed that the posterior density for $\eta$ given $y$ is also a member of the exponential family with canonical parameter given by the observation $y$ and cumulant generating function (CGF) $\lambda(y) = \log(m(y)/m_0(y))$. Differentiating the CGF yields the posterior mean and variance as: \[ {\mathbb E}(\eta \mid y) = \lambda'(y), \quad {\rm Var}(\eta \mid y) = \lambda''(y). \] efron2011tweedie also discusses how robbins1956empirical's famous empirical Bayes Poisson prediction formula ${\mathbb E}(\lambda \mid y) = (y+1) m(y+1)/m(y)$ can be recovered as an approximation of the posterior mean derived from the Tweedie's formula applied to the exponential family representation of the Poisson probability mass function. brown2013poisson showed that this (now, more than) half-century old empirical Bayes formula can perform really well after some adjustments, viz., Rao-Blackwellization and smoothing (an isotonic regression to ensure montonicity). datta2016bayesian propose a global-local shrinkage prior for the Poisson compound decision problem that can also account for inflation of small counts, a phenomenon called `quasi-sparsity'.

polson1991representation extends the Tweedie's formula to a general location model. Let $\hat{\mu}$ and $a$ be the MLE and the maximal ancillary respectively and $y = (\hat{\mu}, a)$. Then the posterior mean ${\mathbb E}(\mu \mid y)$ for the location model under an arbitrary likelihood $f$ and a normal prior $\mu \sim \mathcal{N}(m, \tau^2)$ would be represented by the score function as:

equation[equation omitted — 144 chars of source]

where $p(\hat \mu \mid a ) = \int p(\hat \mu \mid \mu, a) dF(\mu)$. polson1991representation also shows that a similar score representation formula holds for a normal scale-mixture prior $\mu \sim \int \mathcal \varphi(\mu \mid m,\tau^2)\,dH(\tau^2)$, with $\tau^2$ in (ref) replaced by a suitable posterior quantity. polson2019bayesian provides the Tweedie's formula for normal linear regression: \[ {\bf y} = {\bf X}\boldsymbol{\beta} + {\rm e}, \qquad {\rm e} \sim \mathcal{N}(0,\Sigma). \] Let $g(\beta)$ denote the prior density of $\beta$, and define the marginal (prior predictive) density $m(y)=\int f(y \mid \beta)\,g(\beta)\,d\beta.$ Then, assuming $(X^\top\Sigma^{-1}X)^{-1}$ exists, the posterior mean of $\beta$ given $y$ is

equation[equation omitted — 209 chars of source]

The gradient of the prior predictive score $ \nabla_{{\bf y}} \log m({\bf y})$ is the `Bayesian correction' term in the context of linear regression. Clearly, (ref) reduces to (ref) for ${\bf X} = {\bf I}$, which reduces the regression problem to a normal means model.

Finally, west1982aspects provides us with the Tweedie's formula for the known location ($\mu = 0$), and unknown scale parameter $\sigma >0$, so that the likelihood is: $p(y \mid \sigma) = \sigma^{-1} p(\sigma^{-1}y)$. Define $\lambda = \sigma^{-2}$, then Theorem. 2.3.1 of west1982aspects says:

theoremLet the prior for the precision $\lambda>0$ be $\pi(\lambda)=\mathrm{Ga}\!\left(\tfrac{\alpha}{2},\tfrac{\beta}{2}\right)$ with density $\pi(\lambda)\propto \lambda^{\alpha/2-1}\exp\!\{-\tfrac{\beta}{2}\lambda\}$, and let \[ p_y(y)=\int_{0}^{\infty} \lambda^{1/2}\,p(\lambda^{1/2}y)\,\pi(\lambda)\,d\lambda, \qquad g_y(y)=-\frac{\partial}{\partial y}\log p_y(y),\quad G_y(y)=\frac{\partial}{\partial y}g_y(y). \] Then the posterior moments of $\lambda$ are \[ \mathbb{E}[\lambda\mid y]=\beta^{-1}\big\{(\alpha+1)-y\,g_y(y)\big\}, \qquad \operatorname{Var}(\lambda\mid y)=\beta^{-2}\big\{2(\alpha+1)-3y\,g_y(y)-y^{2}G_y(y)\big\}. \]

SURE: Tweedie's formula is also related to the Stein's unbiased risk estimation stein1981estimation framework for achieving a finite sample unbiased estimate of the prediction risk. For the normal means model, ${\bf y} \sim \mathcal{N}(\boldsymbol{\mu}, \sigma^2 I)$, with predictions denoted by $\hat{{\bf y}}$, the SURE for total prediction error is: $\left|\left|\hat{{\bf y}}-{\bf y}\right|\right|^2 + 2 \sigma^2 \sum_{i=1}^{n}\frac{\partial }{\partial y_i}{\hat{y}_i}$. bhadra2019prediction show that the SURE framework can be extended to regression by considering the singular value decomposition of the design matrix, and prove that the gain in predictive accuracy obtained by using the horseshoe prior can be achieved in finite samples, outperforming existing competitors, such as ridge regression and PCR. In a recent work, ghosh2025stein show that taking $h(y)=s_F(y)$ in $\delta(y)= y + h(y)$ yields \[ \mathrm{SURE}(\delta)=1+\frac{1}{n}\sum_{i=1}^n\Big\{ s_F(Y_i)^2 + 2\,\frac{\partial }{\partial y} s_F(Y_i) \Big\}. \] Choosing $F$ to minimize $\mathrm{SURE}(F)$ provides a data-driven approximation to the oracle Bayes rule. ghosh2025stein further provide a unified `$g$-modeling' efron2014two view showing that minimizing Stein's unbiased risk estimate is equivalent (up to a constant) to Hyv\'arinen's score matching hyvarinen2005estimation for learning the prior, and establishes near-parametric convergence rates for regret/Fisher divergence with oracle inequalities under misspecification. For a comprehensive discussion of empirical Bayes and the Tweedie--Eddington formula, we refer readers to the lecture notes by Ignatiadis and Sen ignatiadis_sen_empirical_bayes_2025.

Avoiding Deconvolution

All empirical Bayes methods, thus, share the goal of estimating the mixing measure $F$, irrespective of the loss function used gu2016problem. The deconvolution problem, i.e., estimating $F$, was referred to as `estimating the inestimable' by Robbins koenker2024empirical, and has a long history, as pointed out by ritov2024no. However, the typical reason is to provide a solution to the normal means or compound decision problem robbins1956empirical, i.e., a good estimate for the normal means $(\mu_1 , \ldots , \mu_n)$ something of interest in many different applications. For this problem, deconvolution is often an intermediate goal, but not the primary objective.

We do not attempt a full review of the vast literature on the normal-means problem. Instead, we focus on a specific class of empirical-Bayes estimators called empirical decision rules koenker2014convex. An empirical decision rule attempts to estimate the posterior mean function ${\mathbb E}(\mu_i \mid y_i) $ directly. Such a rule implicitly invokes the existence of some prior distribution $f_0$ but does not attempt to estimate it efron2011tweedie.

To understand this empirical-Bayes approach, return to the simple normal-mean model $(y_i \mid \mu_i) \sim \mathcal{N}(\mu_i,1)$. Absent any prior information that would distinguish the $\mu_i $, a natural assumption is that they arise exchangeably from some prior, $ \mu_i \sim F$. An important result discussed by robbins1956empirical holds that the posterior mean can be written as the unbiased estimate plus a Bayes correction:

equation[equation omitted — 91 chars of source]

where $m(y) = \int \varphi(y \mid \mu) dF(\mu)$ is the marginal density of the data under this prior. The prior does not appear explicitly, but its effect is incorporated in $m(y)$. efron2011tweedie nicely summarizes its appeal: “{The crucial advantage of Tweedie's formula is that it works directly with the marginal density,}" thereby offering the statistician the opportunity to avoid deconvolution entirely.

Nonparametric Maximum Likelihood

Among methods that directly estimate the marginal, brown2009nonparametric advocate for replacing $m(y)$ by a kernel density estimate of it, while the usage of a nonparametric maximum likelihood estimator kiefer1956consistency for $m(y)$ has been investigated by various authors including jiang2009general, saha2020nonparametric, greenshtein2022generalized. koenker2014convex argue that recent advances in convex optimization has paved the way for wider usage and applicability of the Kiefer--Wolfowitz algorithm for the nonparametric empirical Bayes compound decision problems. Ritov ritov2024no shows that the Tweedie plug-in based on the nonparametric MLE of the marginal (equivalently, the Kiefer-Wolfowitz NPMLE for the mixing law) yields a decision rule (based on the Tweedie's formula) that is minimax in both the empirical Bayes and compound-decision formulations, establishes concentration for the NPMLE decision, and links its risk directly to SURE - thereby providing a justification for NPMLE without oracle arguments.

In other words, the key insight here is that, since we observe $m(y)$ directly, a non-parametric approach to estimate $m(y)$ and use that as a plug-in in Tweedie's formula is a natural solution. Moreover, since (ref) only involves the marginal $m(y)$, one can bypass priors entirely and estimate: \[ \hat\mu_i^{\text{TW}} \;=\; y_i + \partial_y \log \hat m(y_i), \] where $\hat m$ is any nonparametric estimate of $m$ (e.g., a kernel density, a Kiefer--Wolfowitz NPMLE Gaussian mixture, or a shape-constrained fit), thus avoiding deconvolution and explicit prior modeling brown2009nonparametric, kiefer1956consistency, koenker2014convex, jiang2009general, saha2020nonparametric, greenshtein2022generalized. Writing $\ell(y)=\log m(y)$ and $s(y)=\ell'(y)$, one may estimate $\ell$ by minimizing a score-matching/SURE objective with a smoothness penalty, \[ \min_{\ell}\;\Big\{\;\frac{1}{n}\sum_{i=1}^{n}\big[s(y_i)^2+2\,s'(y_i)\big] \;+\; \rho\,\mathcal R(\ell)\;\Big\}, \qquad \mathcal R(\ell)=\int\big(\ell''(y)\big)^2\,dy, \] optionally under shape constraints such as convexity of $c(y)=\tfrac12 y^2+\ell(y)$ on a grid (linear inequalities), which yields a convex quadratic program and directly returns the Tweedie plug-in $\hat\mu_i^{\text{TW}}$ ghosh2025stein.

Lindsey's method

efron2011tweedie argues that the Tweedie's formula $y + \frac{\partial }{\partial y} \log m(y)$ requires a smoothly differentiable estimate of the log marginal $\log m(y)$, and describes a Poisson regression technique, called the `Lindsey's method' efron2004estimation, efron2011tweedie and efron2012large. The Lindsey's method assumes that the log-marginal is a $K$-th degree polynomial, for some integer $K$. That is, it estimates the $\mu_i $ by first estimating the marginal density directly using a polynomial or spline on the log scale, e.g.,

equation[equation omitted — 105 chars of source]

The $\beta_k$'s are estimated via a Poisson regression (using usual generalized linear model tools) to the histogram counts after binning the data. Under this form of the marginal density, the posterior mean function is easily computed by taking the log and differentiating. Following on the work of efron2011tweedie, koenker2014convex considered an estimator based on the observation that the posterior mean function ${\mathbb E}(\mu \mid y) $ in Tweedie's formula (ref) is non-decreasing in $y$. This implies that the function \[ c(y)= \frac{1}{2} y^2 + \log m(y) \] is convex, a result which holds not merely for the Gaussian case, but for any exponential-family observational model, regardless of the prior $F$. However, an unconstrained or traditional kernel-based estimate of $m(y)$ does not ensure convexity of $c(y)$. As discussed before, koenker2014convex proposed to estimate $m(y)$ non-parametrically.

A Bayesian view of empirical decision rules

From a Bayesian perspective, a natural question is: for what class of prior distributions could a particular empirical-decision rule arise? In the case of Efron's estimator, we provide an explicit answer in the following theorem: there are no such priors, save the trivial case of a Gaussian prior (and a quadratic log marginal). We state this as the following theorem:

theoremLet $\varphi(y\mid\mu)=\exp\{\mu y-\tfrac12\mu^2\}\,\varphi(y)$ be the $\mathcal{N}(\mu,1)$ exponential-family form, where $\varphi$ is the $\mathcal{N}(0,1)$ density, and let $m(y)=\int \varphi(y\mid\mu)\,dF(\mu)$ be the marginal induced by a prior $F$ on $\mu\in\mathbb{R}$. If $\log m(y)$ is a polynomial of degree $K\ge 3$ on $\mathbb{R}$, then there is no such prior $F$. In other words, for all prior measures $F$ supported on $\mathbb{R}$, and for all nonzero choices of $\beta_1 , \ldots , \beta_K $, $$ \int_{- \infty}^\infty \phi( y \mid \mu ) ~d F(\mu ) \neq \exp \left ( {\beta}_0 + \sum_{k=1}^K {\beta}_k y^k \right ), \quad \text{if} \quad K \geq 3. $$ If $K=2$, then $F$ must be (possibly degenerate) Gaussian.

Hence, Efron's empirical Bayes rule is therefore not a formal Bayes rule, in that it cannot arise from a valid application of Tweedie's formula under any prior. The only exception is $K = 2$, implying that $F$ is Gaussian. Interestingly, when $\log m(y)$ is quadratic ($K=2$), Tweedie's formula gives ${\mathbb E}(\mu_i\mid Y_i = y)=(1-1/V)y$ for a normal marginal $\mathcal{N}(0,V)$, and replacing $1/V$ by the unbiased estimator $(n-2)/\sum_{j=1}^n y_j^2$ yields the James--Stein estimator efron2011tweedie. We note that the argument extends trivially to other smooth bases, such as spline expansions, since these do not define entire functions on the complex plane and therefore cannot yield a valid prior either.

koenker2014convex cast two EB estimators as convex programs: (i) a direct MLE of the marginal density $m$ under the shape constraint that $c(y)=\tfrac12 y^2+\log m(y)$ is convex : this enforces monotonicity of $\delta(y)=y+m'(y)/m(y)=c'(y)$ and yields a piecewise-linear $c$, so $\delta(y)=c'(y)$ is piecewise constant; and (ii) a Kiefer--Wolfowitz NPMLE of the mixing law $F$ via a finite-dimensional convex dual whose solution is atomic (with at most $n$ support points), with interior-point implementations delivering substantial speedups over EM while retaining excellent risk performance. For the koenker2014convex estimator, all Gaussian convolutions do satisfy the convexity of $c$, but convexity is only necessary, not sufficient. Additional conditions are given by the following result adapted from Theorem VIII.6.3 of widder1955convolution.

We note that koenker2014convex themselves emphasize the practical advantages of the weaker convexity constraint, viewing it as a deliberate relaxation rather than a deficiency\footnote{“(S)ince the monotonized Bayes rule imposes a weaker form of convexity constraint, it may have some advantages in misspecified settings where the original Gaussian location mixture assumption fails to be satisfied." koenker2014convex}. Theorem (ref) identifies precisely when such convexity-based rules cease to correspond to any Bayesian model.

theorem[widder1955convolution] Let $ u: \mathbb{R} \times(0,\infty) \mapsto \mathbb{R}$ be $C^{2}$ in $x$, $C^{1}$ in $t$, and bounded for each $t>0$. Then the following are equivalent: \begin{enumerate} • $u$ solves the heat equation $u_t=u_{xx}$ on $\mathbb{R}\times(0,\infty)$ and, for some bounded Borel measure $\mu$, \[ u(x,t)=\frac{1}{\sqrt{4\pi t}}\int_{\mathbb{R}} e^{-(x-y)^2/(4t)}\,d\mu(y)\qquad(t>0). \] • For each fixed $t>0$, $u(\cdot,t)$ is the Weierstrass transform (Gaussian convolution) of a bounded function; in particular $m_t(x):=u(x,t)$ is a Gaussian convolution. \end{enumerate} If, in addition, $\mu$ has a bounded density $f$, then $u(\cdot,t)=\phi_t*f$ with $\phi_t$ the $N(0,2t)$ kernel. Consequently, for a single marginal $m(x)=u(x,t_0)$ (e.g., $t_0=1$), being a Gaussian convolution is stronger than having $c(x):=\tfrac12 x^2+\log m(x)$ convex: $m(\cdot)$ must extend to a bounded solution of the heat equation in a time neighborhood of $t_0$.
remark{One-parameter antecedents.} For one-parameter models, early work showed important rigidity phenomena: for a class where the sample total is sufficient for any sample size, the fiducial distribution of the parameter cannot coincide with any Bayesian posterior under any prior grundy1956fiducial. In one-parameter natural exponential families, sampson1975characterizing shows that the moment generating function (as a function of the natural parameter) uniquely characterizes the family and provides necessary and sufficient conditions for such MGFs. In one-parameter settings where the posterior expectation of the population mean is a linear function of the sample observation, goldstein1975uniqueness proves that the moments of the prior distribution are uniquely determined, giving precise uniqueness relations for such linear Bayes rules. Finally, in a one-parameter exponential family under squared-error loss, goldstein1977contractions shows that the only Bayes estimators whose 1-Lipschitz (contraction) transforms are also Bayes are the linear (affine) estimators in the canonical sufficient statistic.

Proof of Theorem (ref)

proofFirst, note that for any probability prior $F$, the mixture $m(y)=\int\varphi(y-\mu)\,dF(\mu)$ is strictly positive on $\mathbb{R}$; thus $\log m(\cdot)$ is well-defined and finite on $\mathbb{R}$. Next, writing \[ m(y)=\varphi(y)\int e^{\mu y-\mu^2/2}\,dF(\mu)=\varphi(y)\,M_H(y), \] where $H$ is the finite positive measure $dH(\mu)=e^{-\mu^2/2}\,dF(\mu)$ and $M_H(z)=\int e^{z\mu}\,dH(\mu)$ is its moment generating function. Then, for $z=x+it\in\mathbb{C}$, \[ |M_H(z)| \le \int e^{x\mu}e^{-\mu^2/2}\,dF(\mu) = \int \exp\!\Bigl(-\tfrac12(\mu-x)^2+\tfrac12 x^2\Bigr)dF(\mu) \le e^{x^2/2}. \] Hence $M_H$ is entire and of at most Gaussian growth in $\Re z$. Put $M:=H(\mathbb{R})=\int e^{-\mu^2/2}dF(\mu)\in(0,1]$ \footnote{Since \(F\) is a probability measure and \(0<e^{-\mu^{2}/2}\le 1\) for all \(\mu\in\mathbb{R}\), we have: $M:=H(\mathbb{R})=\int_{\mathbb{R}} e^{-\mu^{2}/2}\,dF(\mu)\in(0,1].$} If $\log m(y)$ is a polynomial on $\mathbb{R}$, then \[ \log M_H(y)=\log m(y)-\log\varphi(y)=:P(y) \] is also a polynomial (since, $-\log\varphi(y)=\tfrac12 y^2+\mathrm{const}$ is quadratic). Both $M_H(z)$ and $e^{P(z)}$ are entire and agree for all real $z$; by the identity theorem for entire functions, $M_H(z)\equiv e^{P(z)}$ on $\mathbb{C}$ (see, e.g., Ahlfors1979). Evaluating on the imaginary axis and normalizing yields a characteristic function: \[ \phi_{\widetilde H}(t):=\frac{M_H(it)}{M}=\exp\!\bigl(P(it)-\log M\bigr)=:\exp Q(t), \] so $\log\phi_{\widetilde H}$ is a polynomial. By Marcinkiewicz's theorem, if a characteristic function has the form $\phi(t)=\exp\{Q(t)\}$ with $Q$ a polynomial on $\mathbb{R}$, then $\deg Q\le 2$. Moreover, $\deg Q=0,1$ correspond to a degenerate distribution and $\deg Q=2$ to a (nondegenerate) Gaussian distribution Bryc1995; see also Lukacs1970 and Marcinkiewicz1939. Therefore $\deg P\le 2$, so $\deg\log m(\cdot)\le 2$. In particular, no prior $F$ can produce a marginal with polynomial $\log m$ of degree $K\ge 3$. (Marcinkiewicz's theorem applies here since $\phi_{\tilde H}(t)=M_H(it)/M$ is a bona fide characteristic function.) If $\log M_H(y)=\alpha y^2+\beta y+\gamma$, then \[ m(y)=\varphi(y)\,e^{\alpha y^2+\beta y+\gamma} = C\exp\!\Bigl(-\tfrac12(1-2\alpha)\,y^2+\beta y\Bigr), \] with $C=e^\gamma/\sqrt{2\pi}$. Since $m$ is integrable, $1-2\alpha>0$, so $m$ is a (nondegenerate) Gaussian density. Write $Y=\mu+\varepsilon$ with $\varepsilon\sim N(0,1)$ independent of $\mu\sim F$. Then $Y$ is normal and $\varepsilon$ is nondegenerate normal. By Cram\'er's decomposition theorem, $\mu$ must be (possibly degenerate) normal Bryc1995; see also Cramer1936. Hence, $F$ is Gaussian.
remarkThe argument extends verbatim to any LEF (linear exponential family) with quadratic carrier log-density, since Step 2 only uses that $\log f_0(y)$ is quadratic to pass from $\log m$ polynomial to $\log M_H$ polynomial. For general LEFs this implication is even more restrictive. Background on characteristic functions and cumulants (used implicitly in Marcinkiewicz/Cramer arguments) is collected in Feller1971.

Conclusion

Why care whether an empirical Bayes estimator corresponds to a formal Bayes procedure, that is, whether it can arise from a hierarchical model with an actual prior? Proper Bayes rules bring structural guarantees: admissibility under a proper prior with finite Bayes risk, coherence of the posterior, and, in favorable cases, least-favorable-prior minimaxity. A standard result (see Lehmann and Casella) states: let $R(\theta,\delta)={\mathbb E}_\theta L(\theta,\delta(X))$ be the risk and $r(\delta,\Pi)=\int R(\theta,\delta)\,d\Pi(\theta)$ the Bayes risk for a proper prior $\Pi$; if $\delta_\Pi$ is Bayes with respect to $\Pi$ and $r(\delta_\Pi,\Pi)<\infty$, then $\delta_\Pi$ is admissible. Admissibility itself is a weak filter and does not, by itself, recommend a procedure: as strawderman2021charles notes, “Admissibility is a relatively weak optimality property, and the fact that a procedure is admissible does not give strong evidence that it should be adopted. Typically, there are a very large number of admissible procedures in any given problem, including all procedures which are unique Bayes with respect to any proper prior." Our results mark the boundary for Tweedie-based rules: polynomial $\log m$ of degree greater than two cannot be Bayes, the quadratic case corresponds to a Gaussian prior, and a marginal is a Gaussian convolution only if it extends locally in the smoothing parameter to a bounded solution of the heat equation. We hope to have shown that empirical Bayes and NPMLE plug-ins optimize observable criteria such as SURE and regret, but they do not automatically inherit the guarantees that come with a proper Bayes formulation.

Acknowledgments

The first author (JD) gratefully acknowledges support from the National Science Foundation (NSF CAREER Award DMS-2443282).