This is the first of three short primers on methods that keep turning up in the lessons, written for someone who has done lesson 2 and lesson 3 and wants the idea behind the tools. Each one is introductory, runs on the same real data as the lessons, and takes about twenty minutes.

A generalised linear model, or GLM, is the thing you reach for when you want to know what a neuron responds to and there is more than one candidate. The tuning curve in lesson 2 answered “how does firing depend on contrast?” by averaging. A GLM answers “how does firing depend on contrast, and on what the mouse did, and on how long the experiment has been running, all at once, with each effect measured while the others are held fixed?” It is the workhorse of systems neuroscience, and the good news is that you have already fitted one: logistic regression in lesson 3 is a GLM.

One recipe, three dishes

Every GLM has the same three parts. A weighted sum of the predictors, exactly as in linear regression: so many points per unit of contrast, so many for turning the wheel, plus a baseline. A link, a fixed function that bends that sum into the range the data can actually occupy. And a noise model, a statement of how the data scatter around the prediction. Change the link and the noise model and you get a different member of the family.

Spike counts are the right-hand panel. They cannot be negative, so the prediction should not be either, which the log link guarantees: the model predicts log(rate), and the rate is exp of that. And their noise is not constant. Look at neuron 141 from lesson 2, counting spikes in the 200 ms after the stimulus:

import numpy as np
import matplotlib.pyplot as plt

dat = np.load("steinmetz_session11.zip", allow_pickle=True)
spks = dat["spks"]
contrast = dat["contrast_right"]
contrast_left = dat["contrast_left"]
response = dat["response"]                 # -1 wheel left, 0 no move, 1 wheel right
dt = float(dat["bin_size"])
t = np.arange(spks.shape[2]) * dt - 0.5
post = (t >= 0.05) & (t < 0.25)

y = spks[141][:, post].sum(axis=1)         # neuron 141's spike count on each trial, the thing we will model
levels = np.unique(contrast)
for c in levels:
    print(f"contrast {c}: mean count {y[contrast == c].mean():5.2f}, variance {y[contrast == c].var():5.2f}")
# contrast 0.0: mean count  1.53, variance  4.13
# contrast 0.25: mean count  7.06, variance 13.66
# contrast 0.5: mean count 10.10, variance 24.17
# contrast 1.0: mean count 13.26, variance 21.34

A straight-line fit with constant noise would treat a two-spike miss on a blank trial and a two-spike miss on a full-contrast trial as equally surprising. They are not. The Poisson model knows that, and that is most of why it fits spike data better.

Fitting one in eight lines

Fitting means choosing the weights that make the observed counts most probable. For a Poisson GLM the gradient of that probability has a form so simple it is worth remembering: for each predictor, add up (observed count minus predicted count) times the predictor, over trials. Walk uphill along that gradient and you arrive at the best weights. This is the same loop as lesson 3’s decoder with exp in place of the sigmoid.

def poisson_glm(X, y, lr=0.05, steps=20000):
    """Fit log(rate) = b0 + X @ b by gradient ascent on the Poisson log-likelihood."""
    X1 = np.column_stack([np.ones(len(y)), X])        # a column of ones for the intercept
    w = np.zeros(X1.shape[1])
    for _ in range(steps):
        rate = np.exp(X1 @ w)                          # the inverse link: weights to a positive rate
        w += lr * X1.T @ (y - rate) / len(y)           # the gradient is just (observed - predicted) times X
    return w

def predict(X, w):
    return np.exp(np.column_stack([np.ones(len(X)), X]) @ w)

X_lin = contrast[:, None]                             # one predictor: contrast as a number
w_lin = poisson_glm(X_lin, y)
print("log rate =", w_lin.round(2), "-> predicted counts:", predict(levels[:, None], w_lin).round(2))

X_cat = np.column_stack([(contrast == c) for c in levels[1:]]).astype(float)   # one column per non-zero contrast
w_cat = poisson_glm(X_cat, y)
print("one-hot ->", predict(np.vstack([np.zeros(3), np.eye(3)]), w_cat).round(2))
# log rate = [1.   1.71] -> predicted counts: [ 2.73  4.18  6.4  15.02]
# one-hot -> [ 1.53  7.06 10.1  13.26]

Two fits of the same neuron. The first uses contrast as a single number and predicts 2.7, 4.2, 6.4 and 15 spikes at the four levels. The real means are 1.5, 7.1, 10.1 and 13.3. The fit is poor, and the reason is instructive: with one weight on contrast, the log link forces the rate to grow exponentially with contrast, and this neuron saturates instead. The second fit gives each contrast level its own weight, and recovers the four means exactly. That is not a coincidence: a tuning curve is a GLM with one weight per condition.

The lesson from the orange curve generalises: the link function imposes a shape on how a numeric predictor acts, and if the real relationship has a different shape you must give the model room, with categories, with a square term, or with a transform of the predictor. A GLM is only as good as the predictors you hand it.

Now add everything else

Here is what averaging cannot do. The mouse in this experiment also saw a stimulus on the left on some trials, turned the wheel left or right or not at all, and got tired or practised as the session went on. Any of those might move neuron 141. Put them all in as columns and fit once.

def log_likelihood(y, rate):
    """How probable the observed counts are under predicted Poisson rates (up to a constant)."""
    return np.sum(y * np.log(rate) - rate)

X_full = np.column_stack([
    X_cat,                                   # contrast on the right, one-hot
    contrast_left > 0,                       # was there anything on the left?
    response == -1,                          # did the mouse turn the wheel left?
    response == 1,                           # or right?
    np.arange(len(y)) / len(y),              # how far into the session (0 to 1)
]).astype(float)
names = ["contrast 0.25", "contrast 0.5", "contrast 1.0", "left stimulus", "turned left", "turned right", "time in session"]
w_full = poisson_glm(X_full, y)
for name, b in zip(names, w_full[1:]):
    print(f"{name:16s} x{np.exp(b):.2f}")
print(f"baseline {np.exp(w_full[0]):.2f} spikes")
print("log-likelihood, contrast only:", log_likelihood(y, predict(X_cat, w_cat)).round(1))
print("log-likelihood, full model:    ", log_likelihood(y, predict(X_full, w_full)).round(1))
# contrast 0.25    x3.77
# contrast 0.5     x5.20
# contrast 1.0     x7.24
# left stimulus    x1.12
# turned left      x1.76
# turned right     x1.52
# time in session  x2.24
# baseline 0.69 spikes
# log-likelihood, contrast only: 2386.7
# log-likelihood, full model:     2449.1

Read the multipliers. A full-contrast stimulus multiplies this neuron’s firing by about seven. A stimulus on the left screen, in the other visual field, does almost nothing, which is what a visual cortex neuron with a receptive field on the right should do. Turning the wheel multiplies firing by 1.5 to 1.8 regardless of direction; this is the movement signal that lesson 3 warned about and lesson 2 saw as the second hump, now measured and separated from the stimulus. And the last line is the surprise: the predictor “time in session” multiplies firing by 2.2, meaning the neuron fires more than twice as much at the end of the hour as at the start. No tuning curve would have shown that, because tuning curves average over time. The log-likelihood, our measure of how probable the data are under the model, improves by 62 units, which is a large amount.

rate_full = predict(X_full, w_full)
fig, ax = plt.subplots(figsize=(10, 3.8))
ax.plot(y, lw=0.8, label="observed count")
ax.plot(rate_full, lw=1.2, label="predicted by the full model")
ax.set(xlabel="trial", ylabel="spikes in the window", xlim=(0, 120))
ax.legend()
plt.show()

What a GLM is for, and what it is not

• It separates effects that co-occur. Stimulus and movement are correlated in this task; the GLM gives each its own weight while accounting for the other. This is the main reason the field uses it.
• Its weights have units and meaning. With a log link, exp(weight) is a multiplier on firing, which is a sentence you can put in a paper.
• It is still linear in the weights. The log link does not let it discover that the contrast response saturates; you have to give it categories. Interactions (does movement matter more at high contrast?) also have to be added by hand as a product column.
• Real neurons are over-dispersed. Their variance is a bit more than their mean, as the second figure showed. The weights are still fine; the confidence you put on them should be a little looser than textbook Poisson theory says.
• Correlated predictors are its weakness. If two columns always move together, the model cannot tell which one matters, and the weights become unstable. Check your design before trusting the weights.

Exercises

1. Fit the full model to a neuron in MOs (secondary motor cortex) or MD (thalamus) instead of 141. Which multiplier is largest there?
2. Add an interaction: a column equal to (contrast == 1.0) * (response != 0). Does movement matter more on high-contrast trials?
3. Replace the Poisson fit with ordinary least squares (np.linalg.lstsq on the same columns) and compare the predictions on blank trials. Where does the linear model go wrong?
4. Harder: fit a separate GLM to every VISp neuron and plot the distribution of the “turned left” multiplier across the population. How many visual neurons carry a movement signal?

Next primer: independent component analysis. Data: Steinmetz et al., Nature 2019, CC-BY 4.0, via Neuromatch Academy.

Independent component analysis, ICA, solves a problem that sounds impossible. Several signals were mixed together before you recorded them, you have only the mixtures, and you want the originals back. It is the cocktail party: two people talking, two microphones, each microphone hearing both. In neuroscience the “people” are an eye blink, a heartbeat, a muscle, a patch of cortex, and the “microphones” are electrodes on the scalp, all hearing everything. This primer shows how ICA does it, writes the algorithm in twenty lines, and uses it to pull the eye blinks out of a real EEG recording. It assumes lesson 4, because ICA starts where PCA stops.

The problem, in two channels

Start with a case where we know the answer. Two sources, a sine wave and a square wave. Two sensors, each recording a different blend of the two. Nobody tells the algorithm the blend.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n = 2000
time = np.linspace(0, 8, n)
s1 = np.sin(2 * np.pi * time)                              # a smooth rhythm
s2 = np.sign(np.sin(2 * np.pi * 1.7 * time + 1))           # a square wave at a different rate
S = np.vstack([s1, s2])                                    # 2 sources x n samples
A = np.array([[1.0, 0.6],                                  # how much of each source reaches each sensor
              [0.4, 1.0]])
X = A @ S + 0.05 * rng.standard_normal(S.shape)            # 2 sensors x n samples: mixtures plus a little noise
print("sensor 1 = 1.0 * s1 + 0.6 * s2;  sensor 2 = 0.4 * s1 + 1.0 * s2")
# sensor 1 = 1.0 * s1 + 0.6 * s2;  sensor 2 = 0.4 * s1 + 1.0 * s2

Why PCA cannot do it, and what can

The obvious tool is PCA from lesson 4: find the directions the data vary along most. But PCA has a built-in assumption that gives the wrong answer here. Its components are always at right angles to each other, and the directions the sources mix along are not at right angles unless you are lucky. PCA will decorrelate the sensors, which removes the linear dependence between them, but a sine and a square wave that are uncorrelated can still be thoroughly mixed.

ICA uses a different clue: independence, which is stronger than uncorrelatedness, and a fact about sums. When you add independent signals together the result is always closer to a bell curve than the ingredients were; this is the central limit theorem doing what it does. So a mixture is more Gaussian than its sources, and the way to unmix is to search for the directions along which the data look least Gaussian. A sine wave and a square wave are both very un-bell-shaped, and in the right directions you see that; in the wrong directions you see a blur.

FastICA, the standard algorithm, does exactly this. First it whitens the data with PCA, which handles the scaling and leaves only a rotation to find. Then it rotates, step by step, to maximise a measure of non-Gaussianity, keeping the rows orthogonal so that it finds different sources rather than the same one twice. The non-Gaussianity measure is built from tanh, for reasons that are in the paper and do not matter here.

def whiten(X, k):
    """PCA-whiten: rotate to the top k components and scale each to unit variance."""
    Xc = X - X.mean(axis=1, keepdims=True)
    cov = Xc @ Xc.T / Xc.shape[1]
    variance, vectors = np.linalg.eigh(cov)
    order = np.argsort(variance)[::-1][:k]
    variance, vectors = variance[order], vectors[:, order]
    W = (vectors / np.sqrt(variance)).T                    # k x channels
    unwhiten = vectors * np.sqrt(variance)                 # channels x k: undoes W
    return W @ Xc, unwhiten

def fastica(Z, seed=0, steps=200):
    """Find a rotation of whitened data Z that makes the rows as non-Gaussian as possible."""
    k = Z.shape[0]
    W = np.random.default_rng(seed).standard_normal((k, k))
    for _ in range(steps):
        Y = W @ Z
        g, g_prime = np.tanh(Y), 1 - np.tanh(Y) ** 2       # the non-linearity and its derivative
        W_new = g @ Z.T / Z.shape[1] - g_prime.mean(axis=1)[:, None] * W
        u, _, vt = np.linalg.svd(W_new)
        W_new = u @ vt                                     # keep the rows orthogonal (decorrelated)
        if np.max(np.abs(np.abs((W_new * W).sum(axis=1)) - 1)) < 1e-6:
            break                                          # nothing changed: converged
        W = W_new
    return W

Z, unwhiten = whiten(X, 2)
W = fastica(Z)
recovered = W @ Z
for i in range(2):
    print(f"recovered {i}: correlation with s1 {abs(np.corrcoef(recovered[i], s1)[0, 1]):.2f}, with s2 {abs(np.corrcoef(recovered[i], s2)[0, 1]):.2f}")
# recovered 0: correlation with s1 0.99, with s2 0.00
# recovered 1: correlation with s1 0.03, with s2 1.00

That last sentence is the first thing to remember about ICA in practice: the components come out in no particular order, scaled arbitrarily, and possibly upside down. Nothing in the data says which source is “first” or what its units are. You identify them afterwards by what they look like, which is what the EEG example is about.

A real cocktail party: 64 channels of EEG

The file holds one minute of resting EEG from one person, 64 electrodes at 160 samples per second, from the PhysioNet motor-imagery dataset. Plot four channels. Fp1 sits on the forehead above the left eye.

eeg = np.load("eeg_s001r01.npz")
E = eeg["eeg"].astype(float)                               # 64 channels x 9760 samples, microvolts
fs = float(eeg["fs"])                                      # 160 samples per second
names = list(eeg["channels"])
E = E - E.mean(axis=1, keepdims=True)
te = np.arange(E.shape[1]) / fs
print(E.shape, "channels x samples,", E.shape[1] / fs, "seconds")
# (64, 9760) channels x samples, 61.0 seconds

fig, axes = plt.subplots(4, 1, figsize=(10, 6), sharex=True)
for ax, ch in zip(axes, ["Fp1", "Fz", "Cz", "Oz"]):
    ax.plot(te, E[names.index(ch)], lw=0.7)
    ax.set(ylabel=ch, xlim=(0, 20))
axes[-1].set_xlabel("time (s)")
plt.show()

Blinks are the classic ICA target because they satisfy every assumption: one source, fixed location, mixed into every channel in fixed proportions, and about as non-Gaussian as a signal gets. Run the same two functions on 64 channels, keeping 20 components, and ask which component looks most like the forehead.

Z, unwhiten = whiten(E, 20)                                # keep the 20 strongest components of 64 channels
W = fastica(Z)
sources = W @ Z                                            # 20 independent components x samples
patterns = unwhiten @ W.T                                  # 64 x 20: how much of each component each channel sees

front = (E[names.index("Fp1")] + E[names.index("Fp2")]) / 2
similarity = [abs(np.corrcoef(s, front)[0, 1]) for s in sources]
blink = int(np.argmax(similarity))
print(f"component {blink} looks most like the forehead channels: correlation {similarity[blink]:.2f}")
strongest = np.argsort(-np.abs(patterns[:, blink]))[:6]
print("channels where it is strongest:", [names[i] for i in strongest])
# component 18 looks most like the forehead channels: correlation 0.93
# channels where it is strongest: ['Fp1', 'Fp2', 'Af8', 'Fpz', 'Af7', 'Af3']

The component the algorithm found correlates 0.93 with the forehead channels, and its scalp pattern is a textbook blink: largest at the three electrodes nearest the eyes, falling off smoothly backwards. Nothing told the algorithm where the eyes are. The pattern comes from the mixing proportions it discovered, which are the physics of a dipole under the forehead.

Removing it

Because ICA gives you both the sources and the mixing pattern, you can set one source to zero and rebuild the channels without it. That is artifact removal, and it is what every EEG pipeline does with blinks, heartbeats and muscle.

sources_clean = sources.copy()
sources_clean[blink] = 0                                   # silence the blink component
E_clean = patterns @ sources_clean                         # and rebuild the 64 channels without it

for ch in ["Fp1", "Fz", "Cz"]:
    i = names.index(ch)
    print(f"{ch}: standard deviation {E[i].std():5.1f} uV before, {E_clean[i].std():5.1f} uV after")
# Fp1: standard deviation 110.3 uV before,  46.1 uV after
# Fz: standard deviation  60.8 uV before,  41.6 uV after
# Cz: standard deviation  54.1 uV before,  47.7 uV after

fig, axes = plt.subplots(3, 1, figsize=(10, 5.5), sharex=True)
axes[0].plot(te, sources[blink], lw=0.8)
axes[0].set(ylabel="blink component")
axes[1].plot(te, E[names.index("Fp1")], lw=0.7, label="Fp1 before")
axes[1].plot(te, E_clean[names.index("Fp1")], lw=0.7, label="Fp1 after")
axes[1].legend()
axes[2].plot(te, E[names.index("Cz")], lw=0.7, label="Cz before")
axes[2].plot(te, E_clean[names.index("Cz")], lw=0.7, label="Cz after")
axes[2].legend()
axes[2].set(xlabel="time (s)", xlim=(0, 20))
plt.show()

Fp1’s standard deviation drops from 110 to 46 microvolts, which is the blinks leaving. Cz barely changes, because little of the blink reached it. What remains on Fp1 is the EEG that was always there under the blinks, and that is what you analyse.

What ICA assumes, and when it fails

• Sources are independent. Not just uncorrelated. Two brain regions that co-activate are not independent and ICA may split them strangely.
• At most one source is Gaussian. The method finds non-Gaussianity; two Gaussian sources cannot be told apart by any rotation.
• Mixing is linear and fixed. True for the volume conduction of EEG and MEG, which is why ICA works so well there. Not true if the source moves, which is why a subject who shifts in the chair halfway through gives two blink components.
• You need enough data. The usual rule of thumb is at least 20 to 30 times the square of the number of components. One minute at 160 Hz, 9,760 samples, is enough for 20 components (20 x 400 = 8,000) and nowhere near enough for 64 (which would want 80,000 or more).
• It does not tell you what the components are. Order, sign and scale are arbitrary, and identifying the blink was our job. Automated labelling exists (ICLabel is the common one) but it is a classifier trained on human labels, not part of ICA.

Exercises

1. Run fastica with a different seed. Which components change order? Which change sign? Do any change shape?
2. Keep 10 components instead of 20, then 40. At what point does the blink component stop being clean?
3. Find the component with the highest kurtosis (((s - s.mean())**4).mean() / s.var()**2 - 3) instead of the highest correlation with Fp1. Is it the blink? If not, plot it and its scalp pattern and work out what it is.
4. Harder: the recording also contains a heartbeat somewhere, a sharp spike about once a second with a pattern that is strongest at the back and sides of the head. Find it.

Next primer: Bayesian inference. Data: EEG Motor Movement/Imagery Dataset (Schalk et al. 2004) via PhysioNet (Goldberger et al. 2000), subject 1, run 1, Open Data Commons Attribution licence. FastICA: Hyvärinen and Oja, 2000.

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.

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()

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()

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.