Homework 2

Solution.

The log-likelihood function is

To optimize, we take the derivative with respect to $\lambda$ and set it equal to $0$ to get

Solution.

The MLE for a Bernoulli trial is $\widehat{p} = \mean{X}$, where $X_i = 1$ if the coin lands on heads and $0$ otherwise. In our case, we have the observed value $\mean{X} = \frac{5}{30} = \frac{1}{6}$, so

Solution.

1. We just integrate with respect to $x$:

   $$\begin{aligned} \int_y^\infty \theta y^\theta x^{-\theta-1} \,\diff{x} &= \theta y^\theta \int_y^\infty x^{-\theta-1} \,\diff{x} \\ &= \theta y^\theta \left. \frac{x^{-\theta}}{-\theta} \right\rvert_y^\infty \\ &= \theta y^\theta \p{0 - \frac{y^{-\theta}}{-\theta}} \\ &= 1. \end{aligned}$$

   Note that $\lim_{x\to\infty} x^{-\theta} = 0$ since $\theta > 0$.

2. The log-likelihood function is

   $$\begin{aligned} \ln L\p{\theta} &= \ln \prod_{i=1}^n \theta y^\theta x_i^{-\theta-1} \\ &= n\ln\theta + n\theta \ln y - \p{\theta + 1} \sum_{i=1}^n \ln x_i. \end{aligned}$$

   Setting its derivative equal to $0$, we get

   $$\begin{gathered} \frac{\partial \ln L\p{\theta}}{\partial \theta} = \frac{n}{\theta} + n\ln y - \sum_{i=1}^n \ln x_i = 0 \\ \implies \boxed{\widehat{\theta} = \frac{1}{\frac{1}{n} \sum_{i=1}^n \ln X_i - \ln y}}. \end{gathered}$$

Solution.

The least squares regression line is given by $y = \widehat{\alpha} + \widehat{\beta} x$ with square residuals $\widehat{\sigma}^2$, where

The residuals are relatively small, so this is a good fit.

Python Code

I used numpy since it's faster to type (and runs faster), but you can do the computations without it if you don't want to. However, using numpy makes plotting with matplotlib a bit easier.

from matplotlib import pyplot as plt
import numpy as np

x = np.array([62, 74, 66, 30, 84, 45, 14, 70, 82, 42])
y = np.array([71, 84, 76, 36, 95, 53, 18, 80, 93, 49])

x_bar = np.mean(x)
y_bar = np.mean(y)

beta_hat = np.sum((y - y_bar) * (x - x_bar)) / np.sum((x - x_bar) ** 2)
alpha_hat = y_bar - beta_hat * x_bar
sigma2_hat = np.mean((y - alpha_hat - beta_hat * x) ** 2)

# Or without numpy:
# beta_hat = sum([(y[i] - y_bar) * (x[i] - x_bar) for i in range(n)]) / sum(
#     [(x[i] - x_bar) ** 2 for i in range(n)]
# )
# alpha_hat = y_bar - beta_hat * x_bar
# sigma_hat_2 = sum([y[i] - alpha_hat - beta_hat * x[i] for i in range(n)]) / n

print(f"alpha_hat:  {alpha_hat}")   # 2.9788895969084805
print(f"beta_hat:   {beta_hat}")    # 1.0987892865218194
print(f"sigma2_hat: {sigma2_hat}")  # 0.09218688976135175

fig = plt.figure(figsize=(8, 16))

# Plot sample and regression line
ax = fig.add_subplot(2, 1, 1)
ax.set_title(
    f"$\\alpha$ (intercept): {alpha_hat:>.2f}, $\\beta$ (slope): {beta_hat:>.2f}"
)
ax.set_xlim(0, np.max(x))
ax.set_ylim(0, np.max(y))
ax.scatter(x, y, c="black", label="sample")
ax.axline((0, alpha_hat), slope=beta_hat, label="regression line", zorder=-1)
ax.legend()

# Plot residuals
ax = fig.add_subplot(2, 1, 2)
ax.set_title("residuals")
ax.spines["top"].set_visible(False)
ax.spines["right"].set_visible(False)
ax.spines["bottom"].set_position("zero")
ax.set_ylim(-1, 1)
ax.set_yticks(np.linspace(-1, 1, 9))
ax.scatter(x, y - (alpha_hat + beta_hat * x), c="black", label="residual")
ax.legend()

Solution.

This is a straight calculation. The mean and variance of a Beta distribution are

Note that

so when computing the method of moments estimator, we can just set

instead of using the second moment. This is only slightly easier, and the computations will be similar.

From the first moment,

which we can substitute into

For $\widehat{\beta}$,

so in summary,