---
title: "Homework 05 Solutions"
pagetitle: "Homework 05 Solutions"
body-classes: "lecture-page practice-page"
engine: knitr
knitr:
  opts_chunk:
    jupyter_compat: true
execute:
  enabled: true
  cache: false
  warning: false
  message: false
format:
  html:
    html-math-method: mathml
    page-layout: full
    toc: true
    toc-location: body
    toc-title: "On this page"
    toc-depth: 2
    code-fold: true
    code-summary: "Show the solution code"
    fig-align: center
---

## Question 1: Estimate the bias and variance of KNN

### Original question

::: {.callout-note appearance="simple" icon=false}
In our previous homework, we used repeated simulations to study the behavior of linear regression estimators. We will now use the same idea for KNN. Consider the regression model

$$
Y_i=f(\mathbf X_i)+\epsilon_i,
\qquad
f(\mathbf x)=0.5x_1+\sin(x_2)-0.3x_3^2,
$$

where the three covariates are independent standard normal variables and $\epsilon_i\sim\mathcal N(0,0.1^2)$ independently of the covariates and other observations. Generate $n=400$ training observations. Our goal is to estimate the mean response at

$$
\mathbf x_0=(0.5,0.7,1)^{\mathsf T}.
$$

Use Euclidean distance and consider $k=1,3,5,\ldots,29$.

a. Calculate $f(\mathbf x_0)$. Generate one training dataset and obtain the KNN prediction at $\mathbf x_0$ for each value of $k$.

b. Independently repeat the generation of the training covariates and responses 200 times. Within each repetition, use the same training data for every value of $k$. Store the predictions in a matrix with 200 rows and one column for each value of $k$. Use this matrix to estimate the bias, variance, and mean squared error of the estimated mean response. Use divisor 200 when calculating the empirical variance, and verify numerically that

$$
\text{mean squared error}
=\text{squared bias}+\text{variance}.
$$

Plot squared bias, variance, and mean squared error against $k$ in one figure. Explain the pattern and compare it with the bias-variance trade-off for ridge regression or lasso.
:::

### Solution

**a.** The target is the mean response, whose true value is

$$
f(\mathbf x_0)=0.5(0.5)+\sin(0.7)-0.3(1)^2
\approx 0.594218.
$$

We first generate one training dataset and average the responses of the nearest $k$ observations for each choice of $k$.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q1-one-r
library(FNN)
set.seed(43251)
n = 400
p = 3
k_values = seq(1, 29, 2)
x0 = matrix(c(0.5, 0.7, 1), nrow = 1)
f0 = 0.5*x0[1, 1] + sin(x0[1, 2]) - 0.3*x0[1, 3]^2

X = matrix(rnorm(n*p), nrow = n, ncol = p)
y = 0.5*X[, 1] + sin(X[, 2]) - 0.3*X[, 3]^2 + rnorm(n, sd = 0.1)
pred = numeric(length(k_values))
for (m in seq_along(k_values)) {
  k = k_values[m]
  knn.fit = knn.reg(train = X, test = x0, y = y,
                    k = k, algorithm = "brute")
  pred[m] = knn.fit$pred
}
data.frame(k = k_values, prediction = pred, truth = f0)
```

#### Python

```{python}
#| label: hw5-q1-one-python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsRegressor

np.random.seed(43251)
n = 400
p = 3
k_values = np.arange(1, 30, 2)
x0 = np.array([[0.5, 0.7, 1]])
f0 = 0.5*x0[0, 0] + np.sin(x0[0, 1]) - 0.3*x0[0, 2]**2

X = np.random.normal(size=(n, p))
y = (0.5*X[:, 0] + np.sin(X[:, 1]) - 0.3*X[:, 2]**2
     + np.random.normal(0, 0.1, n))
pred = np.zeros(len(k_values))
for m, k in enumerate(k_values):
    knn_fit = KNeighborsRegressor(n_neighbors=int(k), weights="uniform",
                                  algorithm="brute")
    knn_fit.fit(X, y)
    pred[m] = knn_fit.predict(x0)[0]
print(pd.DataFrame({"k": k_values, "prediction": pred, "truth": f0}))
```

:::

**b.** Each column of the prediction matrix contains 200 estimates of the same target. Its mean minus $f(\mathbf x_0)$ estimates bias; its variance, calculated with divisor 200, estimates sampling variance. Averaging the squared differences from $f(\mathbf x_0)$ estimates MSE. Using the same divisor makes the empirical decomposition an exact identity, apart from rounding.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q1-repeat-r
nsim = 200
pred = matrix(NA, nrow = nsim, ncol = length(k_values))
for (l in 1:nsim) {
  X = matrix(rnorm(n*p), nrow = n, ncol = p)
  y = 0.5*X[, 1] + sin(X[, 2]) - 0.3*X[, 3]^2 + rnorm(n, sd = 0.1)
  for (m in seq_along(k_values)) {
    k = k_values[m]
    knn.fit = knn.reg(train = X, test = x0, y = y,
                      k = k, algorithm = "brute")
    pred[l, m] = knn.fit$pred
  }
}
bias = colMeans(pred) - f0
variance = colMeans(sweep(pred, 2, colMeans(pred), "-")^2)
mse = colMeans((pred - f0)^2)
data.frame(k = k_values, bias = bias, variance = variance, MSE = mse)
max(abs(mse - bias^2 - variance))
```

```{r}
#| label: hw5-q1-plot-r
#| fig-width: 7
#| fig-height: 4.5
#| fig-alt: "Squared bias, variance, and mean squared error at the fixed target as the number of neighbors increases."
matplot(k_values, cbind(bias^2, variance, mse), type = "b", pch = 1,
        lty = 1, col = c("darkorange", "deepskyblue3", "black"),
        xlab = "Number of neighbors k", ylab = "Squared error",
        ylim = c(0, max(mse)))
legend("topright", c("Squared bias", "Variance", "MSE"),
       col = c("darkorange", "deepskyblue3", "black"), lty = 1, pch = 1)
```

#### Python

```{python}
#| label: hw5-q1-repeat-python
nsim = 200
pred = np.full((nsim, len(k_values)), np.nan)
for l in range(nsim):
    X = np.random.normal(size=(n, p))
    y = (0.5*X[:, 0] + np.sin(X[:, 1]) - 0.3*X[:, 2]**2
         + np.random.normal(0, 0.1, n))
    for m, k in enumerate(k_values):
        knn_fit = KNeighborsRegressor(n_neighbors=int(k), weights="uniform",
                                      algorithm="brute")
        knn_fit.fit(X, y)
        pred[l, m] = knn_fit.predict(x0)[0]

bias = pred.mean(axis=0) - f0
variance = pred.var(axis=0, ddof=0)
mse = ((pred - f0)**2).mean(axis=0)
print(pd.DataFrame({"k": k_values, "bias": bias,
                    "variance": variance, "MSE": mse}))
print("Largest decomposition difference:", np.max(np.abs(mse - bias**2 - variance)))
```

```{python}
#| label: hw5-q1-plot-python
#| fig-width: 7
#| fig-height: 4.5
#| fig-alt: "Squared bias, variance, and mean squared error at the fixed target as the number of neighbors increases."
plt.figure(figsize=(7, 4.5))
plt.plot(k_values, bias**2, "o-", color="darkorange", label="Squared bias")
plt.plot(k_values, variance, "o-", color="deepskyblue", label="Variance")
plt.plot(k_values, mse, "o-", color="black", label="MSE")
plt.xlabel("Number of neighbors k")
plt.ylabel("Squared error")
plt.ylim(bottom=0)
plt.legend()
plt.show()
```

:::

Larger $k$ averages more responses, reducing variability, but it also includes observations farther from the target. Their mean responses can differ more from $f(\mathbf x_0)$, introducing bias. Here an intermediate neighborhood gives a lower MSE than either very few or many neighbors. Small irregularities in the curves are expected with 200 repetitions; monotonic changes in each component are not guaranteed.

Increasing $k$ plays a similar smoothing role to increasing the penalty in ridge regression or lasso: predictions become less flexible, with a trade-off between bias and variance. KNN changes which local responses are averaged; ridge and lasso shrink fitted coefficients. In this experiment both the training locations and their responses change across repetitions, so the variance includes changing neighborhoods as well as response noise. It is therefore not simply $0.1^2/k$, which describes the conditional noise variance when the neighbor locations are fixed.

## Question 2: Compare KNN with lasso

### Original question

::: {.callout-note appearance="simple" icon=false}
Use the response model from Question 1, including the noise standard deviation $0.1$, but generate $p=30$ covariates. The response still depends only on the first three. Consider these two settings:

- **Setting 1:** all 30 covariates are independent standard normal variables.
- **Setting 2:** the covariates have a multivariate normal distribution with mean zero and covariance matrix $\boldsymbol\Sigma$, where

  $$
  \Sigma_{jj}=1,
  \qquad
  \Sigma_{jl}=0.8\quad\text{for }j\ne l.
  $$

In each setting, generate 400 training observations and 1,000 independent test observations. Fit 5NN using Euclidean distance on the supplied covariate scale. All covariates already have population variance one.

For lasso, use a linear model containing the 30 covariates and a separately fitted, unpenalized intercept. Follow the Week 4 objective,

$$
L_\lambda(\boldsymbol\beta)
=\frac{1}{2n}\left\lVert
\widetilde{\mathbf y}-\mathbf X\boldsymbol\beta
\right\rVert_2^2
+\lambda\lVert\boldsymbol\beta\rVert_1,
$$

where $\mathbf X\in\mathbb R^{n\times30}$ contains the centered and standardized covariates and $\widetilde{\mathbf y}$ is the centered training response. Select the penalty by ten-fold cross-validation over 41 logarithmically spaced values from $10^{-4}$ to $1$. Estimate all means and scales using only the training portion of each fold, with divisor equal to that portion's sample size. After selecting the penalty, refit on all training observations and apply the fitted transformations unchanged to the test observations. Do not add nonlinear terms to the lasso model.

a. Fit both methods in each setting. Use the same training and test data for the two methods within a setting. Report their test MSE values in one table, together with the selected lasso penalty and the number of nonzero lasso slopes. You only need to generate one dataset for each setting; repeated simulations are not required.

b. Compare the results between the two settings. Explain how correlation can change the usefulness of Euclidean neighborhoods, and why lasso does not use all covariates in the same way as KNN.

c. The response model contains $\sin(X_2)$ and $X_3^2$, but the lasso fit uses only linear terms. Explain how this difference affects the comparison. Would the results justify a general claim that correlation always helps KNN, or that one method is always better? Explain why.
:::

### Solution

**a.** We compare the two methods on the same observations within each setting. KNN uses all 30 covariates when finding neighbors. Lasso can set slopes to zero, but it is restricted to a linear prediction rule. The following code selects the lasso penalty using training data only.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}
#### R

```{r}
#| label: hw5-q2-r
library(FNN)
library(glmnet)
set.seed(43252)

n <- 400
p <- 30
lambda <- 10^seq(0, -4, length.out = 41)
results <- data.frame(
  Setting = c("Independent", "Correlated"),
  KNN_MSE = NA_real_, Lasso_MSE = NA_real_,
  Lambda = NA_real_, Nonzero = NA_integer_
)

for (l in 1:2) {
  Sigma <- diag(p)
  if (l == 2) Sigma <- 0.8 + 0.2 * diag(p)

  X <- matrix(rnorm((n + 1000) * p), ncol = p) %*% chol(Sigma)
  y <- 0.5 * X[, 1] + sin(X[, 2]) - 0.3 * X[, 3]^2 +
    rnorm(n + 1000, sd = 0.1)
  test.X <- X[(n + 1):(n + 1000), ]
  test.y <- y[(n + 1):(n + 1000)]
  X <- X[1:n, ]
  y <- y[1:n]

  knn.fit <- knn.reg(X, test.X, y, k = 5, algorithm = "brute")
  results$KNN_MSE[l] <- mean((test.y - knn.fit$pred)^2)

  # Estimate means and scales separately in each training fold.
  fold <- sample(rep(1:10, each = n / 10))
  cv.mse <- matrix(NA_real_, nrow = 10, ncol = length(lambda))
  for (m in 1:10) {
    X.mean <- colMeans(X[fold != m, ])
    X.fit <- sweep(X[fold != m, ], 2, X.mean, "-")
    X.sd <- sqrt(colMeans(X.fit^2))
    X.fit <- sweep(X.fit, 2, X.sd, "/")
    X.valid <- sweep(X[fold == m, ], 2, X.mean, "-")
    X.valid <- sweep(X.valid, 2, X.sd, "/")
    y.mean <- mean(y[fold != m])

    lasso.fit <- glmnet(
      X.fit, y[fold != m] - y.mean, alpha = 1, lambda = lambda,
      intercept = FALSE, standardize = FALSE
    )
    y.hat <- predict(lasso.fit, newx = X.valid, s = lambda) + y.mean
    cv.mse[m, ] <- colMeans((y[fold == m] - y.hat)^2)
  }
  lambda.min <- lambda[which.min(colMeans(cv.mse))]

  # Refit using all training observations and keep these transformations.
  X.mean <- colMeans(X)
  X.fit <- sweep(X, 2, X.mean, "-")
  X.sd <- sqrt(colMeans(X.fit^2))
  X.fit <- sweep(X.fit, 2, X.sd, "/")
  test.X.fit <- sweep(test.X, 2, X.mean, "-")
  test.X.fit <- sweep(test.X.fit, 2, X.sd, "/")
  y.mean <- mean(y)
  lasso.fit <- glmnet(
    X.fit, y - y.mean, alpha = 1, lambda = lambda,
    intercept = FALSE, standardize = FALSE
  )
  test.pred <- as.vector(
    predict(lasso.fit, newx = test.X.fit, s = lambda.min)
  ) + y.mean

  results$Lasso_MSE[l] <- mean((test.y - test.pred)^2)
  results$Lambda[l] <- lambda.min
  results$Nonzero[l] <- sum(coef(lasso.fit, s = lambda.min)[-1, 1] != 0)
}
knitr::kable(results, digits = 4)
```

#### Python

```{python}
#| label: hw5-q2-python
import numpy as np
import pandas as pd
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import Lasso

np.random.seed(43252)
n = 400
p = 30
lambda_grid = np.logspace(0, -4, 41)
results = pd.DataFrame({
    "Setting": ["Independent", "Correlated"],
    "KNN_MSE": np.nan, "Lasso_MSE": np.nan,
    "Lambda": np.nan, "Nonzero": 0
})

for l in range(2):
    Sigma = np.eye(p)
    if l == 1:
        Sigma = 0.8 + 0.2 * np.eye(p)

    X = np.random.normal(size=(n + 1000, p)) @ np.linalg.cholesky(Sigma).T
    y = (0.5 * X[:, 0] + np.sin(X[:, 1]) - 0.3 * X[:, 2]**2
         + np.random.normal(scale=0.1, size=n + 1000))
    test_X = X[n:, :]
    test_y = y[n:]
    X = X[:n, :]
    y = y[:n]

    knn_fit = KNeighborsRegressor(n_neighbors=5, weights="uniform",
                                  algorithm="brute")
    knn_fit.fit(X, y)
    results.loc[l, "KNN_MSE"] = np.mean((test_y - knn_fit.predict(test_X))**2)

    # Estimate means and scales separately in each training fold.
    fold = np.random.permutation(np.repeat(np.arange(10), n // 10))
    cv_mse = np.zeros((10, len(lambda_grid)))
    for m in range(10):
        X_mean = X[fold != m, :].mean(axis=0)
        X_sd = X[fold != m, :].std(axis=0, ddof=0)
        X_fit = (X[fold != m, :] - X_mean) / X_sd
        X_valid = (X[fold == m, :] - X_mean) / X_sd
        y_mean = y[fold != m].mean()

        for h in range(len(lambda_grid)):
            lasso_fit = Lasso(alpha=lambda_grid[h], fit_intercept=False,
                              tol=1e-8, max_iter=100000)
            lasso_fit.fit(X_fit, y[fold != m] - y_mean)
            y_hat = lasso_fit.predict(X_valid) + y_mean
            cv_mse[m, h] = np.mean((y[fold == m] - y_hat)**2)
    lambda_min = lambda_grid[np.argmin(cv_mse.mean(axis=0))]

    # Refit using all training observations and keep these transformations.
    X_mean = X.mean(axis=0)
    X_sd = X.std(axis=0, ddof=0)
    X_fit = (X - X_mean) / X_sd
    test_X_fit = (test_X - X_mean) / X_sd
    y_mean = y.mean()
    lasso_fit = Lasso(alpha=lambda_min, fit_intercept=False,
                      tol=1e-8, max_iter=100000)
    lasso_fit.fit(X_fit, y - y_mean)
    test_pred = lasso_fit.predict(test_X_fit) + y_mean

    results.loc[l, "Lasso_MSE"] = np.mean((test_y - test_pred)**2)
    results.loc[l, "Lambda"] = lambda_min
    results.loc[l, "Nonzero"] = np.count_nonzero(lasso_fit.coef_)
print(results.round(4).to_string(index=False))
```
:::

The scales use divisor equal to the number of observations being fitted. Centering the response and adding its mean back gives an unpenalized intercept. Here `glmnet`'s `lambda` and scikit-learn's `alpha` both use the penalty in the stated $1/(2n)$ objective; `glmnet`'s `alpha = 1` requests lasso. R and Python generate different observations, so their numerical results need not match.

**b.** In these realizations, lasso has the smaller test MSE with independent covariates, whereas KNN has the smaller test MSE with correlated covariates. With independent covariates, many coordinates that do not affect the response contribute to Euclidean distance. The nearest observations in 30-dimensional Euclidean distance need not be the nearest in the three relevant coordinates.

With strong positive correlation, the covariates share substantial variation. Euclidean neighborhoods can then be more informative about the response-relevant coordinates. Lasso treats covariates differently: it estimates a separate slope for each and can set some slopes to zero. In the correlated setting, a covariate outside the first three can also act as a proxy for a relevant covariate, and several correlated covariates can compete in the fitted model.

**c.** Lasso uses $X_2$ and $X_3$ as linear terms, so it cannot reproduce $\sin(X_2)$ or the quadratic contribution $-0.3X_3^2$. Its intercept can account for the average quadratic contribution, but not its variation across observations. KNN does not impose a linear form, although its accuracy depends on whether it finds informative neighbors. The comparison therefore reflects both neighborhood quality and the restrictions of the fitted lasso model. Two single datasets do not show that correlation always helps KNN or that either method is universally better.

## Question 3: Observed dimension and latent dimension

### Original question

::: {.callout-note appearance="simple" icon=false}
The lecture's handwritten digit example shows that KNN can work with many measured covariates. We will use a simulation to investigate why the structure of those covariates matters.

Let $p=100$ be the number of observed covariates and let $m$ be the number of latent variables. For each observation, generate

$$
\mathbf X_i=\mathbf A\mathbf Z_i+\boldsymbol\eta_i,
\qquad
Y_i=Z_{i1}+Z_{i2}+\epsilon_i,
$$

where $\mathbf A\in\mathbb R^{100\times m}$ has independent entries from Uniform$[-2,2]$. The entries of $\mathbf Z_i\in\mathbb R^m$ and $\boldsymbol\eta_i\in\mathbb R^{100}$, and the scalar $\epsilon_i$, are independent standard normal variables, also independent of $\mathbf A$. Generate $\mathbf A$ once and use it for all training and test observations in that dataset.

Compare $m=3$ and $m=30$. In each setting, generate 200 training observations and 200 independent test observations. Fit KNN using only the 100 observed covariates in $\mathbf X$, with Euclidean distance on their generated scale. The latent variables in $\mathbf Z$ are used to generate the data but are not available to the fitted model. Consider $k=2,10,18,\ldots,82$.

a. Generate one dataset for each value of $m$. Fit KNN for every value of $k$ and calculate test MSE. Within a setting, use the same training and test observations for all values of $k$.

b. Repeat the complete experiment 50 times, generating a new $\mathbf A$ and new observations in each repetition. Plot average test MSE against $k$, with one curve for each value of $m$. These curves compare the specified choices of $k$ in a controlled simulation; they are not a cross-validation procedure for selecting a final model.

c. Both settings have 100 observed covariates, and the response depends on only two latent variables. Explain why their KNN performance can still differ. What role does the noise $\boldsymbol\eta_i$ play, and why are the observed data not confined exactly to an $m$-dimensional linear space?
:::

### Solution

**a.** The observed dimension is 100 in both settings. We generate the training and test observations using the same matrix $\mathbf A$, then use only $\mathbf X$ to fit KNN.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q3-one-r
library(FNN)
set.seed(43254)
p = 100
m_values = c(3, 30)
k_values = seq(2, 82, 8)
test_mse = matrix(NA, nrow = length(k_values), ncol = 2)

for (m in m_values) {
  A = matrix(runif(p*m, -2, 2), nrow = p, ncol = m)
  Z = matrix(rnorm(400*m), nrow = 400, ncol = m)
  X = Z %*% t(A) + matrix(rnorm(400*p), nrow = 400, ncol = p)
  y = Z[, 1] + Z[, 2] + rnorm(400)
  for (k in k_values) {
    knn.fit = knn.reg(train = X[1:200, ], test = X[201:400, ],
                      y = y[1:200], k = k, algorithm = "brute")
    test_mse[match(k, k_values), match(m, m_values)] =
      mean((y[201:400] - knn.fit$pred)^2)
  }
}
data.frame(k = k_values, m3 = test_mse[, 1], m30 = test_mse[, 2])
```

#### Python

```{python}
#| label: hw5-q3-one-python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.neighbors import KNeighborsRegressor

np.random.seed(43254)
p = 100
m_values = [3, 30]
k_values = np.arange(2, 83, 8)
test_mse = np.full((len(k_values), 2), np.nan)

for s, m in enumerate(m_values):
    A = np.random.uniform(-2, 2, size=(p, m))
    Z = np.random.normal(size=(400, m))
    X = Z @ A.T + np.random.normal(size=(400, p))
    y = Z[:, 0] + Z[:, 1] + np.random.normal(size=400)
    for l, k in enumerate(k_values):
        knn_fit = KNeighborsRegressor(n_neighbors=int(k), weights="uniform",
                                      algorithm="brute")
        knn_fit.fit(X[:200], y[:200])
        test_mse[l, s] = np.mean((y[200:] - knn_fit.predict(X[200:]))**2)
print(pd.DataFrame({"k": k_values, "m=3": test_mse[:, 0],
                    "m=30": test_mse[:, 1]}))
```

:::

**b.** The following loop repeats the entire experiment, including generation of $\mathbf A$. The result array stores repetitions, choices of $k$, and the two latent dimensions along its three axes.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q3-repeat-r
nsim = 50
test_mse = array(NA, dim = c(nsim, length(k_values), 2))
for (l in 1:nsim) {
  for (m in m_values) {
    A = matrix(runif(p*m, -2, 2), nrow = p, ncol = m)
    Z = matrix(rnorm(400*m), nrow = 400, ncol = m)
    X = Z %*% t(A) + matrix(rnorm(400*p), nrow = 400, ncol = p)
    y = Z[, 1] + Z[, 2] + rnorm(400)
    for (k in k_values) {
      knn.fit = knn.reg(train = X[1:200, ], test = X[201:400, ],
                        y = y[1:200], k = k, algorithm = "brute")
      test_mse[l, match(k, k_values), match(m, m_values)] =
        mean((y[201:400] - knn.fit$pred)^2)
    }
  }
}
mean_mse = apply(test_mse, c(2, 3), mean)
data.frame(k = k_values, m3 = mean_mse[, 1], m30 = mean_mse[, 2])
```

```{r}
#| label: hw5-q3-plot-r
#| fig-width: 7
#| fig-height: 4.5
#| fig-alt: "Average test MSE across 50 repetitions for latent dimensions 3 and 30, both observed through 100 covariates."
matplot(k_values, mean_mse, type = "b", pch = c(1, 2), lty = 1,
        col = c("deepskyblue3", "darkorange"),
        xlab = "Number of neighbors k", ylab = "Average test MSE")
legend("topright", c("m = 3", "m = 30"),
       col = c("deepskyblue3", "darkorange"), lty = 1, pch = c(1, 2))
```

#### Python

```{python}
#| label: hw5-q3-repeat-python
nsim = 50
test_mse = np.full((nsim, len(k_values), 2), np.nan)
for l in range(nsim):
    for s, m in enumerate(m_values):
        A = np.random.uniform(-2, 2, size=(p, m))
        Z = np.random.normal(size=(400, m))
        X = Z @ A.T + np.random.normal(size=(400, p))
        y = Z[:, 0] + Z[:, 1] + np.random.normal(size=400)
        for h, k in enumerate(k_values):
            knn_fit = KNeighborsRegressor(n_neighbors=int(k), weights="uniform",
                                          algorithm="brute")
            knn_fit.fit(X[:200], y[:200])
            test_mse[l, h, s] = np.mean((y[200:] - knn_fit.predict(X[200:]))**2)

mean_mse = test_mse.mean(axis=0)
print(pd.DataFrame({"k": k_values, "m=3": mean_mse[:, 0],
                    "m=30": mean_mse[:, 1]}))
```

```{python}
#| label: hw5-q3-plot-python
#| fig-width: 7
#| fig-height: 4.5
#| fig-alt: "Average test MSE across 50 repetitions for latent dimensions 3 and 30, both observed through 100 covariates."
plt.figure(figsize=(7, 4.5))
plt.plot(k_values, mean_mse[:, 0], "o-", color="deepskyblue", label="m = 3")
plt.plot(k_values, mean_mse[:, 1], "^-", color="darkorange", label="m = 30")
plt.xlabel("Number of neighbors k")
plt.ylabel("Average test MSE")
plt.legend()
plt.show()
```

:::

The $m=3$ curve has lower average test MSE in this experiment. Averaging a few neighbors reduces noise, while very large neighborhoods blur differences in the regression function. With $m=30$, even a small neighborhood is less local in the latent variables. The two curves need not have their minima at the same $k$. These are simulation averages over independent experiments, not validation scores used to choose a model for a final test set.

**c.** When $m=3$, most of the structured variation in 100 covariates is driven by three latent variables. With $m=30$, distance also responds to many latent directions that do not enter the response. Two observations close in the two relevant latent coordinates may be far apart in those other directions, making Euclidean neighborhoods less useful for predicting $Y$.

The measurement noise $\boldsymbol\eta_i$ perturbs every observed coordinate and can change which observations are nearest. Without it, $\mathbf A\mathbf Z_i$ lies in the column space of $\mathbf A$, of dimension at most $m$. With independent noise in all 100 coordinates, $\mathbf X_i$ is not confined to that space. Thus the simulation concerns a lower-dimensional signal with noise, not data lying exactly in an $m$-dimensional subspace. Keeping the response dependent on two latent variables does not make distances in all observed coordinates equally informative in the two settings.

## Question 4: Classify handwritten digits

### Original question

::: {.callout-note appearance="simple" icon=false}
Use `zip.train` and `zip.test` from the R package `ElemStatLearn`, as in the lecture. The package has been retired from the active CRAN repository. If needed, install the [archived package](https://cran.r-project.org/src/contrib/Archive/ElemStatLearn/ElemStatLearn_2015.6.26.2.tar.gz) into your usual R library, then load the two datasets. Python users can use the course copies of the same [training data](https://teazrq.github.io/stat432rpy/data/knn/zip-train.csv.gz) and [test data](https://teazrq.github.io/stat432rpy/data/knn/zip-test.csv.gz). These compressed CSV files have no headers and are included in the homework ZIP under `data/knn`. Run the Python code from the folder containing `data`.

Each image has $16\times16=256$ pixel values. The first column is the digit label and the remaining columns are the pixel values. Use the first 1,000 rows of `zip.train` and the first 500 rows of `zip.test`. Use the original pixel values without standardization, and do not include the digit label in a distance calculation.

For a test image $\mathbf x_0$ and a training image $\mathbf x_i$, the two distances are

$$
d_E(\mathbf x_0,\mathbf x_i)
=\sqrt{\sum_{j=1}^{256}(x_{0j}-x_{ij})^2},
\qquad
d_M(\mathbf x_0,\mathbf x_i)
=\sum_{j=1}^{256}|x_{0j}-x_{ij}|.
$$

a. Write a 5NN classification function that allows either Euclidean or Manhattan distance. For each test image, calculate its distances to the training images directly, find the five nearest, and return the most frequent digit label. Break distance ties by the smaller training row number and voting ties by the smallest digit. Do not use a built-in distance or KNN function for this part.

b. Apply the classifier to all 500 test images using each distance. For each fit, report the test classification error and a $10\times10$ confusion matrix, with rows representing predicted digits and columns representing true digits. Keep all digit labels from 0 through 9 in the matrix, including any with a zero count. Which digits are most often confused with one another?

c. For Euclidean distance, display three incorrectly classified test images, together with their true and predicted labels. Comment on whether the mistakes are visually understandable. Compare the overall results of the two distance measures. Can they have similar total error while making mistakes on different images? Use your predictions to support your answer.
:::

### Solution

With the specified training and test rows, Euclidean 5NN makes 61 errors out of 500, and Manhattan 5NN makes 69. Both recognize most digits, but some of their errors occur on different images.

**a. Write the classifier.** For each test image, calculate the distances, sort the training rows, and count the labels of the five closest images. Passing either of the distance functions below changes the distance without changing the remaining classification rule.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q4-classifier-r
library(ElemStatLearn)

euclidean_distance <- function(x0, X) {
  sqrt(rowSums(sweep(X, 2, x0, "-")^2))
}

manhattan_distance <- function(x0, X) {
  rowSums(abs(sweep(X, 2, x0, "-")))
}

X <- zip.train[1:1000, -1]
y <- zip.train[1:1000, 1]
X_test <- zip.test[1:500, -1]
y_test <- zip.test[1:500, 1]

knn_classify <- function(X, y, X_test, distance) {
  y_hat <- integer(nrow(X_test))
  for (i in seq_len(nrow(X_test))) {
    d <- distance(X_test[i, ], X)
    neighbors <- order(d, seq_len(nrow(X)))[1:5]
    y_hat[i] <- which.max(tabulate(y[neighbors] + 1, nbins = 10)) - 1
  }
  y_hat
}

y_hat <- cbind(Euclidean = knn_classify(X, y, X_test, euclidean_distance),
               Manhattan = knn_classify(X, y, X_test, manhattan_distance))
```

#### Python

```{python}
#| label: hw5-q4-classifier-python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

zip_train = np.loadtxt("data/knn/zip-train.csv.gz", delimiter=",")
zip_test = np.loadtxt("data/knn/zip-test.csv.gz", delimiter=",")

def euclidean_distance(x0, X):
    return np.sqrt(np.sum((X - x0)**2, axis=1))

def manhattan_distance(x0, X):
    return np.sum(np.abs(X - x0), axis=1)

X = zip_train[:1000, 1:]
y = zip_train[:1000, 0].astype(int)
X_test = zip_test[:500, 1:]
y_test = zip_test[:500, 0].astype(int)

def knn_classify(X, y, X_test, distance):
    y_hat = np.zeros(X_test.shape[0], dtype=int)
    for i in range(X_test.shape[0]):
        d = distance(X_test[i], X)
        neighbors = np.argsort(d, kind="stable")[:5]
        y_hat[i] = np.argmax(np.bincount(y[neighbors], minlength=10))
    return y_hat

y_hat = np.column_stack((knn_classify(X, y, X_test, euclidean_distance),
                          knn_classify(X, y, X_test, manhattan_distance)))
```

:::

**b. Evaluate classification.** Each column of `y_hat` contains 500 predictions. The confusion matrices below have predicted digits in rows and true digits in columns.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q4-error-r
for (l in 1:2) {
  cat(colnames(y_hat)[l], "test error:", mean(y_hat[, l] != y_test), "\n")
  print(table(Predicted = factor(y_hat[, l], levels = 0:9),
              True = factor(y_test, levels = 0:9)))
}
```

#### Python

```{python}
#| label: hw5-q4-error-python
for l in range(2):
    print(["Euclidean", "Manhattan"][l],
          "test error:", np.mean(y_hat[:, l] != y_test))
    print(pd.crosstab(pd.Series(y_hat[:, l], name="Predicted"),
                      pd.Series(y_test, name="True"))
          .reindex(index=range(10), columns=range(10), fill_value=0))
```

:::

The test error is 0.122 for Euclidean distance and 0.138 for Manhattan distance. The most frequent individual confusion is a true 2 predicted as 0: seven images under Euclidean distance and six under Manhattan distance. Other frequent errors include 5 predicted as 0 or 3, 9 predicted as 7, and 4 predicted as 9. These are counts for this test set, whose digit classes have different sample sizes.

**c. Inspect mistakes and compare predictions.** We display the first three Euclidean errors, keeping the lecture's image orientation.

::: {.panel-tabset .sync-code-panels group="homework-05-language"}

#### R

```{r}
#| label: hw5-q4-images-r
#| fig-width: 9
#| fig-height: 3
incorrect <- which(y_hat[, 1] != y_test)[1:3]
old_par <- par(mfrow = c(1, 3), mar = c(1, 1, 3, 1))
for (i in incorrect) {
  image(zip2image(zip.test, i), col = gray(256:0/256), zlim = c(0, 1),
        xlab = "", ylab = "", axes = FALSE,
        main = paste("True:", y_test[i], "Predicted:", y_hat[i, 1]))
}
par(old_par)

cat("Different predictions:", sum(y_hat[, 1] != y_hat[, 2]), "\n")
table(Euclidean_correct = y_hat[, 1] == y_test,
      Manhattan_correct = y_hat[, 2] == y_test)
```

#### Python

```{python}
#| label: hw5-q4-images-python
#| fig-width: 9
#| fig-height: 3
incorrect = np.flatnonzero(y_hat[:, 0] != y_test)[:3]
plt.figure(figsize=(9, 3))
for l in range(3):
    i = incorrect[l]
    plt.subplot(1, 3, l + 1)
    plt.imshow(zip_test[i, 1:].reshape(16, 16), cmap="gray_r",
               vmin=0, vmax=1, interpolation="nearest")
    plt.title(f"True: {y_test[i]} Predicted: {y_hat[i, 0]}")
    plt.axis("off")
plt.tight_layout()
plt.show()

print("Different predictions:", np.sum(y_hat[:, 0] != y_hat[:, 1]))
print(pd.crosstab(pd.Series(y_hat[:, 0] == y_test, name="Euclidean correct"),
                  pd.Series(y_hat[:, 1] == y_test, name="Manhattan correct")))
```

:::

These images are test rows 4, 13, and 14, with true labels 6, 2, and 4, but predicted labels 0, 0, and 7. The 6 has a rounded lower loop, the 2 has an unusually curved outline, and the 4 lacks a typical horizontal crossbar. These shapes make pixel-distance mistakes plausible, even when a person can recognize the digit. Visual recognition also uses information about handwriting that this distance rule does not explicitly represent.

The methods disagree on 22 predicted labels. Euclidean distance is correct and Manhattan distance is wrong on 11 images; the reverse happens on three. On eight additional images, both are wrong but predict different digits. Thus, similar total errors can hide different mistakes. Euclidean distance has the smaller error on this split, but this comparison does not establish that it will always be better.

## Question 5: Create a plot-formatting skill

### Original question

::: {.callout-note appearance="simple" icon=false}
Create a short `SKILL.md` file for improving the formatting of STAT 432 homework plots. You may use AI to help write it.

a. Include a title, a description, and a few formatting rules. The description should limit the skill to STAT 432 homework plot improvements and require an explicit request to use it. Your rules should address margins, titles, axes, labels, colors, and sizes for plots produced in R or Python. Keep them portable: do not refer to specific homework questions.

b. Ask your agent to install the skill. Explicitly call it to improve the plot from Question 1, then compare the original and revised versions. Do you like the changes? You may revise the skill and test it a couple more times, or uninstall it if it remains unhelpful.

c. Append your final `SKILL.md` content to your submission and briefly explain your choices and experience.
:::

### Solution

There is no single correct visual style. The following is one possible `SKILL.md`:

````markdown
---
name: format-homework-plots
description: Use only when the user explicitly calls this skill to improve the formatting of a STAT 432 homework plot.
---

# STAT 432 Plot Formatting

- Preserve the plotted data, calculations, and statistical meaning.
- Leave enough margin for titles, axis labels, tick labels, and legends without clipping.
- Use a short, informative title.
- Label axes clearly, include units when relevant, and use readable tick spacing and labels.
- Use a colorblind-friendly palette and distinguish curves with line styles when helpful.
- Start with a 6.5 by 4.5 inch figure, 12-point labels, and a 14-point title; adjust for multiple panels.
- Apply these rules to the existing R or Python plotting code.
````

The title identifies the skill, and the description limits both its scope and when it should be used. The rules express reusable formatting preferences without referring to an individual homework question. Figure and font sizes provide a starting point that can be revised after viewing the result.

Save the file as `format-homework-plots/SKILL.md` and ask your agent to install it. An explicit test request is: "Use `format-homework-plots` to improve the formatting of my Question 1 plot."

Compare the original and revised plots in your report. Describe the actual changes, whether they make the plot easier to read, and which preferences you revised or kept. Keeping, revising, or uninstalling the skill can all be reasonable outcomes; the explanation should reflect your own test.