Statistical models, estimators and maximum likelihood - Data Science Lecture Notes from Konrad Völkel

In the data science process, the road from the data to a product (be it a report, a prediction or new and interesting data), we model the data. Strictly speaking we model partial aspects, sections of the data (just as the data themselves are only a section of reality). In doing so we may be interested in how the data came about, but we may equally be driven by the wish to be able to say something particular about them, or to do something particular with them. So it helps to be clear about what kind of information one wants to coax out of the data, and then to tune the models to that. Often several approaches are available for computing the same kind of information. The results may well differ, and frequently (in unsupervised learning) we cannot compare the computed results with the “right” result. All that is left then is to compare the information content of the models (to compute the entropy), to examine every model assumption for plausibility, and to develop general criteria for the quality of models. The last of these is what we now set about.

The classical approach – classical in the sense of historically oldest – begins by formalising the models and the aspects one wants to compute from the data with them. One speaks of a statistical model, of a statistic, and of various estimators for a statistic which can then be compared with one another. One criterion is whether an estimator lands near the “right” result at least for very large samples. But hold on – if we do not know the right result, how can we say anything of the sort?

The solution is to restrict ourselves to those statistical models in which we can describe the rough structure of the right result for an estimator well enough that statements such as “in the limit this estimator estimates the expected value of that statistic on the data” become available after all. And then one can talk about how good the values look away from the limit, for somewhat smaller samples. Along the way we meet a notion that matters a great deal to data science and machine learning: the maximum likelihood principle.

A very concrete example first. Given data that come as a vector $x$ of length $n$ of real numbers ( $\texttt{float}$ ), of which we know that they arose as a sample of a one-dimensional normally distributed random variable $X$ with expected value $\mu(X)$ and variance $\sigma^2(X)$ , we want to estimate the standard deviation $\sigma(X)$ . For this the books offer two different formulas, $\sigma(x)$ and $\sigma^\ast(x)$ , which differ only in their denominator:

\sigma(x)^2 = \dfrac{1}{n} \sum_{i=1}^n (x_i-\mu)^2,\qquad \sigma^\ast(x)^2 = \dfrac{1}{n-1} \sum_{i=1}^n (x_i-\mu)^2

(1)

Quite concretely, for a sample of $n = 48$ values that really were drawn from a normal distribution:

import numpy as np

x = np.array([6.5281, 3.8003, 4.9575, 7.4818, 6.7351, 1.0454, 4.9002, 2.6973,
              2.7936, 3.8212, 3.2881, 5.9085, 4.5221, 3.2434, 3.8877, 3.6673,
              5.9882, 2.5897, 3.6261, 1.2918, -2.1060, 4.3072, 4.7289, 1.5157,
              7.5395, 0.0913, 3.0915, 2.6256, 6.0656, 5.9387, 3.3099, 3.7563,
              1.2244, -0.9616, 2.3042, 3.3127, 5.4606, 5.4048, 2.2253, 2.3954,
              0.9029, 0.1600, -0.4125, 6.9016, 1.9807, 2.1239, 0.4944, 4.5550])

print(np.mean(x))            # 3.3689...
print(np.var(x))             # 5.0037...  denominator n
print(np.var(x, ddof=1))     # 5.1102...  denominator n-1
print(np.std(x))             # 2.2369...
print(np.std(x, ddof=1))     # 2.2606...

So the two formulas deliver noticeably different numbers – a good two per cent apart in the variance, and that at $n = 48$ . Which is the right one? And in passing: that we may replace $\mu(X)$ by $\texttt{mean}(x)$ at all is something we have so far only taken to be plausible; the justification is the example “the sample mean is unbiased” in the section Unbiasedness below.

With these formulas one has to be a little careful when implementing them in practice – because what is $\mu$ ? In an application I do not know $\mu(X)$ and must therefore use $\mu(x)$ . But then it turns out that the information in $\{x_i-\mu(x)\}_{i=1}^{n-1}$ suffices to compute $x_n - \mu(x)$ , since by the definition of $\mu(x)$ we have $x_n - \mu(x) = \sum_{i=1}^{n-1} (\mu(x) - x_i)$ . In a sense, then, using $\mu(x)$ in place of $\mu(X)$ leaves us one degree of freedom fewer than $n$ . We shall see in a moment how that affects the properties of these two estimators for $\sigma(X)$ .

Once statistical models have been defined formally, one can also talk very well about confidence intervals and hypothesis tests, which is what the chapter Statistical tests does. Confidence intervals extend the bare estimation of a statistic by the information on what margin of error one attaches to that estimate. Hypothesis tests extend it by arriving at decisions in the end. In discussing hypothesis tests it is unavoidable to speak also about the problematic notion of significance. But one thing at a time!

As a map for the definitions that follow – on the left what one talks about, on the right what it formally becomes:

in everyday terms	formally
a family of distributions from which we want to pick one	statistical model $(\mathcal{X}, P_\theta : \theta \in \Theta)$
the quantity we want to estimate (say $\mu$ or $\sigma$ )	map $\tau \colon \Theta \to \Sigma$
the sample we are holding in our hands	an element $x \in \mathcal{X}$ – or, as long as it has not been drawn yet, a random variable with values in $\mathcal{X}$
the procedure that turns the sample into a number	statistic, or estimator, $T \colon \mathcal{X} \to \Sigma$

The sample space $\mathcal{X}$ here is not the space of a single measured value but the space of all samples: with $n$ real measurements $\mathcal{X} = \mathbb{R}^n$ , so a sample is a single point in a fairly large space.

The crucial move is in the last row. Because the sample is a random variable and the estimation procedure is a map on it, the result of estimating is again a random variable. That lets us speak of its expected value and its variance – and that is exactly what we need in order to be able to say at all what a good estimator is: bias, consistency and mean squared error are all formulated in this language.

Statistical models¶

Remark

A statistical model is therefore not yet the ‘finished model’ with which we want to classify or make predictions in the data science sense, but a whole class of models parametrised by the set $\Theta$ . This set $\Theta$ may be as wild as you like; it should only have at least two elements, so that in the case $\Theta = \{\theta\}$ the notion does not shrink to ‘the sample space with one probability measure on it’. When speaking of expected values (and hence, say, of variances) in the context of a statistical model, one has to be careful, and it is best to carry the $\theta$ along in the notation: $\mathbb{E}_\theta(Y)$ for the expected value of a random variable $Y : \mathcal{X} \to \mathbb{R}$ , because that value depends on $\theta$ .

Definition

A statistical model is called a parametric model if the parameter set $\Theta$ is a subset of some $\mathbb{R}^d$ for a $d \in \mathbb{N}$ . For $d=1$ one also says one-parameter model.

If the sample space $\mathcal{X}$ is discrete (at most countable) and $\mathcal{F} = \mathcal{P}(\mathcal{X})$ , the model is called discrete, and in that case every $P_\theta$ has a discrete density $p_\theta$ , that is $p_\theta(x) = P_\theta(\{x\})$ .

If the sample space $\mathcal{X}$ is a (Borel) measurable subset of $\mathbb{R}^n$ (an open set, say, or the level set of a continuous function), if $\mathcal{F}$ is a restriction of the Borel $\sigma$ -algebra of $\mathbb{R}^n$ , and if moreover every $P_\theta$ has a density function $f_\theta$ , then the statistical model is called a continuous model.

A statistical model that is either discrete or continuous is also called a standard model.

Remark

To a first approximation we can picture the sample space as the set of possible outcomes at sample size 1 (the population $\Omega$ sends its regards). But: one can also relabel the repeated running of an experiment as really-one-experiment, call a sequence of individual outcomes a single outcome, and so model sample sizes greater than 1. For every statistical model the product model can be constructed:

(\mathcal{X}^N, \mathcal{F}^{\otimes N}, P_{\theta}^{\otimes N} : \theta \in \Theta).

(2)

The projections $X_i \colon \mathcal{X}^N \to \mathcal{X}$ are then independent, identically distributed random variables with respect to every distribution $P_\theta$ . The notation $\mathcal{F}^{\otimes N}$ denotes the obvious way of defining a $\sigma$ -algebra on $\mathcal{X}^N$ (you may ignore it, as usual), and $P_\theta^{\otimes N}$ assigns to an event in $\mathcal{X}^N$ of the form $A_1 \times \cdots \times A_N$ the product $\prod_{i=1}^N P_\theta(A_i)$ . With that the $N$ -fold repetition is modelled.

Example

We can also model whether the two throws came up as doubles or not – and so take $\mathcal{X} = \{\text{doubles}, \text{no doubles}\}$ as the sample space. Here we can simply keep the parameter space $\Theta = \Delta_6$ ; that is one possible statistical model. We do see, though, that $P_\theta^{\otimes 2}(k,k) = \theta_k^2$ , so that $P_\theta(\text{doubles}) = \sum_{k=1}^6 \theta_k^2$ , and hence the parameter $s = \sum_{k=1}^6 \theta_k^2$ really suffices and we could imagine a smaller parameter space. But then the question arises which distribution $P_s$ we want to assign to the parameter $s$ . One concretely sensible choice is the distribution $\theta$ with maximal entropy, that is, a distribution as close to the uniform distribution as possible. This does not become unique, because permutations of $\{1,\dots,6\}$ change neither the parameter $s$ nor the entropy $H(\theta)$ . If one nevertheless wants to pick a model that is canonical in some way, one can for instance demand an ordering on the $\theta$ , that is $\theta_k \geq \theta_{k+1}$ . So we see here: a ‘too large’ parameter space is sometimes the simpler model.

Example

Imagine this. We have found some data $x$ in the wild (read: on the web), now believe that all instances $x_i$ are samples of a random variable $X$ with values in $\Omega$ , and so model with $\mathcal{X} = \Omega^n$ . If we notice that the values of $x$ are floats of some kind, we can look at the histogram. Perhaps we do not get much out of it and cannot really tell whether $X$ is a variable uniformly distributed on $[-1,1]$ , or a normal distribution centred at 0, or even an exponential distribution. In reality that could happen, for instance, when $n$ is very small. The parameter space would then be

\Theta = \{\text{uni}\} \times \{0\} \cup \{\text{norm}\} \times \{\sigma > 0\} \cup \{\text{exp}\} \times \{\lambda > 0\}.

(3)

The first component of $\Theta$ is the type of the distribution, $\{\text{uni, norm, exp}\}$ , the second then a possible parameter. Such a model is not a parametric model! We can, though, make one out of it in an unnatural way, by encoding the categorical variable $\{\text{uni, norm, exp}\}$ as the numbers $\{1,2,3\}$ and thus representing $\Theta$ as a subset of $\mathbb{R}^2$ .

Estimators¶

Example

In most statistical models we can consider the map $\sigma \colon \Theta \to \mathbb{R}_{\geq 0}$ that assigns to a parameter $\theta$ the standard deviation $\sigma(P_\theta)$ of the distribution $P_\theta$ (it only fails when the second moment does not exist, that is when $\mathbb{V}(P_\theta)=\infty$ occurs). For it we know the two estimators $\sigma,\ \sigma^\ast \colon \mathcal{X} \to \mathbb{R}_{\geq 0}$ . The two are related by $\sigma^\ast(x) = \sqrt{\frac{n}{n-1}}\,\sigma(x)$ , equivalently $\sigma^\ast(x)^2 = \frac{n}{n-1}\,\sigma(x)^2$ ; this factor is called the Bessel correction (cf. Descriptive statistics). Which of the two estimators is the better one is settled by the next section.

Unbiasedness¶

Definition

An estimator $T \colon \mathcal{X} \to \Sigma$ for $\tau \colon \Theta \to \Sigma$ is called unbiased if $\mathbb{E}_{\theta}( T ) = \tau(\theta)$ holds for all $\theta \in \Theta$ . If it does not, the estimator is called biased, and the difference is the bias $\mathbb{B}_\theta(T) := \mathbb{E}_{\theta}( T ) - \tau(\theta)$ .

The mean squared error of the estimator against the value to be estimated is

\texttt{mse}(T,\theta) := \mathbb{E}_\theta\left ( (T-\tau(\theta))^2 \right ) = \mathbb{V}_\theta(T) + (\mathbb{B}_\theta(T))^2

(4)

A sequence of estimators $T_n \colon \mathcal{X}^n \to \Sigma$ for $\tau$ in $n$ -fold product models is called consistent if for $n\to\infty$ the estimator converges to $\tau$ in probability: $T_n \xrightarrow{P_\theta} \tau(\theta)$ . Convergence in probability means $\forall \epsilon > 0 : \lim_{n\to\infty} P(|T_n - \tau(\theta)| > \epsilon) = 0$ .

The sequence $T_n$ is called asymptotically unbiased if $\lim_{n\to\infty} \mathbb{E}_{\theta}(T_n) = \tau(\theta)$ .

Example

For the variance the estimator $\sigma_N^2 = \frac{1}{N} \sum_{i=1}^N (X_i - \texttt{mean}(X_j))^2$ is not unbiased, since $\mathbb{E}(\sigma_N^2) = \frac{N-1}{N}\sigma^2(X)$ . The bias is therefore $-\frac{\sigma^2(X)}{N}$ , which vanishes as $N\to\infty$ . So the estimator is asymptotically unbiased.

The estimator adjusted by the Bessel correction, $\frac{N}{N-1} \sigma_N^2 = \frac{1}{N-1} \sum_{i=1}^N (X_i - \texttt{mean}(X_j))^2$ , is unbiased. So we say: $\sigma_N^2$ is asymptotically unbiased, the corrected estimator even unbiased (free of bias).

The bias arises solely from the fact that $\mu$ is estimated along with everything else. If the expected value $\mu(X)$ is known, then already

\dfrac{1}{N} \sum_{i=1}^N {\left(X_i - \mu(X)\right)}^2

(5)

is unbiased, since $\mathbb{E}\left((X_i-\mu(X))^2\right) = \sigma^2(X)$ holds for each individual $i$ and the expected value is linear. Only putting $\texttt{mean}(x)$ in place of $\mu(X)$ costs the degree of freedom that the introduction to this chapter spoke of – and it is exactly that one degree of freedom which the factor $\frac{N}{N-1}$ wins back.

Such an estimator therefore has the smallest spread around the value to be estimated among all unbiased estimators.

In general it is very hard to determine uniformly best unbiased estimators (in general there are none). Looking ahead: the Cramér–Rao inequality for regular statistical models bounds the variance of an estimator from below by the Fisher information (a quantity that can also be computed from the relative entropy), and so lets us recognise uniformly best estimators by the fact that they do not exceed this variance bound. For an unbiased estimator $T$ of a real parameter $\theta$ it reads

\mathbb{V}_\theta(T) \geq \frac{1}{I(\theta)}, \qquad I(\theta) := \mathbb{E}_\theta\left[\left(\frac{\partial}{\partial\theta} \log p_\theta(X)\right)^2\right],

(6)

where $I(\theta)$ is the Fisher information. An unbiased estimator that attains this bound is therefore a uniformly best estimator. For distributions from a so-called exponential family (in a certain sense a common generalisation of the normal distribution, the binomial distribution, the gamma distribution and so on), the uniformly best estimators can be described directly in this way.

For many models that matter in practice, though, there is no closed formula for a good estimator at all – the contaminated normal distribution is the example that still awaits us in the course. What one needs then is a principle that yields, instead of a formula, an optimisation problem whose solution can be approximated numerically. That is exactly what the maximum likelihood principle does, and the iterative procedure that goes with it is the EM algorithm.

Maximum likelihood¶

So if we are given data $x$ , the likelihood function for $x$ expresses, for every parameter $\theta$ , how the probability density of $P_\theta$ comes out at the point $x$ . If the density is large, then $x$ is an observation that is probable under the parameter $\theta$ . If the density is small, then $x$ is improbable under $\theta$ . As a reminder: if we take the hypothesis $H = ( \theta$ is the true parameter $)$ and write $D=x$ (for data), then $P(D|H) = L_x(\theta) = f_\theta(x)$ is the likelihood (cf. Bayes’ theorem, version 3, in the chapter Covariance and dependence).

We have already worked out a first example of this without using the name: in the chapter The most important discrete distributions we determined, for $k$ successes in $n$ trials, the $p$ under which that observation is most probable, and obtained $\hat p = \frac{k}{n}$ . In the language of this chapter that is the maximum likelihood estimator in the binomial model $\mathcal{X} = \{0,\dots,n\}$ , $\Theta = [0,1]$ , $P_p = \operatorname{Bin}(n,p)$ .

Remark

So for a family of distributions $P_\theta$ with densities $f_\theta$ we are interested in maximising the likelihood $f_\theta(x)$ of an observation $x$ in $\theta$ . To maximise functions we like to look at the derivative (for which $f$ should be differentiable in $\theta$ ). Already in the example of the univariate normal density $\phi_{\mu,\sigma}$ we see that the derivative $\frac{d}{d\mu} \phi_{\mu,\sigma}(x) =\frac{(x-\mu)}{\sigma^2} \phi_{\mu,\sigma}(x)$ vanishes at $x=\mu$ . So $\hat\mu = x$ is the maximum likelihood estimator for $\mu$ for a sample of size 1 – which one believes without any calculation: if you observe exactly one value, the normal distribution under which that value is most probable is the one whose expected value sits at exactly that point. That this is a maximum and not a minimum can be seen from the fact that $\phi_{\mu,\sigma}(x)$ tends to 0 as $\mu \to \pm\infty$ . In general, though, the likelihood is harder to compute (think of an $n$ -fold product model!), and one can help oneself with a simple trick:

Trick

Instead of maximising the likelihood of $x$ over all possible $\theta$ , one can also take the logarithm of the likelihood, the so-called log-likelihood $\log \rho(x,\theta)$ , and maximise that. Because the logarithm is strictly increasing, likelihood and log-likelihood attain their maxima in the same places.

If, for instance, $\rho$ has the form $e^{A}e^{B}$ (as in a 2-fold product model), then $\log\rho$ has the form $A+B$ and setting derivatives to 0 becomes a great deal easier. Even with this trick, however, no closed analytic solution can be found in general – which is where computer science gets to shine with approximation algorithms, such as the expectation maximization algorithm, which we take up later in the course.

Example

Maximum likelihood for the exponential distribution. We observe $n$ independent waiting times $x_1,\dots,x_n \geq 0$ and model them as $\operatorname{Exp}(\lambda)$ -distributed, with density $f_\lambda(x) = \lambda e^{-\lambda x}$ . Because of the independence the likelihood is the product of the individual densities:

L(\lambda) = \prod_{i=1}^n \lambda e^{-\lambda x_i} = \lambda^n e^{-\lambda \sum_{i=1}^n x_i}.

(9)

Here the trick pays off at once – the product turns into a sum:

\ell(\lambda) := \log L(\lambda) = n\log\lambda - \lambda \sum_{i=1}^n x_i.

(10)

Differentiating and setting equal to 0 gives

\ell'(\lambda) = \frac{n}{\lambda} - \sum_{i=1}^n x_i \overset{!}{=} 0 \qquad\Longleftrightarrow\qquad \hat\lambda = \frac{n}{\sum_{i=1}^n x_i} = \frac{1}{\texttt{mean}(x)}.

(11)

Since $\ell''(\lambda) = -\frac{n}{\lambda^2} < 0$ , this really is a maximum. The result is intuitive again: the exponential distribution has expected value $\frac{1}{\lambda}$ , so we estimate $\frac{1}{\lambda}$ by the sample mean. For the observations $1, 2, 3, 4$ , say, $\texttt{mean}(x) = 2.5$ and hence $\hat\lambda = 0.4$ .

Example

The k-means algorithm (see the chapter k-Means) can be read as an estimator: the statistical model is $\mathcal{X} = \mathbb{R}^{nN}$ , where $n$ is the dimension of the data (so with 4 real variables, as in the Iris dataset, $n=4$ ) and $N$ the sample size (for Iris, $N=150$ ).

The parameter space is $\mathbb{R}^{nk}$ , and the distribution belonging to a parameter $\theta = (c_1,\dots,c_k) \in \mathbb{R}^{nk}$ (with $c_i \in \mathbb{R}^n$ ) is a mixture model of multivariate normal distributions with expected value $c_i$ and the $n\times n$ identity matrix as covariance matrix.

The ‘perfect’ k-means solution, that is the centroids $(c_1,\dots,c_k)$ that the k-means algorithm tries to determine, are those for which the sum of the squared Euclidean distances from every point $x_i$ of a sample $x \in \mathcal{X}$ to its nearest centroid is minimal. These are exactly the centroids for which the likelihood of the data under the corresponding mixture distribution becomes maximal.

That leads us to an obvious generalisation of the k-means procedure: we could allow more complicated covariance matrices and so push the likelihood up further. The problem that then arises, however, is that no closed form can be determined for the log-likelihood either, and that one has to think up a new algorithm for the estimation!