import plotly.io as pio; pio.renderers.default = "notebook_connected"
Lecture 04: Probability and Density Estimation – CS 189, Fall 2026
This is based on a lecture notebook created by Prof. Gonzalez but then rewritten by Claude to add greater clarity around steps.
Accompanying demonstration: the accuracy of a wake-word detector.
This notebook follows the lecture. It fixes three plausible numbers for a wake-word detector, computes the probability that a positive detection is correct, converts that probability into a daily cost, and then estimates the three numbers from data using maximum likelihood.
Notation follows the slides. For each segment of sound there are three random variables:
| $Z \in \{0, 1\}$ | a person says a wake word |
| $\mathbf{x} \in \mathbb{R}^D$ | the sound collected by the microphone |
| $Y \in \{0, 1\}$ | the detector reports a wake word |
Sections 1 through 6 concern $Z$ and $Y$. Section 7 onward models a feature computed from $\mathbf{x}$.
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
rng = np.random.default_rng(189)
# Berkeley palette, so the notebook matches the slides.
BLUE, GOLD, GREEN, ORANGE = '#002675', '#FDB515', '#028842', '#DF5327'
# Plotly abbreviates small numbers as "10u", which collides with the parameter mu.
# These give the log axes explicit powers of ten instead.
DEC = [1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
DECTXT = ['10^-6', '10^-5', '10^-4', '10^-3', '10^-2']
1. Three numbers describe the detector¶
A detector is characterized by two conditional probabilities, the detection rate and the false alarm rate. A third number, the rate at which wake words are actually said, is a property of the environment rather than of the detector.
prior = 0.0001 # p(Z = 1) wake words are rare: 1 segment in 10,000
recall = 0.99 # p(Y = 1 | Z = 1) detection rate, or true positive rate
fpr = 0.001 # p(Y = 1 | Z = 0) false alarm rate
prior, recall, fpr
(0.0001, 0.99, 0.001)
2. Simulating one day¶
A device that scores one segment per second faces $60 \times 60 \times 24 = 86{,}400$ segments in a day. Each segment is a Bernoulli trial.
The sampling procedure below draws $Z$ first and then draws $Y$ from the conditional distribution that matches each value of $Z$. This is the product rule, $p(z, y) = p(y \mid z)\, p(z)$, used as a recipe for generating samples.
N = 86_400
Z = rng.random(N) < prior # was a wake word said?
p_fire = np.where(Z, recall, fpr) # p(Y = 1 | Z), elementwise
Y = rng.random(N) < p_fire # did the detector fire?
Z.sum(), Y.sum()
(np.int64(7), np.int64(73))
Counting the four combinations gives the empirical joint distribution over $(Z, Y)$, that is $\hat{p}(Z = z_i, Y = y_j) = n_{ij}/N$ computed from samples rather than assumed.
counts = pd.crosstab(pd.Series(Z.astype(int), name='Z (said)'),
pd.Series(Y.astype(int), name='Y (detected)'))
counts
| Y (detected) | 0 | 1 |
|---|---|---|
| Z (said) | ||
| 0 | 86327 | 66 |
| 1 | 0 | 7 |
joint = counts / N # the empirical joint distribution
joint
| Y (detected) | 0 | 1 |
|---|---|---|
| Z (said) | ||
| 0 | 0.999155 | 0.000764 |
| 1 | 0.000000 | 0.000081 |
The table satisfies the two properties from lecture. Its four entries sum to one (normalization), and summing across a row gives $p(Z = z_i)$ (the sum rule).
print('sums to 1 :', joint.values.sum())
print('p(Z = 1) marginal:', joint.sum(axis=1)[1], ' vs assumed', prior)
print('p(Y = 1) marginal:', joint.sum(axis=0)[1])
sums to 1 : 1.0 p(Z = 1) marginal: 8.101851851851852e-05 vs assumed 0.0001 p(Y = 1) marginal: 0.0008449074074074075
3. Precision¶
The detection rate and the false alarm rate both condition on $Z$. They describe the behaviour of the detector given the state of the world. The device faces the reverse problem: it has observed $Y = 1$ and requires
$$p(Z = 1 \mid Y = 1)$$
which is a different conditional. Reversing the direction of a conditional is what Bayes' theorem does.
The empirical value is computed first: among the segments where the detector fired, the fraction in which a wake word was actually said.
precision_empirical = Z[Y].mean() # among fired segments, fraction that were real
print(f'fired {Y.sum():,} times today, {Z[Y].sum():,} were real')
print(f'empirical precision: {precision_empirical:.3f}')
fired 73 times today, 7 were real empirical precision: 0.096
# Bayes' theorem, with the denominator expanded by the sum rule.
def precision(prior, recall, fpr):
p_fire_and_real = recall * prior
p_fire = recall * prior + fpr * (1 - prior)
return p_fire_and_real / p_fire
print(f'Bayes: {precision(prior, recall, fpr):.4f}')
print(f'simulated: {precision_empirical:.4f}')
Bayes: 0.0901 simulated: 0.0959
The precision is approximately 9%. A detector that catches 99% of wake words and rarely fires spuriously is nevertheless wrong more than nine times out of ten when it does fire.
The cause is the base rate. Silent segments outnumber wake-word segments by a factor of 10,000, so a 0.1% error rate on the silent segments produces more false alarms than there are true events. Stated in whole numbers, out of 1,000,000 segments:
M = 1_000_000
said = M * prior
pd.DataFrame({
'segments': [said, M - said],
'detector fires': [said * recall, (M - said) * fpr],
}, index=['wake word said', 'no wake word']).astype(int)
| segments | detector fires | |
|---|---|---|
| wake word said | 100 | 99 |
| no wake word | 999900 | 999 |
99 correct detections against 999 false alarms.
4. Comparing the two design parameters¶
The base rate is fixed by the environment. The detection rate and the false alarm rate are design choices. Sweeping each one shows how much precision responds to each.
fpr_grid = np.logspace(-6, -2, 200)
fig = go.Figure()
for r, dash in [(0.9999, 'solid'), (0.90, 'dash'), (0.50, 'dot')]:
fig.add_trace(go.Scatter(x=fpr_grid, y=precision(prior, r, fpr_grid),
name=f'recall = {r}', line=dict(dash=dash, color=BLUE)))
fig.add_vline(x=fpr, line_color=ORANGE,
annotation_text='our detector', annotation_position='top')
fig.update_xaxes(type='log', title='false alarm rate p(Y=1 | Z=0)',
tickmode='array', tickvals=DEC, ticktext=DECTXT, minor_showgrid=False)
fig.update_yaxes(title='precision p(Z=1 | Y=1)', tickformat='.0%')
fig.update_layout(title='Precision as a function of the false alarm rate',
legend_title_text='recall p(Y=1 | Z=1)',
width=850, height=480)
fig
At the operating point, halving the detection rate from 0.99 to 0.50 reduces precision from 9.0% to 4.8%, a factor of two. Reducing the false alarm rate by one decade, from $10^{-3}$ to $10^{-4}$, raises precision from 9.0% to 49.8%. A second decade raises it to 90.8%.
Precision therefore responds roughly linearly to the detection rate and by decades to the false alarm rate. Reducing false alarms is the more effective design change.
5. Expressing precision as a daily cost¶
Each time the detector fires, the device invokes the full speech model, which consumes energy or incurs a cloud charge. Let $C$ denote the cost of one invocation and let $I_t$ be the indicator that segment $t$ is a false alarm. The expected cost over a day is
$$\mathbb{E}\!\left[\sum_{t=1}^{86{,}400} C\, I_t\right] = C \sum_{t=1}^{86{,}400} \mathbb{E}[I_t] = C \sum_{t=1}^{86{,}400} p(Z = 0, Y = 1)$$
The exchange of expectation and summation is linearity of expectation, which requires no independence assumption. Consecutive seconds of audio are strongly dependent, and the result holds regardless.
# Expected number of times the detector fires in a day.
def wakeups_per_day(prior, recall, fpr, seconds=86_400):
return seconds * (recall * prior + fpr * (1 - prior))
false_alarms = 86_400 * fpr * (1 - prior)
total = wakeups_per_day(prior, recall, fpr)
print(f'expected false alarms per day: {false_alarms:.1f}')
print(f'expected wake-ups in total: {total:.1f}')
print(f'expected wanted wake-ups: {total - false_alarms:.1f}')
print(f'observed wake-ups today: {Y.sum()}')
expected false alarms per day: 86.4 expected wake-ups in total: 94.9 expected wanted wake-ups: 8.6 observed wake-ups today: 73
The device performs approximately 86.4 unnecessary invocations per day against approximately 8.6 useful ones. This restates the 9% precision in units that can be compared against a battery budget.
Sweeping the false alarm rate again, now in units of invocations per day rather than probability:
total_grid = wakeups_per_day(prior, recall, fpr_grid)
useful_grid = total_grid * precision(prior, recall, fpr_grid)
fig = go.Figure()
fig.add_trace(go.Scatter(x=fpr_grid, y=total_grid, name='all wake-ups',
line=dict(color=ORANGE)))
fig.add_trace(go.Scatter(x=fpr_grid, y=useful_grid, name='wanted wake-ups',
line=dict(color=GREEN)))
fig.add_vline(x=fpr, line_color=BLUE, annotation_text='our detector')
fig.update_xaxes(type='log', title='false alarm rate p(Y=1 | Z=0)',
tickmode='array', tickvals=DEC, ticktext=DECTXT, minor_showgrid=False)
fig.update_yaxes(type='log', title='expected wake-ups per day',
tickmode='array', tickvals=[10, 100, 1000],
ticktext=['10', '100', '1000'], minor_showgrid=False)
fig.update_layout(title='Expected wake-ups per day over 86,400 segments',
width=850, height=480, template='plotly_white')
fig
The lower curve is constant. The expected number of wanted wake-ups is $86{,}400 \times p(Z{=}1) \times \text{recall} \approx 8.6$ per day, an expression in which the false alarm rate does not appear. The vertical distance between the curves is the wasted cost, and reducing the false alarm rate collapses the upper curve onto the lower one: 94.9 wake-ups per day at $10^{-3}$, 17.2 at $10^{-4}$, and 9.4 at $10^{-5}$.
6. Estimating $p(Z = 1)$¶
The value $p(Z = 1) = 0.0001$ was assumed. Every quantity computed above depends on it, so it must itself be estimated from data.
Each segment is a Bernoulli trial, $Z_n \sim \mathrm{Bern}(\mu)$ with $\mu = p(Z = 1)$. The simulated day provides a labelled sample.
n1 = int(Z.sum()) # segments containing a wake word
n0 = int((~Z).sum()) # segments without
n1, n0
(7, 86393)
The likelihood function¶
For a candidate value of $\mu$, the probability of the observed data is
$$p(\mathcal{D} \mid \mu) = \prod_{n=1}^{N} \mu^{z_n}(1-\mu)^{1-z_n} = \mu^{n_1}(1-\mu)^{n_0}$$
This is a function of $\mu$. The data is fixed; it has already been observed. Evaluating the function on a grid of candidate values gives:
mu_grid = np.linspace(1e-6, 5e-4, 400)
likelihood = mu_grid ** n1 * (1 - mu_grid) ** n0
likelihood.max()
np.float64(2.090026321235871e-32)
The largest value on the grid is of order $10^{-32}$. The likelihood is a product of 86,400 factors, each less than one, so it decays rapidly toward zero. One day of audio is not quite enough to exhaust the floating point range. Thirty days is.
# The same calculation for 30 days of recording.
likelihood_month = mu_grid ** (30 * n1) * (1 - mu_grid) ** (30 * n0)
print('largest value on the grid:', likelihood_month.max())
print('grid points that are exactly 0:', (likelihood_month == 0).sum(), 'out of', len(mu_grid))
largest value on the grid: 0.0 grid points that are exactly 0: 400 out of 400
Every grid point underflows to zero. The function to be maximized has become identically zero in floating point, and no maximum can be located. This is the practical reason for working with the log likelihood,
$$\ln p(\mathcal{D} \mid \mu) = \sum_{n=1}^{N} \ln p(z_n \mid \mu) = n_1 \ln \mu + n_0 \ln(1-\mu)$$
A sum of 86,400 moderate negative numbers is representable. Since $\ln$ is monotonically increasing, the maximum occurs at the same value of $\mu$.
def log_likelihood(mu, n1, n0):
return n1 * np.log(mu) + n0 * np.log1p(-mu)
ll = log_likelihood(mu_grid, n1, n0)
mu_hat_grid = mu_grid[ll.argmax()] # best value on the grid
mu_hat_closed = n1 / (n0 + n1) # the closed form derived in lecture
print(f'grid search : {mu_hat_grid:.6f}')
print(f'n1 / N : {mu_hat_closed:.6f}')
print(f'true value : {prior:.6f}')
grid search : 0.000081 n1 / N : 0.000081 true value : 0.000100
fig = go.Figure()
fig.add_trace(go.Scatter(x=mu_grid, y=ll, line=dict(color=BLUE), name='log likelihood'))
fig.add_vline(x=mu_hat_closed, line_color=GOLD, line_width=3,
annotation_text=' n1/N', annotation_position='top left')
fig.add_vline(x=prior, line_color=GREEN, line_dash='dash',
annotation_text='true mu ', annotation_position='top right')
fig.update_xaxes(title='mu (candidate values of p(Z = 1))',
tickmode='array', tickvals=[0, 1e-4, 2e-4, 3e-4, 4e-4, 5e-4],
ticktext=['0', '0.0001', '0.0002', '0.0003', '0.0004', '0.0005'])
fig.update_yaxes(title='ln p(D | mu)', range=[ll.max() - 40, ll.max() + 3])
fig.update_layout(title='Log likelihood of one day of labelled audio',
width=850, height=480,showlegend=True)
fig
The function has a single maximum, and the closed form $\mu_{ML} = n_1/N$ coincides with it. Maximum likelihood estimation consists of these two steps: write the log likelihood, then maximize it. Here the maximum is available analytically. In later lectures it will not be, and the maximum will be located by gradient ascent instead.
The estimate $\mu_{ML} = n_1/N$ is the observed frequency, which is the empirical distribution $\hat{p}$ defined earlier in the lecture. Maximum likelihood establishes that counting and dividing is the optimal estimate under a Bernoulli model, not merely a convenient one.
Sample size and rare events¶
The estimate above is usable because 86,400 segments happened to contain several wake words. Shorter recordings behave differently. The experiment below is repeated 2,000 times at two sample sizes.
def mle_replicates(N, reps=2000, mu=prior):
n1 = rng.binomial(N, mu, size=reps) # counts, without materialising every segment
return n1 / N
small = mle_replicates(10_000)
large = mle_replicates(1_000_000)
print(f'N = 10,000 : estimate is exactly 0 in {100 * (small == 0).mean():.0f}% of runs')
print(f'N = 1,000,000 : estimate is exactly 0 in {100 * (large == 0).mean():.0f}% of runs')
N = 10,000 : estimate is exactly 0 in 37% of runs N = 1,000,000 : estimate is exactly 0 in 0% of runs
fig = go.Figure()
fig.add_trace(go.Histogram(x=small, name='N = 10,000', marker_color=ORANGE, opacity=0.75, nbinsx=50))
fig.add_trace(go.Histogram(x=large, name='N = 1,000,000', marker_color=BLUE, opacity=0.75))
fig.add_vline(x=prior, line_color=GREEN, line_dash='dash', annotation_text='true mu')
fig.update_xaxes(title='mu_ML', range=[0, 4e-4])
fig.update_yaxes(title='count (out of 2,000 runs)')
fig.update_layout(barmode='overlay', title='Sampling distribution of the Bernoulli MLE',
width=850, height=460)
fig
With ten thousand segments the estimate is frequently exactly zero, and a zero estimate makes the precision calculation of section 3 undefined. The maximum likelihood estimate is correct with respect to the observed data and uninformative about an event that never occurred in it.
The required sample size scales with the rarity of the event rather than with the total quantity of data. An alternative is to introduce a prior distribution over $\mu$ and maximize the posterior, which is MAP estimation, covered later in the course.
7. Maximum a posteriori estimation¶
Section 6 ended with a defect. When the event is rare and the sample is small, the maximum likelihood estimate is frequently exactly zero, and a zero estimate makes the precision calculation of section 3 undefined.
Maximum a posteriori estimation repairs this by treating $\mu$ as a random variable and supplying a prior distribution over it. The estimate maximizes the posterior rather than the likelihood:
$$\mu_{MAP} = \arg\max_\mu \; p(\mathcal{D} \mid \mu)\, p(\mu)$$
# How often does maximum likelihood return exactly zero?
counts_small = rng.binomial(10_000, prior, size=2000)
print(f'N = 10,000: mu_ML is exactly 0 in {100 * (counts_small == 0).mean():.0f}% of runs')
N = 10,000: mu_ML is exactly 0 in 37% of runs
The Beta prior¶
A prior over $\mu$ must live on $[0, 1]$. The Beta distribution does, and it has the same algebraic form as the Bernoulli likelihood:
$$p(\mu \mid a, b) \propto \mu^{a-1}(1-\mu)^{b-1} \qquad p(\mathcal{D} \mid \mu) \propto \mu^{n_1}(1-\mu)^{n_0}$$
Because the two match, the log posterior is the log likelihood with $n_1$ replaced by $n_1 + a - 1$ and $n_0$ by $n_0 + b - 1$. Differentiating and solving, exactly as before, gives
$$\mu_{MAP} = \frac{n_1 + a - 1}{N + a + b - 2}$$
The prior therefore acts as $a - 1$ imaginary successes and $b - 1$ imaginary failures added to the observed counts.
def map_estimate(n1, N, a=2.0, b=2.0):
return (n1 + a - 1) / (N + a + b - 2)
def ml_estimate(n1, N):
return n1 / N
# The failing case from section 6: no wake words observed in 10,000 segments.
print('mu_ML =', ml_estimate(0, 10_000))
print('mu_MAP =', map_estimate(0, 10_000), ' (a = b = 2)')
print('true =', prior)
mu_ML = 0.0 mu_MAP = 9.998000399920016e-05 (a = b = 2) true = 0.0001
With $a = b = 2$ the formula is $(n_1 + 1)/(N + 2)$, one imaginary observation of each outcome. The estimate can no longer be exactly 0 or exactly 1, so the precision calculation is defined again. This is Laplace smoothing, the same device used in naive Bayes classifiers and n-gram language models.
The closeness of this particular estimate to the true rate is a coincidence of the numbers. What the prior guarantees is that the estimate is strictly inside $(0, 1)$, not that it is accurate.
Watching the prior lose influence¶
The prior contributes a fixed number of pseudo-counts while the data contributes $N$ real ones, so the two estimates should converge as $N$ grows.
Ns = np.array([100, 300, 1_000, 3_000, 10_000, 30_000, 100_000, 300_000, 1_000_000])
reps = 400
ml, mp = [], []
for N in Ns:
n1 = rng.binomial(N, prior, size=reps)
ml.append(ml_estimate(n1, N).mean())
mp.append(map_estimate(n1, N).mean())
fig = go.Figure()
fig.add_trace(go.Scatter(x=Ns, y=ml, mode='lines+markers', name='mu_ML',
line=dict(color=ORANGE)))
fig.add_trace(go.Scatter(x=Ns, y=mp, mode='lines+markers', name='mu_MAP (a = b = 2)',
line=dict(color=GREEN)))
fig.add_hline(y=prior, line_dash='dash', line_color=BLUE,
annotation_text='true mu = 0.0001')
fig.update_xaxes(type='log', title='N (segments observed)')
fig.update_yaxes(type='log', title='average estimate over 400 runs')
fig.update_layout(title='The prior matters only while the data is scarce',
width=850, height=470)
fig
At small $N$ the MAP estimate sits well above the true rate, pulled toward the prior mean of $a/(a+b) = 0.5$. The maximum likelihood curve is closer on average but is the one that returns exactly zero, which the logarithmic axis cannot even display. By $N = 10^6$ the two agree.
The prior is a modeling assumption like any other. A badly chosen one biases the estimate, and enough data will eventually overrule it.
The posterior, not just its maximum¶
MAP reports the single most probable parameter value. The posterior itself carries more information, namely how sharply that value is determined.
from scipy.stats import beta as beta_dist
grid = np.linspace(0, 6e-4, 500)
fig = go.Figure()
for N, color in [(1_000, GOLD), (10_000, ORANGE), (100_000, GREEN), (1_000_000, BLUE)]:
n1 = rng.binomial(N, prior)
post = beta_dist(n1 + 2, N - n1 + 2) # Beta(n1 + a, n0 + b)
fig.add_trace(go.Scatter(x=grid, y=post.pdf(grid), name=f'N = {N:,} (n1 = {n1})',
line=dict(color=color)))
fig.add_vline(x=prior, line_dash='dash', annotation_text='true mu')
fig.update_xaxes(title='mu')
fig.update_yaxes(title='posterior density p(mu | D)')
fig.update_layout(title='The posterior concentrates as data accumulates',
width=850, height=470)
fig
Each curve is the posterior after a different amount of data. With a thousand segments the posterior is broad and the data has barely constrained $\mu$; with a million it is a narrow spike. The MAP estimate is the location of each peak, and it discards the width.
Reporting the whole posterior rather than its maximum is Bayesian inference proper. MAP is the point estimate that keeps the prior while remaining as cheap to compute as maximum likelihood.
The same term reappears in Lecture 7 under a different name. Written as a minimization, the MAP objective is $-\ln p(\mathcal{D} \mid w) - \ln p(w)$, which is a loss plus a penalty. The prior is the regularizer.
Summary¶
| Bayes' theorem | converts $p(Y \mid Z)$, which the detector provides, into $p(Z \mid Y)$, which the device requires. The base rate dominates the result. |
| Expectation | converts a precision of 9% into roughly 86.4 unnecessary invocations per day, using linearity and no independence assumption. |
| Likelihood | $p(\mathcal{D} \mid w)$ regarded as a function of $w$; logarithms make it numerically tractable. |
| Maximum likelihood | for a Bernoulli, the observed frequency $n_1/N$. |
| Maximum a posteriori | the same estimate with $a-1$ and $b-1$ pseudo-counts from a Beta prior, which repairs the rare-event case. |
Lecture 5 applies both estimators to continuous data, starting with the Gaussian.