## Nonlinear example

x = np.arange(.05,1,.01)
y = np.sin(1/x) + np.random.normal(0,.25,len(x))
plt.plot(x,y,'.')

"""
Explanation: Model Selection and Bias-Variance Tradeoff
StatML: Lecture 2
Prof. James Sharpnack
LinearRegression.fit
Fit in OLS solves the following, on training set
$$
\hat \beta = (X^\top X)^{-1} X^\top y
$$
where $X,y$ are $n \times p$ and $n$ arrays.

Matrix multiply
Linear solve:
$$
(X^\top X) \hat \beta = X^\top y
$$

O notation

Sequence $A_n = O(B_n)$ if there is an M > 0 such that $A_n < M B_n$ for $n$ large enough.
An algorithm takes $O(A_n)$ time if the number of basic operations (floating point arithmetic) to run is $O(A_n)$. 

The for loop to compute
$$\sum_{i=1}^n i$$
takes $O(n)$ time while computing
$$n (n + 1)/2$$
takes $O(1)$ time.
LinearRegression.fit


Matrix mult complexity
$$
(X^\top X){j,k} = \sum{i=1}^n X_{i,j} X_{i,k} 
$$
sum is $O(n)$, full complexity is $O(p^2 n)$ time.


Matrix solve - Cholesky decomposition - computational complexity O(p^3) time

Total complexity of fit: $O(n p^2 + p^3)$.

LinearRegression.predict
Given $x_{n+1}$,
$$
\hat f(x_{n+1}) = x_{n+1}^\top \hat \beta = \hat \beta_0 + \sum_{j=1}^n \hat \beta_j x_{n+1,j}.
$$
is $O(p)$ time.
Conclusions

Predict is fast: $O(p)$
Fit is slow: $O(np^2)$ (if $n > p$)
$X^\top X$ needs to be invertible (impossible if $p > n$!)
$\hat f$ is linear (restrictive)
End of explanation
"""


## Fit K neighbors regression
knn = neighbors.KNeighborsRegressor(n_neighbors=5)
X = x.reshape(-1,1)
knn.fit(X,y)
y_pred = knn.predict(X)

## 5-nearest neighbors regression
plt.plot(x,y,'.')
plt.plot(x,y_pred)

"""
Explanation: K nearest neighbors
Given a metric (d), the K-nearest neighbors of $x$ in $x_1,\ldots,x_n$ is $x_{j_1},\ldots,x_{j_K}$ such that
$$
d(x,x_{j_1}) \le d(x,x_{j_2}) \le \ldots
$$
- KNN methods fit a local model for predicting $x$ using the KNNs
- Intuition is that $Y|X=x$ is approximately iid with $Y|X=x_{j_k}$
- For regression,
$$
\hat f_K(x) = \frac 1K \sum_{k=1}^K y_{j_k}.
$$
End of explanation
"""


import pandas as pd
import numpy as np
from sklearn import model_selection, linear_model, neighbors, preprocessing, metrics
import matplotlib.pyplot as plt
import seaborn as sns

## Reading in data - using pandas!
datapath = "../../data/"
filename = datapath + 'winequality-red.csv'
wine = pd.read_csv(filename,delimiter=';')

wine.head()

wine.describe()

## Preparing training/test data

y = wine['quality'].values
X = wine.drop(['quality'],axis=1).values
test_size = .33
X_tr, X_te, y_tr, y_te = model_selection.train_test_split(X,y,test_size = .33)

X_tr.shape, X_te.shape # checking

"""
Explanation: K nearest neighbors.fit

Naive implementation, save the data
Better is sort or arrange the data
in 1D sort makes nearest neighbor lookup fast!
sort takes $O(n \log n)$ time
in higher dimensions there is specialized data structions (KDTree)

K nearest neighbors.predict

For each new point find the KNNs
if unsorted then this can take $O(n)$ time for each new point - compare to $O(p)$ for OLS
if sorted then this can take $O(\log n)$ time (different for $>1$D)
Take average of $K$ neighbors is fast

Tuning parameter

K is a tuning parameter
Know how to evaluate - use test error
Treat each K as resulting in a different model

Stop
End of explanation
"""


## Answer to 2.1

## TRAINING - DON'T TOUCH Y_te (better not to touch X_te)

## Preprocessing

stand = preprocessing.StandardScaler()
stand.fit(X_tr)
X_tr = stand.transform(X_tr)

## Fitting OLS and KNN

lr = linear_model.LinearRegression()
lr.fit(X_tr,y_tr)

knns = []
neighbors_K = np.arange(1,241,2)
M = len(neighbors_K)
for k in neighbors_K:
    knn = neighbors.KNeighborsRegressor(n_neighbors=k)
    knn.fit(X_tr,y_tr)
    knns.append(knn)
    
models = [lr] + knns

SEL_COLS = [1,10] # acid volatility, alcohol

class WinePredictor:
    """
    Custom predictor that selects vars 1 and 10 and performs knn
    """
    def __init__(self,n_neighbors=5):
        self.n_neighbors = n_neighbors
        self.knn = neighbors.KNeighborsRegressor(n_neighbors=n_neighbors)
        
    def fit(self, X_tr, y_tr):
        X_sub = X_tr[:,SEL_COLS]
        self.knn.fit(X_sub,y_tr)
        
    def predict(self,X_te):
        X_sub = X_te[:,SEL_COLS]
        return self.knn.predict(X_sub)
    
## Can do this more cleanly with Class inheritance (maybe a future lecture)

## Fitting wine predictor

wine_knns = []

for k in neighbors_K:
    knn = WinePredictor(n_neighbors=k)
    knn.fit(X_tr,y_tr)
    wine_knns.append(knn)
    
models = models + wine_knns

## TESTING - DON'T TOUCH MODELS

## transformations

X_te = stand.transform(X_te)

## prediction

MSEs = []
for m in models:
    y_pred = m.predict(X_te)
    MSEs.append(metrics.mean_squared_error(y_te,y_pred))

plt.plot(neighbors_K,MSEs[1:M+1],label='knn')
plt.plot(neighbors_K,MSEs[M+1:],label='wine pred')
plt.hlines(MSEs[0],0,neighbors_K[-1],label='OLS')
plt.xlabel('k (neighbors)')
plt.ylabel('MSE')
plt.legend()

"""
Explanation: Exercise 2.1

Look at the documentation for StandardScalar notice that there is the fit method to be applied to training data and the transform method for the test data.  Fit and transform an instance on the wine data training Xs.
Fit the following models on the training data: LinearRegression with intercept, KNN for a large range of K.  You can store the trained models in a list.
Create a custom model class (based on the KNeighborsRegressor) that selects columns 1,10 from the dataset and fits KNN only on this subset of columns.  Give it a fit and predict method that has the same input and output as that of KNeighborsRegressor.  Fit this as well for a large range of K.
Use each model to predict the test data and calculate the mean square error.  Plot this as a function of K for each KNN method and compare to linear regression.
End of explanation
"""


## Fit K neighbors regression
x = np.arange(.05,1,.01)
y = np.sin(1/x) + np.random.normal(0,.25,len(x))
knn = neighbors.KNeighborsRegressor(n_neighbors=2)
X = x.reshape(-1,1)
knn.fit(X,y)
y_pred = knn.predict(X)

## 2-nearest neighbors regression
plt.plot(x,y,'.')
plt.plot(x,y_pred)

knn = neighbors.KNeighborsRegressor(n_neighbors=21)
X = x.reshape(-1,1)
knn.fit(X,y)
y_pred = knn.predict(X)

## 21-nearest neighbors regression
plt.plot(x,y,'.')
plt.plot(x,y_pred)

plt.plot(neighbors_K,MSEs[1:M+1],label='knn')
plt.xlabel(r"""k (neighbors)
low bias $\rightarrow$ high bias
high var $\rightarrow$ low var""")
plt.ylabel('MSE')

"""
Explanation: Stop
Bias-Variance tradeoff for square error loss
(STA 251 notes scribed by Jianing Fan, Yuhan Ning)

square error loss
any predictor $f$

Risk is
$$
R(f) = \mathbb E (Y - f(X))^2.
$$
\begin{align}
    & \eta(x) = \mathbb{E}[Y|X=x] \quad \text{minimizes}\, R(f),\
    & R(\eta) = \mathbb{E}(Y-\eta(X))^2=: R^{}.
\end{align*}
\begin{align}
    R(f) &= \mathbb{E}(Y-f(X))^2 = \mathbb{E}(Y-\eta(X)+\eta(X)-f(X))^2\
    &= \mathbb{E}(Y-\eta(X))^2+\mathbb{E}(\eta(X)-f(X))^2+2\mathbb{E}(Y-\eta(X))(\eta(X)-f(X))\
    &=R^{}+\mathbb{E}(\eta(X)-f(X))^2+2\mathbb{E}(Y-\eta(X))(\eta(X)-f(X)),
\end{align*}
\begin{align}
    \mathbb{E}(Y-\eta(X))(\eta(X)-f(X)) &= \mathbb{E}\big[\mathbb{E}[(Y-\eta(X))(\eta(X)-f(X))|X]\big]\
    &= \mathbb{E}\big[\mathbb{E}(Y-\eta(X))|X\big]\
    &= \mathbb{E}\big[\mathbb{E}(Y-\mathbb{E}[Y|X])|X\big]\
    &=0.
\end{align}
$$\mathbb{E}(\hat{f}(X)-\eta(X))^2 = \mathbb{E}\left[ \left[ \hat{f}(X)-\eta(X))^2 | X \right] \right]$$ 
So let's fix the predicted point $X = x$.
Notice that $f$ is random ($\hat{f}$) because it comes from training data.  Let's condition on the training X's ($X_{[n]}$)...
\begin{align}
    \mathbb{E}(\hat{f}(x)-\eta(x))^2 &= \mathbb{E}(\hat{f}(x)-\mathbb{E}[\hat{f}(x)|X_{[n]}])^2+\mathbb{E}(\mathbb{E}[\hat{f}(x)|X_{[n]}]-\eta(x))^2+ (\text{cross-term}=0)\
    &= \mathbb{E}(\hat{f}(x)-\mathbb{E}[\hat{f}(x)|X_{[n]}])^2+\mathbb{E}(\mathbb{E}[\hat{f}(x)|X_{[n]}]-\eta(x))^2\
    &= \mathbb E\left[(\hat{f}(x)-\mathbb E[\hat f(x) | X_{[n]}])^2 \right] + \mathbb{E} \left[ \left( \mathbb E[\hat{f}(x)|X_{[n]}]-\eta(x)\right)^2 \right]\
    &= \mathbb{E}[\text{Var}(\hat{f}(x)|X_{[n]})]+\mathbb{E}[\text{Bias}^2(\hat{f}(x)|X_{[n]})],
\end{align}
where $\text{Bias}(\hat{f}(x)|X_{[n]}) = \mathbb{E}[\hat{f}(x)|X_{[n]}] -\eta(x)$.
$$
R(\hat{f}) = R^* + \mathbb{E}[\text{Bias}^2(\hat{f}(X)|X_{[n]})+\text{Var}(\hat{f}(X)|X_{[n]})].
$$
Risk is irreducible error + bias + variance
Exercise 2.2
Recall KNN regression:
$$
\hat f_K(x) = \frac 1K \sum_{k=1}^K y_{j_k(x)}.
$$
where $j_k(x)$ is the index of the kth nearest neighbor of $x$ in $X_1,\ldots,X_n$. (Consider $x$ fixed.)
Assume homoscedasticity: $\text{Var}(y_i) = \sigma^2$ for all i.
1. Derive a simple expression for $\text{Var}(\hat{f}(x)| X_{[n]})$
2. Derive a simple expression for $\text{Bias}(\hat{f}(x)|X_{[n]})$
3. What do you expect the tradeoff to be between these as $K$ gets large?
$$\text{Var}(\hat{f}(x)| X_{[n]}) = \frac{1}{K^2} \sum_{k=1}^K \text{Var}(y_{j_k(x)}| X_{[n]}) = \frac{\sigma^2}{K}$$
$$\text{Bias}(\hat{f}(x)|X_{[n]}) = \mathbb{E}[\hat{f}(x)|X_{[n]}]-\eta(x) = \frac 1K \sum_{k=1}^K \eta(X_{j_k(x)}) - \eta(x)$$

We can expect that the conditional variance decreases,
The bias may increase because $\eta(X_{j_k(x)})$ may be quite different from $\eta(x)$ as $X_{j_k(x)}$ is far from $x$.

BV tradeoff for KNN

low K
high variance
low bias
high K
low variance
high bias
End of explanation
"""


