MoreRSS

site iconJohn D. CookModify

I have decades of consulting experience helping companies solve complex problems involving applied math, statistics, and data privacy.
Please copy the RSS to your reader, or quickly subscribe to:

Inoreader Feedly Follow Feedbin Local Reader

Rss preview of Blog of John D. Cook

Posterior variance

2026-07-13 03:26:29

A few days ago I wrote a post entitled Does additional data always reduce posterior variance?. In a nutshell, the answer is no, not always.

That led the previous post which looked at posterior means for three Bayesian models, showing how the posterior mean is a weighted average of the prior mean and the mean of the new data. The weights are precisions, which means something different for each model.

For the beta-binomial model, variance may increase when seeing unexpected data (details here), but precision always increases.

For the normal-normal model precision is the reciprocal of variance. Every new data point makes precision go up and posterior variance go down.

The Poisson-gamma model may be the most interesting. As stated in the previous post, if data has a Poisson distribution with parameter λ, and λ has a gamma(α0, β0) prior distribution, then the posterior distribution on λ after observing k events over time t has a gamma(α0 + k, β0 + t) posterior distribution. Therefore the posterior variance is

0 + k) / (β0 + t)².

Note the posterior variance is an increasing function of k and a decreasing function of t. This means that the posterior variance increases every time an event is observed, and it decreases quadratically between observations.

Here’s an illustration. I simulated data from a Poisson process with λ and used a gamma(1, 1) prior on λ. Here’s a plot of the posterior variance.

The post Posterior variance first appeared on John D. Cook.

Posterior mean

2026-07-13 01:48:14

Common sense says that what you believe after seeing new data should be some sort of compromise between what you believed before and what the new data says. You don’t want to ignore previous information or new information.

How much should new data change your prior beliefs? When prior judgment and new information are in conflict, which one should be given the benefit of the doubt?

Bayesian data models provide a framework for making such decisions quantitative and objective. The choice of a data model is somewhat subjective—whether it’s a Bayesian model or not—but given a Bayesian model, the rules for updating the representation of your beliefs are objective. As some put it, you “turn the Bayesian crank.” A likelihood model and a prior on parameters together specify how new data changes the prior distribution into a posterior distribution.

We will make this more concrete with three examples.

Normal-normal model

Suppose that data X has a normal distribution with unknown mean μ and known variance σ², and we assume that a priori μ has a normal distribution with mean μ0 and variance σ0².

After observing x, the posterior distribution on μ also has a normal distribution, but with a different mean and variance. Its mean is somewhere between the prior mean and x. We will ignore the change in the variance for this post.

The posterior mean of μ is

\mu_{\text{post}} = \frac{\dfrac{\mu_0}{\sigma_0^2} + \dfrac{x}{\sigma^2}}{\dfrac{1}{\sigma_0^2} + \dfrac{1}{\sigma^2}}

This equation becomes more understandable when we introduce precisons τ = 1/σ² and τ0 = 1/σ0².

Then we have

\mu_{\text{post}} = \frac{\mu_0 \cdot \tau_0 + x \cdot \tau}{\tau_0 + \tau}

which you can read as saying the posterior mean is the weighted average of the prior mean and x, with the weights given by the precision. Intuitively, you take the weighted mean of your conclusions from previous data and new data, weighting the mean according to how much confidence you have in each.

Beta-binomial model

Now let’s switch over to a different data model. Now assume X is a binary random variable, with probability of success p and probability of failure 1 − p, and we assume p has a beta(a, b) distribution.

After observing s successes and f failures, the posterior mean of the distribution on p becomes

p_{\text{post}} = \frac{a + s}{a + b + s + f}

We can rewrite this as

p_{\text{post}} = \frac{(a + b) \dfrac{a}{a+b} + (s + f) \dfrac{s}{s+f}}{(a + b) + (s +f)}

This says that the posterior mean is the weighted average of the prior mean a/(a + b) and the mean of the data s/n. The weights are the prior effective sample size a + b and the sample size of the new data n. In this example (effective) sample size is playing the role that precision played in the normal-normal model above.

Gamma-Poisson model

Suppose data have a Poisson distribution with parameter λ, and λ has a gamma(α0, β0) prior distribution [1]. And suppose you observe k events over time t. Then the posterior distribution of λ given the data has a gamma(α0 + k, β0 + t) prior distribution and the mean of the posterior distribution is given by

\lambda_{\text{post}} = \frac{\alpha_0 + k}{\beta_0 + t} = \frac{\beta_0 (\alpha_0 / \beta_0) + t (k / t)}{\beta_0 + t}

As before, the posterior mean is a weighted average of the prior mean and new data, and the weights are interpretable as some sort of measure of confidence, namely time. The variable t is directly time and the parameter β0 is sort of an effective time, just as a + b is an effective sample size for the beta distribution.

Common thread

In each example the posterior mean is the weighted average of the prior mean and the mean of the data, with the weights given by a precision. However, precision means something different in each example. In the normal-normal model, precision is the reciprocal of variance, but in the beta-binomial model precision is sample size and in the Poisson-gamma model precision is time.

What all three examples have in common is that they are conjugate models using distributions from the “exponential family” of probability distributions. In technical terms, precision is the multiplicative factor on the sufficient statistic in the exponent of the posterior kernel.

Related posts

[1] There are multiple conventions for parameterizing the gamma distribution. Here we’re using the shape-rate parameterization, where the mean is α/β.

The post Posterior mean first appeared on John D. Cook.

Progress on Gilbreath’s conjecture

2026-07-12 05:30:11

Years ago I wrote about Gilbreath’s conjecture. It’s a simple conjecture; you could explain it to anyone who understands what prime numbers are. See the linked post for a description of the problem.

Gilbreath’s conjecture is simple, but it’s also kinda weird. As I wrote before,

Paul Erdős speculated that Gilbreath’s conjecture is true but it would be 200 years before anyone could prove it. I find Erdős’s conjecture more interesting than Gilbreath’s conjecture.

The conjecture is hard in a way that, say, solving a nasty-looking differential equation is not. Over the last three centuries, mathematics has developed quite a toolbox for solving differential equations. But Gilbreath’s conjecture is just odd enough that it’s not at all clear what kind of tool might be useful in approaching it.

Terence Tao has a new blog post announcing a paper that he and two coauthors wrote on a random model intended to mimic Gilbreath’s calculation on primes. This random model is more sophisticated than the little game Gilbreath was playing, but it’s also much more amenable to analysis by established techniques. Tao’s post gives a heuristic explanation for why Gilbreath’s conjecture is plausible, but then adds

However, it seems well beyond current technology to try to make these heuristics rigorous; even the first step … is far out of reach.

Related posts

The post Progress on Gilbreath’s conjecture first appeared on John D. Cook.

Reproducing a geometry theorem diagram

2026-07-06 22:22:35

I ran across a geometry theorem with the following diagram.

The theorem corresponding to the diagram is interesting, but I found reproducing the diagram more interesting.

The segment AB is a diameter and the line CD is perpendicular to the diameter.

Assume the outer circle is a unit circle. I guessed C = (cos(1), sin(1)) and made the following diagram.

I guessed the value of C by eyeballing it, but in retrospect this would have been a convenient value for the creator of the original diagram to have chosen.

Drawing the blue circle inscribed in the triangle was easy using the equations for the center and radius from this post. Drawing the other two circles, the green and orange circles, was harder. They are also inscribed circles, but not inscribed in a triangle. They’re inscribed in a three-sided figure with two perpendicular sides and a circular arc.

The radius r of the green circle is the distance from the center of the circle to each of its tangent lines. Also, the distance from the origin to the center of the circle must be 1 − r. This is enough information to set up a quadratic equation for r. The same reasoning applies to the orange circle.

The original diagram comes from [1] and the theorem it illustrates says the diameter of the blue circle equals the sum of the radii of the green and orange circles.

Python code

In case you’re interested, here’s the code that created the diagram.

#!/usr/bin/env -S uv run --script

# /// script
# dependencies = ["numpy", "matplotlib"]
# ///

import numpy as np
import matplotlib.pyplot as plt

def connect(A, B, color='gray'):
    plt.plot([A[0], B[0]], [A[1], B[1]], color=color, linewidth=2)

def circle(c, r, color='gray'):
    t = np.linspace(0, 2*np.pi)
    plt.plot(c[0] + r*np.cos(t), c[1] + r*np.sin(t), color=color, linewidth=2)

def quadratic(a, b, c):
    det = b**2 - 4*a*c
    return ((-b - det**0.5)/(2*a), (-b + det**0.5)/(2*a))

A = np.array([-1, 0])
B = np.array([ 1, 0])
C = np.array([np.cos(1), np.sin(1)])
a = np.linalg.norm(B - C)
b = np.linalg.norm(A - C)
c = np.linalg.norm(B - A)
s = (a + b + c)/2

circle([0,0], 1)
connect(A, B,)
connect(A, C)
connect(C, B)
connect(C, C*np.array([1, -1]))

center = (a*A + b*B + c*C)/(2*s)
radius = 0.5*a*b/s
circle(center, radius, 'C0')

Ex = C[0]
roots = quadratic(1, 2 + 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
print(roots)
center = (r + Ex, -r)
circle(center, r, 'C1')

roots = quadratic(1, 2 - 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
center = (Ex - r, -r)
circle(center, r, 'C2')

plt.gca().set_aspect("equal")
plt.axis("off")
plt.show()

[1] Leon Bankoff. A Geometrical Coincidence. Mathematics Magazine, Vol. 37, No. 5 (Nov., 1964), p. 324.

The post Reproducing a geometry theorem diagram first appeared on John D. Cook.

e approximation

2026-07-06 20:22:29

I ran across the approximation

e ≈ 2721/1001

recently. What makes this remarkable is its accuracy relative to the size of the denominator.

You can create a trivial approximation just by truncating a decimal expansion

e ≈ 2718/1000

but this is only good to four significant figures, but 2721/1001 is good to seven, almost eight, significant figures.

e         = 2.71828182… 
2721/1001 = 2.71828171…

The comparison is more impressive in binary.

$ bc -l
>>> obase=2
>>> 2721/1001
10.10110111111000010100…
>>> e(1)
10.10110111111000010101…

The denominator is a 10-bit number but the approximation is accurate to 21 bits.

The post e approximation first appeared on John D. Cook.

Does additional data always reduce posterior variance?

2026-07-04 10:50:36

A discussion over lunch today brought up the fact that additional data does not always decrease the size of a confidence interval. This post will look at this from a Bayesian perspective.

In general, new information reduces your uncertainty regarding whatever you’re estimating. The posterior distribution becomes more concentrated as more data are collected.

That’s what happens “in general” but does it necessarily happen every time you get new data? Conceivably if you get surprising data, data that is very unlikely given your current prior, posterior uncertainty might increase.

Beta-binomial model

To show that this is the case, suppose the probability of success in some binary trial has parameter θ and that θ has a beta prior. You could imagine this prior to be the posterior after having made some number of previous observations. Can a new observation increase the posterior variance in θ? If so, under what conditions?

The variance of a beta(a, b) random variable is

ab / (a + b)²(a + b + 1).

After observing a successful trial, the posterior distribution on θ is beta(a + 1, b). We can calculate the ratio of the posterior variance to the prior variance and ask under what circumstances, if any, the ratio is greater than 1.

If 2ab the posterior variance will be strictly less than the prior variance. This says if the prior mean odds against a success are no more than 2 : 1, observing a success will reduce the variance. (So will observing a failure.) But for any value of b, you can find a small enough value of a that observing a success will increase the variance.

Normal-normal model

Whether an observation can increase the posterior variance depends on the data model. If your data have a normal likelihood function with known variance and a normal prior on the mean θ, the posterior variance is always less than the prior observation, and it reduces by the same amount, independent of the observation x. If x is very unlikely a priori then it will pull the posterior mean toward itself more than an observation that is more concordant with the prior would have, but the change in the posterior variance is the same.

Proof of beta theorem

Here is a proof in Lean 4 of the statement above that if 2ab the posterior variance will be strictly less than the prior variance.

import Mathlib

set_option linter.style.header false

noncomputable def f (a b : ℝ) : ℝ := a * b / ((a + b) ^ 2 * (a + b + 1))

theorem f_ratio_lt_one' (a b : ℝ) (ha : 0 < a) (hb : 0 < b) (hab : b ≤ 2 * a) :
    f (a + 1) b / f a b < 1 := by
  have hs : 0 < a + b := by linarith
  have h2ab : 0 ≤ 2 * a - b := by linarith
  have hprod : 0 ≤ (a + b) * (2 * a - b) := mul_nonneg hs.le h2ab
  -- key polynomial inequality (∗)
  have key : (a + 1) * (a + b) ^ 2 < a * ((a + b + 1) * (a + b + 2)) := by
    nlinarith [hprod, ha]
  -- nonzero facts needed to clear denominators
  have ha' : a ≠ 0 := ne_of_gt ha
  have hb' : b ≠ 0 := ne_of_gt hb
  have hs' : a + b ≠ 0 := ne_of_gt hs
  have hs1' : a + b + 1 ≠ 0 := by positivity
  have hs2' : a + b + 2 ≠ 0 := by positivity
  have ha1' : a + 1 ≠ 0 := by positivity
  -- express the ratio as a single closed-form fraction
  have hratio : f (a + 1) b / f a b
      = ((a + 1) * (a + b) ^ 2) / (a * ((a + b + 1) * (a + b + 2))) := by
    unfold f
    have e : a + 1 + b = a + b + 1 := by ring
    rw [e]
    field_simp
    ring
  rw [hratio, div_lt_one (by positivity)]
  exact key
The post Does additional data always reduce posterior variance? first appeared on John D. Cook.