<jupyter_start><jupyter_text>2.5 Local Methods in High Dimensions **Mean squared error for estimating f(0):***Assume that the relationship between X and Y is: $Y = f(X)$*$$\begin{align} \text{MSE}(x_0) & = E_\tau[f(x_0) - \hat{y_0}]^2\\ & = E_\tau[(f(x_0) - E_\tau(\hat{y_0})) + (E_\tau(\hat{y_0}) - \hat{y_0})]^2\\ & = E_\tau[(E_\tau(\hat{y_0}) - \hat{y_0})^2  + 2(f(x_0) - E_\tau(\hat{y_0}))(E_\tau(\hat{y_0}) - \hat{y_0})+ (f(x_0) - E_\tau(\hat{y_0}))^2]\\ & = E_\tau[(E_\tau(\hat{y_0}) - \hat{y_0})^2] + E_\tau[(E_\tau(\hat{y_0}) - f(x_0))^2]\\ & = E_\tau[\hat{y_0} - E_\tau(\hat{y_0})]^2 + [E_\tau(\hat{y_0}) - f(x_0)]^2\\ & = Var_\tau(\hat{y_0}) + Bias^2(\hat{y_0})\end{align}$$We have broken down the MSE into two components: variance and squared bias. Such decomposition is always possible and is known as *the bias-variance decomposition*.*(2.26) Suppose that the relationship between Y and X is linear with some noise:*$$Y = X^T\beta + \varepsilon $$where $\varepsilon \sim N(0, \sigma^2)$ and we fit the model by least squares to the training data. For a test point $x_0$ we have $\hat{y_0}=x_0^T\hat{\beta}$ which can be written as $\hat{y_0} = x_0^T\beta + \sum_{i=1}^N {l_i(x_0)\varepsilon_i}$ where $l_i(x_0)$ is the $i$th element of $\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}x_0$Proof:$$\begin{equation} \hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y} \\ \hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T(\mathbf{X}\beta + \varepsilon) \\ \hat{\beta}=(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T(\mathbf{X}\beta + \varepsilon) \\ \hat{\beta}=\beta + (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\varepsilon \\\end{equation}$$and by plugging $\hat{B}$ into the linear model:$$\begin{equation} \hat{y_0} = x_0^T(\beta+(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\varepsilon) \\ \hat{y_0} = x_0^T\beta+x_0^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\varepsilon\end{equation}$$we can get $\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}x_0$ from $(x_0^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T)^T$ by using two matrix properties:1. $(\mathbf{AB})^T=\mathbf{B}^T\mathbf{A}^T$2. $(\mathbf{A}^{-1})^T = (\mathbf{A}^T)^{-1}$Under this model the least square estimates are unbiased, so the expected prediction error will be:$$\begin{align}\text{EPE}(x_0) & = E_{y_0|x_0}E_\tau(y_0-\hat{y_0})^2\\& = \text{Var}(y_0|x_0) + Var_\tau(\hat{y_0}) + \text{Bias}^2(\hat{y_0})\\& = \sigma^2 + E_{\tau}x_0^T(\mathbf{X}^T\mathbf{X})^{-1}x_0\sigma^2 + 0^2\end{align}$$*Proof*:$$\begin{align}\text{EPE}(x_0) & = E_{y_0|x_0}E_\tau(y_0-\hat{y_0})^2\\& = E_{y_0|x_0}E_\tau((y_0 - f(x_0)) + (f(x_0) - \hat{y_0}))^2\\& = E_{y_0|x_0}E_\tau(y_0 - f(x_0))^2 + 2E_{y_0|x_0}E_\tau(y_0 - f(x_0))(f(x_0) - \hat{y_0}) + E_{y_0|x_0}E_\tau(f(x_0) - \hat{y_0})^2\\& = U_1 + U_2 + U_3\end{align}$$There are three components $U_1$, $U_2$, $U_3$ and we're going to expand them as well. $U_1 = E_{y_0|x_0}E_\tau(y_0 - f(x_0))^2 = E_{y_0|x_0}(y_0-f(x_0))^2 = \sigma^2$ *Note*: $f(x_0) = E_{y_0|x_0}(y_0) $$U_2 = 2E_{y_0|x_0}E_\tau(y_0 - f(x_0))(f(x_0) - \hat{y_0}) = 0$ *Note*: $E_{y_0|x_0}(y_0-f(x_0)) = 0$*$U_3$*:$$\begin{align}U_3 & = E_{y_0|x_0}E_\tau(f(x_0) - \hat{y_0})^2\\& = E_{y_0|x_0}E_\tau((\hat{y_0} - E_\tau(\hat{y_0})) + (E_\tau(\hat{y_0}) - f(x_0)))^2\\& = E_{y_0|x_0}E_\tau(\hat{y_0} - E_\tau(\hat{y_0}))^2 + 2E_{y_0|x_0}E_\tau[(\hat{y_0} - E_\tau(\hat{y_0}))(E_\tau(\hat{y_0}) - f(x_0))] + E_{y_0|x_0}E_\tau(E_\tau(\hat{y_0}) - f(x_0))^2\\& = E_\tau(\hat{y_0} - E_\tau(\hat{y_0}))^2 + (E_\tau(\hat{y_0}) - f(x_0))^2\\& = \text{Var}_\tau(\hat{y_0}) + \text{Bias}_\tau^2(\hat{y_0}) \end{align}$$Finally if we sum all $U_i$ we get: $$\text{EPE}(x_0) = U_1+U_2+U_3 = \sigma^2 + 0 + (\text{Var}_\tau(\hat{y_0}) + \text{Bias}_\tau^2(\hat{y_0}))$$$E_\tau(\hat{y_0}) = E_\tau(x_0^T\beta + \sum_{i=1}^N {l_i(x_0)\varepsilon_i})=x_0^T\beta + E(\sum_{i=1}^N {l_i(x_0)\varepsilon_i}) = x_0^T\beta + 0$ thus $\text{Bias}_\tau{\hat{y_0}} = 0$ (2.27) and we can find variance: $$\begin{align}\text{Var}_\tau(\hat{y_0}) & = E_\tau(\hat{y_0} - E_\tau(\hat{y_0})) ^ 2\\& = E_\tau(x_0^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\varepsilon)\\& = E_\tau(x_0^T(\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\varepsilon\varepsilon^T\mathbf{X}(\mathbf{X}^T\mathbf{X})^{-1}x_0)\end{align}$$where $\varepsilon\varepsilon^T=\sigma^2\mathbf{I}_n$, so we can simplify further:$$\text{Var}_\tau(\hat{y_0}) = \sigma^2x_0^{T}E_\tau[(\mathbf{X}^T\mathbf{X})^{-1})]x_0$$(2.28) if N is large and $\tau$ were selected at random, and assuming E(X) = 0, then $\mathbf{X}^T\mathbf{X}$->$NCov(\mathbf{X})$.*Proof*: By definition of covariance $\text{Cov}(X) = E[(X-E(X))(X-E(X))^T] = E(XX^T) = \frac{\mathbf{X}^T\mathbf{X}}{N}$and we can derive that:$$\begin{align}E_{x_0}\text{EPE}(x_0) & = E_{x_0}x_0^{T}\text{Cov}^{-1}(X)x_0\sigma^2/N+\sigma^2\\& = \text{trace}[\text{Cov}^{-1}(X)\text{Cov}(x_0)]\sigma^2/N+\sigma^2\\& = \sigma^2(p/N)+\sigma^2\end{align}$$**FIGURE 2.9**<jupyter_code>%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

simulations = 10000
sample_size = 500

def least_square_error(x_0, y_0, train_x, train_y):
    X = np.c_[np.ones((len(train_x), 1)), train_x]
    beta = np.linalg.lstsq(X, train_y, rcond = None)[0]
    return (np.dot(np.array([1, *x_0]), beta) - y_0) ** 2

# 1-nearest neighbor error
def nn_error(x_0, y_0, train_x, train_y):
    X = (train_x * train_x).sum(axis=1)
    return (y_0 - train_y[X.argmin()]) ** 2

cubic_epe_ratio = []
linear_epe_ratio = []

for p in range(1, 11):
    least_square_epe = [0, 0]
    nn_epe = [0, 0]
    for _ in range(simulations):
        error = np.random.standard_normal(sample_size)
        train_x = np.random.uniform(-1, 1, size=(sample_size, p))
        train_y = [train_x[:, 0] + error, 
                   0.5 * (train_x[:, 0] + 1) ** 3 + error]
        x_0 = np.zeros(p)
        y_0 = [np.random.standard_normal(),
               0.5 + np.random.standard_normal()]
        for i in range(2):
            least_square_epe[i] += least_square_error(x_0, y_0[i], 
                                                      train_x, train_y[i])
            nn_epe[i] += nn_error(x_0, y_0[i], 
                                  train_x, train_y[i])
    for i in range(2):
        least_square_epe[i] /= simulations
        nn_epe[i] /= simulations
    linear_epe_ratio.append(nn_epe[0] / least_square_epe[0])
    cubic_epe_ratio.append(nn_epe[1] / least_square_epe[1])

# plot
fig = plt.figure(1, figsize = (8, 8))
axes = fig.add_subplot(1, 1, 1)

axes.set_title("Expected Prediction Error of 1NN vs. OLS")

axes.plot(np.arange(1, 11), linear_epe_ratio, '-o',
          color = 'orange', label = "Linear")

axes.plot(np.arange(1, 11), cubic_epe_ratio, '-o',
          color = 'blue', label = "Cubic")

axes.legend()
axes.set_xlabel("Dimension")
axes.set_ylabel("EPE Ration")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2573765702600058, 'fraction_numerical': 0.04659655273152206, 'mean_word_length': 4.109701492537313, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '361842', 'license': 'MIT', 'n_tokens_mistral': 3244, 'n_tokens_neox': 3119, 'n_words': 471, 'path': 'chapter-02/2.5-local-methods-in-high-dimensions.ipynb', 'repo_name': 'debajitd/the-elements-of-statistical-learning', 'size': '6846'}