31  Point estimation

In this chapter, in the sense of the terminology introduced in Chapter 30, we always assume a parametric product model \[\begin{equation} \mathcal{M} := \{\mathcal{Y},\mathcal{A}, \{\mathbb{P}_\theta| \theta \in \Theta\}\} \end{equation}\] with \(n\)-dimensional sample space (e.g. \(\mathcal{Y} := \mathbb{R}^n\)), \(d\)-dimensional parameter space \(\Theta \subset \mathbb{R}^d\), and given PMF or PDF \(p_\theta\) for all \(\theta \in \Theta\). \(y := (y_1,...,y_n)\) denotes the sample of independent and identically distributed random variables belonging to \(\mathcal{M}\); throughout, therefore, \[\begin{equation} y_1,...,y_n \sim \mathbb{P}_\theta. \end{equation}\] The nature and aim of the point estimation treated here is to give, based on the sample, a guess that is as good as possible for a characteristic quantity of the distribution \(\mathbb{P}_\theta\) of a sample variable. The guess has the same mathematical nature as the corresponding characteristic quantity, for example a scalar value for a scalar parameter. This is not the only possibility for estimation; with confidence intervals we will encounter one way of estimating scalar values by intervals. The characteristic quantities of \(\mathbb{P}_\theta\) to be estimated are often simply the true, but unknown, parameters themselves. We treat this case in detail in Section 31.1. However, many basic results of frequentist point estimation also remain valid when the quantities to be estimated are not the parameters themselves, but functions of them. Examples are the estimation of the expectation, variance, or standard deviation of \(\mathbb{P}_\theta\). We begin, however, with parameter estimation. To estimate the true, but unknown, parameter of an inference model, frequentist inference uses so-called parameter point estimators.

Definition 31.1 \(\mathcal{M} := (\mathcal{Y}, \mathcal{A}, \{\mathbb{P}_\theta|\theta \in \Theta\})\) be a frequentist inference model, let \((\Theta,\mathcal{S})\) be a measurable space, and let \(\hat{\theta} : \mathcal{Y} \to \Theta\) be a mapping. Then \(\hat{\theta}\) is called a parameter point estimator for \(\theta\).

Parameter point estimators are usually also simply called parameter estimators. In the sense of Definition 30.3, parameter point estimators are estimators with \(\tau := \mbox{id}_\Theta\). Parameter point estimators are therefore functions of data and take numerical values in the parameter space. As functions of random variables, parameter estimators are of course also random variables. Often, notation does not distinguish between \(\hat{\theta}\) as a random variable and \(\hat{\theta}(y)\) as the value of this random variable.

Definition 31.1 evidently does not specify how a parameter point estimator is to be constructed or to what extent it may then be a sensible estimator. In what follows, with maximum likelihood estimation, we first discuss a general principle that allows parameter estimators to be determined for a given frequentist inference model and that, as we will see later, guarantees certain desirable properties (Section 31.4). These properties generally concern its qualitative distributional behavior for a fixed sample size or in the limit of an infinitely large sample size. We introduce these properties in Section 31.2 and Section 31.3.

31.1 Maximum likelihood estimation

The basic idea of maximum likelihood estimation is to choose, as a guess for a true, but unknown, parameter value, that parameter value for which the probability of the observed data is maximal. For this, it is first necessary to consider the probability of observed data in a frequentist inference model as a function of the parameter in question. This is enabled and formalized by the likelihood function and its logarithm, the log-likelihood function. We define these terms here for parametric product models.

Definition 31.2 (Likelihood function and log-likelihood function) Let \(\mathcal{M}\) be a parametric product model with PMF or PDF \(p_\theta\). Then the likelihood function is defined as \[\begin{equation} L : \Theta \to [0,\infty[, \theta \mapsto L(\theta) := \prod_{i=1}^n p_\theta(y_i) \end{equation}\] and the log-likelihood function is defined as \[\begin{equation} \ell_n : \Theta \to \mathbb{R}, \theta \mapsto \ell(\theta) := \ln L(\theta). \end{equation}\]

The likelihood function is thus a function of the parameter, and its function values are the values of the joint PMF or PDF of observed data values \(y_1,...,y_n\). In general, there is no reason to assume that a likelihood function integrates to 1 over the parameter space; the likelihood function is therefore generally not a PMF or PDF. The log-likelihood function is simply the logarithm of the likelihood function. A parameter estimator obtained according to the principle of maximum likelihood estimation is now intended to maximize the likelihood function or the log-likelihood function. This leads to the following definition of a maximum likelihood estimator.

Definition 31.3 (Maximum likelihood estimator) Let \(\mathcal{M}\) be a parametric product model with parameter \(\theta \in \Theta\). A maximum likelihood estimator of \(\theta\) is defined as \[\begin{equation} \hat{\theta}^{\mbox{\tiny ML}} : \mathcal{Y} \to \Theta, y \mapsto \hat{\theta}^{\mbox{\tiny ML}}(y) := \mbox{argmax}_{\theta \in \Theta} L(\theta) = \mbox{argmax}_{\theta \in \Theta} \ell(\theta) \end{equation}\]

Note in Definition 31.3 that a maximum point of the log-likelihood function corresponds to the maximum point of the likelihood function, because the logarithm function is monotonically increasing. Working with the log-likelihood function is often easier than working directly with the likelihood function, for example when exponential functions occur in the PMF or PDF of the model. Furthermore, note that Definition 31.2 implies \[\begin{equation} \hat{\theta}^{\mbox{\tiny ML}}(y) = \mbox{argmax}_{\theta \in \Theta} \prod_{i=1}^n p_\theta(y_i) = \mbox{argmax}_{\theta \in \Theta} \sum_{i=1}^n \ln p_\theta(y_i) \end{equation}\] which makes the dependence of a maximum likelihood estimator on the data clear.

With Definition 31.3, maximum likelihood estimation is therefore the problem of determining extrema of a function. For these extrema, differential calculus provides, as is well known, necessary and sufficient conditions (cf. Section 6.2). In its application to obtaining maximum likelihood estimators, one usually contents oneself, due to the functional form of the functions considered, with the fulfillment of the necessary condition. Depending on the nature of the log-likelihood function, methods of analytic or numerical optimization then suggest themselves. In the following classical examples we use an analytic approach. We proceed in the following steps:

  1. Formulation of the log-likelihood function.
  2. Determination of the first derivative of the log-likelihood function and setting it to zero.
  3. Solving for potential maximum points.

In Theorem 31.1 we show that the maximum likelihood estimator for the parameter of the Bernoulli model from Example 30.2 is given by the corresponding sample mean, and in Theorem 31.2 we show that the maximum likelihood estimators for the expectation and variance parameters of the normal distribution model from Example 30.1 are given by the sample mean and a modified sample variance, respectively.

Examples

Theorem 31.1 (Maximum likelihood estimator of the Bernoulli model) Let \(\mathcal{M}\) be the Bernoulli model, that is, \(y_1,...,y_n \sim \mbox{Bern}(\mu)\). Then \[\begin{equation} \hat{\mu}^{\mbox{\tiny ML}} : \{0,1\}^n \to [0,1], y \mapsto \hat{\mu}^{\mbox{\tiny ML}}(y):= \frac{1}{n}\sum_{i=1}^n y_i \end{equation}\] is a maximum likelihood estimator of \(\mu\).

Proof. We first formulate the log-likelihood function. For the likelihood function, \[\begin{equation} L: ]0,1[ \to ]0,1[, \mu \mapsto L(\mu) := \prod_{i=1}^n \mu^{y_i}(1 - \mu)^{1-y_i} = \mu^{\sum_{i=1}^n y_i}(1 - \mu)^{n - \sum_{i=1}^n y_i}. \end{equation}\] Taking logarithms gives \[\begin{equation} \ell: ]0,1[ \to \mathbb{R}, \mu \mapsto \ell(\mu) = \ln \mu \sum_{i=1}^n y_i + \ln (1- \mu) \left(n - \sum_{i=1}^n y_i \right). \end{equation}\] We then evaluate the derivative of the log-likelihood function. We have \[\begin{align} \begin{split} \frac{d}{d\mu} \ell(\mu) & = \frac{d}{d\mu}\left(\ln \mu \sum_{i=1}^n y_i + \ln (1- \mu) \left(n - \sum_{i=1}^n y_i \right)\right) \\ & = \frac{d}{d\mu} \ln \mu \sum_{i=1}^n y_i + \frac{d}{d\mu} \ln (1 - \mu) \left(n - \sum_{i=1}^n y_i \right) \\ & = \frac{1}{\mu}\sum_{i=1}^n y_i - \frac{1}{1-\mu} \left(n - \sum_{i=1}^n y_i \right). \end{split} \end{align}\] Setting this to zero then gives the following as a necessary condition for a maximum likelihood estimator in the Bernoulli model: \[\begin{equation} \frac{1}{\hat{\mu}^{\mbox{\tiny ML}}}\sum_{i=1}^n y_i - \frac{1}{1-\hat{\mu}^{\mbox{\tiny ML}}} \left(n - \sum_{i=1}^n y_i \right) = 0. \end{equation}\] Solving the maximum likelihood equation for \(\hat{\mu}^{\mbox{\tiny ML}}\) gives \[\begin{align} \begin{split} \frac{1}{\hat{\mu}^{\mbox{\tiny ML}}}\sum_{i=1}^n y_i - \frac{1}{1-\hat{\mu}^{\mbox{\tiny ML}}} \left(n - \sum_{i=1}^n y_i \right) & = 0 \\ \Leftrightarrow \hat{\mu}^{\mbox{\tiny ML}}(1 - \hat{\mu}^{\mbox{\tiny ML}})\left(\frac{1}{\hat{\mu}^{\mbox{\tiny ML}}}\sum_{i=1}^n y_i - \frac{1}{1-\hat{\mu}^{\mbox{\tiny ML}}} \left(n - \sum_{i=1}^n y_i \right) \right) & = 0 \\ \Leftrightarrow \sum_{i=1}^n y_i - \hat{\mu}^{\mbox{\tiny ML}} \sum_{i=1}^n y_i - n \hat{\mu}^{\mbox{\tiny ML}} + \hat{\mu}^{\mbox{\tiny ML}}\sum_{i=1}^n y_i & = 0 \\ \Leftrightarrow n \hat{\mu}^{\mbox{\tiny ML}} & = \sum_{i=1}^n y_i \\ \Leftrightarrow \hat{\mu}^{\mbox{\tiny ML}} & = \frac{1}{n} \sum_{i=1}^n y_i. \\ \end{split} \end{align}\] Thus, \(\hat{\mu}^{\mbox{\tiny ML}} = \frac{1}{n}\sum_{i=1}^n y_i\) is a candidate for a maximum likelihood estimator of \(\mu\). This could be verified by considering the second derivative of \(\ell\), but we refrain from doing so here.

Theorem 31.2 (Maximum likelihood estimators of the normal distribution model) Let \(\mathcal{M}\) be the normal distribution model, that is, \(y_1,...,y_n \sim N\left(\mu,\sigma^2\right)\). Then \[\begin{equation} \hat{\mu}^{\mbox{\tiny ML}} : \mathbb{R}^n \to \mathbb{R}, y \mapsto \hat{\mu}^{\mbox{\tiny ML}}(y) := \frac{1}{n}\sum_{i=1}^n y_i \end{equation}\] and \[\begin{equation} \hat{\sigma}^{2^{\mbox{\tiny ML}}} : \mathbb{R}^n \to \mathbb{R}_{\ge 0}, y \mapsto \hat{\sigma}^{2^{\mbox{\tiny ML}}}(y) := \frac{1}{n}\sum_{i=1}^n \left(y_i - \hat{\mu}^{\mbox{\tiny ML}}\right)^2. \end{equation}\] are maximum likelihood estimators for \(\mu\) and \(\sigma^2\), respectively.

Proof. We first formulate the log-likelihood function. For the likelihood function we obtain \[\begin{align} \begin{split} L: \mathbb{R} \times \mathbb{R}_{>0} \to \mathbb{R}_{>0}, (\mu,\sigma^2) \mapsto L(\mu,\sigma^2) := & \prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left(-\frac{1}{2\sigma^2}(y_i-\mu)^2\right) \\ = & \left(2 \pi \sigma^2\right)^{-\frac{n}{2}}\exp\left(-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2\right). \\ \end{split} \end{align}\] Taking logarithms gives \[\begin{equation} \ell: \mathbb{R} \times \mathbb{R}_{>0} \to \mathbb{R}, (\mu,\sigma^2) \mapsto \mathcal{\ell}_n(\mu,\sigma^2) = -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \sigma^2-\frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2. \end{equation}\] Evaluation of the partial derivatives of the log-likelihood function gives \[\begin{equation} \frac{\partial}{\partial{\mu}} \ell(\mu,\sigma^2) = - \frac{\partial}{\partial{\mu}} \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2 = - \frac{1}{2\sigma^2}\sum_{i=1}^n \frac{\partial}{\partial{\mu}} (y_i-\mu)^2 = \frac{1}{\sigma^2}\sum_{i=1}^n (y_i-\mu) \end{equation}\] and \[\begin{align} \begin{split} \frac{\partial}{\partial\sigma^2} \ell(\mu,\sigma^2) = - \frac{n}{2} \frac{\partial}{\partial\sigma^2} \ln \sigma^2 - \frac{\partial}{\partial\sigma^2} \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2 = - \frac{n}{2 \sigma^2} + \frac{1}{2\sigma^4}\sum_{i=1}^n(y_i-\mu)^2. \end{split} \end{align}\] The system of maximum likelihood equations as an expression of the necessary conditions for extrema of the log-likelihood function therefore has the form \[\begin{equation} \sum_{i=1}^n (y_i-\hat{\mu}^{\mbox{\tiny ML}}) = 0 \mbox{ and } - \frac{n}{2 \hat{\sigma}^{2^{\mbox{\tiny ML}}}} + \frac{1}{2\hat{\sigma}^{4^{\mbox{\tiny ML}}}}\sum_{i=1}^n(y_i-\mu)^2 = 0. \end{equation}\] Solving the system of maximum likelihood equations first gives \[\begin{equation} \sum_{i=1}^n (y_i-\hat{\mu}^{\mbox{\tiny ML}}) = 0 \Leftrightarrow \sum_{i=1}^n y_i = n\hat{\mu}^{\mbox{\tiny ML}} \Leftrightarrow \hat{\mu}^{\mbox{\tiny ML}} = \frac{1}{n}\sum_{i=1}^n y_i. \end{equation}\] Thus, \[\begin{equation} \hat{\mu}^{\mbox{\tiny ML}} = \frac{1}{n}\sum_{i=1}^n y_i \end{equation}\] is a potential maximum likelihood estimator of \(\mu\). Substituting this estimator into the second maximum likelihood equation then gives \[\begin{align} \begin{split} - \frac{n}{2 \hat{\sigma}^{2^{\mbox{\tiny ML}}}} + \frac{1}{2\hat{\sigma}^{4^{\mbox{\tiny ML}}}}\sum_{i=1}^n(y_i-\hat{\mu}^{\mbox{\tiny ML}})^2 & = 0 \\ \Leftrightarrow - n\hat{\sigma}^{2^{\mbox{\tiny ML}}} + \sum_{i=1}^n(y_i-\hat{\mu}^{\mbox{\tiny ML}})^2 & = 0 \\ \Leftrightarrow \hat{\sigma}^{2^{\mbox{\tiny ML}}} & = \frac{1}{n} \sum_{i=1}^n(y_i-\hat{\mu}^{\mbox{\tiny ML}})^2. \end{split} \end{align}\] Thus, \[\begin{equation} \hat{\sigma}^{2^{\mbox{\tiny ML}}} = \frac{1}{n}\sum_{i=1}^n\left(y_i-\hat{\mu}^{\mbox{\tiny ML}}\right)^2 \end{equation}\] is a potential maximum likelihood estimator of \(\sigma^2\). Both potential maximum likelihood estimators can be verified by considering the second derivative of \(\ell\), but we refrain from doing so here.

Note in Theorem 31.2 that \(\hat{\mu}^{\mbox{\tiny ML}}\) is identical to the sample mean \(\bar{y}\), but \(\hat{\sigma}^{2^{\mbox{\tiny ML}}}\) does not coincide with the sample variance \(S^2\). In contrast to the sample variance, the maximum likelihood estimator of \(\sigma^2\) contains the multiplicative factor \(\frac{1}{n}\), not, as in the sample variance, the multiplicative factor \(\frac{1}{n-1}\). We will return to this difference in the context of estimator properties.

Example 31.3 (Application example) To conclude this section, we consider Theorem 31.2 in the context of the application example from Example 30.5. There we based the observed dBDI values on the normal distribution model \[\begin{equation} y_1,...,y_n \sim N(\mu,\sigma^2) \end{equation}\] The maximum likelihood estimators for the parameters of this model can then be evaluated, using Theorem 31.2, with the R sample mean and sample variance functions mean() and var(), and taking into account the identity \[\begin{equation} \frac{n-1}{n}s^2 = \frac{n-1}{n}\cdot\frac{1}{n-1}\sum_{i=1}^n (y_i - \bar{y})^2 = \frac{1}{n}\sum_{i=1}^n \left(y_i - \hat{\mu}^{\mbox{\tiny ML}}\right)^2 = \hat{\sigma}^{2^{\mbox{\tiny ML}}} \end{equation}\] as in the following R code.

D           = read.csv("./_data/bdi-ii-dataset.csv")       # read dataset
y           = D$dBDI                                        # select data
mu_hat      = mean(y)                                       # maximum likelihood estimate of expectation parameter
n           = length(y)                                     # number of data points
sigsqr_hat  = ((n-1)/n)*var(y)                              # maximum likelihood estimate of variance parameter
cat("mu_hat     :", mu_hat,"\nsigsqr_hat :", sigsqr_hat)    # output
mu_hat     : 3.166667 
sigsqr_hat : 12.63889

Based on the principle of maximum likelihood estimation and the available \(n = 12\) data points, \[\begin{equation} \hat{\mu}^{\mbox{\tiny ML}} = 3.17 \mbox{ and } \hat{\sigma}^{2^{\mbox{\tiny ML}}} = 12.6 \end{equation}\] are therefore guesses for the true, but unknown, parameters of the model.

31.2 Estimator properties in finite samples

In general, frequentist estimator properties concern the distribution of estimators as a function of the distribution of the data on which they are based. Because data are random in frequentist inference, estimators are random as well. In particular, observed data values are interpreted as realizations of random variables. Estimators as functions of random variables are therefore also random variables, even though, of course, they take only one concrete value when a concrete dataset is available. We distinguish between estimator properties in finite samples and asymptotic estimator properties. The former are the content of this section and concern the properties of an estimator for a fixed sample size \(n\); the latter are the content of Section 31.3 and concern the properties of an estimator in the limiting case \(n \to \infty\) of large sample sizes.

Let \((\Sigma,S)\) first be a measurable space and \(\hat{\tau} : \mathcal{Y} \to \Sigma\) an estimator of \(\tau : \Theta \to \Sigma\) (cf. Definition 30.3). In what follows, in addition to parameter estimators of the form \[\begin{equation} \tau: \Theta \to \Sigma, \tau(\theta) := \theta \end{equation}\] we will repeatedly consider estimators that, in parametric product models, estimate only functions of the parameters, such as the expectation, variance, or standard deviation of the sample variables. Since, by assumption, the distributions of the sample variables \(y_1,...,y_n\) are identical, these are estimators of the form \[\begin{align} \begin{split} \tau: \Theta \to \Sigma,\, \theta \mapsto \tau(\theta) \mbox{ with } \tau(\theta) := \mathbb{E}_\theta(y_1), \tau(\theta) := \mathbb{V}_\theta(y_1) \mbox{ and } \tau(\theta) := \mathbb{S}_\theta(y_1). \end{split} \end{align}\]

In this section we specifically illuminate four aspects of estimator properties in finite samples. In Section 31.2.1 we first deal with the unbiasedness of an estimator. An estimator is called unbiased if its expectation is identical to the true, but unknown, value \(\tau(\theta)\) for all \(\theta \in \Theta\). In Section 31.2.2 we introduce, with the variance and the standard error of an estimator, two measures of the frequentist variability of estimators. With the mean squared error of an estimator \(\hat{\tau}\) as the expectation of the squared deviation of \(\hat{\tau}(y)\) from \(\tau(\theta)\), in Section 31.2.3 we then introduce an estimator property that allows the accuracy and variability of an estimator to be related in the sense of a so-called bias-variance tradeoff. Finally, the Cramér-Rao inequality discussed in Section 31.2.4 gives a lower bound for the variance of unbiased estimators. An unbiased estimator with variance equal to the lower bound given by the Cramér-Rao inequality has the smallest possible variance of all unbiased estimators and is in this sense an optimal estimator.

31.2.1 Unbiasedness

The concept of unbiasedness of an estimator arises in the context of the error and the bias of an estimator as follows.

Definition 31.4 (Error, bias, and unbiasedness) Let \(y = (y_1,...,y_n)\) be a sample of a frequentist inference model and let \(\hat{\tau}\) be an estimator for \(\tau\).

  • The error of \(\hat{\tau}\) is defined as \[\begin{equation} \hat{\tau}(y) - \tau(\theta). \end{equation}\]
  • The bias of \(\hat{\tau}\) is defined as \[\begin{equation} \mbox{B}(\hat{\tau} ) := \mathbb{E}_{\theta}(\hat{\tau} (y)) - \tau(\theta). \end{equation}\]
  • The estimator \(\hat{\tau}\) is called unbiased if \[\begin{equation} \mbox{B}(\hat{\tau} ) = 0\Leftrightarrow \mathbb{E}_{\theta}(\hat{\tau} (y)) = \tau(\theta) \mbox{ for all } \theta \in \Theta \mbox{ and all } n \in \mathbb{N}. \end{equation}\] Otherwise, \(\hat{\tau}\) is called biased.

Note that in Definition 31.4 the error of an estimator depends on the specific realization of the sample \(y\). The bias, in contrast, is the expected error over sample realizations and is therefore independent of a specific realization. For the special case of a parameter point estimator, Definition 31.4 says that it is unbiased if \[\begin{equation} \mathbb{E}_{\theta}(\hat{\theta}) = \theta. \end{equation}\]

As first examples of unbiased estimators, the following theorem considers the sample mean and the sample variance as estimators of the expectation and the variance of a sample variable.

Theorem 31.3 (Unbiasedness of sample mean and sample variance) Let \(y := (y_1,...,y_n)\) be the sample of a parametric product model. Then:

  1. The sample mean \[\begin{equation} \bar{y} := \frac{1}{n}\sum_{i=1}^n y_i \end{equation}\] is an unbiased estimator of the expectation \(\mathbb{E}_\theta(y_1)\).
  2. The sample variance \[\begin{equation} S^2 := \frac{1}{n-1}\sum_{i=1}^n (y_i - \bar{y})^2 \end{equation}\] is an unbiased estimator of the variance \(\mathbb{V}_\theta(y_1)\).

Proof. (1) The unbiasedness of the sample mean follows with Theorem 23.2 from \[\begin{align} \mathbb{E}_\theta(\bar{y}) = \mathbb{E}_\theta \left(\frac{1}{n}\sum_{i=1}^n y_i \right) = \frac{1}{n}\sum_{i=1}^n \mathbb{E}_\theta\left(y_i \right) = \frac{1}{n}\sum_{i=1}^n \mathbb{E}_\theta\left(y_1 \right) = \frac{1}{n} n \mathbb{E}_\theta\left(y_1 \right) = \mathbb{E}_\theta\left(y_1 \right). \end{align}\]

(2) To show the unbiasedness of the sample variance, we first note that with Theorem 25.6 and independence of the random variables considered, \[\begin{equation} \mathbb{V}_\theta(\bar{y}) = \mathbb{V}_\theta\left(\frac{1}{n} \sum_{i=1}^n y_i \right) = \frac{1}{n^2} \sum_{i=1}^n \mathbb{V}_\theta\left(y_i \right) = \frac{1}{n^2} \sum_{i=1}^n \mathbb{V}_\theta\left(y_1 \right) = \frac{1}{n^2} n \mathbb{V}_\theta\left(y_1 \right) = \frac{\mathbb{V}_\theta\left(y_1 \right)}{n}. \end{equation}\] Furthermore, for the term of summed squared deviations in the sample variance, \[\begin{align} \sum_{i=1}^n \left(y_i - \bar{y}\right)^2 = \sum_{i=1}^n (y_i - \mathbb{E}_\theta(y_1))^2 - n(\bar{y} - \mathbb{E}_\theta(y_1))^2, \end{align}\] because \[\begin{align} \begin{split} \sum_{i=1}^n \left(y_i - \bar{y}\right)^2 & = \sum_{i=1}^n \left(y_i - \mathbb{E}_\theta(y_1) - \bar{y} + \mathbb{E}_\theta(y_1) \right)^2 \\ & = \sum_{i=1}^n \left((y_i - \mathbb{E}_\theta(y_1)) - (\bar{y} - \mathbb{E}_\theta(y_1)) \right)^2 \\ & = \sum_{i=1}^n (y_i-\mathbb{E}_\theta(y_1))^2 - 2(\bar{y}-\mathbb{E}_\theta(y_1))\left(\sum_{i=1}^n(y_i-\mathbb{E}_\theta(y_1))\right) + \sum_{i=1}^n (\bar{y}-\mathbb{E}_\theta(y_1))^2 \\ & = \sum_{i=1}^n (y_i-\mathbb{E}_\theta(y_1))^2 - 2(\bar{y}-\mathbb{E}_\theta(y_1))\left(\sum_{i=1}^ny_i- n\mathbb{E}_\theta(y_1)\right) + n(\bar{y}-\mathbb{E}_\theta(y_1))^2 \\ & = \sum_{i=1}^n (y_i-\mathbb{E}_\theta(y_1))^2 - 2(\bar{y}-\mathbb{E}_\theta(y_1))\left(n\left(\frac{1}{n}\sum_{i=1}^ny_i\right)- n\mathbb{E}_\theta(y_1)\right) + n(\bar{y}-\mathbb{E}_\theta(y_1))^2 \\ & = \sum_{i=1}^n (y_i-\mathbb{E}_\theta(y_1))^2 - 2n(\bar{y}-\mathbb{E}_\theta(y_1))^2 + n(\bar{y}-\mathbb{E}_\theta(y_1))^2 \\ & = \sum_{i=1}^n (y_i - \mathbb{E}_\theta(y_1))^2 - n(\bar{y} - \mathbb{E}_\theta(y_1))^2. \end{split} \end{align}\] Together, this gives \[\begin{align} \mathbb{E}_\theta\left((n-1)S^2\right) & = \mathbb{E}_\theta\left(\sum_{i=1}^n \left(y_i - \bar{y}\right)^2 \right) \\ & = \mathbb{E}_\theta\left(\sum_{i=1}^n (y_i - \mathbb{E}_\theta(y_1))^2 - n(\bar{y} - \mathbb{E}_\theta(y_1))^2 \right) \\ & = \sum_{i=1}^n \mathbb{E}_\theta\left((y_i - \mathbb{E}_\theta(y_1))^2\right) - n \mathbb{E}_\theta\left((\bar{y} - \mathbb{E}_\theta(y_1))^2 \right) \\ & = n \mathbb{V}_\theta(y_1) - n \mathbb{V}_\theta(\bar{y}) \\ & = n \mathbb{V}_\theta(y_1) - n \frac{\mathbb{V}_\theta(y_1)}{n} \\ & = n \mathbb{V}_\theta(y_1) - \mathbb{V}_\theta(y_1) \\ & = (n - 1) \mathbb{V}_\theta(y_1). \end{align}\] Finally, \[\begin{equation} \mathbb{E}_\theta\left(S^2\right) = \mathbb{E}_\theta\left(\frac{1}{n-1}(n-1)S^2 \right) = \frac{1}{n-1}\mathbb{E}_\theta\left((n-1)S^2 \right) = \frac{1}{n-1}(n - 1) \mathbb{V}_\theta(y_1) = \mathbb{V}_\theta(y_1) \end{equation}\] and thus the unbiasedness of the sample variance as an estimator of the variance.

Of course, in Theorem 31.3, due to the identical distribution of the sample variables of a parametric product model, the sample mean and the sample variance are also unbiased estimators of the expectation and variance of any sample variable \(y_i\) with \(1 \le i \le n\). Note that in the proof of the unbiasedness of the sample variance, the denominator \(n-1\) in the definition of the sample variance plays a decisive role.

Although the sample variance is an unbiased estimator of the variance of a sample variable of a parametric product model, this does not hold for the sample standard deviation as an estimator of the standard deviation. This is the content of the following theorem.

Theorem 31.4 (Bias of the sample standard deviation) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model. Then the sample standard deviation \[\begin{equation} S := \sqrt{S^2} \end{equation}\] is a biased estimator of the standard deviation \(\mathbb{S}_\theta(y_1)\).

Proof. We first note that \(\sqrt{\cdot}\) is a strictly concave function and \(\sigma^2 > 0\). With Jensen’s inequality \(\mathbb{E}(f(\xi)) < f(\mathbb{E}(\xi))\) for strictly concave functions (cf. Theorem 26.5), it follows that \[\begin{equation} \mathbb{E}_\theta(S) = \mathbb{E}_\theta\left(\sqrt{S^2}\right) < \sqrt{\mathbb{E}_\theta(S^2)} = \sqrt{\mathbb{V}_\theta(y_1)} = \mathbb{S}_\theta(y_1). \end{equation}\]

In general, nonlinear transformations of unbiased estimators often lead to biased estimators, but we will not deepen this point further here.

Example 31.1 (Simulation example) The following R code illustrates the concepts of unbiasedness and bias of the sample mean, sample variance, and sample standard deviation using the example of a parametric product model with sample distribution \[\begin{equation} y_1,...,y_{12} \sim N(1.7,2) \end{equation}\] The expectations of the estimators are estimated by their sample means across many realizations of \(y_1,...,y_{12}\) as a function of the number of realizations (simulations).

# model formulation
set.seed(0)                          # random number generator
mu      = 1.7                         # true, but unknown, expectation parameter
sigsqr  = 2                           # true, but unknown, variance parameter
n       = 12                          # sample size n
nsim    = 5e4                         # number of simulations
y_bar   = rep(NaN,nsim)               # sample mean array
s_sqr   = rep(NaN,nsim)               # sample variance array
s       = rep(NaN,nsim)               # sample standard deviation array

# simulation iterations
for(sim in 1:nsim){

    # sample realization of y_1,...,y_{12}
    y          = rnorm(n,mu,sqrt(sigsqr))

    # expectation, variance, and standard deviation estimators
    y_bar[sim] = mean(y)              # sample mean
    s_sqr[sim] = var(y)               # sample variance
    s[sim]     = sd(y)                # sample standard deviation
}

# expectation estimation
E_hat_y_bar = cumsum(y_bar)/(1:nsim)  # \mathbb{E}(\bar{y}) estimates
E_hat_s_sqr = cumsum(s_sqr)/(1:nsim)  # \mathbb{E}(S^2) estimates
E_hat_s     = cumsum(s)    /(1:nsim)  # \mathbb{E}(S) estimates

Figure 31.1 visualizes the results of the above simulation. Shown are estimates of the expectations of the sample mean, sample variance, and sample standard deviation as a function of the number of realizations of the sample variables \(y_1,...,y_{12}\), as well as the true, but unknown, values of the expectation, variance, and standard deviation of \(y_1,...,y_{12}\). It is noticeable that these estimates are more variable for a small number of realizations. From an estimate based on about 10,000 realizations of \(y_1,...,y_{12}\) onward, the sample means of \(\bar{y}\) and \(S^2\) correspond to their true, but unknown, values in accordance with their unbiasedness. The sample standard deviation, in contrast, in accordance with its bias, still shows a consistently too low estimate of the true, but unknown, standard deviation even as the number of realizations of \(y_1,...,y_{12}\) continues to increase.

Figure 31.1: Simulation of the unbiasedness of the sample mean and sample variance as estimators of the expectation and variance under normally distributed sample variables, and simulation of the bias of the sample standard deviation as an estimator of the standard deviation under normally distributed sample variables.

31.2.2 Variance and standard error

In the previous section we considered the expectation of an estimator. In this section we consider its variance and its standard deviation. We use the following definition.

Definition 31.5 (Variance and standard error) Let \(y = (y_1,...,y_n)\) be the sample of a frequentist inference model and let \(\hat{\tau}\) be an estimator of \(\tau\).

  • The variance of \(\hat{\tau}\) is defined as \[\begin{equation} \mathbb{V}_\theta(\hat{\tau} ) := \mathbb{E}_\theta \left((\hat{\tau} (y) - \mathbb{E}_\theta(\hat{\tau} (y)))^2\right). \end{equation}\]
  • The standard error of \(\hat{\tau}\) is defined as \[\begin{equation} \mbox{SE}(\hat{\tau} ) := \sqrt{\mathbb{V}_\theta(\hat{\tau})}. \end{equation}\]

The variance of an estimator \(\hat{\tau}\) is thus defined as the variance of the random variable \(\hat{\tau}(y)\). The standard error of an estimator \(\hat{\tau}\) is defined as the standard deviation of \(\hat{\tau}(y)\). As a first example of a standard error, we consider the standard error of the sample mean.

Theorem 31.5 (Standard error of the sample mean) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model. Then the is given by \[\begin{equation} \mbox{SE}(\bar{y}) = \frac{\mathbb{S}_\theta(y_1)}{\sqrt{n}}. \end{equation}\].

Proof. With the variance of the sample mean, we obtain \[\begin{equation} \mbox{SE}(\bar{y}) = \sqrt{\mathbb{V}_\theta(\bar{y})} = \sqrt{\frac{\mathbb{V}_\theta(y_1)}{n}} = \frac{\mathbb{S}_\theta(y_1)}{\sqrt{n}}. \end{equation}\]

The standard error of the mean describes the variability of the sample mean. Since the standard deviation \(\mathbb{S}_\theta(y_1)\) is unknown, the standard error \(\mbox{SE}(\bar{y})\) is also unknown and can therefore only be estimated. With the sample standard deviation as a biased estimator of the standard deviation \(\mathbb{S}_\theta(y_1)\), an also biased estimator for the standard error of the sample mean is \[\begin{equation} \hat{\mbox{SE}}(\bar{y}) = \frac{S}{\sqrt{n}}. \end{equation}\]

As a second example, we consider the standard error of the maximum likelihood estimator for the parameter of a Bernoulli model.

Theorem 31.6 (Standard error of the maximum likelihood estimator of the Bernoulli model parameter) Let \(y = (y_1,...,y_n)\) be the sample of a Bernoulli model and let \(\hat{\mu}^{\mbox{\tiny ML}}\) be the maximum likelihood estimator for the Bernoulli model parameter \(\mu\). Then the standard error of \(\hat{\mu}^{\mbox{\tiny ML}}\) is given by \[\begin{equation} \mbox{SE}\left(\hat{\mu}^{\mbox{\tiny ML}}\right) = \sqrt{\frac{\mu(1-\mu)}{n}}. \end{equation}\]

Proof. We have \[\begin{align} \begin{split} \mbox{SE}\left(\hat{\mu}^{\mbox{\tiny ML}}\right) = \sqrt{\mathbb{V}_\mu\left(\hat{\mu}^{\mbox{\tiny ML}}\right)} = \sqrt{\mathbb{V}_\mu\left(\frac{1}{n}\sum_{i=1}^n y_i \right)} = \sqrt{\frac{1}{n^2}\sum_{i=1}^n \mathbb{V}_\mu(y_i)} = \sqrt{\frac{n \mu(1-\mu)}{n^2}} = \sqrt{\frac{\mu(1-\mu)}{n}}, \end{split} \end{align}\] where the third equality follows from the independence of the \(y_i\) and the fourth equality from the variance \[\begin{equation} \mathbb{V}_\mu(y_1) = \mathbb{V}_\mu(y_i) = \mu(1-\mu) \end{equation}\] of the sample variables.

As in the case of the standard error of the sample mean, the standard error of the maximum likelihood estimator of the Bernoulli model parameter is a true, but unknown, value. An estimator for \(\mbox{SE}\left(\hat{\mu}^{\mbox{\tiny ML}}\right)\) is obtained with the maximum likelihood estimator for the Bernoulli model parameter as \[\begin{equation} \hat{\mbox{SE}}\left(\hat{\mu}^{\mbox{\tiny ML}}\right) = \sqrt{\frac{\hat{\mu}^{\mbox{\tiny ML}}(1-\hat{\mu}^{\mbox{\tiny ML}})}{n}}. \end{equation}\]

Example 31.2 (Simulation example) The following R code simulates the distribution of the maximum likelihood estimator for the parameter of a Bernoulli model with true, but unknown, parameter value \(\mu := 0.4\) for the sample sizes \(n = 20\), \(n = 100\), and \(n = 200\). Figure 31.2 visualizes the resulting distributions using histograms. The variability of the estimates, that is, the width of the histogram distributions, evidently depends on the sample size. Larger sample sizes result in lower variability of the estimator. We will deepen this idea in Section 31.3.

# model formulation
mu       = 0.4                                                      # true, but unknown, parameter value
n_all    = c(20,100,200)                                            # sample sizes n
ns       = 1e4                                                      # number of simulations
mu_hat   = matrix(rep(NaN, length(n_all)*ns), nrow = length(n_all)) # maximum likelihood estimator array

# sample size iterations
for(i in seq_along(n_all)){

    # simulation iterations
    for(s in 1:ns){
        y            = rbinom(n_all[i],1,mu)                        # sample realization of y_1,...,y_n
        mu_hat[i,s]  = mean(y)                                      # sample mean
    }
}
Figure 31.2: Simulation of the distribution of the maximum likelihood estimator of a Bernoulli model. The variability of the estimator evidently depends on the sample size \(n\).

31.2.3 Mean squared error

With the unbiasedness and variance of an estimator, we have encountered in the previous two sections two independent criteria for the quality of estimators. The mean squared error of an estimator introduced in this section enables an integrated consideration of the accuracy (unbiasedness) and variability (variance) of an estimator in the sense of its so-called bias-variance decomposition. We first define the mean squared error of an estimator as follows.

Definition 31.6 (Mean squared error) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model and \(\hat{\tau}\) an estimator for \(\tau\). Then the mean squared error of \(\hat{\tau}\) is defined as \[\begin{equation} \mbox{MSE}(\hat{\tau}) := \mathbb{E}_\theta\left((\hat{\tau}(y) - \tau(\theta))^2\right). \end{equation}\]

The mean squared error of \(\hat{\tau}\) is thus the expected squared deviation of \(\hat{\tau}(y)\) from \(\tau(\theta)\). Note that, by contrast, the variance of \(\hat{\tau}\) is the expected squared deviation of \(\hat{\tau}\) from \(\mathbb{E}_\theta(\hat{\tau}(y))\). As seen in Section 31.2.1, \(\mathbb{E}_\theta(\hat{\tau}(y))\) may coincide with \(\tau(\theta)\), so an estimator may be unbiased, but need not be. If one uses the mean squared error as a quality criterion for an estimator, for example by trying to construct an estimator with the smallest possible mean squared error, then one may accept possible slight deviations from unbiasedness in favor of a small estimator variance. For the mean squared error, the following theorem holds.

Theorem 31.7 (Decomposition of the mean squared error) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model, let \(\hat{\tau}\) be an estimator for \(\tau\), and let \(\mbox{MSE}(\hat{\tau})\) be the mean squared error of \(\hat{\tau}\). Then \[\begin{equation} \mbox{MSE}(\hat{\tau}) = \mbox{B}(\hat{\tau})^2 + \mathbb{V}_\theta(\hat{\tau}). \end{equation}\]

Proof. To simplify notation, let \(\tau := \tau(\theta)\), \(\hat{\tau} := \hat{\tau}(y)\), and \(\bar{\tau}_n := \mathbb{E}_\theta(\hat{\tau}(y))\). Then: \[\begin{align} \begin{split} \mathbb{E}_\theta\left((\hat{\tau} - \tau)^2\right) & = \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n + \bar{\tau}_n - \tau)^2\right) \\ & = \mathbb{E}_\theta \left( (\hat{\tau} - \bar{\tau}_n)^2 + 2(\hat{\tau} - \bar{\tau}_n)(\bar{\tau}_n - \tau) + (\bar{\tau}_n - \tau)^2 \right) \\ & = \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)^2\right) + 2\mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)(\bar{\tau}_n - \tau)\right) + \mathbb{E}_\theta\left((\bar{\tau}_n - \tau)^2\right) \\ & = \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)^2\right) + 2\mathbb{E}_\theta\left( \hat{\tau}\bar{\tau}_n - \hat{\tau}\tau - \bar{\tau}_n\bar{\tau}_n + \bar{\tau}_n\tau \right) + \mathbb{E}_\theta\left((\bar{\tau}_n - \tau)^2\right) \\ & = \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)^2\right) + 2\left( \bar{\tau}_n\bar{\tau}_n - \bar{\tau}_n\tau \right) + \mathbb{E}_\theta\left((\bar{\tau}_n - \tau)^2\right) \\ & = \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)^2\right) + 0 + \mathbb{E}_\theta\left((\bar{\tau}_n - \tau)^2\right) \\ & = \mathbb{E}_\theta\left((\bar{\tau}_n - \tau)^2\right) + \mathbb{E}_\theta\left((\hat{\tau} - \bar{\tau}_n)^2\right) \\ & = \mathbb{E}_\theta\left((\mathbb{E}_\theta(\hat{\tau}) - \tau)^2\right) + \mathbb{E}_\theta\left((\hat{\tau} - \mathbb{E}_\theta(\hat{\tau}))^2\right) \\ & = (\mathbb{E}_\theta(\hat{\tau}) - \tau)^2 + \mathbb{V}_\theta(\hat{\tau}) \\ & = \mbox{B}(\hat{\tau})^2 + \mathbb{V}_\theta(\hat{\tau}). \end{split} \end{align}\]

31.2.4 Cramér-Rao inequality

If several unbiased estimators are available, the estimator with the smallest variance most reliably serves its purpose. But because sample realizations of frequentist inference models are generally variable, the variability of unbiased estimators cannot be arbitrarily small either. The Cramér-Rao inequality gives a lower bound for the variance of unbiased estimators. An unbiased estimator with variance equal to this lower bound therefore has the smallest possible variance of all unbiased estimators and is, in this sense, an optimal estimator.

The Cramér-Rao inequality is based on the concept of the so-called Fisher information, which in turn is based on the concept of the score function of a frequentist inference model. We therefore first introduce these two concepts before formulating and proving the Cramér-Rao inequality.

The results presented generally hold only under a series of mathematical assumptions, the so-called Fisher regularity conditions. For a frequentist inference model with PMF or PDF \(p_\theta\) and parameter space \(\Theta\), these consist in assuming that (1) \(\Theta\) is an open set, so the true, but unknown, parameter value cannot lie on a boundary of the parameter space, (2) the subset of \(\Theta\) on which \(p_\theta\) takes nonzero values does not depend on \(\theta\), (3) the model itself is identifiable, that is, PMFs or PDFs with different parameter values are different functions and therefore imply different sample distributions, (4) the likelihood function of the model is twice continuously differentiable, and (5) for the likelihood function integration and differentiation may be interchanged. We thus assume the Fisher regularity conditions to be fulfilled and consider only models with one-dimensional parameter spaces \(\Theta \subseteq \mathbb{R}\). We first define the terms score function and Fisher information as follows.

Definition 31.7 (Score function and Fisher information) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model with one-dimensional parameter \(\theta \in \Theta \subseteq \mathbb{R}\), and let \(\ell\) be the associated log-likelihood function. Then:

  • The first derivative of \(\ell\) is called the score function of the sample and is denoted by \[\begin{equation} S(\theta) := \frac{d}{d\theta}\ell(\theta) \end{equation}\] For \(n = 1\) we write \(S(\theta) := S_1(\theta)\) and call \(S(\theta)\) the score function of a random variable.
  • The negative second derivative of \(\ell\) is called the Fisher information of the sample and is denoted by \[\begin{equation} I(\theta) := -\frac{d^2}{d\theta^2}\ell(\theta) \end{equation}\] For \(n = 1\) we write \(I(\theta) := I_1(\theta)\) and call \(I(\theta)\) the Fisher information of a random variable.

Because likelihood and log-likelihood functions depend on the realization of a sample, they are random functions against the background of a frequentist inference model. Since Fisher information, as a function of the log-likelihood function, is therefore also a random variable, one must distinguish between the observed and the expected values of the Fisher information.

Definition 31.8 (Observed and expected Fisher information) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model with one-dimensional parameter \(\theta \in \Theta \subseteq \mathbb{R}\), let \(\ell\) be the associated log-likelihood function, and let \(\hat{\theta}^{\mbox{\tiny ML}}\) be a maximum likelihood estimator of \(\theta\). Then:

  • The observed Fisher information of the sample is defined as \[\begin{equation} I\left(\hat{\theta}^{\mbox{\tiny ML}}\right) := -\frac{d^2}{d\theta^2}\ell\left(\hat{\theta}^{\mbox{\tiny ML}}\right), \end{equation}\] that is, the observed Fisher information of the sample is the Fisher information at the point of the maximum likelihood estimator \(\hat{\theta}^{\mbox{\tiny ML}}\).
  • The expected Fisher information of the sample is defined as \[\begin{equation} J(\theta) := \mathbb{E}_\theta(I(\theta)). \end{equation}\] For \(n = 1\) we write \(J(\theta) := J_1(\theta)\) and call \(J(\theta)\) the

Before illustrating these concepts using the Bernoulli model (Theorem 31.10) and the normal distribution model (Theorem 31.12 and Theorem 31.11), we introduce, with the additivity of Fisher information in parametric product models (Theorem 31.8) and the expectation and variance of the score function (Theorem 31.9), important properties of these concepts that simplify the following discussion.

Theorem 31.8 (Additivity of Fisher information) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model with parameter \(\theta \in \Theta \subseteq \mathbb{R}\), let \(\ell\) be the associated log-likelihood function, and let \(I(\theta)\) and \(J(\theta)\) be the Fisher information and the expected Fisher information of the sample, respectively. Then \[\begin{equation} I(\theta) = nI_1(\theta) \mbox{ and } J(\theta) = nJ_1(\theta). \end{equation}\]

Proof. We show the result for the expected Fisher information; the result for the observed Fisher information then holds implicitly. With the linearity of derivatives and expectations, \[\begin{align} \begin{split} J(\theta) & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \ell(\theta)\right) \\ & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \ln \left(\prod_{i=1}^n p_\theta(y_i)\right)\right) \\ & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \sum_{i=1}^n \ln p_\theta(y_i)\right) \\ & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \sum_{i=1}^n \ln p_\theta(y_1)\right) \\ & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \ell_1(\theta)n\right) \\ & = n \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2}\ell_1(\theta))\right) \\ & = n J_1(\theta). \end{split} \end{align}\]

According to Theorem 31.8, to compute the observed or expected Fisher information of a sample in parametric product models, it is sufficient to compute the observed or expected Fisher information of one of the random variables in the sample. Further simplifications in determining Fisher information and justifying the Cramér-Rao inequality follow from the identities formulated in the following theorem.

Theorem 31.9 (Expectation and variance of the score function) The expectation of the score function of a random variable is \[\begin{equation} \mathbb{E}_\theta(S(\theta)) = 0 \end{equation}\] and the variance of the score function of a random variable is \[\begin{equation} \mathbb{V}_\theta(S(\theta)) = J(\theta). \end{equation}\]

Proof. We consider only the case in which \(p_\theta\) is a PDF and first show that \(\mathbb{E}_\theta(S(\theta)) = 0\). \[\begin{align} \begin{split} \mathbb{E}_\theta(S(\theta)) & = \int S(\theta)p_\theta(x) \,dx \\ & = \int \frac{d}{d\theta}\ell(\theta)p_\theta(x) \,dx \\ & = \int \frac{d}{d\theta} \ln L(\theta) p_\theta(x) \,dx \\ & = \int \frac{1}{L(\theta)}\frac{d}{d\theta}L(\theta) p_\theta(x) \,dx \\ & = \int \frac{1}{p_\theta(x)}\frac{d}{d\theta}L(\theta) p_\theta(x) \,dx \\ & = \int \frac{d}{d\theta}L(\theta) \,dx \\ & = \frac{d}{d\theta} \int p_\theta(x)\,dx \\ & = \frac{d}{d\theta} 1 \\ & = 0. \end{split} \end{align}\] With the definition of variance, it follows immediately that \(\mathbb{V}_\theta(S(\theta)) = \mathbb{E}_\theta(S(\theta)^2)\). Next we show that \(J(\theta) = \mathbb{E}_\theta(S(\theta)^2)\) and hence \(\mathbb{V}_\theta(S(\theta)) = J(\theta)\). \[\begin{align} \begin{split} J(\theta) & = \mathbb{E}_\theta\left(-\frac{d^2}{d\theta^2} \ln L(\theta)\right) \\ & = \mathbb{E}_\theta\left(-\frac{d}{d\theta} \frac{\frac{d}{d\theta}L(\theta)}{L(\theta)}\right) \\ & = \mathbb{E}_\theta\left(-\frac{\frac{d^2}{d\theta^2} L(\theta) L(\theta) - \frac{d}{d\theta}L(\theta)\frac{d}{d\theta}L(\theta)}{L(\theta)L(\theta)}\right) \\ & = - \mathbb{E}_\theta\left(\frac{\frac{d^2}{d\theta^2} L(\theta)}{L(\theta)}\right) + \mathbb{E}_\theta\left(\frac{\left(\frac{d}{d\theta}L(\theta)\right)^2}{(L(\theta))^2}\right) \\ & = - \int \frac{\frac{d^2}{d\theta^2} L(\theta)}{L(\theta)} p_\theta(x) \,dx + \int \frac{\left(\frac{d}{d\theta}L(\theta)\right)^2}{(L(\theta))^2} p_\theta(x) \,dx \\ & = - \frac{d^2}{d\theta^2} \int p_\theta(x) \,dx + \int \left(\frac{1}{L(\theta)}\frac{d}{d\theta}L(\theta)\right)^2 p_\theta(x) \,dx \\ & = - \frac{d^2}{d\theta^2} 1 + \int \left(\frac{d}{d\theta} \ln L(\theta) \right)^2 p_\theta(x) \,dx \\ & = \mathbb{E}_\theta\left(S(\theta)^2\right). \end{split} \end{align}\]

The expectation of the derivative of the log-likelihood function is thus always zero, and the expected Fisher information is always equal to the variance of the score function. We now compute the score function and the different forms of Fisher information concretely for the familiar frequentist inference models. The following theorem first summarizes the results for the Bernoulli model.

Theorem 31.10 (Score function and Fisher information of the Bernoulli model) Let \(y = (y_1,...,y_n)\) be the sample of a Bernoulli model with parameter \(\mu \in ]0,1[\). Then:

  • The score function of the sample is \[\begin{equation} S : ]0,1[ \to \mathbb{R}, \mu \mapsto S(\mu) := \frac{1}{\mu}\sum_{i=1}^n y_i - \frac{1}{1-\mu} \left(n - \sum_{i=1}^n y_i \right). \end{equation}\]
  • The Fisher information of the sample is \[\begin{equation} I : ]0,1[ \to \mathbb{R}, \mu \mapsto I(\mu) := I(\mu) = \frac{ny}{\mu^2} + \frac{n(1 - y)}{1-\mu}^{2}. \end{equation}\]
  • The observed Fisher information of the sample is \[\begin{equation} I: ]0,1[ \to \mathbb{R}, \hat{\mu}^{\mbox{\tiny ML}} \mapsto I\left(\hat{\mu}^{\mbox{\tiny ML}}\right) := \frac{ny}{\hat{\mu}_{n}^{{\mbox{\tiny ML}}^2}} + \frac{n(1 - y)}{1-\hat{\mu}^{\mbox{\tiny ML}}}. \end{equation}\]
  • The expected Fisher information of the sample is \[\begin{equation} J: ]0,1[ \to \mathbb{R}, \mu \mapsto J(\mu) := \frac{n}{\mu(1-\mu)}. \end{equation}\]

Proof. The score function was already derived in the context of maximum likelihood estimation of \(\mu\). We consider the Fisher information of a single Bernoulli random variable \(y\). \[\begin{align} \begin{split} I(\mu) & := -\frac{d^2}{d\mu^2} \ell_1(\mu) \\ & = -\frac{d^2}{d\mu^2} \ln p_\mu(y) \\ & = -\frac{d^2}{d\mu^2}\left(y \ln \mu + (1 - y) \ln (1-\mu)\right) \\ & = -\frac{\partial}{\partial{\mu}}\left(\frac{\partial}{\partial{\mu}}\left(y \ln \mu + (1 - y) \ln (1-\mu)\right)\right) \\ & = -\frac{\partial}{\partial{\mu}}\left(\frac{y}{\mu} + \frac{(1 - y)}{1-\mu}\right) \\ & = -\left(-\frac{y}{\mu^2} - \frac{(1 - y)}{1-\mu}^{2}\right) \\ & = \frac{y}{\mu^2} + \frac{(1 - y)}{1-\mu}^{2}. \\ \end{split} \end{align}\] Thus the expected Fisher information of the random variable \(y\) is \[\begin{align} \begin{split} J(\mu) & = \mathbb{E}_\mu(I(\mu)) \\ & = \mathbb{E}_\mu \left(\frac{y}{\mu^2} + \frac{(1 - y)}{1-\mu}^{2} \right) \\ & = \frac{\mathbb{E}_\mu(y)}{\mu^2} + \frac{(1 - \mathbb{E}_\mu(y))}{1-\mu}^{2} \\ & = \frac{\mu}{\mu^2} + \frac{(1 - \mu)}{1-\mu}^{2} \\ & = \frac{1}{\mu(1-\mu)}. \\ \end{split} \end{align}\] With the additivity property of Fisher information and the definition of observed Fisher information, it follows immediately that \[\begin{equation} I(\mu) = \frac{ny}{\mu^2} + \frac{n(1 - y)}{1-\mu}^{2} \mbox{ and } J(\mu) = \frac{n}{\mu(1-\mu)}. \end{equation}\]

We consider the score function and Fisher information of the normal distribution model only under the additional assumption of a known variance parameter (Theorem 31.11) or a known expectation parameter (Theorem 31.12).

Theorem 31.11 (Score function and Fisher information of the normal distribution model with known variance parameter) Let \(y = (y_1,...,y_n)\) be the sample of a normal distribution model and assume that the variance parameter \(\sigma^2\) is known. Then:

  • The score function of the sample is \[\begin{equation} S: \mathbb{R} \to \mathbb{R}, \mu \mapsto S(\mu) := \frac{1}{\sigma^2}\sum_{i=1}^n(y_i-\mu). \end{equation}\]
  • The Fisher information of the sample is \[\begin{equation} I: \mathbb{R} \to \mathbb{R}, \mu \mapsto I(\mu) := \frac{n}{\sigma^2}. \end{equation}\]
  • The observed Fisher information of the sample is \[\begin{equation} I(\hat{\mu}^{\mbox{\tiny ML}}_n) = \frac{n}{\sigma^2}. \end{equation}\]
  • The expected Fisher information of the sample is \[\begin{equation} J : \mathbb{R} \to \mathbb{R}, \mu \mapsto J(\mu) := \frac{n}{\sigma^2}. \end{equation}\]

Proof. Recall that the log-likelihood function of a normal distribution model with known variance parameter \(\sigma^2\) is given by \[\begin{equation} \ell: \mathbb{R} \to \mathbb{R}, \mu \mapsto \ell(\mu) := -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \sigma^2 - \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2 \end{equation}\] Thus the score function is \[\begin{align} \begin{split} S(\mu) = \frac{\partial}{\partial{\mu}}\ell(\mu) = \frac{1}{\sigma^2}\sum_{i=1}^n(y_i-\mu) \end{split} \end{align}\] The Fisher information of the sample is \[\begin{align} \begin{split} I(\mu) = -\frac{d^2}{d\mu^2}\ell(\mu) = -\frac{\partial}{\partial{\mu}}S(\mu) = -\frac{1}{\sigma^2}\frac{\partial}{\partial{\mu}}\left(\sum_{i=1}^n y_i - n\mu \right) = \frac{n}{\sigma^2}. \end{split} \end{align}\] The observed Fisher information is the Fisher information at the point of the maximum likelihood estimator \(\hat{\mu}^{\mbox{\tiny ML}}_n\). The expected Fisher information finally is \[\begin{align} \begin{split} J(\mu) = \mathbb{E}_\mu(I(\mu)) = \mathbb{E}_\mu\left(\frac{n}{\sigma^2}\right) = \frac{n}{\sigma^2}. \end{split} \end{align}\]

Theorem 31.12 (Score function and Fisher information of the normal distribution model with known expectation parameter) Let \(y = (y_1,...,y_n)\) be the sample of a normal distribution model and assume that the expectation parameter \(\mu\) is known. Then:

  • The score function of the sample is given by \[\begin{equation} S: \mathbb{R}_{>0} \to \mathbb{R}, \sigma^2 \mapsto S(\sigma^2) := - \frac{n}{2 \sigma^2} + \frac{1}{2\sigma^4}\sum_{i=1}^n(y_i-\mu)^2 \end{equation}\]
  • The Fisher information of the sample is given by \[\begin{equation} I: \mathbb{R}_{>0} \to \mathbb{R}, \sigma^2 \mapsto I(\sigma^2) := \frac{1}{\sigma^6}\sum_{i=1}^n (y_i - \mu)^2 - \frac{n}{2\sigma^4} \end{equation}\]
  • The observed Fisher information of the sample is given by \[\begin{equation} I(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n) = \frac{n}{2\hat{\sigma}_{\mbox{\tiny ML}}^4} \end{equation}\]
  • The expected Fisher information of the sample is given by \[\begin{equation} J: \mathbb{R}_{>0} \to \mathbb{R}, \sigma^2 \mapsto J(\sigma^2) := \frac{n}{2\sigma^4}. \end{equation}\]

Proof. Recall that the log-likelihood function of the sample of a normal distribution model with known expectation parameter \(\mu\) is given by \[\begin{equation} \ell: \mathbb{R}_{>0} \to \mathbb{R}, \sigma^2 \mapsto \ell(\sigma^2) := -\frac{n}{2} \ln 2\pi - \frac{n}{2} \ln \sigma^2 - \frac{1}{2\sigma^2}\sum_{i=1}^n(y_i-\mu)^2. \end{equation}\] The score function is therefore \[\begin{align} \begin{split} S(\sigma^2) = \frac{\partial}{\partial\sigma^2}\ell(\sigma^2) = - \frac{n}{2 \sigma^2} + \frac{1}{2\sigma^4}\sum_{i=1}^n(y_i-\mu)^2. \end{split} \end{align}\] The Fisher information of the sample is \[\begin{align} \begin{split} I(\sigma^2) & = -\frac{\partial}{\partial\sigma^2}S(\sigma^2) \\ & = - \frac{\partial}{\partial\sigma^2}\left(\frac{n}{2\sigma^4} - \frac{1}{\sigma^6}\sum_{i=1}^n (y_i - \mu)^2\right) \\ & = \frac{1}{\sigma^6}\sum_{i=1}^n (y_i - \mu)^2 - \frac{n}{2\sigma^4}. \end{split} \end{align}\] The observed Fisher information is the Fisher information at the point of the maximum likelihood estimator \(\hat{\sigma}^{2^{\mbox{\tiny ML}}}\), \[\begin{align} \begin{split} I(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n) & = \frac{\sum_{i=1}^n (y_i - \mu)^2}{\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^3} - \frac{n}{2\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} \\ & = \frac{\sum_{i=1}^n (y_i - \mu)^2}{\frac{1}{n^3}\left(\sum_{i=1}^n (y_i - \mu)^2 \right)^3} - \frac{n}{2\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} \\ & = \frac{1}{\frac{1}{n^3}\left(\sum_{i=1}^n (y_i - \mu)^2 \right)^2} - \frac{n}{2\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} \\ & = \frac{n}{\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} - \frac{n}{2\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} \\ & = \frac{n}{2\left(\hat{\sigma}^{2\,\mbox{\tiny ML}}_n \right)^2} \\ & = \frac{n}{2\hat{\sigma}^{4\,\mbox{\tiny ML}}_n}. \end{split} \end{align}\] The expected Fisher information finally is \[\begin{align} \begin{split} J(\sigma^2) & = \mathbb{E}_{\sigma^2}(I(\sigma^2)) \\ & = \mathbb{E}_{\sigma^2}\left(\frac{1}{\sigma^6}\sum_{i=1}^n (y_i - \mu)^2 - \frac{n}{2\sigma^4}\right) \\ & = \frac{1}{\sigma^6}\sum_{i=1}^n \mathbb{E}_{\sigma^2}\left((y_i - \mu)^2 \right) - \frac{n}{2\sigma^4} \\ & = \frac{1}{\sigma^6}\sum_{i=1}^n \sigma^2 - \frac{n}{2\sigma^4} \\ & = \frac{n\sigma^2}{\sigma^6} - \frac{n}{2\sigma^4} \\ & = \frac{n}{\sigma^4} - \frac{n}{2\sigma^4} \\ & = \frac{n}{2\sigma^4}. \\ \end{split} \end{align}\]

With the properties of the score function discussed above, we can now formulate and prove the Cramér-Rao inequality.

Theorem 31.13 (Cramér-Rao inequality) Let a frequentist inference model with one-dimensional parameter \(\theta \in \Theta \subseteq \mathbb{R}\), PMF or PDF \(p_\theta\), be given, and let \(\hat{\tau}\) be an unbiased estimator of \(\tau(\theta)\). Then \[\begin{equation} \mathbb{V}_\theta(\hat{\tau}) \ge \frac{\left(\frac{d}{d\theta}\tau(\theta)\right)^2}{J(\theta)}. \end{equation}\] In particular, for a parameter estimator with \(\tau(\theta) := \theta\) and thus \[\begin{equation} \hat{\tau} = \hat{\theta} \mbox{ and } \left(\frac{d}{d\theta}\tau(\theta)\right)^2 = 1, \end{equation}\] we have \[\begin{equation} \mathbb{V}_\theta(\hat{\theta}) \ge \frac{1}{J(\theta)}. \end{equation}\] The right-hand side of the above inequalities is called the Cramér-Rao bound.

Proof. We first note that for the random variables \(S(\theta)\) and \(\hat{\tau}\), with the correlation inequality (Theorem 26.4) and the identity of \(\mathbb{V}_\theta(S(\theta))\) and \(J(\theta)\) (Theorem 31.9), \[\begin{equation} \frac{\mathbb{C}_\theta(S(\theta), \hat{\tau})^2}{\mathbb{V}_\theta(S(\theta))\mathbb{V}_\theta(\hat{\tau})} \le 1 \Leftrightarrow \mathbb{V}_\theta(\hat{\tau}) \ge \frac{\mathbb{C}_\theta(S(\theta),\hat{\tau})^2}{J(\theta)}. \end{equation}\] With the covariance shift theorem (Theorem 25.2), the fact that the expectation of the score function is always zero (Theorem 31.9), and the assumed unbiasedness of \(\hat{\tau}\), we first obtain \[\begin{align} \begin{split} \mathbb{C}_\theta(S(\theta),\hat{\tau}) & = \mathbb{E}_\theta(S(\theta)\hat{\tau}) - \mathbb{E}_\theta(S(\theta))\mathbb{E}_\theta(\hat{\tau}) \\ & = \mathbb{E}_\theta(S(\theta)\hat{\tau}) \\ & = \int S(\theta)\,\hat{\tau}\,p_\theta(x) \,dx \\ & = \int \frac{d}{d\theta} \ln L(\theta)\,\hat{\tau}\,p_\theta(x) \,dx \\ & = \int \frac{\frac{d}{d\theta} L(\theta)}{L(\theta)}\,\hat{\tau}\,p_\theta(x) \,dx \\ & = \int \frac{\frac{d}{d\theta} L(\theta)}{p_\theta(x)}\,\hat{\tau}\,p_\theta(x) \,dx \\ & = \int \frac{d}{d\theta} L(\theta)\, \hat{\tau} \,dx \\ & = \frac{d}{d\theta} \int L(\theta)\, \hat{\tau} \,dx \\ & = \frac{d}{d\theta} \int \hat{\tau}\, p_\theta(x) \,dx \\ & = \frac{d}{d\theta} \mathbb{E}_\theta(\hat{\tau}) \\ & = \frac{d}{d\theta} \tau(\theta). \end{split} \end{align}\] Thus it follows directly that \[\begin{equation} \mathbb{V}_\theta(\hat{\tau}) \ge \frac{\left(\frac{d}{d\theta}\tau(\theta) \right)^2}{J(\theta)}. \end{equation}\]

For parameter estimators, in particular, the variance of an unbiased parameter estimator \(\hat{\theta}\) is therefore always greater than or equal to the reciprocal expected Fisher information \(J(\theta)\). If even \[\begin{equation} \mathbb{V}_\theta(\hat{\theta}) = \frac{1}{J(\theta)} \end{equation}\] holds, then the variance of the parameter estimator is minimal and the estimator is thus shown to be an optimal estimator in the sense of the Cramér-Rao inequality. We return to this idea in Section 31.4.

31.3 Asymptotic estimator properties

In this section we give a brief introduction to asymptotic statistics (Vaart (1998)). Asymptotic statistics is the area of frequentist inference that deals with the behavior of statistics and estimators at large sample sizes. The methods of asymptotic statistics are used, on the one hand, to study qualitative estimator properties, as here, and, on the other hand, to approximate estimator properties at large sample sizes. Since sample sizes can nowadays be quite large, the methods of asymptotic statistics are therefore well motivated for application and can be used in many ways there.

Continuing Section 31.2, in this section we illuminate four asymptotic estimator properties. To emphasize that in this section the properties of an estimator \(\hat{\tau}\) depend on the sample size \(n\), we write \(\hat{\tau}_n\) for an estimator in this section. In Section 31.3.1 we consider the asymptotic unbiasedness of an estimator. An estimator \(\hat{\tau}_n\) for \(\tau\) is called asymptotically unbiased if the expectation of \(\hat{\tau}_n\) for large sample sizes \(n \to \infty\) is identical to the true, but unknown, value \(\tau(\theta)\). In Section 31.3.2 we introduce the concept of consistency of an estimator. Intuitively, an estimator \(\hat{\tau}_n\) for \(\tau\) is called consistent if, for large sample sizes \(n \to \infty\), the probability that \(\hat{\tau}_n(y)\) deviates from the true, but unknown, value \(\tau(\theta)\) becomes arbitrarily small. For large sample sizes, the distributions of estimators often result in normal distributions. In Section 31.3.3 we introduce a corresponding formalization with the concept of an asymptotic normal distribution. An estimator \(\hat{\tau}_n\) for \(\tau\) is then called asymptotically normally distributed if, for large sample sizes \(n \to \infty\), the distribution of \(\hat{\tau}_n\) is given by a normal distribution. Finally, in Section 31.3.4, with asymptotic efficiency, we consider an optimality criterion for estimators at sample sizes tending to infinity. An estimator \(\hat{\tau}_n\) for \(\tau\) is then called asymptotically efficient if, for large sample sizes \(n \to \infty\), the distribution of \(\hat{\tau}_n\) is given by a normal distribution with expectation parameter \(\tau(\theta)\) and variance parameter equal to the Cramér-Rao bound.

31.3.1 Asymptotic unbiasedness

Asymptotic unbiasedness of an estimator refines the concept of an unbiased estimator as follows.

Definition 31.9 (Asymptotic unbiasedness) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\tau}_n\) be an estimator for \(\tau\). \(\hat{\tau}_n\) is called asymptotically unbiased if \[\begin{equation} \lim_{n\to\infty} \mathbb{E}_\theta(\hat{\tau}_n(y)) = \tau(\theta) \mbox{ for all } \theta \in \Theta. \end{equation}\]

Asymptotically unbiased estimators are therefore unbiased only when the sample size tends to infinity. Unbiased estimators are always also asymptotically unbiased, because their unbiasedness is independent of the sample size. As an example of an estimator that is only asymptotically unbiased, we consider the maximum likelihood estimator of the variance parameter of the normal distribution model.

Theorem 31.14 (Asymptotic unbiasedness of the variance parameter estimator) Let \(y = (y_1,...,y_n)\) be the sample of a normal distribution model with variance parameter \(\sigma^2\). Then the maximum likelihood estimator of \(\sigma^2\), \[\begin{equation} \hat{\sigma}_n^{2^{\mbox{\tiny ML}}} := \frac{1}{n}\sum_{i=1}^n \left(y_i - \bar{y}_n \right)^2, \end{equation}\] is not unbiased, but asymptotically unbiased.

Proof. With the unbiasedness of the sample variance we first obtain (cf. Theorem 31.3) \[\begin{equation} \mathbb{E}_{\mu,\sigma^2}\left(\hat{\sigma}_n^{2^{\mbox{\tiny ML}}} \right) = \mathbb{E}_{\mu,\sigma^2}\left(\frac{1}{n}\sum_{i=1}^n \left(y_i - \bar{y}_n \right)^2 \right) = \frac{1}{n}\mathbb{E}_{\mu,\sigma^2}\left(\sum_{i=1}^n \left(y_i - \bar{y}_n \right)^2 \right) = \frac{n-1}{n}\sigma^2. \end{equation}\] Thus \[\begin{equation} \mathbb{E}_{\mu,\sigma^2}\left(\hat{\sigma}_n^{2^{\mbox{\tiny ML}}} \right) \neq \sigma^2 \end{equation}\] and \(\hat{\sigma}_n^{2^{\mbox{\tiny ML}}}\) is not an unbiased estimator of \(\sigma^2\). However, \[\begin{equation} \frac{n-1}{n} \to 1 \mbox{ for } n \to \infty, \end{equation}\] so that \[\begin{equation} \lim_{n \to \infty} \mathbb{E}_{\mu,\sigma^2}\left(\hat{\sigma}_n^{2^{\mbox{\tiny ML}}}\right) = \lim_{n \to \infty} \frac{n-1}{n}\sigma^2 = \sigma^2 \lim_{n \to \infty} \frac{n-1}{n} = \sigma^2 \end{equation}\] and the maximum likelihood estimator of \(\sigma^2\) is therefore asymptotically unbiased.

Example 31.4 (Simulation example) The following R code simulates the behavior of the maximum likelihood estimator for the variance parameter of a normal distribution model with true, but unknown, parameters \(\mu = 1\) and \(\sigma^2 = 2\) as a function of the sample size.

# model formulation
mu        = 1                                               # true, but unknown, expectation parameter
sigsqr    = 2                                               # true, but unknown, variance parameter
n         = seq(1,100, by = 2)                              # sample sizes
ns        = 1e3                                             # number of simulations per sample size
sigsqr_ml = matrix(rep(NaN, length(n)*ns),ncol = length(n)) # \hat{\sigma^2}^{ML} array

# simulation
for(i in seq_along(n)){                                     # sample size iterations
    for(s in 1:ns){                                         # sample realization iterations
        y               = rnorm(n[i], mu, sqrt(sigsqr))     # sample realization
        sigsqr_ml[s,i]  = ((n[i]-1)/n[i])*var(y)            # \hat{\sigma^2}^{ML}
    }
}
E_sigsqr_ml = colMeans(sigsqr_ml)                           # estimator expectation estimate

We visualize the expectation of the estimator estimated on the basis of the above simulations as a function of sample size in Figure 31.3. For small sample sizes, the estimator clearly underestimates the true, but unknown, variance parameter; for larger sample sizes, it does not.

Figure 31.3: Simulation of the expectation of the maximum likelihood estimator for the variance parameter of a normal distribution model. For small sample sizes, the maximum likelihood estimator systematically underestimates the true, but unknown, variance parameter; for larger sample sizes, the estimator is approximately unbiased.

31.3.2 Consistency

The concept of consistency of an estimator generalizes the weak law of large numbers from sample means and expectations (cf. Theorem 27.1) to arbitrary estimators and true, but unknown, parameter values.

Definition 31.10 (Consistency) Let \(y := (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\tau}_n\) be an estimator of \(\tau\). Then a sequence of estimators \(\hat{\tau}_1, \hat{\tau}_2, ...\) is called a consistent sequence of estimators if, for every arbitrarily small \(\epsilon > 0\) and every \(\theta \in \Theta\), \[\begin{equation} \lim_{n\to \infty} \mathbb{P}_\theta\left(|\hat{\tau}_n(y) - \tau(\theta)| \ge \epsilon \right) = 0. \end{equation}\] If \(\hat{\tau}_1,\hat{\tau}_2,...\) is a consistent sequence of estimators, then \(\hat{\tau}_n\) is called a consistent estimator.

The intuition behind Definition 31.10 corresponds to the law of large numbers. For sample sizes with \(n \to \infty\), the probability that \(\hat{\tau}_n(y)\) is arbitrarily close to \(\tau(\theta)\), or equivalently that \(\hat{\tau}_n(y)\) deviates far from \(\tau(\theta)\), becomes small. However, these properties must hold for all possible true, but unknown, parameter values. As will be shown below, the sample mean is a consistent estimator for the expectation parameter of a normal distribution model.

Example 31.5 Analogously to Section 27.1, we demonstrate the meaning of consistency by simulation for \(\mu = 1\) and \(\sigma^2 = 2\) using the following R code.

# model formulation
mu      = 1                                           # true, but unknown, value of \mu
sigsqr  = 2                                           # true, but unknown, value of \sigma^2
n       = seq(1,1e3,by = 10)                          # sample size n
eps     = c(0.15, 0.10, 0.05)                         # \epsilon values
ne      = length(eps)                                 # number of \epsilon values
nn      = length(n)                                   # number of sample sizes
ns      = 1000                                        # number of simulations
E       = array(rep(NaN,nn*ne*ns),dim = c(nn,ne,ns))  # event indicator array

# simulation
for(e in seq_along(eps)){                             # \epsilon iterations
    for(i in seq_along(n)){                           # n iterations
        for(s in 1:ns){                               # simulation iterations

            # sample realizations
            y   = rnorm(n[i], mu, sqrt(sigsqr))
            if(abs(mean(y) - mu) >= eps[e]){          # |y_bar - \mu)| \ge \epsilon
                E[i,e,s] = 1
            } else {                                  # |y_bar - \mu)| < \epsilon
                E[i,e,s] = 0
            }
        }
    }
}

# estimate of \mathbb{P}(|\hat{\tau}_n(y)-\tau(\theta)| \ge \epsilon)
P_hat       = apply(E, c(1,2), mean)

We visualize the estimate of the probability \(\mathbb{P}(|\hat{\tau}_n(y)-\tau(\theta)| \ge \epsilon)\) for this example in Figure 31.4 as a function of sample size and criterion value \(\epsilon\). For larger values of \(\epsilon\), small sample sizes suffice for a small probability of deviation of the estimator from the true, but unknown, parameter value; for smaller values of \(\epsilon\), larger sample sizes are necessary.

Figure 31.4: Simulation of the consistency of the sample mean as an estimator for the expectation parameter of the normal distribution model.

The criteria for consistency of estimators stated in Theorem 31.15 and Theorem 31.16 simplify the proof of consistency of an estimator with the help of the mean squared error (Definition 31.6).

Theorem 31.15 (Mean squared error criterion for consistency) Let \(y := (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\tau}_n\) be an estimator of \(\tau\). If \[\begin{equation} \lim_{n\to \infty} \mbox{MSE}(\hat{\tau}_n) = 0, \end{equation}\] then \(\hat{\tau}_n\) is a consistent estimator.

Proof. With Chebyshev’s inequality, \[\begin{equation} \mathbb{P}_\theta\left(|\hat{\tau}_n(y) - \tau(\theta)| \ge \epsilon \right) \le \frac{\mathbb{E}_\theta\left((\hat{\tau}_n(y) - \tau(\theta))^2\right)}{\epsilon^2} \end{equation}\] Taking limits gives \[\begin{equation} \lim_{n\to \infty}\mathbb{P}_\theta\left(|\hat{\tau}_n(y) - \tau(\theta)| \ge \epsilon \right) \le \frac{1}{\epsilon^2}\lim_{n\to\infty}\mathbb{E}_\theta\left((\hat{\tau}_n(y) - \tau(\theta))^2\right). \end{equation}\] If therefore \(\lim_{n\to\infty}\mathbb{E}_\theta\left((\hat{\tau}_n(y) - \tau(\theta))^2\right) = 0\), then, with \(\mathbb{P}_\theta(|\hat{\tau}_n(y) - \tau(\theta)| \ge \epsilon)\ge 0\), \[\begin{equation} \lim_{n\to \infty}\mathbb{P}_\theta\left(|\hat{\tau}_n(y) - \tau(\theta)| \ge \epsilon \right) = 0. \end{equation}\] Thus \(\hat{\tau}_n\) is a consistent estimator.

Theorem 31.16 (Bias-variance criterion for consistency) Let \(y := (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\tau}_n\) be an estimator of \(\tau\). If \[\begin{equation} \lim_{n\to \infty} \mbox{B}(\hat{\tau}_n) = 0 \mbox{ and } \lim_{n\to \infty} \mathbb{V}_\theta(\hat{\tau}_n) = 0 \end{equation}\] hold, then \(\hat{\tau}_n\) is a consistent estimator.

Proof. If \(n \to \infty\), then \(\mbox{B}(\hat{\tau}_n) \to 0\), and thus also \(\mbox{B}(\hat{\tau}_n)^2 \to 0\). If for \(n \to \infty\) both \(\mbox{B}(\hat{\tau}_n)^2 \to 0\) and \(\mathbb{V}_\theta(\hat{\tau}_n) \to 0\), then also \(\lim_{n\to \infty} \mbox{MSE}(\hat{\theta}_n) = 0\). Hence, by Theorem 31.15, \(\hat{\tau}_n\) is consistent.

As an application of Theorem 31.16, the following theorem proves consistency of the sample mean as an estimator of the expectation under normality.

Theorem 31.17 (Consistency of the expectation estimator under normality) Let \(y\) be the sample of a normal distribution model. Then \(\bar{y}_n\) is a consistent estimator of \(\mathbb{E}(y_1)\).

Proof. With the unbiasedness of the sample mean as an estimator of the expectation, we first have \[\begin{equation} \lim_{n \to \infty} \mbox{B}(\bar{y}_n) = 0. \end{equation}\] Furthermore, with the variance of the sample mean, \[\begin{equation} \lim_{n \to \infty} \mathbb{V}_\theta(\bar{y}_n) = \lim_{n\to \infty} \frac{1}{n}\mathbb{V}(y_1) = 0. \end{equation}\] By the bias-variance criterion, the consistency of \(\bar{y}_n\) as an estimator of \(\mathbb{E}(y_1)\) follows.

Note that Theorem 31.17 is of course already implied by the weak law of large numbers (cf. Theorem 27.1).

31.3.3 Asymptotic normality

In some cases, the frequentist distribution of an estimator approaches a normal distribution at large sample sizes. Such an estimator is called asymptotically normally distributed.

Definition 31.11 (Asymptotically normally distributed estimator) Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\theta}_n\) be a parameter estimator for \(\theta\). Furthermore, let \[\begin{equation} \tilde{\theta} \sim N(\mu,\sigma^2) \end{equation}\] be a normally distributed random variable with expectation parameter \(\mu\) and variance parameter \(\sigma^2\). If \(\hat{\theta}_n\) converges in distribution to \(\tilde{\theta}\), that is, if for the CDFs \(P_n\) and \(P\) of \(\hat{\theta}_n\) and \(\tilde{\theta}\), respectively, \[\begin{equation} \lim_{n\to \infty} P_n(\hat{\theta}_n) = P(\tilde{\theta}), \end{equation}\] then \(\hat{\theta}_n\) is called asymptotically normally distributed and we write \[\begin{equation} \hat{\theta}_n \stackrel{a}{\sim} N(\mu,\sigma^2). \end{equation}\]

As an example of an asymptotically normally distributed estimator, in Section 31.3.4 we consider the maximum likelihood estimator of the Bernoulli model parameter. It is remarkable that this estimator is asymptotically normally distributed, since the sample variables of the Bernoulli model take only the values zero and one.

31.3.4 Asymptotic efficiency

In some cases, in addition to asymptotic normality of an estimator, statements about the expectation parameter and variance parameter of this asymptotic normal distribution can also be justified by the form of the frequentist model. The concept of an efficient estimator formulates such a special case.

Definition 31.12 Let \(y = (y_1,...,y_n)\) be the sample of a parametric product model and let \(\hat{\theta}_n\) be a parameter estimator for \(\theta\). Furthermore, let \(J(\theta)\) be the expected Fisher information of the sample \(y\). If \[\begin{equation} \hat{\theta}_n \stackrel{a}{\sim} N\left(\theta, J(\theta)^{-1}\right), \end{equation}\] then \(\hat{\theta}_n\) is called .

An asymptotically efficient estimator is evidently always also asymptotically normally distributed and unbiased. The variance of the asymptotic distribution of an estimator is also called the asymptotic variance. For an asymptotically efficient estimator, the asymptotic variance is identical to the Cramér-Rao bound and is therefore minimal in the set of unbiased asymptotically normally distributed estimators.

Example 31.6 As an example of an asymptotically efficient estimator, we consider the maximum likelihood estimator of the Bernoulli model parameter. The following R code simulates the frequentist distribution of this estimator and determines the PDF of its asymptotic distribution.

# model formulation
mu          = 0.4                                                       # true, but unknown, parameter value
n_all       = c(1e1,5e1,1e2)                                            # sample size
ns          = 1e4                                                       # number of sample realizations
mu_hat      = matrix(rep(NaN, length(n_all)*ns), nrow = length(n_all))  # maximum likelihood estimator array
mu_hat_r    = 1e3                                                       # maximum likelihood estimator-space resolution
mu_hat_y    = seq(0,1,len = mu_hat_r)                                   # maximum likelihood estimator space
mu_hat_p    = matrix(rep(NaN, length(n_all)*mu_hat_r), nrow = length(n_all)) # maximum likelihood PDF array

# sample size iterations
for(i in seq_along(n_all)){

    # simulation iterations
    for(s in 1:ns){
        y           = rbinom(n_all[i],1,mu)                             # sample realization
        mu_hat[i,s] = mean(y)                                           # maximum likelihood estimator realization
    }
    mu_hat_p[i,] = dnorm(mu_hat_y, mu, sqrt(mu*(1-mu)/n_all[i]))        # PDF of asymptotic distribution
}
Figure 31.5: Simulation of the asymptotic efficiency of the maximum likelihood parameter estimator for the parameter of a Bernoulli model. With increasing sample size, the distribution of the simulated estimator values shown by a histogram becomes similar to the normal distribution formulated in Definition 31.12.

31.4 Properties of maximum likelihood estimators

The maximum likelihood principle for obtaining estimators for parameters of frequentist inference models is justified by the following theorem on the properties of maximum likelihood estimators.

Theorem 31.18 (Properties of maximum likelihood estimators) Let \(y\) be the sample of a parametric product model and let \(\hat{\theta}_n^{\mbox{\tiny{ML}}}\) be a maximum likelihood estimator for \(\theta\). Then \(\hat{\theta}_n^{\mbox{\tiny{ML}}}\)

  1. is not necessarily unbiased, but
  2. is asymptotically unbiased,
  3. is consistent,
  4. is asymptotically normally distributed, and
  5. is asymptotically efficient.

For a proof of Theorem 31.18 we refer to Held & Sabanés Bové (2014), Section 3.4. Thus, if one uses the maximum likelihood principle to obtain an estimator, the estimator obtained is guaranteed to satisfy the estimator quality criteria summarized in Theorem 31.18.

31.5 Bibliographic remarks

The results presented in this chapter go back in a very substantial way to Fisher (1922). Aldrich (1997) provides historical context.

Study questions

  1. State the definition of a parameter point estimator.
  2. Explain the concept of a parameter point estimator.
  3. State the definition of the likelihood function and the log-likelihood function.
  4. State the definition of a maximum likelihood estimator.
  5. Explain the procedure for obtaining maximum likelihood estimators.
  6. State the theorem on the maximum likelihood estimator of the Bernoulli model parameter.
  7. State the theorem on the maximum likelihood estimators of the normal distribution model parameters.
  8. State the definition of unbiasedness of an estimator.
  9. Explain the concept of unbiasedness of an estimator.
  10. State the definition of the variance and standard error of an estimator.
  11. State the definition of asymptotic unbiasedness of an estimator.
  12. Explain the concept of asymptotic unbiasedness of an estimator.
  13. State the definition of consistency of an estimator.
  14. Explain the concept of consistency of an estimator.
  15. State the first three statements of the theorem on the properties of maximum likelihood estimators.

Study question answers

  1. See Definition 31.1.
  2. A parameter point estimator gives, based on a sample, a guess for the corresponding true, but unknown, parameter.
  3. See Definition 31.2.
  4. See Definition 31.3.
  5. To obtain a maximum likelihood estimator, one first formulates the log-likelihood function and then determines the zeros of its derivative as potential maximum points. In classical examples one usually uses analytic optimization for this; in applications to more complex models, numerical optimization is usually used.
  6. See Theorem 31.1.
  7. See Theorem 31.2.
  8. See Definition 31.4.
  9. An estimator is called unbiased if its expectation is identical to the true, but unknown, value it estimates.
  10. See Definition 31.5.
  11. See Definition 31.9.
  12. An estimator is called asymptotically unbiased if, for sample sizes tending to infinity, its expectation is identical to the true, but unknown, value it estimates.
  13. See Definition 31.10.
  14. An estimator is called consistent if, for large sample sizes, the probability that the estimator value deviates from the true, but unknown, value becomes arbitrarily small.
  15. See Theorem 31.18, statements (1), (2), and (3).
Aldrich, J. (1997). R.A. Fisher and the making of maximum likelihood 1912-1922. Statistical Science, 12(3), 162–176. https://doi.org/10.1214/ss/1030037906
Fisher, R. A. (1922). On the mathematical foundations of theoretical statistics. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 222(594-604), 309–368. https://doi.org/10.1098/rsta.1922.0009
Held, L., & Sabanés Bové, D. (2014). Applied Statistical Inference. Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-37887-4
Vaart, A. W. van der. (1998). Asymptotic statistics. Cambridge University Press.