Session 8: Classification Analysis - KNN and Friends

# ---- Setup: the previous slide, condensed, so this cell runs on its own ----
import numpy as np
import pandas as pd
by_hand = pd.DataFrame({"X1": [0, 2, 0, 0, -1, 1], "X2": [3, 0, 1, 1, 0, 1],
                        "X3": [0, 0, 3, 2, 1, 1],
                        "Y": ["Red", "Red", "Red", "Green", "Green", "Red"]}, index=[1, 2, 3, 4, 5, 6])
by_hand["distance"] = np.sqrt(by_hand["X1"]**2 + by_hand["X2"]**2 + by_hand["X3"]**2)
by_hand = by_hand.sort_values("distance")

# ---- Predict ----
# The prediction is the most common colour among the K nearest rows.
print("K = 1 predicts:", by_hand["Y"].iloc[:1].mode()[0])
print("K = 3 predicts:", by_hand["Y"].iloc[:3].mode()[0])

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Smarket: 1,250 trading days of the S&P 500, 2001-2005 (ISLR textbook data)
URL = "https://raw.githubusercontent.com/intro-stat-learning/ISLP/main/ISLP/data/Smarket.csv"
smarket = pd.read_csv(URL)

print(smarket.shape)                         # (rows, columns)
print(smarket.head())
print(smarket["Direction"].value_counts())   # how many Up days, how many Down

# ---- Setup: Steps 1-4, condensed, so this cell runs on its own ----
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
URL = "https://raw.githubusercontent.com/intro-stat-learning/ISLP/main/ISLP/data/Smarket.csv"
smarket = pd.read_csv(URL)
smarket = smarket[smarket["Year"] >= 2002]
predictors = ["Lag2"]
new_day, history = smarket.iloc[[-3]], smarket.iloc[:-3]   # 28 Dec 2005, and every day before
train, test = history[history["Year"] < 2005], history[history["Year"] == 2005]
X_train, y_train = train[predictors], train["Direction"]
X_test, y_test = test[predictors], test["Direction"]
folds = TimeSeriesSplit(n_splits=5)

# ---- Step 6 ----

# K = 1, 3, 5, ..., 125: odd numbers, so a two-class vote can never tie. It stops at 125
# because the first round fits on only 126 days, and K cannot be larger than that.
# "kneighborsclassifier__n_neighbors" = inside the pipeline, the KNN step's n_neighbors.
k_values = list(range(1, 126, 2))
search = GridSearchCV(make_pipeline(StandardScaler(), KNeighborsClassifier()),
                      {"kneighborsclassifier__n_neighbors": k_values},
                      cv=folds)
search.fit(X_train, y_train)

best_k = search.best_params_["kneighborsclassifier__n_neighbors"]
majority = y_train.value_counts(normalize=True).max()   # share of the commonest class
print(f"best K = {best_k}")
print(f"cross-validated accuracy at that K = {search.best_score_:.3f}")
print(f"always predicting the commonest direction = {majority:.3f}   <- the bar to beat")

# ---- Setup: Steps 1-6, condensed, so this cell runs on its own ----
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
URL = "https://raw.githubusercontent.com/intro-stat-learning/ISLP/main/ISLP/data/Smarket.csv"
smarket = pd.read_csv(URL)
smarket = smarket[smarket["Year"] >= 2002]
predictors = ["Lag2"]
new_day, history = smarket.iloc[[-3]], smarket.iloc[:-3]   # 28 Dec 2005, and every day before
train, test = history[history["Year"] < 2005], history[history["Year"] == 2005]
X_train, y_train = train[predictors], train["Direction"]
X_test, y_test = test[predictors], test["Direction"]
folds = TimeSeriesSplit(n_splits=5)
k_values = list(range(1, 126, 2))
search = GridSearchCV(make_pipeline(StandardScaler(), KNeighborsClassifier()),
                      {"kneighborsclassifier__n_neighbors": k_values}, cv=folds).fit(X_train, y_train)
best_k = search.best_params_["kneighborsclassifier__n_neighbors"]
majority = y_train.value_counts(normalize=True).max()

# ---- The curve ----
plt.figure(figsize=(8, 3))
plt.plot(k_values, search.cv_results_["mean_test_score"], "o-", markersize=3,
         label="KNN, cross-validated accuracy")
plt.axhline(majority, color="red", linestyle="--", label="always predict the commonest direction")
plt.axvline(best_k, color="grey", linestyle=":", label=f"chosen K = {best_k}")
plt.xlabel("K (number of neighbours)")
plt.ylabel("accuracy on the hidden blocks")
plt.legend()
plt.tight_layout()
plt.show()          # in a script: close the plot window to let it continue

# Data with a real signal: three iris species from four flower measurements.
# The data ships with scikit-learn, so there is nothing to download.
from sklearn.datasets import load_iris
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

iris = load_iris(as_frame=True)
iris_search = GridSearchCV(make_pipeline(StandardScaler(), KNeighborsClassifier()),
                           {"kneighborsclassifier__n_neighbors": range(1, 16)}, cv=10)
iris_search.fit(iris.data, iris.target)

print("best K:", iris_search.best_params_["kneighborsclassifier__n_neighbors"])
print(f"CV accuracy: {iris_search.best_score_:.3f}")

# Predict a new flower. GridSearchCV has already refitted the best K on all 150 flowers.
import pandas as pd
new_flower = pd.DataFrame({"sepal length (cm)": [6.0], "sepal width (cm)": [2.9],
                           "petal length (cm)": [4.5], "petal width (cm)": [1.5]})
species = iris.target_names[iris_search.predict(new_flower)[0]]
print("new flower:", species)
for name, prob in zip(iris.target_names, iris_search.predict_proba(new_flower)[0]):
    print(f"   Pr({name}) = {prob:.2f}")

# KNN can predict a number too: the average of the K nearest houses' prices.
import pandas as pd
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

URL = "https://www.warin.ca/datalake/courses_data/qmibr/session5/Real_Estate_Sample.csv"
houses = pd.read_csv(URL, index_col=0).dropna()
X = houses[["Living.area", "bedrooms", "bathrooms", "year", "garage"]]
y = houses["Price"]

models = {
    "KNN, K = 5": make_pipeline(StandardScaler(), KNeighborsRegressor(n_neighbors=5)),
    "KNN, K = 25": make_pipeline(StandardScaler(), KNeighborsRegressor(n_neighbors=25)),
    "linear regression": LinearRegression(),
}
for name, model in models.items():
    # scikit-learn returns NEGATIVE errors (it always maximises); the minus undoes it
    rmse = -cross_val_score(model, X, y, cv=5, scoring="neg_root_mean_squared_error").mean()
    print(f"{name:18}  RMSE = {rmse:10,.0f}")
print(f"price standard deviation = {y.std():,.0f}   <- the cost of predicting the average")

# Predict a new house with the model of lowest RMSE above, refitted on all the houses.
best = models["KNN, K = 25"].fit(X, y)
new_house = pd.DataFrame({"Living.area": [2400], "bedrooms": [4], "bathrooms": [2.5],
                          "year": [1970], "garage": [2]})
print(f"new house (2,400 sq ft, 4 bedrooms, 2.5 bathrooms, 1970, 2-car garage): "
      f"{best.predict(new_house)[0]:,.0f}")

# ---- Setup: File 1, Steps 1-6, condensed, so this cell runs on its own ----
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.neighbors import KNeighborsClassifier
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
URL = "https://raw.githubusercontent.com/intro-stat-learning/ISLP/main/ISLP/data/Smarket.csv"
smarket = pd.read_csv(URL)
smarket = smarket[smarket["Year"] >= 2002]
predictors = ["Lag2"]
new_day, history = smarket.iloc[[-3]], smarket.iloc[:-3]   # 28 Dec 2005, and every day before
train, test = history[history["Year"] < 2005], history[history["Year"] == 2005]
X_train, y_train = train[predictors], train["Direction"]
X_test, y_test = test[predictors], test["Direction"]
folds = TimeSeriesSplit(n_splits=5)
k_values = list(range(1, 126, 2))
search = GridSearchCV(make_pipeline(StandardScaler(), KNeighborsClassifier()),
                      {"kneighborsclassifier__n_neighbors": k_values}, cv=folds).fit(X_train, y_train)
best_k = search.best_params_["kneighborsclassifier__n_neighbors"]
majority = y_train.value_counts(normalize=True).max()

# ---- The forest ----
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

# 500 trees, on every core of the machine (n_jobs=-1)
forest = RandomForestClassifier(n_estimators=500, n_jobs=-1, random_state=0)
# cv=folds: the same time-ordered folds as KNN, so the past always predicts the future
forest_accuracy = cross_val_score(forest, X_train, y_train, cv=folds).mean()

print(f"random forest CV accuracy = {forest_accuracy:.3f}")
print(f"KNN CV accuracy           = {search.best_score_:.3f}")
print(f"commonest class           = {majority:.3f}")

# Predict File 2's date, 28 December 2005, with the forest fitted on the training years
forest.fit(X_train, y_train)
new_day = pd.DataFrame({"Lag2": [0.043]})
print("forest's prediction for 28 December 2005:", forest.predict(new_day)[0],
      f"(Pr(Up) = {forest.predict_proba(new_day)[0][1]:.2f}); what happened: Up")