Introduction. Statistics and Estimators

WarningDraft

This set of notes is still being edited: the text may contain typos, inaccuracies, and locally unpolished spots. If you find something, please let the author know (contact information on the course page).

1 Mathematical Statistics: Introduction

1.1 Probability Theory VS Mathematical Statistics

In probability theory, we assumed that we knew the nature of some phenomenon or experiment (the distribution of a random variable), and this knowledge allowed us to study the behavior of various derived quantities. Given the distribution of some random variable \(\xi\), we could, for example, find the distribution of \(\xi^2\), find how the expectations, variances, and distributions of sums of dependent/independent random variables relate to one another, and so on — that is, we went from a known distribution to the properties and behavior of quantities related to it.

In mathematical statistics everything is reversed: the data are the results of experiments (values of a random variable), and we need to understand the nature of the phenomenon being studied (the distribution of the quantity).

Let \(\xi\) be a random variable taking values in the set \(\mathbb {R}\) and having distribution \(\mathbb {P}\)¹. We will call the random variable \(\xi\) an observation. One of the main problems of mathematical statistics, which is essentially what we will be doing, is the following:

Based on a set of observations, determine some information about the distribution \(\mathbb{P}\), which is assumed to be unknown.

As a rule, it is possible to carry out independent experiments in order to determine the nature of the phenomenon under study. Speaking mathematically, the distribution needs to be recovered from a set of independent observations, called a sample. That is, a sample is a vector \(\left(X_{1}, \ldots , X_{n}\right)^T\) (or simply \(X_{1}, \ldots , X_{n}\)) made up of independent random variables distributed according to \(\mathbb {P}\), i.e., in the same way as \(\xi\). Such a vector takes values in the set \(\mathbb {R}^{n}\), and its distribution \(\mathbb {P}^{n}\) is defined on the Borel \(\sigma\)-algebra \(\mathscr {B}\left(\mathbb {R}^{n}\right)\).

There are several approaches to the problem of estimating a distribution. Let us consider the main ones.

1.2 Parametric VS Nonparametric Approaches

In real experiments, as a rule, it is possible to narrow the possible values of \(\mathbb {P}\) somewhat before carrying out the study, i.e., to find some family of distributions \(\mathcal{P}\) to which \(\mathbb {P}\) is guaranteed to belong. The triple \((\mathbb {R}, \mathcal{B}(\mathbb {R}), \mathcal{P})\) is called a probabilistic-statistical model.

In the parametric approach, we assume that the family of distributions \(\mathcal{P}\) to which our unknown distribution \(\mathbb {P}\) belongs consists only of distributions given by some parameter.

ExampleExample 1

Suppose we want to estimate the rate at which customers visit a store. Let \(\xi\) be the time between visits of consecutive customers, that is, how much time passes between the visit of one customer and the next. Let \(\mathbb {P}\) be the distribution of \(\xi\).

It is clear that \(\xi\) is a nonnegative random variable with a continuous distribution, i.e., it is known in advance that \(\mathbb {P} \in \mathcal{P}\), where \(\mathcal{P}\) is the class of all continuous distributions with support \([0, +\infty )\). We need to estimate the density \(f_0\) of the distribution \(\mathbb {P}\).

To estimate \(\mathbb {P}\), the following sample (more precisely, its realization) was collected: \[ \begin{pmatrix} X_1(\omega ) \\ X_2(\omega ) \\ \vdots \\ X_{20}(\omega ) \end{pmatrix} = \left(\begin{smallmatrix} 0.5 \\ 0.8 \\ 1.2 \\ 5.8 \\ 8.5 \\ 7.2 \\ 1.5 \\ 1.8 \\ 2.0 \\ 2.3 \\ 2.6 \\ 3.5 \\ 3.9 \\ 4.5 \\ 5.2 \\ 6.4 \\ 6.8 \\ 7.6 \\ 9.2 \\ 9.8 \end{smallmatrix}\right) \] (time is measured in minutes). Estimate \(\mathbb {P}\) using the parametric and nonparametric approaches.

SolutionSolution

Enum-item1.

Parametric approach. We need to narrow \(\mathcal{P}\), that is, make an additional assumption about the form of \(\mathbb {P}\). It is common to model the time between recurring events using the exponential distribution, since it has the memoryless property. So, we set \(\mathcal{P} = \{ \operatorname {Exp}(\theta ) : \theta > 0\}\), i.e., \(\xi \sim \mathbb {P} = \operatorname {Exp}(\theta )\) for some unknown \(\theta > 0\). Density:

\[ f(x) = f(x,\theta ) = \theta e^{-\theta x} \cdot \; \mathbb {1}_{[0,+\infty )}(x) \]

The problem of estimating the density has been reduced to constructing an estimator of the parameter \(\theta\).

We have \(\mathbb {E}\left[\xi \right] = 1/\theta\), i.e., \(\theta = \frac{1}{\mathbb {E}\left[\xi \right]}\). It is reasonable to estimate the expectation \(\mathbb {E}\left[\xi \right]\) by the sample mean \(\frac{1}{n}\sum_{k=1}^{n}X_k\) (we will discuss this in more detail later). Hence, a natural estimator of \(\theta\) is

\[ \widehat{\theta } = \frac{1}{\frac{1}{n}\sum _{k=1}^{n}X_k}. \]

In our specific case:

\[ \frac{1}{n}\sum _{k=1}^{n}X_k(\omega _0) = \frac{0.5 + 0.8 +\ldots + 9.8}{20} = \frac{91.1}{20} = 4.555, \]

whence \(\widehat{\theta }(\omega_0) \approx 0.220\).

Parametric density estimate:

\[ \widehat{f}^{\mathrm{P}}(x) = \widehat{\theta }\, e^{-\widehat{\theta }\, x} \cdot \; \mathbb {1}_{[0,+\infty )}(x), \]

where P stands for parametric. In our specific case:

\[ \widehat{f}^{\mathrm{P}}(x, \omega _0) \approx 0.220 \cdot e^{-0.220 \cdot x} \cdot \; \mathbb {1}_{[0,+\infty )}(x) \]

Enum-item2.

Nonparametric approach. Here we make no assumptions about the shape of the distribution. Let us use the following method: over each sample point \(X_i\) we inflate a symmetric “bump” of area \(\frac{1}{n}\). The shape of the bump is given by a kernel \(K\) — a function satisfying \(K(u) \geq 0\), \(\int K(u)\, du = 1\), \(K(u) = K(-u)\). Popular kernels:

• the Epanechnikov kernel \(K(u) = \frac{3}{4}(1 - u^2) \cdot \; \mathbb {1}_{[-1,1]}(u)\);

• the Gaussian kernel \(K(u) = \frac{1}{\sqrt{2\pi }} e^{-u^2/2}\);

• the uniform kernel \(K(u) = \frac{1}{2}\cdot \; \mathbb {1}_{[-1,1]}(u)\).

The width of the bump is controlled by the parameter \(h > 0\) (the bandwidth). Preliminary density estimate:

\[ \widehat{f}_0^{\mathrm{NP}}(x) = \sum _{i=1}^{n} \frac{1}{n} \cdot \frac{1}{h} \cdot K\! \left(\frac{x - X_i}{h}\right). \]

where \(\frac{1}{h} \cdot K\! \left(\frac{x - X_i}{h}\right)\) is a bump of width \(h\) centered at \(X_i\). The factor \(\frac{1}{h}\) is included so that the area under the graph remains equal to \(1\).

This is the so-called kernel density estimate, abbreviated KDE (Kernel Density Estimation).

Since the kernel \(K\) is symmetric, part of the area under \(\widehat{f}_0^{\mathrm{NP}}\) sometimes “leaks” onto the negative half-axis, which contradicts the fact that \(\xi \geq 0\). Denote this leaked area by

\[ s = \int _{-\infty }^{0} \widehat{f}_0^{\mathrm{NP}}(x)\, dx. \]

Let us correct the estimate by renormalizing it on \([0, +\infty )\):

\[ \widehat{f}^{\mathrm{NP}}(x) = \begin{cases} \dfrac {1}{1 - s}\, \widehat{f}_0^{\mathrm{NP}}(x), & x \geq 0, \\ 0, & x < 0, \end{cases} \]

where NP stands for non-parametric. Then \(\int_0^{+\infty } \widehat{f}^{\mathrm{NP}}(x)\, dx = 1\).

Density estimation: nonparametric and parametric approaches

Orange points are the observations. Drag observations to move them. Click on an empty spot to add an observation. Double-click or right-click to remove an observation.

Bandwidth \(h\) 1.7 Kernel Epanechnikov Gaussian Uniform

\(\hat{f}_0^{\mathrm{NP}}\) — preliminary KDE (all \(x\))

\(\hat{f}^{\mathrm{NP}}\) — corrected KDE (\(x \ge 0\))

\(\hat{f}^{\mathrm{P}}\) — \(\mathrm{Exp}(\hat{\theta})\), \(\hat{\theta} = \bar{X}\)

Sample size: \(\bar{X} =\) \(\hat{\theta} =\) \(s =\)

As the graph shows, the density estimates differ substantially. The parametric estimate is coarser — it does not account for the fact that the distribution has 2 peaks: near the point \(2\) and near the point \(7\). On the other hand, the parametric estimate is simpler.

1.3 Frequentist VS Bayesian Approaches

Recall that in probability theory there are several ways to understand the word “probability.” A simple example: you flip a coin, it lands on your hand, but you immediately cover it with your other hand. What is the probability that it came up heads?

• Frequentist interpretation of probability. The probability of an event is the limit of the relative frequency of its occurrence under unlimited repetition of the same random experiment:

  \[ \mathbb {P}(A) = \lim _{n \to \infty } \frac{\text{number of successes}}{n}. \]

  To speak of the probability of something, one needs to imagine (at least mentally) a series of repetitions of the experiment. Therefore, probabilistic statements only make sense for random, repeatable events.

  In the coin example, “probability” makes no sense from this point of view: the coin has already landed on your hand, it has already come up somehow. There is no repeatability here: you could open your hand 1000 times, and all 1000 times you would see the same result (though you don’t know it yet, it is nonetheless certain that it will be the same).

  When, within the frequentist approach, one says “a fair coin has probability \(\frac{1}{2}\) of landing heads,” what is meant is that if you flip the coin 1000 times, heads will come up in roughly half of the cases.

• Bayesian interpretation of probability. Probability is a degree of belief (degree of plausibility, degree of confidence) in the truth of some statement, based on available information. Such a probability can be assigned to any statement, including a one-time, non-repeatable event.

  From the Bayesian point of view, the question of the probability that an already-tossed coin landed heads or tails makes perfect sense: you are about \(1/2\) confident that it landed heads.

  Incidentally, this leads to a curious observation: an ordinary person, not versed in mathematical subtleties, when asked about the probability of heads for a coin that has already landed, will most likely confidently answer \(1/2\). Likewise, they would not be puzzled by the question “what is the probability that it will rain tomorrow,” even though there is no repeatability here and the frequentist approach does not apply. That is, the layperson’s interpretation of probability is, in essence, Bayesian.

  One more curious point: since in the Bayesian approach probability is a degree of belief, it can be updated as new, even indirect, information arrives. Suppose that initially (a priori) you had no information about the coin, and you were confident it was fair. Then, after the toss and covering it with your hand, the probability (degree of belief) of heads for you equals \(\frac{1}{2}\). Next, instead of uncovering your hand, you might try to determine the outcome by touch. Suppose it seems to you that you feel heads on top — then the probability (degree of belief) of heads might shift, say, to \(\frac{2}{3}\). Then, suppose it turns out that the coin was lent to you by a friend who says it is unfair and more prone to landing tails. After this new information, you can once again update your probability (degree of belief), shifting it, say, to \(\frac{3}{8}\). The probability obtained after taking new information into account is called the posterior probability. And so on: the Bayesian approach allows one to successively “update” probabilities as data arrive.

Let us consider the difference between the frequentist and Bayesian approaches using the example of estimating a parameter in a parametric model. In the frequentist approach, the parameter is considered a fixed but unknown number — it is not random, and it is formally incorrect to speak of “the probability that the parameter lies in some fixed interval”: either the parameter lies there, or it does not. Randomness in this approach is associated only with the data — with the fact that if the experiment is repeated many times, the sample will come out differently each time.

In the Bayesian approach, the parameter itself is declared to be a random variable. This does not mean that the experiment is “actually” carried out each time with a new parameter — it simply means that our knowledge about the parameter, before and after the experiment, is described by a probability distribution. Before we have seen the data, this knowledge is expressed by a prior distribution (from the Latin a priori — “from what is before”) — this is a distribution that we choose ourselves, based on what we know about the parameter beforehand (or based on convenience, if nothing is known). After the data are obtained, the prior distribution is recomputed into the posterior distribution (a posteriori — “from what is after”) using Bayes’ formula — this is the updated knowledge about the parameter, taking the observations into account.

ExampleExample 2

A coin was tossed \(n = 10\) times, and \(k = 8\) heads came up. Estimate the probability \(p\) of heads: a) using the frequentist approach (point estimate and confidence interval); b) using the Bayesian approach (posterior distribution and confidence interval).

SolutionSolution

a) Frequentist approach.

Maximum likelihood estimate:

\[ \hat{p} = \frac{k}{n} = \frac{8}{10} = 0{,}8. \]

Approximate 95% (Wald) confidence interval:

\[ \hat{p} \pm 1{,}96 \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} = 0{,}8 \pm 1{,}96 \cdot \sqrt{\frac{0{,}8 \cdot 0{,}2}{10}} \approx 0{,}8 \pm 0{,}248, \]

that is, \((0{,}55;\ 1{,}0)\) (after truncation at the boundary \(1\)).

Interpretation: if the experiment (10 tosses) is repeated many times and such an interval is constructed each time, then in 95 out of 100 cases it will cover the true value \(p\). Here \(p\) itself does not change and is not random — what is random is the interval itself (from sample to sample).

b) Bayesian approach.

We need to specify a prior distribution for \(p\). Let us take

\[ p \sim \mathrm{Beta}(\alpha =1,\ \beta =1), \]

that is, the uniform distribution on \([0,1]\) — it reflects a situation where we have no preferences beforehand: any value of \(p\) from 0 to 1 is considered a priori equally likely.

The likelihood of the data under the binomial model:

\[ L(p) \propto p^{k}(1-p)^{n-k}. \]

By Bayes’ formula, the posterior density is proportional to the product of the prior density and the likelihood; for the Beta-Binomial pair this is again a beta distribution:

\[ p \mid \text{data} \sim \mathrm{Beta}(\alpha + k,\ \beta + n - k) = \mathrm{Beta}(9, 3). \]

Posterior expectation:

\[ \mathbb {E}[p \mid \text{data}] = \frac{9}{9+3} = 0{,}75. \]

95% confidence interval (using the \(2{,}5\%\) and \(97{,}5\%\) quantiles of \(\mathrm{Beta}(9,3)\)):

\[ (0{,}48;\ 0{,}94). \]

Interpretation: after we have seen the data and updated our knowledge about \(p\), we can say — with probability 95%, the true value of \(p\) lies in this interval. Here this is a legitimate probabilistic statement about the parameter itself, because in this approach \(p\) is a random variable with posterior distribution \(\mathrm{Beta}(9,3)\).

Bayesian updating: the coin

Sliders for the data \(n\) (number of tosses) and \(k\) (number of heads) and for the prior \(\mathrm{Beta}(\alpha_0,\beta_0)\) parameters — the posterior density, confidence interval, and labels are recomputed instantly.

number of tosses \(n\) = 10

number of heads \(k\) = 8

\(\alpha_0\) (prior) = 1

\(\beta_0\) (prior) = 1

posterior: Beta(9, 3) posterior mean = 0.750 \(\hat p_{\mathrm{MLE}} = k/n\) = 0.800 95% Bayesian CI = (0.482, 0.940)

Comparison. The point estimate in the Bayesian approach (\(0{,}75\)) is shifted slightly toward the center of the prior distribution (\(0{,}5\)) compared to the purely empirical estimate \(0{,}8\) — this is precisely the influence of the prior distribution. If we take a more “confident” prior distribution, for example \(\mathrm{Beta}(50,50)\) (a strong expectation that \(p \approx 0{,}5\), say, from previous experiments with similar coins), the shift in the posterior estimate would be much stronger. The frequentist estimate \(\hat p = k/n\), on the other hand, does not depend at all on the choice of any prior distribution — in this approach, no such object even exists.

In this course we will mainly deal with parametric frequentist statistics. Although we will also touch on nonparametric and Bayesian methods.

Parametric statistics reduces to constructing and studying estimators of parameters, which is what we will be doing in the upcoming seminars.

2 Statistics: Definition, Examples

Definition 1 Suppose we have a sample of size \(n\). Any \((\mathscr {B}\left(\mathbb {R}^{n}\right) \mid \mathcal{E})\)-measurable map \(S: \mathbb {R}^{n} \rightarrow E\) is called a statistic, where \((E, \mathcal{E})\) is an arbitrary measurable space².

Suppose we are working with a parametric model, i.e., the family of distributions is parameterized as \(\mathcal{P}=\left\{ \mathbb {P}_{\theta }, \theta \in \Theta \right\}\). As a rule, the parameter is a number or a vector (as, for example, in the case of the normal distribution, where \(\theta = (\mu , \sigma^2)^T\)), i.e., \(\Theta \subset \mathbb {R}\) (or \(\mathbb {R}^{k}\)), and in this case it is assumed that a Borel sigma-algebra \(\mathscr {B}\left(\Theta \right)\) is defined on \(\Theta\). A measurable (\(\mathscr {B}\left(\mathbb {R}^{n}\right) \mid \mathscr {B}\left(\Theta \right)\)-measurable) map \(S: \mathbb {R}^{n} \to \Theta\) is called an estimator of the parameter \(\theta\).

As a rule, an estimator of the parameter \(\theta\) is denoted \(\widehat{\theta } = \widehat{\theta }(X)\). Note that an estimator of a parameter is a random variable, since it depends on the sample, i.e., on random variables. At the same time, the parameter itself is a nonrandom constant — it is simply unknown to us.

To summarize, if \(X=\left(X_{1}, \ldots , X_{n}\right)\) is a sample, then a statistic \(S(X)\) is a measurable function of the sample. An estimator of the parameter \(\theta \in \Theta\) is a statistic that takes values in \(\Theta\).

In the case of a parametric model, the problem of finding the “true” distribution reduces to finding the “true” value of the parameter. For this, a successful (in some sense) estimator of the parameter is found. In this chapter, we will discuss which estimators are considered “successful.”

Let \(X_{1}, \ldots , X_{n}\) be a sample from an unknown distribution \(\mathbb {P}\).

• Let \(g: \mathbb {R} \to \mathbb {R}\) be a Borel function. The quantity

  \[ \overline{g(X)} = \frac{1}{n} \sum _{i=1}^{n} g(X_i) \]

  is called the sample characteristic of \(g(x)\). For example,

  \[ \overline{X} := \frac{1}{n} \sum _{i=1}^{n} X_{i} \]

  is called the sample mean (here \(g\) is the identity function), and

  \[ \overline{X^k} := \frac{1}{n} \sum _{i=1}^{n} X_{i}^{k} \]

  is the sample \(k\)-th moment. Or, for example, one can speak of the sample logarithm:

  \[ \overline{\ln X} := \frac{1}{n} \sum _{i=1}^{n} \ln (X_{i}) \]

• Suppose the unknown distribution \(\mathbb {P}\) (\(\xi \sim \mathbb {P}\)) has a variance \(\sigma^2 := \operatorname {Var}\left[\xi \right]\), which we will treat as a parameter. How should we estimate it? Recall that \(\operatorname {Var}\left[\xi \right] = \mathbb {E}\left[\xi^2\right] - \left(\mathbb {E}\left[\xi \right]\right)^2\). The statistic

  \[ s^2 = \overline{X^2} - \left(\overline{X}\right)^2 \]

  is called the (biased) sample variance. It will become clear below why this statistic is called biased.

• Let us arrange the values of the sample \(\left(X_{1}, \ldots , X_{n}\right)\) in increasing order:

  \[ X_{(1)} \leq X_{(2)} \leq \ldots \leq X_{(n)} \]

  The resulting collection of statistics is called the order statistics, and its member \(X_{(k)}\) is the \(k\)-th order statistic. In particular, \(X_{(1)}\) is the minimum:

  \[ X_{(1)} = \min _{k \in \left\{ 1,\ldots ,n\right\} }X_k, \]

  and \(X_{(n)}\) is the sample maximum:

  \[ X_{(n)} = \max _{k \in \left\{ 1,\ldots ,n\right\} }X_k, \]

  The sample median is the statistic that splits the order statistics (the sample sorted in increasing order) in half:

  \[ X_{\mathrm{med}} := \begin{cases} X_{(m+1)}, & n = 2m + 1 \\ \frac{1}{2}\left(X_{(m)} + X_{(m+1)}\right), & n = 2m \end{cases} \]

3 Properties of Estimators

3.1 Non-Asymptotic Properties

Let \(X = (X_1, \ldots , X_n)\) be an observation from an unknown distribution \(\mathbb {P} \in \left\{ \mathbb {P}_{\theta }, \theta \in \Theta \right\}\).

Definition 2 An estimator \(\widehat{\theta } = \widehat{\theta }(X)\) is called an unbiased estimator of the parameter \(\theta\) if for every \(\theta \in \Theta\) the equality \[ \mathbb {E}_{\theta }\left[\widehat{\theta }\right] = \theta \] holds. That is, whatever the true parameter \(\theta \in \Theta\) may be, the expectation of the random variable \(\widehat{\theta }\) coincides with \(\theta\). Note here that the expectation, generally speaking, depends on \(\theta\), since the distribution \(\mathbb {P}_\theta\) depends on \(\theta\), and the expectation is a Lebesgue integral with respect to the unknown distribution. This is why \(\mathbb {E}\) carries the subscript \(\theta\).

ExampleExample 3

Let \(X_{1}, \ldots , X_{n}\) be a sample from an unknown distribution \(\mathbb {P} \in \left\{ \mathbb {P}_{\theta }, \theta \in \Theta \right\}\), \(\Theta \subset \mathbb {R}\), where the parameter \(\theta\) equals the expectation of the distribution: \(\mathbb {E}_{\theta }\left[X_1\right] = \theta\) for every \(\theta \in \Theta\). Prove that \(\overline{X}\) is an unbiased estimator of the parameter \(\theta\).

SolutionSolution

By linearity of expectation, \[ \mathbb {E}_{\theta }\left[\overline{X}\right] = \frac{1}{n}\sum _{i=1}^n \mathbb {E}_{\theta }\left[X_i\right] = \mathbb {E}_{\theta }\left[X_1\right] = \theta , \] for every \(\theta \in \Theta\).

ProblemProblem 1

Let \(X_{1}, \ldots , X_{n}\) be a sample from an unknown distribution \(\mathbb {P}_{\sigma^2}\) with parameter \(\sigma^{2}\), whose variance exists and equals \(\sigma^2\). Consider the (biased) sample variance \[ s^2 = \overline{X^2} - \left(\overline{X}\right)^2 \]

1. Prove that

\[ s^2 = \frac{1}{n} \sum _{i=1}^n \left(X_i - \overline{X}\right)^2. \]

2. Prove that

\[ \mathbb {E}_{\sigma ^2}\left[s^2\right] = \frac{n-1}{n} \sigma ^2, \quad \forall \sigma ^2 > 0, \]

i.e., \(s^2\) is indeed a biased estimator.

3. Prove that the estimator

\[ S^2 = \frac{n}{n-1}s^2 = \frac{1}{n-1}\sum _{i=1}^n \left(X_i - \overline{X}\right)^2, \]

called the unbiased sample variance, is indeed an unbiased estimator of \(\sigma^2\).

The following problem shows that unbiased estimators do not exist for all (complicated) parameters.

ProblemProblem 2

Let \(X_{1}, \ldots , X_{n}\) be a sample from the distribution \(\operatorname {Ber}(\theta )\), \(\theta \in \Theta = [0,1]\). For which functions \(\tau (\theta )\) do unbiased estimators exist?

3.2 Asymptotic Properties

In this section we assume that we have an infinitely large sample \(X_1, X_2, X_3, \ldots\) from an unknown distribution \(\mathbb {P}\). One can prove that there exists a probability space \(\left(\mathbb {R}^{\infty }, \mathscr {B}\left(\mathbb {R}^{\infty }\right), \mathbb {P}^{\infty }\right)\) on which one can define an infinite sequence of random variables \(X_{1}, X_{2}, \ldots\) (a sample of infinite size), distributed in the same way as \(\xi\)³. Thereby, a sample of any size \(n\) is defined on this space. Going forward, for simplicity of notation, we will write \((\mathbb {R}, \mathscr {B}\left(\mathbb {R}\right), \mathbb {P})\) instead of \(\left(\mathbb {R}^{\infty }, \mathscr {B}\left(\mathbb {R}^{\infty }\right), \mathbb {P}^{\infty }\right)\), respectively.

Suppose we have a sequence of estimators constructed according to the same principle, where the \(n\)-th estimator in the sequence depends only on the first \(n\) elements of the sample: \(\widehat{\theta }_1 = \widehat{\theta }_1(X_1)\), \(\widehat{\theta }_2 = \widehat{\theta }_2(X_1,X_2)\), \(\ldots\),

\[ \widehat{\theta }_n = \widehat{\theta }_n(X_1, \ldots , X_n), \quad \ldots , \]

For example, the sample mean: \(\widehat{\theta }_n = \frac{1}{n}\sum_{k=1}^{n}X_k\). Or the maximum: \(\widehat{\theta }_n = \max_{1 \leq k \leq n}X_k\). In such cases, it is common to omit the dependence on \(n\) and simply speak of the estimator \(\widehat{\theta }\).

Definition 3 Let \(\widehat{\theta }\) be an estimator (a sequence of estimators) of an unknown parameter \(\theta\).

• \(\widehat{\theta }\) is called a consistent estimator of \(\theta\) if, whatever the true parameter \(\theta\) may be, \(\widehat{\theta }\) converges to it in probability:

  \[ \widehat{\theta } \xrightarrow [n \to \infty ]{\mathbb {P}_{\theta }} \theta \]

  for every \(\theta \in \Theta\).

• \(\widehat{\theta }\) is called a strongly consistent estimator of \(\theta\) if \(\forall \theta \in \Theta\) we have

  \[ \widehat{\theta } \xrightarrow [n \to \infty ]{\mathbb {P}_{\theta }\text{-a.s.}} \theta \]

• \(\widehat{\theta }\) is called an asymptotically normal estimator of \(\theta\) if \(\forall \theta \in \Theta\) we have

  \[ \sqrt{n}\left(\widehat{\theta } - \theta \right) \xrightarrow [n \to \infty ]{d_\theta } \mathscr {N}\left(0, \sigma ^2(\theta )\right) \]

  for some function \(\sigma^{2}(\theta )\), which is called the asymptotic variance of the estimator \(\widehat{\theta }\).

  In the definition above we used the notation \(d_{\theta }\), which denotes convergence in distribution under the condition that the elements of the sample have distribution \(\mathbb {P}_{\theta }\).

All the definitions above carry over to the case where \(\theta\) is a multidimensional parameter. In this case, in the definition of asymptotic normality, the asymptotic variance \(\sigma^2(\theta )\) is replaced by the asymptotic covariance matrix \(\Sigma (\theta )\).

ExampleExample 4

Let a sample \(X_{1}, \ldots , X_{n}\) be given from an unknown distribution \(\mathbb {P} \in \mathcal{P}\). Suppose, moreover, that \(k \in \mathbb {N}\) and \[ 0 < \mathbb {E}_{\mathbb {P}}\left[\left|X_1\right|^{2k}\right] <\infty \] for all \(\mathbb {P} \in \mathcal{P}\), i.e., whatever the unknown distribution may be, it is known in advance that it has a finite \(2k\)-th moment. Prove that the sample \(k\)-th moment \(\overline{X^k}\) is an unbiased, strongly consistent, and asymptotically normal estimator of \(\mathbb {E}_{\mathbb {P}}\left[X_1^{k}\right]\).

SolutionSolution

• Unbiasedness, as in the previous problem, follows from the linearity of expectation.

• Strong consistency follows from the strong law of large numbers (Kolmogorov’s): the quantities \(X_1^k, X_2^k, \ldots\) are independent, identically distributed, and \(\mathbb {E}_{\mathbb {P}}\left[\left|X_1^k\right|\right] = \mathbb {E}_{\mathbb {P}}\left[\left|X_1\right|^{k}\right] < \infty\) (finiteness follows from the finiteness of the \(2k\)-th moment), so

  \[ \overline{X^k} = \frac{1}{n} \sum _{i=1}^{n} X_{i}^{k} \xrightarrow [n \to \infty ]{\text{a.s.}} \mathbb {E}_{\mathbb {P}}\left[X_{1}^{k}\right], \quad \forall \mathbb {P} \in \mathcal{P}. \]

  Almost sure convergence implies convergence in probability, so we obtain (ordinary) consistency along the way.

• Asymptotic normality follows from the central limit theorem: apply the CLT to \(X_1^k, X_2^k, \ldots\). Indeed, \(X_1, X_2, \ldots\) are independent and identically distributed, hence the sequence \(X_1^k, X_2^k, \ldots\) also consists of independent identically distributed random variables, with

  \[ \operatorname {Var}_{\mathbb {P}}\left[X_1^k\right] = \mathbb {E}_{\mathbb {P}}\left[X_1^{2k}\right] - \left(\mathbb {E}_{\mathbb {P}}\left[X_1^k\right]\right)^2 < \infty . \]

  Thus, we can apply the CLT. Applying it:

  \[ \sqrt{n}\left(\overline{X^{k}} - \mathbb {E}_{\mathbb {P}}\left[X_{1}^{k}\right]\right) \stackrel{d}{\rightarrow } \mathscr {N}\left(0, \operatorname {Var}_{\mathbb {P}}\left[X_1^k\right]\right). \]

The main tools for proving that an estimator is consistent or asymptotically normal are the laws of large numbers, the central limit theorem, the theorem on inheritance of convergence, and also the theorem on inheritance of asymptotic normality, formulated below.

ProblemProblem 3

Let \(X_1, \ldots , X_n\) be a sample from the distribution \(\operatorname {U}[0,\theta ]\), \(\theta \in \Theta = (0,+\infty )\). Check the estimator for unbiasedness, consistency, and strong consistency

1. \(2\overline{X}\),

2. \(\overline{X} + \frac{X_{(n)}}{2}\),

3. \((n+1)X_{(1)}\),

4. \(X_{(1)} + X_{(n)}\),

5. \(\frac{n+1}{n}X_{(n)}\).

HintHint

\(X_{(n)} = \max (X_1, \ldots , X_n)\), so \(X_{(n)} \leq t\) holds if and only if \(X_1, X_2\ldots ,X_n \leq t\). We have \[ \begin{align} \mathbb {P}\left(X_{(n)} \leq t\right) & = \mathbb {P}\left(X_1 \leq t, X_2 \leq t, \ldots , X_n \leq t\right) \stackrel{*}{=} \\ & \stackrel{*}{=}\mathbb {P}\left(X_1 \leq t\right) \cdot \ldots \cdot \mathbb {P}\left(X_n \leq t\right) \stackrel{**}{=} \\ & \stackrel{**}{=} \left(\mathbb {P}\left(X_1 \leq t\right)\right)^n, \end{align} \] where \(*\) holds by independence of the sample, and \(**\) by the sample elements \(X_1, \ldots , X_n\) being identically distributed.

ProblemProblem 4

Prove that asymptotic normality implies consistency. That is, let \(\widehat{\theta }(X) = \widehat{\theta }_{n}(X_1, \ldots , X_n)\) be an asymptotically normal estimator of the parameter \(\theta\) with asymptotic variance \(\sigma^{2}(\theta )\). Prove that then \(\widehat{\theta }(X)\) is a consistent estimator of \(\theta\).

HintHint

Apply Slutsky’s theorem from the next chapter.

4 Sample VS Realization of a Sample

Let us say a few words about what this distinction consists of. Imagine that you are conducting a real study — for example, collecting data on the intervals between customer visits. In real life you obtain concrete data, concrete numbers. Numbers, of course, are not random variables. So why does mathematical statistics work with random variables at all?

The point is that we are reasoning as if immediately before conducting the experiment. The experiment has already been planned, and we are about to go collect the data — but have not yet collected it. At this moment, for us, the data are unknown, i.e., they represent random variables. And it is precisely at this moment that we want to understand in advance what properties our future data will have: will, for example, the sample mean approximate the parameter of interest well (will unbiasedness hold, will the asymptotic properties hold), and so on.

In other words, mathematical statistics works with a sample — a collection of random variables

\[ X_1, X_2, X_3, \ldots , X_n, \]

whereas as a result of a real experiment, that is, for a fixed outcome \(\omega_0 \in \Omega\), these random variables turn into concrete numbers — into a realization of the sample

\[ X_1(\omega _0), X_2(\omega _0), X_3(\omega _0), \ldots , X_n(\omega _0). \]

5 Popular Families of Distributions

We will consider two main macro-families (i.e., families of families) of models: group families and exponential families. Together they cover most of the commonly used statistical models.

5.1 Group Families

Recall that a set with an operation \((G, \cdot )\) is called a group if

1. \(G\) is closed under the operation: \(g_1 \cdot g_2 \in G\) for all \(g_1, g_2 \in G\);

2. associativity holds: \((g_1 \cdot g_2) \cdot g_3 = g_1 \cdot (g_2 \cdot g_3)\);

3. an identity exists: there is an element in \(G\), denoted \(e\) or \(\operatorname {id}\), such that \(g \cdot e = e \cdot g = g\) for every \(g \in G\);

4. every element has an inverse: for every element \(g \in G\) there exists an element in \(G\), denoted \(g^{-1}\), such that \(g \cdot g^{-1} = g^{-1} \cdot g = e\).

In statistics, groups of transformations (with the composition operation) are of primary interest. Recall that a transformation is a map of a set into itself. Examples:

• shift transformations on \(\mathbb {R}\):

  \[ \operatorname {TL} := \left\{ x \mapsto a + x \; : \; a \in \mathbb {R}\right\} ; \]

  in the case of \(\mathbb {R}^{m}\):

  \[ \operatorname {TL}_{m} := \left\{ \vec{x} \mapsto \vec{a} + \vec{x} \; : \; \vec{a} \in \mathbb {R}^{m}\right\} \]

• linear maps on \(\mathbb {R}^{m}\):

  \[ \operatorname {GL}_m := \left\{ f: \mathbb {R}^{m} \to \mathbb {R}^{m} \; : \; f \text{ is linear and bijective}\right\} \]

  \(\operatorname {GL}_m\) is called the general linear group. It is, as is well known, isomorphic to the group of nonsingular matrices of order \(m\) with the operation of multiplication. Namely,

  \[ \operatorname {GL}_m = \left\{ \vec{x} \mapsto \mathbf{b}\vec{x}\; : \; \mathbf{b} \in \mathbb {R}^{m \times m} \text{ -- a nonsingular matrix}\right\} \]

  \(\operatorname {GL}_m\) is often identified with the group of matrices.

• A subgroup of the previous group:

  \[ \operatorname {GL}^+_m := \left\{ \vec{x} \mapsto \mathbf{b}\vec{x}\; : \; \mathbf{b} \in \mathbb {R}^{m \times m} , \; \operatorname {det}\mathbf{b} > 0\right\} \]

• All possible compositions of linear transformations and shifts (parallel translations). Such transformations are called affine:

  \[ \operatorname {Aff}_m := \left\{ \vec{x} \mapsto \vec{a} + \mathbf{b}\vec{x} \; : \; \vec{a} \in \mathbb {R}^{m}, \; \mathbf{b} \in \mathbb {R}^{m \times m}, \; \operatorname {det}\mathbf{b} \neq 0\right\} \]

  We will also consider its subgroup

  \[ \operatorname {Aff}^+_m := \left\{ \vec{x} \mapsto \vec{a} + \mathbf{b}\vec{x} \; : \; \vec{a} \in \mathbb {R}^{m}, \; \mathbf{b} \in \mathbb {R}^{m \times m}, \; \operatorname {det}\mathbf{b} > 0\right\} \]

• All possible continuous bijections (i.e., homeomorphisms) on \(\mathbb {R}\). Notation: \(\operatorname {Homeo}(\mathbb {R})\). One can prove that

  \[ \operatorname {Homeo}(\mathbb {R}) = \left\{ g \in C(\mathbb {R}) \; : \; \left[\begin{array}{ll} g \text{ strictly increasing, } & \lim \limits _{x \to \pm \infty }g(x) = \pm \infty \\ g \text{ strictly decreasing, } & \lim \limits _{x \to \pm \infty }g(x) = \mp \infty \\ \end{array}\right.\right\} \]

• A subgroup of the previous group: all possible orientation-preserving homeomorphisms:

  \[ \operatorname {Homeo}_+(\mathbb {R}) = \left\{ g \in C(\mathbb {R}) \; :\; g \text{ strictly increasing, } \lim \limits _{x \to \pm \infty }g(x) = \pm \infty \right\} \]

• A subgroup of the previous group: all possible homeomorphisms preserving orientation and symmetry about zero:

  \[ \operatorname {Homeo}^{\operatorname {odd}}_+(\mathbb {R}) = \left\{ g \in C(\mathbb {R}) \; :\; g \text{ is odd, strictly increasing, } \lim \limits _{x \to \pm \infty }g(x) = \pm \infty \right\} \]

Why are these groups? Associativity for an arbitrary collection of transformations always holds, since it holds for composition. One needs to check the 3 remaining properties: closure under composition, the presence of the identity transformation in the collection, the presence of an inverse transformation. The last point splits into 2: for each transformation one needs to check bijectivity and that the inverse map also lies in the collection.

Let us prove the group properties using one example.

ExampleExample 5

Show that \(\operatorname {Aff}^+_{m}\) is indeed a group.

SolutionSolution

• Closure under composition? Let \(\vec{a}_1, \vec{a}_2 \in \mathbb {R}^{m}\), let \(\mathbf{b}_1, \mathbf{b}_2 \in \mathbb {R}^{m \times m}\), \(\operatorname {det}\left(\mathbf{b}_1\right), \operatorname {det}\left(\mathbf{b}_2\right) > 0\). Let \(f_{i}(\vec{x}) = \vec{a}_i + \mathbf{b}_i\vec{x}\), \(i = 1, 2\). Then

  \[ \begin{align} (f_1 \circ f_2)(\vec{x}) & = f_1 \left(f_2(\vec{x})\right) = \vec{a}_1 + \mathbf{b}_1 \cdot \left(\vec{a}_2 + \mathbf{b}_2 \cdot \vec{x}\right) = \\ & = \vec{a}_1 + \mathbf{b}_1\vec{a}_2 + \mathbf{b}_1\mathbf{b}_2 \cdot \vec{x} \end{align} \]

  This is again a transformation of the form \(\vec{x} \mapsto \vec{a} + \mathbf{b}\vec{x}\) with \(\vec{a} = \vec{a}_1 + \mathbf{b}_1\vec{a}_2\), \(\mathbf{b} = \mathbf{b}_1\mathbf{b}_2\), where \(\operatorname {det}\left(\mathbf{b}_1\mathbf{b}_2\right) = \operatorname {det}\left(\mathbf{b}_1\right) \cdot \operatorname {det}\left(\mathbf{b}_2\right) > 0\), i.e., \(f_1 \circ f_2 \in \operatorname {Aff}^+_m\).

• The presence of an identity transformation? It suffices to take \(\vec{a} = \vec{0}\), \(\mathbf{b} = \mathbf{e}\), where \(\mathbf{e}\) is the \(m\)-dimensional identity matrix. Indeed, \(\operatorname {det}\mathbf{e} = 1 > 0\).

• The presence of an inverse transformation? Let \(f(\vec{x}) = \vec{a} + \mathbf{b}\vec{x}\), \(\operatorname {det}\left(\mathbf{b}\right) > 0\). Then \(\mathbf{b}\) is invertible, and the equation \(\vec{y} = \vec{a} + \mathbf{b}\vec{x}\) is uniquely solvable for \(\vec{x}\):

  \[ \vec{x} = \mathbf{b}^{-1}\left(\vec{y} - \vec{a}\right) = -\mathbf{b}^{-1}\vec{a} + \mathbf{b}^{-1}\vec{y}, \]

  i.e., \(f\) is a bijection and \(f^{-1}(\vec{y}) = -\mathbf{b}^{-1}\vec{a} + \mathbf{b}^{-1}\vec{y}\) is again a transformation of the same form, where \(\operatorname {det}\left(\mathbf{b}^{-1}\right) = \frac{1}{\operatorname {det}\left(\mathbf{b}\right)} > 0\). Hence, \(f^{-1} \in \operatorname {Aff}^+_m\).

Definition 4 Let \(\mu\) be a measure on a measurable space \((\mathcal{X}, \mathcal{F})\), and let \(g: \mathcal{X} \to \mathcal{X}\) be some measurable transformation of \(\mathcal{X}\). Then the pushforward of the measure \(\mu\) is the measure \(\mu \circ g^{-1}\), also denoted \(g_{\sharp }\mu\). In other words,

\[ g_{\sharp }\mu (B) := \left(\mu \circ g^{-1}\right)(B) = \mu \left(g^{-1}(B)\right), \quad B \in \mathcal{F} \]

Definition 5 Let \(G\) be some group of transformations of a set \(\mathcal{X}\). A family \(\mathcal{P}\) of distributions on \(\mathcal{X}\) is called a group family with transformation group \(G\) on \(\mathcal{X}\) if \(\mathcal{P}\) is the orbit of some distribution \(\mathbb {P}_0\) under the action of \(G\): \[ \mathcal{P} = \left\{ g_{\sharp }\mathbb {P}_0 \; : \; g \in G\right\} . \] In particular, \(\mathcal{P}\) is closed under transformations from \(G\) (\(g_{\sharp }\mathbb {P} \in \mathcal{P}\) for all \(g \in G\), \(\mathbb {P} \in \mathcal{P}\)), and any two distributions from \(\mathcal{P}\) are carried into each other by some transformation from \(G\): if \(\mathbb {P}_1 = g_{1\sharp }\mathbb {P}_0\) and \(\mathbb {P}_2 = g_{2\sharp }\mathbb {P}_0\), then \(\mathbb {P}_2 = \left(g_2 g_1^{-1}\right)_{\sharp }\mathbb {P}_1\).

5.1.1 One-Dimensional Location-Scale Families

Let \(U\) be a random variable with a fixed distribution \(F\), having (if it exists) a density \(f\) with respect to Lebesgue measure. In addition to the shift group \(\operatorname {TL}\) already introduced, we will need the group of scale transformations — this is precisely the group \(\operatorname {GL}^+_1\) already familiar to us (the one-dimensional special case of \(\operatorname {GL}^+_m\)):

\[ \operatorname {GL}^+_1 = \left\{ x \mapsto bx \; : \; b > 0\right\} . \]

The group generated by \(\operatorname {TL}\) and \(\operatorname {GL}^+_1\) together is precisely the group already familiar to us

\[ \operatorname {LS} := \left\{ x \mapsto a + bx \; : \; a \in \mathbb {R}, \, b > 0\right\} = \operatorname {Aff}^+_1, \]

so the fact that \(\operatorname {LS}\) and \(\operatorname {GL}^+_1\) are indeed groups has, in essence, already been proved above (these are one-dimensional special cases of the groups \(\operatorname {Aff}^+_m\) and \(\operatorname {GL}^+_m\)).

Definition 6 Let \(U \sim F\).

• A group family with transformation group \(\operatorname {TL}\) (i.e., \(\mathcal{P} = \left\{ g_{\sharp }\mathbb {P}_U \; : \; g \in \operatorname {TL}\right\}\)) is called a location family. If \(X = U + a\), then

  \[ \mathbb {P}\left(X \leq x\right) = F(x - a). \]

• A group family with group \(\operatorname {GL}^+_1\) is called a scale family. If \(X = bU\), \(b>0\), then

  \[ \mathbb {P}\left(X \leq x\right) = F(x/b). \]

• A group family with group \(\operatorname {LS} = \operatorname {Aff}^+_1\) is called a location-scale family. If \(X = a + bU\), then

  \[ \mathbb {P}\left(X \leq x\right) = F\left(\frac{x-a}{b}\right). \]

If \(F\) has a density \(f\), then the density of \(X = a+bU\) equals

\[ \frac{1}{b}f\left(\frac{x-a}{b}\right). \]

Below are listed the most commonly used location-scale families on \(\mathbb {R}\), which will come up repeatedly later on.

Name	Density	Support	Notation
Normal	\(\dfrac{1}{\sqrt{2\pi}\,b}e^{-(x-a)^2/2b^2}\)	\(\mathbb{R}\)	\(\mathcal{N}(a,b^2)\)
Double exponential (Laplace)	\(\dfrac{1}{2b}e^{-|x-a|/b}\)	\(\mathbb{R}\)	\(DE(a,b)\)
Cauchy	\(\dfrac{b}{\pi}\cdot\dfrac{1}{b^2+(x-a)^2}\)	\(\mathbb{R}\)	\(C(a,b)\)
Logistic	\(\dfrac{1}{b}\cdot\dfrac{e^{-(x-a)/b}}{\left(1+e^{-(x-a)/b}\right)^2}\)	\(\mathbb{R}\)	\(L(a,b)\)
Exponential	\(\dfrac{1}{b}e^{-(x-a)/b}\,\mathbb{1}_{[a,+\infty)}(x)\)	\((a, +\infty)\)	\(E(a,b)\)
Uniform	\(\dfrac{1}{b}\,\mathbb{1}_{[a-b/2,\,a+b/2]}(x)\)	\(\left(a-\dfrac{b}{2},\ a+\dfrac{b}{2}\right)\)	\(U\!\left(a-\dfrac{b}{2},\ a+\dfrac{b}{2}\right)\)

(here \(-\infty < a < \infty\) is the location parameter, \(b > 0\) is the scale parameter; the connection to the standard notation for distributions from probability theory is given in the last two columns of the table).

5.1.2 Multidimensional Location-Scale Families

All the families from the previous item generalize naturally to the case where \(U\) is a random vector \(\vec{U} = (U_1, \ldots , U_n)\). If \(U_1,\ldots ,U_n\) are i.i.d. (for example, according to one of the distributions in the table above), one can consider the coordinatewise action of the group \(\operatorname {LS}\): \(\vec X = a \vec{1} + b\vec U\) (i.e., the same element \(g\in \operatorname {LS}\) is applied simultaneously to all coordinates). A broader family results if we allow different shift and scale for each coordinate:

\[ X_i = a_i + b_i U_i, \quad i = 1,\ldots ,n, \]

i.e., if one uses the group

\[ \operatorname {LS}_n := \left\{ \vec x \mapsto \vec a + \mathbf{b}\vec x \; : \; \vec a \in \mathbb {R}^{n}, \; \mathbf{b} = \operatorname {diag}(b_1,\ldots ,b_n), \; b_i > 0\right\} \subset \operatorname {Aff}_n. \]

Families generated not by a coordinatewise action but by a full linear transformation deserve special attention — this leads to the following important example.

ExampleExample 6

The multivariate normal distribution.

SolutionSolution

Let \(\vec U \sim \mathscr {N}\left(\vec0, \mathbf{I}_p\right)\) (i.e., the coordinates \(U_i\) are independent and \(U_i \sim \mathscr {N}\left(0, 1\right)\)). Let us act on \(\vec U\) with the group \(\operatorname {Aff}_p\) already familiar to us:

\[ \vec X = \vec a + \mathbf{b} \vec U, \qquad g(\vec x) = \vec a + \mathbf{b} \vec x \in \operatorname {Aff}_p. \]

The resulting group family is precisely the family of all nondegenerate \(p\)-dimensional normal distributions \(\mathscr {N}\left(\vec a, \Sigma \right)\), \(\Sigma = \mathbf{b} \mathbf{b}^T\): the vector \(\vec a\) ranges over all of \(\mathbb {R}^{p}\), while \(\mathbf{b} \mathbf{b}^T\), as \(\mathbf{b}\) ranges over all nonsingular matrices, ranges precisely over all positive definite matrices.

This group family is structured a little differently from the previous ones: different elements of the group \(\operatorname {Aff}_p\) can generate the same distribution. Namely, \(\vec a_1 + \mathbf{b}_1 \vec U\) and \(\vec a_2 + \mathbf{b}_2 \vec U\) are identically distributed if and only if \(\vec a_1 = \vec a_2\) and \(\mathbf{b}_1 \mathbf{b}_1^T = \mathbf{b}_2 \mathbf{b}_2^T\), i.e., when \(\mathbf{b}_2^{-1}\mathbf{b}_1\) is an orthogonal matrix. Therefore, the same family of distributions can be obtained by restricting to a much smaller subgroup \(G \subset \operatorname {Aff}_p\) — for example, by taking only lower triangular matrices for \(\mathbf{b}\) (this is related to the Cholesky decomposition of the matrix \(\Sigma\)).

One more important example is the linear model (in particular, linear regression, which we will discuss in detail later). Let \(\vec U = (U_1,\ldots ,U_n)\) have a fixed distribution (most often, \(U_i\) i.i.d. \(\sim \mathscr {N}\left(0, 1\right)\)), and let

\[ X_i = a_i + bU_i, \quad i=1,\ldots ,n, \]

where the shift vector \(\vec a = (a_1,\ldots ,a_n)\) is required to lie in some fixed \(s\)-dimensional subspace \(\Omega \subset \mathbb {R}^{n}\) (i.e., \(a_i = \sum_{j=1}^s d_{ij}\beta_j\) for some known matrix \(D=(d_{ij})\) of rank \(s\) and arbitrary \(\beta_j \in \mathbb {R}\)), and the scale \(b>0\) is common to all coordinates. The corresponding group of transformations

\[ \left\{ \vec x \mapsto \vec a + b\vec x \; : \; \vec a \in \Omega , \; b > 0\right\} \]

is a subgroup of \(\operatorname {Aff}_n\). In the normal case, the joint density of \(\vec X\) has the form

\[ \frac{1}{(\sqrt{2\pi }b)^n}\exp \left[-\frac{1}{2b^2}\sum _{i=1}^n (x_i - a_i)^2\right], \qquad \vec a \in \Omega . \]

5.1.3 (*) Nonparametric Group Families

Unlike the previous two items, where the families were effectively parameterized by a small number of numbers (\(a\) and \(b\); the elements \(\vec a\) and \(\mathbf{b}\)), groups of transformations can also generate nonparametric families that are much “poorer” in structure.

ExampleExample 7

The nonparametric i.i.d. family.

SolutionSolution

Let \(U_1,\ldots ,U_n\) be independent and identically distributed according to some fixed continuous distribution whose support is the entire real line (for example, \(\mathscr {N}\left(0, 1\right)\)). Consider the already familiar group \(\operatorname {Homeo}_+(\mathbb {R})\) of strictly increasing homeomorphisms of \(\mathbb {R}\) onto itself, acting identically on each coordinate: \[ X_i = g(U_i), \quad i = 1,\ldots ,n, \qquad g \in \operatorname {Homeo}_+(\mathbb {R}). \] The quantities \(X_1,\ldots ,X_n\) are still i.i.d., with some common distribution function \(F_g\). The corresponding group family \(\left\{ F_g \; : \; g \in \operatorname {Homeo}_+(\mathbb {R})\right\}\) coincides with the class of all continuous distributions with support \(\mathbb {R}\) (i.e., with a continuous strictly increasing distribution function on \(\mathbb {R}\)) — there is no finite-dimensional parameter here anymore.

ExampleExample 8

Symmetric distributions.

SolutionSolution

Let us restrict the group from the previous example to the subgroup of odd functions — this is precisely the group \(\operatorname {Homeo}_+^{\operatorname {odd}}(\mathbb {R})\) already defined above. Then \(X_i = g(U_i)\), \(g\in \operatorname {Homeo}_+^{\operatorname {odd}}(\mathbb {R})\), generates the class of all distributions with support \(\mathbb {R}\) that are symmetric about zero. If we additionally allow a shift, i.e., consider \[ X_i = a + g(U_i), \quad a \in \mathbb {R}, \; g \in \operatorname {Homeo}_+^{\operatorname {odd}}(\mathbb {R}), \] we obtain the class of all distributions with support \(\mathbb {R}\) symmetric about an (unknown) point \(a\). Such families are often used in problems of estimating a center of symmetry, when the shape of the distribution itself is not assumed to be known.

Similarly, one can construct other nonparametric group families. For example, if we take \(U_i\) i.i.d. \(\operatorname {U}[0,1]\), and take as \(g\) all increasing continuous bijections of \((0,1)\) onto itself with \(g(0)=0\), \(g(1)=1\) (this too is a group of transformations), and then add a location-scale shift \(X_i = a + bg(U_i)\), we obtain the class of all continuous distributions whose support is some interval.

A separate, third type of group family arises in survey sampling. Suppose there is a finite population of \(N\) elements with values \(v_1,\ldots ,v_N\), and a random sample of size \(n\) without replacement gives pairs “value-index” \((U_1,J_1),\ldots ,(U_n,J_n)\) (all \(\binom {N}{n}\) possible selections being equally likely). The group of shifts

\[ X_i = U_i + a_{J_i}, \quad i=1,\ldots ,n, \qquad (a_1,\ldots ,a_N) \in \mathbb {R}^{N}, \]

turns this model into a group family: setting \(y_i = v_i + a_i\), the pairs \((X_1,J_1),\ldots ,(X_n,J_n)\) are simply a random sample without replacement from the population \((y_1,1),\ldots ,(y_N,N)\) with arbitrary (unknown) values \(y_i\).

5.2 Exponential Families

Definition 7 Let there be a dominated collection of distributions \(\mathcal{P} = \left\{ \mathbb {P}_{\theta }, \theta \in \Theta \subset \mathbb {R}^{k}\right\}\). Let \(p_\theta\) be a generalized density of the distribution \(\mathbb {P}_{\theta }, \; \theta \in \Theta\). That is, there is some fixed measure \(\mu\) on \(\mathbb {R}\) such that \(\mathbb {P}_\theta (dx) = p_\theta (x) \mu (dx)\). As a rule, \(\mu\) is Lebesgue measure \(\lambda\) (in which case \(\mathcal{P}\) is a collection of continuous distributions) or the counting measure \(\#\) on \(\mathbb {Z}\) (in which case \(\mathcal{P}\) is a collection of discrete distributions).

The collection \(\mathcal{P}\) is said to belong to an exponential family if the “generalized” density of the distribution \(\mathbb {P}_{\theta }\) has the form

\[ p_{\theta }(x)=h(x) \exp \left(\sum _{i=1}^{k} a_{i}(\theta ) u_{i}(x)+v(\theta )\right) . \tag{1}\]

Note that the number of terms of the form \(a_{i}(\theta ) u_{i}(x)\) in the exponent must not exceed the dimension of the parameter set. Moreover, the support of the distribution cannot depend on the parameter.

Note that the normal, exponential, binomial, Poisson, and many other distributions (collections of parameterized distributions) belong to an exponential family. Note that the exponential distribution is merely a special case of the exponential family.

ExampleExample 9

Prove that the normal distribution \(\mathscr {N}\left(\mu , \sigma^2\right), \; \theta = (\mu , \sigma^2)^T \in \Theta = \mathbb {R} \times \mathbb {R}_+ \subset \mathbb {R}^{2}\) is an exponential family.

SolutionSolution

Let us move the factor \(\frac{1}{\sqrt{2\pi \sigma^2}}\) under the exponent (via the logarithm) and expand the square in the exponent: \[ \begin{align} p_{\theta }(x) & = \frac{1}{\sqrt{2 \pi \sigma ^2}}e^{-\frac{(x-\mu )^2}{2\sigma ^2}} = \exp \left(-\frac{(x-\mu )^2}{2\sigma ^2} - \frac{1}{2}\ln (2\pi \sigma ^2)\right) = \\ & = \exp \left(-\frac{x^2 - 2\mu x + \mu ^2}{2\sigma ^2} - \frac{1}{2}\ln (2\pi \sigma ^2)\right) = \\ & = \exp \left(-\frac{1}{2\sigma ^2}x^2 + \frac{\mu }{\sigma ^2}x - \frac{\mu ^2}{2\sigma ^2} - \frac{1}{2}\ln (2\pi \sigma ^2)\right) \end{align} \] Hence \[ \begin{align} \alpha _1(\theta ) & = -\frac{1}{2\sigma ^2} , \; u_1(x) = x^2 \\ \alpha _2(\theta ) & = \frac{\mu }{\sigma ^2} , \; u_2(x) = x \\ v(\theta ) & = - \frac{\mu ^2}{2\sigma ^2} - \frac{1}{2}\ln (2\pi \sigma ^2) \end{align} \]

Name	Density in standard form	Standard parameterization	Standard parameter set	Density in the form from the definition above	Canonical parameterization	Canonical parameter set	\(k\)	Support	Notation
Gamma	\(\dfrac{1}{\Gamma(a)b^a}x^{a-1}e^{-x/b}\)	\((a,b)\)	\((0,+\infty)^2\)	\(\mathbb{1}_{(0,+\infty)}(x)\exp\left[(a-1)\ln x - \dfrac{x}{b} - \ln\Gamma(a) - a\ln b\right]\)	\(\alpha_1 = a-1,\ \alpha_2 = -\dfrac{1}{b}\)	\((-1,+\infty)\times(-\infty,0)\)	2	\((0,+\infty)\)	\(\Gamma(a,b)\)
Chi-squared	\(\dfrac{1}{\Gamma(f/2)2^{f/2}}x^{f/2-1}e^{-x/2}\)	\((f)\)	\((0,+\infty)\)	\(e^{-x/2}\mathbb{1}_{(0,+\infty)}(x)\exp\left[\left(\dfrac{f}{2}-1\right)\ln x - \ln\Gamma\!\left(\dfrac{f}{2}\right) - \dfrac{f}{2}\ln 2\right]\)	\(\alpha_1 = \dfrac{f}{2} - 1\)	\((-1,+\infty)\)	1	\((0,+\infty)\)	\(\chi^2_f\)
Beta	\(\dfrac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}x^{a-1}(1-x)^{b-1}\)	\((a,b)\)	\((0,+\infty)^2\)	\(\mathbb{1}_{(0,1)}(x)\exp\left[(a-1)\ln x + (b-1)\ln(1-x) + \ln\Gamma(a+b) - \ln\Gamma(a) - \ln\Gamma(b)\right]\)	\(\alpha_1 = a-1,\ \alpha_2 = b-1\)	\((-1,+\infty)^2\)	2	\((0,1)\)	\(B(a,b)\)
Bernoulli	\(p^x(1-p)^{1-x}\)	\((p)\)	\((0,1)\)	\(\exp\left[x\ln\dfrac{p}{1-p} + \ln(1-p)\right]\)	\(\alpha_1 = \ln\dfrac{p}{1-p}\)	\(\mathbb{R}\)	1	\(\{0,1\}\)	\(\mathrm{Ber}(p)\)
Binomial	\(\dbinom{n}{x}p^x(1-p)^{n-x}\)	\((p)\); \(n\) known	\((0,1)\)	\(\dbinom{n}{x}\exp\left[x\ln\dfrac{p}{1-p} + n\ln(1-p)\right]\)	\(\alpha_1 = \ln\dfrac{p}{1-p}\)	\(\mathbb{R}\)	1	\(\{0,1,\ldots,n\}\)	\(\mathrm{Bin}(p, n)\)
Poisson	\(\dfrac{\lambda^x}{x!}e^{-\lambda}\)	\((\lambda)\)	\((0,+\infty)\)	\(\dfrac{1}{x!}\exp\left[x\ln\lambda - \lambda\right]\)	\(\alpha_1 = \ln\lambda\)	\(\mathbb{R}\)	1	\(\{0,1,\ldots\}\)	\(\mathrm{Pois}(\lambda)\)
Negative binomial	\(\dbinom{m+x-1}{m-1}p^m q^x\)	\((p)\); \(m\) known	\((0,1)\)	\(\dbinom{m+x-1}{m-1}\exp\left[x\ln(1-p) + m\ln p\right]\)	\(\alpha_1 = \ln(1-p)\)	\((-\infty,0)\)	1	\(\{0,1,\ldots\}\)	\(Nb(p,m)\)

(the first three rows are densities with respect to Lebesgue measure, the last four with respect to counting measure; \(q := 1-p\). The column \(k\) is the dimension of the family, i.e., the number of terms \(\alpha_i(\theta)u_i(x)\) in the definition above; the "canonical parameterization" is precisely the quantities \(\alpha_i(\theta)\) from that definition, and the "canonical parameter set" is the set of values that the vector \((\alpha_1(\theta),\ldots,\alpha_k(\theta))\) traces out as \(\theta\) ranges over the standard parameter set.)

It is often convenient, instead of \(\theta\), to parameterize the family directly by the quantities \(\eta_i := a_i(\theta )\), \(i=1,\ldots ,k\) — these are called the natural (canonical) parameters. In this notation, the density Equation 1 takes the canonical form

\[ p(x \mid \eta ) = h(x)\exp \left(\sum _{i=1}^k \eta _i u_i(x) + v(\eta )\right), \qquad \eta = (\eta _1,\ldots ,\eta _k) \in \mathbb {R}^{k}. \]

The set

\[ \Xi := \left\{ \eta \in \mathbb {R}^{k} \; : \; \int e^{\sum _{i=1}^k \eta _i u_i(x)} h(x)\, \mu (dx) < \infty \right\} \]

is called the natural parameter set: this is precisely the set of those \(\eta\) for which one can choose a normalizing term \(v(\eta )\) so that \(p(x\mid \eta )\) is a density. One can show that \(\Xi\) is always convex.

Note that the representation Equation 1 is not unique: for example, one can multiply \(u_i\) by a constant \(c \neq 0\) while simultaneously dividing \(\eta_i\) by \(c\); more generally, one can perform simultaneous (identical) linear transformations of the collections \((u_1,\ldots ,u_k)\) and \((\eta_1,\ldots ,\eta_k)\).

ExampleExample 10

In \(n\) independent trials with \(s+1\) possible outcomes, the probability of the \(i\)-th outcome equals \(p_i\), \(i=0,1,\ldots ,s\) (\(\sum_{i=0}^s p_i = 1\)). Let \(X_i\) be the number of trials ending in outcome \(i\); the vector \((X_0,\ldots ,X_s)\) has the multinomial distribution \(M(p_0,\ldots ,p_s;n)\): \[ \mathbb {P}\left(X_0=x_0,\ldots ,X_s=x_s\right) = \frac{n!}{x_0!\cdots x_s!}\, p_0^{x_0}\cdots p_s^{x_s}. \] Prove that this is an \(s\)-dimensional (not \((s+1)\)-dimensional!) exponential family.

SolutionSolution

Let us rewrite the probability as \[ \frac{n!}{x_0!\cdots x_s!}\exp \left(x_0\ln p_0 + \ldots + x_s \ln p_s\right). \] Since the trials add up to \(n\), i.e., \(x_0 = n - x_1 - \ldots - x_s\), the exponent can be rewritten so that only \(x_1,\ldots ,x_s\) appear in it: \[ x_0 \ln p_0 + \sum _{i=1}^s x_i \ln p_i = n\ln p_0 + \sum _{i=1}^s x_i \ln \frac{p_i}{p_0}. \] Hence, \[ \mathbb {P}\left(X_0=x_0,\ldots ,X_s=x_s\right) = \frac{n!}{x_0!\cdots x_s!}\exp \left(\sum _{i=1}^s \eta _i x_i + v(\eta )\right), \] where \(\eta_i = \ln (p_i/p_0)\), \(u_i(x)=x_i\), \(i=1,\ldots ,s\), and \(v(\eta ) = n\ln p_0 = -n\ln \left(1+\sum_{i=1}^s e^{\eta_i}\right)\) (the last equality is a consequence of the condition \(\sum_i p_i=1\)). This is indeed an \(s\)-dimensional family: even though there are \(s+1\) outcomes, the linear constraint \(\sum_i x_i = n\) on the statistics “cuts off” one dimension.

Not every representation of the form Equation 1 is economical — this is exactly why the definition includes the caveat that the number of terms “must not exceed the dimension of the parameter set.”

Definition 8  

• Let \(X \sim \mathbb {P}_\theta\). The parameter \(\theta\) is called unidentifiable from \(X\) if there exist \(\theta_1 \neq \theta_2\) such that \(\mathbb {P}_{\theta_1} = \mathbb {P}_{\theta_2}\).

• If the statistics \(u_1,\ldots ,u_k\) or the parameters \(\eta_1(\theta ),\ldots ,\eta_k(\theta )\) are connected by a linear relation, the number of terms in Equation 1 can be reduced; otherwise the representation is called minimal, and the family itself full rank.

It is precisely the linear independence of \(\eta_1(\theta ), \ldots , \eta_k(\theta )\) (i.e., full rank) that is the condition required in the theorem below. If, on the other hand, the parameters \(\eta_i(\theta )\), while formally forming a \(k\)-dimensional family, actually lie, as functions of \(\theta\), on a submanifold of smaller dimension, the representation is not minimal; such a family is called a curved exponential family.

ExampleExample 11

Show that the normal family \(\mathscr {N}\left(\xi , \xi^2\right)\), \(\xi > 0\) (i.e., \(\sigma = \xi\) in the example above), is formally given by two canonical parameters, but is in fact a curved one-dimensional exponential family.

SolutionSolution

Substituting \(\mu =\xi , \sigma^2=\xi^2\) into the example above, we get \[ p_\xi (x) = \exp \left(\frac1\xi x - \frac{1}{2\xi ^2}x^2 - \frac12\right)\cdot \frac{1}{\sqrt{2\pi }\, \xi }. \] Formally this is still an expression of the form Equation 1 with \(k=2\), \(\eta_1(\xi )=1/\xi\), \(\eta_2(\xi )=-1/(2\xi^2)\). However, \(\theta =\xi\) here is a single parameter, and the point \((\eta_1,\eta_2)\), as \(\xi\) varies, traces out not an open subset of the plane but a one-dimensional curve \(\eta_2 = -\eta_1^2/2\). The representation is not minimal — the family is curved, and, in particular, the theorem below does not apply to it.

Theorem 1 (on exponential families) Let \(\theta \in \Theta \subset \mathbb {R}^{k}\), and suppose condition Equation 1 holds for the family of distributions \(\left\{ \mathbb {P}_{\theta }, \theta \in \Theta \right\}\). Suppose, moreover, that the set of all values of the vector \(\left(a_{1}(\theta ), \ldots , a_{k}(\theta )\right)\) for \(\theta \in \Theta\) contains a \(k\)-dimensional ball of the space \(\mathbb {R}^{k}\). Then \(S(X)=\left(\overline{u_{1}(X)}, \ldots , \overline{u_{k}(X)}\right)\) is a complete sufficient statistic.

For the normal distribution, for example, the statistic \((\overline{X}, \overline{X^2})\) will be complete and sufficient. Note also that any bijective function preserves completeness and sufficiency. That is, for example, for the normal distribution the statistic \((\overline{X}, s^2) = (\overline{X}, \overline{X^2} - \overline{X}^2)\) will also be complete and sufficient.

It is worth noting that the condition of the theorem holds if the set \(\Theta\) is “solid,” i.e., contains its interior points, and the functions \(a_{1}(\theta ), \ldots , a_{k}(\theta )\) are linearly independent and smooth in a neighborhood of some interior point \(\theta_{0} \in \Theta\).

5.2.1 Closure Under Sampling

The exponential structure is stable under combining independent observations. Namely, if \(X\) and \(Y\) are independent and their densities are written in canonical form with the same set of natural parameters \(\eta\):

\[ p(x\mid \eta ) = h(x)\exp \left(\sum _{i=1}^k \eta _i T_i(x) + v(\eta )\right), \qquad q(y\mid \eta ) = g(y)\exp \left(\sum _{i=1}^k \eta _i U_i(y) + w(\eta )\right), \]

then the joint distribution of the pair \((X,Y)\) is again a \(k\)-dimensional exponential family with the same \(\eta\), and the statistics add: \(T_i(x) + U_i(y)\). By induction this carries over to any (finite) number of independent terms.

In particular, if \(X_1,\ldots ,X_n\) is a sample (i.i.d.) from the family Equation 1, then the joint density of the sample equals

\[ \exp \left(\sum _{i=1}^k a_i(\theta )\sum _{j=1}^n u_i(X_j) + nv(\theta )\right)\prod _{j=1}^n h(X_j) = \exp \left(n\sum _{i=1}^k a_i(\theta )\overline{u_i(X)} + nv(\theta )\right)\prod _{j=1}^n h(X_j), \]

i.e., the sample also forms a \(k\)-dimensional exponential family — with the same number of parameters \(k\), not depending on the sample size \(n\)! — and its natural statistics are precisely \(\overline{u_1(X)},\ldots ,\overline{u_k(X)}\). These are exactly the statistics that \(S(X)\) is made up of in the theorem above: the theorem states that, under full rank, they are not merely sufficient, but complete as well.

5.2.2 Moments in Canonical Form

The canonical form of writing the density allows one to find the moments of the statistics \(u_i(X)\) without explicit integration — directly from the function \(v(\eta )\).

Proposition 1 (Moments via \(v(\eta )\)) Let the density of \(X\) have the canonical form \(p(x\mid \eta ) = h(x)\exp \left(\sum_{i=1}^k \eta_i u_i(x) + v(\eta )\right)\), and let \(\eta\) be an interior point of the natural parameter set \(\Xi\). Then \[ \mathbb {E}_{\eta }\left[u_j(X)\right] = -\frac{\partial v(\eta )}{\partial \eta _j}, \qquad \operatorname {Cov}\left[u_j(X), u_k(X)\right] = -\frac{\partial ^2 v(\eta )}{\partial \eta _j \partial \eta _k}. \]

This is proved by differentiating the identity \(\int p(x\mid \eta )\, \mu (dx) = 1\) with respect to \(\eta_j\) (we will not discuss the legitimacy of differentiating under the integral sign here):

\[ \int u_j(x)\, p(x\mid \eta )\, \mu (dx) + \frac{\partial v(\eta )}{\partial \eta _j} \underbrace{\int p(x\mid \eta )\, \mu (dx)}_{=1} = 0 \quad \Longrightarrow \quad \mathbb {E}_{\eta }\left[u_j(X)\right] = -\frac{\partial v(\eta )}{\partial \eta _j}. \]

The formula for the covariance is obtained by one more differentiation, this time with respect to \(\eta_k\) (carry this out yourself).

ExampleExample 12

The Poisson distribution \(\operatorname {Pois}\left(\lambda \right)\), \(\lambda > 0\), can be written in canonical form with \(\eta = \ln \lambda\), \(u(x) = x\), \(v(\eta ) = -e^{\eta }\). Using proposition Proposition 1, find \(\mathbb {E}\left[X\right]\) and \(\operatorname {Var}\left[X\right]\) for \(X \sim \operatorname {Pois}\left(\lambda \right)\).

SolutionSolution

Indeed, \(\mathbb {P}\left(X=x\right) = \frac{\lambda^x}{x!}e^{-\lambda } = \frac{1}{x!}\exp \left(x\ln \lambda - \lambda \right)\), from which we see that \(h(x)=1/x!\), \(\eta =\ln \lambda\), \(u(x)=x\), \(v(\eta ) = -\lambda = -e^{\eta }\). By proposition Proposition 1, \[ \mathbb {E}\left[X\right] = -v'(\eta ) = e^{\eta } = \lambda , \qquad \operatorname {Var}\left[X\right] = -v''(\eta ) = e^{\eta } = \lambda , \] which agrees with the well-known formulas \(\mathbb {E}\left[X\right] = \operatorname {Var}\left[X\right] = \lambda\) for the Poisson distribution — but obtained without a single integration or summation of a series.

Below are a few more classical one- and two-parameter exponential families.

Mn-remarkRemark

The family of normal distributions \(\mathscr {N}\left(\xi , \sigma^2\right)\) with fixed \(\sigma\) is simultaneously a (one-dimensional) location family and an exponential family. It turns out that these two examples essentially exhaust all families possessing both properties at once: Dynkin (1951) and Ferguson (1962) showed that the exponential location families are precisely \(\mathscr {N}\left(\xi , \sigma^2\right)\) (in the parameter \(\xi\), with \(\sigma\) fixed) and \(c\ln Y\), where \(Y \sim \Gamma (a,b)\), \(\theta = \ln b\) is the location parameter, and \(c \neq 0\) is an arbitrary fixed constant. Thus exponential and group (location-scale) families — two, at first glance, independent ways of constructing models — turn out to be closely related to one another.

5.3 Other Families

Let us finish with an example showing that not every natural parametric family is a group family or an exponential family.

ExampleExample 13

Let \(X\) have the zero-inflated Poisson distribution (ZIP) with parameters \(\pi \in [0,1)\), \(\lambda >0\): with probability \(\pi\) a “structural” zero is observed, and with probability \(1-\pi\) a value from \(\operatorname {Pois}\left(\lambda \right)\). That is, \[ p_\theta (0) = \pi + (1-\pi )e^{-\lambda }, \qquad p_\theta (k) = (1-\pi )\frac{\lambda ^k}{k!}e^{-\lambda }, \; k=1,2,\ldots , \] where \(\theta = (\pi ,\lambda ) \in \Theta = [0,1) \times (0,+\infty )\). Such a model is often used for count data with an “excess” number of zeros (number of doctor visits, number of insurance claims, etc.): part of the observations are guaranteed zeros (say, part of the insured population is not exposed to the risk at all), and the rest behave like an ordinary Poisson counter.

Does the family \(\mathcal{P}_{\mathrm{ZIP}} = \left\{ \mathrm{ZIP}(\theta ) \; : \; \theta \in \Theta \right\}\) form a group family? If so — indicate the generating group of transformations; if not — prove it.

SolutionSolution

No, it does not. Let us first establish a general fact, which we will then apply.

Let \(\mathcal X\) be a countable set (here \(\mathcal X = \left\{ 0,1,2,\ldots \right\}\)) and let \(g: \mathcal X \to \mathcal X\) be an arbitrary \(1{:}1\) transformation. In the discrete case, no additional restrictions such as continuity or monotonicity are imposed on \(g\) (unlike on \(\mathbb {R}\), any bijection of a countable set with the discrete topology is automatically continuous). If \(Q = g_{\sharp }P\), then by definition \(q(g(x)) = p(x)\) for all \(x \in \mathcal X\), i.e., the density \(q\) is the density \(p\) with the points at which it takes its values rearranged (according to the rule \(g\)). Hence, the multiset of values \(\left\{ p(x) \; : \; x \in \mathcal X\right\}\) does not change under the action of any bijection \(g\), and neither does any symmetric function of it, for example

\[ M(P) := \sup _{x \in \mathcal X} p(x). \]

If \(\mathcal P\) is a group family with group \(G\), i.e., \(\mathcal P = \left\{ g_\sharp \mathbb {P}_0 \; : \; g \in G\right\}\) for some \(\mathbb {P}_0 \in \mathcal P\), then for any \(P_1 = g_{1\sharp }\mathbb {P}_0,\, P_2 = g_{2\sharp }\mathbb {P}_0 \in \mathcal P\) the element \(g := g_2 g_1^{-1} \in G\) carries \(P_1\) into \(P_2\). From this we obtain a useful general principle: for a group family on a countable \(\mathcal X\), the quantity \(M(P)\) must be the same for all \(P \in \mathcal P\).

Let us check this for \(\mathcal P_{\mathrm{ZIP}}\), taking two specific members of the family.

• \(\theta_1 = (0,1)\) (ordinary \(\operatorname {Pois}\left(1\right)\), with no excess zeros): \(p_{\theta_1}(0) = p_{\theta_1}(1) = e^{-1} \approx 0{,}368\), and this is the maximum (for \(k\geq 2\) the values \(p_{\theta_1}(k)\) only decrease). Hence, \(M(P_{\theta_1}) = e^{-1} \approx 0{,}368\).

• \(\theta_2 = (1/2,1)\): \(p_{\theta_2}(0) = \frac12 + \frac12 e^{-1} \approx 0{,}684\), whereas \(p_{\theta_2}(k) = \frac12 p_{\theta_1}(k) \leq \frac12 e^{-1} \approx 0{,}184\) for \(k \geq 1\). Hence, \(M(P_{\theta_2}) = p_{\theta_2}(0) \approx 0{,}684\).

So, \(M(P_{\theta_1}) \neq M(P_{\theta_2})\), even though both \(\theta_1,\theta_2 \in \Theta\). By the principle above, no bijection \(\mathcal X \to \mathcal X\) can carry \(\mathrm{ZIP}(\theta_1)\) into \(\mathrm{ZIP}(\theta_2)\) — and hence, all the more, there is no group of transformations \(G\) that could connect these two distributions. The family \(\mathcal P_{\mathrm{ZIP}}\) is not a group family.

Substantively this is not surprising: the parameter \(\pi\) is not responsible for a symmetric transformation of the support (like a shift or a scaling), but for mixing in the degenerate distribution \(\delta_0\), which qualitatively changes the shape of the distribution rather than merely permuting its values — the mode of \(\operatorname {Pois}\left(\lambda \right)\) can be at any point \(\lfloor \lambda \rfloor\), whereas for a ZIP with large \(\pi\) the entire mass concentrates at the point \(0\) regardless of \(\lambda\).

Footnotes

1. \(\xi\) can also be a random vector taking values in \(\mathbb {R}^{k}\), with an unknown \(k\)-dimensional distribution \(\mathbb {P}\). In this case the model is constructed analogously. Going forward, for simplicity of exposition, the theory is presented for the one-dimensional case.

2. Recall that a map is called measurable if for every set \(A \in \mathcal{E}\) we have \(\left\{ x \; : \; S(x) \in A\right\} \in \mathscr {B}\left(\mathbb {R}^{n}\right)\)

3. One can even prove that a sample of infinite size from an arbitrary distribution can be defined on \(([0,1], \mathscr {B}\left([0,1]\right), \operatorname {U}[0,1])\).