Theorem 1.2.9 and Bonferroni ’s Inequality

If \(P\) is a probability function and \(A\) and \(B\) are any sets in \(B\), then
- a. \(P(B\cap A^c)=P(B)-P(A\cap B)\);
- b. \(P(A\cup B)=P(A)+P(B)-P(A\cap B)\);
- c. iF \(A\subset B\) then \(P(A)\le P(B)\);

Proof: To establish (a) note that for any sets \(A\) and \(B\) we have \[B=\{B\cap A\}\cup\{B\cap A^c\}\] and therefor \[P(B)=P(\{B\cap A\}\cup\{B\cap A^c\})=P(B\cap A)+P(B\cap A^c)\] since \(B\cap A\) and \(B\cap A^c\) are disjoint.

To establish (b), we use the identity \[A\cup B=A\cup \{B\cap A^c\}\] since \(A\) and \(B\cap A^c\) are disjoint, we have \[P(A\cup B)=P(A)+P(B\cap A^c)=P(A)+P(B)-P(A\cap B)\] from (a).

If \(A\subset B\), then \(A\cap B=A\). Therefore, using (a) we have \[0\le P(B\cap A^c)=P(B)-P(A)\] establishing (c).

Formula (b) of Theorem 1.2.9 gives a useful inequality for the probability of an intersection. Since \(P(A\cup B)\le 1\), we have from \[P(A\cup B)=P(A)+P(B\cap A^c)=P(A)+P(B)-P(A\cap B)\le 1\], after some rearranging, \[P(A\cap B)\ge P(A)+P(B)-1\] This inequality is a special case of what is known as Bonferroni ’s Inequality.

Boole’s Inequality

If \(P\) is a probability junction, then

• a. \(P(A)=\sum_{i=1}^{\infty}P(A\cap C_i)\) for any partition \(C_1,C_2,\cdots\);
  
• b. \(P(\bigcup_{i=1}^{\infty}A_i)\le\sum_{i=1}^{\infty}P(A_i)\) for any sets \(A_1,A_2,\cdots\).

Proof: Since \(C_1,C_2,\cdots\) form a partition, we have that \(C_i\cap C_j=\emptyset\) for all \(i \ne j\), and \(S=\bigcup_{i=1}^{\infty}C_i\). Hence, \[A=A\cap S=A\cap\left(\bigcup_{i=1}^{\infty}C_i\right)=\bigcup_{i=1}^{\infty}(A\cap C_i)\] where the last equality follows from the Distributive Law. We therefore have \[P(A)=P\left(\bigcup_{i=1}^{\infty}(A\cap C_i)\right)\] Now, since the \(C_i\) are disjoint, the sets \(A\cap C_i\) are also disjoint, and from the properties of a probability function we have \[P\left(\bigcup_{i=1}^{\infty}(A\cap C_i)\right)=\sum_{i=1}^{\infty}P(A\cap C_i),\] establishing (a).

To establish (b), construct a disjoint collection \(A_1^*,A_2^*,\cdots\), with the property that \(\cup_{i=1}^{\infty}A_i^*=\cup_{i=1}^{\infty}A_i\), we define \(A_i^*\) by \[A_1^*=A_1, A_i^*=A_i\setminus\left(\bigcup_{j=1}^{i-1}A_j\right),\quad i=2,3,\cdots,\] where \(A\setminus B=A\cap B^c\). It should be easy to see that \(\cup_{i=1}^{\infty}A_i^*=\cup_{i=1}^{\infty}A_i\) and we have \[P\left(\bigcup_{i=1}^{\infty}A_i\right)=P\left(\bigcup_{i=1}^{\infty}A_i^*\right)=\sum_{i=1}^{\infty}P(A_i^*)\le\sum_{i=1}^{\infty} P(A_i)\] Since, by construction \(A_i^*\subset A_i\). establishing (b).

Bonferroni ’s Inequality

If we apply Boole’s Inequality to \(A^c\) of Bonferroni ’s Inequality, we have \[P\left(\bigcup_{i=1}^{n}A_i^c\right)\le\sum_{i=1}^{n}P(A_i^c)\] and using the fact that \(\bigcup_{i=1}^{n}A_i^c=\left(\bigcap_{i=1}^{n} A_i\right)^c\) and \(P(A_i^c)=1-P(A_i)\) so we obtain \[1-P\left(\bigcap_{i=1}^{n} A_i\right)=P\left(\bigcup_{i=1}^{n}A_i^c\right)\le\sum_{i=1}^{n}P(A_i^c)= n-\sum_{i=1}^{n}P(A_i)\] This becomes \[P\left(\bigcap_{i=1}^{n} A_i\right)\ge\sum_{i=1}^{n}P(A_i)-(n-1)\] which is a more general version of the Bonferroni Inequality.

Probability integral transformation

Let \(X\) have continuous cdf \(F_{X}(x)\) and define the random variable \(Y\) as \(Y = F_{X}(X)\). Then \(Y\) is uniformly distributed on \((0, 1)\), that is, \(P(Y \le y) = y, 0 < y < 1\).

Let \(F_{X}^{-1}\) be the inverse of the cdf \(F_{X}\). If \(F_{X}\) is strictly increasing, then \(F_{X}^{-1}\) is well defined by \[F_{X}^{-1}(y)=x\iff F_{X}(x)=y.\] However, if \(F_{X}\) is constant on some interval \(x_1\le x \le x_2\), then \(F_{X}^{-1}\) is not well defined by \(F_{X}^{-1}(y)=x\iff F_{X}(x)=y.\) Any \(x\) satisfying \(x_1\le x \le x_2\) satisfies \(F_{X}(x) = y\). This problem is avoided by defining \(F_{X}^{-1}(y)\) for \(0 < y < 1\) by \[F_{X}^{-1}(y) = \inf\; \{x: F_{X}(x) \ge y\},\] Using this definition, we have \(F_{X}^{-1}(y)=x_1\). At the endpoints of the range of \(y\), \(F_{X}^{-1}(y)\) can also be defined. \(F_{X}^{-1}(1)=\infty\) if \(F_X(x) < 1\) for all \(x\) and, for any \(F_X\), \(F_{X}^{-1}(0) = -\infty.\)

Proof: For \(Y = F_X (X)\) we have, for \(0 < y < 1\), \[\begin{align} P(Y\le y)&=P(F_X (X)\le y)\\ &=P(F_X^{-1}[F_X (X)]\le F_X^{-1}(y))\quad\quad(F_X^{-1} \text{ is increasing})\\ &=P(X\le F_X^{-1}(y))\\ &=F_X(F_X^{-1}(y))\\ &=y\\ \end{align}\] At the endpoints we have \(P(Y\le y)=1\) for \(y \ge 1\) and \(P (Y\le y) = 0\) for \(y \le 0\), showing that \(Y\) has a uniform distribution.

Gamma Distribution

\[f(x|\alpha,\beta)=\frac{1}{\Gamma(\alpha)\beta^{\alpha}}x^{\alpha-1}e^{-x/\beta}\quad 0<x<\infty;\quad\alpha,\beta>0\]

Normal Distribution

If \(\mathbf{X} \sim \text{n}(\mathbf{\mu}, \boldsymbol{\Sigma})\) then \(\mathbf{X}\) has pdf \[f_{\mathbf{X}}(x_1, \dots, x_k \vert{} \mathbf{\mu}, \boldsymbol{\Sigma}) = \frac{1}{\sqrt{(2\pi)^k \det\boldsymbol{\Sigma}}} \exp\left(-\frac12 (\mathbf{ x}-\mathbf{\mu})^\top \boldsymbol{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})\right).\]

Chi-Squared Distribution

The chi-squared distribution with \(p\) degrees of freedom has pdf \[\chi_p^2 \sim \frac{1}{\Gamma(p/2)2^{p/2}} x^{(p/2)-1} e^{-x/2}, \quad 0< x< \infty.\]

[Some facts]

1. If \(Z \sim \text{n}(0, 1)\) then \(Z^2 \sim \chi_1^2\)

2. If \(X_1, \dots, X_n\) are independent \(X_i \sim \chi_{p_i}^2\) then \(X_1 + \cdots + X_n \sim \chi_{p_1 + \cdots + p_n}^2\).

Student’s \(t\)-Distribution

Student’s \(t\)-distribution \(T \sim t_p\), a \(t\)-distribution with \(p\) degrees of freedom if it has pdf \[f_T(t) = \frac{\Gamma\left(\frac{p+1}{2}\right)}{\Gamma\left(\frac{p}{2}\right)} \frac{1}{\sqrt{p\pi}} \frac{1}{(1 + t^2/p)^{(p+1)/2}}, \quad t \in \mathbb R{}\]

If \(p = 1\) then this is the Cauchy distribution.

If \(X_1, \dots, X_n\) are a random sample from \(\text{n}(\mu, \sigma^2)\) then \[\frac{\bar{X} - \mu}{S/\sqrt{n}} \sim t_{n-1}.\] This is often taken as the definition. Note that the denominator is independent of the numerator.

[Moments and mgf of \(t\)-distribution] Student’s \(t\) has no mgf because it does not have moments of all orders: \(t_p\) has only \(p-1\) moments. If \(T_p \sim t_p\) then

\[\begin{aligned} \mathbb E{}[T_p] &= 0 \quad p > 1 \\ \text{Var}{}[T_p] &= \frac{p}{p-2} \quad p > 2 \end{aligned}\]

Snedcor’s \(F\)-Distribution

Snedcor’s \(F\)-distribution A random variable \(X \sim F_{p,q}\) has \(F\)-distribution with p and q degrees of freedom if its pdf is \[f_X(x) = \frac{\Gamma\left(\frac{p+q}{2}\right)}{\Gamma\left(\frac p2 \right)\Gamma\left(\frac q2 \right)} (p/q)^{p/2} \frac{x^{(p/2) - 1}}{(1 + px/q)^{(p+q)/2}}, \quad 0 < x < \infty.\]

If \(X_1, \dots, X_n\) is a random sample from \(\text{n}(\mu_X, \sigma_X^2)\) and \(Y_1, \dots, Y_m\) is an independent random sample from \(\text{n}(\mu_Y, \sigma_Y^2)\), then \[\frac{S^2_X/\sigma^2_X}{S_Y^2/\sigma_Y^2} \sim F_{n-1, m-1}.\] This is often taken as the definition.

[Some facts]

1. \(X \sim F_{p,q} \,\, \implies \,\, 1/X \sim F_{q,p}\)

2. \(X \sim t_q \,\, \implies \,\, X^2 \sim F_{1,q}\)

3. \(X \sim F_{p,q} \,\, \implies \,\, \frac{(p/q)X}{1 + (p/q)X} \sim \text{beta}(p/2, q/2)\)

Multinomial Distribution

Multinomial Distribution Let \(m\) and \(n\) be positive integers and let \(p_1, \dots, p_n \in [0, 1]\) satisfy \(\sum_{i=1}^n p_i = 1\). Then the random vector \((X_1, \dots, X_n)\) has multinomial distribution with m trials and cell probabilities \(p_1, \dots, p_n\) if the joint pmf of \((X_1, \dots, X_n)\) is \[f(x_1, \dots , x_n) = \frac{m!}{x_1! \cdots x_n!} p_1^{x_1} \cdots p_n^{x_n} = m! \prod_{i=1}^n \frac{p_i^{x_i}}{x_i!}\] on the set of \((x_1, \dots, x_n)\) such that each \(x_i\) is a nonnegative integer and \(\sum_{i=1}^n x_i = m\).

The marginal distributions have \(X_i \sim \text{binomial}(m, p_i)\).

[Multinomial Theorem] Let \(m\) and \(n\) be positive integers and let \(\mathcal{A}\) be the set of vectors \(\mathbf{x} = (x_1, \dots, x_n)\) such that each \(x_i\) is a nonnegative integer and \(\sum_{i=1}^n x_i = m\). Then for any real numbers \(p_1, \dots, p_n\), \[(p_1 + \cdots + p_n)^m = \sum_{\mathbf{x} \in \mathcal{A}} \frac{m!}{x_1! \cdots x_n!} p_1^{x_1}\cdots p_n^{x_n}.\]

Chebyshev’s inequality

Let \(X\) be a random variable and let \(g(x)\) be a nonnegative function. Then, for any \(r > 0\), \[\text{P}(g(X) \geq r) \leq \frac{\mathbb E{}[g(X)]}{r}.\]

Proof: \[\begin{align} \mathbb E g(X)&=\int_{-\infty}^{\infty}g(x)f_{X}(x)dx\\ &\ge\int_{x:g(x)\ge r}g(x)f_{X}(x)dx\\ &\ge r\int_{x:g(x)\ge r}f_{X}(x)dx\\ &=r P(g(x)\ge r) \end{align}\]

This bound is conservative and almost never attained.

Chebyshev’s inequality (variance version)

The Markov inequality is the special case with \(g = \frac{(X-\mu)^2}{\sigma^2}\), where \(\mu=\mathbb E X\) and \(\sigma^2=\text{Var} X\). For convenience write \(r=t^2\).

\[\text{P}\left(\frac{(X-\mu)^2}{\sigma^2} \geq t^2\right) \leq \frac{1}{t^2}.\]

Let \(X_{\alpha, \beta}\) denote a gamma\((\alpha, \beta)\) random variable with pdf \(f(x \vert{} \alpha, \beta)\), where \(\alpha > 1\). Then for any constants \(a\) and \(b\): \[\text{P}(a < X_{\alpha, \beta} < b) = \beta (f(a \vert{} \alpha, \beta) - f(b \vert{} \alpha, \beta)) + \text{P}(a < X_{\alpha - 1, \beta} < b)\]

[Stein’s Lemma] Let \(X \sim \text{n}(\theta, \sigma^2)\) and let \(g\) be a differentiable function with \(\mathbb E{}[g'(x)] < \infty\). Then \[\mathbb E{}[g(X)(X - \theta)] = \sigma^2 \mathbb E{}[g'(X)]\]

The proof is just integration by parts.

Stein’s lemma is useful for moment calculations

Let \(\chi^2_p\) denote a chi squared distribution with \(p\) degrees of freedom. For any function \(h(x)\), \[\mathbb E{}[h(\chi^2_p)] = p \mathbb E{}\left[\frac{h(\chi_{p+2}^2)}{\chi_{p+2}^2}\right]\] provided the expressions exist.

Let \(g(x)\) be a function that is bounded at \(-1\) and has finite expectation, then

1. If \(X \sim \text{Poisson}(\lambda)\), \[\mathbb E{}[\lambda g(X)] = \mathbb E{}[Xg(X-1)].\]

2. If \(X\sim \text{negative-binomial}(r, p)\), \[\mathbb E{}[(1-p)g(X)] = \mathbb E{}\left[ \frac{X}{r + X - 1}g(X) \right].\]

Beta-binomial hierarchy

One generalization of the binomial distribution is to allow the success probability to vary according to a distribution. A standard model for this situation is \[X|P\sim\text{binomia1}(P),\quad i=1,\cdots,n,\] \[P\sim\text{beta}(\alpha,\beta).\] By iterating the expectation, we calculate the mean of \(X\) as \[\mathbb E X=\mathbb E [\mathbb E (X|P)]=\mathbb E [nP]=n\frac{\alpha}{\alpha+\beta}\]

Conditional variance identity

For any two random variables \(X\) and \(Y\), \[\text{Var} X=\mathbb E (\text{Var} (X|Y))+\text{Var} (\mathbb E (X|Y))\] provided that the expectations exist.

Proof: By definition, we have \[\begin{align} \text{Var} X&=\mathbb E([X-\mathbb E X]^2)\\ &=\mathbb E([X-\mathbb E(X|Y)+\mathbb E(X|Y)-\mathbb E X]^2)\\ &=\mathbb E([X-\mathbb E(X|Y)]^2)+\mathbb E([\mathbb E(X|Y)-\mathbb E X]^2)+2\mathbb E([X-\mathbb E(X|Y)][\mathbb E(X|Y)-\mathbb E X])\\ &=\mathbb E([X-\mathbb E(X|Y)]^2)+\mathbb E([\mathbb E(X|Y)-\mathbb E X]^2)\\ &=\mathbb E (\mathbb E\{[X-\mathbb E(X|Y)]^2|Y\})+\text{Var} (\mathbb E (X|Y))\\ &=\mathbb E (\text{Var} (X|Y))+\text{Var} (\mathbb E (X|Y)) \end{align}\]

Young’s Inequalities

Let \(a\) and \(b\) be any positive numbers and let \(p, q > 1\) satisfy \(1/p + 1/q = 1\), then \[p+q=pq\] \[q(p-1)=p\] \[1-1/q=1/p\] \[\frac1p a^p + \frac1q b^q \geq ab\] with equality if and only if \(a^p = b^q\).

Proof

Put \(\displaystyle t=1/p\) and \(\displaystyle (1-t)=1/q.\) Because the logarithm function is concave,

\[\displaystyle \ln \left(ta^{p}+(1-t)b^{q}\right)\geq t\ln \left(a^{p}\right)+(1-t)\ln \left(b^{q}\right)\\ =tp\ln(a)+(1-t)q\ln(b)\\ =\ln(a)+\ln(b)=\ln(ab)\] with the equality holding if and only if \(\displaystyle a^{p}=b^{q}.\) Young’s inequality follows by exponentiating.

H\(\ddot{o}\)lder’s Inequality Let \(X\) and \(Y\) be any random variables and let \(p, q > 1\) satisfy \(1/p + 1/q = 1\), then \[\lvert{\mathbb E{}[XY]}\lvert \leq \mathbb E{}[\lvert{XY}\lvert] \le \mathbb E{}[\lvert{X}\lvert^p]^{1/p} \mathbb E{}[\lvert{Y}\lvert^q]^{1/q}\] Proof: The first inequality follows from \(- |XY| \leq XY \leq |XY|\). To prove the second inequality, define \[a=\frac{|X|}{(\mathbb E|X|^p)^{1/p}}\quad\text{and}\quad b=\frac{|Y|}{(\mathbb E|Y|^q)^{1/q}}\] then \[\frac1p a^p+ \frac1q b^q=\frac1p\frac{|X|^p}{(\mathbb E|X|^p)}+\frac1q\frac{|Y|^q}{(\mathbb E|Y|^q)}\ge ab=\frac{|XY|}{(\mathbb E|X|^p)^{1/p}(\mathbb E|Y|^q)^{1/q}}\] Now take expectations of both sides. The expectation of the left-hand side is \(1\), then \[1\ge \frac{\mathbb E|XY|}{(\mathbb E|X|^p)^{1/p}(\mathbb E|Y|^q)^{1/q}}\].

• Cauchy-Schwarz Inequality is the special case \(p = q = 2\) \[\mathbb E|XY|\le (\mathbb E|X|^2)^{1/2}(\mathbb E|Y|^2)^{1/2}\]
  Also let \[u(t) = E[(tX - Y)^2]\] Then:\[t^2E[X^2] - 2tE[XY] + E[Y^2] \ge 0\]

This is a quadratic in \(t\). Thus the discriminant must be non-positive. Therefore: \[(\mathbb E[XY])^2 - \mathbb E[X^2]\mathbb E[Y^2] \le 0\] or \[(\mathbb E[XY])^2 \le \mathbb E[X^2]\mathbb E[Y^2]\] Then \[|(\mathbb E[XY])|\le(\mathbb E|X|^2)^{1/2}(\mathbb E|Y|^2)^{1/2}\]

• Covariance inequality \[(\text{Cov}[X, Y])^2 \leq \sigma_X^2 \sigma_Y^2=\text{Var}(X)\text{Var}(Y)\] If \(X\) and \(Y\) have means \(\mu_X\) and \(\mu_Y\) and variances \(\sigma_X^2\) and \(\sigma_Y^2\), respectively, we can apply the Cauchy-Schwarz Inequality to get \[\mathbb E\lvert(X-\mu_X)(Y-\mu_Y)\rvert\le \{\mathbb E(X-\mu_X)^2\}^{1/2}\{\mathbb E(Y-\mu_Y)^2\}^{1/2}\] Squaring both sides and using statistical notation, we have \[(\text{Cov}(X, Y))^2 \leq \sigma_X^2\sigma_Y^2.\]

• \[\mathbb E{}[\lvert{X}\lvert] \leq \mathbb E{}[\lvert{X}\lvert^p]^{1/p}\]

• Liapounov’s Inequality \[\mathbb E{}[\lvert{X}\lvert^r]^{1/r} \leq \mathbb E{}[\lvert{X}\lvert^s]^{1/s}\] where \(1 < r < s < \infty\).

• For numbers \(a_i,b_i,i=1,\cdots,n\), \[\sum_{i=1}^{n}|a_ib_i|\leq\left(\sum_{i=1}^{n}a_i^p\right)^{1/p}\left(\sum_{i=1}^{n}b_i^q\right)^{1/q}\] where \(\frac{1}{p}+\frac{1}{q}=1\)

• \[\frac{1}{n}\left(\sum_{i=1}^{n}|a_i|\right)^2\leq\sum_{i=1}^{n}a_i^2\]

Minkowski’s Inequality
Let \(X\) and \(Y\) be any two random variables. Then for \(1\le p<\infty\) \[(\mathbb E|X+Y|^p)^{1/p}\le (\mathbb E|X|^p)^{1/p}+(\mathbb E|Y|^p)^{1/p}\] Proof: Write \[\mathbb E(|X + Y|^p) = \mathbb E \left(|X + Y||X + Y|^{p-1}\right)\le \mathbb E\left(|X| |X + Y|^{p-1}\right) + \mathbb E\left(|Y||X + Y|^{p-1}\right),\] where we have used the fact that \(|X + Y| \le |X| + |Y|\) (the triangle inequality. Now apply HOlder’s Inequality to each expectation on the right-hand side to get \[\begin{align} \mathbb E(|X + Y|^p) &\leq \mathbb E\left(|X| |X + Y|^{p-1}\right) + \mathbb E\left(|Y||X + Y|^{p-1}\right)\\ &=\left(\mathbb E|X|^p\right)^{1/p} \left(\mathbb E|X + Y|^{q(p-1)}\right)^{1/q} + \left(\mathbb E|Y|^p\right)^{1/p}\left(\mathbb E|X + Y|^{q(p-1)}\right)^{1/q}\\ \end{align}\] where \(q\) satisfies \(1/p + 1/q = 1\). Now divide through by \(\left(\mathbb E|X + Y|^{q(p-1)}\right)^{1/q}\). Then \[\frac{\mathbb E(|X + Y|^p)}{\left(\mathbb E|X + Y|^{q(p-1)}\right)^{1/q}}=\left(\mathbb E(|X + Y|^p\right)^{1/p}\le \left(\mathbb E|X|^p\right)^{1/p} + \left(\mathbb E|Y|^p\right)^{1/p}\]

Functional Inequalities

Convex Function A function \(g(x)\) is convex on a set \(S\) if for all \(x, y \in S\) and \(0< \lambda < 1\) \[g(\lambda x + (1 - \lambda)y) \leq \lambda g(x) + (1 - \lambda)g(y).\] Strictly convex is when the inequality is strict. \(g\) is concave if \(-g\) is convex.

\(g(x)\) is convex on \(S\) if \(g''(x) \geq 0\) \(\forall x \in S\).

Jensen’s Inequality If \(g(x)\) is convex, then for any random variable \(X\) \[\mathbb E{}[g(X)] \leq g(\mathbb E{}[X]).\] Equality holds if and only if, for every line \(a + bx\) that is tangent to \(g(x)\) at \(x = \mathbb E{}[X]\), \(\text{P}{}\{g(X) = a + bX\} = 1\). (So if and only if \(g\) is affine with probability 1.)

• \(x^2\) is convex, \(\mathbb E[X^2] \geq \mathbb E{}[X]^2\)

• \(1/x\) is convex, \(\mathbb E[1/X] \geq 1/ \mathbb E{}[X]\)

• \(\log x\) is concave, \(\mathbb E(\log X) \leq \log (\mathbb E X)\)

• Arithmetic mean: \[a_A=\frac{1}{n}(a_1+a_2+\cdots+a_n)\] Geometric mean: \[a_G=(a_1a_2\cdots a_n)^{1/n}\] Harmonic mean: \[a_H=\frac{1}{\frac{1}{n}\left(\frac{1}{a_1}+\frac{1}{a_2}+\cdots+\frac{1}{a_n}\right)}\] \[a_H\leq a_G\leq a_A\]

Covariance Inequality Let \(X\) be any random variable and \(g(x)\) and \(h(x)\) any functions such that \(\mathbb E g(X), \mathbb E h(X)\) and \(\mathbb E (g(X)h(X))\) exist.
- If \(g(x)\) is a nondecreasing function and \(h(x)\) is a nonincreasing function, then \[\mathbb E \bigg(g(X)h(X)\bigg)\leq\bigg(\mathbb E g(X)\bigg)\bigg(\mathbb E h(X)\bigg)\]

• If \(g(x)\) and \(h(x)\) are both nondecreasing or nonincreasing functions, then \[\mathbb E \bigg(g(X)h(X)\bigg)\geq\bigg(\mathbb E g(X)\bigg)\bigg(\mathbb E h(X)\bigg)\]

Properties of the Sample Mean and Variance

The Derived Distributions: Student’s \(t\) and Snedecor’s \(F\)

Convergence in Probability

Proof is using Chebychev’s inequality.

We have for every \(\epsilon>0\), \[P(|\bar X_n - \mu|\ge\epsilon)=P((\bar X_n - \mu)^2\ge\epsilon^2)\leq\frac{\mathbb E(\bar X_n - \mu)^2}{\epsilon^2}=\frac{\text{Var}\bar X_n}{\epsilon^2}=\frac{\sigma^2}{n\epsilon^2}\] Hence, \[P(|\bar X_n - \mu|<\epsilon)=1-P(|\bar X_n - \mu|\ge\epsilon)\ge 1-\frac{\sigma^2}{n\epsilon^2}\to 1 \text{ as } n\to \infty\]

The property summarized by the WLLN, that a sequence of the “same” sample quantity approaches a constant as \(n \to \infty\), is known as consistency.

Almost Sure Convergence

Convergence in Distribution

The Delta Method

If we are interested in the convergence of some function of a sequence of RVs, rather than the RVs themselves, then we can use the Delta Method (follows from an application of Taylor’s theorem and Slutsky’s theorem).

If \(g'(\theta) = 0\) then we take the next term in the Taylor series.

Proof The Taylor expansion of \(g(Y_n)\) around \(Y_n = 0\) is \[g(Y_n) = g(\theta) + g'(\theta) (Y_n - \theta) + \text{Remainder},\] where the remainder \(\to 0\) as \(Y_n \to \theta\). Since \(Y_n \to \theta\) in probability it follows that the remainder \(\to 0\) in probability. By applying Slutsky’s Theorem to \[\sqrt{n}[g(Y_n) - g(\theta)] = g'(\theta)\sqrt{n}(Y_n-\theta),\] the result now follows.

The second-order Delta Method

If \(g'(\theta) = 0\), we take one more term in the Taylor expansion to get \[g(Y_n) = g(\theta) + g'(\theta)(Y_n -\theta) +\frac{g''(\theta)}{2} (Y_n -\theta)^2 + \text{Remainder}.\] If we do some rearranging (setting \(g' = 0\)), we have \[ g(Y_n) -g(\theta) = \frac{g''(\theta)}{2} (Y_n -\theta)^2 + \text{Remainder}.\] Now recall that the square of a \(n(0, 1)\) is a \(\chi_1^2\) , which implies that \[\frac{n(Y_n-\theta)^2}{\sigma^2}\to\chi_1^2\] in distribution. Therefore,

(Second-order Delta Method) Let \(Y_n\) be a sequence of random variables that satisfies \(\sqrt{n}(Y_n -\theta) \to n(0, \sigma^2)\) in distribution. For a given function \(g\) and a specific value of \(\theta\), suppose that \(g'(\theta)=0\) and \(g''(\theta)\) exists and is not \(0\). Then \[n[g(Y_n) -g(\theta)]\approx n\frac{g''(\theta)}{2} (Y_n -\theta)^2\to \sigma^2\frac{g''(\theta)}{2}\chi_1^2\] in distribution.

Multivariate Delta Method

Note that we must deal with multiple random variables although the ultimate CLT is a univariate one. Suppose the vector-valued random variable \(\mathbf X=(X_1,\cdots, X_p)\) has mean \(\boldsymbol\mu = (\mu_1, \cdots , \mu_p)\) and covariances \(Cov(X_i, X_j) = \sigma_{ij}\), and we observe an independent random sample \(\mathbf X_1, \cdots, \mathbf X_n\) and calculate the means \[\bar X_i = \frac{1}{n}\sum_{k=1}^{n}X_{ik},\quad i=1,\cdots,p\] For a function \(g(\mathbf x) = g(x_1, \cdots, x_p)\) we can write \[g(\bar x_1, \cdots, \bar x_p) = g(\mu_1, \cdots, \mu_p) + \sum_{k=1}^{p}g'_k(\mathbf x)(\bar x_k - \mu_k),\] and we then have the following theorem:

Multivariate Delta Method Let \(\mathbf X_1, \cdots, \mathbf X_n\) be a random sample with \(E(X_{ij})=\mu_i\) and \(Cov(X_{ik}, X_{jk}) = \sigma_{ij}\) . For a given function \(g\) with continuous first partial derivatives and a specific value of \(\boldsymbol\mu = (\mu_1, \cdots , \mu_p)\) for which \[\tau^2=\sum\sum\sigma_{ij}\frac{\partial g(\boldsymbol\mu)}{\partial{\mu_i}}\frac{\partial g(\boldsymbol\mu)}{\partial{\mu_j}}>0,\] \[\sqrt{n}\big[g(\bar x_1,\cdots,\bar x_p)-g(\mu_1,\cdots,\mu_p)\big]\to n(0,\tau^2)\;\;\text {in distribution}\]

Direct Method

Probability integral transformation Let \(X\) have continuous cdf \(F_X(x)\) and define the random variable \(Y\) as \(Y = F_X (X)\). Then \(Y\) is uniformly distributed on \((0, 1)\), that is, \(P(Y \leq y) = y, 0 < y < 1\).

Proof For \(Y = F_X (X)\) we have, for \(0 < y < 1\), \[\begin{align} P(Y\leq y)&=P(F_X (X)\leq y)\\ &=P(F_X^{-1} [F_X (X)]\leq F_X^{-1}(y))\quad (F_X^{-1} \text{ is increasing})\\ &=P(X\leq F_X^{-1}(y))\\ &=F_X(F_X^{-1}(y))\\ &=y \end{align}\] At the endpoints we have \(P (Y \leq y) = 1\) for \(y \ge 1\) and \(P (Y \leq y)=0\) for \(y \leq 0\), showing that \(Y\) has a uniform distribution.

Let \(X\) have cdf \(F_X(x)\), let \(Y=g(X)\), and let \(\mathcal X\) and \(\mathcal Y\) be defined as \(\mathcal X \{x: f_X(x) > 0\}\) and \(\mathcal Y = \{y: y = g(x)\},\text{ for some }x \in \mathcal X\}\).

• 1. If \(g\) is an increasing function on \(\mathcal X\), \[F_Y (y) =\int_{x\in \mathcal X:x\leq g^{-1}(y)}f_X(x)dx=\int_{-\infty}^{g^{-1}(y)}f_X(x)dx= F_X (g^{-1}(y)) \text{ for }y \in \mathcal Y\]
• 2. If \(g\) is a decreasing function on \(\mathcal X\) and \(X\) is a continuous random variable, \[F_Y (y)=\int_{g^{-1}(y)}^{\infty}f_X(x)dx=1 - F_X(g^{-1} (y)) \text{ for } y \in \mathcal Y\]

Indirect Methods

The Accept/Reject Algorithm

• It is typical to call \(V\) the candidate density and \(Y\) the target density.

• One would normally try to choose a candidate density with heavier tails than the target density (e.g. Cauchy and normal) to ensure that the tails of the target are well represented. If the target has heavy tails, however, it can be hard to find a candidate that results in finite \(M\). In this case people turn to MCMC methods.

• Note that \(\text{P }(\texttt{terminate}) = 1/M\). The number of trials to generate one sample of \(Y\) is therefore \(\text{geometric}(1/M)\), with \(M\) the expected number of trials.

• The intuition behind this algorithm is that if we consider placing the density of a random variable \(Y\) in a box (2d for simplicity) with coordinates \((v,u)\), we express the cdf of \(Y\) using \(V, U \sim \text{uniform}(0, 1)\) \[\text{P }(Y \leq y) = \text{P }(V \leq y \vert{} U \leq \frac1c f_Y(V))\] where \(c = \sup_y f_Y(y)\). In the actual algorithm we take \(U \sim \text{uniform}(0,1)\) and \(V\) to be an RV that has common support with \(Y\).

Definition: sample percentile

The notation \(\{b\}\), when appearing in a subscript, is defined to be the number \(b\) rounded to the nearest integer in the usual way. More precisely, if \(i\) is an integer and \(i - .5 \le b < i + .5\), then \(\{b\} = i\).

The \((100p)th\) sample percentile is \(X_{(\{np\})}\) if \(\frac{1}{2n}<p<.5\) and \(X_{(n+1-\{n(1-p)\})}\) if \(.5<p<1-\frac{1}{2n}\)

Factorization Theorem

Let \(f(\mathbf{x} \vert{} \theta)\) denote the pmf/pdf of a sample \(\mathbf{X}\). A statistic \(T(\mathbf{X})\) is a sufficient statistic for \(\theta\) if and only if there exists functions \(g(T(\mathbf{X})\vert{} \theta)\) and \(h(\mathbf{x})\) such that \(\forall \mathbf{x} \in \mathcal X{}\), \(\forall \theta \in \Theta\) \[\begin{equation} f(\mathbf{x} \vert{} \theta) = g(T(\mathbf{x}) \vert{} \theta)h(\mathbf{x}). \end{equation}\]

This theorem shows that the identity is a sufficient statistic. It is straightforward to show from this that any bijection of a sufficient statistic is a sufficient statistic.

Proof Suppose \(T(\mathbf X)\) is a sufficient statistic. Choose \(g(t\lvert \theta) = P_{\theta}(T(\mathbf X) = t)\) and \(h(\mathbf x) =P\big(\mathbf X = \mathbf x\lvert T(\mathbf X) = T(\mathbf x)\big)\). Because \(T(\mathbf X)\) is sufficient, the conditional probability defining \(h(\mathbf x)\) does not depend on \(\theta\). Thus this choice of \(h(\mathbf x)\) and \(g(t\lvert \theta)\) is legitimate, and for this choice we have \[\begin{align} f(\mathbf x\lvert \theta) &= P_{\theta}(\mathbf X = \mathbf x)\\ &= P_{\theta}\bigg(\mathbf X = \mathbf x \text{ and } T(\mathbf X) = T(\mathbf x)\bigg)\\ &= P_{\theta}(T(\mathbf X) = T(\mathbf x)) P(\mathbf X=\mathbf x\lvert T(\mathbf X) = T(\mathbf x))\quad \text{(sufficiency)}\\ &= g(T(\mathbf x)\lvert \theta)h(\mathbf x). \end{align}\]

So factorization has been exhibited. We also see from the last two lines above that \[P_{\theta}\big(T(\mathbf X) = T(\mathbf x)\big) =g(T(\mathbf x)\lvert \theta),\] so \(g(T(\mathbf x)\lvert \theta)\) is the pmf of \(T(\mathbf X).\)

Now assume the factorization exists. Let \(q(t\lvert \theta)\) be the pmf of \(T(\mathbf X)\). To show that \(T(\mathbf X)\) is sufficient we examine the ratio \(f(\mathbf x\lvert \theta)/q(T(\mathbf x)\lvert \theta)\). Define \(A_{T(\mathbf x)} =\{\mathbf y: T(\mathbf y) = T(\mathbf x)\}\). Then \[\begin{align} \frac{f(\mathbf x\lvert \theta)}{q(T(\mathbf x)\lvert \theta)}&=\frac{g(T(\mathbf x)\lvert \theta)h(\mathbf x)}{q(T(\mathbf x)\lvert \theta)}\\ &=\frac{g(T(\mathbf x)\lvert \theta)h(\mathbf x)}{\sum_{A_{T(\mathbf x)}}g(T(\mathbf y)\lvert \theta)h(\mathbf y)}\quad \text{(definition of the pmf of T)}\\ &=\frac{g(T(\mathbf x)\lvert \theta)h(\mathbf x)}{g(T(\mathbf y)\lvert \theta)\sum_{A_{T(\mathbf x)}}h(\mathbf y)}\quad \text{(since T is constant on} \;A_{T(\mathbf x)})\\ &=\frac{g(T(\mathbf x)\lvert \theta)h(\mathbf x)}{g(T(\mathbf x)\lvert \theta)\sum_{A_{T(\mathbf x)}}h(\mathbf y)}\\ &=\frac{h(\mathbf x)}{\sum_{A_{T(\mathbf x)}}h(\mathbf y)}\\ \end{align}\]

Since the ratio does not depend on \(\theta\), \(T(\mathbf X)\) is a sufficient statistic for \(\theta\).

Let \(X_1, \dots, X_n\) be iid observations from a pmf/pdf \(f(x \vert{} \boldsymbol\theta)\) from an exponential family \[f(x|\boldsymbol{\theta}) = h(x)c(\boldsymbol{\theta})\exp\left( \sum_{i=1}^k w_i(\boldsymbol{\theta}) t_i(x) \right),\] where \(\boldsymbol{\theta} = (\theta_1, \dots, \theta_d)\), \(d \leq k\). Then \(T(\mathbf{X})\) defined by \[T(\mathbf{X}) = \left(\sum_{j=1}^n t_1(\mathbf{X}_j),\cdots,\sum_{j=1}^n t_k(\mathbf{X}_j)\right)\] is a sufficient statistic for \(\boldsymbol{\theta}\).

Minimal sufficient statistic A sufficient statistic \(T(\mathbf{X})\) is called a minimal sufficient statistic if, for any other sufficient statistic \(T'(\mathbf{X})\), \(T(\mathbf{x})\) is a function of \(T'(\mathbf{x})\).

By ‘function of’ we mean that if \(T'(\mathbf{x}) = T'(\mathbf{y})\) then \(T(\mathbf{x}) = T(\mathbf{y})\), \(T\) varies with respect to \(X\) only as it varies with \(T'\). This means that each tile of the partition \(\{B_{t'}:t'\in \mathcal T'\}\) of the sample space according to the image of \(T'\) is a subset of some tile in the partition \(\{A_{t}:t\in \mathcal T\}\) according to \(T\). This means that minimal sufficient statistics provide the coarsest possible tiling of the sample space and thus are the sufficient statistics that provide the greatest data reduction.

Theorem Let \(f(\mathbf{x} \vert{} \theta)\) be the pmf/pdf of a sample \(\mathbf{X}\). Suppose there exists a function \(T(\mathbf{X})\) such that \(\forall \mathbf{x}, \mathbf{y} \in \mathcal X{}\) the ratio \(f(\mathbf{x}\vert{}\theta)/f(\mathbf{y}\vert{}\theta)\) is constant as a function of \(\theta\) if an only if \(T(\mathbf{x}) = T(\mathbf{y})\). Then \(T(\mathbf{X})\) is a minimal sufficient statistic for \(\theta\).

Proof To simplify the proof, we assume \(f(\mathbf x\lvert \theta) > 0\) for all \(\mathbf x \in \mathcal X\) and \(\theta\).

First we show that \(T(\mathbf X)\) is a sufficient statistic. Let \(\mathcal T = \{t : t=T(\mathbf x)\; \text{for some }\; \mathbf x \in \mathcal X\}\) be the image of \(\mathcal X\) under \(T(\mathbf x)\). Define the partition sets induced by \(T(\mathbf x)\) as \(A_t = \{\mathbf x: T(\mathbf x) = t\}\). For each \(A_t\), choose and fix one element \(\mathbf x_t \in A_t\). For any \(\mathbf x \in \mathcal X\), \(\mathbf x_{T(\mathbf x)}\) is the fixed element that is in the same set, \(A_t\) , as \(\mathbf x\). Since \(\mathbf x\) and \(\mathbf x_{T(\mathbf x)}\) are in the same set \(A_t\) , \(T(\mathbf x) = T(\mathbf x_{T(\mathbf x)})\) and, hence, \(f(\mathbf x\lvert\theta)/ f(\mathbf x_{T(\mathbf x)}\lvert\theta)\) is constant as a function of \(\theta\). Thus, we can define a function on \(\mathcal X\) by \(h(\mathbf x) = f(\mathbf x\lvert\theta)/ f(\mathbf x_{T(\mathbf x)}\lvert\theta)\) and \(h\) does not depend on \(\theta\). Define a function on \(\mathcal T\) by \(g(t\lvert\theta) = f(\mathbf x_t\lvert\theta)\). Then it can be seen that \[f(\mathbf x_t\lvert\theta)=\frac{f(\mathbf x_{T(\mathbf x)}\lvert\theta)f(\mathbf x\lvert\theta)}{f(\mathbf x_{T(\mathbf x)}\lvert\theta)} = g(T(\mathbf x)\lvert\theta)h(\mathbf x)\]
and, by the Factorization Theorem, \(T(\mathbf X)\) is a sufficient statistic for \(\theta\).

Now to show that \(T(\mathbf X)\) is minimal, let \(T'(\mathbf X)\) be any other sufficient statistic. By the Factorization Theorem, there exist functions \(g'\) and \(h'\) such that \(f(\mathbf x\lvert\theta) =g'(T'(\mathbf x)\lvert\theta)h'(\mathbf x)\). Let \(\mathbf x\) and \(\mathbf y\) be any two sample points with \(T'(\mathbf x) = T'(\mathbf y)\). Then \[\frac{f(\mathbf x\lvert\theta)}{f(\mathbf y\lvert\theta)}=\frac{g'(T'(\mathbf x)\lvert\theta)h'(\mathbf x)}{g'(T'(\mathbf y)\lvert\theta)h'(\mathbf y)}=\frac{h'(\mathbf x)}{h'(\mathbf y)}\] Since this ratio does not depend on \(\theta\), the assumptions of the theorem imply that \(T(\mathbf x) = T(\mathbf y)\). Thus, \(T(\mathbf x)\) is a function of \(T'(\mathbf x)\) and \(T(\mathbf x)\) is minimal.

Necessary Statistic A statistic is necessary if it can be written as a function of every sufficient statistic.

A statistic is minimal sufficient if and only if it is a necessary and sufficient statistic.

Method of Moments

The method of moments is performed using a random sample \(X_1, \dots, X_n\) from a population with unknown parameters \(\theta_1, \dots, \theta_k\) by computing the first \(k\) empirical moments \(m_k = \frac{1}{n}\sum_{i=1}^nX_i^k\) and matching them with the corresponding moments of the population \(\mu_k' = \mathbb E[X^t]\).

\[\begin{aligned} m_1 &= \mu_1'(\theta_1, \dots, \theta_k)\\ m_2 &= \mu_2'(\theta_1, \dots, \theta_k)\\ &\vdots \\ m_k &= \mu_k' (\theta_1, \dots, \theta_k)\\ \end{aligned}\]

From this you get \(k\) simultaneous equations that you can use to solve for the parameters of the population.

Maximum Likelihood Estimators

Maximum Likelihood Estimator For each sample point \(\mathbf x\), let \(\hat{\boldsymbol\theta}(\mathbf x)\) be a parameter value at which the likelihood function \[L(\boldsymbol\theta \vert{} \mathbf x)=\prod_{i=1}^{n}f(x_i\lvert \theta_1,\cdots,\theta_k)\] attains its maximum as a function of \(\boldsymbol\theta\), with \(\mathbf x\) held fixed. A maximum likelihood estimator (MLE) of the parameter \(\boldsymbol\theta\) based on a sample \(\mathbf X\) is \(\hat{\boldsymbol\theta}(\mathbf X)\).

Suppose we want to find the MLE for some function of the parameter \(\tau(\theta)\).

Induced Likelihood Given some function of the parameter \(\tau(\theta)\), we define the induced likelihood function \(L^*\) by \[L^*(\eta \vert{} \mathbf x) = \sup_{\theta : \tau(\theta) = \eta} L(\theta \vert{} \mathbf x).\] The value \(\hat{\eta}\) that maximises \(L^*(\eta \vert{} \mathbf x)\) will be called the MLE of \(\eta = \tau(\theta)\). It can be seen that the maxima of \(L^*\) and \(L\) coincide.

Invariance Property of MLEs IF \(\hat{\theta}\) is the MLE of \(\theta\), then for any function \(\tau(\theta)\), the MLE of \(\tau(\theta)\) is \(\tau(\hat{\theta})\).

1. The MLE can be an unstable function of the data.

2. When verifying maxima for multi-dimensional problems try to avoid going down the hessian route, which can be tedious.

To use two-variate calculus to verify that a function \(H(\theta_1,\theta_2)\) has a local maximum at \((\hat\theta_1,\hat\theta_2)\), it must be shown that the following three conditions hold.

• 1. The first-order partial derivatives are \(0\), \[\frac{\partial}{\partial\theta_1}H(\theta_1,\theta_2)|_{\theta_1=\hat\theta_1,\theta_2=\hat\theta_2}=0\] and \[\frac{\partial}{\partial\theta_2}H(\theta_1,\theta_2)|_{\theta_1=\hat\theta_1,\theta_2=\hat\theta_2}=0\]
• 2. At least one second-order partial derivative is negative, \[\frac{\partial^2}{\partial\theta_1^2}H(\theta_1,\theta_2)|_{\theta_1=\hat\theta_1,\theta_2=\hat\theta_2}<0\] or \[\frac{\partial^2}{\partial\theta_2^2}H(\theta_1,\theta_2)|_{\theta_1=\hat\theta_1,\theta_2=\hat\theta_2}<0\]
• 3. The Jacobian of the second-order partial derivatives is positive, \[\begin{vmatrix} \frac{\partial^2}{\partial\theta_1^2}H(\theta_1,\theta_2)&\frac{\partial^2}{\partial\theta_1\partial\theta_2}H(\theta_1,\theta_2)\\ \frac{\partial^2}{\partial\theta_1\partial\theta_2}H(\theta_1,\theta_2)&\frac{\partial^2}{\partial\theta_2^2}H(\theta_1,\theta_2) \end{vmatrix}_{\theta_1=\hat\theta_1,\theta_2=\hat\theta_2}\\ =\frac{\partial^2}{\partial\theta_1^2}H(\theta_1,\theta_2)\frac{\partial^2}{\partial\theta_2^2}H(\theta_1,\theta_2)-\left(\frac{\partial^2}{\partial\theta_1\partial\theta_2}H(\theta_1,\theta_2)\right)^2 >0\]
     Since \[f(\theta) \approx f(\theta_0) + \nabla f(\theta_0)^T (\theta - \theta_0) + \frac12 (\theta - \theta_0)^T \nabla^2 f(\theta_0) (\theta - \theta_0).\] If \(\theta_0\) is a local extremum for \(f\), then \(\nabla f(\theta_0) = 0\), so if \(\nabla^2 f(\theta_0)\) is negative definite, then \((\theta - \theta_0)^T \nabla^2 f(\theta_0) (\theta - \theta_0) \leq 0\) for all \(\theta\), so \(\theta_0\) is a local maximum. Then the eigenvalues of \(\nabla^2 f(\theta_0)\) should be both negative and the determinant should be positive.

Bayes Estimators

If we denote the prior distribution by \(\pi(\theta)\) and the sampling distribution by \(f(\mathbf x \vert{} \theta)\), then the posterior distribution of \(\theta\) given the sample \(\mathbf x\) is \[\pi(\theta \vert{} \mathbf x) = \frac{f(\mathbf x\vert{}\theta)\pi(\theta)}{m(\mathbf x)},\] where \(m(\mathbf x)\) is the marginal distribution of the sample \(\mathbf x\) \[m(\mathbf x) = \int f(\mathbf x\vert{}\theta) \pi(\theta) d\theta.\]

Conjugate Priors Let \(\mathcal{F}\) denote the class of pdfs/pmfs \(f(\mathbf x \vert{} \theta)\) (indexed by \(\theta\)). A class \(\Pi\) of prior distribution is a conjugate family for \(\mathcal{F}\) if, \(\forall f \in \mathcal{F}\), \(\forall \, \text{priors} \in \mathcal{F}\) and \(\forall \mathbf x \in \mathcal X\), the posterior distribution is in \(\Pi\).

Note that this relation is not said to be symmetric.

[Some Examples]

• Normal is self-conjugate as a family.

• Beta distribution is conjugate to binomial.

• Gamma distribution is conjugate to Poisson.
  Let \(X_1,\cdots, X_n\) be iid \(Poisson(\lambda)\), and let \(\lambda\) have a \(gamma(\alpha,\beta)\) distribution, the conjugate family for the Poisson.

For \(n\) observations, \(Y =\sum_i X_i\sim Poisson(n\lambda)\). - a. The marginal pmf of \(Y\) is \[\begin{align} m(y)&=\int_{-\infty}^{\infty}f(Y|\lambda)f(\lambda)d\lambda\\ &=\int_{-\infty}^{\infty}\frac{(n\lambda)^ye^{-n\lambda}}{y!}\frac{1}{\Gamma(\alpha)\beta^{\alpha}}\lambda^{\alpha-1}e^{-\lambda/\beta}d\lambda\\ &=\frac{n^y}{y!\Gamma(\alpha)\beta^{\alpha}}\int_{-\infty}^{\infty}\lambda^{(y+\alpha)-1}e^{-n\lambda-\lambda/\beta}d\lambda\\ &=\frac{n^y}{y!\Gamma(\alpha)\beta^{\alpha}}\int_{-\infty}^{\infty}\lambda^{(y+\alpha)-1}e^{-\frac{\lambda}{\beta/(n\beta+1)}}d\lambda\\ &=\frac{n^y}{y!\Gamma(\alpha)\beta^{\alpha}}\Gamma(y+\alpha)\left(\frac{\beta}{n\beta+1}\right)^{y+\alpha}\\ \end{align}\] Thus \[\pi(\lambda|y)=\frac{f(y|\lambda)\pi(\lambda)}{m(y)}=\frac{\lambda^{(y+\alpha)-1}e^{-\frac{\lambda}{\beta/(n\beta+1)}}}{\Gamma(y+\alpha)\left(\frac{\beta}{n\beta+1}\right)^{y+\alpha}}\\ \sim gamma\left(y+\alpha, \frac{\beta}{n\beta+1}\right)\]

The EM (Expectation-Maximization) Algorithm

In general, we can move in either direction, from the complete-data problem to the incomplete-data problem or the reverse. If \(\mathbf Y = (Y_1, \cdots , Y_n)\) are the incomplete data, and \(\mathbf X = (X_1, \cdots, X_m)\) are the augmented data, making \((\mathbf Y, \mathbf X)\) the complete data, the densities \(g(\cdot\lvert \theta)\) of \(\mathbf Y\) and \(f(\cdot\lvert\theta)\) of \((\mathbf Y, \mathbf X)\) have the relationship \[g(\mathbf y\lvert\theta) = \int f(\mathbf y, \mathbf x\lvert\theta) d\mathbf x\] with sums replacing integrals in the discrete case.
If we turn these into likelihoods, \[L(\theta\lvert\mathbf y) = g(\mathbf y\lvert\theta)\] is the incomplete-data likelihood and \[L(\theta\lvert\mathbf y, \mathbf x)=f(\mathbf y,\mathbf x\lvert\theta)\] is the complete-data likelihood.

Mean Squared Error

The mean squared error (MSE) of an estimator \(W\) of a parameter \(\theta\) is defined by \(\mathbb E_{\theta}[(W - \theta)^2]\).

Bias-Variance Decomposition \[\mathbb E_{\theta}[(W - \theta)^2] = \text{Var}_{\theta} [W] + (\mathbb E_{\theta} [W] - \theta)^2 = \text{Var}_{\theta} [W] + (\text{Bias}_{\theta} [W])^2\]

Bias The bias of a point estimator \(W\) of a parameter \(\theta\) is given by \[\text{Bias}_{\theta} [W] = \mathbb E_{\theta} [W] - \theta.\] An estimator whose bias is \(0\) is called an unbiased estimator and has \(\mathbb E_{\theta}[W] = \theta,\quad\forall \theta\).

• Clearly, if an estimator is unbiased, then its MSE is equal to its variance.

• The MSE makes sense for location parameters but not so much for scale parameters, since it is symmetric and scale parameters have a natural floor at \(0\).

• The MSE may be a function of the thing you’re trying to estimate. So which estimator you choose as being the ‘best’ may depend on the range you expect the parameter to lie within.

The MLE (biased) estimator of the variance of a normal distribution is \[\hat \sigma^2 = \frac{1}{N} \sum_{i=1}^{N}(X_i - \bar X)^2,\] where \(\bar X\) is the sample mean and \(X_i \sim^{iid} \mathcal{N}(\mu,\sigma^2)\). Since \[\frac{1}{\sigma ^2} \sum_{i=1}^{N}(X_i-\bar X)^2 \sim \chi^2_{N},\] we have that \[\begin{align} \text{Var}\big(\frac{N\hat \sigma^2}{\sigma^2}\big) = \text{Var}(\chi^2_{N-1}) = 2(N-1) \\ \implies \text{Var}(\hat \sigma ^2) = \frac{2(N-1) \sigma^4}{N^2} \end{align} \] The variance of this estimator is equal to \(\frac{2(N-1) \sigma^4}{N^2}\). The MSE is \[\text{MSE}=\text{Var}(\hat \sigma ^2)+\text{Bias} =\frac{2\sigma ^4}{N}+ \mathbb E [\hat \sigma ^2] - \sigma ^2=\text{Var}(\hat \sigma ^2)+\text{Bias}^2 \\=\text{Var}(\hat \sigma ^2)+(\mathbb E(\hat \sigma ^2)-\sigma ^2)^2=\frac{2(N-1) \sigma^4}{N^2}+ \left(\frac{N-1}{N}\sigma ^2 - \sigma ^2\right)^2\\ =\frac{2N-1}{N^2}\sigma ^4.\] where \[\begin{align} \mathbb E(\hat \sigma ^2)&=\mathbb E\left(\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})^2\right)\\ &=\frac{1}{N}\mathbb E\left(\sum_{i=1}^{N}x_i^2-2\sum_{i=1}^{N}x_i\bar{x}+\sum_{i=1}^{N}\bar{x}^2\right)\\ &=\frac{1}{N}\mathbb E\left(\sum_{i=1}^{N}x_i^2-N\bar{x}^2\right)\\ &=\frac{1}{N}\mathbb E\left(\sum_{i=1}^{N}x_i^2\right)-\mathbb E\left(\bar{x}^2\right)\\ &=\mathbb E\left(x^2\right)-\mathbb E\left(\bar{x}^2\right)\\ &=(\sigma_x^2+\mu^2)-(\sigma_{\bar x}^2+\mu^2)\\ &=\sigma_x^2-\sigma_{\bar x}^2\\ &=\sigma_x^2-\text{Var}(\bar x)\\ &=\sigma_x^2-\text{Var}\left(\frac{1}{N}\sum_{i=1}^{N}x_i\right)\\ &=\sigma_x^2-\frac{1}{N^2}\text{Var}\left(\sum_{i=1}^{N}x_i\right)\\ &=\sigma_x^2-\frac{1}{N^2}\sum_{i=1}^{N}\text{Var}\left(x\right)\\ &=\sigma_x^2-\frac{1}{N}\text{Var}\left(x\right)\\ &=\sigma_x^2-\frac{1}{N}\sigma_x^2\\ &=\frac{N-1}{N}\sigma_x^2\\ \end{align}\]

A unbiased estimator of the variance of a normal distribution is \[S^2 = \frac{1}{N-1} \sum_{i=1}^{N}(X_i - \bar X)^2,\] since \(\mathbb E S^2=\sigma^2\). And since \[\frac{(N-1) S^2}{\sigma^2} \sim \chi_{n-1}^2,\] \[\begin{align} \text{Var}\big(\frac{(N-1)S^2}{\sigma^2}\big) = \text{Var}(\chi^2_{N-1}) = 2(N-1) \\ \implies \text{Var}(S^2) = \frac{2(N-1) \sigma^4}{(N-1)^2}=\frac{2\sigma^4}{N-1} \end{align} \] The variance of this estimator is equal to \(\frac{2\sigma^4}{N-1}\), which is also the MSE.

The variance (or MSE) in the unbiased case \(\frac{2\sigma^4}{N-1}\) is greater than the variance in the biased case \(\frac{2(N-1) \sigma^4}{N^2}\). The MSE of the unbiased estimator \(\frac{2\sigma^4}{N-1}\) is larger than the MSE of the biased estimator \(\frac{2N-1}{N^2}\sigma ^4\).

\[ \begin{array}{llll} \text{Estimator} & \text{MLE} & \text{best-unbiased} & \text{Bayes} \\ \hline \text{Normal}(\mu) & \bar X & \bar X & \frac{\tau^2}{\tau^2+\sigma^2}\bar X+\frac{\sigma^2}{\tau^2+\sigma^2}\mu \\ \text{Normal}(\sigma^2) & \frac{1}{n}\sum(X_i-\bar X)=\frac{n-1}{n}S^2 & S^2=\frac{1}{n-1}\sum(X_i-\bar X) & \frac{\sigma^2\tau^2}{\sigma^2+\tau^2} \\ \text{Bernoulli}(p) & \bar X & & \\ \text{Bernoulli}(p(1-p)) & & & \\ \text{Binomial}(np) & \frac{\sum_{i=1}^n X_i}{n} & & \frac{(\sum_{i=1}^n X_i) +\alpha}{\alpha+\beta+n} \\ \text{Binomial}(np(1-p)) & & & \\ \end{array} \]

Best Unbiased Estimators

[Best Unbiased Estimator] An estimator \(W^*\) is best unbiased estimator of of \(\tau(\theta)\) if it satisfies \(\mathbb E[W^*] = \tau(\theta) \,\, \forall \theta\), and for any other estimator with \(\mathbb E[W] = \tau(\theta) \,\, \forall \theta\) we have \(\text{Var}[W*] \leq \text{Var}[W] \,\, \forall \theta\). We also call \(W^*\) the uniform minimum variance unbiased estimator (UMVUE) of \(\tau(\theta)\).

Suppose that we are trying to estimate \(\theta\) and consider the class of estimators \[\mathcal{C}_\tau = \{W: \mathbb E[W] = \tau(\theta)\}.\] All estimators in this class have the same bias, so we can compare their MSEs by comparing their variances alone. (So the best estimator in this class is just the minimum variance one.) This means that the considerations of this chapter can be applied to classes like \(\mathcal{C}_\tau\), even if \(\tau(\theta) \neq \theta\).

The best unbiased estimator, if it exists, could be hard to find. The following lower bound at least gives a stopping criterion to our search.

Fisher Information

Fisher Information The quantity \[\mathbb E_{\theta}\left[ \left(\frac{\partial}{\partial \theta} \log f(\mathbf X \vert{} \theta) \right)^2 \right]\] is called the Fisher Information or Information number of the sample \(\mathbf X\). This terminology reflects the fact that the information number gives a bound on the variance of the best unbiased estimator of \(\theta\). As the information number gets bigger and we have more information about \(\theta\), we have a smaller bound on the variance of the best unbiased estimator.

If \(f(x \vert{} \theta)\) satisfies \[\frac{d}{d \theta} \mathbb E_{\theta}\left[ \frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right] = \int \frac{\partial}{\partial \theta} \left[ \left( \frac{\partial}{\partial \theta} \log f(x \vert{} \theta) \right) f(x \vert{} \theta) \right] dx\] (as is true for an exponential family), the partial derivative with respect to \(\theta\) of the natural logarithm of the likelihood function is called the score. Under certain regularity conditions, if \(\theta\) is the true parameter (i.e. \(X\) is actually distributed as \(f(X; \theta)\)), it can be shown that the expected value (the first moment) of the score, evaluated at the true parameter value \(\theta\) , is \(0\)

\[\begin{align} \mathbb E_{\theta}\left[\frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right]&=\int_{\mathbb R}^{}\left[\frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right]f(X \vert{} \theta) dx\\ &=\int_{\mathbb R}^{}\frac{\frac{\partial}{\partial \theta} f(X \vert{} \theta)}{f(X \vert{} \theta)} f(X \vert{} \theta) dx\\ &=\frac{\partial}{\partial \theta}\int_{\mathbb R}^{} f(X \vert{} \theta) dx\\ &=\frac{\partial}{\partial \theta}\cdot 1\\ &=0\\ \end{align}\] The variance of the score is defined to be the Fisher information: \[\mathcal I_{\theta}=\mathbb E_{\theta}\left[ \left(\frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right)^2 \right]=\int_{\mathbb R}\left(\frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right)^2f(X \vert{} \theta)dx\] Since \[\begin{align} \frac{\partial^2}{\partial \theta^2} \log f(X \vert{} \theta)&=\frac{\partial}{\partial \theta}\frac{\frac{\partial}{\partial \theta} f(X \vert{} \theta)}{f(X \vert{} \theta)}\\ &=\frac{\left[\frac{\partial^2}{\partial \theta^2} f(X \vert{} \theta)\right]f(X \vert{} \theta)-\left[\frac{\partial}{\partial \theta} f(X \vert{} \theta)\right]^2}{(f(X \vert{} \theta))^2}\\ &=\frac{\frac{\partial^2}{\partial \theta^2} f(X \vert{} \theta)}{f(X \vert{} \theta)}-\frac{\left[\frac{\partial}{\partial \theta} f(X \vert{} \theta)\right]^2}{(f(X \vert{} \theta))^2}\\ &=\frac{\frac{\partial^2}{\partial \theta^2} f(X \vert{} \theta)}{f(X \vert{} \theta)}-\left(\frac{\partial}{\partial \theta} \log f(X \vert{} \theta)\right)^2\\ \end{align}\] and \[\mathbb E_{\theta}\left[\frac{\frac{\partial^2}{\partial \theta^2} f(X \vert{} \theta)}{f(X \vert{} \theta)}\right]=\frac{\partial^2}{\partial \theta^2}\int_{\mathbb R}f(X \vert{} \theta)dx=0\] then the Fisher information can be written \[\begin{equation} \mathbb E_{\theta}\left[ \left(\frac{\partial}{\partial \theta} \log f(X \vert{} \theta) \right)^2 \right] =\mathbb E_{\theta}\left[ \left(\frac{\frac{\partial}{\partial \theta} f(X \vert{} \theta)}{f(X \vert{} \theta)}\right)^2 \right]= - \mathbb E_{\theta}\left[ \frac{\partial^2}{\partial \theta^2} \log f(X \vert{} \theta) \right] \tag{1} \end{equation}\]

Even if the Cramér-Rao bound is applicable, it may not be sharp – there may not be an estimator that attains this bound.

• For Bernoulli RVs, the Fisher information is \[\begin{align} I_X(p)&=E_p\left[\left(\frac{d}{dp}\log\left(p^X(1-p)^{1-X}\right)\right)^2\right]\\ &=E_p\left[\left(\frac{d}{dp}\left[X\log p+(1-X)\log (1-p)\right]\right)^2\right]\\ &=E_p\left[\left(\frac Xp-\frac{1-X}{1-p}\right)^2\right]\\ &=E_p\left[\left(\frac Xp\right)^2-2\left(\frac Xp\frac{1-X}{1-p}\right)+\left(\frac{1-X}{1-p}\right)^2\right]\\ &=E_p\left(\frac{X^2}{p^2}\right)-2E_p\left(\frac{X(1-X)}{p(1-p)}\right)+E_p\left(\frac{(1-X)^2}{(1-p)^2}\right)\\ &=\frac{1}{p}-0+\frac{1}{(1-p)}\\ &=\frac{1}{p(1-p)}\\ \end{align}\]

Since \[E_p(X^2)=0^2\cdot p_X(0)+1^2\cdot p_X(1)=0^2(1-p)+1^2p=p,\] \[E_p((1-X)^2)=1-p,\]

• For Exponential RVs \(X \sim \exp(\beta)\), If the density is \[f(X|\beta) = \frac1\beta e^{-\frac1\beta X}\] The Fisher information is \[\begin{align*} \mathcal I_{\beta} & = n\mathbb E_{\beta}\left( \left(\frac{\partial \log f(X| \beta)}{\partial \beta}\right)^2\right) \\ & = n\int_0^\infty \left(\frac{\partial \log f(X|\beta)}{\partial \beta}\right)^2 \, f(X|\beta) \, dx \\ &=n\int_0^\infty \left(\frac{\partial}{\partial \beta}\left(-\frac{1}{\beta}x+\log\frac{1}{\beta}\right)\right)^2 \, \frac{1}{\beta}e^{-\frac{1}{\beta}x} \, dx \\ & = n\int_0^\infty \left(\frac{x}{\beta^2} -\frac{1}{\beta}\right)^2 \, \frac{1}{\beta}e^{-\frac{1}{\beta}x} \, dx \\ & = n\int_0^\infty \left(\frac{x^2}{\beta^4} -2\frac{x}{\beta^3}+\frac{1}{\beta^2}\right) \, \frac{1}{\beta}e^{-\frac{1}{\beta}x} \, dx \\ &=n\left(\frac{2}{\beta^2}-\frac{2}{\beta^2}+\frac{1}{\beta^2}\right)\\ &=\frac{n}{\beta^2}\\ \end{align*}\]

• For Exponential RVs \(X \sim \exp(\lambda)\), If the density is \[f(X|\lambda) = \lambda e^{-\lambda X}\] The Fisher information is \[\begin{align*} \mathcal I_{\lambda} & = n\mathbb E_{\lambda}\left( \left(\frac{\partial \log f(X| \lambda)}{\partial \lambda}\right)^2\right) \\ & = n\int_0^\infty \left(\frac{\partial \log f(X|\lambda)}{\partial \lambda}\right)^2 \, f(X|\lambda) \, dx \\ &=n\int_0^\infty \left(\frac{\partial}{\partial \lambda}\left(-\lambda x+\log\lambda\right)\right)^2 \, \lambda e^{-\lambda X}\, dx \\ & = n\int_0^\infty \left(-x +\frac{1}{\lambda}\right)^2 \, \lambda e^{-\lambda X} \, dx \\ & = n\int_0^\infty \left(x^2 -2\frac{x}{\lambda}+\frac{1}{\lambda^2}\right) \, \lambda e^{-\lambda X} \, dx \\ &=n\left(\frac{2}{\lambda^2}-\frac{2}{\lambda^2}+\frac{1}{\lambda^2}\right)\\ &=\frac{n}{\lambda^2}\\ \end{align*}\]

\[ \begin{array}{llll} & \text{Fisher Information} & \text{Cramér-Rao lower bound} \\ \hline \text{Exponential}(\beta) & \frac{n}{\beta^2} & \frac{\beta^2}{n} \\ \text{Bernoulli}(p) & \frac{1}{p(1-p)} & p(1-p) \\ \text{Normal}(\sigma^2) & \frac{n}{2\sigma^4} & \frac{2\sigma^4}{n} \\ \end{array} \]

The Matrix form

When there are \(N\) parameters, so that \(\theta\) is an \(N \times 1\) vector \[\displaystyle\boldsymbol\theta =\begin{bmatrix} \theta _{1}&\theta _{2}&\dots &\theta _{N} \end{bmatrix}^{T}\] then the Fisher information takes the form of an \(N \times N\) matrix. This matrix is called the Fisher information matrix (FIM) and has typical element \[[\mathcal I(\theta)]_{i,j}=\mathbb E_{\theta}\left[\left(\frac{\partial}{\partial\theta_i}\log f(X|\theta)\right)\left(\frac{\partial}{\partial\theta_j}\log f(X|\theta)\right)\right]\] The FIM is a \(N \times N\) positive semidefinite matrix. If it is positive definite, then it defines a Riemannian metric on the \(N\)-dimensional parameter space.

Under certain regularity conditions, the Fisher information matrix may also be written as \[[\mathcal I(\theta)]_{i,j}=-\mathbb E_{\theta}\left[\frac{\partial}{\partial\theta_i}\frac{\partial}{\partial\theta_j}\log f(X|\theta)\right]\]

[Attainment] Let \(X_1, \dots, X_n\) be iid \(X \sim f(x \vert{} \theta)\), where \(f(x \vert{} \theta)\) satisfies the conditions of the Cramér-Rao Theorem. Let \(L(\theta \vert{} \mathbf x) = \prod_{i=1}^n f(x_i \vert{} \theta)\) denote the likelihood function. If \(W(\mathbf X)=W(X_1,\cdots,X_n)\) is any unbiased estimator of \(\tau(\theta)\) then \(W(\mathbf X)\) attains the Cramér-Rao lower bound if and only if exist \(a(\theta)\) such that \[a(\theta) [W(\mathbf x) - \tau(\theta)] = \frac{\partial}{\partial \theta} \log L(\theta \vert{} \mathbf x).\]

Proof: Recalling that \(\mathbb E_{\theta}W=\tau(\theta)\), \(\mathbb E_{\theta}\left(\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)\right)=0\), the Cramer-Rao Inequality, can be written as \[\begin{align} \left[\text{Cov}_{\theta}\left(W(\mathbf X),\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)\right)\right]^2&=\left[\mathbb E_{\theta}\bigg(W(\mathbf X)-\tau(\theta)\bigg)\left(\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)-0\right)\right]^2\\ &=\left[\mathbb E_{\theta}\left(W(\mathbf X)\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)\right)\right]^2\\ &\le\text{Var}_{\theta}\left[W(\mathbf X)\right]\text{Var}_{\theta}\left[\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)\right]\end{align}\] We can have equality if and only if \[W(\mathbf X)-\tau(\theta)\] is proportional to \(\frac{\partial}{\partial\theta}\log\prod_{i=1}^{n}f(X_i|\theta)\).

Sufficiency and Unbiasedness

Rao-Blackwell Let \(W\) be any unbiased estimator of \(\tau(\theta)\) and let \(T\) be a sufficient statistic for \(\theta\). Define \(\phi(T) = \mathbb E[W \vert{} T]\). Then \(\mathbb E_{\theta}[\phi(T)] = \tau(\theta)\) and \(\text{Var}_{\theta}[\phi(T)] \leq \text{Var}_{\theta}[W] \,\, \forall \theta\). That is, \(\phi(T)\) is a uniformly better unbiased estimator of \(\tau(\theta)\).

Proof: Since \[\mathbb E X=\mathbb E[\mathbb E(X|Y)]\] \[\text{Var} X=\text{Var}[\mathbb E(X|Y)]+\mathbb E[\text{Var} (X|Y)]\] then we have \[\tau(\theta)=\mathbb E_{\theta}W=\mathbb E_{\theta}[\mathbb E(W|T)]=\mathbb E_{\theta}\phi(T)\] so \(\phi(T)\) is unbiased for \(\tau(\theta)\). Also \[\begin{align} \text{Var}_{\theta}W&=\text{Var}_{\theta}[\mathbb E(W|T)]+\mathbb E_{\theta}[\text{Var} (W|T)]\\ &=\text{Var}_{\theta}[\phi(T)]+\mathbb E_{\theta}[\text{Var} (W|T)]\\ &\ge\text{Var}_{\theta}[\phi(T)] \end{align}\] Hence \(\phi(T)\) is uniformly better than \(W\), and it only remains to show that \(\phi(T)\) is indeed an estimator. That is, we must show that \(\phi(T) = \mathbb E[W \vert{} T]\) is a function of only the sample and, in particular, is independent of \(\theta\). But it follows from the definition of sufficiency, and the fact that \(W\) is a function only of the sample, that the distribution of \(W|T\) is independent of \(\theta\). Hence \(\phi(T)\) is a uniformly better unbiased estimator of \(\tau(\theta)\).

• Conditioning any unbiased estimator on a sufficient statistic will result in a uniform improvement, so we need consider only statistics that are functions of a sufficient statistic when looking for best unbiased estimators.

• The proof doesn’t require that the statistic we condition on is sufficient, but if it isn’t then the resulting quantity will probably depend on the parameter we are trying to estimate and not be an estimator..

Uniqueness of best unbiased estimator

If \(W\) is a best unbiased estimator of \(\tau(\theta)\) then \(W\) is unique.

Proof: Suppose \(W'\) is another best unbiased estimator, and consider the estimator \(W^*=\frac{1}{2}(W + W')\). Note that \(\mathbb E_{\theta}W^* = \tau(\theta)\) and
\[\begin{align} \text{Var}_{\theta}W^*&=\text{Var}_{\theta}\left(\frac{1}{2}W+\frac{1}{2}W'\right)\\ &=\frac{1}{4}\text{Var}_{\theta}W+\frac{1}{4}\text{Var}_{\theta}W'+\frac{1}{2}\text{Cov}_{\theta}(W,W')\\ &\le\frac{1}{4}\text{Var}_{\theta}W+\frac{1}{4}\text{Var}_{\theta}W'+\frac{1}{2}\left[\text{Var}_{\theta}(W)\text{Var}_{\theta}(W')\right]^{1/2}\\ &=\text{Var}_{\theta}W \quad\quad(\text{Var}_{\theta}W=\text{Var}_{\theta}W')\\ \end{align}\] But if the above inequality is strict, then the best unbiasedness of \(W\) is contradicted, so we must have equality for all \(\theta\). Since the inequality is an application of Cauchy-Schwarz, we can have equality only if \(W' = a(\theta) W + b(\theta)\). Now using properties of covariance, we have \[\begin{align} \text{Cov}_{\theta}(W + W')&=\text{Cov}_{\theta}(W + a(\theta) W + b(\theta))\\ &=\text{Cov}_{\theta}(W + a(\theta) W)\\ &=a(\theta)\text{Var}_{\theta}(W)\\ \end{align}\] But \[\text{Cov}_{\theta}(W + W')=\text{Var}_{\theta}(W)\] since we had equality in last Cauchy-Schwarz, hence \(a(\theta)=1\) and since \(\mathbb E_{\theta}W^* = \tau(\theta)\), we must have \(b(\theta)=0\) and \(W=W'\).

The following theorem is mostly useful to show that a given estimator isn’t best unbiased.

If \(\mathbb E_{\theta}[W] = \tau(\theta)\), \(W\) is the best unbiased estimator of \(\tau(\theta)\) if and only if \(W\) is uncorrelated with all unbiased estimators of \(0\).

Proof: If \(W\) satisfies \(\mathbb E_{\theta}[W] = \tau(\theta)\), and we have another estimator, \(U\), that satisfies \(\mathbb E_{\theta}U=0\) for all \(\theta\), that is, \(U\) is an unbiased estimator of \(0\). The estimator \[\phi_a=W+aU,\] where \(a\) is a constant, satisfies \(\mathbb E_{\theta}\phi_a = \tau(\theta)\) and hence is also an unbiased estimator of \(\tau(\theta)\). The variance of \(\phi_a\) is \[\text{Var}_{\theta}\phi_a=\text{Var}_{\theta}(W+aU)\\ =\text{Var}_{\theta}W+a^2\text{Var}_{\theta}U+2a\text{Cov}_{\theta}(W,U)\] Now, if for some \(\theta=\theta_0,\text{Cov}_{\theta_0}(W,U)<0\), then we can male \(2a\text{Cov}_{\theta_0}(W,U)+a^2\text{Var}_{\theta}U<0\) by choosing \(a\in(0, -2\text{Cov}_{\theta_0}(W,U)/\text{Var}_{\theta_0}U)\). Hence, \(\phi_a\) will be better than \(W\) at \(\theta=\theta_0\) and \(W\) can’t be best unbiased. A similar argument will show that if \(\text{Cov}_{\theta_0}(W,U)>0\) for any \(\theta_0\), \(W\) also cannot be best unbiased.

If \(W\) is best unbiased, the above argument shows that \(W\) must satisfy \(\text{Cov}_{\theta}(W,U)=0\) for all \(\theta\), for any \(U\) satisfying \(\mathbb E_{\theta}U=0\). Hence the necessity is established.

Suppose now that we have an unbiased estimator \(W\) that is uncorrelated with all unbiased estimators of \(0\). Let \(W'\) be any other estimator satisfying \(\mathbb E_{\theta}W'=\mathbb E_{\theta}W=\tau(\theta)\) we will show that \(W\) is better than \(W'\). Write \[W'=W+(W'-W),\] and calculate \[\begin{align} \text{Var}_{\theta}W'&=\text{Var}_{\theta}W+\text{Var}_{\theta}(W'-W)+2\text{Cov}_{\theta}(W,W'-W)\\ &=\text{Var}_{\theta}W+\text{Var}_{\theta}(W'-W) \end{align}\] where the last equality is true because \(W'-W\) is an unbiased estimator of \(0\) and is uncorrelated with \(W\) by assumption. Since \(\text{Var}_{\theta}(W'-W)\ge 0\), then \(\text{Var}_{\theta}W'\ge \text{Var}_{\theta}W\). Since \(W'\) is arbitrary, it follows that \(W\) is the best unbiased estimator of \(\tau(\theta)\).

The idea comes from considering \(\phi_a = W + aU\) where \(\mathbb E[U] = 0\), then considering the variance.

Unbiased estimator of \(0\) Note that an unbiased estimator of \(0\) is simply random noise (there is no information in an estimator of \(0\)). If an estimator can be improved by adding noise, then it is probably defective.

We are now in a position such that, if we can characterise all of the unbiased estimators of \(0\) then we can check if a given estimator is best unbiased. In general this is not easy and requires conditions on the distribution. However, if a distribution is complete then it admits no unbiased estimators of \(0\) other than \(0\) itself (by definition), so we will be done.

Note that due to the Rao-Blackwell theorem, only the distribution of the sufficient statistic needs to be complete (not the underlying population distribution).

completeness and best unbiasedness

Let \(T\) be a complete sufficient statistic for a parameter \(\theta\) and let \(\phi(T)\) be any estimator based only on \(T\). Then \(\phi(T)\) is the unique best unbiased estimator of its expected value.

Lehmann-Scheffé Unbiased estimators based on complete sufficient statistics are unique.

Loss Function Optimality

Likelihood Ratio Tests

Likelihood ratio test (LRT) statistic for testing \(H_0: \theta \in \Theta_0\) versus \(H_1: \theta \in \Theta^c_0\) is:

\[\lambda(\mathbf x) = \frac{\underset{\Theta_0}{\sup}L(\theta|\mathbf x)}{\underset{\Theta}{\sup}L(\theta|\mathbf x)}\] A likelihood ratio test (LRT) is any test that has a rejection region of the form \(\{\mathbf x : \lambda(\mathbf x)\le c\}\) , where \(c\) is any number satisfying \(0\le c \le 1\). \[L(\theta|x_1,..., x_n) = f(\mathbf x|\theta) = \prod^n_{i=1}f(x_i|\theta)\] the numerator is the maximum probability of the observed sample, the maximum being computed over parameters in the null hypothesis. The denominator is maximum probability of observed sample over all possible parameters. Small likelihood implies rejecting \(H_0\).

Suppose \(\hat\theta\) an MLE of \(\theta\), exists; \(\hat\theta\) is obtained by doing an unrestricted maximization of \(L(\theta|\mathbf x)\). We can also consider the MLE of \(\theta\), call it \(\hat\theta_0\), obtained by doing a restricted maximization, assuming \(\Theta_0\) is the parameter space. That is, \(\hat\theta_0=\hat\theta_0(\mathbf x)\) is the value of \(\theta\in\Theta_0\) that maximizes \(L(\theta|\mathbf x)\). Then, the LRT statistic is \[\lambda(\mathbf x)=\frac{L(\hat\theta_0|\mathbf x)}{L(\hat\theta|\mathbf x)}\]

Likelihood Ratio Tests using the sufficiency statistic

If \(T(\mathbf X)\) is a sufficient statistic for \(\theta\) and \(\lambda^*(t)\) and \(\lambda(\mathbf x)\) are the LRT statistics based on \(T\) and \(\mathbf X\), respectively, then \(\lambda^*(T(\mathbf x))=\lambda(\mathbf x)\) for every \(\mathbf x\) in the sample space.

Proof: From the Factorization Theorem, the pdf or pmf of \(\mathbf X\) can be written as \(f(\mathbf x|\theta) = g(T(\mathbf x)|\theta)h(\mathbf x)\), where \(g(t|\theta)\) is the pdf or pmf of \(T\) and \(h(\mathbf x)\) does not depend on \(\theta\). Thus \[\begin{align} \lambda(\mathbf x) &= \frac{\underset{\Theta_0}{\sup}L(\theta|\mathbf x)}{\underset{\Theta}{\sup}L(\theta|\mathbf x)}\\ &=\frac{\underset{\Theta_0}{\sup}f(\mathbf x|\theta)}{\underset{\Theta}{\sup}f(\mathbf x|\theta)}\\ &=\frac{\underset{\Theta_0}{\sup}g(T(\mathbf x)|\theta)h(\mathbf x)}{\underset{\Theta}{\sup}g(T(\mathbf x)|\theta)h(\mathbf x)}\\ &=\frac{\underset{\Theta_0}{\sup}g(T(\mathbf x)|\theta)}{\underset{\Theta}{\sup}g(T(\mathbf x)|\theta)}\\ &= \frac{\underset{\Theta_0}{\sup}L^*(\theta|T(\mathbf x))}{\underset{\Theta}{\sup}L^*(\theta|T(\mathbf x))}\\ &=\lambda^*(T(\mathbf x))\\ \end{align}\]

Bayesian tests

Normal Bayesian test: accept \(H_0\) if and only if \[\frac{1}{2} \leq P(\theta \in \Theta_0|\mathbf X) = P(\theta \leq \theta_0|\mathbf X)\]

Union-intersection and intersection-union tests

Union-intersection test: \[H_0: \theta \in \underset{\gamma \in \Gamma}\bigcap \Theta_{\gamma}\] intersection-union test: \[H_0: \theta \in \underset{\gamma \in \Gamma}\bigcup \Theta_{\gamma}\]

Error Probabilities and the Power Function

• Type I Error (FP): If \(\theta \in \Theta_0\) but the hypothesis test incorrectly decides to reject \(H_0\),
• Type II Error (FN): If \(\theta \in \Theta^c_0\) but the test decides to accept \(H_0\),

Power function: The power function of a hypothesis test with rejection region \(R\) is the function of \(\theta\) defined by \[\beta(\theta) = P_{\theta}(\mathbf X \in R)\] Ideally; - \(\beta(\theta) = 0\) for all \(\theta \in \Theta_0\) - \(\beta(\theta) = 1\) for all \(\theta \in \Theta^c_0\)

Let \(X_1,\cdots,X_n\) be a random sample from a \(n(\theta, \sigma^2)\) population, \(\sigma^2\) known. An LRT of \(H_0 : \theta\le\theta_0\) versus \(H_1 : \theta > \theta_0\) is a test that rejects \(H_0\) if \[\frac{\bar{X} - \theta_0}{\sigma/ \sqrt{n}} > c\] The constant \(c\) can be any positive number. The power function of this test is: \[\begin{align} \beta(\theta) &= P_{\theta}\left(\frac{\bar{X} - \theta_0}{\sigma/ \sqrt{n}} > c\right)\\ &= P_{\theta}\left(\frac{\bar{X} - \theta}{\sigma/ \sqrt{n}} > c+\frac{\theta_0 - \theta}{\sigma/ \sqrt{n}}\right)\\ &=P\left(Z > c + \frac{\theta_0 - \theta}{\sigma/ \sqrt{n}}\right) \end{align}\] where \(Z\) is a standard normal random variable, since \[\frac{\bar{X} - \theta}{\sigma/ \sqrt{n}}\sim n(0,1)\] As \(\theta\) increases from \(-\infty\to \infty\), it is easy to see that this normal probability increases from \(0\) to \(1\). Therefore, it follows that \(\beta(\theta)\) is an increasing function of \(\theta\), with \[\underset{\theta\to-\infty}{\lim}\beta(\theta)=0, \;\;\underset{\theta\to+\infty}{\lim}\beta(\theta)=1,\] and \[\beta(\theta_0)=\alpha\quad\text{if}\quad P(Z>c)=\alpha\] \(\alpha\) is the Type I Error probability.

• For \(0\le\alpha\le 1\), a test with power function \((\beta(\theta)\) is a size \(\alpha\) test if \(\sup_{\theta\in\Theta_0}\beta(\theta) = \alpha\).

• For \(0\le\alpha\le 1\), a test with power function \((\beta(\theta)\) is a level \(\alpha\) test if \(\sup_{\theta\in\Theta_0}\beta(\theta) \le \alpha\).

• Size of LRT a size \(\alpha\) LRT is constructed by choosing \(c\) such that \(\sup_{\theta\in\Theta_0}P_{\theta}(\lambda(\mathbf X)\le c) = \alpha\). How that \(c\) is determined depends on the particular problem. For example, \(\Theta_0\) consists of the single point \(\theta=\theta_0\) and \(\sqrt{n}(\bar{X} - \theta_0) \sim n(0,1)\) if \(\theta=\theta_0\). So the test \[\text{reject }H_0\text{ if }|\bar{X}-\theta_0|\ge z_{\alpha/2}/\sqrt{n}\] where \(z_{\alpha/2}\) satisfies \(P(Z\ge z_{\alpha/2})=\alpha/2\) with \(Z\sim n(0,1)\) is a size \(\alpha\) LRT.

• Size of union-intersection test Finding constants \(t_L\) and \(t_U\) such that \[\underset{\theta\in\Theta_0}{\sup}P_{\theta}\left(\frac{\bar{X}-\mu_0}{\sqrt{S^2/n}}\ge t_L\quad\text{ or }\quad\frac{\bar{X}-\mu_0}{\sqrt{S^2/n}}\le t_U\right)=\alpha\] But for any \((\mu,\sigma^2)=\theta\in\Theta_0,\;\mu=\mu_0\) and thus \[\frac{\bar{X}-\mu_0}{\sqrt{S^2/n}}=\frac{\bar{X}-\mu}{\sqrt{S^2/n}}\] has a Student’s \(t\) distribution with \(n-1\) degree of freedom. So any choice of \(t_U=t_{n-1,1-\alpha_1}\) and \(t_L=t_{n-1,\alpha_2}\) with \(\alpha_1+\alpha_2=\alpha\) will yield a test with Type I Error probability of exactly \(\alpha\) for all \(\theta\in\Theta_0\).

• Unbiased power function A test with power function \(\beta(\theta)\) is unbiased if \(\beta(\theta')\ge\beta(\theta'')\) for every \(\theta'\in\Theta_0^c\) and \(\theta''\in\Theta_0\)

Most Powerful Tests (UMP)

Let \(\mathcal C\) be a class of tests for testing \(H_0:\theta\in\Theta_0\) versus \(H_1:\theta\in\Theta_0^c\). A test in class \(\mathcal C\), with power function \(\beta(\theta)\), is a uniformly most powerful (UMP) class \(\mathcal C\) test if \(\beta(\theta)\ge \beta'(\theta)\) for every \(\theta\in\Theta_0^c\) and every \(\beta'(\theta)\) that is a power function of a test in class \(\mathcal C\).

Sizes of Union-Intersection and Intersection-Union Tests

Union-Intersection tests ( UIT) and Intersection-Union tests (IUT) can often be bounded above by the sizes of some other tests. First consider UITs. In this situation, we are testing a null hypothesis of the form \(H_0:\theta\in\Theta_0\), where \(\Theta_0=\underset{\gamma\in\Gamma}{\bigcap}\Theta_{\gamma}\). To be specific, let \(\lambda_{\gamma}(\mathbf x)\) be the LRT statistic for testing \[H_{0\gamma}:\theta\in\Theta_{\gamma}\] versus \[H_{1\gamma}:\theta\in\Theta_{\gamma}^c\] and let \(\lambda(\mathbf x)\) be the LRT statistic for testing \(H_{0}:\theta\in\Theta_{0}\) versus \(H_{1}:\theta\in\Theta_{0}^c\) Then we have the following relationships between the overall LRT and the UIT based on \(\lambda_{\gamma}(\mathbf x)\).

\(p\)-Values

Definition A \(p\)-value \(p(\mathbf X)\) is a test statistic satisfying \(0\le p(\mathbf x)\le 1\) for every sample point \(\mathbf x\). Small values of \(p(\mathbf X)\) give evidence that \(H_1\) is true. A \(p\)-value is valid if, for every \(\theta\in\Theta_0\) and every \(0\le\alpha\le1\), \[P_{\theta}(p(\mathbf X)\le\alpha)\le\alpha\]

Loss Function Optimality

Inverting a Test Statistic

Pivotal Quantities

Definition: Pivotal Quantities A random variable \(Q(\mathbf X, \theta) = Q(X_1, \cdots, X_n, \theta)\) is a pivotal quantity (or pivot) if the distribution of \(Q(\mathbf X, \theta)\) is independent of all parameters. That is, if \(\mathbf X \sim F(\mathbf x|\theta)\), then \(Q(\mathbf X, \theta)\) has the same distribution for all values of \(\theta\).

Pivoting the CDF

Size and Coverage Probability

Test-Related Optimality

\(\mathbf X\sim f(\mathbf x|\theta)\), and we construct a \(1 - \alpha\) confidence set for \(\theta\), \(C(\mathbf x)\), by inverting an acceptance region, \(A(\theta)\). The probability of coverage of \(C(\mathbf x)\), that is, the probability of true coverage, is the function of \(\theta\) given by \(P_{\theta}(\theta \in C(\mathbf X))\). The Probability of false coverage is the function of \(\theta\) and \(\theta'\) defined by \[P_{\theta}(\theta' \in C(\mathbf X)),\theta\ne \theta', \text{ if }C(\mathbf X))=[L(\mathbf X), U(\mathbf X)],\] \[P_{\theta}(\theta' \in C(\mathbf X)),\theta'< \theta, \text{ if }C(\mathbf X))=[L(\mathbf X), \infty),\] \[P_{\theta}(\theta' \in C(\mathbf X)),\theta'> \theta, \text{ if }C(\mathbf X))=(-\infty, U(\mathbf X)],\] the probability of covering \(\theta'\) when \(\theta\) is the true parameter.

A \(1-\alpha\) confidence set that minimizes the probability of false coverage over a class of \(1-\alpha\) confidence sets is called a uniformly most accurate (UMA) confidence set.

Consistency

Efficiency

Calculations and Comparisons

If an MLE is asymptotically efficient, the asymptotic variance in Consistency of MLEs Theorem is the Delta Method variance (without the \(1/n\) in term). Thus, we can use the Cramer-Rao Lower Bound as an approximation to the true variance of the MLE. Suppose that \(X_1,\cdots, X_n\) are iid \(f(x|\theta)\), \(\hat\theta\) is the MLE of \(\theta\), and \(I_n(\theta)=\mathbb E_{\theta}\left(\frac{\partial}{\partial\theta}\log L(\theta|\mathbf X)\right)^2\) is the information number of the sample. From the Delta Method and asymptotic efficiency of MLEs, the variance of \(h(\hat\theta)\) can be approximated by
\[\begin{align} \text{Var} (h(\hat\theta)|\theta)&\approx\frac{[h(\hat\theta)]^2}{I_n(\theta)}\\ &=\frac{[h(\hat\theta)]^2}{\mathbb E_{\theta}\left(\frac{\partial}{\partial\theta}\log L(\theta|\mathbf X)\right)^2}\\ &=\frac{[h(\hat\theta)]^2}{-\mathbb E_{\theta}\left(\frac{\partial^2}{\partial\theta^2}\log L(\theta|\mathbf X)\right)}\\ &\approx\frac{[h(\hat\theta)]^2}{-\frac{\partial^2}{\partial\theta^2}\log L(\theta|\mathbf X)|_{\theta=\hat{\theta}}}\\ \end{align}\] See Fisher information (1).

The denominator is \(\hat I_n(\hat\theta)\), the observed information number. To estimate \(\text{Var}_{\theta}h(\hat\theta)\), first we approximate \(\text{Var}_{\theta}h(\hat\theta)\); then we estimate the resulting approximation, usually by substituting \(\hat\theta\) for \(\theta\). The resulting estimate can be denoted by \(\text{Var}_{\hat\theta}h(\hat\theta)\) or \(\widehat{\text{Var}}_{\theta}h(\hat\theta).\)

It follows from Theorem of Consistency of MLEs that \(-\frac{1}{n}\frac{\partial^2}{\partial\theta^2}\log L(\theta|\mathbf X)|_{\theta=\hat\theta}\) is a consistent estimator of \(I(\theta)\), so it follows that \(\text{Var}_{\hat\theta}h(\hat\theta)\) is a consistent estimator of \(\text{Var}_{\theta}h(\hat\theta)\).

Bootstrap Standard Errors

We can estimate the variance of the sample mean \(\bar{X}\) by \[\text{Var}^*(\bar{X})=\frac{1}{n^n-1}\sum_{i=1}^{n^n}(\bar{x}^*_i-\bar{\bar{x}^*})^2\] where \[\bar{\bar{x}^*}=\frac{1}{n^n}\sum_{i=1}^{n^n}\bar{x}^*_i\] the mean of the resamples.

The variance formula is applicable to virtually any estimator. Thus, for any estimator \(\hat\theta(\mathbf x)=\hat\theta\), we can write \[\text{Var}^*(\hat\theta)=\frac{1}{n^n-1}\sum_{i=1}^{n^n}(\hat\theta^*_i-\bar{\hat\theta^*})^2\] where \(\hat\theta^*_i\) is the estimator calculated from the \(i\)th resample and \[\bar{\hat\theta^*}=\frac{1}{n^n}\sum_{i=1}^{n^n}\hat\theta^*_i\] the mean of the resampled values.

The Mean and the Median

M-Estimators

Asymptotic Distribution of LRTs

Likelihood ratio method of test construction \(H_0:\theta\in \Theta_0\): \[\lambda(\mathbf x)=\frac{\underset{\Theta_0}{\sup}L(\theta|\mathbf x)}{\underset{\Theta}{\sup}L(\theta|\mathbf x)}\] and an explicit form for the rejection region, \(\{\mathbf x:\lambda(\mathbf x)\le c\}\). To define a level \(\alpha\) test, the constant \(c\) must be chosen so that \[\underset{\theta\in\Theta_0}{\sup}P_{\theta}(\lambda(\mathbf X)\le c)\le\alpha.\]

Model and Distribution Assumptions

The Classic ANOVA Hypothesis

• Hypothesises:
  • \(H_0: \theta_1 = \theta_2 = ... = \theta_k\)
  • \(H_1: \theta_i \neq \theta_j\) for some \(i, j\)

Inferences Regarding Linear Combinations of Means

Working under the oneway ANOVA assumptions, we have that \[Y_{ij}\sim n(\theta_i,\sigma^2), i=1,\cdots,k;\; j=1,\cdots,n_i.\] Therefore, \[\bar{Y}_{i\cdot}=\frac{1}{n_i}\sum_{j=1}^{n_i}Y_{ij}\sim n(\theta_i,\sigma^2/n_i),\;\; i=1,\cdots,k.\] \[\bar{\bar{Y}}=\frac{1}{N}\sum_{i=1}^{k}\sum_{j=1}^{n_i}Y_{ij},\] For any constants \(\mathbf a=(a_1,\cdots,a_k),\;\; \sum_{i=1}^{k}a_i\bar{Y}_{i\cdot}\) is also normal with \[\mathbb E\left(\sum_{i=1}^{k}a_i\bar{Y}_{i\cdot}\right)=\sum_{i=1}^{k}a_i\theta_{i}\quad\text{and}\quad\text{Var}\left(\sum_{i=1}^{k}a_i\bar{Y}_{i\cdot}\right)=\sigma^2\sum_{i=1}^{k}\frac{a_i^2}{n_i}\] and therefore

\[\frac{\sum_{i=1}^{k}a_i\bar{Y}_{i\cdot}-\sum_{i=1}^{k}a_i\theta_{i}}{\sqrt{\sigma^2\sum_{i=1}^{k}\frac{a_i^2}{n_i}}}\sim n(0,1).\]

In each population, if we denote the sample variance by \(S_i^2\), that is, \[S_i^2=\frac{1}{n_i-1}\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y}_{i\cdot})^2,\;\; i=1,\cdots,k,\] then \(S_i^2\) is an estimate of \(\sigma^2\) and \((n_i-1)S_i^2/\sigma^2\sim \chi^2_{n_i-1}\),

Furthermore, under the ANOVA assumptions, since each \(S_i^2\) estimates the same \(\sigma^2\), we can improve the estimators by combining them. We thus use the pooled estimator of \(\sigma^2\), \(S_p^2\)

• Pooled variance: \[S^{2}_{p} = \frac{1}{N-k}\sum^{k}_{i=1}(n_{i}-1)S^{2}_{i} = \frac{1}{N-k}\sum^{k}_{i=1}\sum^{n_{i}}_{j=1}(Y_{ij}-\bar{Y_{i\cdot}})^{2}\] \[(N-k)S^{2}_{p}/\sigma^2\sim\chi^2_{N-k}\] Also, \(S^{2}_{p}\) is independent of each \(\bar{Y_{i\cdot}}\) and thus \[\frac{\sum_{i=1}^{k}a_i\bar{Y_{i\cdot}}-\sum_{i=1}^{k}a_i\theta_i}{\sqrt{S^2_p\sum_{i=1}^{k}a_i^2/n_i}}\sim t_{N-k}\] Student’s \(t\) with \(N - k\) degrees of freedom.

• t statistic:

\[H_0:\sum_{i=1}^{k}a_i\theta_i=0\quad\text{versus}\quad H_1:\sum_{i=1}^{k}a_i\theta_i\ne0\] at level \(\alpha\), we would reject \(H_0\) if \[\left|\frac{\sum_{i=1}^{k}a_i\bar{Y_{i\cdot}}}{\sqrt{S^2_p\sum_{i=1}^{k}a_i^2/n_i}}\right|> t_{N-k,\alpha/2}\]

Furthermore, a pivot that can be inverted to give an interval estimator of \(\sum a_i\theta_i\) with probability \(1-\alpha\) \[\sum_{i=1}^{k}a_i\bar{Y_{i\cdot}}-t_{N-k,\alpha/2}\sqrt{S^2_p\sum_{i=1}^{k}a_i^2/n_i}\le\sum_{i=1}^{k}a_i\theta_i\le\sum_{i=1}^{k}a_i\bar{Y_{i\cdot}}+t_{N-k,\alpha/2}\sqrt{S^2_p\sum_{i=1}^{k}a_i^2/n_i}\]

To compare treatments \(1\) and \(2\), take \(\mathbf a =(1,-1,0,\cdots,0)\). Then to test \(H_0:\theta_1=\theta_2\quad\text{versus}\quad H_1:\theta_1\ne\theta_2\), we would reject \(H_0\) if \[\left|\frac{\bar{Y_{1\cdot}}\bar{Y_{2\cdot}}}{\sqrt{S^{2}_{p}(\frac{1}{n_1} + \frac{1}{n_2}})}\right| > t_{N-k, \alpha/2}\]

The ANOVA F Test

From Theorem of contrast, the ANOVA hypothesis test can be written \[H_0:\sum_{i=1}^{k}a_i\theta_i=0\;\;\forall\mathbf a\in\mathcal A\quad\text{versus}\quad H_1:\sum_{i=1}^{k}a_i\theta_i\ne0\;\;\text{for some }\mathbf a\in\mathcal A.\] where \[\mathcal A=\left\{\mathbf a=(a_1,\cdots,a_k):\sum_{i=1}^{k}a_i=0\right\}.\] To see this more clearly as a union-intersection test, define, for each \(\mathbf a\), the set \[\Theta_{\mathbf a}=\left\{\boldsymbol\theta=(\theta_1,\cdots,\theta_k):\sum_{i=1}^{k}a_i\theta_i=0\right\}\] Then we have \[\theta\in\left\{\theta:\theta_1=\theta_2= \cdots= \theta_k\right\}\iff\theta\in\Theta_{\mathbf a}\;\;\forall\mathbf a\in\mathcal A\iff\theta\in\underset{\mathbf a\in\mathcal A}{\bigcap}\Theta_{\mathbf a},\] showing that the ANOVA null can be written as an intersection.

By the Union-intersection and intersection-union tests, we would reject \(H_0:\theta\in\underset{\mathbf a\in\mathcal A}{\bigcap}\Theta_{\mathbf a}\) if we can reject \[H_{0_{\mathbf a}}:\theta\in\Theta_{\mathbf a}\quad\text{versus}\quad H_{1_{\mathbf a}}:\theta\notin\Theta_{\mathbf a};\quad \forall \mathbf a\] We test \(H_{0_{\mathbf a}}:\theta\in\Theta_{\mathbf a}\) with \(t\) statistic of \[T_{\mathbf a}=\left|\frac{\sum_{i=1}^{k}a_i\bar{Y_{i\cdot}}-\sum_{i=1}^{k}a_i\theta_i}{\sqrt{S^2_p\sum_{i=1}^{k}a_i^2/n_i}}\right|\] We then reject \(H_{0_{\mathbf a}}\) if \(T_{\mathbf a} > k\) for some constant \(k\).

From the union-intersection methodology, it follows that if we could reject for any \(\mathbf a\), we could reject for the \(\mathbf a\) that maximizes \(T_{\mathbf a}\). Thus, the union-intersection test of the ANOVA null is to reject \(H_{0_{\mathbf a}}\) if \(\sup_{\mathbf a} T_{\mathbf a} > k\), where \(k\) is chosen so that \(P_{H_0}(\sup_{\mathbf a} T_{\mathbf a} > k) =\alpha\).

Simultaneous Estimation of Contrasts

Partitioning Sums of Squares

population regression function

\[Y_i = \beta_0 + \beta_1 x_i+\epsilon_i\]

\[E[Y_i] = \beta_0 + \beta_1 x_i\]

Assumes: \(\mathbb E[\epsilon_i] = 0\)

• Another form: \(\mathbb E[Y_i|x_i] \approx \beta_0 + \beta_1 x_i\)

Least Squares

• parameters solutions:
  • \(S_{xx}=\sum(x_i-\bar{x})^2\quad S_{xy}=\sum(x_i-\bar{x})(y_i-\bar{y})\)
  • \(\beta_0 = \bar{y} - \beta_1\bar{x}\)
  • \(\beta_1 = \frac{S_{xy}}{S_{xx}}\)

Best Linear Unbiased Estimators

• Assume the form: \[Y_i = \beta_0 + \beta_1 x_i + \epsilon_i \text{ for } i = 1,\cdots, n\]

• \(\epsilon_i,..., \epsilon_n\) are uncorrelated random variables with: \[E[\epsilon_i] = 0 \text{ and } \text{Var}[\epsilon_i] = \sigma^2\]

  • \(\epsilon_i\) are random errors, therefore \(Y_i\) are uncorrelated
  • An estimator is a linear estimator if it is of the form \(\sum d_iY_i\), where \(d_1,\cdots,d_n\) are known, fixed constants.
    
  • An unbiased estimator of the slope \(\beta_1\) must satisfy \(\mathbb E\sum d_iY_i=\beta_1\)
  • This equality is true for all \(\beta_0\) and \(\beta_1\) if and only if \(\sum_{i=1}^{n} d_i=0\) and \(\sum_{i=1}^{n}d_ix_i=1\).
  • An estimator is the best linear unbiased estimator (BLUE) if it is the linear unbiased estimator with the smallest variance, which must satisfy \(\beta_1=\mathbb E\sum_{i=1}^{n}d_iY_i\).

Because \(Y_1,\cdots, Y_n\) are uncorrelated with equal variance \(\sigma^2\), the variance of any linear estimator is given by \[\text{Var}\sum_{i=1}^{n}d_iY_i=\sum_{i=1}^{n}d_i^2\sigma^2=\sigma^2\sum_{i=1}^{n}d_i^2\] The best linear unbiased estimator of \(\beta_1\) is, therefore, defined by constants \(d_1,\cdots,d_n\) that satisfy \(\sum_{i=1}^{n} d_i=0\) and \(\sum_{i=1}^{n}d_ix_i=1\) and have the minimum value of \(\sum_{i=1}^{n}d_i^2\).

The minimizing values of the constants \(d_1,\cdots,d_n\) can now be found by using Lemma of maximum of multiplied constants: For \(\mathcal A = \left\{\mathbf a = (a_1,\cdots,a_k): \sum a_i = 0\right\}\), \[\underset{\mathbf a\in\mathcal A}{\max}\left\{\frac{\left(\sum_{i=1}^{k}a_iv_i\right)^2}{\sum_{i=1}^{k}a_i^2/c_i}\right\}=\sum_{i=1}^{k}c_i(v_i-\bar{v}_c)^2,\] where \(\bar{v}_c=\sum c_iv_i/\sum c_i\) The maximum is attained at any \(\mathbf a\) of the form \(a_i = K c_i (v_i-\bar{v}_c)\), where \(K\) is a nonzero constant.

To apply the lemma to our minimization problem \[\max\frac{\left(\sum_{i=1}^{n}d_ix_i\right)^2}{\sum_{i=1}^{n}d_i^2}\] where \(\bar{x}=\sum x_i/n\). The maximum is attained at \(d_i=K(x_i-\bar{x}), \;i=1,\cdots,n\).

Then we have \[\sum_{i=1}^{n}d_ix_i=\sum_{i=1}^{n}K(x_i-\bar{x})x_i=KS_{xx}\] The constraint \(\sum_{i=1}^{n}d_ix_i=1\) is satisfied if \(K=\frac{1}{S_{xx}}\). Therefore \[d_i=\frac{(x_i-\bar{x})}{S_{xx}}\quad i=1,\cdots,n\] Finally, \[\frac{\left(\sum_{i=1}^{n}d_ix_i\right)^2}{\sum_{i=1}^{n}d_i^2}=\frac{1}{\sum_{i=1}^{n}d_i^2}\] and \[\max\frac{\left(\sum_{i=1}^{n}d_ix_i\right)^2}{\sum_{i=1}^{n}d_i^2}=\min\sum_{i=1}^{n}d_i^2, \\ \text{under constrain} \quad\sum_{i=1}^{n} d_i=0, \quad\sum_{i=1}^{n}d_ix_i=1\] with \(d_i=(x_i-\bar x)/S_{xx}\) and the linear unbiased estimator defined by these \(d_i\)s, \[\beta_1=\mathbb E\sum_{i=1}^{n}d_iY_i=\sum_{i=1}^{n}\frac{(x_i-\bar{x})}{S_{xx}}y_i=\frac{S_{xy}}{S_{xx}}\]

The variance of \(\beta_1\) is \[\text{Var}\beta_1=\sigma^2\sum_{i=1}^{n}d_i^2=\sigma^2\sum_{i=1}^{n}(x_i-\bar x)^2/S_{xx}^2=\frac{\sigma^2}{\sum_{i=1}^{n}(x_i-\bar x)^2}\]

Since \(X_1 ,\cdots, X_n\) are values chosen by the experimenter, they can be chosen to make \(S_{xx}\) large and the variance of the estimator small that is, the experimenter can design the experiment to make the estimator more precise.

Models and Distribution Assumptions

Estimation and Testing with Normal Errors

\[Y_i\sim\mathrm{n}(\beta_0+\beta_1x_i,\sigma^2),\quad i=1,\cdots,n,\] \[Y_i=\beta_0+\beta_1x_i+\epsilon_i,\quad i=1,\cdots,n,\]

The maximum likelihood estimates of the three parameters: \(\hat\beta_1=\frac{S_{xy}}{S_{xx}}\) which is linear unbiased estimator of \(\beta_1\).

\(\hat\beta_0=\bar{y}-\hat\beta_1\bar{x}\) which is linear unbiased estimator of \(\beta_0\).

\[\hat\sigma^2=\frac{1}{n}\sum_{i=1}^{n}(y_i-\hat\beta_0-\hat\beta_1x_i)^2\] which is the RSS, evaluated at the least squares line, divided by the sample size, which is not an unbiased estimator of \(\sigma^2\)

Estimation and Prediction at a Specified \(x = x_0\)

Simultaneous Estimation and Confidence Bands

Errors in variables model (or measurement error model)

Let \(\mathbb E X_i=\xi_i\), errors in variables model: \[Y_i=\beta_0+\beta_1\xi_i+\epsilon_i,\quad \epsilon_i\sim\mathrm{n}(0,\sigma^2_{\epsilon}),\] \[X_i=\xi_i+\delta_i,\quad \delta_i\sim\mathrm{n}(0,\sigma^2_{\delta}),\] Note that the assumption of normality, although common, is not necessary. Other distributions can be used.

Functional and Structural Relationships

A Least Squares Solution: orthogonal least squares

One way to take account of the fact that the \(x\)s also have error in their measurement is to perform orthogonal least squares, that is, to find the line that minimizes orthogonal (perpendicular to the line) distances rather than vertical distances.

For a particular data point \((x', y')\), the point \((\hat x', \hat y')\) on a line \(y = a + bx\) that is closest when we measure distance orthogonally is given by \[\hat x'=\frac{by'+x'-ab}{1+b^2},\quad\text{and}\quad \hat y'=a+\frac{b}{1+b^2}(by'+x'-ab).\] and the norm line has slope \(-\frac{1}{b}\)

The total least squares problem is to minimize, over all \(a\) and \(b\), the quantity \[\begin{align} &\sum_{i=1}^{n}\left((x_i-\hat x_i)^2+(y_i-\hat y_i)^2\right)\\ &=\frac{1}{1+b^2}\sum_{i=1}^{n}\left(y_i-(a+bx_i)\right)^2\\ \end{align}\]

If \(b\) is fixed, the minimizing choice of \(a\) in the sum is \(a=\bar{y} - b\bar {x}\). Then \[\begin{align} &\sum_{i=1}^{n}\left((x_i-\hat x_i)^2+(y_i-\hat y_i)^2\right)\\ &=\frac{1}{1+b^2}\sum_{i=1}^{n}\left(y_i-(a+bx_i)\right)^2\\ &=\frac{1}{1+b^2}\sum_{i=1}^{n}\left(y_i-(\bar{y} - b\bar {x}+bx_i)\right)^2\\ &=\frac{1}{1+b^2}\sum_{i=1}^{n}\left(y_i-\bar{y}) - b(x_i-\bar {x})\right)^2\\ &=\frac{1}{1+b^2}\left[S_{yy}-2bS_{xy}+b^2S_{xx}\right]. \end{align}\]

With \[S_{xx}=\sum(x_i-\bar{x})^2,\quad S_{yy}=\sum(y_i-\bar{y})^2,\quad S_{xy}=\sum(x_i-\bar{x})(y_i-\bar{y}),\]

The minimum is achieved at \(b=\frac{-(S_{xx}-S_{yy})+\sqrt{(S_{xx}-S_{yy})^2+4S_{xy}^2}}{2S_{xy}}\)

Maximum Likelihood Estimation

The Model

\[\log\left(\frac{\pi_i}{1-\pi_i}\right)=\beta_0+\beta_1x_i\] with \[\pi_i=\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}\]