---
title: $N$-gram models
bibliography: ../references.bib
jupyter: python3
---
::: {.callout-note}
The executable cells on this page use Python 3.14 syntax, matching the project configuration. The Quarto kernel is named `python3`; before rendering with execution, configure that kernelspec to point to a Python 3.14 environment.
:::
In the previous section, I defined a language model as a specification of a probability measure on $\Sigma^*$. For a distribution over variable-length strings, the **termination-inclusive chain rule** (TCR) gives the exact decomposition
$$
p(w_1 w_2 \ldots w_n)
=\prod_{i=1}^{n+1}p(w_i\mid w_1\ldots w_{i-1}),
\qquad w_{n+1}=\langle\text{/s}\rangle.
$$
The final factor is the probability of stopping after $w_n$. Omitting it would instead give the probability of the observed symbols without asserting termination—for instance, a quantity conditioned on a separately specified length.
The trouble with using this decomposition directly is that each conditional $p(w_i \mid w_1 \ldots w_{i-1})$ depends on the entire prefix $w_1 \ldots w_{i-1}$. As the prefix gets longer, the number of possible prefixes grows exponentially: there are $|\Sigma|^{i-1}$ possible prefixes of length $i - 1$. Estimating a separate distribution for every possible prefix would require an enormous amount of data and memory, and for most prefixes we'd never observe enough examples to get a reliable estimate.
An $N$-gram model deals with this by making an independence assumption: the distribution of $w_i$ depends only on the preceding $N - 1$ symbols, not on the entire prefix. That is:
$$p(w_i \mid w_1 \ldots w_{i-1}) \approx p(w_i \mid w_{i-N+1} \ldots w_{i-1})$$
Substituting this approximation into the exact chain rule decomposition gives us:
$$p(w_1 w_2 \ldots w_n) \approx \prod_{i=1}^{n+1} p(w_i \mid w_{i-N+1} \ldots w_{i-1})$$
where we pad the beginning of the string with $N-1$ copies of a start symbol $\langle\text{s}\rangle$ and append an end symbol $\langle\text{/s}\rangle$ at position $n + 1$.^[The start and end symbols ensure that the model can learn which phones tend to begin and end words. Without them, the model wouldn't distinguish between a phone that appears word-initially and one that appears word-medially.]
Each conditional distribution $p(w_i \mid w_{i-N+1} \ldots w_{i-1})$ is a categorical distribution over $\Sigma$ (plus the end symbol), parameterized by a vector $\boldsymbol{\theta}^{(w_{i-N+1} \ldots w_{i-1})}$. So the full model is parameterized by a collection $\boldsymbol{\Theta}$ of such vectors—one for each possible context of length $N - 1$.
## Maximum likelihood estimation
Given a corpus of strings, we want to find the parameters $\boldsymbol{\Theta}$ that maximize the likelihood of the observed data. Let $\mathcal A=\Sigma\cup\{\langle\text{/s}\rangle\}$ be the continuation alphabet and let $K=|\mathcal A|$. For a unigram model ($N = 1$), where there is no conditioning context, this amounts to:
$$\hat{\boldsymbol{\Theta}}_\text{MLE} = \arg_{\boldsymbol{\Theta}}\max \prod_{j=1}^{K} \theta_j^{c_j}$$
where $c_j$ counts continuation symbol $j$, including one end symbol for every corpus string. The solution is to set each $\theta_j$ proportional to its count:
$$\hat{\theta}_j = \frac{c_j}{\sum_k c_k}$$
But why do the relative counts maximize the likelihood? We can show this with a **simplex-constrained likelihood argument** (SLA). Let $C=\sum_j c_j>0$. Because the logarithm is increasing, the parameter vector that maximizes the likelihood also maximizes the log-likelihood
$$
\ell(\boldsymbol\theta)=\sum_{j=1}^K c_j\log\theta_j
$$
subject to $\sum_j\theta_j=1$ and $\theta_j\geq0$. First suppose $c_j>0$ for every $j$, so the optimum is in the interior of the simplex. Introduce a Lagrange multiplier $\mu$ for the normalization constraint:
$$
\mathcal{L}(\boldsymbol\theta,\mu)
=\sum_{j=1}^K c_j\log\theta_j
+\mu\left(1-\sum_{j=1}^K\theta_j\right).
$$
At a stationary point, differentiating with respect to an arbitrary $\theta_j$ gives
$$
\frac{\partial\mathcal{L}}{\partial\theta_j}
=\frac{c_j}{\theta_j}-\mu=0,
$$
so $\theta_j=c_j/\mu$. Summing this identity over $j$ and imposing $\sum_j\theta_j=1$ yields $1=C/\mu$, hence $\mu=C$. Substitution gives $\hat\theta_j=c_j/C$, the claimed relative-frequency estimate.
The stationary point is a global maximum, not merely a local candidate. For positive counts, the Hessian of the log-likelihood is diagonal:
$$
\frac{\partial^2\ell}{\partial\theta_j^2}
=-\frac{c_j}{\theta_j^2}<0,
\qquad
\frac{\partial^2\ell}{\partial\theta_j\partial\theta_k}=0
\quad(j\neq k).
$$
Thus, $\ell$ is strictly concave on the interior of the simplex. For a concave objective subject to an affine equality constraint, the first-order Lagrange conditions are sufficient for global maximality; strict concavity makes that maximizer unique. Hence, the stationary point is the unique global maximum when every count is positive.
What if some counts are zero? Let $P=\{j\mid c_j>0\}$. Any vector with $\theta_j=0$ for some $j\in P$ has likelihood zero and cannot beat a vector that assigns every member of $P$ positive mass. Now consider a feasible vector with finite log-likelihood that assigns total mass $\delta=\sum_{j\notin P}\theta_j>0$ to zero-count symbols. Renormalize the positive-count coordinates by setting $\theta'_j=\theta_j/(1-\delta)$ for $j\in P$ and $\theta'_j=0$ otherwise. The change in log-likelihood is
$$
\begin{aligned}
\ell(\boldsymbol\theta')-\ell(\boldsymbol\theta)
&=\sum_{j\in P}c_j
\log\left(\frac{\theta_j}{1-\delta}\right)
-\sum_{j\in P}c_j\log\theta_j\\
&=-C\log(1-\delta)>0.
\end{aligned}
$$
Hence no maximizing vector can assign positive mass to a zero-count symbol. On the simplex face supported by $P$, every remaining count is positive, so the strict-concavity argument applies and gives the unique maximizer $\hat\theta_j=c_j/C$ for $j\in P$, with $\hat\theta_j=0$ elsewhere. If $C=0$, every likelihood factor is $\theta_j^0=1$; all parameter vectors have the same likelihood, and the MLE is not identified. Here and below, the likelihood uses the convention $\theta^0=1$, including $0^0=1$. The corresponding log-likelihood convention is $0\log\theta=0$; if $c_j>0$ and $\theta_j=0$, however, that term is $-\infty$.
More generally, for any $N$-gram model and any context with a positive continuation count:
$$\hat{\theta}_{w_{i_1} \ldots w_{i_N}} = \frac{c_{i_1 i_2 \ldots i_N}}{\sum_j c_{i_1 i_2 \ldots j}}$$
where $c_{i_1 i_2 \ldots i_N}$ is the number of times the sequence $w_{i_1} \ldots w_{i_N}$ occurs in the corpus. In words: the estimated probability of a continuation symbol, which may be $\langle\text{/s}\rangle$, given a context is the proportion of times that symbol follows the context in the training data. If the denominator is zero, the corpus supplies no likelihood information about that conditional vector, so its unsmoothed MLE is not identified.
To see why the unigram calculation applies, fix an arbitrary context $h=w_{i_1}\ldots w_{i_{N-1}}$. The corpus likelihood groups into one factor $\prod_x(\theta_x^{(h)})^{c(h,x)}$ for each context, and the parameters for that factor satisfy their own constraint $\sum_x\theta_x^{(h)}=1$. Write $C(h)=\sum_{x'}c(h,x')$. When $C(h)>0$, applying the SLA to this factor gives
$$
\hat\theta_x^{(h)}
=\frac{c(h,x)}{\sum_{x'}c(h,x')}.
$$
Because the factors use disjoint parameter vectors, maximizing each identified factor separately maximizes their product. We thus obtain the general $N$-gram estimate for observed contexts while isolating the unidentified $C(h)=0$ case.
## Smoothing
One issue with the MLE is that any $N$-gram that never occurs in the training data will be assigned probability zero. This is a problem because when we compute the probability of a string, we multiply conditionals together—and a single zero sends the entire product to zero. In log space, this looks like a single $-\infty$ sending the sum to $-\infty$.
To see why this matters concretely, consider a bigram model trained on a typical English lexicon. The sequence $\text{ŋk}$ probably appears many times (as in *think*, *bank*), so the model learns a reasonable estimate of $p(\text{k} \mid \text{ŋ})$. But the sequence $\text{ŋm}$ may never appear in the training data at all, so the MLE assigns $p(\text{m} \mid \text{ŋ}) = 0$. Now suppose we want to evaluate a nonword that happens to contain $\text{ŋm}$. Under the MLE, the entire string gets probability zero—not because it's maximally bad, but because the model has never seen that one transition. A single gap in the training data wipes out all the information from the other transitions in the string.
### Add-$\lambda$ smoothing
The fix is to ensure that no $N$-gram probability is ever exactly zero. [Add-$\lambda$ smoothing](https://en.wikipedia.org/wiki/Additive_smoothing) (also called Laplace smoothing when $\lambda = 1$) does this by adding a pseudocount $\lambda > 0$ to every $N$-gram count before normalizing:
$$\hat{\theta}_{w_{i_1} \ldots w_{i_N}} = \frac{c_{i_1 i_2 \ldots i_N} + \lambda}{\sum_j \left[c_{i_1 i_2 \ldots j} + \lambda\right]}$$
The denominator does the required normalization. Fix a context $h$, let $K$ be the number of possible continuation symbols, and write $C(h)=\sum_x c(h,x)$. Then
$$
\sum_x\hat\theta_x^{(h)}
=\sum_x\frac{c(h,x)+\lambda}{C(h)+K\lambda}
=\frac{C(h)+K\lambda}{C(h)+K\lambda}
=1.
$$
When $\lambda>0$, every numerator and the denominator are positive, including for a context with no observed continuation. Thus, smoothing produces a well-defined categorical distribution for every context represented by the model.
The pseudocount works as if we had observed each possible $N$-gram $\lambda$ additional times on top of whatever the data actually contains. This has two effects. First, no probability is ever zero: even an $N$-gram that never occurred in the training data gets a small positive probability proportional to $\lambda$. Second, the probabilities of frequently observed $N$-grams are slightly reduced to make room, since the denominator grows by $K\lambda$, where $K$ includes every symbol that may follow the context, including the end symbol.
The parameter $\lambda$ controls the strength of this effect, and it is worth understanding its two extremes:
- When $\lambda = 0$, we recover the MLE for every observed context. The model trusts the data completely: observed counts determine probabilities, and an unseen continuation of an observed context gets probability zero. An entirely unseen context has $C(h)=0$, so its unsmoothed conditional distribution is not identified by the data.
- As $\lambda \to \infty$, divide the numerator and denominator by $\lambda$ to obtain $[c(h,x)/\lambda+1]/[C(h)/\lambda+K]$. The count terms tend to zero, so every conditional probability tends to $1/K$. The model ignores the relative counts and approaches a uniform continuation distribution.
Between these extremes, $\lambda$ interpolates between fitting the data closely and hedging toward a uniform distribution. The following plot illustrates this on a toy example: a unigram distribution over 5 phones with counts $(10, 5, 3, 1, 0)$. As $\lambda$ increases, the distribution flattens from the MLE toward uniform.
```{python}
#| code-fold: true
#| code-summary: Effect of smoothing on a toy distribution
import numpy as np
import matplotlib.pyplot as plt
counts = np.array([10, 5, 3, 1, 0])
phones = ['a', 'b', 'c', 'd', 'e']
lambdas = [0, 0.1, 0.5, 1, 5, 50]
fig, axes = plt.subplots(1, len(lambdas), figsize=(16, 2.8), sharey=True)
for ax, lam in zip(axes, lambdas):
smoothed = (counts + lam) / (counts.sum() + lam * len(counts))
ax.bar(phones, smoothed, color='steelblue', edgecolor='white')
ax.set_title(f'λ = {lam}', fontsize=10)
ax.set_ylim(0, 0.6)
ax.tick_params(labelsize=8)
if lam == 0:
ax.set_ylabel('Probability', fontsize=9)
fig.tight_layout()
plt.show()
```
A small $\lambda$ stays close to the observed frequencies while giving unseen $N$-grams just enough probability to avoid zeros. A large $\lambda$ washes out the distinctions the model learned from the data. The right value of $\lambda$ depends on how much data we have and how much we trust the model's specific estimates—which is exactly the kind of question that [cross-validation](#evaluating-a-predictor) can answer.^[For those interested: add-$\lambda$ smoothing corresponds to placing a symmetric [Dirichlet prior](https://en.wikipedia.org/wiki/Dirichlet_distribution) $\text{Dir}(\boldsymbol{1} + \lambda)$ on each conditional distribution and computing a [maximum a posteriori](https://en.wikipedia.org/wiki/Maximum_a_posteriori_estimation) (MAP) estimate rather than an MLE. The pseudocounts are the prior's way of encoding a belief about what the distribution looks like before seeing any data. When $\lambda = 0$—a uniform Dirichlet $\text{Dir}(\boldsymbol{1})$—the MAP estimate reduces to the MLE. See @hastie_elements_2009[Ch. 8] for a broader discussion of regularization, of which smoothing is a special case.]
Add-$\lambda$ smoothing is the simplest member of a large family of smoothing techniques for $N$-gram models, all of which address the same core problem—redistributing probability mass from observed $N$-grams to unobserved ones—but differ in how they decide how much mass to redistribute and where to put it. [Good-Turing smoothing](https://en.wikipedia.org/wiki/Good%E2%80%93Turing_frequency_estimation) [@good_population_1953] estimates the total probability of unseen events from the number of events seen exactly once. [Jelinek-Mercer smoothing](https://en.wikipedia.org/wiki/Jelinek%E2%80%93Mercer_smoothing) [@jelinek_interpolated_1980] interpolates between higher-order and lower-order $N$-gram estimates, so that a trigram model can fall back to a bigram estimate when the trigram count is sparse. [Katz backoff](https://en.wikipedia.org/wiki/Katz%27s_back-off_model) [@katz_estimation_1987] takes a similar approach but uses the lower-order model only when the higher-order count falls below a threshold. [Kneser-Ney smoothing](https://en.wikipedia.org/wiki/Kneser%E2%80%93Ney_smoothing) [@kneser_improved_1995] refines this further by basing the lower-order distribution not on raw counts but on the number of distinct contexts in which a word appears, which gives better estimates for words that are common only in a few fixed phrases. @chen_empirical_1996 provide a thorough empirical comparison of these methods. We use add-$\lambda$ smoothing here because it is easy to implement and sufficient for our purposes, but the more sophisticated methods are standard in large-scale language modeling.
We can thus view an $N$-gram model as having two hyperparameters: the order $N$, which determines how much context the model conditions on, and the smoothing parameter $\lambda$, which determines how much it trusts the observed counts versus a uniform fallback.
## Fitting $N$-gram models to a lexicon
To see how well $N$-gram models predict phonotactic acceptability, we'll fit models to the [CMU Pronouncing Dictionary](https://web.archive.org/web/20210216211828/http://www.speech.cs.cmu.edu/cgi-bin/cmudict) and evaluate them against the acceptability judgments from @daland_explaining_2011.
We begin by loading the dictionary, which uses the [ARPABET](https://en.wikipedia.org/wiki/ARPABET) representation, and converting it to IPA.
```{python}
#| code-fold: true
#| code-summary: Load CMU Pronouncing Dictionary and convert to IPA
import re
import numpy as np
import pandas as pd
import nltk
nltk.download('cmudict', quiet=True)
from nltk.corpus import cmudict
arpabet_to_phoneme: dict[str, str] = {
'AA': 'ɑ', 'AE': 'æ', 'AH': 'ʌ', 'AO': 'ɔ', 'AW': 'aʊ', 'AY': 'aɪ',
'B': 'b', 'CH': 'tʃ', 'D': 'd', 'DH': 'ð', 'EH': 'ɛ', 'ER': 'ɝ',
'EY': 'eɪ', 'F': 'f', 'G': 'g', 'HH': 'h', 'IH': 'ɪ', 'IY': 'i',
'JH': 'dʒ', 'K': 'k', 'L': 'l', 'M': 'm', 'N': 'n', 'NG': 'ŋ',
'OW': 'oʊ', 'OY': 'ɔɪ', 'P': 'p', 'R': 'ɹ', 'S': 's', 'SH': 'ʃ',
'T': 't', 'TH': 'θ', 'UH': 'ʊ', 'UW': 'u', 'V': 'v', 'W': 'w',
'Y': 'j', 'Z': 'z', 'ZH': 'ʒ'
}
def arpabet_phone_to_ipa(phone: str) -> str:
"""Convert one stress-marked ARPABET phone to IPA."""
if phone == 'AH0':
return 'ə'
return arpabet_to_phoneme[re.sub(r'\d', '', phone)]
entries: dict[str, list[str]] = {
word: [arpabet_phone_to_ipa(phone) for phone in phones]
for word, phones in cmudict.entries()
if len(set(word)) > 1
if len(phones) > 1
if not re.findall(r'[^a-z]', word)
}
len(entries)
```
Next, we define the $N$-gram model. The implementation stores log-probabilities to avoid numerical underflow.
```{python}
#| code-fold: true
#| code-summary: Define `NgramModel`
from collections.abc import Iterable, Sequence
from collections import defaultdict, Counter
from itertools import product
from math import isfinite
from numbers import Real
from typing import Self
alphabet: set[str] = {
phone for phones in entries.values() for phone in phones
}
class NgramModel:
"""An N-gram language model with add-lambda smoothing.
Parameters
----------
alphabet : set[str]
The set of symbols in the language.
n : int
The order of the model (default 2 for bigram).
lam : float
The smoothing pseudocount (default 0 for MLE).
"""
def __init__(
self,
alphabet: set[str],
n: int = 2,
lam: float = 0.,
) -> None:
if isinstance(n, bool) or not isinstance(n, int):
raise TypeError("n must be an integer")
if n < 1:
raise ValueError("n must be at least 1")
if isinstance(lam, bool) or not isinstance(lam, Real):
raise TypeError("lam must be a real number")
if not isfinite(lam) or lam < 0:
raise ValueError("lam must be finite and nonnegative")
self._alphabet = set(alphabet)
self._n = n
self._lam = float(lam)
self._logprob: dict[tuple[str, ...], dict[str, float]] = {}
def predict(self, string: Sequence[str]) -> float:
"""Compute the log-probability of a string.
Parameters
----------
string : Sequence[str]
The phone sequence to score.
Returns
-------
float
The log-probability of the string under the model.
"""
padded = ['<s>'] * (self._n - 1) + list(string) + ['</s>']
logprob = 0.
for i in range(self._n - 1, len(padded)):
context = tuple(padded[i - self._n + 1:i])
token = padded[i]
if self._n == 1:
context = ()
if token in ('<s>',):
continue
if context in self._logprob and token in self._logprob[context]:
logprob += self._logprob[context][token]
else:
logprob += -np.inf
return logprob
def fit(self, lexicon: Iterable[Sequence[str]]) -> Self:
"""Fit the model to a lexicon of phone sequences.
Parameters
----------
lexicon : Iterable[Sequence[str]]
An iterable of phone sequences to train on.
Returns
-------
NgramModel
The fitted model (self).
"""
vocab = self._alphabet | {'</s>'}
freq: dict[tuple[str, ...], Counter[str]] = {}
# Contexts that start with <s> padding
for k in range(1, self._n):
for ctx in product(*[self._alphabet] * (k - 1)):
full_ctx = tuple(['<s>'] * (self._n - k)) + ctx
freq[full_ctx] = Counter({v: 0 for v in vocab})
# All other contexts
for ctx in product(*[self._alphabet] * (self._n - 1)):
freq[ctx] = Counter({v: 0 for v in vocab})
# Unigram context
if self._n == 1:
freq[()] = Counter({v: 0 for v in vocab})
# Count N-grams
for word in lexicon:
padded = ['<s>'] * (self._n - 1) + list(word) + ['</s>']
for i in range(self._n - 1, len(padded)):
context = tuple(padded[i - self._n + 1:i])
token = padded[i]
if self._n == 1:
context = ()
if token == '<s>':
continue
if context in freq:
freq[context][token] += 1
# Compute log-probabilities with add-λ smoothing
self._logprob = {}
for ctx, counts in freq.items():
denominator = sum(counts.values()) + len(counts) * self._lam
# With λ = 0, an unseen context has no identified MLE.
if denominator == 0:
continue
self._logprob[ctx] = {
tok: (
np.log(count + self._lam) - np.log(denominator)
if count + self._lam > 0
else -np.inf
)
for tok, count in counts.items()
}
return self
```
The implementation preserves the normalization argument above by using the same denominator for every continuation of a fixed context. With $\lambda>0$, every enumerated context receives all $K$ continuation probabilities. With $\lambda=0$, a context with positive total count receives its MLE probabilities, including $-\infty$ log-probability for zero-count continuations. A context with zero total count is omitted because its MLE is unidentified; `predict` then assigns $-\infty$ when a string requires that context. Thus, the code distinguishes an undefined unsmoothed conditional estimate from a `NaN` produced by the indeterminate expression $\log 0-\log 0$.
The following checks exercise the parameter boundary and the termination count. A unigram trained on the one-symbol string `a` sees two continuation tokens—`a` and `</s>`—so its two fitted probabilities are each $1/2$.
```{python}
#| code-fold: true
#| code-summary: Check N-gram parameter and termination invariants
toy_unigram = NgramModel({'a'}, n=1, lam=0).fit([['a']])
assert np.isclose(np.exp(toy_unigram._logprob[()]['a']), 0.5)
assert np.isclose(np.exp(toy_unigram._logprob[()]['</s>']), 0.5)
for invalid_n in (0, -1, 1.5, True):
try:
NgramModel({'a'}, n=invalid_n)
except (TypeError, ValueError):
pass
else:
raise AssertionError(f"accepted invalid n={invalid_n!r}")
for invalid_lam in (-0.1, np.inf, np.nan, True):
try:
NgramModel({'a'}, lam=invalid_lam)
except (TypeError, ValueError):
pass
else:
raise AssertionError(f"accepted invalid lam={invalid_lam!r}")
```
We'll fit models with $N \in \{1, 2, 3, 4\}$ and a range of smoothing values $\lambda \in \{0, 0.1, 0.5, 1, 2, 5, 10\}$.
```{python}
models = defaultdict(dict)
for n in range(1, 5):
for lam in [0., 0.1, 0.5, 1., 2., 5., 10.]:
models[n][lam] = NgramModel(alphabet, n=n, lam=lam).fit(entries.values())
```
## A few test cases
To get an initial sense for how these models behave, we can look at the probabilities they assign to a few test words. These are constructed to have the same vowel and coda but different onsets, ranging from well-attested ($\text{blɪk}$) to unattested ($\text{bnɪk}$, $\text{bzɪk}$, $\text{bsɪk}$).
```{python}
test_words = [
['b', 'l', 'ɪ', 'k'],
['b', 'n', 'ɪ', 'k'],
['b', 'z', 'ɪ', 'k'],
['b', 's', 'ɪ', 'k'],
]
for n in [2, 3]:
for lam in [0., 1.]:
print(f'N={n}, λ={lam}:')
for w in test_words:
print(f' {"".join(w):6s} {models[n][lam].predict(w):8.2f}')
print()
```
Notice that the models with $N > 1$ consistently assign higher probability to $\text{blɪk}$ (which has an attested onset cluster) than to $\text{bnɪk}$, $\text{bzɪk}$, or $\text{bsɪk}$ (which do not). The unigram model doesn't make this distinction, since it doesn't condition on previous phones at all.
## Predicting acceptability
Now let's ask a more systematic question: how well do these models predict the acceptability judgments from @daland_explaining_2011? We load the acceptability data and convert the ARPABET representations to IPA.
```{python}
#| code-fold: true
#| code-summary: Load acceptability data and convert to IPA
acceptability = pd.read_csv('Daland_etal_2011__AverageScores.csv')
acceptability['phono_cmu'] = (
acceptability.phono_cmu
.str.replace(r'Coda$', '', regex=True)
.str.split()
)
def map_to_ipa(phones: Sequence[str]) -> list[str]:
"""Convert ARPABET phone list to IPA symbols.
Parameters
----------
phones : Sequence[str]
ARPABET phone symbols, possibly with stress digits.
Returns
-------
list[str]
The corresponding IPA symbols.
"""
return [arpabet_phone_to_ipa(phone) for phone in phones]
acceptability['phono_ipa'] = acceptability.phono_cmu.map(map_to_ipa)
```
The approach is straightforward. For each model, we compute the log-probability of each nonword in the dataset, then ask how well that log-probability predicts the human acceptability rating. We test the hypothesis that the relationship is linear:
$$a_\mathbf{w} \sim \mathcal{N}(m \log p(\mathbf{w} \mid \boldsymbol\Theta) + b, \sigma^2)$$
Because all the nonwords in the dataset happen to be the same length, we don't need to worry about length normalization here—though in general, we would need to, since longer strings tend to have lower log-probabilities simply because more terms are being multiplied together.
### Evaluating a predictor
Before we run the evaluation, it's worth spelling out the machinery we'll use, since the same ideas will recur throughout these notes.
#### Linear regression
The model above says that a nonword's acceptability rating is, in expectation, a linear function of its log-probability. We defined [conditional expectation](some-more-useful-definitions.qmd#conditional-expectation) earlier as $\mathbb{E}[X \mid Y = y] = \sum_x x \cdot p(x \mid y)$; the claim here is that $\mathbb{E}[a_\mathbf{w} \mid \log p(\mathbf{w} \mid \boldsymbol\Theta) = x] = mx + b$ for some slope $m$ and intercept $b$. The Gaussian noise term $\sigma^2$ captures everything about the rating that the log-probability doesn't predict.
To *fit* this model—to find $m$ and $b$—we use [ordinary least squares](https://en.wikipedia.org/wiki/Ordinary_least_squares): choose the values of $m$ and $b$ that minimize the sum of squared residuals across the observed data:
$$\hat{m}, \hat{b} = \underset{m, b}{\text{argmin}} \sum_{i=1}^{n} \left(a_i - (m x_i + b)\right)^2$$
where $a_i$ is the observed rating for nonword $i$ and $x_i$ is its log-probability. This has a closed-form solution (the [normal equations](https://en.wikipedia.org/wiki/Ordinary_least_squares#Derivation_of_simple_linear_regression_estimators)), but the conceptual point is more important than the algebra: we are choosing the line through our data that makes the observed ratings least surprising under the Gaussian noise model stated above. See @hastie_elements_2009[Ch. 3] for a thorough treatment of [linear regression](https://en.wikipedia.org/wiki/Linear_regression) and its extensions.
#### Coefficient of determination
Once we've fitted the line, we need a way to measure how well it predicts. The [coefficient of determination](https://en.wikipedia.org/wiki/Coefficient_of_determination) $R^2$ does this by comparing the model's errors to the simplest possible baseline—predicting the mean rating $\bar{a}$ for every nonword:
$$R^2 = 1 - \frac{\sum_{i=1}^{n} (a_i - \hat{a}_i)^2}{\sum_{i=1}^{n} (a_i - \bar{a})^2}$$
The denominator is the total variance in the ratings. The numerator is the variance the model leaves unexplained. So $R^2$ is the fraction of variance *explained* by the predictor. If $R^2 = 1$, every prediction is exact; if $R^2 = 0$, the model does no better than predicting the mean.^[An $R^2$ of, say, 0.3 does not mean the model is "bad"—it means 30% of the variance in the response is predictable from the predictor. Whether that is a lot or a little depends on the domain. In psycholinguistics, where individual responses are noisy, $R^2$ values in the 0.1–0.4 range for item-level predictions are common.]
For simple linear regression with a single predictor, $R^2$ is exactly the square of the [Pearson correlation](some-more-useful-definitions.qmd#covariance-and-correlation) between the predictor and the response: $R^2 = \text{corr}(X, Y)^2$. This connection to the correlation we defined earlier is useful, but it breaks down when we evaluate on new data—where, as we'll see, $R^2$ can actually go negative.
#### Cross-validation
If we compute $R^2$ on the same data we used to fit $m$ and $b$, we get an optimistic estimate of how well that fitted calibration generalizes. The regression is free to latch onto idiosyncratic properties of the judgment data, inflating its apparent performance relative to performance on new judgment items. This is [overfitting](https://en.wikipedia.org/wiki/Overfitting). The $N$-gram scores used here have already been fitted to a separate lexicon, so the resampling below evaluates the regression that maps those fixed scores to judgments; it does not refit the $N$-gram counts or measure their sensitivity to a new lexical sample.
[Cross-validation](https://en.wikipedia.org/wiki/Cross-validation_(statistics)) guards against this [@stone_cross-validatory_1974]. The idea is simple: partition the judgment items into $K$ non-overlapping subsets called *folds*. For each fold, fit the linear regression on the other $K - 1$ folds (the *training set*) and evaluate $R^2$ on the held-out fold (the *test set*). Each regression fit thus uses approximately the fraction $(K-1)/K$ of the items for training. Rotate through all $K$ folds and report the average $R^2$ across them.^[The choice of $K = 5$ is a common default that balances bias and variance in the estimate. Smaller $K$ uses a smaller training fraction and places more items in each test fold; larger $K$ (up to leave-one-out, where $K = n$) uses a larger training fraction but leaves fewer evaluation items per fold and increases computation. See @hastie_elements_2009[Sec. 7.10] for analysis of this tradeoff.]
Because the regression is always evaluated on judgment items it hasn't seen, the cross-validated $R^2$ is a more honest estimate of its predictive calibration. It can go negative: under the held-out-fold definition used by `scikit-learn`, that means the predictions are worse than the baseline that assigns every item the held-out fold's mean. Comparing many choices of $N$ and $\lambda$ on these same folds is still model selection; a final performance estimate for the selected setting would require nested cross-validation or a separate test set. See @hastie_elements_2009[Ch. 7] for an extended discussion of cross-validation and related model-selection methods.
#### Descriptive fold-resampling bands
Cross-validation gives us $K$ scores, one per held-out fold. The code below resamples those observed scores and plots percentiles of their means. We call these **fold-resampling bands** (FRBs), not confidence intervals. The five fold scores are not independent because their training sets overlap, and five values do not provide a reliable empirical approximation to a sampling distribution. Thus, an FRB summarizes how the displayed mean changes when the observed folds receive different multiplicities; it is not a calibrated 95% interval for performance on new items. An item-level bootstrap that repeats the entire fitting and evaluation procedure, or repeated cross-validation with a prespecified inferential target, would be needed for that stronger interpretation.
With these tools in hand, we evaluate each model.
```{python}
#| code-fold: true
#| code-summary: Cross-validation
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
acceptability_shuffled = acceptability.sort_values(
'likert_rating'
).sample(frac=1., random_state=40393, replace=False)
cv_stats = []
all_logprobs = []
rng = np.random.default_rng(40393)
for n in range(1, 5):
for lam in [0.1, 0.5, 1., 2., 5., 10.]:
logprob = acceptability_shuffled.phono_ipa.map(models[n][lam].predict)
assert np.isfinite(logprob).all()
logprob_df = pd.DataFrame({
'logprob': logprob.values,
'likert_rating': acceptability_shuffled.likert_rating.values
})
logprob_df = logprob_df[~logprob_df.logprob.isnull()]
logprob_df['n'] = n
logprob_df['lam'] = lam
all_logprobs.append(logprob_df)
scores = cross_val_score(
LinearRegression(),
logprob_df.logprob.values[:, None],
logprob_df.likert_rating.values,
cv=5
)
stats = np.quantile(
[np.mean(rng.choice(scores, size=5)) for _ in range(2_000)],
[0.5, 0.025, 0.975]
)
cv_stats.append([n, lam] + list(stats))
all_logprobs = pd.concat(all_logprobs, axis=0)
cv_stats = pd.DataFrame(
cv_stats,
columns=['n', 'lam', 'r2', 'band_lo', 'band_hi']
)
```
The unsmoothed models are omitted from this regression analysis. They assign probability zero, and hence log-probability $-\infty$, to some test strings. Replacing those values with a data-dependent finite cutoff would define neither the original model nor a normalized floor model. Every evaluated model instead has $\lambda>0$, so every enumerated transition has finite log-probability.
The following plot shows the relationship between log-probability and acceptability rating for each combination of $N$ and $\lambda$.
```{python}
#| code-fold: true
#| code-summary: Scatter plots of log-probability vs. acceptability
import matplotlib.pyplot as plt
fig, axes = plt.subplots(4, 6, figsize=(18, 12), sharex=False, sharey=True)
lam_values = [0.1, 0.5, 1., 2., 5., 10.]
for row, n in enumerate(range(1, 5)):
for col, lam in enumerate(lam_values):
ax = axes[row, col]
subset = all_logprobs[(all_logprobs.n == n) & (all_logprobs.lam == lam)]
ax.scatter(subset.logprob, subset.likert_rating, alpha=0.3, s=8)
if len(subset) > 0:
m, b = np.polyfit(subset.logprob, subset.likert_rating, 1)
x_range = np.linspace(subset.logprob.min(), subset.logprob.max(), 50)
ax.plot(x_range, m * x_range + b, color='C1', linewidth=1.5)
if row == 0:
ax.set_title(f'λ={lam}', fontsize=10)
if col == 0:
ax.set_ylabel(f'N={n}', fontsize=10)
ax.tick_params(labelsize=7)
fig.supxlabel('Log-probability', fontsize=12)
fig.supylabel('Likert rating', fontsize=12)
fig.tight_layout()
plt.show()
```
And the following plot shows the cross-validated $R^2$ as a function of $\lambda$ for each value of $N$, with descriptive FRBs.
```{python}
#| code-fold: true
#| code-summary: Cross-validation performance
fig, axes = plt.subplots(1, 4, figsize=(14, 3.5), sharey=True)
for i, n in enumerate(range(1, 5)):
ax = axes[i]
subset = cv_stats[cv_stats.n == n]
ax.plot(subset.lam, subset.r2, marker='o', markersize=4)
ax.fill_between(subset.lam, subset.band_lo, subset.band_hi, alpha=0.15)
ax.set_title(f'N = {n}', fontsize=11)
ax.set_xlabel('λ', fontsize=10)
if i == 0:
ax.set_ylabel('$R^2$ (5-fold CV)', fontsize=10)
ax.tick_params(labelsize=8)
fig.tight_layout()
plt.show()
```
```{python}
cv_stats.sort_values('r2', ascending=False).head(10)
```
The main pattern is that the unigram model ($N = 1$) does poorly regardless of $\lambda$, which makes sense: it has no way of capturing the fact that some sequences of phones are more acceptable than others, since it treats each phone independently. The bigram and trigram models do substantially better, and the best-performing models tend to use moderate smoothing. Too little smoothing and the model overfits to the training data; too much and it washes out the distinctions it learned.
A simple bigram or trigram model trained on a dictionary can predict a reasonable amount of variance in human acceptability judgments. This suggests that a good deal of phonotactic knowledge can be captured by tracking short-range sequential regularities in the lexicon, consistent with the lexicalist approach that @daland_explaining_2011 investigate.