M1.4 Denoising as a Posterior Mean: Tweedie's Formula
Builds on3 notes
The scale reads 72.3
A randomly selected student’s weight follows a known class distribution, approximately Gaussian with mean 62 kg and standard deviation 8 kg. An unbiased scale adds independent Gaussian noise of standard deviation 1 kg.
The displayed reading is 72.3 kg.
What is the best squared-error estimate of the actual weight?
Symmetric noise makes readings unbiased given a true value. Reversing the conditioning still requires the prior: different clean origins occur at different rates.
Could the observed-reading distribution provide the posterior mean without explicitly evaluating every posterior weight?
DiscussionQ1
What is random, which loss defines the best estimate, and how does the population enter?
Show the reasoning
Zero-mean noise identifies , not the reversed expectation we need.
Random selection gives a clean value with prior . The reading is , with independent standard-Gaussian .
Squared-error prediction selects , the average clean origin among readings compatible with .
Near the high reading, lighter compatible origins may be more common than heavier ones. Posterior reweighting then shifts the estimate toward the population center.
Stand on the slope of a density
Imagine people randomly displaced left or right. At the observed position, their possible origins must be weighted by both displacement probability and original population density.
For the Gaussian class example, more compatible people originate on the lighter side, shifting the posterior center below the displayed reading.
The exact correction uses the observed density , not the clean density . Noise has already smoothed the population before this slope is calculated.
Figure: the correction is the observed score multiplied by noise variance.
Connect three tools to derive Tweedie’s identity
Assume independent Gaussian additive noise with , existing conditional means, and conditions permitting differentiation under the integral. The observed density is
Differentiate in three steps.
First, differentiate under the integral over clean origins:
Second, differentiate the Gaussian kernel. Its derivative is the kernel times :
Third, divide by the marginal and use Bayes’ rule. The normalized integrand becomes posterior weighting:
Rearranging gives Tweedie’s formula:
For isotropic vector noise, replace the scalar derivative by a gradient.
The estimate equals observation plus noise variance times the observed log-density slope. Its direction, magnitude, and noise dependence each have a separate role.
DetailsWhy the observed marginal and its score are useful
The identity needs the observed marginal score without separately reconstructing the clean prior. This connects to empirical Bayes [1,2]. Estimating a score from data still introduces error; a noisy histogram is not an exact score.
A positive constant normalizer cancels from the log gradient. Shape accuracy, particularly derivative accuracy, remains necessary.
Calculate the correction
For and unit noise, the observed variance is 65 and its score is . Thus
With noise standard deviation 3, observed variance becomes 73. The correction is kg, giving about 71.03 kg.
Zero marginal score makes the correction zero at that point, including a symmetric center. Increasing noise also changes the observed score, so the total correction need not scale solely as .
The identity assumes the correct independent Gaussian noise model and scale, the relevant source population, and squared-error loss for posterior-mean optimality. Using another population’s observed density can give the wrong inference.
Simulate noisy observations
import numpy as np
rng = np.random.default_rng(4); n = 2_000_000
def run(noise, sigma=1.0, ys=(66.0, 72.3, 78.0)):
X = rng.normal(62, 8, n)
if noise == "gauss": eps = rng.normal(0, 1, n)
else: eps = np.where(rng.random(n) < 0.95, rng.normal(0, 1, n), rng.normal(0, 10, n))
Y = X + sigma*eps
emp = [X[np.abs(Y-y)<0.1].mean() for y in ys]
hist, edges = np.histogram(Y, bins=350, range=(30, 100), density=True)
c = 0.5*(edges[1:] + edges[:-1])
dlogp = np.gradient(np.log(hist + 1e-12), c)
estimate = np.array(ys) + sigma**2*np.interp(ys, c, dlogp)
print(noise, sigma, "conditional bins", np.round(emp, 2), "score estimate", np.round(estimate, 2))
run("gauss", 1.0); run("gauss", 3.0); run("glitch", 1.0)
Gaussian cases should approximately match posterior bin averages, with finite-sample and histogram-derivative error. The glitch mixture violates the single-Gaussian noise assumption, so the same correction need not match its true posterior mean.
Check your understanding
References
- Robbins, H. An Empirical Bayes Approach to Statistics. Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Vol. 1, 1956.
- Efron, B. Tweedie’s Formula and Selection Bias. Journal of the American Statistical Association 2011.
- Stein, C. M. Estimation of the Mean of a Multivariate Normal Distribution. Annals of Statistics 1981.