Methods of Estimation — A Guided Tour of Modern Statistics

Notes / Parameter Estimation

2.2 Methods of Estimation

Now we introduce loss function to describe the quality of estimators quantitatively. Consider \(\theta \) as the decision rule, and \(R(\theta ,\delta ) = {\mathbb{E}} _\theta L(\theta , \delta (X))\) as the risk function.

2.2.1 Unbiased Estimation

An estimator is unbiased if \({\mathbb{E}} _\theta [\delta (X)] = g(\theta )\) (to estimate). Although uniformly best estimator does not exist, quite often, we can find a unbiased estimator with uniformly minimum risk, that is, an unbiased \(\delta \) satisfying \(R(\theta , \delta ) \le R(\theta , \delta ')\), \(\forall \theta \) and any other unbiased estimators \(\delta '\). Such an estimator is called a uniformly minimum risk unbiased estimator (UMRUE).

Consider \(L(\theta , d) = (\theta -d)^2\), then an UMRUE becomes a uniformly minimum variance unbiased estimator (UMVUE). \begin{align*} \mathbb {E}_\theta \left [(\theta - \delta (X))^2\right ] &= \left (\mathbb {E}_\theta [\delta (X)] - \theta \right )^2 + \mathbb {E}_\theta \left \{\left (\delta (X) - \mathbb {E}_\theta [\delta (X)]\right )^2\right \} \\ &= \text {Bias}^2 + \text {Variance}. \end{align*}

Example 2.16. Let \(X_1, \ldots , X_n \stackrel {iid}{\sim } \mathcal {N}(\mu , \sigma ^2)\). First, note that if \(\sigma ^2\) is known, \(\bar {X}\) is a complete sufficient statistic for \(\mu \) and hence also the UMVUE. Consider the case when \(\theta = (\mu , \sigma ^2)\) is unknown.

(a): The UMVUE for \(\mu \) is \(\bar {X}\).
(b): The UMVUE for \(\sigma ^2\) is \(\sum _{i=1}^{n} \frac {(X_i - \bar {X})^2}{n-1}\).
(c): What is the UMVUE for \(\sigma \)? First, note that \(X_i - \bar {X} \sim \mathcal {N}(0, \frac {n-1}{n} \sigma ^2)\), and hence \(\mathbb {E}[|X_i - \bar {X}|] = \sigma \sqrt {\frac {2}{\pi }} \sqrt {\frac {n-1}{n}}\). This implies \[ \frac {\sqrt {\pi n}}{\sqrt {2(n-1)}} |X_i - \bar {X}| \] is unbiased for \(\sigma \). At this point we could Rao-Blackwellize, but the math is messy. Instead, we will try to stumble upon the solution. Let \[ S^2 = \sum _{i=1}^{n} (X_i - \bar {X})^2. \] We know that \[ S^2 \sim \sigma ^2 \chi ^2_{n-1}. \] Thus, \[ \mathbb {E}(S) = \sigma \mathbb {E}(\chi _{n-1}). \] Which in turn implies that \[ \frac {\mathbb {E}(S)}{\mathbb {E}(\chi _{n-1})} = \sigma , \] meaning \(\frac {S}{\mathbb {E}(\chi _{n-1})}\) is unbiased for \(\sigma \) and hence UMVU.
(d): What is the UMVUE for \(\mu ^2\)? Taking the expectation of the UMVUE for \(\mu \) and squaring it yields \[ \mathbb {E}(\bar {X}^2) = \mu ^2 + \sigma ^2 / n. \] So, \[ \delta _n(X) = \bar {X}^2 - \frac {S^2}{n(n-1)} \] is the UMVUE. Note that \(\delta _n(X)\) may be negative even though it estimates a non-negative quantity. Indeed, \(\delta _n\) is inadmissible and dominated by the biased estimator \(\max (0, \delta _n(X))\).

If we know some prior knowledge of the estimators (eg: Equivariant, Invariant), we can design better estimators.

2.2.2 Bayes Estimators

Our optimality goal, given a measure \(\Lambda \), is to find an estimator \(\delta _\Lambda \) which minimizes the average risk, \[ r(\Lambda , \delta ) = \int R(\theta , \delta ) \, d\Lambda (\theta ). \] If \(\Lambda \) is a probability distribution on \(\Omega \), we call \(\Lambda \) the prior distribution. The estimator \(\delta _\Lambda \), if it exists, is called the Bayes estimator with respect to \(\Lambda \), and the minimized average risk \(r(\Lambda , \delta _\Lambda )\) is called the Bayes risk.

Example 2.18. Suppose that \(X \sim \text {Bin}(n, \theta )\) given \(\Theta = \theta \) and that \(\Theta \) has prior distribution \(\text {Beta}(a, b)\). The prior density is given by \[ \pi (\theta ) = \frac {\Gamma (a+b)}{\Gamma (a)\Gamma (b)} \theta ^{a-1} (1-\theta )^{b-1} \mathbb {I}(0 < \theta < 1) \] The likelihood (model density) is given by \[ f(x \mid \theta ) = \binom {n}{x} \theta ^x (1-\theta )^{n-x} \] The marginal density is given by \[ f(x) = \int f(x \mid \theta ) \pi (\theta ) \, d\theta . \] The posterior density may be calculated using Bayes rule which states that \[ \text {posterior} = \frac {\text {joint}}{\text {marginal}} = \frac {\text {prior} \cdot \text {likelihood}}{\text {marginal}}. \] In our notation, the posterior density is given by the formula \begin{align*} \pi (\theta \mid x) &= \frac {\pi (\theta ) f(x \mid \theta )}{f(x)} \\ &= \frac {\pi (\theta ) f(x \mid \theta )}{\int \pi (\theta ') f(x \mid \theta ') \, d\theta '} \end{align*}

\begin{align*} \pi (\theta \mid x) &\propto \binom {n}{x} \theta ^x (1-\theta )^{n-x} \frac {\Gamma (a+b)}{\Gamma (a)\Gamma (b)} \theta ^{a-1} (1-\theta )^{b-1} \\ &\propto \theta ^{x+a-1} (1-\theta )^{n-x+b-1} \sim \text {Beta}(x+a, n-x+b). \end{align*}

Hence, the Bayes estimator of \(\theta \) under the squared error loss is given by \[ \mathbb {E}[\Theta \mid X = x] = \frac {x+a}{n+a+b}. \] This posterior mean may be expressed as \[ \frac {X+a}{n+a+b} = \frac {n}{n+a+b} \left (\frac {X}{n}\right ) + \frac {a+b}{n+a+b} \left (\frac {a}{a+b}\right ). \] Hence, the Bayes estimate is a convex combination of the sample proportion \(X/n\) (which is the UMVUE) and the prior mean \(a/(a+b)\). Thus, the Bayes estimate modifies the sample estimate in light of prior information by “shrinking” the sample estimate towards the prior mean. (This is a commonly recurring property of Bayes estimators.) In addition, as the sample size \(n\) tends to infinity, the weight of the prior mean tends to zero, the empirical evidence increasingly outweighs the prior information, and the posterior mean becomes less distinguishable from the sample proportion.

2.2.3 Minimax Estimators

In minimax estimation, we collapse our risk function by looking at the worse-case risk. Given \(X \sim \mathbb {P}_\theta \), where \(\theta \in \Omega \), and a loss function \(L(\theta , d)\), we want to find an estimator \(\delta \) that minimizes the maximum risk: \[ \sup _{\theta \in \Omega } R(\theta , \delta ). \] Any such \(\delta \) is called a minimax estimator.

Example 2.28. Suppose \(X \sim \text {Binom}(n, \theta )\) for some \(\theta \in (0,1)\) and that we use the squared error loss function. Is the sample proportion \(\frac {X}{n}\) minimax? The risk of this estimator is \[ R\left (\theta , \frac {X}{n}\right ) = \frac {\theta (1-\theta )}{n}. \] The graph of \(R\left (\theta , \frac {X}{n}\right )\) versus \(\theta \) looks like the following:

The risk has a unique maximum at \(\theta = \frac {1}{2}\), so the worst-case risk is \[ \sup _{\theta \in \Omega } R\left (\theta , \frac {X}{n}\right ) = R\left (\frac {1}{2}, \frac {X}{n}\right ) = \frac {1}{4n}. \]

We can use the Corollary 2.27 to find a minimax estimator and then compare the risk of the minimax estimator with that of \(\frac {X}{n}\). To find a minimax estimator, we will search for a prior such that the Bayes estimator has constant risk.

Recall the following useful fact. Under the prior distribution \(\text {Beta}(a, b)\), the Bayes estimator under the squared error loss is \[ \delta _{a,b}(X) = \frac {X + a}{n + a + b}. \] For any \(a\) and \(b\), \begin{align*} R(\theta , \delta _{a,b}) &= \mathbb {E}_\theta \left [\left (\frac {X + a}{n + a + b} - \theta \right )^2\right ] \\ &= \frac {1}{(n + a + b)^2} \mathbb {E}_\theta \left [(X + a - (n + a + b)\theta )^2\right ] \\ &= \frac {1}{(n + a + b)^2} \mathbb {E}_\theta \left [(X - n\theta - a(\theta - 1) - \theta b)^2\right ] \\ &= \frac {1}{(n + a + b)^2} \left (n\theta (1 - \theta ) + (a(\theta - 1) + \theta b)^2\right ). \end{align*}

This is a quadratic function of \(\theta \). To eliminate the \(\theta \) dependence in \(R(\theta , \delta _{a,b})\), we need the coefficients of the linear and quadratic terms to equal zero. The coefficient of \(\theta ^2\) is \[ -n + (a + b)^2, \] so we need \(a + b = \sqrt {n}\) (since \(a, b > 0\)). The coefficient of \(\theta \) is \[ n - 2a(a + b) = n - 2a\sqrt {n}, \] so we need \(a = b = \frac {\sqrt {n}}{2}\). With these choices of \(a\) and \(b\), the risk of \(R(\theta , \delta _{a,b})\) is constant, which implies that \(\text {Beta}\left (\frac {\sqrt {n}}{2}, \frac {\sqrt {n}}{2}\right )\) is a least favorable prior with constant risk. Then our Bayes estimator \[ \delta _{\frac {\sqrt {n}}{2}, \frac {\sqrt {n}}{2}}(X) = \frac {X + \frac {\sqrt {n}}{2}}{n + \sqrt {n}}, \] is minimax with constant risk of \[ \frac {1}{4(\sqrt {n} + 1)^2}. \] Since the worst-case risk of \(\frac {X}{n}\) is \(\frac {1}{4n} > \frac {1}{4(\sqrt {n}+1)^2}\), we can conclude that \(\frac {X}{n}\) is not minimax.

2.2.4 James Stein Estimator

Example 2.29 (James Stein Estimator). Let \(X_1, X_2, \ldots , X_p\) be independent with \(X_i \sim \mathcal {N}(\theta _i, \sigma ^2)\) for \(1 \leq i \leq p\). For the sake of simplicity, say \(\sigma ^2 = 1\). Now our goal is to estimate \(\theta = (\theta _1, \theta _2, \ldots , \theta _p)\) under the loss function: \[ L(\theta , d) = \sum _{i=1}^{p} (d_i - \theta _i)^2 \]

A natural estimator for \(\theta \) is \(X = (X_1, X_2, \ldots , X_p)\). It can be shown that \(X\) is the UMRUE, the maximum likelihood estimator, a generalized Bayes estimator, and a minimax estimator for \(\theta \). So, it would be natural to think that \(X\) is admissible. However, counterintuitively, it turns out that this is not the case when \(p \geq 3\).

When \(p \geq 3\), \(X\) is dominated by the James-Stein estimator (and that too, strictly dominated): \[ \delta (X) = (\delta _1(X), \delta _2(X), \ldots , \delta _p(X)) \text { where} \] \[ \delta _i(X) = \left (1 - \frac {p-2}{\|X\|_2^2}\right ) X_i. \] It turns out that the James-Stein estimator is not itself admissible because it is dominated by the positive part James-Stein estimator \[ \delta _i(X) = \max \left ( 1 - \frac {p-2}{\|X\|_2^2}, 0 \right ) X_i. \]

Remark 2.30. Intuitively, the problem with the estimate \(X\) is that \(\|X\|_2^2\) is typically much larger than \(\|\theta \|_2^2\): \[ \mathbb {E}[\|X\|_2^2] = E\left [\sum _{i=1}^{p} X_j^2\right ] = p + \sum _{i=1}^{p} \theta _i^2 = p + \|\theta \|_2^2 \] where \(p\) is actually \(\sigma ^2 p = p\) in this case. So, we may view the J-S estimator as a method for correcting the bias in the size of \(X\). It achieves this by shrinking each coordinate of \(X\) toward 0.

1.

Suppose \(\theta _i \stackrel {iid}{\sim } \mathcal {N}(0, A)\) then the Bayes estimator for \(\theta _i\) is \[ \delta _{A,i}(X) = \frac {X_i}{1 + \frac {1}{A}} = \left (1 - \frac {1}{A+1}\right ) X_i \]

2.

In this step we must choose \(A\). Marginalizing over \(\theta \), we see that \(X\) has the distribution, \[ X_i \stackrel {iid}{\sim } \mathcal {N}(0, A+1) \] We will use \(X\) and the knowledge of this marginal distribution to find an estimate of \(\frac {1}{A+1}\). One could, in principle, use any estimate of \(A\), and it is common to use a maximum likelihood estimate, but here we will used an unbiased estimate.

It can then be shown that \[ \mathbb {E}\left [\frac {1}{\|X\|_2^2}\right ] = \frac {1}{(p-2)(A+1)} \] (\(\frac {1}{A+1}\|X\|_2^2\) follows a \(\chi _n^2\) distribution). So \[ 1 - \frac {p-2}{\|X\|_2^2} \] must be UMVU for \(1 - \frac {1}{A+1}\).

If we plug this estimator into our Bayes estimator we obtain the J-S estimator: \[ \delta (X_i) = \left (1 - \frac {p-2}{\|X\|_2^2}\right ) X_i. \]

2.2.5 How to Derive an Estimator

Example 2.31 (Moment Estimation). Let \(X_1, \ldots , X_n\) be i.i.d. random variables from \(P_\theta \), \(\theta \in \Theta \subset \mathbb {R}^k\), and assume \(\mathbb {E}|X_1|^k < \infty \).

Let \(\mu _j = \mathbb {E} X_1^j\) be the \(j\)th moment of \(P\) and then \[ \hat {\mu }_j = \frac {1}{n} \sum _{i=1}^{n} X_i^j \] becomes the \(j\)th sample moment, which is an unbiased estimator of \(\mu _j\), \(j = 1, \ldots , k\).
Typically, \begin{equation} \mu _j = h_j(\theta ), \quad j = 1, \ldots , k \end{equation} for some functions \(h_j\) on \(\mathbb {R}^k\). By substituting \(\mu _j\)’s on the left-hand side of (5) by the sample moments \(\hat {\mu }_j\), we obtain a moment estimator \(\hat {\theta }\), i.e., \(\hat {\theta }\) satisfies \[ \hat {\mu }_j = h_j(\hat {\theta }), \quad j = 1, \ldots , k, \] which is a sample analogue of (5). This method of deriving estimators is called the method of moments.
Let \(\hat {\mu } = (\hat {\mu }_1, \ldots , \hat {\mu }_k)\) and \(h = (h_1, \ldots , h_k)\). Then \(\hat {\mu } = h(\hat {\theta })\). If the inverse function \(h^{-1}\) exists, then the unique moment estimator of \(\theta \) is \(\hat {\theta } = h^{-1}(\hat {\mu })\).
When \(h^{-1}\) does not exist (i.e., \(h\) is not one-to-one), any solution of \(\hat {\mu } = h(\hat {\theta })\) is a moment estimator of \(\theta \). If possible, we always choose a solution \(\hat {\theta }\) in the parameter space \(\Theta \). In some cases, however, a moment estimator may not exist.

Example 2.32 (Maximum Likelihood Estimation). Let \(X \in \mathcal {X}\) be a sample with a PDF \(f_{\boldsymbol {\theta }}\) w.r.t. a \(\sigma \)-finite measure, where \(\boldsymbol {\theta } \in \Theta \subset \mathbb {R}^k\).

For each \(x \in \mathcal {X}\), \(f_{\boldsymbol {\theta }}(x)\) considered as a function of \(\boldsymbol {\theta }\) is called the likelihood function and denoted by \(\ell (\theta )\).
Let \(\bar {\Theta }\) be the closure of \(\Theta \). A \(\hat {\theta } \in \bar {\Theta }\) satisfying \(\ell (\hat {\theta }) = \max _{\theta \in \bar {\Theta }} \ell (\theta )\) is called a maximum likelihood estimate (MLE) of \(\theta \). If \(\hat {\theta }\) is a measurable function of \(X\), then \(\hat {\theta }\) is called a maximum likelihood estimator (MLE) of \(\theta \).