AOS chapter10 Parametric Inference

10. Parametric Inference

Parametric models are of the form

\[ \mathfrak{F} = \bigg\{ f(x; \theta) : \; \theta \in \Theta \bigg\} \]

where \(\Theta \subset \mathbb{R}^k\) is the parameter space and \(\theta = (\theta_1, \dots, \theta_k)\) is the parameter. The problem of inference then reduces to the problem of estimating parameter \(\theta\).

10.13 Exercises

Exercise 10.13.1. Let \(X_1, \dots, X_n \sim \text{Gamma}(\alpha, \beta)\). Find the method of moments estimator for \(\alpha\) and \(\beta\).

Solution.

Let \(X \sim \text{Gamma}(\alpha, \beta)\). The PDF is

\[ f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x} \quad \text{for } x > 0 \]

We have:

\[ \begin{align} \mathbb{E}(X) &= \int x f_X(x) dx \\ &= \int_0^\infty x \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x} dx \\ &= \frac{\alpha}{\beta} \int_0^\infty\frac{\beta^{\alpha + 1}}{\Gamma(\alpha + 1)} x^\alpha e^{-\beta x} dx \\ &= \frac{\alpha}{\beta} \end{align} \] We also have:

\[ \begin{align} \mathbb{E}(X^2) &= \int x^2 f_X(x) dx \\ &= \int_0^\infty x^2 \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x} dx \\ &= \frac{\alpha (\alpha + 1)}{\beta^2} \int_0^\infty\frac{\beta^{\alpha + 2}}{\Gamma(\alpha + 2)} x^{\alpha + 1} e^{-\beta x} dx \\ &= \frac{\alpha (\alpha + 1)}{\beta^2} \end{align} \] Therefore,

\[ \mathbb{V}(X) = \mathbb{E}(X^2) - \mathbb{E}(X)^2 = \frac{\alpha (\alpha + 1)}{\beta^2} - \frac{\alpha^2}{\beta^2} = \frac{\alpha}{\beta^2} \]

The first two moments are:

\[ \alpha_1 = \mathbb{E}(X) = \frac{\alpha}{\beta}\\ \alpha_2 = \mathbb{E}(X^2) = \mathbb{V}(X) + \mathbb{E}(X)^2 = \frac{\alpha}{\beta^2} + \frac{\alpha^2}{\beta^2} = \frac{\alpha(\alpha + 1)}{\beta^2} \]

We have the sample moments:

\[\hat{\alpha}_1 = \frac{1}{n}\sum_{i=1}^n X_i \quad \quad \hat{\alpha}_2 = \frac{1}{n}\sum_{i=1}^n X_i^2\]

Equating these we get:

\[ \hat{\alpha}_1 = \frac{\hat{\alpha}_n}{\hat{\beta}_n} \quad \quad \hat{\alpha}_2 = \frac{\hat{\alpha}_n(\hat{\alpha}_n + 1)}{\hat{\beta}_n^2} \]

Solving these we get the method of moments estimators for \(\alpha\) and \(\beta\):

\[ \hat{\alpha}_n = \frac{\hat{\alpha}_1^2}{\hat{\alpha}_2 - \hat{\alpha}_1^2} \quad \quad \hat{\beta}_n = \frac{\hat{\alpha}_1}{\hat{\alpha}_2 - \hat{\alpha}_1^2} \]

Exercise 10.13.2. Let \(X_1, \dots, X_n \sim \text{Uniform}(a, b)\) where \(a, b\) are unknown parameters and \(a < b\).

(a) Find the method of moments estimators for \(a\) and \(b\).

(b) Find the MLE \(\hat{a}\) and \(\hat{b}\).

(c) Let \(\tau = \int x dF(x)\). Find the MLE of \(\tau\).

(d) Let \(\hat{\tau}\) be the MLE from the previous item. Let \(\tilde{\tau}\) be the nonparametric plug-in estimator of \(\tau = \int x dF(x)\). Suppose that \(a = 1\), \(b = 3\) and \(n = 10\). Find the MSE of \(\hat{\tau}\) by simulation. Find the MSE of \(\tilde{\tau}\) analytically. Compare.

Solution.

(a)

Let \(X \sim \text{Uniform}(a, b)\). Then:

• \[\mathbb{E}(X) = \int_a^b x \frac{1}{b - a} dx = \frac{a + b}{2}\]
• \[\mathbb{E}(X^2) = \int_a^b x^2 \frac{1}{b - a} dx = \frac{a^2 + ab + b^2}{3}\]
• \[\mathbb{V}(X) = \mathbb{E}(X^2) - \mathbb{E}(X)^2 = \frac{a^2 + ab + b^2}{3} - \frac{a^2 + 2ab + b^2}{4} = \frac{(b - a)^2}{12}\] The first two moments are:

\[ \alpha_1 = \mathbb{E}(X) = \frac{a + b}{2} \\ \alpha_2 = \mathbb{E}(X^2) = \mathbb{V}(X) + \mathbb{E}(X)^2 = \frac{(b - a)^2}{12} + \frac{(a + b)^2}{4} = \frac{a^2 + ab + b^2}{3} \]

We have the sample moments:

\[\hat{\alpha}_1 = \frac{1}{n}\sum_{i=1}^n X_i \quad \quad \hat{\alpha}_2 = \frac{1}{n}\sum_{i=1}^n X_i^2\]

Equating these we get:

\[ \hat{\alpha}_1 = \frac{\hat{a} + \hat{b}}{2} \quad \quad \hat{\alpha}_2 = \frac{(\hat{b} - \hat{a})^2}{12} + \frac{(\hat{a} + \hat{b})^2}{4} \]

Solving these we get the method of moment estimators for \(a\) and \(b\):

\[ \hat{a} = \hat{\alpha}_1 - \sqrt{3}(\hat{\alpha}_1^2 - \hat{\alpha}_2) \quad \quad \hat{b} = \hat{\alpha}_1 + \sqrt{3}(\hat{\alpha}_1^2 - \hat{\alpha}_2) \]

(b)

The probability density function for each \(X_i\) is

\[f(x; (a, b)) = \begin{cases} (b - a)^{-1} & \text{if } a \leq x \leq b \\ 0 & \text{otherwise} \end{cases}\]

The likelihood function is

\[\mathcal{L}_n(a, b) = \prod_{i=1}^n f(X_i; (a, b)) = \begin{cases} (b-a)^{-n} & \text{if } a \leq X_i \leq b \text{ for all } X_i\\ 0 & \text{otherwise} \end{cases} \]

The parameters that maximize the likelihood function make the \(b - a\) as small as possible – that is, we should pick the maximum \(a\) and the minimum \(b\) for which the likelihood function is non-zero. So the MLEs are:

\[\hat{a} = \min \{X_1, \dots, X_n \} \quad \quad \hat{b} = \max \{X_1, \dots, X_n \}\]

(c)

\(\tau = \int x dF(x) = \mathbb{E}(x) = (a + b)/2\), so since the MLE is equivariant, the MLE of \(\tau\) is

\[\hat{\tau} = \frac{\hat{a} + \hat{b}}{2} = \frac{\min \{X_1, \dots, X_n\} + \max\{X_1, \dots, X_n\}}{2}\]

(d)

import numpy as np

a = 1
b = 3
n = 10

X = np.random.uniform(low=a, high=b, size=n)

tau_hat = (X.min() + X.max()) / 2

# Nonparametric bootstrap to find MSE of tau_hat
B = 10000
t_boot = np.empty(B)
for i in range(B):
    xx = np.random.choice(X, n, replace=True)
    t_boot[i] = (xx.min() + xx.max()) / 2
    
se = t_boot.std()
print("MSE for tau_hat: \t %.3f" % se)

MSE for tau_hat:     0.091

Analytically, we have:

\[ \begin{align} \mathbb{V}(\tilde{\tau}) &= \mathbb{E}(\tilde{\tau}^2) - (\mathbb{E}(\tilde{\tau}))^2 \\ &= \frac{1}{n^2}\left(\mathbb{E}\left[ \left(\sum_{i=1}^n X_i\right)^2\right] - \left(\mathbb{E}\left[\sum_{i=1}^n X_i \right]\right)^2\right) \\ &= \frac{1}{n^2}\left( \mathbb{E}\left[ \sum_{i=1}^n X_i^2 + \sum_{i=1}^n \sum_{j=1, j \neq i}^n X_i X_j \right] - \left(n \frac{a + b}{2}\right)^2\right) \\ &= \frac{1}{n^2}\left( \sum_{i=1}^n \mathbb{E}[X_i^2] + \sum_{i=1}^n \sum_{j=1, j \neq i}^n \mathbb{E}[X_i]\mathbb{E}[X_j] - \left(n \frac{a + b}{2}\right)^2\right) \\ &= \frac{1}{n^2}\left( n \frac{a^2 + ab + b^2}{3} + n(n-1) \left(\frac{a+b}{2}\right)^2 - n^2\left(\frac{a + b}{2}\right)^2\right) \\ &= \frac{1}{n^2}\left( n \frac{a^2 + ab + b^2}{3} - n \left(\frac{a+b}{2}\right)^2 \right) \\ &= \frac{1}{n} \left( \frac{a^2 + ab + b^2}{3} - \frac{a^2 + 2ab + b^2}{4}\right) \\ &= \frac{1}{n} \frac{(b - a)^2}{12} \end{align} \]

Therefore,

\[\text{se}(\tilde{\tau}) = \sqrt{\frac{1}{n} \frac{(b - a)^2}{12}}\]

se_tau_tilde = np.sqrt((1/n) * ((b - a)**2 / 12))

print("MSE for tau_tilde: \t %.3f" % se_tau_tilde)

MSE for tau_tilde:   0.183

Exercise 10.13.3. Let \(X_1, \dots, X_n \sim N(\mu, \sigma^2)\). Let \(\tau\) be the 0.95 percentile, i.e. \(\mathbb{P}(X < \tau) = 0.95\).

(a) Find the MLE of \(\tau\).

(b) Find an expression for an approximate \(1 - \alpha\) confidence interval for \(\tau\).

(c) Suppose the data are:

\[ \begin{matrix} 3.23 & -2.50 & 1.88 & -0.68 & 4.43 & 0.17 \\ 1.03 & -0.07 & -0.01 & 0.76 & 1.76 & 3.18 \\ 0.33 & -0.31 & 0.30 & -0.61 & 1.52 & 5.43 \\ 1.54 & 2.28 & 0.42 & 2.33 & -1.03 & 4.00 \\ 0.39 \end{matrix} \]

Find the MLE \(\hat{\tau}\). Find the standard error using the delta method. Find the standard error using the parametric bootstrap.

Solution.

(a)

Let \(Z \sim N(0, 1)\), so \((X - \mu) / \sigma \sim Z\). We have:

\[ \begin{align} \mathbb{P}(X < \tau) &= 0.95 \\ \mathbb{P}\left(\frac{X - \mu}{\sigma} < \frac{\tau - \mu}{\sigma}\right) &= 0.95 \\ \mathbb{P}\left(Z < \frac{\tau - \mu}{\sigma}\right) &= 0.95 \\ \frac{\tau - \mu}{\sigma} &= z_{5\%} \\ \tau &= \mu + z_{5\%} \sigma \end{align} \]

Since the MLE is equivariant, \(\hat{\tau} = \hat{\mu} + z_{5\%} \hat{\sigma}\), where \(\hat{\mu}, \hat{\sigma}\) are the MLEs for the Normal distribution parameters:

\[ \hat{\mu} = n^{-1} \sum_{i=1}^n X_i \quad \quad \hat{\sigma} = \sqrt{n^{-1} \sum_{i=1}^n (X_i - \overline{X})^2} \]

(b)

Let’s use the multiparameter delta method.

We have \(\tau = g(\mu, \sigma) = \mu + z_{5\%} \sigma\), so

\[ \nabla g = \begin{bmatrix} \partial g / \partial \mu \\ \partial g / \partial \sigma \end{bmatrix} = \begin{bmatrix} 1 \\ z_{5\%} \end{bmatrix}\].

Let \(X\) be a Normal trial. The log-likelihood is \[\ell(\sigma)=-n\log\sigma-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(x_i-\mu)^2\] Let \[\ell'(\sigma)=-n\frac{1}{\sigma}+\frac{1}{\sigma^3}\sum_{i=1}^{n}(x_i-\mu)^2=0\] get \[\hat{\sigma}=\sqrt{\frac{\sum_{i=1}^{n}(x_i-\mu)^2}{n}}\] To get the Fisher information , \[\log f(X;\sigma)=-\log\sigma-\frac{(X-\mu)^2}{2\sigma^2}\] with \[\frac{\partial^2\log f(X;\sigma)}{\partial^2\sigma}=\frac{1}{\sigma^2}-\frac{3(X-\mu)^2}{\sigma^4}\] and hence

\[I(\mu)=\frac{1}{\sigma^2}\] \[ \begin{align} I(\sigma) & = -\mathbb{E} \left( \frac{1}{\sigma^2}-\frac{3(X-\mu)^2}{\sigma^4} \right) \\ &=-\frac{1}{\sigma^2}+\frac{3\sigma^2}{\sigma^4}\\ &=\frac{2}{\sigma^2} \end{align} \] The Fisher Information Matrix for the Normal process is

\[ I_n(\mu, \sigma) = \begin{bmatrix} n / \sigma^2 & 0 \\ 0 & 2n / \sigma^2 \end{bmatrix} \]

then its inverse is

\[ J_n = I_n^{-1}(\mu, \sigma) = \frac{1}{n} \begin{bmatrix} \sigma^2 & 0 \\ 0 & \sigma^2/2 \end{bmatrix} \]

and the standard error estimate for our new parameter variable is

\[\hat{\text{se}}(\hat{\tau}) = \sqrt{(\hat{\nabla} g)^T \hat{J}_n (\hat{\nabla} g)} = \hat{\sigma} \sqrt{n^{-1}(1 + z_{5\%}^2 / 2)}\]

A \(1 - \alpha\) confidence interval for \(\hat{\tau}\), then, is

\[ C_n = \left( \hat{\mu} + \hat{\sigma}\left( z_{5\%} - z_{\alpha / 2} \sqrt{n^{-1}(1 + z_{5\%}^2 / 2)} \right), \; \hat{\mu} + \hat{\sigma}\left( z_{5\%} + z_{\alpha / 2} \sqrt{n^{-1}(1 + z_{5\%}^2 / 2)} \right) \right) \]

(c)

import numpy as np
from scipy.stats import norm

z_05 = norm.ppf(0.95)
z_025 = norm.ppf(0.975)

X = np.array([
    3.23, -2.50,  1.88, -0.68,  4.43, 0.17,
    1.03, -0.07, -0.01,  0.76,  1.76, 3.18,
    0.33, -0.31,  0.30, -0.61,  1.52, 5.43,
    1.54,  2.28,  0.42,  2.33, -1.03, 4.00,
    0.39   
])

# Estimate the MLE tau_hat

n = len(X)
mu_hat = X.mean()
sigma_hat = X.std()
tau_hat = mu_hat + z_05 * sigma_hat

print("Estimated tau: %.3f" % tau_hat)

Estimated tau: 4.180

# Confidence interval using delta method

se_tau_hat = sigma_hat * np.sqrt((1/n) * (1 + z_05 * z_05 / 2))
confidence_interval = (tau_hat - z_025 * se_tau_hat, tau_hat + z_025 * se_tau_hat)

print("Estimated tau (delta method, 95%% confidence interval): \t (%.3f, %.3f)" % confidence_interval)

Estimated tau (delta method, 95% confidence interval):   (3.088, 5.273)

# Confidence interval using parametric bootstrap

n = len(X)
mu_hat = X.mean()
sigma_hat = X.std()
tau_hat = mu_hat + z_05 * sigma_hat

B = 10000
t_boot = np.empty(B)
for i in range(B):
    xx = norm.rvs(loc=mu_hat, scale=sigma_hat, size=n)
    t_boot[i] = np.quantile(xx, 0.95)
    
se_tau_hat_bootstrap = t_boot.std()
confidence_interval = (tau_hat - z_025 * se_tau_hat_bootstrap, tau_hat + z_025 * se_tau_hat_bootstrap)

print("Estimated tau (parametric bootstrap, 95%% confidence interval): \t (%.3f, %.3f)" % confidence_interval)

Estimated tau (parametric bootstrap, 95% confidence interval):   (2.898, 5.463)

Exercise 10.13.4 Let \(X_1, \dots, X_n \sim \text{Uniform}(0, \theta)\). Show that the MLE is consistent.

Hint: Let \(Y = \max \{ X_1, \dots, X_n \}\). For any c, \(\mathbb{P}(Y < c) = \mathbb{P}(X_1 < c, X_2 < c, \dots, X_n < c) = \mathbb{P}(X_1 < c)\mathbb{P}(X_2 < c)\dots\mathbb{P}(X_n < c)\).

Solution.

The probability density function is

\[ f(x, \theta) = \mathbb{P}(Y < x) = \prod_{i = 1}^n \mathbb{P}(X_i < x) = f_{\text{Uniform}(0, \theta)}(x)^n \]

The probability density function for the original distribution is

\[ f_{\text{Uniform}(0, \theta)}(x) = \begin{cases} \theta^{-1} & \text{if } 0 \leq x \leq \theta \\ 0 & \text{otherwise} \end{cases} \]

so

\[ f(x, \theta) = \begin{cases} \theta^{-n} & \text{if } 0 \leq x \leq \theta \\ 0 & \text{otherwise} \end{cases} \]

The likelihood is maximized when \(\theta\) is as small as possible while keeping all samples within the first case, so \(\hat{\theta}_n = \max \{X_1, \dots, X_n \}\).

For a given \(\epsilon > 0\), we have

\[\mathbb{P}(\hat{\theta}_n < \theta - \epsilon) = \prod_{i=1}^n \mathbb{P}(X_i < \theta - \epsilon) = \left(1 - \frac{\epsilon}{\theta} \right)^n\]

which goes to 0 as \(n \rightarrow \infty\), so \(\lim _{n \rightarrow \infty} \hat{\theta}_n = \theta\), and thus the MLE is consistent.

Exercise 10.13.5. Let \(X_1, \dots, X_n \sim \text{Poisson}(\lambda)\). Find the method of moments estimator, the maximum likelihood estimator, and the Fisher information \(I(\lambda)\).

Solution.

The first moment is:

\[\mathbb{E}(X) = \lambda\]

We have the sample moment:

\[\hat{\alpha}_1 = \frac{1}{n} \sum_{i=1}^n X_i\]

Equating those, the method of moments estimator for \(\hat{\lambda}\) is:

\[\hat{\lambda} = \hat{\alpha_1} = \frac{1}{n} \sum_{i=1}^n X_i\]

The likelihood function is

\[\mathcal{L}_n(\lambda) = \prod_{i=1}^n f(X_i; \lambda) = \prod_{i=1}^n \frac{\lambda^{X_i}e^{-\lambda}}{(X_i)!}\]

so the log likelihood function is

\[\ell_n(\lambda) = \log \mathcal{L}_n(\lambda) = \sum_{i=1}^n (\log(\lambda^{X_i}e^{-\lambda}) - \log X_i!)\\ = \sum_{i=1}^n (X_i \log \lambda - \lambda - \log X_i!) = -n \lambda + (\log \lambda) \sum_{i=1}^n X_i - \sum_{i=1}^n \log X_i! \]

To find the MLE, we differentiate this equation with respect to 0 and equate it to 0:

\[ \frac{\partial \ell_n(\lambda)}{\partial \lambda} = 0 \\ -n + \frac{\sum_{i=1}^n X_i}{\hat{\lambda}} = 0 \\ \hat{\lambda} = \frac{1}{n} \sum_{i=1}^n X_i \]

The score function is:

\[ s(X; \lambda) = \frac{\partial \log f(X; \lambda)}{\partial \lambda} = \frac{X}{\lambda} - 1\]

and the Fisher information is:

\[ I_n(\lambda) = \sum_{i=1}^n \mathbb{V}\left( s(X_i; \lambda) \right) = \sum_{i=1}^n \mathbb{V} \left( \frac{X_i}{\lambda} - 1 \right) = \frac{1}{\lambda^2} \sum_{i=1}^n \mathbb{V}(X_i) = \frac{n}{\lambda}\]

In particular, \(I(\lambda) = I_1(\lambda) = 1 / \lambda\).

Exercise 10.13.6. Let \(X_1, \dots, X_n \sim N(\theta, 1)\). Define

\[Y_i = \begin{cases} 1 & \text{if } X_i > 0 \\ 0 & \text{if } X_i \leq 0 \end{cases}\]

Let \(\psi = \mathbb{P}(Y_1 = 1)\).

(a) Find the maximum likelihood estimate \(\hat{\psi}\) of \(\psi\).

(b) Find an approximate 95% confidence interval for \(\psi\).

(c) Define \(\overline{\psi} = (1 / n) \sum_i Y_i\). Show that \(\overline{\psi}\) is a consistent estimator of \(\psi\).

(d) Compute the asymptotic relative efficiency of \(\overline{\psi}\) to \(\hat{\psi}\). Hint: Use the delta method to get the standard error of the MLE. Then compute the standard error (i.e. the standard deviation) of \(\overline{\psi}\).

(e) Suppose that the data are not really normal. Show that \(\psi\) is not consistent. What, if anything, does \(\hat{\psi}\) converge to?

Solution.

Note that, from the definition, \(Y_1, \dots, Y_n \sim \text{Bernoulli}(\Phi(\theta))\), where \(\Phi\) is the CDF for the normal distribution. Let \(p = \Phi(\theta)\).

(a) We have \(\psi = \mathbb{P}(Y_1 = 1) = p\), so the MLE is \(\hat{\psi} = \hat{p} = \Phi(\hat{\theta}) = \Phi(\overline{X})\), where \(\overline{X} = n^{-1} \sum_{i=1}^n X_i\).

(b) Let \(g(\theta) = \Phi(\theta)\). Then \(g'(\theta) = \phi(\theta)\), where \(\phi\) is the standard normal PDF. By the delta method, \(\hat{\text{se}}(\hat{\psi}) = |g'(\hat{\theta})| \hat{\text{se}}(\hat{\theta}) = \phi(\overline{X}) n^{-1/2}\).

Then, an approximate 95% confidence interval is

\[ C_n = \left(\Phi(\overline{X}) \left(1 - \frac{z_{2.5\%}}{\sqrt{n}}\right), \; \Phi(\overline{X}) \left(1 + \frac{z_{2.5\%}}{\sqrt{n}}\right) \right)\]

(c) \(\overline{\psi}\) has mean \(p\), so consistency follows from the law of large numbers.

(d) We have \(\mathbb{V}(Y_1) = \psi (1 - \psi)\), since \(Y_1\) follows a Bernoulli distribution, so \(\mathbb{V}(\overline{\psi}) = \mathbb{V}(Y_1) / n = \psi (1 - \psi) / n\).

From (b), \(\mathbb{V}{\hat{\psi}} = \phi(\theta) / n\).

Therefore, the asymptotic relative efficiency is

\[\frac{\psi(1 - \psi)}{\phi(\theta)} = \frac{\Phi(\theta)(1 - \Phi(\theta))}{\phi(\theta)}\]

(e) By the law of large numbers, we still have that \(\overline{X}\) converges to \(\mathbb{E}(Y_1) = \mathbb{P}(Y_1 = 1) \cdot 1 + \mathbb{P}(Y_1 = 0)\cdot 0 = \mathbb{P}(Y_1 = 1) = 1 - F_X(0) = \mu_Y\). Then \(\hat{\psi} = \Phi(\overline{X})\) converges to \(\Phi(\mu_Y)\). But the true value of \(\psi\) is \(\mathbb{P}(Y_1 = 1) = 1 - F_X(0)\).

But for an arbitrary distribution \(1 - F_X(0) \neq \Phi(1 - F_X(0))\).

Exercise 10.13.7. (Comparing two treatments). \(n_1\) people are given treatment 1 and \(n_2\) people are given treatment 2. Let \(X_1\) be the number of people on treatment 1 who respond favorably to the treatment and let \(X_2\) be the number of people on treatment 2 who respond favorably. Assume that \(X_1 \sim \text{Binomial}(n_1, p_1)\), \(X_2 \sim \text{Binomial}(n_2, p_2)\). Let \(\psi = p_1 - p_2\).

(a) Find the MLE of \(\psi\).

(b) Find the Fisher Information Matrix \(I(p_1, p_2)\).

(c) Use the multiparameter delta method to find the asymptotic standard error of \(\hat{\psi}\).

(d) Suppose that \(n_1 = n_2 = 200\), \(X_1 = 160\) and \(X_2 = 148\). Find \(\hat{\psi}\). Find an approximate 90% confidence interval for \(\psi\) using (i) the delta method and (ii) the parametric bootstrap.

Solution.

(a) The MLE is equivariant, so

\[\hat{\psi} = \hat{p_1} - \hat{p_2} = \frac{X_1}{n_1} - \frac{X_2}{n_2}\]

(b)

The probability mass function is

\[f((x_1, x_2); \psi) = f_1(x_1; p_1) f_2(x_2; p_2) = \binom{n_1}{x_1} p_1^{x_1} (1 - p_1)^{n_1 - x_1} \binom{n_2}{x_2} p_2^{x_2} (1 - p_2)^{n_2 - x_2} \]

The log likelihood is

\[ \begin{align} \ell_n &= \log f((x_1, x_2); \psi) \\ &= \sum_{i=1}^2 \left[\log \binom{n_i}{x_i} + x_i \log p_i + (n_i - x_i) \log (1 - p_i)\right] \end{align}\]

Calculating the partial derivatives and their expectations,

\[ \begin{align} H_{11} & = \frac{\partial^2 \ell_n}{\partial p_1^2} = \frac{\partial}{\partial p_1} \left( \frac{x_1}{p_1} - \frac{n_1 - x_1}{1 - p_1}\right) = -\frac{x_1}{p_1^2} - \frac{n_1 - x_1}{(1 - p_1)^2} \\ \mathbb{E}[H_{11}] &= -\frac{\mathbb{E}[x_1]}{p_1^2} - \frac{\mathbb{E}[n - x_1]}{(1 - p_1)^2} = -\frac{n_1 p_1}{p_1^2} - \frac{n_1(1 - p_1)}{(1 - p_1)^2} = -\frac{n_1}{p_1} - \frac{n_1}{1 - p_1} = -\frac{n_1}{p_1(1 - p_1)} \end{align} \]

\[ \begin{align} H_{22} &= -\frac{x_2}{p_2^2} - \frac{n_2 - x_2}{(1 - p_2)^2} \\ \mathbb{E}[H_{22}] &= -\frac{n_2}{p_2(1 - p_2)} \end{align} \]

\[ H_{12} = \frac{\partial^2 \ell_n}{\partial p_1 \partial p_2} = 0\] \[ H_{21} = 0\]

So the Fisher Information Matrix is:

\[ I(p_1, p_2) = \begin{bmatrix} \frac{n_1}{p_1(1 - p_1)} & 0\\ 0 & \frac{n_2}{p_2(1 - p_2)} \end{bmatrix}\]

(c) Using the multiparameter delta method, \(g(\psi) = p_1 - p_2\), so

\[ \nabla g = \begin{bmatrix} \partial g / \partial p_1 \\ \partial g / \partial p_2 \end{bmatrix} = \begin{bmatrix} 1 \\ -1 \end{bmatrix}\]

The inverse of the Fisher Information Matrix is

\[J(p_1, p_2) = I(p_1, p_2)^{-1} = \begin{bmatrix} \frac{p_1(1 - p_1)}{n_1} & 0 \\ 0 & \frac{p_2(1 - p_2)}{n_2} \end{bmatrix}\]

Then the asymptotic standard error of \(\hat{\psi}\) is:

\[\hat{\text{se}}(\hat{\psi}) = \sqrt{(\hat{\nabla} g)^T \hat{J}_n (\hat{\nabla} g)} = \sqrt{\frac{p_1(1 - p_1)}{n_1} + \frac{p_2(1 - p_2)}{n_2}}\]

(d)

import numpy as np
from scipy.stats import norm, binom

n = 200
X1 = 160
X2 = 148

p1_hat = X1 / n
p2_hat = X2 / n
psi_hat = p1_hat - p2_hat

print("Estimated psi: \t %.3f" % psi_hat)

Estimated psi:   0.060

# Confidence using delta method

z = norm.ppf(.95)

se_delta = np.sqrt(p1_hat * (1 - p1_hat)/n + p2_hat * (1 - p2_hat) / n)
confidence_delta = (psi_hat - z * se_delta, psi_hat + z * se_delta)

print("90%% confidence interval (delta method): \t %.3f, %.3f" % confidence_delta)

90% confidence interval (delta method):      -0.009, 0.129

# Confidence using parametric bootstrap

B = 1000
xx1 = binom.rvs(n, p1_hat, size=B)
xx2 = binom.rvs(n, p2_hat, size=B)
t_boot = xx1 / n - xx2 / n

se_bootstrap = t_boot.std()
confidence_delta = (psi_hat - z * se_bootstrap, psi_hat + z * se_bootstrap)

print("90%% confidence interval (parametric bootstrap): \t %.3f, %.3f" % confidence_delta)

90% confidence interval (parametric bootstrap):      -0.010, 0.130

Exercise 10.13.8. Find the Fisher information matrix for Example 10.29:

Let \(X_1, \dots, X_n \sim N(\mu, \sigma^2)\).

Solution The log likelihood is:

\[\ell_n = \sum_i \log f(x; (\mu, \sigma)) = n \left[ \log \left( \frac{1}{\sigma \sqrt{2 \pi}} \right) + \left( -\frac{1}{2} \left(\frac{x - \mu}{\sigma} \right)^2\right) \right]\]

From this,

\[H_{11} = \frac{\partial^2 \ell_n}{\partial \mu^2} = -\frac{n}{\sigma^2}\] \[H_{22} = \frac{\partial^2 \ell_n}{\partial \sigma^2} = -\frac{n}{\sigma^2} - \frac{n}{\sigma^2} = -\frac{2n}{\sigma^2}\] \[H_{12} = H_{21} = \frac{\partial^2 \ell_n}{\partial \mu \partial \sigma} = 0\]

So the Fisher Information Matrix is

\[I(\mu, \sigma) = -\begin{bmatrix} \mathbb{E}[H_{11}] & \mathbb{E}[H_{12}] \\ \mathbb{E}[H_{21}] & \mathbb{E}[H_{22}] \end{bmatrix} = \begin{bmatrix} \frac{n}{\sigma^2} & 0 \\ 0 & \frac{2n}{\sigma^2} \end{bmatrix}\]

Exercises 10.13.9 and 10.13.10. See final exercises from chapter 9.