<commit_before>"""
==================================================
Statistical comparison of models using grid search
==================================================

This example illustrates how to statistically compare the performance of models
trained and evaluated using :class:`~sklearn.model_selection.GridSearchCV`.
"""

# %%
# We will start by simulating moon shaped data (where the ideal separation
# between classes is non-linear), adding to it a moderate degree of noise.
# Datapoints will belong to one of two possible classes to be predicted by two
# features. We will simulate 50 samples for each class:

print(__doc__)
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_moons

X, y = make_moons(noise=0.352, random_state=1, n_samples=100)

sns.scatterplot(
    x=X[:, 0], y=X[:, 1], hue=y,
    marker='o', s=25, edgecolor='k', legend=False
).set_title("Data")
plt.show()

# %%
# We will compare the performance of :class:`~sklearn.svm.SVC` estimators that
# vary on their `kernel` parameter, to decide which choice of this
# hyper-parameter predicts our simulated data best.
# We will evaluate the performance of the models using
# :class:`~sklearn.model_selection.RepeatedStratifiedKFold`, repeating 10 times
# a 10-fold stratified cross validation using a different randomization of the
# data in each repetition. The performance will be evaluated using
# :class:`~sklearn.metrics.roc_auc_score`.

from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold
from sklearn.svm import SVC

param_grid = [
    {'kernel': ['linear']},
    {'kernel': ['poly'], 'degree': [2, 3]},
    {'kernel': ['rbf']}
]

svc = SVC(random_state=0)

cv = RepeatedStratifiedKFold(
    n_splits=10, n_repeats=10, random_state=0
)

search = GridSearchCV(
    estimator=svc, param_grid=param_grid,
    scoring='roc_auc', cv=cv
)
search.fit(X, y)

# %%
# We can now inspect the results of our search, sorted by their
# `mean_test_score`:

import pandas as pd

results_df = pd.DataFrame(search.cv_results_)
results_df = results_df.sort_values(by=['rank_test_score'])
results_df = (
    results_df
    .set_index(results_df["params"].apply(
        lambda x: "_".join(str(val) for val in x.values()))
    )
    .rename_axis('kernel')
)
results_df[
    ['params', 'rank_test_score', 'mean_test_score', 'std_test_score']
]

# %%
# We can see that the estimator using the `'rbf'` kernel performed best,
# closely followed by `'linear'`. Both estimators with a `'poly'` kernel
# performed worse, with the one using a two-degree polynomial achieving a much
# lower perfomance than all other models.
#
# Usually, the analysis just ends here, but half the story is missing. The
# output of :class:`~sklearn.model_selection.GridSearchCV` does not provide
# information on the certainty of the differences between the models.
# We don't know if these are **statistically** significant.
# To evaluate this, we need to conduct a statistical test.
# Specifically, to contrast the performance of two models we should
# statistically compare their AUC scores. There are 100 samples (AUC
# scores) for each model as we repreated 10 times a 10-fold cross-validation.
#
# However, the scores of the models are not independent: all models are
# evaluated on the **same** 100 partitions, increasing the correlation
# between the performance of the models.
# Since some partitions of the data can make the distinction of the classes
# particularly easy or hard to find for all models, the models scores will
# co-vary.
#
# Let's inspect this partition effect by plotting the performance of all models
# in each fold, and calculating the correlation between models across folds:

# create df of model scores ordered by perfomance
model_scores = results_df.filter(regex=r'split\d*_test_score')

# plot 30 examples of dependency between cv fold and AUC scores
fig, ax = plt.subplots()
sns.lineplot(
    data=model_scores.transpose().iloc[:30],
    dashes=False, palette='Set1', marker='o', alpha=.5, ax=ax
)
ax.set_xlabel("CV test fold", size=12, labelpad=10)
ax.set_ylabel("Model AUC", size=12)
ax.tick_params(bottom=True, labelbottom=False)
plt.show()

# print correlation of AUC scores across folds
print(f"Correlation of models:\n {model_scores.transpose().corr()}")

# %%
# We can observe that the performance of the models highly depends on the fold.
#
# As a consequence, if we assume independence between samples we will be
# underestimating the variance computed in our statistical tests, increasing
# the number of false positive errors (i.e. detecting a significant difference
# between models when such does not exist) [1]_.
#
# Several variance-corrected statistical tests have been developed for these
# cases. In this example we will show how to implement one of them (the so
# called Nadeau and Bengio's corrected t-test) under two different statistical
# frameworks: frequentist and Bayesian.

# %%
# Comparing two models: frequentist approach
# ------------------------------------------
#
# We can start by asking: "Is the first model significantly better than the
# second model (when ranked by `mean_test_score`)?"
#
# To answer this question using a frequentist approach we could
# run a paired t-test and compute the p-value. This is also known as
# Diebold-Mariano test in the forecast literature [5]_.
# Many variants of such a t-test have been developed to account for the
# 'non-independence of samples problem'
# described in the previous section. We will use the one proven to obtain the
# highest replicability scores (which rate how similar the performance of a
# model is when evaluating it on different random partitions of the same
# dataset) while mantaining a low rate of false postitives and false negatives:
# the Nadeau and Bengio's corrected t-test [2]_ that uses a 10 times repeated
# 10-fold cross validation [3]_.
#
# This corrected paired t-test is computed as:
#
# .. math::
#    t=\frac{\frac{1}{k \cdot r}\sum_{i=1}^{k}\sum_{j=1}^{r}x_{ij}}
#    {\sqrt{(\frac{1}{k \cdot r}+\frac{n_{test}}{n_{train}})\hat{\sigma}^2}}
#
# where :math:`k` is the number of folds,
# :math:`r` the number of repetitions in the cross-validation,
# :math:`x` is the difference in performance of the models,
# :math:`n_{test}` is the number of samples used for testing,
# :math:`n_{train}` is the number of samples used for training,
# and :math:`\hat{\sigma}^2` represents the variance of the observed
# differences.
#
# Let's implement a corrected right-tailed paired t-test to evaluate if the
# performance of the first model is significantly better than that of the
# second model. Our null hypothesis is that the second model performs at least
# as good as the first model.

import numpy as np
from scipy.stats import t


def corrected_std(differences, n_train, n_test):
    """Corrects standard deviation using Nadeau and Bengio's approach.

    Parameters
    ----------
    differences : ndarray of shape (n_samples, 1)
        Vector containing the differences in the score metrics of two models.
    n_train : int
        Number of samples in the training set.
    n_test : int
        Number of samples in the testing set.

    Returns
    -------
    corrected_std : int
        Variance-corrected standard deviation of the set of differences.
    """
    n = n_train + n_test
    corrected_var = (
        np.var(differences, ddof=1) * ((1 / n) + (n_test / n_train))
    )
    corrected_std = np.sqrt(corrected_var)
    return corrected_std


def compute_corrected_ttest(differences, df, n_train, n_test):
    """Computes right-tailed paired t-test with corrected variance.

    Parameters
    ----------
    differences : array-like of shape (n_samples, 1)
        Vector containing the differences in the score metrics of two models.
    df : int
        Degrees of freedom.
    n_train : int
        Number of samples in the training set.
    n_test : int
        Number of samples in the testing set.

    Returns
    -------
    t_stat : float
        Variance-corrected t-statistic.<commit_msg>EXA Fix issue in Statistical comparison of models using grid search (#19966)<commit_after>"""
==================================================
Statistical comparison of models using grid search
==================================================

This example illustrates how to statistically compare the performance of models
trained and evaluated using :class:`~sklearn.model_selection.GridSearchCV`.
"""

# %%
# We will start by simulating moon shaped data (where the ideal separation
# between classes is non-linear), adding to it a moderate degree of noise.
# Datapoints will belong to one of two possible classes to be predicted by two
# features. We will simulate 50 samples for each class:

print(__doc__)
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_moons

X, y = make_moons(noise=0.352, random_state=1, n_samples=100)

sns.scatterplot(
    x=X[:, 0], y=X[:, 1], hue=y,
    marker='o', s=25, edgecolor='k', legend=False
).set_title("Data")
plt.show()

# %%
# We will compare the performance of :class:`~sklearn.svm.SVC` estimators that
# vary on their `kernel` parameter, to decide which choice of this
# hyper-parameter predicts our simulated data best.
# We will evaluate the performance of the models using
# :class:`~sklearn.model_selection.RepeatedStratifiedKFold`, repeating 10 times
# a 10-fold stratified cross validation using a different randomization of the
# data in each repetition. The performance will be evaluated using
# :class:`~sklearn.metrics.roc_auc_score`.

from sklearn.model_selection import GridSearchCV, RepeatedStratifiedKFold
from sklearn.svm import SVC

param_grid = [
    {'kernel': ['linear']},
    {'kernel': ['poly'], 'degree': [2, 3]},
    {'kernel': ['rbf']}
]

svc = SVC(random_state=0)

cv = RepeatedStratifiedKFold(
    n_splits=10, n_repeats=10, random_state=0
)

search = GridSearchCV(
    estimator=svc, param_grid=param_grid,
    scoring='roc_auc', cv=cv
)
search.fit(X, y)

# %%
# We can now inspect the results of our search, sorted by their
# `mean_test_score`:

import pandas as pd

results_df = pd.DataFrame(search.cv_results_)
results_df = results_df.sort_values(by=['rank_test_score'])
results_df = (
    results_df
    .set_index(results_df["params"].apply(
        lambda x: "_".join(str(val) for val in x.values()))
    )
    .rename_axis('kernel')
)
results_df[
    ['params', 'rank_test_score', 'mean_test_score', 'std_test_score']
]

# %%
# We can see that the estimator using the `'rbf'` kernel performed best,
# closely followed by `'linear'`. Both estimators with a `'poly'` kernel
# performed worse, with the one using a two-degree polynomial achieving a much
# lower perfomance than all other models.
#
# Usually, the analysis just ends here, but half the story is missing. The
# output of :class:`~sklearn.model_selection.GridSearchCV` does not provide
# information on the certainty of the differences between the models.
# We don't know if these are **statistically** significant.
# To evaluate this, we need to conduct a statistical test.
# Specifically, to contrast the performance of two models we should
# statistically compare their AUC scores. There are 100 samples (AUC
# scores) for each model as we repreated 10 times a 10-fold cross-validation.
#
# However, the scores of the models are not independent: all models are
# evaluated on the **same** 100 partitions, increasing the correlation
# between the performance of the models.
# Since some partitions of the data can make the distinction of the classes
# particularly easy or hard to find for all models, the models scores will
# co-vary.
#
# Let's inspect this partition effect by plotting the performance of all models
# in each fold, and calculating the correlation between models across folds:

# create df of model scores ordered by perfomance
model_scores = results_df.filter(regex=r'split\d*_test_score')

# plot 30 examples of dependency between cv fold and AUC scores
fig, ax = plt.subplots()
sns.lineplot(
    data=model_scores.transpose().iloc[:30],
    dashes=False, palette='Set1', marker='o', alpha=.5, ax=ax
)
ax.set_xlabel("CV test fold", size=12, labelpad=10)
ax.set_ylabel("Model AUC", size=12)
ax.tick_params(bottom=True, labelbottom=False)
plt.show()

# print correlation of AUC scores across folds
print(f"Correlation of models:\n {model_scores.transpose().corr()}")

# %%
# We can observe that the performance of the models highly depends on the fold.
#
# As a consequence, if we assume independence between samples we will be
# underestimating the variance computed in our statistical tests, increasing
# the number of false positive errors (i.e. detecting a significant difference
# between models when such does not exist) [1]_.
#
# Several variance-corrected statistical tests have been developed for these
# cases. In this example we will show how to implement one of them (the so
# called Nadeau and Bengio's corrected t-test) under two different statistical
# frameworks: frequentist and Bayesian.

# %%
# Comparing two models: frequentist approach
# ------------------------------------------
#
# We can start by asking: "Is the first model significantly better than the
# second model (when ranked by `mean_test_score`)?"
#
# To answer this question using a frequentist approach we could
# run a paired t-test and compute the p-value. This is also known as
# Diebold-Mariano test in the forecast literature [5]_.
# Many variants of such a t-test have been developed to account for the
# 'non-independence of samples problem'
# described in the previous section. We will use the one proven to obtain the
# highest replicability scores (which rate how similar the performance of a
# model is when evaluating it on different random partitions of the same
# dataset) while mantaining a low rate of false postitives and false negatives:
# the Nadeau and Bengio's corrected t-test [2]_ that uses a 10 times repeated
# 10-fold cross validation [3]_.
#
# This corrected paired t-test is computed as:
#
# .. math::
#    t=\frac{\frac{1}{k \cdot r}\sum_{i=1}^{k}\sum_{j=1}^{r}x_{ij}}
#    {\sqrt{(\frac{1}{k \cdot r}+\frac{n_{test}}{n_{train}})\hat{\sigma}^2}}
#
# where :math:`k` is the number of folds,
# :math:`r` the number of repetitions in the cross-validation,
# :math:`x` is the difference in performance of the models,
# :math:`n_{test}` is the number of samples used for testing,
# :math:`n_{train}` is the number of samples used for training,
# and :math:`\hat{\sigma}^2` represents the variance of the observed
# differences.
#
# Let's implement a corrected right-tailed paired t-test to evaluate if the
# performance of the first model is significantly better than that of the
# second model. Our null hypothesis is that the second model performs at least
# as good as the first model.

import numpy as np
from scipy.stats import t


def corrected_std(differences, n_train, n_test):
    """Corrects standard deviation using Nadeau and Bengio's approach.

    Parameters
    ----------
    differences : ndarray of shape (n_samples, 1)
        Vector containing the differences in the score metrics of two models.
    n_train : int
        Number of samples in the training set.
    n_test : int
        Number of samples in the testing set.

    Returns
    -------
    corrected_std : int
        Variance-corrected standard deviation of the set of differences.
    """
    kr = len(differences)
    corrected_var = (
        np.var(differences, ddof=1) * (1 / kr + n_test / n_train)
    )
    corrected_std = np.sqrt(corrected_var)
    return corrected_std


def compute_corrected_ttest(differences, df, n_train, n_test):
    """Computes right-tailed paired t-test with corrected variance.

    Parameters
    ----------
    differences : array-like of shape (n_samples, 1)
        Vector containing the differences in the score metrics of two models.
    df : int
        Degrees of freedom.
    n_train : int
        Number of samples in the training set.
    n_test : int
        Number of samples in the testing set.

    Returns
    -------
    t_stat : float
        Variance-corrected t-statistic.
{'fraction_non_alphanumeric': 0.10952469250262686, 'fraction_numerical': 0.00846776685827307, 'mean_word_length': 4.075282308657465, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3948792', 'n_tokens_mistral': 4762, 'n_tokens_neox': 4499, 'n_words': 1923}