A Tutorial on Fisher Information

TL;DR

This tutorial explains the concept of Fisher information across frequentist, Bayesian, and minimum description length paradigms, highlighting its roles in hypothesis testing, prior definition, and model complexity measurement.

Contribution

It provides a clear, unified explanation of Fisher information's applications in three different statistical paradigms, enhancing understanding for researchers.

Findings

Clarifies Fisher information's role in hypothesis testing and confidence intervals

Shows how Fisher information defines default priors in Bayesian analysis

Demonstrates Fisher information as a measure of model complexity in MDL

Abstract

In many statistical applications that concern mathematical psychologists, the concept of Fisher information plays an important role. In this tutorial we clarify the concept of Fisher information as it manifests itself across three different statistical paradigms. First, in the frequentist paradigm, Fisher information is used to construct hypothesis tests and confidence intervals using maximum likelihood estimators; second, in the Bayesian paradigm, Fisher information is used to define a default prior; lastly, in the minimum description length paradigm, Fisher information is used to measure model complexity.

Tables

Table 1. Table 1: The description lengths for three observations x obs n = [ y L , y M , y R ] subscript superscript 𝑥 𝑛 obs subscript 𝑦 𝐿 subscript 𝑦 𝑀 subscript 𝑦 𝑅 x^{n}_{\textnormal{obs}}=[y_{L},y_{M},y_{R}] , where y L , y M , y R subscript 𝑦 𝐿 subscript 𝑦 𝑀 subscript 𝑦 𝑅 y_{L},y_{M},y_{R} are the number of observed responses L , M 𝐿 𝑀 L,M and R 𝑅 R respectively.

$x_{\text{obs}}^n=[y_L,y_M,y_R]$	${\text{FIA}}_{{ℳ}_1}\left(x_{\text{obs}}^n\right)$	${\text{FIA}}_{{ℳ}_2}\left(x_{\text{obs}}^n\right)$	Preferred model
$x_{\text{obs},1}^n=[12,1,17]$	42	26	${ℳ}_2$
$x_{\text{obs},2}^n=[14,10,6]$	34	34	tie
$x_{\text{obs},3}^n=[12,16,2]$	29	32	${ℳ}_1$

Equations

$$f(x_i|\theta )={\theta }^{x_i}(1-\theta {)}^{1-x_i},\text{ where }{x}_i\in \mathcal{X}=\{0,1\}.$$

$$f(x^n|\theta )=f(x_1|\theta )\times \dots \times f(x_n|\theta )=\prod _{i=1}^{n}f(x_i|\theta ).$$

$$f(x^n|\theta )=f(x_1|\theta )\times \dots \times f(x_n|\theta )={\theta }^{{\sum }_{i=1}^n{x}_i}(1-\theta {)}^{n-{\sum }_{i=1}^n{x}_i},\text{ where }{x}^n\in {\mathcal{X}}^n.$$

$$f(y|\theta )=(\genfrac{}{}{0ex}{}{n}{y}){\theta }^y(1-\theta {)}^{n-y},\text{ where }y\in \mathcal{Y}=\{0,1,\dots ,n\},$$

$$\displaystyle I_{X}(\theta)=-E\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(X\,|\,\theta)\Big{)}=\begin{cases}-\sum_{x\in\mathcal{X}}\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(x\,|\,\theta)\Big{)}p_{\theta}(x)&\text{if }X\text{ is discrete,}\\ -\int_{\mathcal{X}}\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(x\,|\,\theta)\Big{)}p_{\theta}(x)\textnormal{d}x&\text{if }X\text{ is continuous.}\end{cases}$$

$$\displaystyle I_{X}(\theta)=\begin{cases}\sum_{x\in\mathcal{X}}\Big{(}\tfrac{\textnormal{d}}{\textnormal{d}\theta}\log f(x\,|\,\theta)\Big{)}^{2}p_{\theta}(x)&\text{if }X\text{ is discrete,}\\ \int_{\mathcal{X}}\Big{(}\tfrac{\textnormal{d}}{\textnormal{d}\theta}\log f(x\,|\,\theta)\Big{)}^{2}p_{\theta}(x)\textnormal{d}x&\text{if }X\text{ is continuous.}\end{cases}$$

$$I_{X^n}(\theta )\ge I_T(\theta ),$$

$$I_{X^n}(\theta )=n{I}_X(\theta )=n\frac{1}{\theta (1-\theta )},$$

$$\sqrt{n}(\hat{\theta }-{\theta }^{\ast })\stackrel{D}{\to }\mathcal{N}(0,I_X^{-1}({\theta }^{\ast })),\text{ as n }\to \infty .$$

$$\displaystyle(\hat{\theta}-\theta^{*})\overset{D}{\approx}\mathcal{N}\Big{(}0,1/\big{(}nI_{X}(\theta^{*})\big{)}\Big{)}.$$

$$\left({\theta }^{\ast }-1.96\sqrt{\frac{1}{n}{I}_X^{-1}({\theta }^{\ast })},{\theta }^{\ast }+1.96\sqrt{\frac{1}{n}{I}_X^{-1}({\theta }^{\ast })}\right),$$

$$\displaystyle f(x_{i}\,|\,\theta)=\tfrac{1}{2b}\exp\big{(}-\tfrac{|x_{i}-\theta|}{b}\big{)},$$

$$g(\theta |X^n=x_{\text{obs}}^n)=\frac{f(x_{\text{obs}}^n|\theta )g(\theta )}{{\int }_{\Theta }f(x_{\text{obs}}^n|\theta )g(\theta )\text{d}\theta }.$$

$$\text{posterior}=\frac{\text{likelihood}\times \text{prior}}{\text{marginal likelihood}}.$$

$$\displaystyle P\Big{(}\theta^{*}\in J_{\theta}\Big{)}=\int_{J_{\theta}}g(\theta)\textnormal{d}\theta=\frac{1}{V_{\Theta}}\int_{\theta_{a}}^{\theta_{b}}1\textnormal{d}\theta=\frac{\theta_{b}-\theta_{a}}{V_{\Theta}}.$$

$$g_J(\theta )=\frac{1}{V}\sqrt{I_X(\theta )},\text{ where }V={\int }_{\Theta }\sqrt{I_X(\theta )}\text{d}\theta ,$$

$$g_J(\varphi )=\frac{3{\varphi }^2}{V\sqrt{{\pi }^6-{\varphi }^6}},\text{ where }V=\pi ,$$

$$g_J(\theta )=\frac{1}{V\sqrt{\theta (1-\theta )}},\text{ where }V=\pi .$$

$$V={\int }_{{\mathcal{M}}_{\Theta }}1\text{d}{m}_{\theta }(X)={\int }_{\Theta }\sqrt{I_X(\theta )}\text{d}\theta .$$

$$\text{AIC}=$$

$$\text{BIC}=$$

$$\text{FIA}=$$

$$p_{\text{NML}}(x_{\text{obs}}^n|{\mathcal{M}}_j)=\frac{f_j(x_{\text{obs}}^n|{\hat{\theta }}_j(x_{\text{obs}}^n))}{{\sum }_{x^n\in {\mathcal{X}}^n}{f}_j(x^n|{\hat{\theta }}_j(x^n))},$$

$$\displaystyle-\log p_{\text{NML}}(x^{n}_{\textnormal{obs}}\,|\,\mathcal{M}_{j})=-\log f_{j}\big{(}x^{n}_{\textnormal{obs}}\,|\,\hat{\theta}_{j}(x^{n}_{\textnormal{obs}})\big{)}+\underbrace{\log\sum f_{j}(x^{n}\,|\,\hat{\theta}_{j}(x^{n}))}_{\text{Model complexity}},$$

$$\text{FIA}=$$

$$\displaystyle f_{1}(X\,|\,\vartheta)=\big{[}(1-\vartheta)^{2},2\vartheta(1-\vartheta),\vartheta^{2}\big{]}.$$

$$\displaystyle f_{2}(X\,|\,\alpha)=\big{[}(1-\alpha)/2,\alpha,(1-\alpha)/2\big{]}.$$

$${\hat{\theta }}_1=\hat{\vartheta }=(Y_M+2Y_R)/(2n)\text{ for }{\mathcal{M}}_1\text{, and }{\hat{\theta }}_2=\hat{\alpha }=Y_M/n\text{ for }{\mathcal{M}}_2,$$

$$V_{{\mathcal{M}}_1}$$

$$V_{{\mathcal{M}}_2}$$

$$I_X(\vec{\theta }{)}_{i,j}=$$

$$=$$

Full text

62-01,

62B10,

62F03,

62F12,

62F15,

62B10,

Confidence intervals,

hypothesis testing,

Jeffreys’s prior,

minimum description length,

model complexity,

model selection,

statistical modeling,

\startlocaldefs\endlocaldefs

, , ,

and label=u1, url]www.alexander-ly.com/ label=u2, url]https://jasp-stats.org/

t1This work was supported by the starting grant “Bayes or Bust” awarded by the European Research Council (283876). Correspondence concerning this article may be addressed to Alexander Ly, email address: a.ly@uva.nl. The authors would like to thank Jay Myung, Trisha Van Zandt, and three anonymous reviewers for their comments on an earlier version of this paper. The discussions with Helen Steingroever, Jean-Bernard Salomond, Fabian Dablander, Nishant Mehta, Alexander Etz, Quentin Gronau and Sacha Epskamp led to great improvements of the manuscript. Moreover, the first author is grateful to Chris Klaassen, Bas Kleijn and Henk Pijls for their patience and enthusiasm with which they taught, and answered questions from a not very docile student.

1. 1 Introduction
2. 2 The Role of Fisher Information in Frequentist Statistics
3. 3 The Role of Fisher Information in Bayesian Statistics
4. 4 The Role of Fisher Information in Minimum Description Length
5. 5 Concluding Comments
6. A Generalization to Vector-Valued Parameters: The Fisher Information Matrix
7. B Frequentist Statistics based on Asymptotic Normality
8. C Bayesian use of the Fisher-Rao Metric: The Jeffreys’s Prior
9. D MDL: Coding Theoretical Background
10. E Regularity conditions

1 Introduction

Mathematical psychologists develop and apply quantitative models in order to describe human behavior and understand latent psychological processes. Examples of such models include Stevens’ law of psychophysics that describes the relation between the objective physical intensity of a stimulus and its subjectively experienced intensity (Stevens, 1957); Ratcliff’s diffusion model of decision making that measures the various processes that drive behavior in speeded response time tasks (Ratcliff, 1978); and multinomial processing tree models that decompose performance in memory tasks into the contribution of separate latent mechanisms (Batchelder and Riefer, 1980; Chechile, 1973).

When applying their models to data, mathematical psychologists may operate from within different statistical paradigms and focus on different substantive questions. For instance, working within the classical or frequentist paradigm a researcher may wish to test certain hypotheses or decide upon the number of trials to be presented to participants in order to estimate their latent abilities. Working within the Bayesian paradigm a researcher may wish to know how to determine a suitable default prior on the parameters of a model. Working within the minimum description length (MDL) paradigm a researcher may wish to compare rival models and quantify their complexity. Despite the diversity of these paradigms and purposes, they are connected through the concept of Fisher information.

Fisher information plays a pivotal role throughout statistical modeling, but an accessible introduction for mathematical psychologists is lacking. The goal of this tutorial is to fill this gap and illustrate the use of Fisher information in the three statistical paradigms mentioned above: frequentist, Bayesian, and MDL. This work builds directly upon the Journal of Mathematical Psychology tutorial article by Myung (2003) on maximum likelihood estimation. The intended target group for this tutorial are graduate students and researchers with an affinity for cognitive modeling and mathematical statistics.

To keep this tutorial self-contained we start by describing our notation and key concepts. We then provide the definition of Fisher information and show how it can be calculated. The ensuing sections exemplify the use of Fisher information for different purposes. Section 2 shows how Fisher information can be used in frequentist statistics to construct confidence intervals and hypothesis tests from maximum likelihood estimators (MLEs). Section 3 shows how Fisher information can be used in Bayesian statistics to define a default prior on model parameters. In Section 4 we clarify how Fisher information can be used to measure model complexity within the MDL framework of inference.

1.1 Notation and key concepts

Before defining Fisher information it is necessary to discuss a series of fundamental concepts such as the nature of statistical models, probability mass functions, and statistical independence. Readers familiar with these concepts may safely skip to the next section.

A statistical model is typically defined through a function $f(x_{i}\,|\,\theta)$ that represents how a parameter $\theta$ is functionally related to potential outcomes $x_{i}$ of a random variable $X_{i}$. For ease of exposition, we take $\theta$ to be one-dimensional throughout this text. The generalization to vector-valued $\theta$ can be found in Appendix A, see also Myung and Navarro (2005).

As a concrete example, $\theta$ may represent a participant’s intelligence, $X_{i}$ a participant’s (future) performance on the $i$th item of an IQ test, $x_{i}=1$ the potential outcome of a correct response, and $x_{i}=0$ the potential outcome of an incorrect response on the $i$th item. Similarly, $X_{i}$ is the $i$th trial in a coin flip experiment with two potential outcomes: heads, $x_{i}=1$, or tails, $x_{i}=0$. Thus, we have the binary outcome space $\mathcal{X}=\{0,1\}$. The coin flip model is also known as the Bernoulli distribution $f(x_{i}\,|\,\theta)$ that relates the coin’s propensity $\theta\in(0,1)$ to land heads to the potential outcomes as

[TABLE]

Formally, if $\theta$ is known, fixing it in the functional relationship $f$ yields a function $p_{\theta}(x_{i})=f(x_{i}\,|\,\theta)$ of the potential outcomes $x_{i}$. This $p_{\theta}(x_{i})$ is referred to as a probability density function (pdf) when $X_{i}$ has outcomes in a continuous interval, whereas it is known as a probability mass function (pmf) when $X_{i}$ has discrete outcomes. The pmf $p_{\theta}(x_{i})=P(X_{i}=x_{i}\,|\,\theta)$ can be thought of as a data generative device as it specifies how $\theta$ defines the chance with which $X_{i}$ takes on a potential outcome $x_{i}$. As this holds for any outcome $x_{i}$ of $X_{i}$, we say that $X_{i}$ is distributed according to $p_{\theta}(x_{i})$. For brevity, we do not further distinguish the continuous from the discrete case, and refer to $p_{\theta}(x_{i})$ simply as a pmf.

For example, when the coin’s true propensity is $\theta^{*}=0.3$, replacing $\theta$ by $\theta^{*}$ in the Bernoulli distribution yields the pmf $p_{0.3}(x_{i})=0.3^{x_{i}}0.7^{1-x_{i}}$, a function of all possible outcomes of $X_{i}$. A subsequent replacement $x_{i}=0$ in the pmf $p_{0.3}(0)=0.7$ tells us that this coin generates the outcome [math] with 70% chance.

In general, experiments consist of $n$ trials yielding a potential set of outcomes $x^{n}=(x_{1},\ldots,x_{n})$ of the random vector $X^{n}=(X_{1},\ldots,X_{n})$. These $n$ random variables are typically assumed to be independent and identically distributed (iid). Identically distributed implies that each of these $n$ random variables is governed by one and the same $\theta$, while independence implies that the joint distribution of all these $n$ random variables simultaneously is given by a product, that is,

[TABLE]

As before, when $\theta$ is known, fixing it in this relationship $f(x^{n}\,|\,\theta)$ yields the (joint) pmf of $X^{n}$ as $p_{\theta}(x^{n})=p_{\theta}(x_{1})\times\ldots\times p_{\theta}(x_{n})=\prod_{i=1}^{n}p_{\theta}(x_{i})$.

In psychology the iid assumption is typically evoked when experimental data are analyzed in which participants have been confronted with a sequence of $n$ items of roughly equal difficulty. When the participant can be either correct or incorrect on each trial, the participant’s performance $X^{n}$ can then be related to an $n$-trial coin flip experiment governed by one single $\theta$ over all $n$ trials. The random vector $X^{n}$ has $2^{n}$ potential outcomes $x^{n}$. For instance, when $n=10$, we have $2^{n}=1{,}024$ possible outcomes and we write $\mathcal{X}^{n}$ for the collection of all these potential outcomes. The chance of observing a potential outcome $x^{n}$ is determined by the coin’s propensity $\theta$ as follows

[TABLE]

When the coin’s true propensity $\theta$ is $\theta^{*}=0.6$, replacing $\theta$ by $\theta^{*}$ in Eq. (1.3) yields the joint pmf $p_{0.6}(x^{n})=f(x^{n}\,|\,\theta=0.6)=0.6^{\sum_{i=1}^{n}x_{i}}0.4^{n-\sum_{i=1}^{n}x_{i}}$. The pmf with a particular outcome entered, say, $x^{n}=(1,1,1,1,1,1,1,0,0,0)$ reveals that the coin with $\theta^{*}=0.6$ generates this particular outcome with 0.18% chance.

1.2 Definition of Fisher information

In practice, the true value of $\theta$ is not known and has to be inferred from the observed data. The first step typically entails the creation of a data summary. For example, suppose once more that $X^{n}$ refers to an $n$-trial coin flip experiment and suppose that we observed $x^{n}_{\textnormal{obs}}=(1,0,0,1,1,1,1,0,1,1)$. To simplify matters, we only record the number of heads as $Y=\sum_{i=1}^{n}X_{i}$, which is a function of the data. Applying our function to the specific observations yields the realization $y_{\textnormal{obs}}=Y(x^{n}_{\textnormal{obs}})=7$. Since the coin flips $X^{n}$ are governed by $\theta$, so is a function of $X^{n}$; indeed, $\theta$ relates to the potential outcomes $y$ of $Y$ as follows

[TABLE]

where ${n\choose y}={n!\over y!(n-y)!}$ enumerates the possible sequences of length $n$ that consist of $y$ heads and $n-y$ tails. For instance, when flipping a coin $n=10$ times, there are 120 possible sequences of zeroes and ones that contain $y=7$ heads and $n-y=3$ tails. The distribution $f(y\,|\,\theta)$ is known as the binomial distribution.

The summary statistic $Y$ has $n+1$ possible outcomes, whereas $X^{n}$ has $2^{n}$. For instance, when $n=10$ the statistic $Y$ has only $11$ possible outcomes, whereas $X^{n}$ has $1{,}024$. This reduction results from the fact that the statistic $Y$ ignores the order with which the data are collected. Observe that the conditional probability of the raw data given $Y=y$ is equal to $P(X^{n}\,|\,Y=y,\theta)=1/{n\choose y}$ and that it does not depend on $\theta$. This means that after we observe $Y=y$ the conditional probability of $X^{n}$ is independent of $\theta$, even though each of the distributions of $X^{n}$ and $Y$ separately do depend on $\theta$. We, therefore, conclude that there is no information about $\theta$ left in $X^{n}$ after observing $Y=y$ (Fisher, 1920; Stigler, 1973).

More generally, we call a function of the data, say, $T=t(X^{n})$ a statistic. A statistic is referred to as sufficient for the parameter $\theta$, if the expression $P(X^{n}\,|\,T=t,\theta)$ does not depend on $\theta$ itself. To quantify the amount of information about the parameter $\theta$ in a sufficient statistic $T$ and the raw data, Fisher introduced the following measure.

Definition 1.1** (Fisher information).**

The Fisher information $I_{X}(\theta)$ of a random variable $X$ about $\theta$ is defined as111Under mild regularity conditions Fisher information is equivalently defined as

\displaystyle I_{X}(\theta)=-E\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(X\,|\,\theta)\Big{)}=\begin{cases}-\sum_{x\in\mathcal{X}}\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(x\,|\,\theta)\Big{)}p_{\theta}(x)&\text{if }X\text{ is discrete,}\\ -\int_{\mathcal{X}}\Big{(}\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(x\,|\,\theta)\Big{)}p_{\theta}(x)\textnormal{d}x&\text{if }X\text{ is continuous.}\end{cases}

(1.5)

where $\tfrac{\textnormal{d}^{2}}{\textnormal{d}\theta^{2}}\log f(x\,|\,\theta)$ denotes the second derivate of the logarithm of $f$ with respect $\theta$.

[TABLE]

The derivative $\tfrac{\textnormal{d}}{\textnormal{d}\theta}\log f(x\,|\,\theta)$ is known as the score function, a function of $x$, and describes how sensitive the model (i.e., the functional form $f$) is to changes in $\theta$ at a particular $\theta$. The Fisher information measures the overall sensitivity of the functional relationship $f$ to changes of $\theta$ by weighting the sensitivity at each potential outcome $x$ with respect to the chance defined by $p_{\theta}(x)=f(x\,|\,\theta)$. The weighting with respect to $p_{\theta}(x)$ implies that the Fisher information about $\theta$ is an expectation.

Similarly, Fisher information $I_{X^{n}}(\theta)$ within the random vector $X^{n}$ about $\theta$ is calculated by replacing $f(x\,|\,\theta)$ with $f(x^{n}\,|\,\theta)$, thus, $p_{\theta}(x)$ with $p_{\theta}(x^{n})$ in the definition. Moreover, under the assumption that the random vector $X^{n}$ consists of $n$ iid trials of $X$ it can be shown that $I_{X^{n}}(\theta)=nI_{X}(\theta)$, which is why $I_{X}(\theta)$ is also known as the unit Fisher information.222Note the abuse of notation – we dropped the subscript $i$ for the $i$th random variable $X_{i}$ and denote it simply by $X$ instead. Intuitively, an experiment consisting of $n=10$ trials is expected to be twice as informative about $\theta$ compared to an experiment consisting of only $n=5$ trials. $\hfill\diamond$

Intuitively, we cannot expect an arbitrary summary statistic $T$ to extract more information about $\theta$ than what is already provided by the raw data. Fisher information adheres to this rule, as it can be shown that

[TABLE]

with equality if and only if $T$ is a sufficient statistic for $\theta$.

Example 1.1** (name=The information about $\theta$ within the raw data and a summary statistic).**

A direct calculation with a Bernoulli distributed random vector $X^{n}$ shows that the Fisher information about $\theta$ within an $n$-trial coin flip experiment is given by

[TABLE]

where $I_{X}(\theta)=\tfrac{1}{\theta(1-\theta)}$ is the Fisher information of $\theta$ within a single trial. As shown in Fig. 1, the unit Fisher information $I_{X}(\theta)$ depends on $\theta$.

Similarly, we can calculate the Fisher information about $\theta$ within the summary statistic $Y$ by using the binomial model instead. This yields $I_{Y}(\theta)={n\over\theta(1-\theta)}$. Hence, $I_{X^{n}}(\theta)=I_{Y}(\theta)$ for any value of $\theta$. In other words, the expected information in $Y$ about $\theta$ is the same as the expected information about $\theta$ in $X^{n}$, regardless of the value of $\theta$. $\hfill\diamond$

Observe that the information in the raw data $X^{n}$ and the statistic $Y$ are equal for every $\theta$, and specifically also for its unknown true value $\theta^{*}$. That is, there is no statistical information about $\theta$ lost when we use a sufficient statistic $Y$ instead of the raw data $X^{n}$. This is particular useful when the data set $X^{n}$ is large and can be replaced by single number $Y$.

2 The Role of Fisher Information in Frequentist Statistics

Recall that $\theta$ is unknown in practice and to infer its value we might: (1) provide a best guess in terms of a point estimate; (2) postulate its value and test whether this value aligns with the data, or (3) derive a confidence interval. In the frequentist framework, each of these inferential tools is related to the Fisher information and exploits the data generative interpretation of a pmf. Recall that given a model $f(x^{n}\,|\,\theta)$ and a known $\theta$, we can view the resulting pmf $p_{\theta}(x^{n})$ as a recipe that reveals how $\theta$ defines the chances with which $X^{n}$ takes on the potential outcomes $x^{n}$.

This data generative view is central to Fisher’s conceptualization of the maximum likelihood estimator (MLE; Fisher, 1912; Fisher, 1922; Fisher, 1925; LeCam, 1990; Myung, 2003). For instance, the binomial model implies that a coin with a hypothetical propensity $\theta=0.5$ will generate the outcome $y=7$ heads out of $n=10$ trials with 11.7% chance, whereas a hypothetical propensity of $\theta=0.7$ will generate the same outcome $y=7$ with 26.7% chance. Fisher concluded that an actual observation $y_{\textnormal{obs}}=7$ out of $n=10$ is therefore more likely to be generated from a coin with a hypothetical propensity of $\theta=0.7$ than from a coin with a hypothetical propensity of $\theta=0.5$. Fig. 2 shows that for this specific observation $y_{\textnormal{obs}}=7$, the hypothetical value $\theta=0.7$ is the maximum likelihood estimate; the number $\hat{\theta}_{\textnormal{obs}}=0.7$.

This estimate is a realization of the maximum likelihood estimator (MLE); in this case, the MLE is the function $\hat{\theta}=\tfrac{1}{n}\sum_{i=1}^{n}X_{i}=\tfrac{1}{n}Y$, i.e., the sample mean. Note that the MLE is a statistic, that is, a function of the data.

2.1 Using Fisher information to design an experiment

Since $X^{n}$ depends on $\theta$ so will a function of $X^{n}$, in particular, the MLE $\hat{\theta}$. The distribution of the potential outcomes of the MLE $\hat{\theta}$ is known as the sampling distribution of the estimator and denoted as $f(\hat{\theta}_{\textnormal{obs}}\,|\,\theta)$. As before, when $\theta^{*}$ is assumed to be known, fixing it in $f(\hat{\theta}_{\textnormal{obs}}\,|\,\theta)$ yields the pmf $p_{\theta^{*}}(\hat{\theta}_{\textnormal{obs}})$, a function of the potential outcomes of $\hat{\theta}$. This function $f$ between the parameter $\theta$ and the potential outcomes of the MLE $\hat{\theta}$ is typically hard to describe, but for $n$ large enough it can be characterized by the Fisher information.

For iid data and under general conditions,333Basically, when the Fisher information exists for all parameter values. For details see the advanced accounts provided by Bickel et al. (1993), Hájek (1970), Inagaki (1970), LeCam (1970) and Appendix E. the difference between the true $\theta^{*}$ and the MLE converges in distribution to a normal distribution, that is,

[TABLE]

Hence, for large enough $n$, the “error” is approximately normally distributed444Note that $\hat{\theta}$ is random, while the true value $\theta^{*}$ is fixed. As such, the error $\hat{\theta}-\theta^{*}$ and the rescaled error $\sqrt{n}(\hat{\theta}-\theta^{*})$ are also random. We used $\overset{D}{\rightarrow}$ in Eq. (2.1) to convey that the distribution of the left-hand side goes to the distribution on the right-hand side. Similarly, $\overset{D}{\approx}$ in Eq. (2.2) implies that the distribution of the left-hand side is approximately equal to the distribution given on the right-hand side. Hence, for finite $n$ there will be an error due to using the normal distribution as an approximation to the true sampling distribution. This approximation error is ignored in the constructions given below, see Appendix B.1 for a more thorough discussion.

[TABLE]

This means that the MLE $\hat{\theta}$ generates potential estimates $\hat{\theta}_{\textnormal{obs}}$ around the true value $\theta^{*}$ with a standard error given by the inverse of the square root of the Fisher information at the true value $\theta^{*}$, i.e., $1/\sqrt{nI_{X}(\theta^{*})}$, whenever $n$ is large enough. Note that the chances with which the estimates of $\hat{\theta}$ are generated depend on the true value $\theta^{*}$ and the sample size $n$. Observe that the standard error decreases when the unit information $I_{X}(\theta^{*})$ is high or when $n$ is large. As experimenters we do not have control over the true value $\theta^{*}$, but we can affect the data generating process by choosing the number of trials $n$. Larger values of $n$ increase the amount of information in $X^{n}$, heightening the chances of the MLE producing an estimate $\hat{\theta}_{\textnormal{obs}}$ that is close to the true value $\theta^{*}$. The following example shows how this can be made precise.

Example 2.1** (name=Designing a binomial experiment with the Fisher information, label=designBin).**

Recall that the potential outcomes of a normal distribution fall within one standard error of the population mean with 68% chance. Hence, when we choose $n$ such that $1/\sqrt{nI_{X}(\theta^{*})}=0.1$ we design an experiment that allows the MLE to generate estimates within $0.1$ distance of the true value with 68% chance. To overcome the problem that $\theta^{*}$ is not known, we solve the problem for the worst case scenario. For the Bernoulli model this is given by $\theta=1/2$, the least informative case, see Fig. 1. As such, we have ${1/\sqrt{nI_{X}(\theta^{*})}}\leq{1/\sqrt{nI_{X}(1/2)}}={1/(2\sqrt{n})}=0.1$, where the last equality is the target requirement and is solved by $n=25$.

This leads to the following interpretation. After simulating $k=100$ data sets $x^{n}_{\textnormal{obs},1},\ldots,x^{n}_{\textnormal{obs},k}$ each with $n=25$ trials, we can apply to each of these data sets the MLE yielding $k$ estimates $\hat{\theta}_{\textnormal{obs},1},\ldots,\hat{\theta}_{\textnormal{obs},k}$. The sampling distribution implies that at least $68$ of these $k=100$ estimate are expected to be at most $0.1$ distance away from the true $\theta^{*}$. $\hfill\diamond$

2.2 Using Fisher information to construct a null hypothesis test

The (asymptotic) normal approximation to the sampling distribution of the MLE can also be used to construct a null hypothesis test. When we postulate that the true value equals some hypothesized value of interest, say, $\theta^{*}=\theta_{0}$, a simple plugin then allows us to construct a prediction interval based on our knowledge of the normal distribution. More precisely, the potential outcomes $x^{n}$ with $n$ large enough and generated according to $p_{\theta^{*}}(x^{n})$ leads to potential estimates $\hat{\theta}_{\textnormal{obs}}$ that fall within the range

[TABLE]

with (approximately) 95% chance. This 95%-prediction interval Eq. (2.3) allows us to construct a point null hypothesis test based on a pre-experimental postulate $\theta^{*}=\theta_{0}$.

Example 2.2** (name=A null hypothesis test for a binomial experiment, label=bernoulliH0Ex).**

Under the null hypothesis $H_{0}:\theta^{*}=\theta_{0}=0.5$, we predict that an outcome of the MLE based on $n=10$ trials will lie between $(0.19,0.81)$ with 95% chance. This interval follows from replacing $\theta^{*}$ by $\theta_{0}$ in the 95%-prediction interval Eq. (2.3). The data generative view implies that if we simulate $k=100$ data sets each with the same $\theta^{*}=0.5$ and $n=10$, we would then have $k$ estimates $\hat{\theta}_{\textnormal{obs},1},\ldots,\hat{\theta}_{\textnormal{obs},k}$ of which five are expected to be outside this 95% interval $(0.19,0.81)$. Fisher, therefore, classified an outcome of the MLE that is smaller than 0.19 or larger than 0.81 as extreme under the null and would then reject the postulate $H_{0}:\theta_{0}=0.5$ at a significance level of $.05$. $\hfill\diamond$

The normal approximation to the sampling distribution of the MLE and the resulting null hypothesis test is particularly useful when the exact sampling distribution of the MLE is unavailable or hard to compute.

Example 2.3** (name=An MLE null hypothesis test for the Laplace model, label=laplaceEx).**

Suppose that we have $n$ iid samples from the Laplace distribution

[TABLE]

where $\theta$ denotes the population mean and the population variance is given by $2b^{2}$. It can be shown that the MLE for this model is the sample median, $\hat{\theta}=\hat{M}$, and the unit Fisher information is $I_{X}(\theta)=b^{-2}$. The exact sampling distribution of the MLE is unwieldy (Kotz, Kozubowski and Podgorski, 2001) and not presented here. Asymptotic normality of the MLE is practical, as it allows us to discard the unwieldy exact sampling distribution and, instead, base our inference on a more tractable (approximate) normal distribution with a mean equal to the true value $\theta^{*}$ and a variance equal to ${b^{2}/n}$. For $n=100$, $b=1$ and repeated sampling under the hypothesis $H_{0}:\theta^{*}=\theta_{0}$, approximately 95% of the estimates (the observed sample medians) are expected to fall in the range $(\theta_{0}-0.196,\theta_{0}+0.196)$. $\hfill\diamond$

2.3 Using Fisher information to compute confidence intervals

An alternative to both point estimation and null hypothesis testing is interval estimation. In particular, a 95%-confidence interval can be obtained by replacing in the prediction interval Eq. (2.3) the unknown true value $\theta^{*}$ by an estimate $\hat{\theta}_{\textnormal{obs}}$. Recall that a simulation with $k=100$ data sets each with $n$ trials leads to $\hat{\theta}_{\textnormal{obs},1},\ldots,\hat{\theta}_{\textnormal{obs},k}$ estimates, and each estimate leads to a different 95%-confidence interval. It is then expected that $95$ of these $k=100$ intervals encapsulate the true value $\theta^{*}$.555But see Brown, Cai and DasGupta (2001). Note that these intervals are centred around different points whenever the estimates differ and that their lengths differ, as the Fisher information depends on $\theta$.

Example 2.4** (name=An MLE confidence interval for the Bernoulli model, label=bernoulliCIEx).**

When we observe $y_{\textnormal{obs},1}=7$ heads in $n=10$ trials, the MLE then produces the estimate $\hat{\theta}_{\textnormal{obs},1}=0.7$. Replacing $\theta^{\ast}$ in the prediction interval Eq. (2.3) with $\theta^{*}=\hat{\theta}_{\textnormal{obs},1}$ yields an approximate 95%-confidence interval $(0.42,0.98)$ of length $0.57$. On the other hand, had we instead observed $y_{\textnormal{obs},2}=6$ heads, the MLE would then yield $\hat{\theta}_{\textnormal{obs},2}=0.6$ resulting in the interval $(0.29,0.90)$ of length $0.61$. $\hfill\diamond$

In sum, Fisher information can be used to approximate the sampling distribution of the MLE when $n$ is large enough. Knowledge of the Fisher information can be used to choose $n$ such that the MLE produces an estimate close to the true value, construct a null hypothesis test, and compute confidence intervals.

3 The Role of Fisher Information in Bayesian Statistics

This section outlines how Fisher information can be used to define the Jeffreys’s prior, a default prior commonly used for estimation problems and for nuisance parameters in a Bayesian hypothesis test (e.g., Bayarri et al., 2012; Dawid, 2011; Gronau, Ly and Wagenmakers, 2017; Jeffreys, 1961; Liang et al., 2008; Li and Clyde, 2015; Ly, Verhagen and Wagenmakers, 2016a, b; Ly, Marsman and Wagenmakers, in press; Ly et al., 2017a; Robert, 2016). To illustrate the desirability of the Jeffreys’s prior we first show how the naive use of a uniform prior may have undesirable consequences, as the uniform prior depends on the representation of the inference problem, that is, on how the model is parameterized. This dependence is commonly referred to as lack of invariance: different parameterizations of the same model result in different posteriors and, hence, different conclusions. We visualize the representation problem using simple geometry and show how the geometrical interpretation of Fisher information leads to the Jeffreys’s prior that is parameterization-invariant.

3.1 Bayesian updating

Bayesian analysis centers on the observations $x^{n}_{\textnormal{obs}}$ for which a generative model $f$ is proposed that functionally relates the observed data to an unobserved parameter $\theta$. Given the observations $x^{n}_{\textnormal{obs}}$, the functional relationship $f$ is inverted using Bayes’ rule to infer the relative plausibility of the values of $\theta$. This is done by replacing the potential outcome part $x^{n}$ in $f$ by the actual observations yielding a likelihood function $f(x^{n}_{\textnormal{obs}}\,|\,\theta)$, which is a function of $\theta$. In other words, $x^{n}_{\textnormal{obs}}$ is known, thus, fixed, and the true $\theta$ is unknown, therefore, free to vary. The candidate set of possible values for the true $\theta$ is denoted by $\Theta$ and referred to as the parameter space. Our knowledge about $\theta$ is formalized by a distribution $g(\theta)$ over the parameter space $\Theta$. This distribution is known as the prior on $\theta$, as it is set before any datum is observed. We can use Bayes’ theorem to calculate the posterior distribution over the parameter space $\Theta$ given the data that were actually observed as follows

[TABLE]

This expression is often verbalized as

[TABLE]

The posterior distribution is a combination of what we knew before we saw the data (i.e., the information in the prior), and what we have learned from the observations in terms of the likelihood (e.g., Lee and Wagenmakers, 2013). Note that the integral is now over $\theta$ and not over the potential outcomes.

3.2 Failure of the uniform distribution on the parameter as a noninformative prior

When little is known about the parameter $\theta$ that governs the outcomes of $X^{n}$, it may seem reasonable to express this ignorance with a uniform prior distribution $g(\theta)$, as no parameter value of $\theta$ is then favored over another. This leads to the following type of inference:

Example 3.1** (Uniform prior on $\theta$).**

Before data collection, $\theta$ is assigned a uniform prior, that is, $g(\theta)={1/V_{\Theta}}$ with a normalizing constant of $V_{\Theta}=1$ as shown in the left panel of Fig. 3.

Suppose that we observe coin flip data $x^{n}_{\textnormal{obs}}$ with $y_{\textnormal{obs}}=7$ heads out of $n=10$ trials. To relate these observations to the coin’s propensity $\theta$ we use the Bernoulli distribution as our $f(x^{n}\,|\,\theta)$. A replacement of $x^{n}$ by the data actually observed yields the likelihood function $f(x^{n}_{\textnormal{obs}}\,|\,\theta)=\theta^{7}(1-\theta)^{3}$, which is a function of $\theta$. Bayes’ theorem now allows us to update our prior to the posterior that is plotted in the right panel of Fig. 3. $\hfill\diamond$

Note that a uniform prior on $\theta$ has the length, more generally, volume, of the parameter space as the normalizing constant; in this case, $V_{\Theta}=1$, which equals the length of the interval $\Theta=(0,1)$. Furthermore, a uniform prior can be characterized as the prior that gives equal probability to all sub-intervals of equal length. Thus, the probability of finding the true value $\theta^{*}$ within a sub-interval $J_{\theta}=(\theta_{a},\theta_{b})\subset\Theta=(0,1)$ is given by the relative length of $J_{\theta}$ with respect to the length of the parameter space, that is,

[TABLE]

Hence, before any datum is observed, the uniform prior expresses the belief $P(\theta^{*}\in J_{\theta})=0.20$ of finding the true value $\theta^{*}$ within the interval $J_{\theta}=(0.6,0.8)$. After observing $x^{n}_{\textnormal{obs}}$ with $y_{\textnormal{obs}}=7$ out of $n=10$, this prior is updated to the posterior belief of $P(\theta^{*}\in J_{\theta}\,|\,x^{n}_{\textnormal{obs}})=0.54$, see the shaded areas in Fig. 3.

Although intuitively appealing, it can be unwise to choose the uniform distribution by default, as the results are highly dependent on how the model is parameterized. In what follows, we show how a different parameterization leads to different posteriors and, consequently, different conclusions.

Example 3.2** (Different representations, different conclusions).**

The propensity of a coin landing heads up is related to the angle $\phi$ with which that coin is bent. Suppose that the relation between the angle $\phi$ and the propensity $\theta$ is given by the function \theta=h(\phi)=\frac{1}{2}+\frac{1}{2}\big{(}\tfrac{\phi}{\pi}\big{)}^{3}, chosen here for mathematical convenience.666Another example involves the logit formulation of the Bernoulli model, that is, in terms of $\phi=\log(\tfrac{\theta}{1-\theta})$, where $\Phi=\mathbb{R}$. This logit formulation is the basic building block in item response theory. We did not discuss this example as the uniform prior on the logit cannot be normalized and, therefore, not easily represented in the plots. When $\phi$ is positive the tail side of the coin is bent inwards, which increases the coin’s chances to land heads. As the function $\theta=h(\phi)$ also admits an inverse function $h^{-1}(\theta)=\phi$, we have an equivalent formulation of the problem in Example 3.1, but now described in terms of the angle $\phi$ instead of the propensity $\theta$.

As before, in order to obtain a posterior distribution, Bayes’ theorem requires that we specify a prior distribution. As the problem is formulated in terms of $\phi$, one may believe that a noninformative choice is to assign a uniform prior $\tilde{g}(\phi)$ on $\phi$, as this means that no value of $\phi$ is favored over another. A uniform prior on $\phi$ is in this case given by $\tilde{g}(\phi)=1/V_{\Phi}$ with a normalizing constant $V_{\Phi}=2\pi$, because the parameter $\phi$ takes on values in the interval $\Phi=(-\pi,\pi)$. This uniform distribution expresses the belief that the true $\phi^{*}$ can be found in any of the intervals $(-1.0\pi,-0.8\pi),(-0.8\pi,-0.6\pi),\ldots,(0.8\pi,1.0\pi)$ with 10% probability, because each of these intervals is 10% of the total length, see the top-left panel of Fig. 4.

For the same data as before, the posterior calculated from Bayes’ theorem is given in top-right panel of Fig. 4. As the problem in terms of the angle $\phi$ is equivalent to that of $\theta=h(\phi)$ we can use the function $h$ to translate the posterior in terms of $\phi$ to a posterior on $\theta$, see the bottom-right panel of Fig. 4. This posterior on $\theta$ is noticeably different from the posterior on $\theta$ shown in Figure 3.

Specifically, the uniform prior on $\phi$ corresponds to the prior belief $\tilde{P}(\theta^{*}\in J_{\theta})=0.13$ of finding the true value $\theta^{*}$ within the interval $J_{\theta}=(0.6,0.8)$. After observing $x^{n}_{\textnormal{obs}}$ with $y_{\textnormal{obs}}=7$ out of $n=10$, this prior is updated to the posterior belief of $\tilde{P}(\theta^{*}\in J_{\theta}\,|\,x^{n}_{\textnormal{obs}})=0.29$,777The tilde makes explicit that the prior and posterior are derived from the uniform prior $\tilde{g}(\phi)$ on $\phi$. see the shaded areas in Fig. 4. Crucially, the earlier analysis that assigned a uniform prior to the propensity $\theta$ yielded a posterior probability $P(\theta^{*}\in J_{\theta}\,|\,x^{n}_{\textnormal{obs}})=0.54$, which is markedly different from the current analysis that assigns a uniform prior to the angle $\phi$.

The same posterior on $\theta$ is obtained when the prior on $\phi$ is first translated into a prior on $\theta$ (bottom-left panel) and then updated to a posterior with Bayes’ theorem. Regardless of the stage at which the transformation is applied, the resulting posterior on $\theta$ differs substantially from the result plotted in the right panel of Fig. 3. $\hfill\diamond$

Thus, the uniform prior distribution is not a panacea for the quantification of prior ignorance, as the conclusions depend on how the problem is parameterized. In particular, a uniform prior on the coin’s angle $\tilde{g}(\phi)={1/V_{\Phi}}$ yields a highly informative prior in terms of the coin’s propensity $\theta$. This lack of invariance caused Karl Pearson, Ronald Fisher and Jerzy Neyman to reject 19th century Bayesian statistics that was based on the uniform prior championed by Pierre-Simon Laplace. This rejection resulted in, what is now known as, frequentist statistics, see also Hald (2008), Lehmann (2011), and Stigler (1986).

3.3 A default prior by Jeffreys’s rule

Unlike the other fathers of modern statistical thoughts, Harold Jeffreys continued to study Bayesian statistics based on formal logic and his philosophical convictions of scientific inference (see, e.g., Aldrich, 2005; Etz and Wagenmakers, 2017; Jeffreys, 1961; Ly, Verhagen and Wagenmakers, 2016a, b; Robert, Chopin and Rousseau, 2009; Wrinch and Jeffreys, 1919, 1921, 1923). Jeffreys concluded that the uniform prior is unsuitable as a default prior due to its dependence on the parameterization. As an alternative, Jeffreys (1946) proposed the following prior based on Fisher information

[TABLE]

which is known as the prior derived from Jeffreys’s rule or the Jeffreys’s prior in short. The Jeffreys’s prior is parameterization-invariant, which implies that it leads to the same posteriors regardless of how the model is represented.

Example 3.3** (Jeffreys’s prior).**

The Jeffreys’s prior of the Bernoulli model in terms of $\phi$ is

[TABLE]

which is plotted in the top-left panel of Fig. 5.

The corresponding posterior is plotted in the top-right panel, which we transformed into a posterior in terms of $\theta$ using the function $\theta=h(\phi)$ shown in the bottom-right panel.888The subscript $J$ makes explicit that the prior and posterior are based on the prior derived from Jeffreys’s rule, i.e., $g_{J}(\theta)$ on $\theta$, or equivalently, $g_{J}(\phi)$ on $\phi$.

Similarly, we could have started with the Jeffreys’s prior in terms of $\theta$ instead, that is,

[TABLE]

The Jeffreys’s prior and posterior on $\theta$ are plotted in the bottom-left and the bottom-right panel of Fig. 5, respectively. The Jeffreys’s prior on $\theta$ corresponds to the prior belief $P_{J}(\theta^{*}\in J_{\theta})=0.14$ of finding the true value $\theta^{*}$ within the interval $J_{\theta}=(0.6,0.8)$. After observing $x^{n}_{\textnormal{obs}}$ with $y_{\textnormal{obs}}=7$ out of $n=10$, this prior is updated to the posterior belief of $P_{J}(\theta^{*}\in J_{\theta}\,|\,x^{n}_{\textnormal{obs}})=0.53$, see the shaded areas in Fig. 5. The posterior is identical to the one obtained from the previously described updating procedure that starts with the Jeffreys’s prior on $\phi$ instead of on $\theta$. $\hfill\diamond$

This example shows that the Jeffreys’s prior leads to the same posterior knowledge regardless of how we as researcher represent the problem. Hence, the same conclusions about $\theta$ are drawn regardless of whether we (1) use Jeffreys’s rule to construct a prior on $\theta$ and update with the observed data, or (2) use Jeffreys’s rule to construct a prior on $\phi$, update to a posterior distribution on $\phi$, which is then transformed to a posterior on $\theta$.

3.4 Geometrical properties of Fisher information

In the remainder of this section we make intuitive that the Jeffreys’s prior is in fact uniform in the model space. We elaborate on what is meant by model space and how this can be viewed geometrically. This geometric approach illustrates (1) the role of Fisher information in the definition of the Jeffreys’s prior, (2) the interpretation of the shaded area, and (3) why the normalizing constant is $V=\pi$, regardless of the chosen parameterization.

3.4.1 The model space $\mathcal{M}$

Before we describe the geometry of statistical models, recall that at a pmf can be thought of as a data generating device of $X$, as the pmf specifies the chances with which $X$ takes on the potential outcomes [math] and $1$. Each such pmf has to fulfil two conditions: (i) the chances have to be non-negative, that is, $0\leq p(x)=P(X=x)$ for every possible outcome $x$ of $X$, and (ii) to explicitly convey that there are $w=2$ outcomes, and none more, the chances have to sum to one, that is, $p(0)+p(1)=1$. We call the largest set of functions that adhere to conditions (i) and (ii) the complete set of pmfs $\mathcal{P}$.

As any pmf from $\mathcal{P}$ defines $w=2$ chances, we can represent such a pmf as a vector in $w$ dimensions. To simplify notation, we write $p(X)$ for all $w$ chances simultaneously, hence, $p(X)$ is the vector $p(X)=[p(0),p(1)]$ when $w=2$. The two chances with which a pmf $p(X)$ generates outcomes of $X$ can be simultaneously represented in the plane with $p(0)=P(X=0)$ on the horizontal axis and $p(1)=P(X=1)$ on the vertical axis. In the most extreme case, we have the pmf $p(X)=[1,0]$ or $p(X)=[0,1]$. These two extremes are linked by a straight line in the left panel of Fig. 6.

Any pmf –and the true pmf $p^{*}(X)$ of $X$ in particular– can be uniquely identified with a vector on the line and vice versa. For instance, the pmf $p_{e}(X)=[1/2,1/2]$ (i.e., the two outcomes are generated with the same chance) is depicted as the dot on the line.

This vector representation allows us to associate to each pmf of $X$ a norm, that is, a length. Our intuitive notion of length is based on the Euclidean norm and entails taking the root of the sums of squares. For instance, we can associate to the pmf $p_{e}(X)$ the length $\|p_{e}(X)\|_{2}=\sqrt{(1/2)^{2}+(1/2)^{2}}=1/\sqrt{2}\approx 0.71$. On the other hand, the length of the pmf that states that $X=1$ is generated with 100% chance has length one. Note that by eye, we conclude that $p_{e}(X)$, the arrow pointing to the dot in the left panel in Fig. 6 is indeed much shorter than the arrow pointing to extreme pmf $p(X)=[0,1]$.

This mismatch in lengths can be avoided when we represent each pmf $p(X)$ by two times its square root instead (Kass, 1989), that is, by $m(X)=2\sqrt{p(X)}=[2\sqrt{p(0)},2\sqrt{p(1)}]$.999The factor two is used to avoid a scaling of a quarter, though, its precise value is not essential for the ideas conveyed here. To simplify matters, we also call $m(X)$ a pmf. A pmf that is identified as the vector $m(X)$ is now two units away from the origin, that is, $\|m(X)\|_{2}=\sqrt{m(0)^{2}+m(1)^{2}}=\sqrt{4(p(0)+p(1))}=2$. For instance, the pmf $p_{e}(X)$ is now represented as $m_{e}(X)\approx[1.41,1.41]$. The model space $\mathcal{M}$ is collection of all transformed pmfs and represented as the surface of (the positive part of) a circle, see the right panel of Fig. 6.101010Hence, the model space $\mathcal{M}$ is the collection of all functions on $\mathcal{X}$ such that (i) $m(x)\geq 0$ for every outcome $x$ of $X$, and (ii) $\sqrt{m(0)^{2}+m(1)^{2}}=2$. This vector representation of all the pmfs on $X$ has the advantage that it also induces an inner product, which allows one to project one vector onto another, see Rudin (1991, p. 4), van der Vaart (1998, p. 94) and Appendix E. By representing the set of all possible pmfs of $X$ as vectors $m(X)=2\sqrt{p(X)}$ that reside on the sphere $\mathcal{M}$, we adopted our intuitive notion of distance. As a result, we can now, by simply looking at the figures, clarify that a uniform prior on the parameter space may lead to a very informative prior in the model space $\mathcal{M}$.

3.4.2 Uniform on the parameter space versus uniform on the model space

As $\mathcal{M}$ represents the largest set of pmfs, any model defines a subset of $\mathcal{M}$. Recall that the function $f(x\,|\,\theta)$ represents how we believe a parameter $\theta$ is functionally related to an outcome $x$ of $X$. For each $\theta$ this parameterization yelds a pmf $p_{\theta}(X)$ and, thus, also $m_{\theta}(X)=2\sqrt{p_{\theta}(X)}$. We denote the resulting set of vectors $m_{\theta}(X)$ so created by $\mathcal{M}_{\Theta}$. For instance, the Bernoulli model $f(x\,|\,\theta)=\theta^{x}(1-\theta)^{1-x}$ consists of pmfs given by $p_{\theta}(X)=[f(0\,|\,\theta),f(1\,|\,\theta)]=[1-\theta,\theta]$, which we represent as the vectors $m_{\theta}(X)=[2\sqrt{1-\theta},2\sqrt{\theta}]$. Doing this for every $\theta$ in the parameter space $\Theta$ yields the candidate set of pmfs $\mathcal{M}_{\Theta}$. In this case, we obtain a saturated model, since $\mathcal{M}_{\Theta}=\mathcal{M}$, see the left panel in Fig. 7, where the right most square on the curve corresponds to $m_{0}(X)=[2,0]$. By following the curve in an anti-clockwise manner we encounter squares that represent the pmfs $m_{\theta}(X)$ corresponding to $\theta=0.1,0.2,\ldots,1.0$ respectively.

In the right panel of Fig. 7 the same procedure is repeated, but this time in terms of $\phi$ at $\phi=-1.0\pi,-0.8\pi,\ldots,1.0\pi$. Indeed, filling in the gaps shows that the Bernoulli model in terms of $\theta$ and $\phi$ fully overlap with the largest set of possible pmfs, thus, $\mathcal{M}_{\Theta}=\mathcal{M}=\mathcal{M}_{\Phi}$. Fig. 7 makes precise what is meant when we say that the models $\mathcal{M}_{\Theta}$ and $\mathcal{M}_{\Phi}$ are equivalent; the two models define the same candidate set of pmfs that we believe to be viable data generating devices for $X$.

However, $\theta$ and $\phi$ represent $\mathcal{M}$ in a substantially different manner. As the representation $m(X)=2\sqrt{p(X)}$ respects our natural notion of distance, we conclude, by eye, that a uniform division of $\theta$s with distance, say, $\textnormal{d}\theta=0.1$ does not lead to a uniform partition of the model. More extremely, a uniform division of $\phi$ with distance $\textnormal{d}\phi=0.2\pi$ (10% of the length of the parameter space) also does not lead to a uniform partition of the model. In particular, even though the intervals $(-\pi,-0.8\pi)$ and $(-0.2\pi,0)$ are of equal length in the parameter space $\Phi$, they do not have an equal displacement in the model $\mathcal{M}_{\Phi}$. In effect, the right panel of Fig. 7 shows that the 10% probability that the uniform prior on $\phi$ assigns to $\phi^{*}\in(-\pi,-0.8\pi)$ in parameter space is redistributed over a larger arc length of the model $\mathcal{M}_{\Phi}$ compared to the 10% assigned to $\phi^{*}\in(-0.2\pi,0)$. Thus, a uniform distribution on $\phi$ favors the pmfs $m_{\phi}(X)$ with $\phi$ close to zero. Note that this effect is cancelled by the Jeffreys’s prior, as it puts more mass on the end points compared to $\phi=0$, see the top-left panel of Fig. 5. Similarly, the left panel of Fig. 7 shows that the uniform prior $g(\theta)$ also fails to yield an equiprobable assessment of the pmfs in model space. Again, the Jeffreys’s prior in terms of $\theta$ compensates for the fact that the interval $(0,0.1)$ as compared to $(0.5,0.6)$ in $\Theta$ is more spread out in model space. However, it does so less severely compared to the Jeffreys’s prior on $\phi$. To illustrate, we added additional tick marks on the horizontal axis of the priors in the left panels of Fig. 5. The tick mark at $\phi=-2.8$ and $\theta=0.15$ both indicate the 25% quantiles of their respective Jeffreys’s priors. Hence, the Jeffreys’s prior allocates more mass to the boundaries of $\phi$ than to the boundaries of $\theta$ to compensate for the difference in geometry, see Fig. 7. More generally, the Jeffreys’s prior uses Fisher information to convert the geometry of the model to the parameter space.

Note that because the Jeffreys’s prior is specified using the Fisher information, it takes the functional relationship $f(x\,|\,\theta)$ into account. The functional relationship makes precise how the parameter is linked to the data and, thus, gives meaning and context to the parameter. On the other hand, a prior on $\phi$ specified without taking the functional relationship $f(x\,|\,\phi)$ into account is a prior that neglects the context of the problem. For instance, the right panel of Fig. 7 shows that this neglect with a uniform prior on $\phi$ results in having the geometry of $\Phi=(-\pi,\pi)$ forced onto the model $\mathcal{M}_{\Phi}$.

3.5 Uniform prior on the model

Fig. 7 shows that neither a uniform prior on $\theta$, nor a uniform prior on $\phi$ yields a uniform prior on the model. Alternatively, we can begin with a uniform prior on the model $\mathcal{M}$ and convert this into priors on the parameter spaces $\Theta$ and $\Phi$. This uniform prior on the model translated to the parameters is exactly the Jeffreys’s prior.

Recall that a prior on a space $S$ is uniform, if it has the following two defining features: (i) the prior is proportional to one, and (ii) a normalizing constant given by $V_{S}=\int_{S}1\textnormal{d}s$ that equals the length, more generally, volume of $S$. For instance, a replacement of $s$ by $\phi$ and $S$ by $\Phi=(-\pi,\pi)$ yields the uniform prior on the angles with the normalizing constant $V_{\Phi}=\int_{\Phi}1\textnormal{d}\phi=2\pi$. Similarly, a replacement of $s$ by the pmf $m_{\theta}(X)$ and $S$ by the function space $\mathcal{M}_{\Theta}$ yields a uniform prior on the model $\mathcal{M}_{\Theta}$. The normalizing constant then becomes a daunting looking integral in terms of displacements $\textnormal{d}m_{\theta}(X)$ between functions in model space $\mathcal{M}_{\Theta}$. Fortunately, it can be shown, see Appendix C, that $V$ simplifies to

[TABLE]

Thus, $V$ can be computed in terms of $\theta$ by multiplying the distances $\textnormal{d}\theta$ in $\Theta$ by the root of the Fisher information. Heuristically, this means that the root of the Fisher information translates displacements $\textnormal{d}m_{\theta}(X)$ in the model $\mathcal{M}_{\Theta}$ to distances $\sqrt{I_{X}(\theta)}\textnormal{d}\theta$ in the parameter space $\Theta$.

Recall from Example 3.3 that regardless of the parameterization, the normalizing constant of the Jeffreys’s prior was $\pi$. To verify that this is indeed the length of the model, we use the fact that the circumference of a quarter circle with radius $r=2$ can also be calculated as $V=(2\pi r)/4=\pi$.

Given that the Jeffreys’s prior corresponds to a uniform prior on the model, we deduce that the shaded area in the bottom-left panel of Fig. 5 with $P_{J}(\theta^{*}\in J_{\theta})=0.14$, implies that the model interval J_{m}=\Big{(}m_{0.6}(X),m_{0.8}(X)\Big{)}, the shaded area in the left panel of Fig. 7, accounts for 14% of the model’s length. After updating the Jeffreys’s prior with the observations $x^{n}_{\textnormal{obs}}$ consisting of $y_{\textnormal{obs}}=7$ out of $n=10$ the probability of finding the true data generating pmf $m^{*}(X)$ in this interval of pmfs $J_{m}$ is increased to 53%.

In conclusion, we verified that the Jeffreys’s prior is a prior that leads to the same conclusion regardless of how we parameterize the problem. This parameterization-invariance property is a direct result of shifting our focus from finding the true parameter value within the parameter space to the proper formulation of the estimation problem –as discovering the true data generating pmf $m_{\theta^{*}}(X)=2\sqrt{p_{\theta^{*}}(X)}$ in $\mathcal{M}_{\Theta}$ and by expressing our prior ignorance as a uniform prior on the model $\mathcal{M}_{\Theta}$.

4 The Role of Fisher Information in Minimum Description Length

In this section we graphically show how Fisher information is used as a measure of model complexity and its role in model selection within the minimum description length framework (MDL; de Rooij and Grünwald, 2011; Grünwald, Myung and Pitt, 2005; Grünwald, 2007; Myung, Forster and Browne, 2000; Myung, Navarro and Pitt, 2006; Pitt, Myung and Zhang, 2002).

The primary aim of a model selection procedure is to select a single model from a set of competing models, say, models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$, that best suits the observed data $x^{n}_{\textnormal{obs}}$. Many model selection procedures have been proposed in the literature, but the most popular methods are those based on penalized maximum likelihood criteria, such as the Akaike information criterion (AIC; Akaike, 1974; Burnham and Anderson, 2002), the Bayesian information criterion (BIC; Raftery, 1995; Schwarz, 1978), and the Fisher information approximation (FIA; Grünwald, 2007; Rissanen, 1996). These criteria are defined as follows

[TABLE]

where $n$ denotes the sample size, $d_{j}$ the number of free parameters, $\hat{\theta}_{j}$ the MLE, $I_{\mathcal{M}_{j}}(\theta_{j})$ the unit Fisher information, and $f_{j}$ the functional relationship between the potential outcome $x^{n}$ and the parameters $\theta_{j}$ within model $\mathcal{M}_{j}$.111111For vector-valued parameters $\theta_{j}$, we have a Fisher information matrix and $\text{det }I_{\mathcal{M}_{j}}(\theta_{j})$ refers to the determinant of this matrix. This determinant is always non-negative, because the Fisher information matrix is always a positive semidefinite symmetric matrix. Intuitively, volumes and areas cannot be negative (Appendix C.3.3). Hence, except for the observations $x^{n}_{\textnormal{obs}}$, all quantities in the formulas depend on the model $\mathcal{M}_{j}$. We made this explicit using a subscript $j$ to indicate that the quantity, say, $\theta_{j}$ belongs to model $\mathcal{M}_{j}$.121212For the sake of clarity, we will use different notations for the parameters within the different models. We introduce two models in this section: the model $\mathcal{M}_{1}$ with parameter $\theta_{1}=\vartheta$ which we pit against the model $\mathcal{M}_{2}$ with parameter $\theta_{2}=\alpha$. For all three criteria, the model yielding the lowest criterion value is perceived as the model that generalizes best (Myung and Pitt, in press).

Each of the three model selection criteria tries to strike a balance between model fit and model complexity. Model fit is expressed by the goodness-of-fit terms, which involves replacing the potential outcomes $x^{n}$ and the unknown parameter $\theta_{j}$ of the functional relationships $f_{j}$ by the actually observed data $x^{n}_{\textnormal{obs}}$, as in the Bayesian setting, and the maximum likelihood estimate $\hat{\theta}_{j}(x^{n}_{\textnormal{obs}})$, as in the frequentist setting.

The positive terms in the criteria account for model complexity. A penalization of model complexity is necessary, because the support in the data cannot be assessed by solely considering goodness-of-fit, as the ability to fit observations increases with model complexity (e.g., Roberts and Pashler, 2000). As a result, the more complex model necessarily leads to better fits but may in fact overfit the data. The overly complex model then captures idiosyncratic noise rather than general structure, resulting in poor model generalizability (Myung, Forster and Browne, 2000; Wagenmakers and Waldorp, 2006).

The focus in this section is to make intuitive how FIA acknowledges the trade-off between goodness-of-fit and model complexity in a principled manner by graphically illustrating this model selection procedure, see also Balasubramanian (1996), Kass (1989), Myung, Balasubramanian and Pitt (2000), and Rissanen (1996). We exemplify the concepts with simple multinomial processing tree (MPT) models (e.g., Batchelder and Riefer, 1999; Klauer and Kellen, 2011; Wu, Myung and Batchelder, 2010). For a more detailed treatment of the subject we refer to Appendix D, de Rooij and Grünwald (2011), Grünwald (2007), Myung, Navarro and Pitt (2006), and the references therein.

4.0.1 The description length of a model

Recall that each model specifies a functional relationship $f_{j}$ between the potential outcomes of $X$ and the parameters $\theta_{j}$. This $f_{j}$ is used to define a so-called normalized maximum likelihood (NML) code. For the $j$th model its NML code is defined as

[TABLE]

where the sum in the denominator is over all possible outcomes $x^{n}$ in $\mathcal{X}^{n}$, and where $\hat{\theta}_{j}$ refers to the MLE within model $\mathcal{M}_{j}$. The NML code is a relative goodness-of-fit measure, as it compares the observed goodness-of-fit term against the sum of all possible goodness-of-fit terms. Note that the actual observations $x^{n}_{\textnormal{obs}}$ only affect the numerator, by a plugin of $x^{n}_{\textnormal{obs}}$ and its associated maximum likelihood estimate $\hat{\theta}(x^{n}_{\textnormal{obs}})$ into the functional relationship $f_{j}$ belonging to model $\mathcal{M}_{j}$. The sum in the denominator consists of the same plugins, but for every possible realization of $X^{n}$.131313As before, for continuous data, the sum is replaced by an integral. Hence, the denominator can be interpreted as a measure of the model’s collective goodness-of-fit or the model’s fit capacity. Consequently, for every set of observations $x^{n}_{\textnormal{obs}}$, the NML code outputs a number between zero and one that can be transformed into a non-negative number by taking the negative logarithm as141414Quite deceivingly the minus sign actually makes this definition positive, as $-\log(y)=\log(1/y)\geq 0$ if $0\leq y\leq 1$.

[TABLE]

which is called the description length of model $\mathcal{M}_{j}$. Within the MDL framework, the model with the shortest description length is the model that best describes the observed data $x^{n}_{\textnormal{obs}}$.

The model complexity term is typically hard to compute, but Rissanen (1996) showed that it can be well-approximated by the dimensionality and the geometrical complexity terms. That is,

[TABLE]

is an approximation of the description length of model $\mathcal{M}_{j}$. The determinant is simply the absolute value when the number of free parameters $d_{j}$ is equal to one. Furthermore, the integral in the geometrical complexity term coincides with the normalizing constant of the Jeffreys’s prior, which represented the volume of the model. In other words, a model’s fit capacity is proportional to its volume in model space as one would expect.

In sum, within the MDL philosophy, a model is selected if it yields the shortest description length, as this model uses the functional relationship $f_{j}$ that best extracts the regularities from $x^{n}_{\textnormal{obs}}$. As the description length is often hard to compute, we approximate it with FIA instead (Heck, Moshagen and Erdfelder, 2014). To do so, we have to characterize (1) all possible outcomes of $X$, (2) propose at least two models which will be pitted against each other, and (3) identify the model characteristics: the MLE $\hat{\theta}_{j}$ corresponding to $\mathcal{M}_{j}$, and its volume $V_{\mathcal{M}_{j}}$. In the remainder of this section we show that FIA selects the model that is closest to the data with an additional penalty for model complexity.

4.1 A new running example and the geometry of a random variable with $w=3$ outcomes

To graphically illustrate the model selection procedure underlying MDL we introduce a random variable $X$ that has $w=3$ number of potential outcomes.

Example 4.1** (A psychological task with three outcomes).**

In the training phase of a source-memory task, the participant is presented with two lists of words on a computer screen. List $\mathcal{L}$ is projected on the left-hand side and list $\mathcal{R}$ is projected on the right-hand side. In the test phase, the participant is presented with two words, side by side, that can stem from either list, thus, $ll,lr,rl,rr$. At each trial, the participant is asked to categorize these pairs as either:

• •

$L$* meaning both words come from the left list, i.e., $``ll"$,*

• •

$M$* meaning the words are mixed, i.e., $``lr"$ or $``rl"$,*

• •

$R$* meaning both words come from the right list, i.e., $``rr"$.*

For simplicity we assume that the participant will be presented with $n$ test pairs $X^{n}$ of equal difficulty. $\hfill\diamond$

For the graphical illustration of this new running example, we generalize the ideas presented in Section 3.4.1 from $w=2$ to $w=3$. Recall that a pmf of $X$ with $w$ number of outcomes can be written as a $w$-dimensional vector. For the task described above we know that a data generating pmf defines the three chances $p(X)=[p(L),p(M),p(R)]$ with which $X$ generates the outcomes $[L,M,R]$ respectively.151515As before we write $p(X)=[p(L),p(M),p(R)]$ with a capital $X$ to denote all the $w$ number of chances simultaneously and we used the shorthand notation $p(L)=p(X=L)$, $p(M)=p(X=M)$ and $p(R)=p(X=R)$. As chances cannot be negative, (i) we require that $0\leq p(x)=P(X=x)$ for every outcome $x$ in $\mathcal{X}$, and (ii) to explicitly convey that there are $w=3$ outcomes, and none more, these $w=3$ chances have to sum to one, that is, $\sum_{x\in\mathcal{X}}p(x)=1$. We call the largest set of functions that adhere to conditions (i) and (ii) the complete set of pmfs $\mathcal{P}$. The three chances with which a pmf $p(X)$ generates outcomes of $X$ can be simultaneously represented in three-dimensional space with $p(L)=P(X=L)$ on the left most axis, $p(M)=P(X=M)$ on the right most axis and $p(R)=P(X=R)$ on the vertical axis as shown in the left panel of Fig. 8.161616This is the three-dimensional generalization of Fig. 6.

In the most extreme case, we have the pmf $p(X)=[1,0,0]$, $p(X)=[0,1,0]$ or $p(X)=[0,0,1]$, which correspond to the corners of the triangle indicated by $pL,pM$ and $pR$, respectively. These three extremes are linked by a triangular plane in the left panel of Fig. 8. Any pmf –and the true pmf $p^{*}(X)$ in particular– can be uniquely identified with a vector on the triangular plane and vice versa. For instance, a possible true pmf of $X$ is $p_{e}(X)=[1/3,1/3,1/3]$ (i.e., the outcomes $L,M$ and $R$ are generated with the same chance) depicted as a (red) dot on the simplex.

This vector representation allows us to associate to each pmf of $X$ the Euclidean norm. For instance, the representation in the left panel of Fig. 8 leads to an extreme pmf $p(X)=[1,0,0]$ that is one unit long, while $p_{e}(X)=[1/3,1/3,1/3]$ is only $\sqrt{(1/3)^{2}+(1/3)^{2}+(1/3)^{2}}\approx 0.58$ units away from the origin. As before, we can avoid this mismatch in lengths by considering the vectors $m(X)=2\sqrt{p(X)}$, instead. Any pmf that is identified as $m(X)$ is now two units away from the origin. The model space $\mathcal{M}$ is the collection of all transformed pmfs and represented as the surface of (the positive part of) the sphere in the right panel of Fig. 8. By representing the set of all possible pmfs of $X$ as $m(X)=2\sqrt{p(X)}$, we adopted our intuitive notion of distance. As a result, the selection mechanism underlying MDL can be made intuitive by simply looking at the forthcoming plots.

4.2 The individual-word and the only-mixed strategy

To ease the exposition, we assume that both words presented to the participant come from the right list $\mathcal{R}$, thus, $``rr"$ for the two models introduced below.

As model $\mathcal{M}_{1}$ we take the so-called individual-word strategy. Within this model $\mathcal{M}_{1}$, the parameter is $\theta_{1}=\vartheta$, which we interpret as the participant’s “right-list recognition ability”. With chance $\vartheta$ the participant then correctly recognizes that the first word originates from the right list and repeats this procedure for the second word, after which the participant categorizes the word pair as $L,M$, or $R$, see the left panel of Fig. 9 for a schematic description of this strategy as a processing tree. Fixing the participant’s “right-list recognition ability” $\vartheta$ yields the following pmf

[TABLE]

For instance, when the participant’s true ability is $\vartheta^{*}=0.9$, the three outcomes $[L,M,R]$ are then generated with the following three chances $f_{1}(X\,|\,0.9)=[0.01,0.18,0.81]$, which is plotted as a circle in Fig. 8. On the other hand, when $\vartheta^{*}=0.5$ the participant’s generating pmf is then $f_{1}(X\,|\,\vartheta=0.5)=[0.25,0.5,0.25]$, which is depicted as the cross in model space $\mathcal{M}$. The set of pmfs so defined forms a curve that goes through both the cross and the circle, see the left panel of Fig. 10.

As a competing model $\mathcal{M}_{2}$, we take the so-called only-mixed strategy. For the task described in Example 4.1, we might pose that participants from a certain clinical group are only capable of recognizing mixed word pairs and that they are unable to distinguish the pairs $``rr"$ from $``ll"$ resulting in a random guess between the responses $L$ and $R$, see the right panel of Fig. 9 for the processing tree. Within this model $\mathcal{M}_{2}$ the parameter is $\theta_{2}=\alpha$, which is interpreted as the participant’s “mixed-list differentiability skill” and fixing it yields the following pmf

[TABLE]

For instance, when the participant’s true differentiability is $\alpha^{*}=1/3$, the three outcomes $[L,M,R]$ are then generated with the equal chances $f_{2}(X\,|\,1/3)=[1/3,1/3,1/3]$, which, as before, is plotted as the dot in Fig. 10. On the other hand, when $\alpha^{*}=0.5$ the participant’s generating pmf is then given by $f_{2}(X\,|\,\alpha=0.5)=[0.25,0.5,0.25]$, i.e., the cross. The set of pmfs so defined forms a curve that goes through both the dot and the cross, see the left panel of Fig. 10.

The plots show that the models $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are neither saturated nor nested, as the two models define proper subsets of $\mathcal{M}$ and only overlap at the cross. Furthermore, the plots also show that $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are both one-dimensional, as each model is represented as a line in model space. Hence, the dimensionality terms in all three information criteria are the same. Moreover, AIC and BIC will only discriminate these two models based on goodness-of-fit alone. This particular model comparison, thus, allows us to highlight the role Fisher information plays in the MDL model selection philosophy.

4.3 Model characteristics

4.3.1 The maximum likelihood estimators

For FIA we need to compute the goodness-of-fit terms, thus, we need to identify the MLEs for the parameters within each model. For the models at hand, the MLEs are

[TABLE]

where $Y_{L},Y_{M}$ and $Y_{R}=n-Y_{L}-Y_{M}$ are the number of $L,M$ and $R$ responses in the data consisting of $n$ trials.

Estimation is a within model operation and it can be viewed as projecting the so-called empirical (i.e., observed) pmf corresponding to the data onto the model. For iid data with $w=3$ outcomes the empirical pmf corresponding to $x^{n}_{\textnormal{obs}}$ is defined as $\hat{p}_{\textnormal{obs}}(X)=[y_{L}/n,y_{M}/n,y_{R}/n]$. Hence, the empirical pmf gives the relative occurrence of each outcome in the sample. For instance, the observations $x^{n}_{\textnormal{obs}}$ consisting of $[y_{L}=3,y_{M}=3,y_{R}=3]$ responses correspond to the observed pmf $\hat{p}_{\textnormal{obs}}(X)=[1/3,1/3,1/3]$, i.e., the dot in Fig. 10. Note that this observed pmf $\hat{p}_{\textnormal{obs}}(X)$ does not reside on the curve of $\mathcal{M}_{1}$.

Nonetheless, when we use the MLE $\hat{\vartheta}$ of $\mathcal{M}_{1}$, we as researchers bestow the participant with a “right-list recognition ability” $\vartheta$ and implicitly assume that she used the individual-word strategy to generate the observations. In other words, we only consider the pmfs on the curve of $\mathcal{M}_{1}$ as viable explanations of how the participant generated her responses. For the data at hand, we have the estimate $\hat{\vartheta}_{\textnormal{obs}}=0.5$. If we were to generalize the observations $x^{n}_{\textnormal{obs}}$ under $\mathcal{M}_{1}$, we would then plug this estimate into the functional relationship $f_{1}$ resulting in the predictive pmf $f_{1}(X\,|\,\hat{\vartheta}_{\textnormal{obs}})=[0.25,0.5,0.25]$. Hence, even though the number of $L,M$ and $R$ responses were equal in the observations $x^{n}_{\textnormal{obs}}$, under $\mathcal{M}_{1}$ we expect that this participant will answer with twice as many $M$ responses compared to the $L$ and $R$ responses in a next set of test items. Thus, for predictions, part of the data is ignored and considered as noise.

Geometrically, the generalization $f_{1}(X\,|\,\hat{\vartheta}_{\textnormal{obs}})$ is a result of projecting the observed pmf $\hat{p}_{\textnormal{obs}}(X)$, i.e., the dot, onto the cross that does reside on the curve of $\mathcal{M}_{1}$.171717This resulting pmf $f_{1}(X\,|\,\hat{\vartheta}_{\textnormal{obs}})$ is also known as the Kullback-Leibler projection of the empirical pmf $\hat{p}_{\textnormal{obs}}(X)$ onto the model $\mathcal{M}_{1}$. White (1982) used this projection to study the behavior of the MLE under model misspecification. Observe that amongst all pmfs on $\mathcal{M}_{1}$, the projected pmf is closest to the empirical pmf $\hat{p}_{\textnormal{obs}}(X)$. Under $\mathcal{M}_{1}$ the projected pmf $f_{1}(X\,|\,\hat{\vartheta}_{\textnormal{obs}})$, i.e., the cross, is perceived as structural, while any deviations from the curve of $\mathcal{M}_{1}$ is labeled as noise. When generalizing the observations, we ignore noise. Hence, by estimating the parameter $\vartheta$, we implicitly restrict our predictions to only those pmfs that are defined by $\mathcal{M}_{1}$. Moreover, evaluating the prediction at $x^{n}_{\textnormal{obs}}$ and, subsequently, taking the negative logarithm yields the goodness-of-fit term; in this case, $-\log f_{1}(x^{n}_{\textnormal{obs}}\,|\,\hat{\vartheta}_{\textnormal{obs}}=0.5)=10.4$.

Which part of the data is perceived as structural or as noise depends on the model. For instance, when we use the MLE $\hat{\alpha}$, we restrict our predictions to the pmfs of $\mathcal{M}_{2}$. For the data at hand, we get $\hat{\alpha}_{\textnormal{obs}}=1/3$ and the plugin yields $f_{2}(X\,|\,\hat{\alpha}_{\textnormal{obs}})=[1/3,1/3,1/3]$. Again, amongst all pmfs on $\mathcal{M}_{2}$, the projected pmf is closest to the empirical pmf $\hat{p}_{\textnormal{obs}}(X)$. In this case, the generalization under $\mathcal{M}_{2}$ coincides with the observed pmf $\hat{p}_{\textnormal{obs}}(X)$. Hence, under $\mathcal{M}_{2}$ there is no noise, as the empirical pmf $\hat{p}_{\textnormal{obs}}(X)$ was already on the model. Geometrically, this means that $\mathcal{M}_{2}$ is closer to the empirical pmf than $\mathcal{M}_{1}$, which results in a lower goodness-of-fit term $-\log f_{2}(x^{n}_{\textnormal{obs}}\,|\,\hat{\alpha}_{\textnormal{obs}}=1/3)=9.9$.

This geometric interpretation allows us to make intuitive that data sets with the same goodness-of-fit terms will be as far from $\mathcal{M}_{1}$ as from $\mathcal{M}_{2}$. Equivalently, $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ identify the same amount of noise within $x^{n}_{\textnormal{obs}}$, when the two models fit the observations equally well. For instance, Fig. 10 shows that observations $x^{n}_{\textnormal{obs}}$ with an empirical pmf $\hat{p}_{\textnormal{obs}}(X)=[0.25,0.5,0.25]$ are equally far from $\mathcal{M}_{1}$ as from $\mathcal{M}_{2}$. Note that the closest pmf on $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ are both equal to the empirical pmf, as $f_{1}(X\,|\,\hat{\vartheta}_{\textnormal{obs}}=0.5)=\hat{p}_{\textnormal{obs}}(X)=f_{2}(X\,|\,\hat{\alpha}_{\textnormal{obs}}=1/2)$. As a result, the two goodness-of-fit terms will be equal to each other.

In sum, goodness-of-fit measures a model’s proximity to the observed data. Consequently, models that take up more volume in model space will be able to be closer to a larger number of data sets. In particular, when, say, $\mathcal{M}_{3}$ is nested within $\mathcal{M}_{4}$, this means that the distance between $\hat{p}_{\textnormal{obs}}(X)$ and $\mathcal{M}_{3}$ (noise) is at least the distance between $\hat{p}_{\textnormal{obs}}(X)$ and $\mathcal{M}_{4}$. Equivalently, for any data set, $\mathcal{M}_{4}$ will automatically label more of the observations as structural. Models that excessively identify parts of the observations as structural are known to overfit the data. Overfitting has an adverse effect on generalizability, especially when $n$ is small, as $\hat{p}_{\textnormal{obs}}(X)$ is then dominated by sampling error. In effect, the more voluminous model will then use this sampling error, rather than the structure, for its predictions. To guard ourselves from overfitting, thus, bad generalizability, the information criteria AIC, BIC and FIA all penalize for model complexity. AIC and BIC only do this via the dimensionality terms, while FIA also take the models’ volumes into account.

4.3.2 Geometrical complexity

For both models the dimensionality term is given by $\tfrac{1}{2}\log(\tfrac{n}{2\pi})$. Recall that the geometrical complexity term is the logarithm of the model’s volume, which for the individual-word and the only-mixed strategy are given by

[TABLE]

respectively. Hence, the individual-word strategy is a more complex model, because it has a larger volume, thus, capacity to fit data compared to the only-mixed strategy. After taking logs, we see that the individual-word strategy incurs an additional penalty of $1/2\log(2)$ compared to the only-mixed strategy.

4.4 Model selection based on the minimum description length principle

With all model characteristics at hand, we only need observations to illustrate that MDL model selection boils down to selecting the model that is closest to the observations with an additional penalty for model complexity.

Table 1 shows three data sets $x^{n}_{\textnormal{obs},1},x^{n}_{\textnormal{obs},2},x^{n}_{\textnormal{obs},3}$ with $n=30$ observations. The three associated empirical pmfs are plotted as the top, middle and lower rectangles in the right panel of Fig. 10, respectively. Table 1 also shows the approximation of each model’s description length using FIA. Note that the first observed pmf, the top rectangle in Fig. 10, is closer to $\mathcal{M}_{2}$ than to $\mathcal{M}_{1}$, while the third empirical pmf, the lower rectangle, is closer to $\mathcal{M}_{1}$. Of particular interest is the middle rectangle, which lies on an additional black curve that we refer to as a non-decision curve; observations that correspond to an empirical pmf that lies on this curve are described equally well by $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$. For this specific comparison, we have the following decision rule: FIA selects $\mathcal{M}_{2}$ as the preferred model whenever the observations correspond to an empirical pmf between the two non-decision curves, otherwise, FIA selects $\mathcal{M}_{1}$. Fig. 10 shows that FIA, indeed, selects the model that is closest to the data except in the area where the two models overlap –observations consisting of $n=30$ trials with an empirical pmf near the cross are considered better described by the simpler model $\mathcal{M}_{2}$. Hence, this yields an incorrect decision even when the empirical pmf is exactly equal to the true data generating pmf that is given by, say, $f_{1}(X\,|\,\vartheta=0.51)$. This automatic preference for the simpler model, however, decreases as $n$ increases.

The left and right panel of Fig. 11 show the non-decision curves when $n=120$ and $n$ (extremely) large, respectively. As a result of moving non-decision bounds, the data set $x^{n}_{\textnormal{obs},4}=[56,40,24]$ that has the same observed pmf as $x^{n}_{\textnormal{obs},2}$, i.e., the middle rectangle, will now be better described by model $\mathcal{M}_{1}$.

For (extremely) large $n$, the additional penalty due to $\mathcal{M}_{1}$ being more voluptuous than $\mathcal{M}_{2}$ becomes irrelevant and the sphere is then separated into quadrants: observations corresponding to an empirical pmf in the top-left or bottom-right quadrant are better suited to the only-mixed strategy, while the top-right and bottom-left quadrants indicate a preference for the individual-word strategy $\mathcal{M}_{1}$. Note that pmfs on the non-decision curves in the right panel of Fig. 11 are as far apart from $\mathcal{M}_{1}$ as from $\mathcal{M}_{2}$, which agrees with our geometric interpretation of goodness-of-fit as a measure of the model’s proximity to the data. This quadrant division is only based on the two models’ goodness-of-fit terms and yields the same selection as one would get from BIC (e.g., Rissanen, 1996). For large $n$, FIA, thus, selects the model that is closest to the empirical pmf. This behavior is desirable, because asymptotically the empirical pmf is not distinguishable from the true data generating pmf. As such, the model that is closest to the empirical pmf will then also be closest to the true pmf. Hence, FIA asymptotically selects the model that is closest to the true pmf. As a result, the projected pmf within the closest model is then expected to yield the best predictions amongst the competing models.

4.5 Fisher information and generalizability

Model selection by MDL is sometimes perceived as a formalization of Occam’s razor (e.g., Balasubramanian, 1996; Grünwald, 1998), a principle that states that the most parsimonious model should be chosen when the models under consideration fit the observed data equally well. This preference for the parsimonious model is based on the belief that the simpler model is better at predicting new (as yet unseen) data coming from the same source, as was shown by Pitt, Myung and Zhang (2002) with simulated data.

To make intuitive why the more parsimonious model, on average, leads to better predictions, we assume, for simplicity, that the true data generating pmf is given by $f(X\,|\,\theta^{*})$, thus, the existence of a true parameter value $\theta^{*}$. As the observations are expected to be contaminated with sampling error, we also expect an estimation error, i.e., a distance $\textnormal{d}\theta$ between the maximum likelihood estimate $\hat{\theta}_{\textnormal{obs}}$ and the true $\theta^{*}$. Recall that in the construction of Jeffreys’s prior Fisher information was used to convert displacement in model space to distances on parameter space. Conversely, Fisher information transforms the estimation error in parameter space to a generalization error in model space. Moreover, the larger the Fisher information at $\theta^{*}$ is, the more it will expand the estimation error into a displacement between the prediction $f(X\,|\,\hat{\theta}_{\textnormal{obs}})$ and the true pmf $f(X\,|\,\theta^{*})$. Thus, a larger Fisher information at $\theta^{*}$ will push the prediction further from the true pmf resulting in a bad generalization. Smaller models have, on average, a smaller Fisher information at $\theta^{*}$ and will therefore lead to more stable predictions that are closer to the true data generating pmf. Note that the generalization scheme based on the MLE plugin $f(X\,|\,\hat{\theta}_{\textnormal{obs}})$ ignores the error at each generalization step. The Bayesian counterpart, on the other hand, does take these errors into account, see Dawid (2011), Ly et al. (2017b), Marsman, Ly and Wagenmakers (2016) and see van Erven, Grünwald and De Rooij (2012), Grünwald and Mehta (2016), van der Pas and Grünwald (2014), Wagenmakers, Grünwald and Steyvers (2006) for a prequential view of generalizability.

5 Concluding Comments

Fisher information is a central statistical concept that is of considerable relevance for mathematical psychologists. We illustrated the use of Fisher information in three different statistical paradigms: in the frequentist paradigm, Fisher information was used to construct hypothesis tests and confidence intervals; in the Bayesian paradigm, Fisher information was used to specify a default, parameterization-invariant prior distribution; lastly, in the paradigm of information theory, data compression, and minimum description length, Fisher information was used to measure model complexity. Note that these three paradigms highlight three uses of the functional relationship $f$ between potential observations $x^{n}$ and the parameters $\theta$. Firstly, in the frequentist setting, the second argument was fixed at a supposedly known parameter value $\theta_{0}$ or $\hat{\theta}_{\textnormal{obs}}$ resulting in a probability mass function, a function of the potential outcomes $f(\cdot\,\,|\,\theta_{0})$. Secondly, in the Bayesian setting, the first argument was fixed at the observed data resulting in a likelihood function, a function of the parameters $f(x_{\textnormal{obs}}\,|\,\cdot)$. Lastly, in the information geometric setting both arguments were free to vary, i.e., $f(\cdot\,|\,\cdot)$ and plugged in by the observed data and the maximum likelihood estimate.

To ease the exposition we only considered Fisher information of one-dimensional parameters. The generalization of the concepts introduced here to vector valued $\theta$ can be found in the appendix. A complete treatment of all the uses of Fisher information throughout statistics would require a book (e.g., Frieden, 2004) rather than a tutorial article. Due to the vastness of the subject, the present account is by no means comprehensive. Our goal was to use concrete examples to provide more insight about Fisher information, something that may benefit psychologists who propose, develop, and compare mathematical models for psychological processes. Other uses of Fisher information are in the detection of model misspecification (Golden, 1995; Golden, 2000; Waldorp, Huizenga and Grasman, 2005; Waldorp, 2009; Waldorp, Christoffels and van de Ven, 2011; White, 1982), in the reconciliation of frequentist and Bayesian estimation methods through the Bernstein-von Mises theorem (Bickel and Kleijn, 2012; Rivoirard and Rousseau, 2012; van der Vaart, 1998; Yang and Le Cam, 2000), in statistical decision theory (e.g., Berger, 1985; Hájek, 1972; Korostelev and Korosteleva, 2011; Ray and Schmidt-Hieber, 2016; Wald, 1949), in the specification of objective priors for more complex models (e.g., Ghosal, Ghosh and Ramamoorthi, 1997; Grazian and Robert, 2015; Kleijn and Zhao, 2017), and computational statistics and generalized MCMC sampling in particular (e.g., Banterle et al., 2015; Girolami and Calderhead, 2011; Grazian and Liseo, 2014; Gronau et al., 2017).

In sum, Fisher information is a key concept in statistical modeling. We hope to have provided an accessible and concrete tutorial article that explains the concept and some of its uses for applications that are of particular interest to mathematical psychologists.

Appendix A Generalization to Vector-Valued Parameters: The Fisher Information Matrix

Let $X$ be a random variable, $\vec{\theta}=(\theta_{1},\ldots,\theta_{d})$ a vector of parameters, and $f$ a functional relationship that relates $\vec{\theta}$ to the potential outcomes $x$ of $X$. As before, it is assumed that by fixing $\vec{\theta}$ in $f$ we get the pmf $p_{\vec{\theta}}(x)=f(x\,|\,\vec{\theta})$, which is a function of $x$. The pmf $p_{\vec{\theta}}(x)$ fully determines the chances with which $X$ takes on the events in the outcome space $\mathcal{X}$. The Fisher information of the vector $\vec{\theta}\in\mathbb{R}^{d}$ is a positive semidefinite symmetric matrix of dimension $d\times d$ with the entry at the $i$th row and $j$th column given by

[TABLE]

where $l(x\,|\,\vec{\theta})=\log f(x\,|\,\vec{\theta})$ is the log-likelihood function, $\tfrac{\partial}{\partial\theta_{i}}l(x\,|\,\vec{\theta})$ is the score function, that is, the partial derivative with respect to the $i$th component of the vector $\vec{\theta}$ and the dot is short-hand notation for the vector of the partial derivatives with respect to $\theta=(\theta_{1},\ldots,\theta_{d})$. Thus, $\dot{l}(x\,|\,\vec{\theta})$ is a $d\times 1$ column vector of score functions, while $\dot{l}^{T}(x\,|\,\vec{\theta})$ is a $1\times d$ row vector of score functions at the outcome $x$. The partial derivative is evaluated at $\vec{\theta}$, the same $\vec{\theta}$ that is used in the pmf $p_{\vec{\theta}}(x)$ for the weighting. In Appendix E it is shown that the score functions are expected to be zero, which explains why $I_{X}(\vec{\theta})$ is a covariance matrix.

Under mild regularity conditions the $i,j$th entry of the Fisher information matrix can be equivalently calculated via the negative expectation of the second order partial derivates, that is,

[TABLE]

Note that the sum (thus, integral in the continuous case) is with respect to the outcomes $x$ of $X$.

Example A.1** (Fisher information for normally distributed random variables).**

When $X$ is normally distributed, i.e., $X\sim\mathcal{N}(\mu,\sigma^{2})$, it has the following probability density function (pdf)

[TABLE]

where the parameters are collected into the vector $\vec{\theta}={\mu\choose\sigma}$, with $\mu\in\mathbb{R}$ and $\sigma>0$. The score vector at a specific $\theta={\mu\choose\sigma}$ is the following vector of functions of $x$

[TABLE]

The unit Fisher information matrix $I_{X}(\vec{\theta})$ is a $2\times 2$ symmetric positive semidefinite matrix, consisting of expectations of partial derivatives. Equivalently, $I_{X}(\vec{\theta})$ can be calculated using the second order partials derivatives

[TABLE]

The off-diagonal elements are in general not zero. If the $i,j$th entry is zero we say that $\theta_{i}$ and $\theta_{j}$ are orthogonal to each other, see Appendix C.3.3 below. $\hfill\diamond$

For iid trials $X^{n}=(X_{1},\ldots,X_{n})$ with $X\sim p_{\theta}(x)$, the Fisher information matrix for $X^{n}$ is given by $I_{X^{n}}(\vec{\theta})=nI_{X}(\vec{\theta})$. Thus, for vector-valued parameters $\vec{\theta}$ the Fisher information matrix remains additive.

In the remainder of the text, we simply use $\theta$ for both one-dimensional and vector-valued parameters. Similarly, depending on the context it should be clear whether $I_{X}(\theta)$ is a number or a matrix.

Appendix B Frequentist Statistics based on Asymptotic Normality

The construction of the hypothesis tests and confidence intervals in the frequentist section were all based on the MLE being asymptotically normal.

B.1 Asymptotic normality of the MLE for vector-valued parameters

For so-called regular parametric models, see Appendix E, the MLE for vector-valued parameters $\theta$ converges in distribution to a multivariate normal distribution, that is,

[TABLE]

where $\mathcal{N}_{d}$ is a $d$-dimensional multivariate normal distribution, and $I_{X}^{-1}(\theta^{*})$ the inverse Fisher information matrix at the true value $\theta^{*}$. For $n$ large enough, we can, thus, approximate the sampling distribution of the “error” of the MLE by a normal distribution, thus,

[TABLE]

In practice, we fix $n$ and replace the true sampling distribution by this normal distribution. Hence, we incur an approximation error that is only negligible whenever $n$ is large enough. What constitutes $n$ large enough depends on the true data generating pmf $p^{*}(x)$ that is unknown in practice. In other words, the hypothesis tests and confidence intervals given in the main text based on the replacement of the true sampling distribution by this normal distribution might not be appropriate. In particular, this means that a hypothesis tests at a significance level of 5% based on the asymptotic normal distribution, instead of the true sampling distribution, might actually yield a type 1 error rate of, say, 42%. Similarly, as a result of the approximation error, a 95%-confidence interval might only encapsulate the true parameter in, say, 20% of the time that we repeat the experiment.

B.2 Asymptotic normality of the MLE and the central limit theorem

Asymptotic normality of the MLE can be thought of as a refinement of the central limit theorem. The (Lindeberg-Lévy) CLT is a general statement about the sampling distribution of the sample mean estimator $\bar{X}=\tfrac{1}{n}\sum_{i=1}^{n}X_{i}$ based on iid trials of $X$ with common population mean $\theta=E(X)$ and variance $\text{Var}(X)<\infty$. More specifically, the CLT states that, with a proper scaling, the sample mean $\bar{X}$ centred around the true $\theta^{*}$ will converge in distribution to a normal distribution, that is, \sqrt{n}(\bar{X}-\theta^{*})\overset{D}{\rightarrow}\mathcal{N}\big{(}0,\text{Var}(X)\big{)}. In practice, we replace the true sampling distribution by this normal distribution at fixed $n$ and hope that $n$ is large enough. Hence, for fixed $n$ we then suppose that the “error” is distributed as $(\bar{X}-\theta^{*})\overset{D}{\approx}\mathcal{N}(0,\tfrac{1}{n}\text{Var}(X))$ and we ignore the approximation error. In particular, when we know that the population variance is $\text{Var}(X)=1$, we then know that we require an experiment with $n=100$ samples for $\bar{X}$ to generate estimates within $0.196$ distance from $\theta$ with approximately 95% chance, that is, $P(|\bar{X}-\theta|\leq 0.196)\approx 0.95$.181818As before, chance refers to the relative frequency, that is, when we repeat the experiment $k=200$ times, each with $n=100$, we get $k$ number of estimates and approximately 95% of these $k$ number of estimates are then expected to be within $0.196$ distance away from the true population mean $\theta^{*}$. This calculation was based on our knowledge of the normal distribution $\mathcal{N}(0,0.01)$, which has its 97.5% quantile at 0.196. In the examples below we re-use this calculation by matching the asymptotic variances to $0.01$.191919Technically, an asymptotic variance is free of $n$, but we mean the approximate variance at finite $n$. For the CLT this means $\tfrac{1}{n}\sigma^{2}$. The 95% statement only holds approximately, because we do not know whether $n=100$ is large enough for the CLT to hold, i.e., this probability could be well below 23%. Note that the CLT holds under very general conditions; the population mean and variance both need to exist, i.e., be finite. The distributional form of $X$ is irrelevant for the statement of the CLT.

On the other hand, to even compute the MLE we not only require that the population quantities to exists and be finite, but we also need to know the functional relationship $f$ that relates these parameters to the outcomes of $X$. When we assume more (and nature adheres to these additional conditions), we know more, and are then able to give stronger statements. We give three examples.

Example B.1** (Asymptotic normality of the MLE vs the CLT: The Gaussian distribution).**

If $X$ has a Gaussian (normal) distribution, i.e., $X\sim\mathcal{N}(\theta,\sigma^{2})$, with $\sigma^{2}$ known, then the MLE is the sample mean and the unit Fisher information is $I_{X}(\theta)=1/\sigma^{2}$. Asymptotic normality of the MLE leads to the same statement as the CLT, that is, $\sqrt{n}(\hat{\theta}-\theta^{*})\overset{D}{\rightarrow}\mathcal{N}(0,\sigma^{2})$. Hence, asymptotically we do not gain anything by going from the CLT to asymptotic normality of the MLE. The additional knowledge of $f(x\,|\,\theta)$ being normal does, however, allow us to come to the rare conclusion that the normal approximation holds exactly for every finite $n$, thus, $(\hat{\theta}-\theta^{*})\overset{D}{=}\mathcal{N}(0,\tfrac{1}{n}\sigma^{2})$. In all other cases, whenever $X\not\sim\mathcal{N}(\theta,\sigma^{2})$, we always have an approximation.202020This is a direct result of Cramér’s theorem that states that whenever $X$ is independent of $Y$ and $Z=X+Y$ with $Z$ a normal distribution, then $X$ and $Y$ themselves are necessarily normally distributed. Thus, whenever $\sigma^{2}=1$ and $n=100$ we know that $P(|\hat{\theta}-\theta^{*}|\leq 0.196)=0.95$ holds exactly. $\hfill\diamond$

Example B.2** (Asymptotic normality of the MLE vs the CLT: The Laplace distribution).**

If $X$ has a Laplace distribution with scale $b$, i.e., $X\sim\text{Laplace}(\theta,b)$, then its population mean and variance are $\theta=E(X)$ and $2b^{2}=\text{Var}(X)$, respectively.

In this case, the MLE is the sample median $\hat{M}$ and the unit Fisher information is $I_{X}(\theta)=1/b^{2}$. Asymptotic normality of the MLE implies that we can approximate the sampling distribution by the normal distribution, that is, $(\hat{\theta}-\theta^{*})\overset{D}{\approx}\mathcal{N}(0,\tfrac{1}{n}b^{2})$, when $n$ is large enough. Given that the population variance is $\text{Var}(X)=1$, we know that $b=1/\sqrt{2}$, yielding a variance of $\tfrac{1}{2n}$ in our normal approximation to the sampling distribution. Matching this variance to $0.01$ shows that we now require only $n=50$ samples for the estimator to generate estimates within 0.196 distance away from the true value $\theta^{*}$ with 95% chance. As before, the validity of this statement only holds approximately, i.e., whenever the normal approximation to the sampling distribution of the MLE at $n=50$ is not too bad.

Hence, the additional knowledge of $f(x\,|\,\theta)$ being Laplace allows us to use an estimator, i.e., the MLE, that has a lower asymptotic variance. Exploiting this knowledge allowed us to design an experiment with twice as few participants. $\hfill\diamond$

Example B.3** (Asymptotic normality of the MLE vs the CLT: The Cauchy distribution).**

If $X$ has a Cauchy distribution centred around $\theta$ with scale $1$, i.e., $X\sim\text{Cauchy}(\theta,1)$, then $X$ does not have a finite population variance, nor a finite population mean. As such, the CLT cannot be used. Even worse, Fisher (1922) showed that the sample mean as an estimator for $\theta$ is in this case useless, as the sampling distribution of the sample mean is a Cauchy distribution that does not depend on $n$, namely, $\bar{X}\sim\text{Cauchy}(\theta,1)$. As such, using the first observation alone to estimate $\theta$ is as good as combining the information of $n=100$ samples in the sample mean estimator. Hence, after seeing the first observation no additional information about $\theta$ is gained using the sample mean $\bar{X}$, not even if we increase $n$.

The sample median estimator $\hat{M}$ performs better. Again, Fisher (1922) already knew that for $n$ large enough that $(\hat{M}-\theta^{*})\overset{D}{\approx}\mathcal{N}(0,\tfrac{1}{n}\tfrac{\pi^{2}}{2})$. The MLE is even better, but unfortunately, in this case, it cannot be given as an explicit function of the data.212121Given observations $x^{n}_{\textnormal{obs}}$ the maximum likelihood estimate $\hat{\theta}_{\textnormal{obs}}$ is the number for which the score function $\dot{l}(x^{n}_{\textnormal{obs}}\,|\,\theta)=\sum_{i=1}^{n}\tfrac{2(x_{\textnormal{obs},i}-\theta)}{1+(x_{\textnormal{obs},i}-\theta)^{2}}$ is zero. This optimization cannot be solved analytically and there are $2n$ solutions to this equation. The Fisher information can be given explicitly, namely, $I_{X}(\theta)=1/2$. Asymptotic normality of the MLE implies that $(\hat{\theta}-\theta^{*})\overset{D}{\approx}\mathcal{N}(0,\tfrac{1}{n}2)$, when $n$ is large enough. Matching the variances in the approximation based on the normal distribution to $0.01$ shows that we require $n=25\pi^{2}\approx 247$ for the sample median and $n=200$ samples for the MLE to generate estimates within 0.196 distance away from the true value of value $\theta^{*}$ with approximate 95% chance. $\hfill\diamond$

B.3 Efficiency of the MLE: The Hájek-LeCam convolution theorem and the Cramér-Fréchet-Rao information lower bound

The previous examples showed that the MLE is an estimator that leads to a smaller sample size requirement, because it is the estimator with the lower asymptotic variance. This lower asymptotic variance is a result of the MLE making explicit use of the functional relationship between the samples $x^{n}_{\textnormal{obs}}$ and the target $\theta$ in the population. Given any such $f$, one might wonder whether the MLE is the estimator with the lowest possible asymptotic variance. The answer is affirmative, whenever we restrict ourselves to the broad class of so-called regular estimators.

A regular estimator $T_{n}=t_{n}(X_{n})$ is a function of the data that has a limiting distribution that does not change too much, whenever we change the parameters in the neighborhood of the true value $\theta^{*}$, see van der Vaart (1998, p. 115) for a precise definition. The Hájek-LeCam convolution theorem characterizes the aforementioned limiting distribution as a convolution, i.e., a sum of, the independent statistics $\Delta_{\theta^{*}}$ and $Z_{\theta^{*}}$. That is, for any regular estimator $T_{n}$ and every possible true value $\theta^{*}$ we have

[TABLE]

where Z_{\theta^{*}}\sim\mathcal{N}\big{(}0,I^{-1}_{X}(\theta^{*})\big{)} and where $\Delta_{\theta^{*}}$ has an arbitrary distribution. By independence, the variance of the asymptotic distribution is simply the sum of the variances. As the variance of $\Delta_{\theta^{*}}$ cannot be negative, we know that the asymptotic variance of any regular estimator $T_{n}$ is bounded from below, that is, $\text{Var}(\Delta_{\theta^{*}})+I_{X}^{-1}(\theta^{*})\geq I^{-1}_{X}(\theta^{*})$.

The MLE is a regular estimator with $\Delta_{\theta^{*}}$ equal to the fixed true value $\theta^{*}$, thus, $\text{Var}(\Delta_{\theta^{*}})=0$. As such, the MLE has an asymptotic variance $I^{-1}_{X}(\theta^{*})$ that is equal to the lower bound given above. Hence, amongst the broad class of regular estimators, the MLE performs best. This result was already foreshadowed by Fisher (1922), though it took another 50 years before this statement was made mathematically rigorous (Hájek, 1970; Inagaki, 1970; LeCam, 1970; van der Vaart, 2002; Yang, 1999), see also Ghosh (1985) for a beautiful review.

We stress that the normal approximation to the true sampling distribution only holds when $n$ is large enough. In practice, $n$ is relatively small and the replacement of the true sampling distribution by the normal approximation can, thus, lead to confidence intervals and hypothesis tests that perform poorly (Brown, Cai and DasGupta, 2001). This can be very detrimental, especially, when we are dealing with hard decisions such as the rejection or non-rejection of a hypothesis.

A simpler version of the Hájek-LeCam convolution theorem is known as the Cramér-Fréchet-Rao information lower bound (Cramér, 1946; Fréchet, 1943; Rao, 1945), which also holds for finite $n$. This theorem states that the variance of an unbiased estimator $T_{n}$ cannot be lower than the inverse Fisher information, that is, $n\text{Var}(T_{n})\geq I_{X}^{-1}(\theta^{*})$. We call an estimator $T_{n}=t(X^{n})$ unbiased if for every possible true value $\theta^{*}$ and at each fixed $n$, its expectation is equal to the true value, that is, $E(T_{n})=\theta^{*}$. Hence, this lower bound shows that Fisher information is not only a concept that is useful for large samples.

Unfortunately, the class of unbiased estimators is rather restrictive (in general, it does not include the MLE) and the lower bound cannot be attained whenever the parameter is of more than one dimensions (Wijsman, 1973). Consequently, for vector-valued parameters $\theta$, this information lower bound does not inform us, whether we should stop our search for a better estimator.

The Hájek-LeCam convolution theorem implies that for $n$ large enough the MLE $\hat{\theta}$ is the best performing statistic. For the MLE to be superior, however, the data do need to be generated as specified by the functional relationship $f$. In reality, we do not know whether the data are indeed generated as specified by $f$, which is why we should also try to empirically test such an assumption. For instance, we might believe that the data are normally distributed, while in fact they were generated according to a Cauchy distribution. This incorrect assumption implies that we should use the sample mean, but Example B.3 showed the futility of such estimator. Model misspecification, in addition to hard decisions based on the normal approximation, might be the main culprit of the crisis of replicability. Hence, more research on the detection of model misspecification is desirable and expected (e.g., Grünwald, 2016; Grünwald and van Ommen, 2014; van Ommen et al., 2016).

Appendix C Bayesian use of the Fisher-Rao Metric: The Jeffreys’s Prior

We make intuitive that the Jeffreys’s prior is a uniform prior on the model $\mathcal{M}_{\Theta}$, i.e.,

[TABLE]

where J_{m}=\big{(}m_{\theta_{a}}(X),m_{\theta_{b}}(X)\big{)} is an interval of pmfs in model space. To do so, we explain why the differential $\textnormal{d}m_{\theta}(X)$, a displacement in model space, is converted into $\sqrt{I_{X}(\theta)}\textnormal{d}\theta$ in parameter space. The elaboration below boils down to an explanation of arc length computations using integration by substitution.

C.1 Tangent vectors

First note that we swapped the area of integration by substituting the interval J_{m}=\big{(}m_{\theta_{a}}(X),m_{\theta_{b}}(X)\big{)} consisting of pmfs in function space $\mathcal{M}_{\Theta}$ by the interval $(\theta_{a},\theta_{b})$ in parameter space. This is made possible by the parameter functional $\nu$ with domain $\mathcal{M}_{\Theta}$ and range $\Theta$ that uniquely assigns to any (transformed) pmf $m_{a}(X)\in\mathcal{M}_{\Theta}$ a parameter value $\theta_{a}\in\Theta$. In this case, we have $\theta_{a}=\nu(m_{a}(X))=(\tfrac{1}{2}m_{a}(1))^{2}$. Uniqueness of the assignment implies that the resulting parameter values $\theta_{a}$ and $\theta_{b}$ in $\Theta$ differ from each other whenever $m_{a}(X)$ and $m_{b}(X)$ in $\mathcal{M}_{\Theta}$ differ from each other.

For example, the map $\nu:\mathcal{M}_{\Theta}\rightarrow\Theta$ implies that in the left panel of Fig. 12 the third square from the left with coordinates $m_{a}(X)=[0.89,1.79]$ can be labeled by $\theta_{a}=0.8\approx(\tfrac{1}{2}(1.79))^{2}$, while the second square from the left with coordinates $m_{b}(X)=[0.63,1.90]$ can be labeled by $\theta_{b}=0.9\approx(\tfrac{1}{2}(1.90))^{2}$.

To calculate the arc length of the curve $J_{m}$ consisting of functions in $\mathcal{M}_{\Theta}$, we first approximate $J_{m}$ by a finite sum of tangent vectors, i.e., straight lines. The approximation of the arc length is the sum of the length of these straight lines. The associated approximation error goes to zero, when we increase the number of tangent vectors and change the sum into an integral sign, as in the usual definition of an integral. First we discuss tangent vectors.

In the left panel in Fig. 12, we depicted the tangent vector at $m_{\theta_{a}}(X)$ as the full arrow. This full arrow is constructed from its components: one broken arrow that is parallel to the horizontal axis associated with the outcome $x=0$, and one broken arrow that is parallel to the vertical axis associated with the outcome $x=1$. The arrows parallel to the axes are derived by first fixing $X=x$ followed by a Taylor expansion of the parameterization $\theta\mapsto m_{\theta}(x)$ at $\theta_{a}$. The Taylor expansion is derived by differentiating with respect to $\theta$ at $\theta_{a}$ yielding the following “linear” function of the distance $\textnormal{d}\theta=|\theta_{b}-\theta_{a}|$ in parameter space,

[TABLE]

where the slope, a function of $x$, $A_{\theta_{a}}(x)$ at $m_{\theta_{a}}(x)$ in the direction of $x$ is given by

[TABLE]

and with an “intercept” $B_{\theta_{a}}(x)=o(\textnormal{d}\theta)$ that goes fast to zero whenever $\textnormal{d}\theta\rightarrow 0$. Thus, for $\textnormal{d}\theta$ small, the intercept $B_{\theta_{a}}(x)$ is practically zero. Hence, we approximate the displacement between $m_{\theta_{a}}(x)$ and $m_{\theta_{b}}(x)$ by a straight line.

Example C.1** (Tangent vectors).**

In the right panel of Fig. 12 the right most triangle is given by $m_{\phi_{a}}(X)=[1.25,1.56]$, while the triangle in the middle refers to $m_{\phi_{b}}(X)=[0.99,1.74]$. Using the functional $\tilde{\nu}$, i.e., the inverse of the parameterization, $\phi\mapsto 2\sqrt{f(x\,|\,\phi)}$, where f(x\,|\,\phi)=\big{(}\frac{1}{2}+\frac{1}{2}\big{(}\tfrac{\phi}{\pi}\big{)}^{3}\big{)}^{x}\big{(}\frac{1}{2}-\frac{1}{2}\big{(}\tfrac{\phi}{\pi}\big{)}^{3}\big{)}^{1-x}, we find that these two pmfs correspond to $\phi_{a}=0.6\pi$ and $\phi_{b}=0.8\pi$.

The tangent vector at $m_{\phi_{a}}(X)$ is constructed from its components. For the horizontal displacement, we fill in $x=0$ in $\log f(x\,|\,\phi)$ followed by the derivation with respect to $\phi$ at $\phi_{a}$ and a multiplication by $m_{\phi_{a}}(x)$ resulting in

[TABLE]

where $\textnormal{d}\phi=|\phi_{b}-\phi_{a}|$ is the distance in parameter space $\Phi$. The minus sign indicates that the displacement along the horizontal axis is from right to left. Filling in $\textnormal{d}\phi=|\phi_{b}-\phi_{a}|=0.2\pi$ and $\phi_{a}=0.6\pi$ yields a horizontal displacement of $0.17$ at $m_{\phi_{a}}(0)$ from right to left in model space. Similarly, the vertical displacement in terms of $\phi$ is calculated by first filling in $x=1$ and leads to

[TABLE]

By filling in $\textnormal{d}\phi=0.2$ and $\phi_{a}=0.6\pi$, we see that a change of $\textnormal{d}\phi=0.2\pi$ at $\phi_{a}=0.6\pi$ in the parameter space corresponds to a vertical displacement of $0.14$ at $m_{\phi_{a}}(1)$ from bottom to top in model space. Note that the axes in Fig. 12 are scaled differently.

The combined displacement $\frac{\textnormal{d}m_{\phi_{a}}(X)}{\textnormal{d}\phi}\textnormal{d}\phi$ at $m_{\phi_{a}}(X)$ is the sum of the two broken arrows and plotted as a full arrow in the right panel of Fig. 12. $\hfill\diamond$

The length of the tangent vector $\frac{\textnormal{d}m_{\theta_{a}}(X)}{\textnormal{d}\theta}$ at the vector $m_{\theta_{a}}(X)$ is calculated by taking the root of the sum of its squared component, the natural measure of distance we adopted above and this yields

[TABLE]

The second equality follows from the definition of ${\textnormal{d}m_{\theta_{a}}(X)\over\textnormal{d}\theta}$, i.e., Eq. (C.3), and the last equality is due to the definition of Fisher information.

Example C.2** (name=Length of the tangent vectors, continue=tangentVectors).**

The length of the tangent vector in the right panel of Fig. 12 can be calculated as the root of the sums of squares of its components, that is, $\|\tfrac{\textnormal{d}m_{\phi_{a}}(X)}{\textnormal{d}\phi}\textnormal{d}\phi\|_{2}=\sqrt{(-0.14)^{2}+0.17^{2}}=0.22$. Alternatively, we can first calculate the square root of the Fisher information at $\phi_{a}=0.6\pi$, i.e.,

[TABLE]

and a multiplication with $\textnormal{d}\phi=0.2\pi$ results in $\|\tfrac{\textnormal{d}m_{\phi_{a}}(X)}{\textnormal{d}\phi}\|_{2}\textnormal{d}\phi=0.22$. $\hfill\diamond$

More generally, to approximate the length between pmfs $m_{\theta_{a}}(X)$ and $m_{\theta_{b}}(X)$, we first identify $\nu(m_{\theta_{a}}(X))=\theta_{a}$ and multiply this with the distance $\textnormal{d}\theta=|\theta_{a}-\nu(m_{\theta_{b}}(X))|$ in parameter space, i.e.,

[TABLE]

In other words, the root of the Fisher information converts a small distance $\textnormal{d}\theta$ at $\theta_{a}$ to a displacement in model space at $m_{\theta_{a}}(X)$.

C.2 The Fisher-Rao metric

By virtue of the parameter functional $\nu$, we send an interval of pmfs J_{m}=\Big{(}m_{\theta_{a}}(X),m_{\theta_{b}}(X)\Big{)} in the function space $\mathcal{M}_{\Theta}$ to the interval $(\theta_{a},\theta_{b})$ in the parameter space $\Theta$. In addition, with the conversion of $\textnormal{d}m_{\theta}(X)=\sqrt{I_{X}(\theta)}\,\textnormal{d}\theta$ we integrate by substitution, that is,

[TABLE]

In particular, choosing $J_{\theta}=\mathcal{M}_{\Theta}$ yields the normalizing constant $V=\int_{0}^{1}\sqrt{I_{X}(\theta)}\textnormal{d}\theta$. The interpretation of $V$ as being the total length of $\mathcal{M}_{\Theta}$ is due to the use of $\textnormal{d}m_{\theta}(X)$ as the metric, a measure of distance, in model space. To honour Calyampudi Radhakrishna Rao’s (1945) contribution to the theory, this metric is also known as the Fisher-Rao metric (e.g., Amari et al., 1987; Atkinson and Mitchell, 1981; Burbea, 1984; Burbea and Rao, 1982, 1984; Dawid, 1977; Efron, 1975; Kass and Vos, 2011).

C.3 Fisher-Rao metric for vector-valued parameters

C.3.1 The parameter functional $\nu:\mathcal{P}\rightarrow B$ and the categorical distribution

For random variables with $w$ number of outcomes, the largest set of pmfs $\mathcal{P}$ is the collection of functions $p$ on $\mathcal{X}$ such that (i) $0\leq p(x)=P(X=x)$ for every outcome $x$ in $\mathcal{X}$, and (ii) to explicitly convey that there are $w$ outcomes, and none more, these $w$ chances have to sum to one, that is, $\sum_{x\in\mathcal{X}}p(x)=1$. The complete set of pmfs $\mathcal{P}$ can be parameterized using the functional $\nu$ that assigns to each $w$-dimensional pmf $p(X)$ a parameter $\beta\in\mathbb{R}^{w-1}$.

For instance, given a pmf $p(X)=[p(L),p(M),p(R)]$ we typically use the functional $\nu:\mathcal{P}\rightarrow\mathbb{R}^{2}$ that takes the first two coordinates, that is, $\nu(p(X))=\beta={\beta_{1}\choose\beta_{2}}$, where $\beta_{1}=p(L)$ and $\beta_{2}=p(M)$. The range of this functional $\nu$ is the parameter space $B=[0,1]\times[0,\beta_{1}]$. Conversely, the inverse of the functional $\nu$ is the parameterization $\beta\mapsto p_{\beta}(X)=[\beta_{1},\beta_{2},1-\beta_{1}-\beta_{2}]$, where (i’) $0\leq\beta_{1},\beta_{2}$ and (ii’) $\beta_{1}+\beta_{2}\leq 1$. The restrictions (i’) and (ii’) imply that the parameterization has domain $B$ and the largest set of pmfs $\mathcal{P}$ as its range. By virtue of the functional $\nu$ and its inverse, that is, the parameterization $\beta\mapsto p_{\beta}(X)$, we conclude that the parameter space $B$ and the complete set of pmfs $\mathcal{P}$ are isomorphic. This means that each pmf $p(X)\in\mathcal{P}$ can be uniquely identified with a parameter $\beta\in B$ and vice versa. The inverse of $\nu$ implies that the parameters $\beta\in B$ are functionally related to the potential outcomes $x$ of $X$ as

[TABLE]

where $x_{L},x_{M}$ and $x_{R}$ are the number of $L,M$ and $R$ responses in one trial –we either have $x=[x_{L},x_{M},x_{R}]=[1,0,0]$, $x=[0,1,0]$, or $x=[0,0,1]$. The model $f(x\,|\,\beta)$ can be regarded as the generalization of the Bernoulli model to $w=3$ categories. In effect, the parameters $\beta_{1}$ and $\beta_{2}$ can be interpreted as a participant’s propensity of choosing $L$ and $M$, respectively. If $X^{n}$ consists of $n$ iid categorical random variables with the outcomes $[L,M,R]$, the joint pmf of $X^{n}$ is then

[TABLE]

where $y_{L},y_{M}$ and $y_{R}=n-y_{L}-y_{M}$ are the number of $L,M$ and $R$ responses in $n$ trials. As before, the representation of the pmfs as the vectors $m_{\beta}(X)=[2\sqrt{\beta_{1}},2\sqrt{\beta_{2}},2\sqrt{1-\beta_{1}-\beta_{2}}]$ form the surface of (the positive part of) the sphere of radius two, thus, $\mathcal{M}=\mathcal{M}_{B}$, see Fig. 13. The extreme pmfs indicated by $mL,mM$ and $mR$ in the figure are indexed by the parameter values $\beta=(1,0)$, $\beta=(0,1)$ and $\beta=(0,0)$, respectively.

C.3.2 The stick-breaking parameterization of the categorical distribution

Alternatively, we could also have used a “stick-breaking” parameter functional $\tilde{\nu}$ that sends each pmf in $\mathcal{P}$ to the vector of parameters $\tilde{\nu}(p(X))={\gamma_{1}\choose\gamma_{2}}$, where $\gamma_{1}=p_{L}$ and $\gamma_{2}=p_{M}/(1-p_{L})$.222222This only works if $p_{L}<1$. When $p(x_{1})=1$, we simply set $\gamma_{2}=0$, thus, $\gamma=(1,0)$. Again the parameter $\gamma={\gamma_{1}\choose\gamma_{2}}$ is only a label, but this time the range of $\tilde{\nu}$ is the parameter space $\Gamma=[0,1]\times[0,1]$. The functional relationship $f$ associated to $\gamma$ is given by

[TABLE]

For each $\gamma$ we can transform the pmf into the vector

[TABLE]

and write $\mathcal{M}_{\Gamma}$ for the collection of vectors so defined. As before, this collection coincides with the full model, i.e., $\mathcal{M}_{\Gamma}=\mathcal{M}$. In other words, by virtue of the functional $\tilde{\nu}$ and its inverse $\gamma\mapsto p_{\gamma}(x)=f(x\,|\,\gamma)$ we conclude that the parameter space $\Gamma$ and the complete set of pmfs $\mathcal{M}$ are isomorphic. Because $\mathcal{M}=\mathcal{M}_{B}$ this means that we also have an isomorphism between the parameter space $B$ and $\Gamma$ via $\mathcal{M}$, even though $B$ is a strict subset of $\Gamma$. Note that this equivalence goes via parameterization $\beta\mapsto m_{\beta}(X)$ and the functional $\tilde{\nu}$.

C.3.3 Multidimensional Jeffreys’s prior via the Fisher information matrix and orthogonal parameters

The multidimensional Jeffreys’s prior is parameterization-invariant and has as normalization constant $V=\int\sqrt{\det I_{X}(\theta)}\textnormal{d}\theta$, where $\det I_{X}(\theta)$ is the determinant of the Fisher information matrix.

In the previous subsection we argued that the categorical distribution in terms of $\beta$ or parameterized with $\gamma$ are equivalent to each other, that is, $\mathcal{M}_{B}=\mathcal{M}=\mathcal{M}_{\Gamma}$. However, these two parameterizations describe the model space $\mathcal{M}$ quite differently. In this subsection we use the Fisher information to show that the parameterization in terms of $\gamma$ is sometimes preferred over $\beta$.

The complete model $\mathcal{M}$ is easier described by $\gamma$, because the parameters are orthogonal. We say that two parameters are orthogonal to each other whenever the corresponding off-diagonal entries in the Fisher information matrix are zero. The Fisher information matrices in terms of $\beta$ and $\gamma$ are

[TABLE]

respectively.

The left panel of Fig. 13 shows the tangent vectors at $p_{\beta^{*}}(X)=[1/3,1/3,1/3]$ in model space, where $\beta^{*}=(1/3,1/3)$. The green tangent vector corresponds to $\frac{\partial m_{\beta^{*}}(X)}{\partial\beta_{1}}$, thus, with $\beta_{2}=1/3$ fixed and $\beta_{1}$ free to vary, while the red tangent vector corresponds to $\frac{\partial m_{\beta^{*}}(X)}{\partial\beta_{2}}$, thus, with $\beta_{1}=1/3$ and $\beta_{2}$ free to vary. The area of the diamond spanned by these two tangent vectors is $\sqrt{\det I(\beta^{*})}\textnormal{d}\beta_{1}\textnormal{d}\beta_{2}$, where we have taken $\textnormal{d}\beta_{1}=0.1$ and $\textnormal{d}\beta_{2}=0.1$.

The right panel of Fig. 13 shows the tangent vectors at the same point $p_{\gamma^{*}}(X)=[1/3,1/3,1/3]$, where $\gamma^{*}=(1/3,1/2)$. The green tangent vector corresponds to $\frac{\partial m_{\gamma^{*}}(X)}{\partial\gamma_{1}}$, thus, with $\gamma_{2}=1/2$ fixed and $\gamma_{1}$ free to vary, while the red tangent vector corresponds to $\frac{\partial m_{\gamma^{*}}(X)}{\partial\gamma_{2}}$, thus, with $\gamma_{1}=1/3$ and $\gamma_{2}$ free to vary. By glancing over the plots, we see that the two tangent vectors are indeed orthogonal. The area of the rectangle spanned by these these two tangent vectors is $\sqrt{\det I(\gamma^{*})}\textnormal{d}\gamma_{1}\textnormal{d}\gamma_{2}$, where we have taken $\textnormal{d}\gamma_{1}=\textnormal{d}\gamma_{2}=0.1$.

There are now two ways to calculate the normalizing constant of the Jeffreys’s prior, the area, more generally volume, of the model $\mathcal{M}$. In terms of $\beta$ this leads to

[TABLE]

Observe that the inner integral depends on the value of $\beta_{1}$ from the outer integral. This coupling is reflected by the non-zero off-diagonal term of the Fisher information matrix $I_{X}(\beta)$ corresponding to $\beta_{1}$ and $\beta_{2}$. On the other hand, orthogonality implies that the two parameters can be treated independently of each other. That is, knowing and fixing $\gamma_{1}$ and changing $\gamma_{2}$ will not affect $m_{\gamma}(X)$ via $\gamma_{1}$. This means that the double integral decouples

[TABLE]

Using standard geometry we verify that this is indeed the area of $\mathcal{M}$, as an eighth of the surface area of a sphere of radius two is given by $\tfrac{1}{8}4\pi 2^{2}=2\pi$.

Orthogonality is relevant in Bayesian analysis, as it provides an argument to choose a prior on a vector-valued parameter that factorizes (e.g., Berger, Pericchi and Varshavsky, 1998; Huzurbazar, 1950, 1956; Jeffreys, 1961; Kass and Vaidyanathan, 1992; Ly, Verhagen and Wagenmakers, 2016a, b), see also Cox and Reid (1987); Mitchell (1962).

By taking a random variable $X$ with $w=3$ outcomes, we were able to visualize the geometry of model space. For more general $X$ these plots get more complicated and perhaps even impossible to draw. Nonetheless, the ideas conveyed here extend, even to continuous $X$, whenever the model adheres to the regularity conditions given in Appendix E.

Appendix D MDL: Coding Theoretical Background

D.1 Coding theory, code length and log-loss

A coding system translates words, i.e., outcomes of a random variable $X$, into code words with code lengths that behave like a pmf. Code lengths can be measured with a logarithm, which motivates the adoption of log-loss, defined below, as the decision criterion within the MDL paradigm. The coding theoretical terminologies introduced here are illustrated using the random variable $X$ with $w=3$ potential outcomes.

D.1.1 Kraft-McMillan inequality: From code lengths of a specific coding system to a pmf

For the source-memory task we encoded the outcomes as $L,M$ and $R$, but when we communicate a participant’s responses $x^{n}_{\textnormal{obs}}$ to a collaborator over the internet, we have to encode the observations $x^{n}_{\textnormal{obs}}$ as zeroes and ones. For instance, we might use a coding system $\tilde{C}$ with code words $\tilde{C}(X=L)=00$, $\tilde{C}(X=M)=01$ and $\tilde{C}(X=R)=10$. This coding system $\tilde{C}$ will transform any set of responses $x^{n}_{\textnormal{obs}}$ into a code string $\tilde{C}(x^{n}_{\textnormal{obs}})$ consisting of $2n$ bits. Alternatively, we can use a coding system $C$ with code words $C(X=L)=10$, $C(X=M)=0$ and $C(X=R)=11$, instead. Depending on the actual observations $x^{n}_{\textnormal{obs}}$, this coding system outputs code strings $C(x^{n}_{\textnormal{obs}})$ with varying code lengths that range from $n$ to $2n$ bits. For example, if a participant responded with $x^{n}_{\textnormal{obs}}=(M,R,M,L,L,M,M,M)$ in $n=8$ trials, the coding system $C$ would then output the 11-bit long code string $C(x^{n}_{\textnormal{obs}})=01101010000$. In contrast, the first coding system $\tilde{C}$ will always output a $16$-bit long code string when $n=8$. Shorter code strings are desirable as they will lead to a smaller load on the communication network and they are less likely to be intercepted by “competing” researchers.

Note that the shorter code length $C(x^{n}_{\textnormal{obs}})=01101010000$ of 11-bits is a result of having code words of unequal lengths. The fact that one of the code word is shorter does not interfere with the decoding, since no code word is a prefix of another code word. As such, we refer to $C$ as a prefix (free) coding system. This implies that the 11-bit long code string $C(x^{n}_{\textnormal{obs}})$ is self-punctuated and that it can be uniquely deciphered by simply reading the code string from left to right resulting in the retrieval of $x^{n}_{\textnormal{obs}}$. Note that the code lengths of $C$ inherit the randomness of the data. In particular, the coding system $C$ produces a shorter code string with high chance, if the participant generates the outcome $M$ with high chance. In the extreme case, the coding system $C$ produces the 8-bits long code string $C(x^{n})=00000000$ with 100% (respectively, 0%) chance, if the participant generates the outcome $M$ with 100% (respectively, 0%) chance. More generally, Kraft and McMillan (Kraft, 1949; McMillan, 1956) showed that any uniquely decipherable (prefix) coding system from the outcome space $\mathcal{X}$ with $w$ outcomes to an alphabet with $D$ elements must satisfy the inequality

[TABLE]

where $l_{i}$ is the code length of the outcome $w$. In our example, we have taken $D=2$ and code length of 2, 1 and 2 bits for the response $L,M$ and $R$ respectively. Indeed, $2^{-2}+2^{-1}+2^{-2}=1$. Hence, code lengths behave like the logarithm (with base $D$) of a pmf.

D.1.2 Shannon-Fano algorithm: From a pmf to a coding system with specific code lengths

Given a data generating pmf $p^{*}(X)$, we can use the so-called Shannon-Fano algorithm (e.g., Cover and Thomas, 2006, Ch. 5) to construct a prefix coding system $C^{*}$. The idea behind this algorithm is to give the outcome $x$ that is generated with the highest chance the shortest code length. To do so, we encode the outcome $x$ as a code word $C^{*}(x)$ that consists of $-\log_{2}p^{*}(x)$ bits.232323When we use the logarithm with base two, $\log_{2}(y)$, we get the code length in bits, while the natural logarithm, $\log(y)$, yields the code length in nats. Any result in terms of the natural logarithm can be equivalently described in terms of the logarithm with base two, as $\log(y)=\log(2)\log_{2}(y)$.

For instance, when a participant generates the outcomes $[L,M,R]$ according to the chances $p^{*}(X)=[0.25,0.5,0.25]$ the Shannon-Fano algorithm prescribes that we should encode the outcome $L$ with $-\log_{2}(0.25)=2$, $M$ with $-\log_{2}(0.5)=1$ and $R$ with 2 bits; the coding system $C$ given above.242424Due to rounding, the Shannon-Fano algorithm actually produces code words $C(x)$ that are at most one bit larger than the ideal code length $-\log_{2}p^{*}(x)$. We avoid further discussions on rounding. Moreover, in the following we consider the natural logarithm instead. The Shannon-Fano algorithm works similarly for any other given pmf $p_{\beta}(X)$. Hence, the Kraft-McMillan inequality and its inverse, i.e., the Shannon-Fano algorithm imply that pmfs and uniquely decipherable coding systems are equivalent to each other. As such we have an additional interpretation of a pmf. To distinguish the different uses, we write $f(X\,|\,\beta)$ when we view the pmf as a coding system, while we retain the notation $p_{\beta}(X)$ when we view the pmf as a data generating device. In the remainder of this section we will not explicitly construct any other coding system, as the coding system itself is irrelevant for the discussion at hand –only the code lengths matter.

D.1.3 Entropy, cross entropy, log-loss

With the true data generating pmf $p^{*}(X)$ at hand, thus, also the true coding system $f(X\,|\,\beta^{*})$, we can calculate the (population) average code length per trial

[TABLE]

Whenever we use the logarithm with base 2, we refer to this quantity $H(p^{*}(X))$ as the Shannon entropy.252525Shannon denoted this quantity with an $H$ to refer to the capital Greek letter for eta. It seems that John von Neumann convinced Claude Shannon to call this quantity entropy rather than information (Tribus and McIrvine, 1971). If the true pmf is $p^{*}(X)=[0.25,0.5,0.25]$ we have an average code length of $1.5$ bits per trail whenever we use the true coding system $f(X\,|\,\beta^{*})$. Thus, we expect to use 12 bits to encode observations consisting of $n=8$ trials.

As coding theorists, we have no control over the true data generating pmf $p^{*}(X)$, but we can choose the coding system $f(X\,|\,\beta)$ to encode the observations. The (population) average code length per trial is given by

[TABLE]

The quantity $H(p^{*}(X)\|\,\beta)$ is also known as the cross entropy from the true pmf $p^{*}(X)$ to the postulated $f(X\,|\,\beta)$.262626Observe that the entropy $H(p^{*}(X))$ is the just the cross entropy from the true $p^{*}(X)$ to the true coding system $f(X\,|\,\beta^{*})$. For instance, when we use the pmf $f(X\,|\,\beta)=[0.01,0.18,0.81]$ to encode data that are generated according to $p^{*}(X)=[0.25,0.5,0.25]$, we will use $2.97$ bits on average per trial. Clearly, this is much more than the $1.5$ bits per trial that we get from using the true coding system $f(X\,|\,\beta^{*})$.

More generally, Shannon (1948) showed that the cross entropy can never be smaller than the entropy, i.e., $H(p^{*}(X))\leq H(p^{*}(X)\|\,\beta)$. In other words, we always get a larger average code length, whenever we use the wrong coding system $f(X\,|\,\beta)$. To see why this holds, we decompose the cross entropy as a sum of the entropy and the Kullback-Leibler divergence,272727The KL-divergence is also known as the relative entropy. and show that the latter cannot be negative. This decomposition follows from the definition of cross entropy and a subsequent addition and subtraction of the entropy resulting in

[TABLE]

where $D(p^{*}(X)\|\beta)$ defines the Kullback-Leibler divergence from the true pmf $p^{*}(X)$ to the postulated coding system $f(X\,|\,\beta)$. Using the so-called Jensen’s inequality it can be shown that the KL-divergence is non-negative and that it is only zero whenever $f(X\,|\,\beta)=p^{*}(X)$. Thus, the cross entropy can never be smaller than the entropy. Consequently, to minimize the load on the communication network, we have to minimize the cross entropy with respect to the parameter $\beta$. Unfortunately, however, we cannot do this in practice, because the cross entropy is a population quantity based on the unknown true pmf $p^{*}(X)$. Instead, we do the next best thing by replacing the true $p^{*}(X)$ in Eq. (D.3) by the empirical pmf that gives the relative occurrences of the outcomes in the sample rather than in the population. Hence, for any postulated $f(X\,|\,\beta)$, with $\beta$ fixed, we approximate the population average defined in Eq. (D.3) by the sample average

[TABLE]

We call the quantity $H(x^{n}_{\textnormal{obs}}\|\,\beta)$ the log-loss from the observed data $x^{n}_{\textnormal{obs}}$, i.e., the empirical pmf $\hat{p}_{\textnormal{obs}}(X)$, to the coding system $f(X\,|\,\beta)$.

D.2 Data compression and statistical inference

The entropy inequality $H(p^{*}(X))\leq H(p^{*}(X)\|\beta)$ implies that the coding theorist’s goal of finding the coding system $f(X\,|\,\beta)$ with the shortest average code length is in fact equivalent to the statistical goal of finding the true data generating process $p^{*}(X)$. The coding theorist’s best guess is the coding system $f(X\,|\,\beta)$ that minimizes the log-loss from $x^{n}_{\textnormal{obs}}$ to the model $\mathcal{M}_{B}$. Note that minimizing the negative log-likelihood is the same as maximizing the likelihood. Hence, the log-loss is minimized by the coding system associated with the MLE, thus, the predictive pmf $f(X\,|\,\hat{\beta}_{\textnormal{obs}})$. Furthermore, the cross entropy decomposition shows that minimization of the log-loss is equivalent to minimization of the KL-divergence from the observations $x^{n}_{\textnormal{obs}}$ to the model $\mathcal{M}_{B}$. The advantage of having the optimization problem formulated in terms of KL-divergence is that it has a known lower bound, namely, zero. Moreover, whenever the KL-divergence from $x^{n}_{\textnormal{obs}}$ to the code $f(X\,|\,\hat{\beta}_{\textnormal{obs}})$ is larger than zero, we then know that the empirical pmf associated to the observations does not reside on the model. In particular, Section 4.3.1 showed that the MLE plugin, $f(X\,|\,\hat{\beta}_{\textnormal{obs}})$ is the pmf on the model that is closest to the data. This geometric interpretation is due to the fact that we retrieve the Fisher-Rao metric, when we take the second derivative of the KL-divergence with respect to $\beta$ (Kullback and Leibler, 1951). This connection between the KL-divergence and Fisher information is exploited in Ghosal, Ghosh and Ramamoorthi (1997) to generalize the Jeffreys’s prior to nonparametric models, see also van Erven and Harremos (2014) for the relationship between KL-divergence and the broader class of divergence measures developed by Rényi (1961), see also Campbell (1965).

Appendix E Regularity conditions

A more mathematically rigorous exposition of the subject would have had this section as the starting point, rather than the last section of the appendix. The regularity conditions given below can be seen as a summary, and guidelines for model builders. If we as scientists construct models such that these conditions are met, we can then use the results presented in the main text. We first give a more general notion of statistical models, then state the regularity conditions followed by a brief discussion on these conditions.

The goal of statistical inference is to find the true probability measure $P^{*}$ that governs the chances with which $X$ takes on its events. A model $\mathscr{P}_{\Theta}$ defines a subset of $\mathscr{P}$, the largest collection of all possible probability measures. We as model builders choose $\mathscr{P}_{\Theta}$ and perceive each probability measure $P$ within $\mathscr{P}_{\Theta}$ as a possible explanation of how the events of $X$ were or will be generated. When $P^{*}\in\mathscr{P}_{\Theta}$ we have a well-specified model and when $P^{*}\nin\mathscr{P}_{\Theta}$, we say that the model is misspecified.

By taking $\mathscr{P}_{\Theta}$ to be equal to the largest possible collection $\mathscr{P}$, we will not be misspecified. Unfortunately, this choice is not helpful as the complete set is hard to track and leads to uninterpretable inferences. Instead, we typically construct the candidate set $\mathscr{P}_{\Theta}$ using a parameterization that sends a label $\theta\in\Theta$ to a probability measure $P_{\theta}$. For instance, we might take the label $\theta={\mu\choose\sigma^{2}}$ from the parameter space $\Theta=\mathbb{R}\times(0,\infty)$ and interpret these two numbers as the population mean and variance of a normal probability $P_{\theta}$. This distributional choice is typical in psychology, because it allows for very tractable inference with parameters that are generally overinterpreted. Unfortunately, the normal distribution comes with rather stringent assumptions resulting in a high risk of misspecification. More specifically, the normal distribution is far too ideal, as it supposes that the population is nicely symmetrically centred at its population mean and outliers are practically not expected due to its tail behavior.

Statistical modeling is concerned with the intelligent construction of the candidate set $\mathscr{P}_{\Theta}$ such that it encapsulates the true probability measure $P^{*}$. In other words, the restriction of $\mathscr{P}$ to $\mathscr{P}_{\Theta}$ in a meaningful manner. Consequently, the goal of statistical inference is to give an informed guess $\tilde{P}$ within $\mathscr{P}_{\Theta}$ for $P^{*}$ based on the data. This guess should give us insights to how the data were generated and how yet unseen data will be generated. Hence, the goal is not to find the parameters as they are mere labels. Of course parameters can be helpful, but they should not be the goal of inference.

Note that our general description of a model as a candidate set $\mathscr{P}_{\Theta}$ does not involve any structure –thus, the members of $\mathscr{P}_{\Theta}$ do not need to be related to each other in any sense. We use the parameterization to transfer the structure of our labels $\Theta$ to a structure on $\mathscr{P}_{\Theta}$. To do so, we require that $\Theta$ is a nice open subset of $\mathbb{R}^{d}$. Furthermore, we require that each label defines a member $P_{\theta}$ of $\mathscr{P}_{\Theta}$ unambiguously. This means that if $\theta^{*}$ and $\theta$ differ from each other that the resulting pair of probability measure $P_{\theta^{*}}$ and $P_{\theta}$ also differ from each other. Equivalently, we call a parameterization identifiable whenever $\theta^{*}=\theta$ leads to $P_{\theta^{*}}=P_{\theta}$. Conversely, identifiability implies that when we know everything about $P_{\theta}$, we can then also use the inverse of the parameterization to pinpoint the unique $\theta$ that corresponds to $P_{\theta}$. We write $\nu:\mathcal{P}_{\Theta}\rightarrow\Theta$ for the functional that attaches to each probability measure $P$ a label $\theta$. For instance, $\nu$ could be defined on the family of normal distribution such that $P\mapsto\nu(P)={E_{P}(X)\choose\text{Var}_{P}(X)}={\mu\choose\sigma^{2}}$. In this case we have $\nu(\mathscr{P}_{\Theta})=\Theta$ and, therefore, a one-to-one correspondence between the probability measures $P_{\theta}\in\mathscr{P}_{\Theta}$ and the parameters $\theta\in\Theta$.

By virtue of the parameterization and its inverse $\nu$, we can now transfer additional structure from $\Theta$ to $\mathscr{P}_{\Theta}$. We assume that each probability measure $P_{\theta}$ that is defined on the events of $X$ can be identified with a probability density function (pdf) $p_{\theta}(x)$ that is defined on the outcomes of $X$. For this assumption, we require that the set $\mathscr{P}_{\Theta}$ is dominated by a so-called countably additive measure $\lambda$. When $X$ is continuous, we usually take for $\lambda$ the Lebesgue measure that assigns to each interval of the form $(a,b)$ a length of $b-a$. Domination allows us to express the probability of $X$ falling in the range $(a,b)$ under $P_{\theta}$ by the “area under the curve of $p_{\theta}(x)$”, that is, P_{\theta}\big{(}X\in(a,b)\big{)}=\int_{a}^{b}p_{\theta}(x)\textnormal{d}x. For discrete variables $X$ taking values in $\mathcal{X}=\{x_{1},x_{2},x_{3},\ldots\}$, we take $\lambda$ to be the counting measure. Consequently, the probability of observing the event $X\in A$ where $A=\{a=x_{1},x_{2},\ldots,b=x_{k}\}$ is calculated by summing the pmf at each outcome, that is, $P_{\theta}(X\in A)=\sum_{x=a}^{x=b}p_{\theta}(x)$. Thus, we represent $\mathscr{P}_{\Theta}$ as the set $\mathcal{P}_{\Theta}=\{p_{\theta}(x)\,:\,\theta\in\Theta,P_{\theta}(x)=\int_{-\infty}^{x}p_{\theta}(y)\textnormal{d}y\text{ for all x}\in\mathcal{X}\}$ in function space. With this representation of $\mathcal{P}_{\Theta}$ in function space, the parameterization is now essentially the functional relationship $f$ that pushes each $\theta$ in $\Theta$ to a pdf $p_{\theta}(x)$. If we choose $f$ to be regular, we can then also transfer additional topological structure from $\Theta$ to $\mathcal{P}_{\Theta}$.

Definition E.1** (Regular parametric model).**

We call the model $\mathcal{P}_{\Theta}$ a regular parametric model, if the parameterization $\theta\mapsto p_{\theta}(x)=f(x\,|\,\theta)$, that is, the functional relationship $f$, satisfies the following conditions

1. (i)

its domain $\Theta$ is an open subset of $\mathbb{R}^{d}$, 2. (ii)

at each possible true value $\theta^{*}\in\Theta$, the spherical representation $\theta\mapsto m_{\theta}(x)=2\sqrt{p_{\theta}(x)}=2\sqrt{f(x\,|\,\theta)}$ is so-called Fréchet differentiable in $L_{2}(\lambda)$. The tangent function, i.e., the “derivative” in function space, at $m_{\theta^{*}}(x)$ is then given by

[TABLE]

where $\dot{l}(x\,|\,\theta^{*})$ is a $d$-dimensional vector of score functions in $L_{2}(P_{\theta^{*}})$, 3. (iii)

the Fisher information matrix $I_{X}(\theta)$ is non-singular, 4. (iv)

the map $\theta\mapsto\dot{l}(x\,|\,\theta)m_{\theta}(x)$ is continuous from $\Theta$ to $L_{2}^{d}(\lambda)$.

Note that (ii) allows us to generalize the geometrical concepts discussed in Appendix C.3 to more general random variables $X$. $\hfill\diamond$

We provide some intuition. Condition (i) implies that $\Theta$ inherits the topological structure of $\mathbb{R}^{d}$. In particular, we have an inner product on $\mathbb{R}^{d}$ that allows us to project vectors onto each other, a norm that allows us to measure the length of a vector, and the Euclidean metric that allows us to measure the distance between two vectors by taking the square root of the sums of squares, that is, \|\theta^{*}-\theta\|_{2}=\sqrt{\sum_{i=1}^{d}\big{(}\theta^{*}_{i}-\theta_{i}\big{)}^{2}}. For $d=1$ this norm is just the absolute value, which is why we previously denoted this as $|\theta^{*}-\theta|$.

Condition (ii) implies that the measurement of distances in $\mathbb{R}^{d}$ generalizes to the measurement of distance in function space $L_{2}(\lambda)$. Intuitively, we perceive functions as vectors and say that a function $h$ is a member of $L_{2}(\lambda)$, if it has a finite norm (length), i.e., $\|h(x)\|_{L_{2}(\lambda)}<\infty$, meaning

[TABLE]

As visualized in the main text, by considering $\mathcal{M}_{\Theta}=\{m_{\theta}(x)=\sqrt{p_{\theta}(x)}\,|\,p_{\theta}\in\mathcal{P}_{\theta}\}$ we relate $\Theta$ to a subset of the sphere with radius two in the function space $L_{2}(\lambda)$. In particular, Section 4 showed that whenever the parameter is one-dimensional, thus, a line, that the resulting collection $\mathcal{M}_{\Theta}$ also defines a line in model space. Similarly, Appendix C.3 showed that whenever the parameter space is a subset of $[0,1]\times[0,1]$ that the resulting $\mathcal{M}_{\Theta}$ also forms a plain.

Fréchet differentiability at $\theta^{*}$ is formalized as

[TABLE]

This implies that the linearization term $\tfrac{1}{2}(\theta-\theta^{*})^{T}\dot{l}(x\,|\,\theta^{*})m_{\theta^{*}}(x)$ is a good approximation to the “error” $m_{\theta}(x)-m_{\theta^{*}}(x)$ in the model $\mathcal{M}_{\Theta}$, whenever $\theta$ is close to $\theta^{*}$ given that the score functions $\dot{l}(x\,|\,\theta^{*})$ do not blow up. More specifically, this means that each component of $\dot{l}(x\,|\,\theta^{*})$ has a finite norm. We say that the component $\tfrac{\partial}{\partial\theta_{i}}l(x\,|\,\theta^{*})$ is in $L_{2}(P_{\theta^{*}})$, if $\|\tfrac{\partial}{\partial\theta_{i}}l(x\,|\,\theta^{*})\|_{L_{2}(P_{\theta^{*}})}<\infty$, meaning

[TABLE]

This condition is visualized in Fig. 12 and Fig. 13 by tangent vectors with finite lengths. Under $P_{\theta^{*}}$, each component $i=1,\ldots,d$ of the tangent vector is expected to be zero, that is,

[TABLE]

This condition follows from the chain rule applied to the logarithm and an exchange of the order of integration with respect to $x$, and derivation with respect to $\theta_{i}$, as

[TABLE]

Note that if $\int\tfrac{\partial}{\partial\theta_{i}}p_{\theta^{*}}(x)\textnormal{d}x>0$, then a small change at $\theta^{*}$ will lead to a function $p_{\theta^{*}+\textnormal{d}\theta}(x)$ that does not integrate to one and, therefore, not a pdf.

Condition (iii) implies that the model does not collapse to a lower dimension. For instance, when the parameter space is a plain the resulting model $\mathcal{M}_{\Theta}$ cannot be line. Lastly, condition (iv) implies that the tangent functions change smoothly as we move from $m_{\theta^{*}}(x)$ to $m_{\theta}(x)$ on the sphere in $L_{2}(\lambda)$, where $\theta$ is a parameter value in the neighborhood of $\theta^{*}$.

The following conditions are stronger, thus, less general, but avoid Fréchet differentiability and are typically easier to check.

Lemma E.1**.**

Let $\Theta\subset\mathbb{R}^{d}$ be open. At each possible true value $\theta^{*}\in\Theta$, we assume that $p_{\theta}(x)$ is continuously differentiable in $\theta$ for $\lambda$-almost all $x$ with tangent vector $\dot{p}_{\theta^{*}}(x)$. We define the score function at $x$ as

[TABLE]

where $1_{[p_{\theta^{*}}>0]}(x)$ is the indicator function

[TABLE]

The parameterization $\theta\mapsto P_{\theta}$ is regular, if the norm of the score vector Eq. (E.7) is finite in quadratic mean, that is, $\dot{l}(X\,|\,\theta^{*})\in L_{2}(P_{\theta^{*}})$, and if the corresponding Fisher information matrix based on the score functions Eq. (E.7) is non-singular and continuous in $\theta$. $\hfill\diamond$

There are many better sources than the current manuscript on this topic that are mathematically much more rigorous and better written. For instance, Bickel et al. (1993) give a proof of the lemma above and many more beautiful, but sometimes rather (agonizingly) technically challenging, results. For a more accessible, but no less elegant, exposition of the theory we highly recommend van der Vaart (1998).