Primer: Bayesian inference

On By @thinquanautIn Primers

Bayesian inference is a rule for changing your mind. You start with a belief about something you cannot see, a neuron’s true firing rate, say, or which stimulus was on the screen. Data arrive. The rule tells you exactly how much to move your belief, and in which direction, and it leaves you with a new belief that is ready for the next piece of data. This primer builds the rule from one line of arithmetic, applies it to a neuron from the lessons, and then uses it to decode the stimulus the way lesson 3 did, but with an answer that includes how sure it is. It assumes lessons 2 and 3.

The rule

Three ingredients, each a curve over the thing you want to know. The prior: how plausible each possible value was before the data. The likelihood: for each possible value, how probable the data you actually saw would have been. The posterior: how plausible each value is now. Bayes’ rule says the posterior is the prior times the likelihood, rescaled so the total is one. That is the whole thing. Everything else is bookkeeping.

Bayes' rule as three pictures, for the question "what is neuron 141's firing rate on blank trials?" after one trial in which it fired twice in 0.2 seconds. Left: a flat prior. Middle: the likelihood, which peaks at 10 Hz (two spikes in a fifth of a second) but is wide, because two spikes is not much evidence. Right: their product. With a flat prior the posterior is just the likelihood, rescaled.
Bayes’ rule as three pictures, for the question “what is neuron 141’s firing rate on blank trials?” after one trial in which it fired twice in 0.2 seconds. Left: a flat prior. Middle: the likelihood, which peaks at 10 Hz (two spikes in a fifth of a second) but is wide, because two spikes is not much evidence. Right: their product. With a flat prior the posterior is just the likelihood, rescaled.

The computational trick that makes this a two-line program is the grid. Rather than deriving formulas, lay out every value you are willing to consider, here 800 firing rates from 0 to 40 Hz, and compute the three curves as arrays. Multiplying two arrays and dividing by the sum is Bayes’ rule.

import numpy as np
import matplotlib.pyplot as plt
from math import factorial

dat = np.load("steinmetz_session11.zip", allow_pickle=True)
spks = dat["spks"]
area = dat["brain_area"]
contrast = dat["contrast_right"]
dt = float(dat["bin_size"])
t = np.arange(spks.shape[2]) * dt - 0.5
post = (t >= 0.05) & (t < 0.25)
window = post.sum() * dt                                   # 0.2 seconds

counts = spks[141][:, post].sum(axis=1)[contrast == 0]     # neuron 141 on the 167 blank-screen trials
print("first ten counts:", counts[:10])

rates = np.linspace(0.05, 40, 800)                         # every firing rate we are willing to consider, in Hz
prior = np.ones_like(rates)                                # flat: no opinion before seeing data
prior /= prior.sum()

def likelihood(count, rate):
    """Probability of seeing `count` spikes in the window if the true rate is `rate` (Poisson)."""
    expected = rate * window
    return expected ** count * np.exp(-expected) / factorial(count)

like = likelihood(counts[0], rates)                        # one trial's evidence, evaluated at every rate
posterior = prior * like
posterior /= posterior.sum()
print(f"after one trial with {counts[0]} spikes: most probable rate {rates[posterior.argmax()]:.1f} Hz")
# first ten counts: [0 2 1 0 0 0 0 1 0 0]
# after one trial with 0 spikes: most probable rate 0.1 Hz

The first blank trial had zero spikes, so after one trial the most probable rate is as low as the grid allows. One trial is weak evidence; the posterior still spreads up past 18 Hz. The interesting part is what happens next.

Updating, trial by trial

The rule is designed to be applied repeatedly: today’s posterior is tomorrow’s prior. Feed in the 167 blank trials one at a time and watch the belief sharpen. The interval function reads off the range that holds 95 percent of the probability, which is the Bayesian answer to “how sure are you?”.

def update(prior, count):
    posterior = prior * likelihood(count, rates)
    return posterior / posterior.sum()

def interval(p, mass=0.95):
    """The narrowest range of rates holding `mass` of the probability."""
    cdf = np.cumsum(p)
    return rates[np.searchsorted(cdf, (1 - mass) / 2)], rates[np.searchsorted(cdf, 1 - (1 - mass) / 2)]

belief = prior.copy()
snapshots = {}
for i, c in enumerate(counts):
    belief = update(belief, c)                              # today's posterior is tomorrow's prior
    if i + 1 in (1, 5, 20, len(counts)):
        snapshots[i + 1] = belief.copy()
        lo, hi = interval(belief)
        print(f"after {i + 1:3d} trials: mean {np.sum(rates * belief):4.1f} Hz, 95% interval {lo:4.1f} to {hi:4.1f} Hz")
# after   1 trials: mean  5.0 Hz, 95% interval  0.2 to 18.4 Hz
# after   5 trials: mean  4.0 Hz, 95% interval  1.1 to  8.8 Hz
# after  20 trials: mean  2.2 Hz, 95% interval  1.1 to  4.0 Hz
# after 167 trials: mean  7.7 Hz, 95% interval  6.8 to  8.7 Hz

fig, ax = plt.subplots(figsize=(10, 3.8))
for n, p in snapshots.items():
    ax.plot(rates, p / p.max(), label=f"after {n} trials")
ax.set(xlabel="firing rate (Hz)", ylabel="posterior (scaled)", xlim=(0, 20))
ax.legend()
plt.show()
The posterior over neuron 141's blank-trial firing rate after 1, 5, 20 and 167 trials. It narrows as evidence accumulates, from a 95% interval 18 Hz wide to one 2 Hz wide. But notice that it also moves: after 20 trials it is confident the rate is around 2 Hz, and after 167 it is confident the rate is around 8. Both were right about the data they had seen.
The posterior over neuron 141’s blank-trial firing rate after 1, 5, 20 and 167 trials. It narrows as evidence accumulates, from a 95% interval 18 Hz wide to one 2 Hz wide. But notice that it also moves: after 20 trials it is confident the rate is around 2 Hz, and after 167 it is confident the rate is around 8. Both were right about the data they had seen.

Two things to take from this figure. The first is the point of the method: the width of the belief tracks the amount of evidence, automatically, with no separate theory of error bars. The second is a warning that every Bayesian analysis carries. After 20 trials the posterior said, with 95 percent confidence, that the rate was between 1.1 and 4.0 Hz. After 167 trials it said, with the same confidence, 6.8 to 8.7. Both cannot be right about the same quantity, and the resolution is that there is no such quantity: this neuron’s firing rate rose through the session, as the GLM primer found, and our model assumed a single constant rate. The posterior was exactly as confident as the model allowed, and the model was wrong. Bayesian inference tells you what to believe given the model. Checking the model is still your job.

Decoding with uncertainty

Now the question from lesson 3, in Bayesian form: given the spike counts of the 66 visual cortex neurons on one trial, what was the contrast? The unknown is now one of four values rather than a rate on a grid, and the data are 66 counts rather than one. The recipe does not change. The prior is how often each contrast occurred in the training trials. The likelihood of the counts, treating each neuron as an independent Poisson counter with the mean rate it showed at that contrast in training, is a product of 66 terms per contrast, which we compute as a sum of logs so nothing underflows. Multiply, rescale, and the posterior is four probabilities.

visp = np.where(area == "VISp")[0]
X = spks[visp][:, :, post].sum(axis=2).T                   # trials x neurons, as in lesson 3
levels = np.unique(contrast)
rng = np.random.default_rng(0)
order = rng.permutation(len(contrast))
train, test = order[:240], order[240:]

expected = np.array([X[train][contrast[train] == c].mean(axis=0) + 0.1 for c in levels])   # 4 x neurons
prior_c = np.array([np.mean(contrast[train] == c) for c in levels])                        # how common each contrast is

def posterior_over_contrast(counts):
    """P(contrast | counts) for one trial, treating neurons as independent Poisson counters."""
    log_like = (counts * np.log(expected) - expected).sum(axis=1)     # one number per contrast level
    log_post = log_like + np.log(prior_c)
    p = np.exp(log_post - log_post.max())                             # subtract the max before exp, for numerical safety
    return p / p.sum()

P = np.array([posterior_over_contrast(X[i]) for i in test])          # 100 test trials x 4
guess = levels[P.argmax(axis=1)]
truth = contrast[test]
print(f"four-way accuracy: {np.mean(guess == truth):.0%}")
print(f"stimulus-or-blank accuracy: {np.mean((guess > 0) == (truth > 0)):.0%}")
print(f"average confidence when right: {P.max(axis=1)[guess == truth].mean():.2f}, when wrong: {P.max(axis=1)[guess != truth].mean():.2f}")
# four-way accuracy: 78%
# stimulus-or-blank accuracy: 100%
# average confidence when right: 0.98, when wrong: 0.93

fig, axes = plt.subplots(1, 4, figsize=(11, 3.2), sharey=True)
for ax, i in zip(axes, [0, 1, 2, 3]):
    ax.bar([str(c) for c in levels], P[i])
    ax.set(title=f"trial {test[i]}: truth {truth[i]}", xlabel="contrast")
axes[0].set_ylabel("posterior probability")
plt.show()
Left: the posterior over the four contrasts on four test trials. Most trials look like the first two: one bar, and it is right. The third is the least confident of the 100 test trials, torn between 0.25 and 0.5. The fourth is confidently wrong. Right: the confusion matrix. Blank versus stimulus is perfect; the three stimulus contrasts blur into each other.
Left: the posterior over the four contrasts on four test trials. Most trials look like the first two: one bar, and it is right. The third is the least confident of the 100 test trials, torn between 0.25 and 0.5. The fourth is confidently wrong. Right: the confusion matrix. Blank versus stimulus is perfect; the three stimulus contrasts blur into each other.

Three things to read off. First, 78 percent four-way accuracy and 100 percent on stimulus-versus-blank, from a model with no training loop at all, just means and Bayes’ rule; this is the “naive Bayes” classifier, and it is often a strong baseline. Second, the posterior gives you something lesson 3’s decoder did not: a confidence per trial, and the confusion matrix shows exactly where it is uncertain, between 0.25 and 0.5 contrast, where the neurons’ responses genuinely overlap. Third, the warning again. The decoder’s average confidence is 0.98 when right and 0.93 when wrong. A decoder that is wrong 22 percent of the time should not be 93 percent sure when it is wrong. It is overconfident because we assumed the 66 neurons were independent and they are not; they share noise, so 66 neurons carry less evidence than 66 independent ones would. The arithmetic was correct. The assumption was generous, and the posterior inherited its generosity.

What to remember

  • Posterior equals prior times likelihood, rescaled. On a grid, that is two lines of NumPy, and it is the same two lines whether the unknown is a rate, a stimulus, or a set of model weights.
  • The width of the posterior is your uncertainty. It narrows with evidence and nothing else; there is no separate recipe for error bars.
  • Priors matter least when data are plentiful and most when they are scarce, which is exactly when you most need to state them honestly. With 167 trials, the flat prior above was irrelevant.
  • The posterior is only as good as the model. Constant-rate and independent-neuron assumptions produced confident, precise, wrong answers above. Bayes guarantees coherence, not truth.
  • Lesson 3’s decoder is a cousin. Logistic regression with a penalty on large weights is the same as finding the most probable weights under a Gaussian prior; most of machine learning can be read this way.

Exercises

  1. Replace the flat prior with one that believes rates near 2 Hz (lesson 1’s typical neuron), for instance prior = np.exp(-(rates - 2)**2 / 8). How many trials does it take for the data to overrule it?
  2. Split the blank trials into first and second halves and run the update separately on each. Do the two posteriors overlap?
  3. Shuffle the contrast labels of the training trials and rerun the decoder. What accuracy and what confidence do you get, and why is the second number the more alarming one?
  4. Harder: the decoder is overconfident because it treats neurons as independent. Use only the 10 neurons with the largest contrast effect and check whether confidence-when-wrong falls. Then think about why fewer neurons can give better-calibrated answers.

This completes the three primers; the next numbered lesson picks up the Signals track. Data: Steinmetz et al., Nature 2019, CC-BY 4.0, via Neuromatch Academy.

Keep reading