Bayesian Statistics

In 1978 Ward Casscells, Arno Schoenberger and Thomas Grayboys stopped sixty house officers, students and attending physicians in the corridors of four Harvard teaching hospitals and asked each of them a single question. A disease affects one person in a thousand. A test for it produces a positive result in five percent of the people who do not have it. A patient tests positive, and you know nothing else about him. What is the chance he has the disease?

Twenty-seven of the sixty answered 95 percent. Eleven gave the answer that follows from the numbers they had just been handed, which is about 2 percent. The distance between 95 and 2 is not an arithmetic slip. It is a whole missing quantity, and this lesson is about that quantity and the rule that forces you to use it.

Counting people instead of manipulating fractions

Take a hundred thousand men and assume the test never misses a genuine case. One in a thousand has the disease, so 100 men do, and all 100 test positive. The other 99,900 are healthy, and five percent of them, which is 4,995 men, test positive anyway. That is 5,095 positive results in total, of which 100 belong to someone who is actually ill.

100 / 5,095 = 0.0196, or about 2 percent.

Nothing about that computation is subtle. The false positives outnumber the true positives fifty to one because there are a thousand times more healthy men than sick ones, and a five percent error rate applied to a very large group produces more positives than a perfect detection rate applied to a very small one. The physicians who said 95 percent were not confused about percentages. They had answered a different question: given that this man is ill, how often does the test say so? That is P(positive | disease). The question asked was P(disease | positive). The two are related, but they are not equal, and the thing that relates them is the base rate.

The rule, and where each piece comes from

Bayes' theorem is two lines of algebra. The probability that two things both happen can be written in either order:

P(H and D) = P(H | D) P(D) = P(D | H) P(H)

Divide by P(D) and you have it:

P(H | D) = P(D | H) P(H) / P(D)

Read the four pieces out loud, because their names are used constantly for the rest of this course. P(H) is the prior: what you believed about the hypothesis before the data arrived. P(D | H) is the likelihood: how probable this particular data would be if the hypothesis were true. P(D) is the marginal likelihood or evidence, the probability of seeing this data at all, averaged over every hypothesis you allow. And P(H | D) is the posterior, the revised belief. The denominator is almost always computed by splitting the world into cases:

P(D) = P(D | H) P(H) + P(D | not H) P(not H)

Thomas Bayes never published this. Richard Price found the essay among Bayes' papers after his death and read it to the Royal Society in December 1763, two years after Bayes died; Pierre-Simon Laplace rediscovered and generalised the result in the 1770s and applied it to birth records, court verdicts and the mass of Saturn. For most of the nineteenth century the rule was simply called inverse probability.

Key idea: Bayes' theorem is not a Bayesian doctrine. It is a theorem of probability that everyone agrees with. What is distinctively Bayesian, and what the rest of this course is about, is applying it to an unknown parameter rather than to an event.

Mammography: the example that has been run on physicians for forty years

David Eddy put a harder version to a hundred physicians in the early 1980s. His figures, rounded the way Gerd Gigerenzer and Ulrich Hoffrage later used them, are: one percent of women in the screened age band have breast cancer; the mammogram detects 80 percent of the cancers that are there; and it reports a positive result for 9.6 percent of the women who do not have cancer. A woman gets a positive mammogram. What is the probability she has cancer?

Ninety-five of the hundred physicians answered around 75 percent. Count ten thousand women instead.

Of the 1,030 women who get a positive result, 80 have cancer. So 80 / 1,030 = 0.0777, near enough 7.8 percent. Carrying the decimals rather than the rounded counts gives 0.01 times 0.80 = 0.00800 over 0.00800 + 0.99 times 0.096 = 0.10304, which is 0.07764. The physicians were high by a factor of ten.

Odds and likelihood ratios: the form you should actually use

Dividing by the marginal likelihood is a nuisance. Write the same rule for the hypothesis and its negation and divide one by the other, and P(D) cancels:

posterior odds = prior odds times likelihood ratio

The likelihood ratio of a positive mammogram is 0.80 / 0.096 = 8.333. The prior odds are 0.01 / 0.99 = 1/99 = 0.010101. So the posterior odds are 0.010101 times 8.333 = 0.08418, and converting back, 0.08418 / 1.08418 = 0.0776. Same answer, no denominator.

The payoff is that evidence now multiplies. Suppose the woman is recalled and a second, conditionally independent reading also comes back positive.

The last row uses the negative likelihood ratio, (1 - 0.80) / (1 - 0.096) = 0.20 / 0.904 = 0.2212. A negative reading is real evidence, but it divides the odds by only 4.5; it does not clear anyone. Note also that the conditional independence assumed here is generous. Two radiologists reading the same dense breast tissue make correlated errors, and if you treat correlated readings as independent the odds climb far too fast. That failure mode returns in Lesson 12, where correlated units are modelled properly.

If you take base-ten logarithms the multiplication becomes addition: each positive mammogram adds log(8.333) = 0.921 to the log odds, and each negative subtracts 0.655. Evidence, in this representation, is a quantity you can accumulate on a piece of paper, which is precisely how a sequential clinical workup feels from the inside.

Change the base rate and the same result means something else

Hold the test fixed at a likelihood ratio of 8.333 and vary only who is being screened.

One printed result, five different meanings, spanning a factor of fifty-eight. This is why the long argument over the age at which routine mammography should begin is not an argument about the machine. It is an argument about the base rate: breast cancer incidence in women in their early forties is several times lower than in women in their early sixties, so the identical positive report carries several times less weight, and the downstream biopsies fall on a population where far fewer of them find anything.

The point: a test result is not a verdict; it is a multiplier applied to whatever you believed beforehand. Refusing to state a prior does not remove the prior from the calculation, it only hides which one you used.

Common misconceptions

• "The test is 80 percent accurate, so a positive result is 80 percent likely to be right." Eighty percent is P(positive | cancer), a property of the instrument measured on known cases. The number a patient wants is P(cancer | positive), which additionally depends on how many people like her have cancer. Swapping the two is common enough to have a name: the prosecutor's fallacy.
• "The base rate is background information about other people, whereas the test is about this patient." The base rate is about this patient, in the only sense probability can be: it says which population she was drawn from before anyone looked at her. Discard it and you are not being more personal, you are silently substituting a prevalence of 50 percent.
• "A negative result rules the disease out." With a sensitivity of 80 percent the negative likelihood ratio is 0.22. Starting from a 10 percent prevalence, a negative mammogram leaves the probability at 2.4 percent, which is not zero and is roughly the prevalence at which screening was thought worth doing in the first place.
• "Bayes' theorem is what makes an analysis Bayesian." Every statistician uses this theorem; it follows from the axioms of probability. The Bayesian move is to put a probability distribution over an unknown quantity that is not random, such as a treatment effect or a physical constant, and that move is what Lessons 5 and 7 examine.

Where this leaves us

Bayes' theorem says posterior is proportional to likelihood times prior, and in odds form the messy denominator disappears: multiply the prior odds by the likelihood ratio. On the Casscells question that gives 2 percent rather than 95. On Eddy's mammogram it gives 7.8 percent rather than 75, rising to 41.2 percent on a second positive reading and falling to 1.8 percent if the second reading is negative. Change the base rate from 0.1 percent to 10 percent and the same positive report means anything from 0.83 percent to 48.1 percent.

So far every hypothesis has had exactly two states, disease or no disease. The next lesson replaces that with a continuous unknown, an unknown proportion, and asks what a dataset actually says about it before any prior is attached. The answer is the likelihood function, and it turns out to say strictly less, and strictly more, than a p-value does.

Sources

1. Wikipedia contributors. (n.d.). Bayes' theorem. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Base rate fallacy. Wikipedia. en.wikipedia.org
3. Casscells, W., Schoenberger, A., & Grayboys, T. B. (1978). Interpretation by physicians of clinical laboratory results. New England Journal of Medicine, 299(18), 999-1001.
4. Eddy, D. M. (1982). Probabilistic reasoning in clinical medicine: Problems and opportunities. In D. Kahneman, P. Slovic, & A. Tversky (Eds.), Judgment under Uncertainty: Heuristics and Biases (pp. 249-267). Cambridge University Press.
5. Gigerenzer, G., & Hoffrage, U. (1995). How to improve Bayesian reasoning without instruction: Frequency formats. Psychological Review, 102(4), 684-704.

Two laboratories send you the same line of data: nine heads in twelve tosses of the same coin. The first laboratory ran exactly twelve tosses, because twelve was the number written on the protocol before anyone touched the coin. The second kept tossing until the third tail appeared, which happened to be toss number twelve. Nothing about the coin differed. Nothing about the sequence of outcomes differed. You are asked, in both cases, whether the coin is fair.

The standard one-sided p-value for the first laboratory is 0.0730. For the second it is 0.0327. At the conventional threshold one result is not significant and the other is, and the only thing that changed was a sentence in a protocol describing an intention that was never acted on.

The function that holds what the data say

Start with the object that does not change between the two laboratories. The likelihood function is the probability of the observed data, written as a function of the unknown parameter with the data held fixed. For nine heads in twelve tosses under the fixed-sample design:

L(theta) = C(12,9) theta9 (1 - theta)3 = 220 theta9 (1 - theta)3

Read the variable carefully. In P(D | theta) as a probability, theta is fixed and D varies. In L(theta) the same expression is read the other way round: D is the nine heads you actually got, and theta ranges over every value it could have. That switch is the whole idea, and it is why a likelihood is not a probability distribution over theta. It does not integrate to one, and its absolute height means nothing. Only ratios of it carry information.

Here is the function, with the constant 220 dropped and every value divided by the largest so the peak reads 1.

The peak sits at theta = 9/12 = 0.75, the maximum likelihood estimate, and you can get it without calculus by noticing that the table is symmetric in a rough way around 0.75. With calculus: differentiate the log, 9/theta - 3/(1 - theta) = 0, giving theta = 9/12.

The ratio L(0.75)/L(0.50) = 0.0011732/0.0002441 = 4.81 says the data are 4.81 times more probable under a coin that lands heads three quarters of the time than under a fair one. That is a statement about evidence with no p-value, no prior and no threshold in it, and it is the same number for both laboratories.

If you keep every value of theta whose relative likelihood is at least 1/8, you get the interval from 0.461 to 0.936. Allan Birnbaum and later Richard Royall used exactly this construction as a way of reporting evidence: a 1/8 likelihood interval is roughly, though not exactly, comparable to a conventional 95 percent interval. Note that 0.50 sits inside it. The data are not strong.

The upshot: the likelihood function is the entire evidential content of the dataset about theta under the assumed model. Everything you do afterwards, Bayesian or not, is a way of summarising this curve or combining it with something else.

Curvature is precision

The width of the peak is not decoration. Take the log-likelihood l(theta) = 9 log theta + 3 log(1 - theta) and differentiate twice:

l''(theta) = -9/theta2 - 3/(1 - theta)2, and at theta = 0.75 this is -16 - 48 = -64.

The observed information is I = 64, and the classical standard error is 1/sqrt(64) = 0.125, which is exactly the familiar sqrt(p(1-p)/n) = sqrt(0.75 times 0.25/12) = 0.125. A sharply curved log-likelihood is a precise estimate; a flat one is a vague one. Keep that picture: in Lesson 4 the same curvature reappears as a precision, and precisions are what get added when you combine a prior with data.

Two designs, two p-values, one likelihood

Now the arithmetic that started the lesson. Under the first laboratory's design the number of tosses is fixed at 12 and the number of heads is random. The one-sided p-value asks how often a fair coin gives nine or more heads in twelve tosses:

P(X >= 9) = [C(12,9) + C(12,10) + C(12,11) + C(12,12)]/212 = (220 + 66 + 12 + 1)/4096 = 299/4096 = 0.0730

Under the second laboratory's design the number of tails is fixed at 3 and the number of tosses is random. The data are extreme if it took a long time to accumulate three tails, so the p-value asks how often a fair coin needs twelve or more tosses to produce three tails. That happens exactly when the first eleven tosses contain at most two tails:

P(N >= 12) = [C(11,0) + C(11,1) + C(11,2)]/211 = (1 + 11 + 55)/2048 = 67/2048 = 0.0327

Both are correct calculations of what they claim to calculate. Now write down the two likelihood functions. Fixed sample: 220 theta9(1 - theta)3. Stop at the third tail, using the negative binomial mass function: C(11,2) theta9(1 - theta)3 = 55 theta9(1 - theta)3.

The two differ by a factor of four that does not depend on theta at all. As functions of theta they have the same shape, the same peak, the same curvature, the same ratios everywhere. A ratio of likelihoods, which is the only thing a likelihood is good for, is identical under the two designs.

The likelihood principle, and what it costs each side

The likelihood principle says that if two experiments produce likelihood functions proportional to one another, they carry the same evidence about the parameter. Allan Birnbaum proved in 1962 that it follows from two assumptions most statisticians already accept: sufficiency, which says a sufficient statistic carries everything relevant, and conditionality, which says that if you flip a coin to decide which of two experiments to run, you should analyse the one you actually ran.

A Bayesian analysis obeys the principle without trying. The posterior is proportional to likelihood times prior, and a constant that does not involve theta is absorbed into the normalising denominator. Put a uniform prior on theta and both laboratories get the same Beta(10, 4) posterior, with mean 10/14 = 0.714 and P(theta > 0.5) = 0.954.

A p-value cannot obey it. Look again at what the two calculations sum over. The fixed-sample p-value adds up the probabilities of ten, eleven and twelve heads, outcomes that did not happen. The sequential p-value adds up the probabilities of runs of thirteen, fourteen, fifteen tosses, which also did not happen. Different designs make different things possible, so the set of unobserved outcomes differs, so the tail area differs. This is not a defect in anyone's arithmetic. It is a consequence of measuring evidence by comparing what you saw against what you might have seen.

Remember: the disagreement between 0.0730 and 0.0327 is not about the coin. It is about whether an inference may depend on data that were never collected. The likelihood principle says no; frequentist error control says it must, because error rates are defined over repetitions that include those outcomes.

Where this bites in practice

Stopping rules are not a philosopher's example. A clinical trial with interim analyses is a stopping rule, and frequentist practice compensates with alpha spending functions, O'Brien-Fleming boundaries and adjusted final p-values, all of which exist because the sampling distribution depends on when you looked. A Bayesian analysis of the same trial uses the same posterior whether the data arrived in one batch or in six, which is convenient and is also the thing critics point at: if peeking costs nothing, what stops an investigator from watching the posterior and stopping the moment it looks good?

The honest answer is that optional stopping does not bias a Bayesian posterior, but it does change what a reported selected result means. If you stop the first time the posterior probability exceeds 0.95, you will publish more than five percent of the time when nothing is happening. The posterior you report is still correctly computed from the data in hand; what has gone wrong is the selection of which posteriors get reported. Keeping the whole sequence, or modelling the selection, fixes it. Pretending the stopping rule is irrelevant to publication decisions does not.

Common misconceptions

• "The likelihood function is the probability that theta takes each value." It is not a distribution over theta. The area under the curve in the table above is about 0.00009, not 1, and doubling every entry would change nothing. To turn a likelihood into a probability statement about theta you must supply a prior, which is exactly the step frequentists decline to take.
• "The two p-values differ because the second design gives more information." Both designs delivered the identical string of twelve outcomes. The difference sits entirely in which unobserved datasets are counted as more extreme.
• "A maximum likelihood estimate is the most probable value of the parameter." It is the value that makes the observed data most probable, which is a different sentence. With a uniform prior the two happen to coincide for the mode; with any other prior they do not.
• "Since Bayesian analysis ignores the stopping rule, you can peek as often as you like with no consequences." The posterior is unaffected, but the decision to publish when it crosses a threshold is a selection effect that inflates the rate of confident wrong claims. Ignoring the stopping rule is not the same as ignoring the selection rule.

What to carry forward

The likelihood function L(theta) = theta9(1 - theta)3 peaks at 0.75, and its curvature there, l'' = -64, gives a standard error of 0.125. Ratios of the curve are the evidence: 4.81 to 1 for a three-quarters coin against a fair one, and a 1/8 likelihood interval from 0.461 to 0.936 that comfortably contains 0.5. The same twelve tosses give a p-value of 0.0730 under a fixed-sample design and 0.0327 under a stop-at-the-third-tail design, because the two designs make different unobserved outcomes possible, while the two likelihoods differ only by the constants 220 and 55. Birnbaum's likelihood principle says the evidence is therefore the same, Bayesian updating obeys it automatically, and tail-area testing cannot.

The likelihood is half of Bayes' theorem. The next lesson supplies the other half and multiplies them together on twenty-three space shuttle flights, where the arithmetic turns out to be a matter of adding two numbers to two other numbers.

Sources

1. Wikipedia contributors. (n.d.). Likelihood function. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Likelihood principle. Wikipedia. en.wikipedia.org
3. Birnbaum, A. (1962). On the foundations of statistical inference. Journal of the American Statistical Association, 57(298), 269-306.
4. Berger, J. O., & Wolpert, R. L. (1988). The Likelihood Principle (2nd ed.). Institute of Mathematical Statistics.
5. Royall, R. (1997). Statistical Evidence: A Likelihood Paradigm. Chapman and Hall.

Seven of the twenty-three space shuttle flights before Challenger came back with damaged O-rings. The boosters were recovered from the Atlantic, floated back to Cape Canaveral, taken apart and inspected joint by joint, and on seven of those flights at least one of the six field-joint O-rings showed erosion or blow-by. Here are the seven.

The other sixteen flights came back clean. The question this lesson answers is narrow and answerable: given 7 out of 23, what should you believe about the per-flight probability of O-ring distress, and how sure should you be? By the end you will have a full probability distribution over that rate, obtained by adding two numbers to two other numbers.

A distribution over a proportion

Call the unknown per-flight probability of distress theta. The data are 23 flights with 7 distressed, and if the flights are treated as independent trials with a common theta the likelihood is binomial:

L(theta) = C(23,7) theta7 (1 - theta)16

Lesson 2 stopped here. The maximum is at 7/23 = 0.3043, and the curvature gives a standard error of sqrt(0.3043 times 0.6957/23) = 0.0960. To get a probability distribution over theta rather than a curve with no scale, you need a prior, and the prior has to be a distribution on the interval from 0 to 1.

The beta distribution is the natural candidate. Its density is

p(theta) proportional to thetaa-1 (1 - theta)b-1, for 0 < theta < 1, with mean a/(a + b) and, for a, b > 1, mode (a - 1)/(a + b - 2).

Beta(1, 1) is flat. Beta(2, 18) is a hump with mean 0.1 and most of its mass below 0.25. Beta(0.5, 0.5), the Jeffreys prior for a proportion, is U-shaped, piling mass near 0 and near 1. The parameters have a reading that makes the rest of the lesson easy: a Beta(a, b) prior behaves like a previous dataset containing a - 1 successes and b - 1 failures.

Why the algebra collapses

Multiply the prior by the likelihood and keep only what depends on theta:

p(theta | y) proportional to thetaa-1(1 - theta)b-1 times thetay(1 - theta)n-y = thetaa+y-1(1 - theta)b+n-y-1

That is a beta density again, with parameters a + y and b + n - y. No integral was evaluated. The marginal likelihood in the denominator of Bayes' theorem is whatever number makes this integrate to one, and because the answer is a known family you can look that constant up rather than compute it. A prior family that reproduces itself under the likelihood in this way is called conjugate.

What matters here: conjugate updating is bookkeeping. Add the successes to the first parameter, add the failures to the second, and you are done. Every conjugate model in the next two lessons has the same shape.

Running it on the shuttle

Start with Beta(1, 1), the flat prior, which says nothing more than that theta lies between 0 and 1. Add 7 successes and 16 failures:

Beta(1 + 7, 1 + 16) = Beta(8, 17)

Now read the posterior off. Mean 8/25 = 0.3200. Mode 7/25 - careful, the mode is (8 - 1)/(25 - 2) = 7/23 = 0.3043, which is exactly the maximum likelihood estimate, as it must be when the prior is flat. Standard deviation sqrt(ab/((a+b)2(a+b+1))) = sqrt(136/(625 times 26)) = 0.0915. The 2.5th and 97.5th percentiles are 0.1563 and 0.5109, so the 95 percent central credible interval is 0.156 to 0.511.

That interval is enormous, and it should be. Twenty-three flights is not much data about a rate near a third; the interval spans a factor of three and a quarter. Two further quantities come from the same distribution with no extra modelling. The probability that the true rate exceeds one half is 0.0320. The probability that it is below 0.10, which is roughly what the flight-readiness reviews were treating as acceptable, is 0.0017.

Change the prior and see what moves.

The first two barely differ; they disagree about a half-observation. The third moves the posterior mean down by a third of its value, because Beta(2, 18) is worth 20 prior observations and 20 is comparable to 23. That is the honest statement about prior sensitivity: a prior matters exactly to the extent that its equivalent sample size is comparable to the real one.

The posterior mean is an average of what you thought and what you saw

Write the posterior mean out and split it:

(a + y)/(a + b + n) = [(a + b)/(a + b + n)] times [a/(a + b)] + [n/(a + b + n)] times [y/n]

The posterior mean is a weighted average of the prior mean and the sample proportion, with weights proportional to the prior sample size a + b and the actual sample size n. Check it on the flat prior: (2/25) times 0.5 + (23/25) times 0.3043 = 0.0400 + 0.2800 = 0.3200. Check it on the engineering prior: (20/43) times 0.100 + (23/43) times 0.3043 = 0.0465 + 0.1628 = 0.2093. In that second case the prior carries 47 percent of the weight, which is what it means to bring 20 pseudo-observations to a 23-observation problem.

This weighted-average reading is the single most portable idea in the course. Lesson 4 shows the same structure for counts and for normal means, where the weights become precisions.

The order of the flights does not matter, and one flight at a time gives the same answer

Take the first ten flights in the table, of which three had distress, and update: Beta(1, 1) becomes Beta(4, 8), with mean 0.333 and standard deviation 0.131. Now feed in the remaining thirteen flights, four of which had distress: Beta(4 + 4, 8 + 9) = Beta(8, 17). Identical to doing it in one batch. Yesterday's posterior is today's prior, and because addition is commutative the order of arrival is irrelevant. This is what people mean when they say Bayesian inference is sequential by construction, and it is why an interim analysis needs no correction, as Lesson 2 argued.

Predicting the next flight

The question a flight-readiness review actually asks is not about theta. It is: will the next one come back damaged? Average the Bernoulli probability over the posterior:

P(next flight distressed) = integral theta p(theta | y) d theta = E[theta | y] = 8/25 = 0.32

The posterior predictive probability for a single new trial is just the posterior mean. This is also Laplace's rule of succession, (y + 1)/(n + 2), published in 1774 and used by him to compute the probability that the sun would rise tomorrow. For two future flights the answer is not 0.322: the flights are conditionally independent given theta but not marginally independent, because a damaged flight raises your estimate of theta. The correct figure, from the beta-binomial predictive, is (8 times 9)/(25 times 26) = 0.1108, against 0.1024 for the naive product.

In short: a posterior is not the answer to a decision problem; it is the input. The predictive distribution is what you push forward into tomorrow, and it is wider than plugging in a point estimate, because it carries the uncertainty in theta with it.

What this model gets wrong

Look back at the temperature column. Every flight is being treated as an exchangeable draw with the same theta, and the coldest launch in the table, STS 51-C at 53 degrees, had two distressed O-rings. On 28 January 1986 the air temperature at Pad 39B was 36 degrees Fahrenheit at launch, and the Rogers Commission put the temperature of the right-hand booster's field joint near 28. The single-theta model has nothing to say about that, because it threw the temperature away.

It is worse than that. The chart the engineers took into the teleconference on the night of 27 January showed only the seven flights with damage, because those were the interesting ones. Plot only those seven and temperature looks unrelated to damage. Plot all twenty-three, including the sixteen clean flights clustered at warm temperatures, and the relationship is visible. Conditioning on the outcome destroyed the evidence. Lesson 14 puts temperature back in as a predictor and Lesson 16 asks whether the model with temperature is better supported than the model without it.

Common misconceptions

• "A flat prior means no prior." Beta(1, 1) is worth two pseudo-observations, one success and one failure, and it drags the posterior mean from 0.3043 up to 0.3200. That is small but not nothing, and Lesson 5 shows a case with zero successes where the choice among flat priors changes the answer by a factor of two.
• "The posterior mean is the estimate and the rest is detail." With a 95 percent interval running from 0.156 to 0.511, reporting 0.32 alone hides everything the analysis learned. The whole point of carrying a distribution is that its width is part of the answer.
• "Conjugacy is what makes the analysis Bayesian." Conjugacy is a computational convenience that saves an integral. Lessons 9 to 11 do the same inference on models with no conjugate form at all, using simulation, and nothing philosophical changes.
• "Seven out of twenty-three is 30 percent, so the shuttle had a 30 percent chance of failing." Theta here is the probability that an O-ring shows erosion or blow-by on a recovered booster, not the probability of losing the vehicle. Most distressed flights returned safely. Naming the event precisely is half of applied probability.

Putting it together

A Beta(a, b) prior and a binomial likelihood give a Beta(a + y, b + n - y) posterior, so updating is addition. On 7 distressed flights out of 23 with a flat prior that is Beta(8, 17): posterior mean 0.320, mode 0.3043, standard deviation 0.0915, 95 percent credible interval 0.156 to 0.511, P(theta > 0.5) = 0.032 and P(theta < 0.10) = 0.0017. The posterior mean is the weighted average (2/25)(0.5) + (23/25)(0.3043), and swapping in a Beta(2, 18) prior, worth 20 pseudo-observations, pulls it down to 0.209 because 20 is comparable to 23. Updating ten flights and then thirteen gives the same Beta(8, 17) as updating all twenty-three at once. The predictive probability for the next flight is the posterior mean, 0.32, and for two flights it is 0.1108 rather than 0.32 squared.

Counts and continuous measurements have their own conjugate pairs, and the next lesson works both: Prussian cavalrymen killed by horses, and Michelson's hundred attempts to measure the speed of light.

Sources

1. Wikipedia contributors. (n.d.). Conjugate prior. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Space Shuttle Challenger disaster. Wikipedia. en.wikipedia.org
3. Dalal, S. R., Fowlkes, E. B., & Hoadley, B. (1989). Risk analysis of the space shuttle: Pre-Challenger prediction of failure. Journal of the American Statistical Association, 84(408), 945-957.
4. Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.), Chapter 2. Chapman and Hall/CRC.

Ladislaus Bortkiewicz counted 122 dead soldiers. Working through the casualty registers of the Prussian army, he pulled out every death caused by the kick of a horse in ten cavalry corps over the twenty years from 1875 to 1894, and tabulated them one corps-year at a time. That is 200 corps-years and 122 deaths, and the table he published in 1898 in Das Gesetz der kleinen Zahlen is still the standard first example of a rare-event count.

This lesson runs two conjugate updates side by side, one on those counts and one on a hundred measurements of the speed of light, and shows that both come out as the same arithmetic: a weighted average in which the weights are precisions.

Counts: gamma is conjugate to Poisson

Model each corps-year as Poisson with rate lambda. For n independent corps-years with total count S = sum yi the likelihood, dropping factorials that do not involve lambda, is

L(lambda) proportional to lambdaS e-n lambda

A gamma density in the shape-rate parameterisation is p(lambda) proportional to lambdaalpha-1 e-beta lambda, with mean alpha/beta and variance alpha/beta2. Multiply:

p(lambda | y) proportional to lambdaalpha + S - 1 e-(beta + n) lambda = Gamma(alpha + S, beta + n)

Add the total count to the shape and the number of exposure units to the rate. The gamma parameters read as pseudo-data in the same way the beta parameters did: alpha is a prior count of events and beta is a prior number of exposure units.

Take Gamma(2, 2) as the prior. It says: I have seen something like 2 deaths in 2 corps-years, so my prior mean is 1.0 death per corps-year, with prior standard deviation sqrt(2)/2 = 0.707. That is a weak and slightly pessimistic belief. Update it:

Gamma(2 + 122, 2 + 200) = Gamma(124, 202)

Posterior mean 124/202 = 0.6139. Posterior standard deviation sqrt(124)/202 = 0.0551. The 2.5th and 97.5th percentiles are 0.5106 and 0.7265, so the rate is somewhere between about half a death and three quarters of a death per corps-year, and the interval is narrow because 200 corps-years is a lot of exposure.

Now split the mean:

(alpha + S)/(beta + n) = [beta/(beta + n)] times [alpha/beta] + [n/(beta + n)] times [S/n]

= (2/202) times 1.000 + (200/202) times 0.610 = 0.0099 + 0.6040 = 0.6139

The prior carries one percent of the weight, because it is worth 2 corps-years against 200 real ones. Compare two other priors: the Jeffreys prior for a Poisson rate is proportional to lambda-1/2, a Gamma(0.5, 0) limit, giving posterior Gamma(122.5, 200) with mean 0.6125; the flat prior Gamma(1, 0) gives Gamma(123, 200) with mean 0.6150. All three answers agree to two decimal places, which is what a prior looks like when the data outnumber it a hundred to one.

Why this matters: arguments about priors are arguments about sample size in disguise. Ask how many observations the prior is worth and compare it with how many you have; if the ratio is 1 to 100 the argument is not worth having, and if it is 20 to 23, as in Lesson 3, it is the whole analysis.

Prediction for a rare count

What is the probability that a given corps loses nobody to a horse next year? Averaging the Poisson mass at zero over the gamma posterior gives a negative binomial predictive, and for zero events in one unit of exposure it comes out as

P(0) = [beta*/(beta* + 1)]alpha* = (202/203)124 = 0.5421

against 0.5412 if you simply plug in the posterior mean, e-0.6139. The difference is 0.0009, which is negligible here and will not be negligible in Lesson 12, where each group has only a handful of observations and the plug-in answer is badly overconfident.

Measurements: normal is conjugate to normal

In the summer of 1879, at the United States Naval Academy at Annapolis, Albert Michelson ran a rotating-mirror apparatus 100 times and recorded 100 determinations of the speed of light. The standard coding of that dataset reports each run in kilometres per second with 299,000 subtracted. The mean of the hundred runs is 852.4 and the sample standard deviation is 79.01, so the standard error of the mean is 7.90.

Treat sigma as known at 79.0 for the moment; Lesson 6 comes back to what changes when it is not. With n observations of mean ybar from N(mu, sigma2) and a prior mu ~ N(mu0, tau02), define precisions as reciprocal variances:

prior precision = 1/tau02, data precision = n/sigma2.

Then the posterior is normal with

posterior precision = 1/tau02 + n/sigma2, and posterior mean = [ (1/tau02) mu0 + (n/sigma2) ybar ] / posterior precision.

Precisions add. That single sentence is the whole normal-normal model, and it is worth more than the formula because it tells you what to expect before you compute: information accumulates, never decreases, and the posterior is always at least as sharp as the prior.

Suppose a physicist in 1879 comes to the data believing the speed of light is near 299,800 km/s with a standard deviation of 20, that is mu0 = 800 and tau0 = 20 in the coded units. Then:

Check the weighted average directly: 0.135 times 800 + 0.865 times 852.4 = 108.0 + 737.3 = 845.3. The posterior standard deviation, 7.348, is smaller than the data-only standard error of 7.90 and smaller than the prior standard deviation of 20, as it has to be once you accept that precisions add.

Where this analysis goes wrong, and why that is instructive

The speed of light in vacuum has been fixed by definition since 1983 at 299,792.458 km/s, which is 792.458 in these coded units. Michelson's mean of 852.4 is 59.9 too high, and the standard error of his mean is 7.90, so he is 7.6 standard errors off. The posterior above puts a 95 percent interval at 830.9 to 859.7, and the truth is not in it.

Nothing in the arithmetic failed. The model assumed the 100 runs were independent draws around the true value, and they were not: they shared a systematic error, from the measured distance between the mirrors and from the reduction to vacuum, which no amount of repetition removes. Averaging 100 runs shrinks random error by a factor of 10 and shrinks systematic error by a factor of 1. Stephen Stigler used this dataset in 1977 precisely to make that point, and it is the reason experimental physicists report two uncertainties.

Bottom line: a posterior is a statement conditional on the model, and a confident wrong posterior almost always means the model, not the arithmetic. In Lesson 15 you will learn how to catch this kind of failure from the data itself, using posterior predictive checks.

The three conjugate pairs side by side

The pattern is not a coincidence. All three likelihoods belong to the exponential family, where the log-likelihood is linear in a sufficient statistic, and for any such likelihood a conjugate prior exists and updating amounts to adding the observed sufficient statistics to prior pseudo-counts. That is the whole supply of models you can do in closed form, which is why Module 4 exists.

Common misconceptions

• "Variances add, so the posterior variance is the prior variance plus the data variance." Precisions add, which means variances combine harmonically. With a prior variance of 400 and a data variance of 62.4, the posterior variance is 54.0, smaller than either. Adding variances would make the posterior less certain than the prior, which is impossible after seeing data.
• "A narrow posterior means an accurate one." Michelson's posterior standard deviation is 7.35 and his error is 59.9. Width measures how much the model thinks it has learned, not whether the model is right.
• "The gamma prior on a rate has to be proper, so Gamma(1, 0) is not allowed." Gamma(1, 0) is improper, since it integrates to infinity, but it produces a proper posterior as soon as any exposure is observed. Improper priors are usable when they yield proper posteriors, and dangerous when nobody checks. Lesson 5 shows a case where nobody checked.
• "Bortkiewicz's data show that horse kicks are random." They show that the count in a corps-year is well described by a Poisson distribution with a rate near 0.61, which is an empirical claim about a table, not a metaphysical one about horses. The fit is close, but the value of 4 in a single corps-year is a mild strain on a model that expects 0.6 such years in 200.

The short version

Gamma is conjugate to Poisson: add the event total to the shape and the exposure to the rate. Bortkiewicz's 122 deaths in 200 corps-years update a Gamma(2, 2) prior to Gamma(124, 202), mean 0.6139, standard deviation 0.0551, 95 percent interval 0.511 to 0.727, and that mean is (2/202)(1.000) + (200/202)(0.610). Normal is conjugate to normal with known variance, and precisions add: Michelson's 100 runs with mean 852.4 and standard error 7.90, against a prior of 800 with standard deviation 20, give a posterior mean of 845.33 with standard deviation 7.348, the data carrying 86.5 percent of the weight. That posterior excludes the true value of 792.458 by a wide margin, because 100 repetitions do nothing to a systematic error. All three conjugate models in the course share one sentence: the posterior mean is a weighted average of the prior mean and the data mean, with weights given by how much each is worth.

Every calculation so far has taken the prior as given. The next lesson asks where it comes from, and shows three priors that are all described as uninformative giving three different answers to the same question.

Sources

1. Wikipedia contributors. (n.d.). Poisson distribution. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Ladislaus Bortkiewicz. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Albert A. Michelson. Wikipedia. en.wikipedia.org
4. Bortkiewicz, L. von. (1898). Das Gesetz der kleinen Zahlen. B. G. Teubner.
5. Stigler, S. M. (1977). Do robust estimators work with real data? Annals of Statistics, 5(6), 1055-1098.

Ten patients received the device. None of them had a serious adverse event. The maximum likelihood estimate of the event rate is therefore exactly zero, which is the one number everybody in the room knows is wrong, because nobody in the room would accept a device that is guaranteed never to harm anyone on the strength of ten patients.

So you go Bayesian, and you decide to be neutral about it, because the trial is small and you do not want to be accused of putting your thumb on the scale. You reach for an uninformative prior. There are at least three standard candidates, and this lesson works out what each of them says.

Three answers, and one non-answer

The likelihood for zero successes in ten Bernoulli trials is (1 - theta)10, a curve that starts at 1 and slides monotonically to 0. It has no interior peak, so the maximum likelihood estimate sits on the boundary at 0.

Two answers differ by a factor of 1.83, and the third does not exist. The Haldane prior, proportional to 1/(theta(1 - theta)), is improper, and with zero successes it stays improper after the data: the integral of theta-1(1 - theta)9 diverges at the origin. There is no posterior mean to report because there is no posterior. That failure is not a curiosity. It is what happens whenever an improper prior meets data that fail to pin down the parameter from one side, and nothing in the software will tell you.

For comparison, a frequentist analysing the same ten patients has an exact one-sided 95 percent upper confidence bound of 1 - 0.051/10 = 0.2589, and the familiar rule of three gives 3/10 = 0.30. The uniform-prior answer of 0.2384 is close to both; the Jeffreys answer of 0.1708 is noticeably more optimistic.

The core of it: there is no such thing as no prior. There are only priors whose content you have not inspected, and with ten observations and zero events, the content is most of the answer.

Why nothing can be uninformative about everything at once

Here is the argument that ends the debate, and it takes three lines. Suppose you put a uniform prior on theta, the event rate, because you claim to know nothing. Then consider psi = theta2, which is an equally legitimate parameter, and about which you also claim to know nothing. Change variables:

p(psi) = p(theta) |d theta/d psi| = 1 times 1/(2 sqrt(psi))

That density is infinite at 0 and falls away; under it, P(psi < 0.25) = 0.5, whereas a uniform prior on psi would put that probability at 0.25. So being ignorant about theta means having strong opinions about theta squared, and being ignorant about theta squared means having opinions about theta. Ignorance cannot be flat in every coordinate system, because flatness is not a property of a belief, it is a property of a belief plus a coordinate system.

The same thing happens on the scale statisticians actually use for proportions. Put a uniform prior on the log-odds x = log(theta/(1 - theta)), and the implied density on theta is 1/(theta(1 - theta)), the Haldane prior, which piles infinite mass at both ends. Go the other way, uniform on theta, and the implied density on the log-odds is the standard logistic density, peaked at x = 0, that is, at theta = 0.5. Someone who calls the uniform prior neutral is asserting on the log-odds scale that the truth is probably near a coin flip.

Jeffreys' repair

Harold Jeffreys proposed a rule with the invariance the flat prior lacks: take the prior proportional to the square root of the Fisher information.

p(theta) proportional to sqrt(I(theta)), where I(theta) = -E[d2 log L/d theta2].

Because the Fisher information transforms with the square of the Jacobian, the square root transforms with the Jacobian itself, which is exactly what a density must do. Derive the Jeffreys prior once in a coordinate system and you get the same beliefs in every other coordinate system, automatically.

The third row explains something you have already met: with a flat prior on mu the normal-normal posterior collapses to N(ybar, sigma2/n), so the 95 percent credible interval is numerically identical to the classical confidence interval. Many of the agreements between the two traditions are of this kind, and Lesson 6 shows where the agreement stops.

Jeffreys' rule is not a solution to everything. In more than one dimension it can misbehave: for a normal with both mu and sigma unknown, the multivariate Jeffreys rule gives 1/sigma2, while Jeffreys himself recommended 1/sigma, obtained by treating the two parameters separately. When a rule's own author overrides it, treat it as a useful default rather than a law.

A vocabulary that is actually usable

Drop the word uninformative. Four categories do more work.

• Informative. Built from a named source: a meta-analysis, a previous trial, a physical constraint. Beta(2, 18) in Lesson 3 was informative and worth 20 flights, and it should have been defended as such.
• Weakly informative. Deliberately vaguer than your real belief, but tight enough to rule out the absurd. It contributes almost nothing when the data are plentiful and stabilises the answer when they are not.
• Reference or objective. Chosen by a formal rule such as Jeffreys' so that the answer is reproducible by anyone with the same model. Reproducible is not the same as neutral.
• Improper. Not a probability distribution at all. Sometimes harmless, sometimes fatal, as the Haldane case above shows and as Lesson 12 shows again for variance components.

To see what a prior really claims, simulate from it and look at the data it produces. This is a prior predictive check, and it is quick. Take a logistic regression coefficient on a standardised predictor and give it a N(0, 10²) prior, the kind of thing routinely described as vague. A 95 percent prior interval for the coefficient runs from -19.6 to +19.6, so the prior asserts that a one standard deviation change in the predictor plausibly multiplies the odds by anything between e-19.6 = 3.1 x 10-9 and e19.6 = 3.3 x 108. No effect in medicine, economics or ecology is that large. Replace it with N(0, 1) and the same interval runs from an odds ratio of 0.14 to 7.1, which covers every real effect anyone has ever measured while excluding the impossible.

So what?: a prior that puts most of its mass on impossible parameter values is not neutral, it is wrong, and it will bite hardest exactly when the data are thin and you needed help most.

Back to the ten patients

With the vocabulary in hand, the trial is easy. You are not ignorant about device-related adverse events; comparable devices run at 1 to 5 percent. A Beta(1, 30) prior encodes a prior mean of 0.032 and is worth 31 pseudo-patients, which is too strong for a first-in-human study. Beta(0.5, 10) is worth about 10 and gives a posterior of Beta(0.5, 20), mean 0.0244, 95th percentile 0.0927. Beta(1, 4), worth 5 pseudo-patients with prior mean 0.2, gives Beta(1, 14), mean 0.0667, 95th percentile 0.1926. Report the sensitivity, name the source of each prior, and let the reader see that the honest conclusion from ten patients is that the rate is probably under 20 percent and could easily be 5.

Common misconceptions

• "A flat prior lets the data speak for themselves." Flat on which scale? Flat on theta is peaked at 0.5 on the log-odds scale, and flat on the log-odds is infinitely peaked at 0 and 1 on the theta scale. There is no scale-free flatness.
• "An improper prior is fine as long as the software runs." The software will happily return samples that look like a posterior when no posterior exists. The Haldane case with zero successes is the simplest example; a flat prior on a hierarchical standard deviation is the one that actually appears in published work.
• "The Jeffreys prior is the objectively correct prior." It is the prior invariant under reparameterisation, which is a real and valuable property and not the same as correctness. It also produces answers Jeffreys himself declined to use in more than one dimension.
• "With enough data the prior does not matter, so the choice is academic." True asymptotically, by the Bernstein-von Mises theorem, and irrelevant to a trial with ten patients, a hierarchical model with four groups, or any model with more parameters than observations. The prior stops mattering exactly when you no longer needed it.

What to remember

Zero events in ten patients gives a posterior mean of 0.0833 under a uniform prior, 0.0455 under Jeffreys, and nothing at all under Haldane, whose posterior is improper. A uniform prior on theta implies a density of 1/(2 sqrt(psi)) on psi = theta2 and the logistic density on the log-odds, so flatness is a property of a coordinate system rather than of a state of belief. The Jeffreys prior, proportional to the square root of the Fisher information, is invariant under reparameterisation and gives Beta(0.5, 0.5) for a proportion, 1/sigma for a scale and a flat prior for a normal mean. Replace the word uninformative with informative, weakly informative, reference or improper, and check any prior by simulating from it: a N(0, 10²) prior on a log-odds coefficient asserts odds ratios from 3 x 10-9 to 3 x 108, which is not vagueness but nonsense.

You now have a posterior and a defensible prior. The next question is what the interval around that posterior means, and whether the confidence interval printed next to it in the same table means the same thing. It does not, and there is a two-observation example where the difference is total.

Sources

1. Wikipedia contributors. (n.d.). Prior probability. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Jeffreys prior. Wikipedia. en.wikipedia.org
3. Gelman, A., Jakulin, A., Pittau, M. G., & Su, Y.-S. (2008). A weakly informative default prior distribution for logistic and other regression models. Annals of Applied Statistics, 2(4), 1360-1383. arxiv.org
4. Jeffreys, H. (1961). Theory of Probability (3rd ed.). Oxford University Press.
5. Kass, R. E., & Wasserman, L. (1996). The selection of prior distributions by formal rules. Journal of the American Statistical Association, 91(435), 1343-1370.

Two numbers come off a gauge: 3.1 and 3.9. The gauge is honest but crude, and its reading error is uniform on plus or minus one half of a unit around the true value theta. You take the interval running from the smaller reading to the larger, from 3.1 to 3.9, and you look up its confidence level. It is exactly 50 percent.

Now think about what theta can be. The reading 3.9 came from somewhere within half a unit of theta, so theta > 3.4. The reading 3.1 came from within half a unit of theta as well, so theta < 3.6. The true value is in the interval from 3.4 to 3.6, and that whole interval sits inside 3.1 to 3.9. Your 50 percent confidence interval contains theta with probability one, and you know it, and the confidence level is still 50 percent.

Two definitions that do different jobs

A confidence interval is a promise about a procedure. A credible interval is a statement about a parameter. The reason nobody notices the difference most of the time is that in the models people teach first, the two come out numerically identical.

Where they agree, to three decimal places

Take the shuttle posterior from Lesson 3, Beta(8, 17) from 7 distressed flights in 23.

The Wilson interval and the flat-prior credible interval agree to within four thousandths at both ends, which is not a coincidence: the Wilson construction is the score interval, and the score statistic is the derivative of the log-likelihood, the same curve the posterior is built from. The same thing happened in Lesson 5: with a flat prior on a normal mean, the credible interval is exactly ybar +/- 1.96 sigma/sqrt(n), endpoint for endpoint.

Numerical agreement is not conceptual agreement. Both intervals for the shuttle run from about 0.16 to about 0.51. The Bayesian says the probability that theta is in there is 0.95. The frequentist says that if you repeated the twenty-three flights many times, 95 percent of the intervals you would construct would contain the true theta, and declines to say anything about this one. Most of the time nobody is harmed by conflating them.

Welch's example, worked

Now the case where it matters. Let X1 and X2 be independent and uniform on the interval from theta - 1/2 to theta + 1/2. Write L for the smaller observation, U for the larger, and d = U - L for the range.

The interval (L, U) contains theta exactly when one observation falls below theta and the other above, and since each falls above with probability 1/2 independently, that happens with probability 2 times (1/2)(1/2) = 1/2. It is a valid 50 percent confidence interval, with coverage exactly 0.5 for every theta.

But condition on the range. Given d, theta must satisfy U - 1/2 < theta < L + 1/2, an interval of length 1 - d, and within it theta is uniformly distributed under a flat prior. The interval (L, U) covers the part of that range where theta also lies between L and U, giving a conditional coverage of d/(1 - d) when d < 1/2 and 1 when d >= 1/2.

Average d/(1 - d), truncated at 1, over the distribution of the range, which has density 2(1 - d), and you get 0.25 + 0.25 = 0.50. The procedure has exactly the coverage it advertises. It achieves that average by being certain when the observations are far apart and nearly useless when they are close together, and the data tell you which case you are in.

The Bayesian analysis with a flat prior on theta is trivial and never makes this mistake. The likelihood is 1 when theta is within half a unit of both observations and 0 otherwise, so the posterior is uniform on (U - 1/2, L + 1/2). For 3.1 and 3.9 that is uniform on (3.4, 3.6), and a 50 percent credible interval is any subinterval of length 0.1, say (3.45, 3.55). For 3.10 and 3.15 the posterior is uniform on (2.65, 3.60), and a 50 percent credible interval has length 0.475. The credible interval widens when the data are uninformative, which is what an interval is supposed to do.

Worth holding on to: coverage is a property averaged over datasets you did not get. When the data contain a statistic that tells you how informative they are, the average can be a poor description of the case in hand.

The general form of the problem

The range d in Welch's example is an ancillary statistic: its distribution does not depend on theta, so it carries no information about theta by itself, and yet it determines how much information the rest of the data carry. David Cox made the point in 1958 with a cleaner story. You measure a quantity with one of two instruments, chosen by a fair coin: instrument A has standard deviation 1, instrument B has standard deviation 100. The coin comes up heads and you use A.

The unconditional 95 percent requirement can be met by a procedure with 99 percent coverage on A-days and 91 percent on B-days, since 0.5 times 0.99 + 0.5 times 0.91 = 0.95. That procedure is perfectly valid by the unconditional criterion, and on the day you actually used instrument A it reports a level of 95 when the truthful figure is 99. Fisher's conditionality principle says: condition on the ancillary, analyse the experiment you performed. Bayesian analysis does this without being told, because the posterior is computed from the likelihood of the data in hand, and the likelihood already knows which instrument was used.

Two smaller cases are worth having in your pocket. The Wald interval for a proportion with zero successes in ten trials is 0 +/- 1.96 sqrt(0), the single point 0, which has zero coverage whenever theta is not zero. And in variance-component models an interval built by inverting a test can come out empty, which is a valid 95 percent procedure and a useless report.

Two ways to cut a posterior

Once you have a posterior there is still a choice. An equal-tailed interval puts 2.5 percent in each tail and is invariant under monotone reparameterisation: transform to the log-odds and the endpoints transform with it. A highest posterior density interval collects the shortest region holding 95 percent, so every point inside has higher density than every point outside, and it is the natural choice for a skewed posterior or one pressed against a boundary. It is not invariant under reparameterisation, since the shortest interval for theta is not the shortest interval for log theta.

For Beta(8, 17) the two differ by two thousandths in width, 0.3525 against 0.3546, and nobody would care. For a posterior that piles up at zero, such as the Beta(1, 11) of Lesson 5, the equal-tailed interval starts at 0.0023 and excludes the most probable value, while the HPD interval starts at 0. That is when the choice matters, and when it matters you should say which one you used.

Common misconceptions

• "There is a 95 percent probability that the parameter is in this confidence interval." Under the frequentist definition theta is a fixed constant, so it is either in the interval or not, and the probability is 0 or 1. The 95 percent describes the procedure across datasets. This is the single most common misstatement in applied science, and it is a Bayesian sentence attached to a frequentist object.
• "The two kinds of interval always differ, so you can tell them apart by their numbers." For a normal mean with a flat prior they are identical to every decimal place, and for a binomial proportion the Wilson interval and the flat-prior credible interval agree to three decimals. Identical numbers, different claims.
• "Guaranteed coverage means the interval is trustworthy on any given day." Welch's interval has exactly 50 percent coverage while being certain on some datasets and 5 percent reliable on others. Guaranteed coverage is a guarantee about an average.
• "An HPD interval is always better than an equal-tailed one." It is shorter by construction, but it is not preserved under reparameterisation, so an HPD interval for a rate and an HPD interval for its logarithm are inconsistent statements. For a roughly symmetric posterior the difference is not worth the trouble.

Pulling it together

A confidence interval is a statement about a procedure over repeated datasets; a credible interval is a probability statement about the parameter given these data. Numerically they often coincide: for Beta(8, 17) the flat-prior credible interval is 0.156 to 0.511 and the Wilson interval is 0.156 to 0.509, and for a normal mean with a flat prior they agree exactly. Welch's two uniform observations separate them completely: the interval from the smaller to the larger reading has 50 percent coverage for every theta, yet is certain to contain theta whenever the two readings differ by more than half a unit, and has conditional coverage d/(1 - d), which is 0.053 when the readings differ by 0.05. The range is ancillary, and conditioning on it, which Bayesian analysis does automatically, is what turns an average into a statement about the data you have. Equal-tailed intervals are reparameterisation-invariant, HPD intervals are shortest, and the difference matters only when the posterior is skewed or against a boundary.

That is the technical case. The next module puts it in front of the people who spent forty years arguing about it, and computes the exact point at which a p-value of 0.05 becomes evidence in favour of the null hypothesis.

Sources

1. Wikipedia contributors. (n.d.). Credible interval. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Confidence interval. Wikipedia. en.wikipedia.org
3. Welch, B. L. (1939). On confidence limits and sufficiency, with particular reference to parameters of location. Annals of Mathematical Statistics, 10(1), 58-69.
4. Cox, D. R. (1958). Some problems connected with statistical inference. Annals of Mathematical Statistics, 29(2), 357-372.
5. Brown, L. D., Cai, T. T., & DasGupta, A. (2001). Interval estimation for a binomial proportion. Statistical Science, 16(2), 101-133.

In March 1935 Jerzy Neyman read a paper on the design of agricultural experiments to the Royal Statistical Society, and Ronald Fisher, who had been asked to propose the vote of thanks, used his reply to say that he had hoped for a paper on a subject the author understood. Fisher was 44 and had invented most of the field. Neyman was 41, had been at University College London for a year, and had just published, with Egon Pearson, a framework that treated Fisher's significance tests as a special case of something more general. They never reconciled. Harold Jeffreys, meanwhile, was a Cambridge geophysicist who had spent 1932 to 1934 arguing with Fisher in the Proceedings of the Royal Society about whether probability could describe a state of belief, and who in 1939 published a book laying out how to do statistics if it could.

Three positions, three sets of problems that produced them. This lesson lays them out and then computes the single number where two of them openly disagree.

What each of them was actually trying to do

Fisher spent fourteen years at Rothamsted Experimental Station analysing field trials. A field trial happens once. You want to know what this experiment, on this soil, in this season, has shown, and Fisher's significance test is built to answer that: compute how improbable the data would be if nothing were happening, and if the answer is small enough, treat the null as discredited. He insisted that a p-value is evidence about the experiment in hand, not a rule for behaviour, and he never accepted a fixed threshold as anything but a convenience. He also spent thirty years trying to get probability statements about parameters without a prior, through fiducial inference, and the attempt failed: for anything beyond the simplest location problem the fiducial distribution is not a probability distribution and different routes to it give different answers.

Neyman and Egon Pearson came from a different problem. If you are inspecting batches of manufactured parts, or screening thousands of compounds, the long run is not a fiction, it is your Tuesday. Their 1933 framework asks you to fix in advance a rule for choosing between two hypotheses and to characterise it by two error rates: how often it rejects a true null, and how often it fails to reject a false one. Neyman called the result inductive behaviour, deliberately refusing to call it inference. Their machinery gives power calculations, sample size planning and, from 1937, confidence intervals, and none of it claims to say anything about the case in front of you.

Jeffreys was trying to decide whether the Earth's core is liquid, whether the continents move, and how the surface waves from an earthquake in Japan constrain the elastic constants of the mantle. In that work you have a handful of specific physical hypotheses and one planet, and the question that matters is which hypothesis is probably true. He wanted a number for that, so he built the apparatus this course teaches: priors chosen by invariance rules, posteriors, and Bayes factors for comparing hypotheses.

Jeffreys' objection to significance testing is worth stating precisely, because it is sharper than the usual complaint. A p-value measures the data against a tail region: everything at least as extreme as what happened. So a hypothesis can be discredited partly because of results it did not predict and which also did not occur. Fisher's objection to Jeffreys was equally specific: a lump of prior probability on an exact point value is a claim about the world that nobody can justify, and the answer depends on how much of it you put there.

Key idea: both objections are correct, and both are about the same structural fact. A p-value must integrate over unobserved data; a Bayes factor must integrate over unobserved parameter values. Neither can avoid supplying something that the data do not contain.

Lindley's paradox, worked at three sample sizes

Dennis Lindley published the sharpest version of the disagreement in 1957. Set it up concretely. You observe n independent measurements from N(theta, 1) and want to test H0: theta = 0. The alternative H1 says theta is unknown, with a prior N(0, 1). Give the two hypotheses equal prior probability.

Suppose the data come in at exactly z = ybar sqrt(n) = 1.96, so the two-sided p-value is exactly 0.05 no matter what n is. Under H0 the sample mean is N(0, 1/n). Under H1, averaging over the prior, it is N(0, 1 + 1/n). The Bayes factor in favour of the null is the ratio of those two densities at the observed ybar, which reduces to

BF01 = sqrt(1 + n) exp[ -(z2/2) n/(1 + n) ]

Read the last row. A result that any journal would report as significant at the five percent level is, on the same data and the same model, evidence of 146 to 1 in favour of the null hypothesis it supposedly rejects. The crossover, where the two verdicts change places, is at n = 41.6 for this prior.

Why it happens, in one sentence and then in detail

Holding the p-value fixed while increasing n means holding ybar sqrt(n) fixed, which forces the observed effect to shrink like 1/sqrt(n). At n = 10,000 the sample mean is 0.0196. The null predicted 0 and the alternative, with its N(0, 1) prior, expected something of order 1. An observation of 0.0196 is very close to what the null predicted and nowhere near what the alternative predicted, so of course it favours the null.

Split the formula and the mechanism is visible. The factor exp[-(z2/2) n/(1+n)] tends to exp(-1.92) = 0.146 and stops moving; that is the part that penalises the null for the data being 1.96 standard errors away. The factor sqrt(1 + n) grows without limit; that is the Ockham penalty on the alternative for spreading its prior mass over a region sqrt(n) times wider than the data can resolve. A hypothesis that permits many outcomes must give each of them a small probability, and it pays for that whenever the outcome that occurs is one the sharper hypothesis also allowed.

That also shows what the paradox depends on. Change the alternative's prior standard deviation from 1 to 0.25 and at n = 10,000 the Bayes factor falls from 14.65 to 3.68; change it to 4 and the factor rises to 58.6. The number is a joint statement about the data and about how vague you allowed the alternative to be, and it does not converge to a data-only answer as n grows. Lesson 16 makes that dependence the whole subject.

In short: the paradox is not a mistake by either side. It shows that a p-value and a posterior probability answer different questions, and that the gap between them widens with sample size rather than closing.

Where they agree, which is most of the time

None of this affects estimation much. With a flat prior on a normal mean the credible interval and the confidence interval have the same endpoints, as Lesson 6 showed. The Bernstein-von Mises theorem says that under regularity conditions the posterior converges to a normal distribution centred at the maximum likelihood estimate with variance equal to the inverse Fisher information, so in large samples with a fixed number of parameters the two traditions produce the same intervals and differ only in what they claim about them.

The disagreement is concentrated in a short list: testing a point null, especially with a large sample; small samples, where the prior does real work; problems with many parameters, where Lessons 12 and 13 show the frequentist estimator can be beaten outright; and sequential or multiple testing, where the two traditions do not even agree on which quantity needs adjusting. James Berger argued in 2003 that a conditional frequentist test, one that reports an error probability conditioned on the strength of the observed evidence, reproduces the Jeffreys answer and would have satisfied Fisher's demand for a data-specific statement. That reconciliation exists on paper and has never been adopted in practice.

Common misconceptions

• "A p-value of 0.05 means there is a 5 percent chance the null is true." The table above gives the actual posterior probability of the null at exactly p = 0.05 under a specific model: 0.33 at n = 5, 0.60 at n = 100, 0.94 at n = 10,000. The p-value is P(data at least this extreme | null); the quantity people want is P(null | data), and Lesson 1 already showed those two can differ by orders of magnitude.
• "Fisher was a frequentist and Jeffreys was a Bayesian, so Fisher and Neyman were on the same side." Fisher and Neyman disagreed with each other at least as violently as either disagreed with Jeffreys. Fisher's significance test with no alternative, no power and no fixed threshold is not the Neyman-Pearson procedure, and the hybrid taught in most first courses is a merger neither man would have endorsed.
• "Lindley's paradox shows Bayesian methods are biased toward the null." It shows that a point null with a lump of prior probability, compared with a diffuse alternative, is favoured when the observed effect is tiny. Replace the point null with a small interval around zero, or make the alternative sharper, and the numbers change accordingly. The paradox is about how you framed the comparison.
• "With enough data the two traditions must converge." They converge for estimation. For testing a point null they diverge, which is exactly what the sqrt(1 + n) factor says.

The takeaway

Fisher built significance tests for one-off field trials and wanted a statement about the experiment in hand; Neyman and Pearson built decision rules with controlled error rates for repeated industrial sampling; Jeffreys built Bayes factors because a geophysicist needs the probability that a specific hypothesis about the Earth is true. Lindley's 1957 example makes the gap arithmetic: with z = 1.96 fixed so that p = 0.05 always, and a N(0, 1) prior on the alternative, the Bayes factor for the null runs 0.494 at n = 5, crosses 1 at n = 41.6, and reaches 14.65 at n = 10,000 and 146 at a million. The mechanism is the Ockham factor sqrt(1 + n), the price the alternative pays for spreading its prior over a range the data have long since narrowed. Change the alternative's prior standard deviation to 0.25 or 4 and the factor at n = 10,000 becomes 3.68 or 58.6, which is the warning that Lesson 16 will collect.

A posterior probability of 0.936 for the null is still not a decision. The next lesson supplies the missing ingredient, which is a statement of what it costs to be wrong in each direction, and works a treatment decision to a threshold of one sixteenth.

Sources

1. Wikipedia contributors. (n.d.). Lindley's paradox. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Ronald Fisher. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Jerzy Neyman. Wikipedia. en.wikipedia.org
4. Lindley, D. V. (1957). A statistical paradox. Biometrika, 44(1-2), 187-192.
5. Berger, J. O. (2003). Could Fisher, Jeffreys and Neyman have agreed on testing? Statistical Science, 18(1), 1-32.

The posterior probability of cancer after one positive mammogram is 0.0776. That number does not tell anyone whether to biopsy. Two people can agree on it exactly and still disagree about what to do, and their disagreement will not be about probability at all. It will be about what a missed cancer costs relative to an unnecessary biopsy, and until somebody writes that ratio down, the argument cannot be settled.

This lesson supplies the missing half. A decision problem needs three things: a posterior, a set of available actions, and a loss function L(a, theta) giving the cost of taking action a when the truth is theta. The rule is then mechanical. Compute the expected loss of each action under the posterior, and take the action with the smallest one.

E[L(a)] = integral L(a, theta) p(theta | data) d theta, then choose a minimising it.

Three losses, three estimates

Start with the simplest decision: report a single number for theta. What number depends entirely on how you will be penalised for being off.

• Squared loss, L = (a - theta)2. Expand the expectation: E[(a - theta)2] = (a - E[theta])2 + Var[theta]. The second term does not involve a, so the minimum is at a = E[theta], the posterior mean.
• Absolute loss, L = |a - theta|. Differentiating gives P(theta < a) - P(theta > a) = 0, so the minimum is at the posterior median.
• Zero-one loss, L = 0 if |a - theta| < eps and 1 otherwise. As eps shrinks the minimiser is the posterior mode, which for a flat prior is the maximum likelihood estimate.

Check it on the shuttle posterior, Beta(8, 17), whose mean is 0.3200, median 0.3151 and mode 0.3043.

Each column has its minimum where the theory says it should, and the differences are tiny because this posterior is nearly symmetric. Take a skewed posterior, say Beta(1, 11) from Lesson 5 with mean 0.0833, median 0.0611 and mode 0, and the three answers differ by more than a factor of the smallest.

Bottom line: the posterior mean is not the estimate. It is the estimate under squared loss, which is a choice about how you will be penalised, made most often by not making it.

Treat or wait, and the threshold that falls out

Now a two-action problem with a binary state. The patient either has cancer, with posterior probability p, or does not. The actions are to treat now or to wait and rescan in six months.

The units are arbitrary; only the ratio matters, and the ratio here says a missed cancer is fifteen times as bad as an unnecessary intervention. Now compute:

E[L(treat)] = p times 0 + (1 - p) times 1 = 1 - p

E[L(wait)] = p times 15 + (1 - p) times 0 = 15p

Treat when 15p > 1 - p, that is when 16p > 1, that is when p > 1/16 = 0.0625. In general the threshold is cFP/(cFP + cFN), and it is worth memorising in the form: the threshold is one over one plus the loss ratio.

At p = 0.0776 the woman is above the threshold, so the decision is to treat, at an expected loss of 1 - 0.0776 = 0.9224 against 1.164 for waiting. Notice how close 0.0776 is to 0.0625. A small change in the loss ratio, from 15 to 12, moves the threshold to 0.0769 and the decision flips. That sensitivity is a genuine feature of the problem and not a defect of the method: it is telling you that this patient is a genuinely marginal case and that the answer turns on a value judgement, which is exactly what a clinician would say.

What is a second mammogram worth?

Suppose you can order a second, conditionally independent reading before deciding. It has the same operating characteristics: sensitivity 0.80, false positive rate 0.096. Is it worth having? Answer by computing the expected loss with the test and comparing it with the expected loss without.

Start from p = 0.0776. The probability the second reading is positive is

0.0776 times 0.80 + 0.9224 times 0.096 = 0.06208 + 0.08855 = 0.15063

E[loss with the test] = 0.15063 times 0.58786 + 0.84937 times 0.27408 = 0.08855 + 0.23280 = 0.32135

Value of the test = 0.92240 - 0.32135 = 0.60105

The second reading is worth 0.601 loss units, which in the units of this table is worth about three fifths of an unnecessary intervention. If the test costs less than that, in money, radiation, delay and anxiety converted to the same scale, take it. If it costs more, do not.

Compare that with perfect information. An oracle that told you the truth would let you treat every genuine case and leave every healthy woman alone, for an expected loss of exactly 0. So the expected value of perfect information is 0.9224 - 0 = 0.9224, and the second mammogram captures 0.601/0.9224 = 65 percent of everything an oracle could offer. That ceiling is useful: it tells you at once that no further imaging, however good, can be worth more than 0.92 units here, so a proposed test costing 1.5 units can be refused without knowing anything about its accuracy.

The upshot: the value of a test is not its accuracy. It is the reduction in expected loss it produces, which is zero whenever the test cannot change the decision. A test with a likelihood ratio of 50 is worthless to a patient already at p = 0.9 if the treatment threshold is 0.0625, because she will be treated either way.

Why every good rule turns out to be a Bayes rule

A decision rule is inadmissible if some other rule is at least as good for every value of theta and strictly better somewhere. Abraham Wald proved in the late 1940s that, under regularity conditions, the class of admissible rules is essentially the class of Bayes rules and their limits. That is a striking result to arrive at from the frequentist side: it says that if you insist only on not being dominated, you end up using a procedure that behaves as though it had a prior, whether or not you admit to one.

The rival criterion is minimax, which chooses the rule with the smallest worst-case risk over theta. Minimax is the right answer when nature really is adversarial, in cryptography or in some engineering safety margins, and it is usually too pessimistic elsewhere: it optimises against a value of theta that your posterior says is very unlikely. A minimax rule is often a Bayes rule for the least favourable prior, which is another way of saying the same thing.

Common misconceptions

• "Decision theory means putting a monetary value on a human life." It means writing down a ratio you are already acting on. A clinic that biopsies above 6 percent has chosen a loss ratio of 15 whether or not anyone in it says so. The choice is between an explicit ratio someone can argue with and an implicit one nobody can audit.
• "The best estimate is the posterior mean." Only under squared loss. Under absolute loss it is the median, under zero-one loss the mode, and in the treat-or-wait problem the best action is not an estimate of theta at all.
• "A more accurate test is always worth ordering." A test is worth exactly the reduction in expected loss it buys. If no possible result would change the action, that reduction is zero however good the test is.
• "Minimax is safer than a Bayes rule because it protects against the worst case." It protects against a worst case chosen by an adversary who knows your rule. When theta is a patient's disease status rather than an opponent's move, that protection is bought by performing worse in every situation you are actually likely to face.

Summing up

Expected loss under the posterior turns a distribution into an action. Squared loss selects the posterior mean, absolute loss the median, zero-one loss the mode, and on Beta(8, 17) the three candidates give expected squared losses of 0.008369, 0.008393 and 0.008614, with the minimum where the theory places it. In a two-action problem the loss ratio fixes a threshold: with a missed case fifteen times worse than a false alarm the threshold is 1/16 = 0.0625, so a posterior of 0.0776 says treat, at an expected loss of 0.9224 against 1.164 for waiting. A second conditionally independent mammogram splits the future into a positive branch with probability 0.15063 and posterior 0.412, and a negative branch with probability 0.84937 and posterior 0.0183, for an expected loss of 0.32135 and a value of 0.601 loss units, which is 65 percent of the 0.9224 an oracle would be worth. Wald's complete class theorem says that the admissible rules are essentially the Bayes rules, so a decision maker who only insists on not being dominated has adopted a prior without saying so.

Every posterior so far has come from a conjugate model or a two-state problem, and every integral has had a closed form. Almost no real model does. The next module replaces exact integration with sampling, and starts by asking how many samples you need before the third decimal place means anything.

Sources

1. Wikipedia contributors. (n.d.). Bayes estimator. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Loss function. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Admissible decision rule. Wikipedia. en.wikipedia.org
4. Berger, J. O. (1985). Statistical Decision Theory and Bayesian Analysis (2nd ed.). Springer.
5. Robert, C. P. (2007). The Bayesian Choice (2nd ed.), Chapter 2. Springer.

Ten thousand draws from the shuttle posterior put the probability that the O-ring rate exceeds one half at 0.0320. Another ten thousand draws, from the same distribution with a different seed, put it at 0.0335. The exact answer, from the incomplete beta function, is 0.031957. Neither run is wrong, and the difference between them is the subject of this lesson: how much of the number you just printed is signal.

Everything you want is an expectation

Look at what a posterior summary actually is. The mean is E[theta]. The variance is E[theta2] - E[theta]2. The probability that theta exceeds one half is E[1(theta > 0.5)], the expectation of an indicator. A quantile is the value at which one of those expectations hits a target. The predictive probability for the next flight is E[theta]. The expected loss of an action is E[L(a, theta)].

All of these have the form E[g(theta)] = integral g(theta) p(theta | y) d theta, and the law of large numbers says that if you can draw theta(1), ..., theta(S) from the posterior, the average of g over those draws converges to the integral. That is the entire idea of Monte Carlo integration, and its virtue is that the cost does not grow with the dimension of theta, unlike every deterministic quadrature rule, which is why a hierarchical model with 300 parameters is tractable at all.

The error bar on a simulation

With independent draws the central limit theorem gives the error directly:

MCSE = sd(g(theta))/sqrt(S)

For a probability, g is an indicator with standard deviation sqrt(p(1 - p)), so MCSE = sqrt(p(1-p)/S). For the shuttle posterior, with p = 0.03196 and posterior standard deviation 0.09148:

Now the two runs at the top make sense. At S = 10,000 the standard error on the tail probability is 0.00176, and the two answers 0.0320 and 0.0335 differ by 0.0015, which is less than one standard error. Both are correct to two decimal places and neither is correct to three.

The square root is unforgiving. Every extra correct digit costs a factor of one hundred in draws, so a run of 1,000 that gives you two digits needs 10 million to give you four. A working rule: report a posterior summary to the precision at which the MCSE is smaller than about half the last digit you print, and print the MCSE alongside so a reader can check.

Remember: a posterior mean printed to five decimal places from 1,000 draws is a claim about the random number generator, not about the parameter.

Getting the draws: rejection sampling

Suppose you can evaluate the target density f(theta) up to a constant but cannot sample from it. Rejection sampling needs a proposal density q you can sample from and a constant M with f(theta) <= M q(theta) everywhere. Draw theta from q, draw u uniform on the interval from 0 to 1, and keep theta if u <= f(theta)/(M q(theta)). The kept draws are exact draws from f, with no approximation at all.

Work it on Beta(8, 17) with a uniform proposal. The density peaks at the mode 7/23 = 0.30435, and its height there is 4.28076, so take M = 4.28076. The acceptance probability is 1/M = 0.2336, which means about 4.3 uniform draws per accepted sample, and you would need roughly 42,800 proposals to accumulate 10,000 posterior draws.

That is fine in one dimension. It is fatal in twenty. If the target is a 20-dimensional normal and the proposal is a normal with 1.2 times the standard deviation in each coordinate, the required M is 1.220 = 38.3, so you accept one proposal in 38, and at 1.5 times the standard deviation the figure is 1.520 = 3325. Acceptance rates fall geometrically in the dimension, and this is precisely the wall that Markov chain methods were invented to get around.

Importance sampling, and the moment it collapses

A second idea keeps every draw. Rewrite the expectation against a proposal you can sample from:

Ep[g] = integral g(theta) [p(theta)/q(theta)] q(theta) d theta = Eq[g(theta) w(theta)], with w = p/q.

Since the posterior is usually known only up to its normalising constant, use the self-normalised form:

estimate = sum wi g(thetai) / sum wi

The health of the sampler is measured by the effective sample size:

ESS = (sum wi)2/sum wi2

If every weight is equal, ESS equals the number of draws. If one weight dominates, ESS falls to 1: you have paid for S draws and are using one of them.

Here is that failure, made with five draws you can check by hand. The target is the shuttle posterior, proportional to theta7(1 - theta)16. Suppose you take the proposal to be the Beta(1, 99) prior of Lesson 5, the one that insists O-ring failures are about a one-in-a-hundred event. Because the proposal is the prior, the weights are just the likelihood. Five draws from Beta(1, 99) came out at 0.031, 0.019, 0.008, 0.004 and 0.002.

sum w = 1.039673 and sum w2 = 1.001565, so ESS = 1.0396732/1.001565 = 1.079. The importance estimate of the posterior mean is 0.0305, against a true value of 0.3200. The estimator is off by a factor of ten and would still be off by a factor of ten with five million draws, because the proposal essentially never visits the region where the posterior lives.

Compare a sane proposal. Uniform draws at 0.31, 0.52, 0.18, 0.63 and 0.09 give normalised weights of 0.673, 0.076, 0.237, 0.005 and 0.010, an ESS of 1.94 from five draws, and an estimate of 0.2944 against 0.3200. Still crude, but converging, and it would keep improving with more draws.

What matters here: importance sampling fails silently. The estimate has a number of digits and no complaint attached, and the only way you learn that it is nonsense is by computing the effective sample size or looking at the largest normalised weight.

Diagnosing weight collapse

Three checks, in increasing order of usefulness. First, the largest normalised weight: 0.96 in the table above is a flashing light, since one draw is carrying the whole estimate. Second, the effective sample size relative to S: an ESS below about one percent of S means the answer is unreliable. Third, and best, Pareto smoothed importance sampling, which fits a generalised Pareto distribution to the largest weights and reports its shape parameter k. When k exceeds 0.7 the weight distribution has an infinite variance and the estimate has no usable error bar. That diagnostic returns in Lesson 16, where leave-one-out cross-validation is computed with importance weights and the k value flags the observations the model cannot predict.

Common misconceptions

• "More draws always fixes it." More draws fix Monte Carlo error, which falls as one over the square root of S. They do not fix a proposal that visits the wrong region: the Beta(1, 99) sampler above converges to the right answer only after enough draws to visit the region the posterior occupies, and under that proposal a draw above 0.2 turns up about once in four billion.
• "Monte Carlo error is the same as posterior uncertainty." They are unrelated quantities that happen to share the word error. The posterior standard deviation of 0.0915 is what the data leave unknown; the MCSE of 0.0029 is how badly you have approximated your own posterior. Only the second one is fixed by buying a bigger computer.
• "Rejection sampling is a general-purpose method." It is exact and it is unusable in more than a few dimensions, because the acceptance rate falls geometrically with dimension. That is why the rest of this module is about Markov chains.
• "An effective sample size of 1.08 just means the estimate is a bit noisy." It means the estimate is one draw from the proposal, dressed up. The reported figure of 0.0305 is not a noisy version of 0.32; it is a completely different number, and no error bar computed from those weights will tell you so.

What you now know

Every posterior summary is an expectation, and the average of g over posterior draws estimates it with a Monte Carlo standard error of sd(g)/sqrt(S). For the shuttle posterior that is 0.00176 on a tail probability at 10,000 draws, which is why 0.0320 and 0.0335 are the same answer, and each additional correct digit costs a hundredfold increase in S. Rejection sampling from Beta(8, 17) with a uniform proposal accepts 1/4.28076 = 23.4 percent of proposals and is exact, but its acceptance rate collapses geometrically with dimension. Importance sampling reuses every draw with weight p/q and reports its own health through ESS = (sum w)2/sum w2: with the Beta(1, 99) prior as the proposal, five weights give an ESS of 1.079 and an estimate of 0.0305 for a posterior mean of 0.3200. Check the largest normalised weight, the ESS, and where available the Pareto k, because none of these methods complains when it fails.

Independent draws are the problem. The next lesson gives them up, constructs a chain in which each draw depends on the last, and shows why that trade is worth making, one Metropolis step at a time.

Sources

1. Wikipedia contributors. (n.d.). Monte Carlo method. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Importance sampling. Wikipedia. en.wikipedia.org
3. Vehtari, A., Simpson, D., Gelman, A., Yao, Y., & Gabry, J. (2015). Pareto smoothed importance sampling. arXiv. arxiv.org
4. Robert, C. P., & Casella, G. (2004). Monte Carlo Statistical Methods (2nd ed.), Chapters 2-3. Springer.

In 1953 five people at Los Alamos published a four-page paper in the Journal of Chemical Physics about the equation of state of a two-dimensional gas of hard disks. They had 224 disks in a square box and a machine called MANIAC, and they could not evaluate the configuration integral, so they invented a way to avoid it: move one disk by a small random displacement, always accept a move that lowers the energy, and accept a move that raises it with probability exp(-delta E/kT). Run that long enough and the configurations you visit are distributed according to the Boltzmann distribution you could not integrate. W. K. Hastings generalised the rule in 1970 in Biometrika.

What a Markov chain is being asked to do

You want draws from p(theta | y), and you can evaluate it only up to a constant, because the marginal likelihood in the denominator is the integral you cannot do. Markov chain Monte Carlo gives up on independent draws and builds a chain in which each state depends only on the previous one, arranged so that its stationary distribution is the posterior. The states you visit are then, marginally, draws from the target. They are correlated, which costs effective sample size, and Lesson 11 measures that cost.

Metropolis-Hastings, in five lines

Given the current state theta and a proposal density q(theta' | theta):

1. Draw a candidate theta' from q(theta' | theta).
2. Compute r = [pi(theta') q(theta | theta')] / [pi(theta) q(theta' | theta)], pi being the unnormalised target.
3. Draw u uniform on the interval from 0 to 1.
4. If u < min(1, r), move to theta'; otherwise stay, and record theta again as the next draw.

Two things deserve emphasis. The normalising constant of pi cancels in the ratio, which is the whole point. And a rejected step still produces a draw: the current state is written down a second time.

Eight steps, worked

Take a posterior with no closed form. You observe one measurement y = 3 from N(theta, 1), and your prior is a standard Cauchy, heavy-tailed and not easily bullied by a single observation. The unnormalised posterior is

pi(theta) = [1/(1 + theta2)] exp[-(3 - theta)2/2], so log pi(theta) = -log(1 + theta2) - (3 - theta)2/2.

Use a symmetric random walk proposal, theta' = theta + Uniform(-1, 1). Symmetry makes the proposal densities cancel, so r = pi(theta')/pi(theta), the original 1953 form. Start at theta = 0. The first uniform in each row sets the step, as 2(u1 - 0.5); the second decides acceptance.

Check one row. Step 4 proposes 0.76 from 1.06. log pi(0.76) = -log(1.5776) - (2.24)2/2 = -0.4558 - 2.5088 = -2.9647, and log pi(1.06) = -log(2.1236) - (1.94)2/2 = -0.7532 - 1.8818 = -2.6349. The difference is -0.3298, so r = exp(-0.3298) = 0.719. The uniform came up at 0.67, below 0.719, so a downhill move was accepted. A chain that only ever moved uphill would climb to the mode, stop, and sample nothing.

Six of eight steps were accepted and the chain has walked from 0 to 2.04. The posterior has mean 2.285 and mode 2.260, so eight steps have taken it most of the way from an arbitrary start to the region that matters, and none of those eight draws is a posterior sample. That initial stretch is the burn-in, or warm-up, and it is discarded.

So what?: the accept-reject rule does one job. It makes the chain spend time in each region in proportion to the target's density there, using only ratios, which is why the impossible constant never appears.

Why the rule is the rule

The condition that guarantees the target is stationary is detailed balance: for every pair of states, the flow from theta to theta' must equal the flow back.

pi(theta) q(theta' | theta) A(theta, theta') = pi(theta') q(theta | theta') A(theta', theta)

Substitute A = min(1, r). If r < 1 then A(theta, theta') = r while the reverse ratio exceeds 1 and A(theta', theta) = 1, so the left side collapses to pi(theta') q(theta|theta'), which is the right side. The case r >= 1 is the same computation read backwards.

Step size: the only knob, and it matters

Make the proposal too small and nearly everything is accepted, but consecutive draws are almost identical and the effective sample size collapses. Make it too large and almost everything is rejected, so the chain sits still. Roberts, Gelman and Gilks showed in 1997 that for a random walk Metropolis on a high-dimensional target the asymptotically optimal acceptance rate is 0.234, and in one dimension the optimum is around 0.44. An acceptance rate of 0.95 or of 0.02 both mean the same thing: retune.

Gibbs sampling: every step accepted

When you can sample from each full conditional, Gibbs sampling updates one coordinate at a time, drawing it from its conditional given the current values of the others. It is Metropolis-Hastings with the full conditional as the proposal, so the acceptance ratio is always exactly 1.

Take a bivariate normal with zero means, unit variances and correlation 0.8. The conditionals are x | y ~ N(0.8y, 0.36) and y | x ~ N(0.8x, 0.36), standard deviation 0.6. Start at (3.0, 3.0) and use standard normal draws -0.5, 0.8, 0.3, -1.2, 0.6, -0.1 in order.

Three sweeps have carried the chain from (3.0, 3.0) to (1.01, 0.74). Gibbs never rejects, which sounds like pure gain and is not: every move is parallel to an axis, so on a correlated target the chain zigzags along the ridge. For a bivariate normal the lag-one autocorrelation of a Gibbs chain is rho2 per sweep, so at rho = 0.8 it is 0.64 and the integrated autocorrelation time is (1 + 0.64)/(1 - 0.64) = 4.56. At rho = 0.99 those become 0.9801 and 99.5, so a million sweeps buy about ten thousand independent draws. Reparameterising to break the correlation beats buying computing time.

Hamiltonian Monte Carlo: proposals with a sense of direction

Random walk proposals are blind. Hamiltonian Monte Carlo uses the gradient of the log posterior. Define a potential energy U(theta) = -log pi(theta), invent a momentum vector p with kinetic energy pTM-1p/2, and note that the joint density proportional to exp(-U - K) has the target as its theta-marginal. Draw a fresh momentum, then follow Hamilton's equations, which conserve the total energy and so move a long way at essentially constant density.

The trajectory is simulated with the leapfrog integrator, which is reversible and volume-preserving, both required for the Metropolis correction to be valid. One leapfrog step of size eps is:

p <- p - (eps/2) grad U(theta); then theta <- theta + eps M-1 p; then p <- p - (eps/2) grad U(theta).

Run L of these, negate the momentum, and accept with probability min(1, exp(Hstart - Hend)). Because the integrator nearly conserves H, that probability stays near 1 even for a proposal far from the start. The payoff is the scaling: to hold acceptance fixed, random walk Metropolis needs a step size shrinking like d-1, while HMC needs only d-1/4. The No-U-Turn Sampler of Hoffman and Gelman picks L by extending the trajectory until it doubles back, and that is what Stan runs by default. The price is that gradients are required, so discrete parameters must be summed out, and when the geometry is bad the energy error explodes and the step is recorded as a divergence, which Lesson 11 takes up.

Worth holding on to: all three samplers target the same posterior and differ only in how efficiently they explore it. Choosing between them is an engineering decision about geometry, not a statistical one about your model.

Common misconceptions

• "Rejected proposals should not be recorded." They must be. Staying put is how the chain gives extra weight to a high-density region; dropping duplicates changes the stationary distribution.
• "A high acceptance rate means the sampler is working well." A rate near 1 usually means the proposal is far too small, so the chain explores slowly while looking healthy. Aim near 0.44 in one dimension and 0.234 in many.
• "Gibbs has no tuning parameters, so it is the safe choice." Its efficiency is fixed by the correlation structure you wrote down. At rho = 0.99 it needs a hundred sweeps per independent draw, and no tuning will help; reparameterising will.
• "HMC is just Metropolis with a cleverer proposal, so it gives a different answer." It gives the same answer, satisfying the same detailed balance condition with the same target, and gives it for far less work.

Recap

Metropolis-Hastings proposes a move, computes r = pi(theta')q(theta|theta')/[pi(theta)q(theta'|theta)], and accepts with probability min(1, r), so the unknown normalising constant cancels. Eight hand steps on the Cauchy-prior posterior for y = 3 accepted six proposals and carried the chain from 0.00 to 2.04, step 4 taking a downhill move because r = 0.719 and the uniform came up 0.67. The rule satisfies detailed balance, which with irreducibility and aperiodicity gives convergence. Step size is the whole tuning problem, optimal acceptance being near 0.234 in high dimensions and 0.44 in one. Gibbs always accepts: three sweeps at rho = 0.8 moved (3.00, 3.00) to (1.01, 0.74), and its autocorrelation time of 4.56 becomes 99.5 at rho = 0.99. Hamiltonian Monte Carlo follows the gradient instead, scaling as d-1/4 rather than d-1, at the cost of requiring differentiability.

You now have chains. The next lesson asks how you know when to believe one.

Sources

1. Wikipedia contributors. (n.d.). Metropolis-Hastings algorithm. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Gibbs sampling. Wikipedia. en.wikipedia.org
3. Neal, R. M. (2011). MCMC using Hamiltonian dynamics. arXiv. arxiv.org
4. Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler. arXiv. arxiv.org
5. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6), 1087-1092.

Four chains, two thousand draws each, ninety seconds of compute, and the summary table reports R-hat = 2.36 for the parameter you care about. The posterior mean printed next to it is 5.50 with a standard deviation of 3.86. Both numbers are meaningless, and this lesson works out how you know that, what the diagnostic is actually measuring, and what to do next.

The symptom, before any statistic

Plot the draws against iteration number. A healthy chain looks like a fat hairy caterpillar: it moves up and down within a fixed band, and different chains lie on top of one another. Here is what came back instead, with two chains of ten draws shown so the arithmetic stays hand-sized.

Chain A wanders around 3, chain B wanders around 8, and neither ever visits the other's territory. Pooling them produces a posterior that is bimodal for no reason the model supplied. Every summary computed from the pooled draws, mean, interval, tail probability, is a statement about how long each chain happened to be run.

R-hat, computed by hand

Gelman and Rubin's 1992 diagnostic makes that comparison quantitative by asking a simple question: is the variation between chains larger than the variation within them? With m chains of n draws, chain means xbarj and grand mean xbar, define

B = [n/(m - 1)] sum (xbarj - xbar)2 and W = (1/m) sum sj2

Then form the overestimate of the target variance that mixes them, and take its ratio to W:

var+ = [(n - 1)/n] W + (1/n) B, and R-hat = sqrt(var+/W).

Run it on the table. Chain A has mean 3.0; its deviations are 2.2, -1.9, 1.5, -0.9, 0.7, -2.2, 1.9, -1.5, 0.9, -0.7, whose squares sum to 24.0, so sA2 = 24/9 = 2.6667. Chain B has mean 8.0, and its squared deviations also sum to 24.0, so sB2 = 2.6667. Therefore W = 2.6667.

The grand mean is 5.5, so sum (xbarj - xbar)2 = (3.0 - 5.5)2 + (8.0 - 5.5)2 = 6.25 + 6.25 = 12.5, and B = (10/1) times 12.5 = 125.0.

var+ = 0.9 times 2.6667 + 0.1 times 125.0 = 2.4000 + 12.5000 = 14.9000

R-hat = sqrt(14.9000/2.6667) = sqrt(5.5875) = 2.364

Read that number as a ratio of standard deviations. It says the posterior spread estimated from the pooled chains is 2.36 times the spread each chain individually explores, so if the chains ran forever the reported interval would shrink by roughly that factor. The classical advice was to accept R-hat below 1.1; the current recommendation, after Vehtari and colleagues re-examined the statistic in 2021, is 1.01, because 1.05 already corresponds to a meaningful bias in tail quantiles.

The core of it: R-hat compares chains with each other, so it can only detect problems that make chains differ. Four chains started from the same point, or one long chain, can be badly wrong with an R-hat of 1.00.

What went wrong, and the run after the fix

In this case the parameter was a group-level standard deviation with a flat prior on the interval from 0 to infinity, and the chains had separated because two of them were stuck near the boundary. Putting a half-normal prior with scale 5 on it removed the boundary attraction, and the reparameterisation described below removed the rest. Here is the same parameter after the fix.

Chain A now has mean 4.8 with squared deviations summing to 8.80, so sA2 = 8.80/9 = 0.9778; chain B has mean 5.2 with squared deviations summing to 5.60, so sB2 = 5.60/9 = 0.6222. Then W = (0.9778 + 0.6222)/2 = 0.8000. The grand mean is 5.0, so B = 10 times [(4.8 - 5.0)2 + (5.2 - 5.0)2] = 10 times 0.08 = 0.80. Then var+ = 0.9(0.80) + 0.1(0.80) = 0.80, and R-hat = sqrt(0.80/0.80) = 1.000. The chains now disagree exactly as much as a single chain disagrees with itself, which is what convergence looks like.

Effective sample size: how much of your sample is real

R-hat says the chains agree. It says nothing about how many independent draws they represent. For a chain with autocorrelation rhot at lag t,

ESS = N/(1 + 2 sumt>=1 rhot), and the denominator is the integrated autocorrelation time.

Suppose four chains of 1,000 draws each, so N = 4,000, with autocorrelations decaying geometrically at 0.85 per lag.

Those ten values sum to 4.551, giving a truncated estimate 1 + 2(4.551) = 10.10 and ESS = 4000/10.10 = 396. Carrying the geometric series to infinity, sum = 0.85/0.15 = 5.667, so the true autocorrelation time is 1 + 2(5.667) = 12.33 and ESS = 4000/12.33 = 324. Truncating too early flatters the sampler, which is why practical estimators use Geyer's initial positive sequence rule, stopping when consecutive pairs of autocorrelations sum to a negative number.

Two rules to carry. Report bulk-ESS, for the centre of the distribution, and tail-ESS, for the 5th and 95th percentiles, because a chain can be fine in the middle and hopeless in the tails. And treat 400, which is 100 per chain across four chains, as a floor rather than a target: the Monte Carlo error on a posterior mean with ESS 324 is sd/sqrt(324) = sd/18, so a posterior standard deviation of 3.9 leaves 0.21 of pure simulation noise on your estimate.

Divergences, and the geometry that causes them

Hamiltonian Monte Carlo reports a diagnostic no other sampler has. When the leapfrog integrator hits curvature too sharp for its step size, the simulated energy blows up instead of being conserved, and the transition is flagged as a divergence. A run with divergences is not merely slow; the sampler is systematically failing to enter a region of the posterior, and the bias does not shrink with more iterations.

The standard cause is a funnel. Take the hierarchical structure of the next lesson: a group-level scale tau with log tau free to roam, and group effects thetaj ~ N(mu, tau2). When tau is 10 the theta values need steps of order 10; when tau is 0.01 they need steps of order 0.01. There is no single step size that works in both places, so the sampler either crawls in the wide part or diverges in the neck.

The repair is algebraic, not computational. Write thetaj = mu + tau zj with zj ~ N(0, 1), and sample the z values instead. Their geometry no longer depends on tau, the funnel is gone, and divergences usually vanish with it. This is called the non-centred parameterisation, and it is the single most useful trick in applied hierarchical modelling. It is worth noting that when the data are strongly informative about each group, the centred version samples better, so the choice is empirical.

The checklist

• Run at least four chains from dispersed starting points, and discard the warm-up.
• Look at the trace plots and the rank plots before any number.
• Require R-hat below 1.01 on every parameter, not just the ones you care about.
• Require bulk-ESS and tail-ESS above 400, and report the Monte Carlo standard error.
• Require zero divergences. One divergence in four thousand is not a rounding error.
• Check that trajectories are not saturating the maximum tree depth, which signals inefficiency rather than bias.

Why this matters: none of these diagnostics can prove a chain has converged. They are all tests that a broken chain can fail, and passing them means only that the specific failures they detect are absent.

Common misconceptions

• "R-hat is close to 1, so the chains have converged." R-hat detects disagreement between chains. Four chains initialised identically, or a target with a mode so isolated that no chain ever finds it, produce R-hat near 1 while missing most of the posterior.
• "More iterations will fix the divergences." A divergence is a region the sampler cannot integrate through at the current step size. Running longer collects more biased draws. Reducing the step size sometimes helps; reparameterising usually does.
• "Thinning the chain removes autocorrelation and improves the estimate." Thinning throws away draws. Keeping every draw always gives an estimate at least as good; thinning is a memory-saving measure, not a statistical improvement.
• "An effective sample size of 324 from 4,000 draws means the sampler failed." It means each independent draw costs about twelve iterations, which for a correlated posterior is ordinary. It fails when 324 is not enough for the precision you need, which is a question about the Monte Carlo standard error rather than about the ratio.

Looking back

Two chains that wander around 3 and around 8 give W = 2.6667, B = 125.0, var+ = 14.9 and R-hat = 2.364, which says the pooled spread is 2.36 times what a single chain explores. After fixing the prior and the parameterisation the same computation gives W = 0.80, B = 0.80, var+ = 0.80 and R-hat = 1.000. Effective sample size measures a different failure: with autocorrelations of 0.85t, the integrated autocorrelation time is 12.33 and 4,000 draws are worth 324 independent ones, while truncating the sum at lag 10 would have flattered it to 396. Divergences in Hamiltonian Monte Carlo mark curvature the integrator cannot follow, usually a funnel created by a hierarchical scale, and the cure is the non-centred parameterisation thetaj = mu + tau zj rather than more iterations. Every one of these is a test a broken chain can fail, and none is proof that a chain is right.

The funnel that caused the divergences is not a computational curiosity. It is the shape of the model that eight American schools produced in 1981, and that model is the subject of the next lesson.

Sources

1. Wikipedia contributors. (n.d.). Effective sample size. Wikipedia. en.wikipedia.org
2. Stan Development Team. (n.d.). MCMC sampling. Stan Reference Manual. mc-stan.org
3. Vehtari, A., Gelman, A., Simpson, D., et al. (2021). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667-718. arxiv.org
4. Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472.
5. Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv. arxiv.org

In 1981 Donald Rubin published a table with eight rows. The Educational Testing Service had run a randomised experiment in eight American high schools to see whether a short coaching programme raised scores on the verbal section of the SAT, and each school had reported an estimated effect together with its standard error. The schools were labelled A to H, and the numbers are in points on a scale where the whole test ran from 200 to 800.

School A gained 28 points. School C lost 3. If you were the superintendent, would you buy school A's programme? Every standard error in that table is large enough that the differences between schools are within noise, and yet the schools are not identical either. This lesson is about the model that lives between those two statements.

Two answers, both wrong

The first answer is no pooling: take each school at face value. School A's effect is 28 with a standard error of 15, so a 95 percent interval runs from -1 to 57. The estimate is unbiased for school A and it is also, obviously, mostly noise. If you funded programmes by picking the largest estimate you would be running a lottery in which the prize goes to whichever school had the smallest sample and the luckiest draw.

The second is complete pooling: assume every school has the same true effect and combine the eight estimates by precision weighting, exactly as in Lesson 4. The inverse variances sum to 0.060312, so the pooled estimate is

mu-hat = [sum yj/sigmaj2]/[sum 1/sigmaj2] = 0.463533/0.060312 = 7.686, with standard error 1/sqrt(0.060312) = 4.07.

That is a much sharper number, 7.7 points with an interval from -0.3 to 15.7, and it is bought by asserting something nobody believes: that the coaching programme, the teachers and the students at school A are interchangeable with those at school C. If that assertion is wrong the pooled estimate is biased for every individual school.

The model that sits between them

A hierarchical model refuses to choose. Let each school have its own true effect, and let those effects themselves be draws from a common distribution:

yj ~ N(thetaj, sigmaj2) for the observed estimate, and thetaj ~ N(mu, tau2) for the true school effects.

The parameter tau is the one that does the work. It measures how much schools genuinely differ. Set tau = 0 and every school has the same effect, which is complete pooling. Let tau go to infinity and the schools tell you nothing about each other, which is no pooling. Anything between is partial pooling, and, crucially, tau is estimated from the data rather than assumed.

Conditional on mu and tau, the posterior for one school is exactly the normal-normal update of Lesson 4, with the group distribution playing the part of the prior:

E[thetaj | y] = [yj/sigmaj2 + mu/tau2]/[1/sigmaj2 + 1/tau2], with posterior standard deviation 1/sqrt(1/sigmaj2 + 1/tau2).

Define the shrinkage factor Bj = sigmaj2/(sigmaj2 + tau2), the fraction of the distance from yj back to mu that the estimate travels. Then E[thetaj] = (1 - Bj) yj + Bj mu. A school with a large standard error, or a world in which schools are similar, gets pulled hard toward the centre.

Key idea: partial pooling is not a compromise between two positions. It is what you get by writing down a model in which schools are similar but not identical, and then doing the arithmetic.

Turning the dial

Hold mu = 8, close to the pooled estimate, and compute the posterior means at four values of tau.

Check the entry for school A at tau = 10. The precisions are 1/225 = 0.004444 from the data and 1/100 = 0.01 from the group distribution, so

E[thetaA] = (28 times 0.004444 + 8 times 0.01)/(0.004444 + 0.01) = (0.124444 + 0.080000)/0.014444 = 14.15

and the shrinkage factor is 225/(225 + 100) = 0.692, so the estimate moves 69 percent of the way from 28 down to 8. Notice school B, whose estimate of 8 coincides with mu, does not move at all at any value of tau: shrinkage pulls toward the centre, and something already at the centre stays put. Notice also that the posterior standard deviation for school A falls from 15 to 8.32: borrowing strength from the other seven schools buys precision as well as moving the point estimate.

Which column is right? That is what estimating tau decides. In the full Bayesian analysis, with a flat prior on tau over the positive half-line, the marginal posterior for tau is highest at zero and falls away with a long right tail, putting most of its mass below 15. Averaging over that posterior, school A's estimated effect comes out close to 10 rather than 28, with a posterior standard deviation near 8. The data are telling you that the eight schools are much more alike than the raw table suggests, and that a 28-point estimate from one school with a standard error of 15 is mostly luck.

Why this is not a fudge

The obvious objection is that shrinking school A's estimate is biased for school A, and it is. The reply is that the total squared error over the eight schools is smaller, often much smaller, and Lesson 13 makes that claim precise and computes the factor. The subtler defence is about selection: you are looking at school A precisely because it had the biggest number, and the biggest of eight noisy numbers is systematically too big. Partial pooling corrects that automatically, without anyone having to remember to apply a multiplicity adjustment.

The same structure covers a very long list of problems, and it is worth recognising it in the wild: several hospitals with different surgical mortality rates, several laboratories measuring the same constant, several regions in an epidemiological map, several batches of a manufactured part, several respondents each answering several questions. Wherever the units are similar but not identical, the hierarchical model is the honest description, and the alternative is a choice between ignoring the differences and ignoring the similarity.

The computational catch, which you have already met

Small values of tau squeeze all eight thetas toward mu, so the region of the posterior near tau = 0 is a narrow neck, and larger tau opens out into a wide bowl. That is the funnel from Lesson 11, and the eight schools model is the example on which it was first widely diagnosed. Fitting the centred form in Stan produces divergences that concentrate at small tau, and the posterior for tau is biased away from zero because the sampler cannot get into the neck. Rewriting the model as

thetaj = mu + tau zj with zj ~ N(0, 1)

gives an identical model with a geometry the sampler can handle, and the divergences disappear. With eight groups and one observation each, the non-centred form is the right choice; with a hundred observations per group it would not be.

One more warning. A flat prior on tau over the positive half-line is proper here only because there are eight groups; with three or fewer the posterior can be improper, and a half-normal or half-Cauchy prior on tau is the standard remedy.

Common misconceptions

• "Partial pooling averages the two extremes." It does not interpolate between the no-pooling and complete-pooling answers by a fixed rule. Each school gets its own weight from its own standard error, which is why school C at 16 shrinks further than school E at 9 from the same distance.
• "Shrinkage biases the estimates, so it should be avoided in a regulated setting." Each individual estimate is biased toward the centre and the collection of estimates is better, in total squared error and in calibration. Which one you should optimise is a decision-theory question, and Lesson 8 gave you the machinery for it.
• "The hierarchical model assumes all schools are the same." That is complete pooling, the special case tau = 0. The hierarchical model estimates how different they are, and if the data say they differ a great deal, tau is large and almost no shrinkage occurs.
• "You need many groups for a hierarchical model to be worth it." Eight is enough to be useful and few enough that the prior on tau matters. Below about five groups tau is very poorly identified and the prior is doing most of the work, which is a reason to state it clearly rather than a reason to abandon the model.

Where this leaves us

Eight schools produced coaching effects from -3 to 28 with standard errors from 9 to 18. Complete pooling gives 7.686 with a standard error of 4.07 and asserts the schools are interchangeable; no pooling reports 28 for school A with a standard error of 15 and asserts they have nothing to do with one another. The hierarchical model yj ~ N(thetaj, sigmaj2) with thetaj ~ N(mu, tau2) estimates the between-school standard deviation and shrinks each school by Bj = sigmaj2/(sigmaj2 + tau2). At mu = 8 and tau = 10, school A moves from 28 to 14.15 with its standard deviation falling from 15 to 8.32, while school B, already at the centre, does not move at all. The full analysis puts the between-school standard deviation low and school A's effect near 10. The price is a funnel-shaped posterior, cured by the non-centred parameterisation.

Shrinking a set of estimates toward their common mean sounds like a Bayesian preference. The next lesson shows that it is a theorem, proved in 1961 by two frequentists, and works it on six baseball players whose later performance we can check.

Sources

1. Wikipedia contributors. (n.d.). Bayesian hierarchical modeling. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Multilevel model. Wikipedia. en.wikipedia.org
3. Stan Development Team. (n.d.). Stan User's Guide. mc-stan.org
4. Rubin, D. B. (1981). Estimation in parallel randomized experiments. Journal of Educational Statistics, 6(4), 377-401.
5. Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.), Chapter 5. Chapman and Hall/CRC.

Roberto Clemente's first 45 at-bats of the 1970 season produced 18 hits, a batting average of .400. Fifteen other major league players had also come to bat exactly 45 times by the same date, and Bradley Efron and Carl Morris used all eighteen of them to ask a question with a checkable answer: given a player's first 45 at-bats, what will he hit over the rest of the season? The rest of the season eventually happened, so every prediction can be scored.

Take six of the eighteen, spanning the range.

The obvious prediction is the observed average: predict .400 for Clemente and .156 for Alvis. Nobody who follows baseball would do that, because nobody hits .400 over a season and nobody who is playing in the major leagues hits .156. But the obvious estimator is unbiased, it is the maximum likelihood estimate, and until 1956 there was a general belief that in a problem like this you could not do better. That belief was false.

Building the alternative

Write ybar for the mean of the six averages: (0.400 + 0.356 + 0.289 + 0.244 + 0.222 + 0.156)/6 = 0.27778. Take the sampling variance of a single average as sigma2 = pbar(1 - pbar)/45 = 0.27778 times 0.72222/45 = 0.004458, so sigma = 0.0668, the standard error attached to any one player's 45 at-bats.

Now measure how far apart the six observed averages are. The sum of squared deviations from their mean is

S = (0.12222)2 + (0.07822)2 + (0.01122)2 + (-0.03378)2 + (-0.05578)2 + (-0.12178)2 = 0.040247

The James-Stein estimator, shrinking toward the common mean, is

theta-hatj = ybar + c (yj - ybar), with c = 1 - (k - 3) sigma2/S.

Here k = 6, so c = 1 - 3 times 0.004458/0.040247 = 1 - 0.33231 = 0.66769. Every player keeps two thirds of his distance from the group mean.

The James-Stein predictions have 3.31 times less total squared error. Efron and Morris report a ratio of about 3.5 across the full set of eighteen players. And look at the individual rows: James-Stein is better for all six here, though it does not have to be. It is guaranteed only in total.

The point: shrinkage did not require knowing anything about baseball. The only inputs were the six numbers, their spread, and the standard error of a proportion from 45 trials.

Stein's result, stated carefully

In 1956 Charles Stein proved that for estimating the mean vector of a multivariate normal with known variance, under total squared error loss, the sample mean is inadmissible in three or more dimensions. In 1961 Willard James and Stein produced the explicit estimator that beats it. Inadmissible means what Lesson 8 defined: there exists another estimator whose expected loss is never larger, for any true mean vector, and strictly smaller somewhere. For the James-Stein estimator, strictly smaller is everywhere.

In one or two dimensions the sample mean is admissible and no such improvement exists. The transition at three is not a technicality about proofs; it is where the estimator's variance reduction begins to outweigh the bias it introduces, and it is why Stein's example was received as a paradox rather than a lemma.

The genuinely uncomfortable part is that the components need not have anything to do with each other. Estimate the speed of light, the wheat yield in Nebraska and the mean weight of a Hungarian schoolchild, all standardised, and shrinking the three toward their common average lowers the expected total squared error. Nothing in the theorem asks the quantities to be related. What the theorem optimises is the sum of the three errors, and if that sum is not the quantity you care about, the result is a curiosity. If you are choosing which of eighteen players to sign, or which of eight schools to fund, the sum is exactly what you care about.

The same estimator, derived from a hierarchical model

Now put Lesson 12's model on it: yj ~ N(thetaj, sigma2) and thetaj ~ N(mu, tau2). The posterior mean is (1 - B) yj + B mu with B = sigma2/(sigma2 + tau2), and a Bayesian would put a prior on tau and integrate. Empirical Bayes instead estimates tau from the data and plugs it in.

Marginally, integrating out the thetas, yj ~ N(mu, sigma2 + tau2), so S follows a chi-square distribution on k - 1 degrees of freedom scaled by sigma2 + tau2. The obvious plug-in is tau-hat2 = S/(k - 1) - sigma2, which here gives 0.040247/5 - 0.004458 = 0.003591 and a shrinkage factor of 0.004458/0.008049 = 0.554.

But that plug-in is not what you want. The quantity to estimate is B itself, and since E[1/chi-squarek-1] = 1/(k - 3), the estimator

B-hat = (k - 3) sigma2/S

is exactly unbiased for B. Substitute the numbers: 3 times 0.004458/0.040247 = 0.3323, so 1 - B-hat = 0.6677. That is the James-Stein constant, arrived at from a completely different direction. James-Stein is empirical Bayes with the shrinkage factor estimated without bias, and the difference between 0.554 and 0.332 is the whole gap between a careless plug-in and a careful one.

One more detail. If the six averages had been very close together, S would be small and c = 1 - (k-3)sigma2/S could come out negative, which would flip every estimate to the far side of the mean. The positive-part James-Stein estimator sets c = max(0, c), and it dominates the original, so there is no reason ever to use the unmodified version.

What empirical Bayes leaves out

Plugging in a single value of tau treats it as known. It is not: with six groups, tau-hat2 = 0.003591 has enormous sampling variability, and pretending otherwise makes every interval too narrow. Carl Morris quantified the correction in 1983, and full Bayesian analysis handles it by integrating over the posterior for tau, which is exactly what Lesson 12 described and what a Stan model does without being asked.

So the ordering is: no pooling ignores the group structure; empirical Bayes uses it but treats the estimated group variance as certain; full Bayes propagates the uncertainty. The point estimates from the last two are usually close. The intervals are not, and it is intervals that get quoted.

Bottom line: James-Stein is not a Bayesian result and does not need a prior, and the hierarchical model that reproduces it is not a trick for getting a smaller number. They are two derivations of the same estimator, and the agreement is why shrinkage survived a dispute that settled almost nothing else.

Common misconceptions

• "James-Stein beats the sample mean for every component." It beats it in total expected squared error. An individual component can be worse, and a player whose true average really is .400 is estimated worse by shrinking. In this dataset all six happen to improve, which is luck on top of a guarantee about the sum.
• "Shrinkage only makes sense when the quantities are genuinely related." The theorem holds for arbitrary unrelated components. Relatedness matters for whether the total squared error is the loss you should be minimising, not for whether the estimator dominates.
• "Empirical Bayes and full Bayes give the same answer." They typically give similar point estimates and different intervals, because empirical Bayes conditions on an estimated tau as though it were known. With few groups the difference is large.
• "The estimator works because .400 is implausible for a batting average." No outside knowledge of baseball enters the formula. The shrinkage constant comes from comparing the observed spread of the six averages, 0.040247, with the spread 45 at-bats would produce by chance alone, 0.004458 per player.

What to carry forward

Six players batting in 45 at-bats give averages from .156 to .400 with a common mean of .27778 and a per-player sampling variance of 0.004458. Their squared deviations sum to S = 0.040247, so the James-Stein constant is c = 1 - 3(0.004458)/0.040247 = 0.6677 and each estimate keeps two thirds of its distance from the mean. Total squared error against the rest of the season falls from 0.014239 to 0.004296, a factor of 3.31, with Efron and Morris reporting about 3.5 over all eighteen players. Stein showed in 1956 that the sample mean is inadmissible in three or more dimensions and James and Stein produced the estimator in 1961. The same constant appears as 1 - (k-3)sigma2/S, the unbiased empirical Bayes estimate of the hierarchical shrinkage factor, against the careless plug-in value of 0.554. Empirical Bayes gets the point estimates and understates the intervals, because it treats an estimated tau as known.

So far the parameters being shrunk have been group means. The next lesson shrinks regression coefficients instead, and computes the exact coefficient size at which shrinkage stops being a good idea.

Sources

1. Wikipedia contributors. (n.d.). James-Stein estimator. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Stein's example. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Empirical Bayes method. Wikipedia. en.wikipedia.org
4. Efron, B., & Morris, C. (1975). Data analysis using Stein's estimator and its generalizations. Journal of the American Statistical Association, 70(350), 311-319.
5. Morris, C. N. (1983). Parametric empirical Bayes inference: Theory and applications. Journal of the American Statistical Association, 78(381), 47-55.

Halve a regression coefficient whose true value is 1.6 standard errors from zero and your estimate improves. Halve one whose true value is 1.8 standard errors from zero and your estimate gets worse. The break-even point is exactly the square root of three, and this lesson derives it, then shows why the choice of prior is really a choice about what happens to the coefficients on the far side of that line.

The cleanest possible setting

Take an orthonormal design with the error standard deviation known and equal to 1, so that the least squares estimate of each coefficient is an independent observation of it:

yj = betaj + ej, ej ~ N(0, 1).

Put a normal prior betaj ~ N(0, tau2) on each coefficient. The normal-normal update from Lesson 4 gives a posterior mean of c yj with

c = tau2/(1 + tau2).

Now compute the mean squared error of the estimator c y for a coefficient whose true value is beta. The bias is (c - 1)beta and the variance is c2, so

MSE(c y) = (1 - c)2 beta2 + c2

The unshrunk estimator has MSE = 1. Shrinking is worth it exactly when (1-c)2beta2 + c2 < 1, that is when

beta2 < (1 - c2)/(1 - c)2 = (1 + c)/(1 - c)

With tau2 = 1, so that the estimate is halved, the threshold is sqrt(3) = 1.732. Below it shrinkage helps and above it hurts, and the damage grows quadratically: at beta = 5 the halved estimator has an MSE of 0.25 times 25 + 0.25 = 6.5 against 1 for doing nothing.

What matters here: a Gaussian prior applies the same fractional shrinkage to every coefficient, no matter how large. That is fine when all the coefficients are small and expensive when one of them is not.

Every penalty is a prior

Two familiar estimators are Bayesian in disguise. Ridge regression minimises

||y - X b||2 + lambda ||b||2

and if you write down the log posterior for a normal likelihood with a N(0, sigma2/lambda) prior on each coefficient, you get exactly that expression with the sign flipped. So ridge is the posterior mode, and because a normal posterior is symmetric it is also the posterior mean.

The lasso minimises ||y - X b||2 + lambda ||b||1, and the absolute-value penalty is the log of a Laplace, or double exponential, density. In the orthonormal case its solution is soft thresholding: bj = sign(yj) max(|yj| - lambda/2, 0), which sets small coefficients exactly to zero.

Here is where the disguise matters. The lasso solution is the posterior mode, and the mode of a Laplace posterior can be exactly zero while the posterior mean is never exactly zero. So the Bayesian lasso of Park and Casella does not perform variable selection: run it and every coefficient comes back with a posterior distribution straddling zero. If sparsity is what you wanted, the Bayesian answer is that the data support a small effect, not no effect, and reporting a zero is a decision made under a particular loss rather than an inference.

The horseshoe: shrink almost everything, and leave the rest alone

The problem with both is the tails. A normal prior has tails that die like exp(-beta2), so a genuinely large coefficient is dragged down hard; a Laplace prior has tails dying like exp(-|beta|), which is better but still applies a constant absolute shrinkage to everything. What you usually want in a wide regression is different treatment for different coefficients: crush the noise, leave the signal.

The horseshoe prior of Carvalho, Polson and Scott does this with a scale mixture:

betaj | lambdaj, tau ~ N(0, tau2 lambdaj2), with lambdaj ~ Cauchy+(0, 1).

Each coefficient gets its own local scale, and the half-Cauchy on that scale has both a spike near zero and a very heavy tail. Write the shrinkage weight kappaj = 1/(1 + tau2lambdaj2), so the posterior mean is (1 - kappaj) yj. With tau = 1 the implied prior on kappa is Beta(1/2, 1/2), which is U-shaped with mass piled at 0 and at 1. The name comes from that shape: a coefficient is either shrunk almost completely or barely shrunk at all, and the data decide which.

Four estimators, one dataset

Ten orthonormal coefficients, error standard deviation 1. The truth is that only the first is real, at beta1 = 3, and the other nine are exactly zero. The least squares estimates come back as

3.2, -0.4, 0.6, 1.1, -0.9, 0.3, -1.4, 0.7, -0.2, 0.5

Read the first column down. Least squares gets the real coefficient nearly right and pays for it by chasing nine pieces of noise. Ridge cuts the noise in half and cuts the signal in half too, which costs 1.96 on that one coefficient alone and leaves the total at 3.30. Lasso zeroes six of the nine noise coefficients and still drags the real one down by a full unit, since soft thresholding subtracts the same amount from everything. The horseshoe crushes the noise to a few hundredths and leaves the signal essentially untouched, for a total error two hundred times smaller than least squares.

That is the whole argument for heavy-tailed shrinkage priors, and it is worth being clear about the conditions. The horseshoe wins here because the truth is sparse. If all ten coefficients had been around 0.5, ridge would win, because a normal prior is the correct description of that world. Choosing a shrinkage prior is a claim about how the coefficients are distributed, and it should be defended as one.

Setting the global scale

The parameter tau controls overall shrinkage and is the hard part. Estimating it from the data with a half-Cauchy prior works, but in problems where the number of predictors exceeds the number of observations the posterior for tau is poorly identified and the heavy tails let weakly identified coefficients run away. Piironen and Vehtari's regularised horseshoe adds a slab: coefficients that escape the shrinkage are still bounded by a wide normal. It also gives a way to set tau from a guess at how many coefficients are non-zero, so a substantive statement, "I expect about five of these forty predictors to matter", becomes a prior scale.

Remember: the lambda in a penalised regression and the tau in a Bayesian shrinkage prior are the same number wearing different clothes, and cross-validating lambda is empirical Bayes with the model comparison hidden inside the loop.

Common misconceptions

• "Shrinkage always improves an estimate." Halving an estimate helps only while the true coefficient is under sqrt(3) standard errors from zero, and the penalty above that grows with the square of the coefficient. The break-even point moves with the prior scale, as the table shows.
• "The Bayesian lasso gives sparse estimates." The lasso point estimate is a posterior mode and can be exactly zero. The posterior mean under a Laplace prior is never exactly zero, and neither is any posterior interval. Sparsity in a Bayesian analysis is a decision, not an inference.
• "Ridge and lasso are frequentist methods, so the prior is a metaphor." The penalties are exactly the negative log densities of a normal and a Laplace prior, and the estimates are exactly posterior modes. Nothing about the correspondence is loose.
• "A heavier-tailed prior is always safer." Heavy tails mean weakly identified coefficients are barely constrained, which is why the plain horseshoe struggles when predictors outnumber observations and why the regularised version exists.

The short version

With an orthonormal design and unit error variance, a N(0, tau2) prior gives a posterior mean of c y with c = tau2/(1 + tau2), and the estimator beats the raw one exactly when beta2 < (1 + c)/(1 - c): 1.414 standard errors at c = 1/3, 1.732 at c = 1/2, 3.000 at c = 0.8. Ridge is the posterior mode under a normal prior and lasso the posterior mode under a Laplace prior, which is why the lasso zeroes coefficients while the Bayesian posterior mean never does. The horseshoe gives each coefficient its own half-Cauchy scale, so the implied shrinkage weight has a Beta(1/2, 1/2) distribution with mass at both ends. On ten coefficients with one true value of 3 and nine zeros, total squared error ran 5.41 for least squares, 3.30 for ridge, 0.81 for lasso and 0.03 for the horseshoe, with ridge estimating the true 3 as 1.60 and the horseshoe as 3.15. The ordering is a consequence of the truth being sparse; a different truth reverses it.

Every model in the last three lessons has been fitted and reported. None of them has been checked against the data that produced it. The next lesson does that, and kills two of them.

Sources

1. Wikipedia contributors. (n.d.). Ridge regression. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Lasso (statistics). Wikipedia. en.wikipedia.org
3. Piironen, J., & Vehtari, A. (2017). Sparsity information and regularization in the horseshoe and other shrinkage priors. arXiv. arxiv.org
4. Park, T., & Casella, G. (2008). The Bayesian lasso. Journal of the American Statistical Association, 103(482), 681-686.
5. Carvalho, C. M., Polson, N. G., & Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2), 465-480.

Simon Newcomb's sixth measurement, made in 1882, came back as minus 44. He was timing how long light took to travel 7,442 metres between a mirror at the Naval Observatory and one at the base of the Washington Monument, and recording each run as a deviation from 24,800 nanoseconds. Sixty-four of his 66 runs fall between 16 and 40. Two do not: one at -2 and one at -44.

Fit a normal model with a flat prior on the mean. The sample mean is 26.212, the sample standard deviation is 10.745, and the posterior for mu is centred at 26.21 with a standard error of 1.323, giving a 95 percent interval of 23.62 to 28.81. The value implied by the modern speed of light is 33.0, which is outside that interval, so something is wrong. This lesson is about how you would have known that from the data alone, without a modern value to check against.

The mechanism

A posterior predictive check asks one question: if this model were true, could it have produced data that look like the data I have? Operationally:

1. Draw a parameter value from the posterior.
2. Simulate a replicated dataset yrep of the same size from the model at that parameter value.
3. Repeat a few thousand times.
4. Pick a test quantity T, compute T(y) on the real data and T(yrep) on each replicate, and see where the real value falls.

The summary is the Bayesian p-value, P(T(yrep) >= T(y)), and a value near 0 or near 1 means the model cannot produce the feature of the data that T measures.

Check one: the mean, which passes

Take T = the sample mean. Simulate 66 observations from N(mu, 10.7452) with mu drawn from its posterior, and take the mean of each replicate. The replicate means cluster around 26.21 with a spread of about 1.3, and the observed 26.212 sits in the middle, giving a Bayesian p-value close to 0.5.

This proves nothing. The model was fitted by matching the mean, so a check on the mean is guaranteed to pass. A useful test quantity is one that the fitting procedure did not target.

Check two: the minimum, which fails by nine orders of magnitude

Take T = the smallest observation. The expected minimum of 66 standard normal draws is about -2.35, so replicated minima should cluster around 26.21 - 2.35 times 10.745 = 0.96, and a spread of a few units either side. The observed minimum is -44.

Put a number on it. For a single draw, P(y < -44) = Phi((-44 - 26.212)/10.745) = Phi(-6.534) = 3.2 x 10-11. For the minimum of 66 independent draws to fall below -44,

P = 1 - (1 - 3.2 x 10-11)66 = 2.1 x 10-9

about one chance in 470 million, and that figure ignores the extra uncertainty in mu, which changes it by nothing that matters. In a simulation of a thousand replicates, none would come close, so the Bayesian p-value is reported as less than 0.001.

The upshot: the normal model is not a slightly imperfect description of Newcomb's data. It is a model under which the dataset in front of you effectively could not have occurred, and every posterior interval it produces inherits that.

Two responses are available, and they are not equivalent. The first is to replace the normal likelihood with a Student t on, say, four degrees of freedom, which accommodates occasional extreme values without letting them drag the centre. The second is to treat -44 and -2 as recording errors, which is what most later analysts have concluded, and to say so in the paper rather than silently deleting them. Refitting the normal model without the two outliers moves the estimate of mu to 27.75, closer to the 33.0 implied by the modern value but still short of it, so the outliers were not the only problem: Newcomb also had a systematic error, as Michelson did in Lesson 4.

A second model, killed a different way

A clinic records the number of unplanned readmissions per patient over a year for 100 patients. The counts total 150, so the mean is 1.50, and a Poisson model is the obvious first choice. Fit it and check.

The mean is matched exactly, as it must be. Everything else is wrong. Take T = the number of zeros. Under Poisson(1.5), P(0) = e-1.5 = 0.22313, so a replicate should contain 22.31 zeros with a standard deviation of sqrt(100 times 0.22313 times 0.77687) = 4.16. The observed count is 45, which is 5.45 standard deviations above, and the Bayesian p-value is again below 0.001.

Take T = the number of patients with five or more readmissions. Under the model P(X >= 5) = 0.01858, so 1.86 expected; 7 observed. Both tails are too heavy and the middle is too light. The variance-to-mean ratio makes the same point in one number: the sample variance is 3.020 against a mean of 1.500, a ratio of 2.01, where a Poisson distribution requires 1.

The diagnosis names the repair. Excess zeros point to a zero-inflated model, in which a fraction of patients are structurally never readmitted and the rest follow a Poisson distribution. Overdispersion in both tails points to a negative binomial, which is the gamma-Poisson mixture of Lesson 4 with the rate varying between patients. In this dataset both features are present, which is common and is why the two models are often compared rather than chosen between. Lesson 16 does that comparison properly.

Which test quantities to use

A check is only as good as its test quantity, and the useful ones are the features the fitting procedure ignored.

A test quantity may also depend on the parameters, T(y, theta), in which case it is called a discrepancy measure and both the observed and the replicated versions move from draw to draw. That is how you check, for instance, whether the residuals are larger than the model's own error scale for a particular subgroup.

The honest limitation

A posterior predictive p-value uses the data twice: once to form the posterior and again to compute T(y). As a result it is not uniformly distributed even when the model is exactly right; it is conservative, clustering toward 0.5, so the check has less power than a frequentist p-value with the same nominal meaning. Bayarri and Berger proposed partial predictive and conditional predictive alternatives that repair the distribution at some cost in complexity, and the practical advice from the people who use these checks most is to treat the number as a rough summary and the plot as the real output.

Which leads to the discipline that makes this work: decide what would surprise you before you look. If you cannot name a feature the model would have to reproduce, you will find yourself looking at a dozen quantities and reporting the one that came out worst, which is a search, not a check.

So what?: a model that fits every summary statistic you fitted it to has told you nothing. A check is informative exactly to the degree that it could have failed.

Common misconceptions

• "A Bayesian p-value near 0.5 means the model is correct." It means the model reproduces that one feature. Newcomb's data give 0.5 on the mean and under 0.001 on the minimum from the same fit.
• "Posterior predictive checks are the Bayesian version of a significance test." They are exploratory. Because the data are used twice, the p-value is conservative and not calibrated, and no threshold should be attached to it.
• "If the check fails, delete the outliers." Sometimes that is right, and it must be stated. Refitting Newcomb's data without the two extreme values moves mu from 26.21 to 27.75, which reduces but does not remove the discrepancy with 33.0, so deleting them would have hidden a second problem rather than fixing the first.
• "Comparing observed and expected counts is just a chi-square test." The table is a display, not the inference. The predictive distribution of the number of zeros under the fitted model carries parameter uncertainty and can be computed for any test quantity, including ones with no classical test attached, such as the longest run of identical values.

Putting it together

A posterior predictive check simulates replicated datasets from the fitted model and asks whether the observed data look like them under a chosen test quantity. On Newcomb's 66 measurements with mean 26.212 and standard deviation 10.745, a check on the mean gives a Bayesian p-value near 0.5 and proves nothing, because the model was fitted to the mean. A check on the minimum gives a p-value below 0.001: the observed -44 is 6.53 standard deviations below the centre, and the minimum of 66 normal draws falls that low with probability 2.1 x 10-9. On 100 patients with 150 readmissions, a Poisson model predicts 22.31 zeros against 45 observed, a discrepancy of 5.45 standard deviations, and 1.86 patients with five or more against 7 observed, with a variance-to-mean ratio of 2.01 where the model demands 1. Excess zeros point to zero inflation, heavy tails to a negative binomial. Useful test quantities are the ones the fit ignored, and the p-value is conservative because the data are used twice.

You now have two candidate repairs for the count model and no way to choose between them. The next lesson supplies three ways, computes them by hand, and shows that one of the three can be moved by a factor of forty thousand without touching the data.

Sources

1. Wikipedia contributors. (n.d.). Posterior predictive distribution. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Simon Newcomb. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Zero-inflated model. Wikipedia. en.wikipedia.org
4. Gelman, A., Meng, X.-L., & Stern, H. (1996). Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6(4), 733-760.
5. Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2013). Bayesian Data Analysis (3rd ed.), Chapter 6. Chapman and Hall/CRC.

Seven of the twenty-three pre-Challenger flights showed O-ring distress. Test the managerial claim that the per-flight rate is one in ten against the alternative that it is something else, and the Bayes factor in favour of one in ten comes out at 1,442. Test the same claim against the same data with a different description of what else it might be, and the Bayes factor comes out at 0.036. Nothing about the flights changed. This lesson computes both numbers, explains the swing, and then gives the tool you should reach for when the swing makes a Bayes factor unusable.

The Bayes factor, and the integral inside it

Bayes' theorem applied to two hypotheses gives

posterior odds = prior odds times BF01, where BF01 = p(y | H0)/p(y | H1).

Each term is a marginal likelihood: the probability of the data under a hypothesis, with every parameter that hypothesis contains integrated out against its prior. For a point null there is nothing to integrate:

p(y | H0) proportional to 0.17 times 0.916 = 1.85302 x 10-8

dropping the binomial coefficient, which appears in both hypotheses and cancels. For an alternative with theta ~ Beta(a, b) the integral is a beta function ratio:

p(y | H1) proportional to B(a + 7, b + 16)/B(a, b)

Now vary only the alternative's prior.

The evidence runs from 1,442 to 1 for the one-in-ten claim down to 28 to 1 against it, a swing of 39,700, driven by a choice nobody in the argument would call the main issue. The last row is a floor: no prior can push BF01 below 0.0255, since the best any alternative can do is put all its mass at the maximum likelihood estimate.

Why the prior refuses to wash out

A posterior becomes insensitive to its prior as the data accumulate. A marginal likelihood does not, and the reason is structural. Apply Laplace's approximation to the integral:

p(y | H1) approximately p(theta-hat) L(theta-hat) sqrt(2 pi/(n I))

The prior enters as its density at the maximum likelihood estimate, and that factor never goes away. A prior that is ten times more spread out halves nothing and divides the marginal likelihood by roughly ten, so the Bayes factor moves by a factor of ten regardless of sample size. Doubling the width of a prior you did not think carefully about doubles your evidence for the null.

That is Lindley's paradox in a different costume, and it explains the top row of the table. Beta(1, 99) puts almost all its mass below 0.05, where 7 of 23 is wildly improbable, so the alternative predicted the data terribly and the point null wins by three orders of magnitude. The Bayes factor is doing what it claims: comparing how well each hypothesis, taken whole with its prior, predicted the data before they arrived.

Why this matters: a Bayes factor is a comparison of two complete probabilistic claims, prior included. If you would not defend the prior in public, you cannot defend the Bayes factor, and no amount of data will rescue it.

When a Bayes factor is the right tool

Jeffreys built the machinery for a setting where the priors are genuinely specified: two physical models of the Earth's interior, each making quantitative predictions with parameters constrained by physics. Modern uses that survive the objection look similar. A genetic linkage analysis whose prior on the recombination fraction comes from the biology. A physics experiment where the alternative predicts a signal of a size the theory fixes. A trial registered in advance with an agreed prior on the effect size.

Where it fails is the routine case: a default prior on a parameter nobody has thought about, in a model chosen for convenience. Kass and Raftery's 1995 review remains the standard warning, and Jeffreys' interpretive scale, in which 1 to 3 is barely worth mentioning, 3 to 10 substantial, 10 to 30 strong, 30 to 100 very strong and above 100 decisive, is a rough guide attached to numbers you can defend.

The predictive alternative, computed

A different question: which model predicts a new observation better? Its answer does not depend on integrating over the prior, because it is asked of the posterior. Define the expected log pointwise predictive density, elpd, as the sum over observations of the log predictive density each would receive from the model fitted to everything else.

Two estimators. WAIC computes, from the posterior draws for the fitted model,

lppd = sumi log(average over draws of p(yi | theta)) and pWAIC = sumi variance over draws of log p(yi | theta),

and reports elpdWAIC = lppd - pWAIC. Leave-one-out cross-validation computes the same target directly, refitting without each observation, which Pareto smoothed importance sampling approximates from the single fit at no extra cost.

Here are four observations and two candidate models: A, a normal likelihood, and B, a Student t. The fourth observation is an outlier.

So WAIC = -2 elpd gives 18.46 for model A and 15.10 for model B, and B is better by 1.68 on the elpd scale. Before believing that, compute the standard error of the difference from the pointwise differences -0.05, -0.05, 0.03, 1.75. Their standard deviation is 0.8874, so

SE = sqrt(4) times 0.8874 = 1.775

The difference is 1.68 with a standard error of 1.78. It is not distinguishable from zero, and reading down the column shows why: three observations are a hair in favour of A and the entire difference comes from observation 4. With four observations you have not learned which model predicts better; you have learned that they differ on one point.

Leave-one-out tells the same story with a warning attached. The pointwise elpdloo values are -1.25, -1.35, -1.85, -4.88 for model A, summing to -9.33 with ploo = 2.10, and -1.30, -1.40, -1.82, -3.05 for model B, summing to -7.57 with ploo = 1.17. The difference is 1.76 with a standard error of 1.86. The Pareto k diagnostic from Lesson 9 comes back at 0.78 for observation 4 under model A, above the 0.7 threshold, so the importance sampling approximation is unreliable for exactly the observation that carries the comparison. That is the diagnostic doing its job: it names the point to refit exactly rather than trust.

Worth holding on to: report the difference in elpd together with its standard error, and look at the pointwise contributions. A model comparison in which one observation supplies the entire difference is a statement about that observation.

Two criteria, two questions

The two disagree systematically because they optimise different things, and neither is a mistake. If a true model is in your candidate set and you want to identify it, Bayes factors are the coherent tool. If your models are all wrong and you want the one that forecasts best, which is the ordinary situation, cross-validation is. When two models are close, the best answer is often neither: stacking combines their predictive distributions with weights chosen to maximise the cross-validated score, and usually beats picking one.

Common misconceptions

• "A Bayes factor is prior-free because the priors cancel." The prior on the hypotheses cancels; the priors on the parameters inside each hypothesis do not, and are the whole source of the swing from 1,442 to 0.036.
• "With enough data the prior stops mattering for a Bayes factor too." Laplace's approximation shows the prior entering as its density at the maximum likelihood estimate, a factor that does not shrink with n. Widening a prior tenfold multiplies the Bayes factor for the null by roughly ten at any n.
• "A lower WAIC means a better model." It means a better estimated out-of-sample predictive fit, on this sample, with an uncertainty attached. A difference of 3.36 in WAIC with a standard error of 3.55 on the same scale is not a result.
• "Cross-validation is not Bayesian." The quantity being estimated, the expected log predictive density of a new observation under the posterior predictive distribution, is defined entirely within the Bayesian model. What cross-validation abandons is the claim that one of the candidates is true.

What to remember

A Bayes factor compares marginal likelihoods, each of which integrates the likelihood against the prior. For 7 distressed flights in 23 against a point null at 0.10, the alternative's prior moves the answer from 1,441.8 under Beta(1, 99) to 0.0363 under Beta(7, 16), with a floor of 0.0255 from a point prior at the maximum likelihood estimate: a swing of nearly forty thousand. Laplace's approximation shows why, since the prior enters as its density at the estimate and that factor does not shrink with n. Bayes factors belong where both priors are genuinely specified and defensible. For predictive comparison, WAIC subtracts an effective-parameter penalty from the log pointwise predictive density, and leave-one-out estimates the same quantity with a Pareto k diagnostic attached: on four observations the two models differ by 1.68 in elpd with a standard error of 1.78, and all of it comes from one point whose k value of 0.78 says the approximation there cannot be trusted. Bayes factors are consistent, cross-validation is efficient, and stacking beats choosing when the candidates are close.

That is the last technique in the course. What remains is the question a colleague will ask when you propose to spend three weeks on a hierarchical model, and the last lesson answers it without a sales pitch.

Sources

1. Wikipedia contributors. (n.d.). Bayes factor. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Cross-validation (statistics). Wikipedia. en.wikipedia.org
3. Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27, 1413-1432. arxiv.org
4. Kass, R. E., & Raftery, A. E. (1995). Bayes factors. Journal of the American Statistical Association, 90(430), 773-795.
5. Watanabe, S. (2010). Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11, 3571-3594.

In October 2015 the American Statistical Association gathered about two dozen statisticians for a two-day meeting to try to write down what a p-value is not. The statement they published in March 2016 runs to six principles, and it is notable for what it does not say. It does not recommend replacing p-values with Bayes factors, or with posterior probabilities, or with anything else. The closest it comes is a sentence noting that other approaches exist. Twenty-odd people who had spent careers on this could agree that the dominant practice was being misused and could not agree on what should replace it.

That is the honest state of the argument, and this lesson closes the course by setting out both sides of the ledger with the numbers you have already computed.

What you get

Probability statements about the thing you asked about. The Beta(8, 17) posterior says the probability that the O-ring rate lies between 0.156 and 0.511 is 0.95, and the probability it exceeds one half is 0.032. A confidence interval cannot say either, as Lesson 6 showed, and the sentence people actually write under a confidence interval is the Bayesian one.

Nuisance parameters handled by integration. Unknown variance, unknown baseline, unknown measurement error: you integrate them out and the uncertainty propagates automatically. There is no separate theory of approximate degrees of freedom.

Sequential analysis without corrections. Lesson 3 updated ten flights and then thirteen and got the same Beta(8, 17) as one batch of twenty-three. Lesson 2 showed why: the stopping rule drops out of the posterior. Group sequential boundaries and alpha spending functions exist to repair something a posterior does not break.

Regularisation as a consequence rather than a bolt-on. Partial pooling in Lesson 12 moved school A from 28 to 14.15 because the model said schools are similar, not because anyone added a penalty. Ridge and lasso turned out in Lesson 14 to be exactly posterior modes under normal and Laplace priors. James-Stein turned out in Lesson 13 to be empirical Bayes with an unbiased shrinkage factor.

Prior information used when data are scarce. Ten patients with zero events support no sensible frequentist point estimate, and the choice among priors in Lesson 5 was visible, arguable and quantified in pseudo-observations.

Predictive distributions that carry the uncertainty. Two future flights both showing damage has probability 0.1108, not 0.32 squared, because the posterior for theta travels with the prediction.

Decisions. Lesson 8 turned 0.0776 into an action through a loss ratio, and computed that a second mammogram is worth 0.601 loss units against a ceiling of 0.9224.

Arbitrary model complexity without new theory. Once you can sample from a posterior, a measurement error model with a mixture likelihood and three levels of hierarchy needs no new distribution theory, only more compute and better diagnostics.

What it costs

Priors, chosen and defended. Beta(2, 18) in Lesson 3 was worth 20 flights against 23 real ones and moved the posterior mean from 0.320 to 0.209. Nothing in the arithmetic warns you. The cost is real and it is largest exactly when the data are thin, which is when you most wanted help.

Computation, measured in hours and in attention. A conjugate update is a line of arithmetic; a hierarchical model with a non-centred parameterisation is a fit you must diagnose. Lesson 11's checklist has six items and each can fail. A frequentist alternative often runs in a second and needs no diagnostics because it has no chain.

Bayes factors that can be moved by a factor of forty thousand. Lesson 16 did exactly that to the O-ring data without touching a single observation, and Lesson 7 showed the same mechanism producing a posterior probability of 0.936 for a null at p = 0.05. If you want to compare hypotheses rather than predict, you have taken on a commitment that a frequentist test does not require.

No error-rate guarantee. A Bayesian procedure under a wrong model can have dreadful long-run properties, and nothing in the posterior reveals it. Regulators who license drugs, and engineers who certify aircraft, often want a procedure with a stated error rate under repetition, and that is a legitimate requirement the framework does not meet by itself.

A confidently wrong answer looks exactly like a confidently right one. Michelson's posterior in Lesson 4 had a standard deviation of 7.35 and an error of 59.9. Newcomb's normal model in Lesson 15 gave a tidy interval for a dataset it could not have generated. Posterior width measures what the model thinks it has learned, and the model is not consulted about whether it is true.

Reporting and review. An honest Bayesian paper carries the priors, their sources, their equivalent sample sizes, a sensitivity analysis, convergence diagnostics and a predictive check. That is a longer methods section, and the number of reviewers equipped to audit it is smaller.

In short: the framework buys coherence and pays for it in specification. Every number the posterior gives you is conditional on choices you made, and the discipline is not avoiding those choices but making them visible.

Where the two traditions agree, with the numbers

Six of the most common tasks in applied statistics, and in all six the numbers match or nearly match. That is why a great deal of practice can proceed without anyone declaring an allegiance, and why the fight is not about most of the work.

Where they do not, and what would settle it

Testing a point null in a large sample. Lindley's paradox is not an empirical disagreement and no experiment settles it. The two procedures answer different questions: one reports how improbable the data are under the null, the other how the data shift the odds between two complete hypotheses. Both answers are correct to their own question, and the lesson is to state which question you asked.

Optional stopping. Here there are factual claims and a simulation settles them. Run a million simulated trials that stop the first time a posterior probability crosses 0.95 and you will find two things: each reported posterior is correctly computed from the data in hand, and the fraction of published confident claims under a true null far exceeds five percent. Both sides are right about different objects, and Lesson 2 named them: the posterior is unaffected, the selection is not.

Many parameters. This one was settled, in 1956, against the intuition of almost everyone. Stein proved the sample mean is inadmissible in three or more dimensions, and Lesson 13 computed the factor of 3.31 on six batting averages. No argument remains about the mathematics; the remaining argument is about whether total squared error is your loss.

A procedure for deciding

Reach for a Bayesian analysis when at least one of these holds: the units are few and structured, so partial pooling buys real precision; there is genuine prior information from a named source; the output must be a probability attached to a hypothesis or a parameter; the decision has asymmetric losses that somebody will have to state; the model is non-standard enough that no frequentist theory exists for it; or the data arrive over time and you want to report continuously.

Do not bother when the sample is large, the model is standard, the prior would be swamped, and the answer will agree to three decimals with a procedure your readers already trust. Lesson 4's horse kicks are a case in point: three different priors gave 0.6139, 0.6125 and 0.6150, and the maximum likelihood estimate was 0.610. Spending an afternoon on that choice would have been a waste of an afternoon.

The core of it: the question is not which framework is true. It is which one answers the question you were asked, at a cost you can pay, in a form the person asking can check.

Common misconceptions

• "Bayesian methods are subjective and frequentist methods are objective." A frequentist analysis chooses a model, a test statistic, an alpha level, a stopping rule and a set of covariates, and every one of those is a judgement. The Bayesian difference is that one more judgement is written in the notation instead of the discussion section.
• "Bayesian methods are now standard, so the dispute is over." Posterior computation is standard; the interpretation is not. The 2016 ASA statement declined to endorse a replacement, and most regulatory and clinical practice still runs on error rates.
• "Since the two approaches usually agree, the choice does not matter." They agree on estimation in large samples and disagree on point null testing, on small samples, on many-parameter problems and on sequential designs. Those four are where the contested results in most literatures live.
• "A Bayesian analysis is safer because it uses all the information." It uses all the information in the model, including whatever the prior asserts. Lesson 15 showed two models producing clean posteriors for data they could not have generated, and no amount of Bayesian machinery substitutes for checking the model against the data.

Looking back

Seventeen lessons, and a single habit: never state a posterior you have not computed. Bayes' theorem turned a positive mammogram at one percent prevalence into 7.8 percent and then 41.2 percent on a second reading. One likelihood function gave two p-values, 0.0730 and 0.0327, under two stopping rules. Twenty-three shuttle flights became Beta(8, 17) with a 95 percent interval of 0.156 to 0.511; 122 horse-kick deaths became Gamma(124, 202) with a mean of 0.6139; Michelson's hundred runs gave a posterior mean of 845.33 that excluded the truth. Three flat priors gave three answers to zero successes in ten trials, and a 50 percent confidence interval turned out to be certain. A p-value of 0.05 favoured the null 14.65 to 1 at n = 10,000; a loss ratio of fifteen set a threshold at one sixteenth and made a second mammogram worth 0.601. Five importance weights gave an effective sample size of 1.08; eight Metropolis steps walked from 0 to 2.04; an R-hat of 2.364 became 1.000. School A fell from 28 to 14.15, six batting averages improved by a factor of 3.31, a horseshoe prior cut total squared error from 5.41 to 0.03, a normal model failed on a minimum of -44, and one Bayes factor moved by a factor of 39,700 on the prior alone.

What you should take away is not a preference. It is the ability to write down a model, state a prior and defend it, compute the posterior and check it, turn it into a decision, and say honestly what in the answer came from the data and what came from you.

Sources

1. Wikipedia contributors. (n.d.). Bayesian inference. Wikipedia. en.wikipedia.org
2. Wikipedia contributors. (n.d.). Frequentist inference. Wikipedia. en.wikipedia.org
3. Wikipedia contributors. (n.d.). Bayesian statistics. Wikipedia. en.wikipedia.org
4. Wasserstein, R. L., & Lazar, N. A. (2016). The ASA statement on p-values: Context, process, and purpose. The American Statistician, 70(2), 129-133.
5. Efron, B. (1986). Why isn't everyone a Bayesian? The American Statistician, 40(1), 1-5.