🧮 Solutions: Modelling

import numpy as np
from sklearn.datasets import fetch_california_housing
from mlxtend.feature_selection import SequentialFeatureSelector
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, KFold, cross_val_score


# 1) Load the California housing dataset
data = fetch_california_housing(as_frame=True)

X = data.data
y = data.target


# 2) Split the data and evaluate a baseline model with all features
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

print(X_train.shape, X_test.shape)
print(y_train.shape, y_test.shape)

cv = KFold(n_splits=5, shuffle=True, random_state=42)  # shuffle, since the rows are ordered
scores = cross_val_score(LinearRegression(), X_train, y_train, cv=cv, scoring='r2')

print("R² scores from each fold:", np.round(scores, 3))
print(f"Average R² score: {np.mean(scores):.3f}")


# 3) Apply a forward stepwise selection with the same folds
sfs_forward = SequentialFeatureSelector(
    estimator=LinearRegression(),
    k_features="parsimonious",
    forward=True,
    floating=False,
    scoring='r2',
    cv=cv,
    verbose=0)

sfs_forward.fit(X_train, y_train)

print(f">> Forward SFS:")
print(f"   Best CV R²      : {sfs_forward.k_score_:.3f}")
print(f"   Optimal # feats : {len(sfs_forward.k_feature_idx_)}")
print(f"   Feature names   : {sfs_forward.k_feature_names_}")


# 4) Evaluate both models on the test set
selected_features = list(sfs_forward.k_feature_names_)

full_model = LinearRegression().fit(X_train, y_train)
selected_model = LinearRegression().fit(X_train[selected_features], y_train)

print(f"Test R² (all features):      {full_model.score(X_test, y_test):.4f}")
print(f"Test R² (selected features): {selected_model.score(X_test[selected_features], y_test):.4f}")

import pandas as pd
import numpy as np
import statsmodels.api as sm 

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LassoCV
from sklearn.model_selection import train_test_split

# Data related processing
hitters = sm.datasets.get_rdataset("Hitters", "ISLR").data
hitters_subset = hitters[["Salary", "AtBat", "Runs","RBI", "CHits", "CAtBat", "CRuns", "CWalks", "Assists", "Hits", "HmRun", "Years", "Errors", "Walks"]].copy()

# 1) Prepare the data
hitters_subset = hitters_subset.drop(columns=["CRuns", "CAtBat"]) # Remove highly correlated features (see previous session)
hitters_subset.dropna(inplace=True) # drop rows containing missing data

y = hitters_subset["Salary"]
X = hitters_subset.drop(columns=["Salary"])
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

scaler = StandardScaler() # Scale predictors to mean=0 and std=1
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)


# 2) Find the optimal lambda with LassoCV
lambda_range = np.logspace(-2, 3, 100)

lasso_cv = LassoCV(alphas=lambda_range, cv=5, max_iter=10000)
lasso_cv.fit(X_train_scaled, y_train) 


# 3) Print the optimal lambda, the training and test R², and the coefficients
print(f"Optimal lambda: {lasso_cv.alpha_:.2f}")
print(f"Training R²:    {lasso_cv.score(X_train_scaled, y_train):.3f}")
print(f"Test R²:        {lasso_cv.score(X_test_scaled, y_test):.3f}\n")

coef_table = pd.DataFrame({
    'Predictor': X_train.columns,
    'Beta': lasso_cv.coef_
})

coef_table = coef_table.reindex(coef_table['Beta'].abs().sort_values(ascending=False).index)
print(coef_table)

from sklearn.datasets import load_diabetes

# Load the data as a DataFrame
diabetes = load_diabetes(as_frame=True)
X = diabetes.data
y = diabetes.target

# Inspect the data and check for missing values
print(X.isna().sum().sum(), "missing values")
X.describe().round(3)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# Standardize the data and fit the PCA
X_std = StandardScaler().fit_transform(X)
pca = PCA().fit(X_std)

evr = pca.explained_variance_ratio_
components = np.arange(1, len(evr) + 1)
print("Explained variance ratio:", evr.round(3))
print("Cumulative:              ", evr.cumsum().round(3))

# Plot the explained variance ratio and the cumulative explained variance
fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].bar(components, evr)
ax[0].set(xlabel="Component", ylabel="Explained variance ratio", title="Scree plot", xticks=components)
ax[1].plot(components, evr.cumsum(), marker="o")
ax[1].axhline(0.8, color="gray", linestyle="--")
ax[1].set(xlabel="Number of components", ylabel="Cumulative explained variance", title="Cumulative", xticks=components, ylim=(0, 1.05));

import pandas as pd
import seaborn as sns

# Create a DataFrame with the loadings of the first three components
loadings = pd.DataFrame(pca.components_[:3], columns=X.columns, index=["PC1", "PC2", "PC3"])

# Plot the loadings as a heatmap
fig, ax = plt.subplots(figsize=(10, 3))
sns.heatmap(loadings, annot=True, fmt=".2f", cmap="coolwarm", center=0, ax=ax);

# Project the data into the PCA space
X_pca = pca.transform(X_std)

# Plot PC1 against PC2, colored by disease progression
fig, ax = plt.subplots()
points = ax.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap="viridis", alpha=0.7)
ax.set(xlabel="PC1", ylabel="PC2")
fig.colorbar(points, label="Disease progression");

sns.heatmap(X.corr(), annot=True, fmt=".2f", cmap="coolwarm", center=0);

from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split, cross_val_score, KFold

# Split the data and set up the cross-validation
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
cv = KFold(n_splits=10, shuffle=True, random_state=42)

# Loop over the number of components, create the pipeline, evaluate it with cross_val_score()
n_components = np.arange(1, X.shape[1] + 1)
pcr_scores = []

for n in n_components:
    model = make_pipeline(StandardScaler(), PCA(n_components=n), LinearRegression())
    scores = cross_val_score(model, X_train, y_train, cv=cv, scoring="r2")
    pcr_scores.append(scores.mean())

print("PCR CV R²:", np.round(pcr_scores, 3))

from sklearn.cross_decomposition import PLSRegression

# Loop over the number of components, create the PLS model, evaluate it with cross_val_score()
pls_scores = []

for n in n_components:
    pls = PLSRegression(n_components=n)
    scores = cross_val_score(pls, X_train, y_train, cv=cv, scoring="r2")
    pls_scores.append(scores.mean())

print("PLS CV R²:", np.round(pls_scores, 3))

# Compute the baseline
lr_score = cross_val_score(LinearRegression(), X_train, y_train, cv=cv, scoring="r2").mean()
print(f"Linear regression CV R²: {lr_score:.3f}")

# Plot PCR, PLS, and the baseline together
fig, ax = plt.subplots()
ax.plot(n_components, pcr_scores, marker="o", label="PCR")
ax.plot(n_components, pls_scores, marker="s", label="PLS")
ax.axhline(lr_score, color="gray", linestyle="--", label="Linear regression")
ax.set(xlabel="Number of components", ylabel="Cross-validated R²", title="PCR vs PLS", xticks=n_components)
ax.legend();

# Fit the scaler and the PCA on the training data and compute the scores
pca_pipe = make_pipeline(StandardScaler(), PCA()).fit(X_train)
scores_train = pca_pipe.transform(X_train)

# Compute the correlation of each component with y_train
correlations = [np.corrcoef(scores_train[:, i], y_train)[0, 1] for i in range(scores_train.shape[1])]

pd.DataFrame({
    "Explained variance": pca_pipe[-1].explained_variance_ratio_,
    "Correlation with y": correlations,
}, index=[f"PC{i}" for i in n_components]).round(3)

# Fit the final models on the training data and evaluate them on the test set
pcr_final = make_pipeline(StandardScaler(), PCA(n_components=4), LinearRegression()).fit(X_train, y_train)
pls_final = PLSRegression(n_components=2).fit(X_train, y_train)
lr_final = LinearRegression().fit(X_train, y_train)

print(f"Test R² PCR (4 components): {pcr_final.score(X_test, y_test):.3f}")
print(f"Test R² PLS (2 components): {pls_final.score(X_test, y_test):.3f}")
print(f"Test R² linear regression:  {lr_final.score(X_test, y_test):.3f}")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, confusion_matrix, classification_report
from ucimlrepo import fetch_ucirepo

# Fetch dataset 
breast_cancer_wisconsin_diagnostic = fetch_ucirepo(id=17) 
  
# Get data (as pandas dataframes) 
X = breast_cancer_wisconsin_diagnostic.data.features 
y = breast_cancer_wisconsin_diagnostic.data.targets

# Convert y to a 1D array (this is the required input for the logistic regression model)
y = np.ravel(y)

breast_cancer_wisconsin_diagnostic.variables

X.head()

X.describe()

y

# 1. Split the data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)

# 2. Baseline model with three features
features_baseline = ["radius1", "texture1", "area1"]
X_train_base = X_train[features_baseline]
X_test_base = X_test[features_baseline]

baseline_model = LogisticRegression(max_iter=10000)
baseline_model.fit(X_train_base, y_train)

intercept = baseline_model.intercept_
coef = baseline_model.coef_

print("Intercept:", intercept)
print("Coefficients:", coef)

# 3. Evaluate the baseline model
# Accuracy on the test set with baseline features
print("Test accuracy (baseline):", baseline_model.score(X_test_base, y_test))

# Predictions on the test set
y_pred_test_base = baseline_model.predict(X_test_base)

# Confusion matrix on the test set
conf = confusion_matrix(y_test, y_pred_test_base)
print("Confusion matrix (test):\n", conf)

# Classification report on the test set
report = classification_report(y_test, y_pred_test_base, target_names=["Benign", "Malignant"])
print(report)

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# 4. Full model with all features
full_model = Pipeline(steps=[
    ("scaler", StandardScaler()),
    ("logreg", LogisticRegression(max_iter=10000))
])

full_model.fit(X_train, y_train)

# 5. Evaluate the full model
y_pred_test = full_model.predict(X_test)

print("Test accuracy:", accuracy_score(y_test, y_pred_test))
print("Confusion matrix:\n", confusion_matrix(y_test, y_pred_test))
print(classification_report(y_test, y_pred_test, target_names=["Benign", "Malignant"]))

# 6. Plot a custom confusion matrix
conf = confusion_matrix(y_test, y_pred_test)

fig, ax = plt.subplots(figsize=(8, 8))
cax = ax.imshow(conf, cmap='Blues')

# labels
ax.xaxis.set(ticks=(0, 1), ticklabels=('Predicted Benign', 'Predicted Malignant'))
ax.yaxis.set(ticks=(0, 1), ticklabels=('Actual Benign', 'Actual Malignant'))

# Annotate the confusion matrix
for i in range(2):
    for j in range(2):
        ax.text(j, i, conf[i, j], ha='center', va='center', color='red', fontsize=16)

ax.set_title('Confusion Matrix', fontsize=15)

plt.colorbar(cax)
plt.show()

# 6. Use scikit-learn
from sklearn.metrics import ConfusionMatrixDisplay

# Display the confusion matrix using scikit-learn
fig, ax = plt.subplots(figsize=(6, 6))
disp = ConfusionMatrixDisplay(confusion_matrix=conf,
                              display_labels=["Benign", "Malignant"])

disp.plot(cmap='Blues', ax=ax, colorbar=True)
plt.title("Confusion Matrix")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.metrics import classification_report
from sklearn.model_selection import train_test_split
from matplotlib.colors import ListedColormap

# 1. Load data
iris = load_iris()
X = iris.data[:, :2]
y = iris.target
target_names = iris.target_names

# 2. Split into train/test
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=123)

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# 3. TODO: Fit a LDA model and print the classification report
lda = LinearDiscriminantAnalysis()
lda.fit(X_train, y_train)

print(classification_report(y_test, lda.predict(X_test)))

from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis

# 3. TODO: Fit a QDA model and print the classification report
qda = QuadraticDiscriminantAnalysis()
qda.fit(X_train, y_train)

print(classification_report(y_test, qda.predict(X_test)))

from sklearn.naive_bayes import GaussianNB

# 3. TODO: Fit a Gaussian Naive Bayes model and print the classification report
gnb = GaussianNB()
gnb.fit(X_train, y_train)

print(classification_report(y_test, gnb.predict(X_test)))

# 4. Plot the decision boundaries for all 3 classifiers

# Plotting function
def plot_decision_boundary(model, X, y, title, ax):
    h = .02
    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
    xx, yy = np.meshgrid(np.arange(x_min, x_max, h),
                         np.arange(y_min, y_max, h))
    Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)

    cmap_light = ListedColormap(['#FFAAAA', '#AAFFAA', '#AAAAFF'])
    cmap_bold = ListedColormap(['#FF0000', '#00FF00', '#0000FF'])

    ax.contourf(xx, yy, Z, cmap=cmap_light, alpha=0.2)
    scatter = ax.scatter(X[:, 0], X[:, 1], c=y, cmap=cmap_bold, s=30)
    ax.set_xlim(xx.min(), xx.max())
    ax.set_ylim(yy.min(), yy.max())
    ax.set_title(title)
    ax.set_xlabel('Sepal length')
    ax.set_ylabel('Sepal width')

# Create plots for all 3 classifiers
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
plot_decision_boundary(lda, X_train, y_train, "LDA Decision Boundary", axes[0])
plot_decision_boundary(qda, X_train, y_train, "QDA Decision Boundary", axes[1])
plot_decision_boundary(gnb, X_train, y_train, "Naïve Bayes Decision Boundary", axes[2])
plt.tight_layout()

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

wage = sm.datasets.get_rdataset("Wage", "ISLR").data

fig, ax = plt.subplots()
ax.scatter(wage["age"], wage["wage"], alpha=0.15)
ax.set(xlabel="Age", ylabel="Wage", title="Wage vs Age");

wage.head()

from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression

poly = PolynomialFeatures(degree=3, include_bias=False)
X_poly = poly.fit_transform(wage[["age"]])
poly_mod = LinearRegression().fit(X_poly, wage["wage"])

import patsy
import statsmodels.api as sm

cuts = [35, 50, 65]

step_basis = patsy.dmatrix(
    "bs(age, knots=cuts, degree=0, include_intercept=False)",
    {"age": wage["age"], "cuts": cuts},
    return_type="dataframe",
    )
step_mod = sm.OLS(wage["wage"], step_basis).fit()

spline_basis = patsy.dmatrix(
    "bs(age, df=6, degree=3, include_intercept=False)",
    {"age": wage["age"]},
    return_type="dataframe",
    )
spline_mod = sm.OLS(wage["wage"], spline_basis).fit()

from statsmodels.nonparametric.smoothers_lowess import lowess

low_mod = lowess(wage["wage"], wage["age"], frac=0.4)

age_grid = np.linspace(wage["age"].min(), wage["age"].max(), 300)

# Polynomial predictions
poly_pred = poly_mod.predict(poly.transform(age_grid.reshape(-1, 1)))

# Stepwise predictions
step_grid = patsy.dmatrix(
    "bs(x, knots=cuts, degree=0, include_intercept=False)",
    {"x": age_grid, "cuts": cuts},
    return_type="dataframe",
)
step_pred = step_mod.predict(step_grid)

# Spline predictions
spline_grid = patsy.dmatrix(
    "bs(x, df=6, degree=3, include_intercept=False)",
    {"x": age_grid},
    return_type="dataframe",
)
spline_pred = spline_mod.predict(spline_grid)

fig, ax = plt.subplots()
ax.scatter(wage["age"], wage["wage"], alpha=0.15, label="Data")
ax.plot(age_grid, poly_pred, label="Polynomial (deg=3)", linewidth=2)
ax.plot(age_grid, step_pred, label="Step function", linewidth=2)
ax.plot(age_grid, spline_pred, label="Spline (df=6)", linewidth=2)
ax.plot(low_mod[:,0], low_mod[:,1], label="LOWESS (frac=0.4)", linewidth=2)
ax.set(xlabel="Age", ylabel="Wage")
plt.legend();

import matplotlib.pyplot as plt
from sklearn.datasets import load_diabetes
from statsmodels.gam.api import GLMGam, BSplines

# 1. Get bmi as x and the target as y
data = load_diabetes(as_frame=True)
x = data.data[['bmi']]
y = data.target

# 2. Define possible parameters
df_values = [4, 6, 12]
degree_values = [2, 3]

# 3. PLot partial effect for each combination of df and degree
fig, axes = plt.subplots(len(df_values), len(degree_values), figsize=(15, 10), sharey=True)

for i, df_val in enumerate(df_values):
    for j, deg_val in enumerate(degree_values):
        bs = BSplines(x, df=df_val, degree=deg_val)
        gam = GLMGam(y, smoother=bs)
        res = gam.fit()

        res.plot_partial(0, cpr=True, ax=axes[i, j])
        axes[i, j].set_title(f'B-spline: df={df_val}, degree={deg_val}')
        axes[i, j].set_xlabel('BMI')
        axes[i, j].set_ylabel('Effect')

plt.tight_layout()
plt.show()

import pandas as pd
from ISLP import load_data
from statsmodels.gam.api import GLMGam, BSplines

# Load data
Wage = load_data('Wage')

# Continuous features
spline_features = ['age', 'year']
X_spline = Wage[spline_features]

# Categorical features — one-hot encode
categoricals = ['education', 'jobclass', 'maritl']
X_cat = pd.get_dummies(Wage[categoricals], drop_first=True)

# Outcome
y = Wage['wage']

# Create BSpline basis
bs = BSplines(X_spline, df=[6]*len(spline_features), degree=[3]*len(spline_features))

# Fit GAM
gam = GLMGam(y, exog=X_cat, smoother=bs)
res = gam.fit()

print(res.summary())

from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeClassifier, plot_tree
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt

# 1) Load and inspect data
diab = fetch_openml("diabetes", version=1, as_frame=True)
X = diab.data
y = diab.target

print(X.head(), "\n")
print(y.head())

# 2) Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=0)

# 3) Fit tree
clf = DecisionTreeClassifier(max_depth=3, random_state=0)
clf.fit(X_train, y_train)

print("\nTrain accuracy:", accuracy_score(y_train, clf.predict(X_train)))
print("Test accuracy: ", accuracy_score(y_test,  clf.predict(X_test)))

# 5) Plot tree
plt.figure(figsize=(14,7))
plot_tree(clf, feature_names=X.columns, class_names=["neg","pos"], filled=True, rounded=True);

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import accuracy_score, classification_report

# 1) Set up Random Forest + parameter grid
rf = RandomForestClassifier(random_state=0)

param_grid = {
    'n_estimators':      [50, 100, 200],
    'max_depth':         [None, 10, 20],
    'min_samples_split': [2, 5, 10],
    'max_features':      ['sqrt', 'log2', 0.5]
}

# 2) Fit on training data
grid = GridSearchCV(rf, param_grid=param_grid, cv=5, scoring='accuracy', verbose=1)
grid.fit(X_train, y_train)

# 3) Print best hyperparameters
print("Best parameters:", grid.best_params_)
print(f"CV accuracy: {grid.best_score_:.3f}")

# 4) Evaluate on the held‐out test set
best_rf = grid.best_estimator_
y_pred = best_rf.predict(X_test)

print(f"\nTest accuracy: {accuracy_score(y_test, y_pred):.3f}\n")
print("Classification Report:")
print(classification_report(y_test, y_pred, target_names=['neg','pos']))

import numpy as np
import matplotlib.pyplot as plt
from sklearn.ensemble import GradientBoostingRegressor

rng = np.random.RandomState(1)

# Data: sin curve + noise
X = np.linspace(0, 6, 160).reshape(-1, 1)
y = np.sin(X).ravel() + rng.normal(scale=0.25, size=len(X))

# Sklearn model 
gb_sklearn = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1, max_depth=2, random_state=0)
gb_sklearn.fit(X, y)
y_pred_sklearn = gb_sklearn.predict(X)

# Initial plot
fig, ax = plt.subplots()
ax.scatter(X, y, s=25, alpha=0.3, color="black")
ax.plot(X, y_pred_sklearn, label="sklearn model", linewidth=2)
ax.set(xlabel="x", ylabel="y")
ax.legend();

from sklearn.tree import DecisionTreeRegressor

def custom_gb_mse(X, y, M=50, max_depth=2, learning_rate=0.1, random_state=0):
    rng = np.random.RandomState(random_state)
    
    # Initialise the model with a constant predictor.
    f0 = y.mean()
    F = np.full_like(y, fill_value=f0, dtype=float)

    learners = []
    for _ in range(M):
        # Compute the pseudo-residuals for the squared error loss.
        # Hint: This requires you to calculate the negative gradient of the loss with respect to the current prediction F.
        r = y - F

        # Fit a weak learner (shallow regression tree) to the residuals. Each tree learns how to correct the current model's mistakes.
        tree = DecisionTreeRegressor(max_depth=max_depth, random_state=rng.randint(0, 10_000))
        tree.fit(X, r)

        # Update the current model by adding the new tree's predictions, scaled by the learning rate.
        F = F + learning_rate * tree.predict(X)

        learners.append(tree)

    return f0, learners

def predict_gb(X_new, f0, learners, learning_rate):
    # Compute predictions of the full boosted model. Start from the constant baseline f0 and add the contributions.
    pred = np.full(X_new.shape[0], f0, dtype=float)
    
    for tree in learners:
        pred += learning_rate * tree.predict(X_new)
    
    return pred

# Train the model
f0, learners = custom_gb_mse(X, y, M=200, max_depth=2, learning_rate=0.1, random_state=0) 
y_pred_custom = predict_gb(X, f0, learners, learning_rate=0.1)

# Plot comparison with sklearn
fig, ax = plt.subplots()
ax.scatter(X, y, s=25, alpha=0.3, color="black")
ax.plot(X, y_pred_sklearn, label="sklearn model", lw=2)
ax.plot(X, y_pred_custom, label="custom model", lw=2, ls="--")
ax.set(xlabel="x", ylabel="y")
ax.legend();

import seaborn as sns
df = sns.load_dataset("fmri")
df

df_wide = df.pivot_table(
    index=["subject","timepoint","event"],
    columns="region",
    values="signal"
).reset_index()
df_wide.columns.name = None

X = df_wide[["frontal","parietal"]] 
y = df_wide["event"].map({"cue":0,"stim":1})

print("\nFeatures:")
print(X.head())
print("\nTarget:")
print(y.head())

import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.metrics import classification_report

# 1. Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.3, stratify=y, random_state=42)

# 2. Scale the features after splitting (important to avoid data leakage)
scaler = StandardScaler()
X_train_sc = scaler.fit_transform(X_train)
X_test_sc  = scaler.transform(X_test)

# 3. Fit the SVC models and compare the classification reports
clf_lin = SVC(kernel="linear", random_state=42)
clf_lin.fit(X_train_sc, y_train)
y_pred_lin = clf_lin.predict(X_test_sc)
print("Linear SVC\n", classification_report(y_test, y_pred_lin))

clf_rbf = SVC(kernel="rbf", random_state=42)
clf_rbf.fit(X_train_sc, y_train)
y_pred_rbf = clf_rbf.predict(X_test_sc)
print("RBF SVC\n", classification_report(y_test, y_pred_rbf))

import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

def plot_svc_decision_function(model, ax=None):
    """Plot the decision boundary for a trained 2D SVC model."""
    # Set up grid
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    xx, yy = np.meshgrid(np.linspace(*xlim, 100), np.linspace(*ylim, 100))
    grid = np.c_[xx.ravel(), yy.ravel()]
    decision_values = model.decision_function(grid).reshape(xx.shape)
    ax.contour(xx, yy, decision_values, levels=[0], linestyles=['-'], colors='k')

# Plot
fig, ax = plt.subplots(1,2, figsize=(12, 6))

legend_elements = [
    Line2D([0], [0], marker='o', linestyle='None', markersize=8, label='Cue', markerfacecolor="#0173B2", markeredgecolor='None'),
    Line2D([0], [0], marker='o', linestyle='None', markersize=8, label='Stim', markerfacecolor="#DE8F05", markeredgecolor='None'),
    Line2D([0], [0], color='k', linestyle='-', label='Decision boundary')]

# Linear SVC
sns.scatterplot(x = X_train_sc[:, 0], y = X_train_sc[:, 1], hue = y_train.map({0:"cue",1:"stim"}), palette = ["#0173B2", "#DE8F05"], s = 60, ax = ax[0], legend=None)
ax[0].set(xlabel = "Frontal signal (scaled)", ylabel = "Parietal signal (scaled)", title  = "Linear SVC Decision Boundary")
plot_svc_decision_function(clf_lin, ax=ax[0])
ax[0].legend(handles=legend_elements, loc="upper left", handlelength=1)

# RBF SVC
sns.scatterplot(x = X_train_sc[:, 0], y = X_train_sc[:, 1], hue = y_train.map({0:"cue",1:"stim"}), palette = ["#0173B2", "#DE8F05"], s = 60, ax = ax[1], legend=None)
ax[1].set(xlabel = "Frontal signal (scaled)", ylabel = "Parietal signal (scaled)", title  = "RBF SVC Decision Boundary")
plot_svc_decision_function(clf_rbf, ax=ax[1])
ax[1].legend(handles=legend_elements, loc="upper left", handlelength=1);

from sklearn.model_selection import GridSearchCV
from sklearn.svm import SVC
import numpy as np

param_grid = [
    # RBF kernel
    {
        "kernel": ["rbf"],
        "C": np.logspace(-2, 2, 10),
        "gamma": np.logspace(-3, 1, 10)
    },
    # Linear kernel
    {
        "kernel": ["linear"],
        "C": np.logspace(-2, 2, 10)
    },
]

grid = GridSearchCV(
    SVC(),
    param_grid,
    cv=5,
    scoring="accuracy",
    n_jobs=-1
)

grid.fit(X_train_sc, y_train)

print("Best params:", grid.best_params_)
print("CV accuracy:", grid.best_score_)
print("Test accuracy:", grid.score(X_test_sc, y_test))