scikitlearn_cross-validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "25dfd6d1",
   "metadata": {},
   "source": [
    "# 1) Cross-validation\n",
    "\n",
    "author: Mat Gilson, https://github.com/MatthieuGilson\n",
    "\n",
    "This notebook shows a key concept at the core of all machine learning procedures. It aims to quantify the generalizability of a trained classifier to unseen data. Here we also see how to get a baseline reference in terms of accuracy. \n",
    "\n",
    "See also the documentation of scikit-learn library (https://scikit-learn.org/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "707a0699",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import librairies\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.model_selection import ShuffleSplit, StratifiedShuffleSplit, LeaveOneOut\n",
    "\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "\n",
    "font = {'family' : 'DejaVu Sans',\n",
    "        'weight' : 'regular',\n",
    "        'size'   : 18}\n",
    "mpl.rc('font', **font)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "219520bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create synthetic dataset where 2 classes of s0+s1 samples of m-dimensional inputs with controlled contrast\n",
    "def gen_inputs(m,        # input dimensionality\n",
    "               s0,       # number of samples for class 0\n",
    "               s1,       # number of samples for class 1\n",
    "               scaling): # scaling factor to separate classes\n",
    "\n",
    "    # labels\n",
    "    lbl = np.zeros([s0+s1], dtype=int)\n",
    "    # inputs\n",
    "    X = np.zeros([s0+s1,m])\n",
    "\n",
    "    # create s0 and s1 samples for the 2 classes\n",
    "    for i in range(s0+s1):\n",
    "        # label\n",
    "        lbl[i] = int(i<s0)\n",
    "        # inputs are random noise plus a shift\n",
    "        for j in range(m):\n",
    "            # positive/negative shift for 1st/2nd class\n",
    "            if i<s0:\n",
    "                a = -scaling\n",
    "            else:\n",
    "                a = scaling\n",
    "            # the shift linearly depends on the feature index j\n",
    "            X[i,j] = a*j/m + np.random.randn()\n",
    "            \n",
    "    return X, lbl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "506698d5",
   "metadata": {},
   "source": [
    "Let's have a first look at the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bc5284f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# input properties\n",
    "m = 10 # input dimensionality\n",
    "s0 = 100 # number of samples for class 0\n",
    "s1 = 100 # number of samples for class 1\n",
    "scaling = 1.0 # class contrast\n",
    "\n",
    "# generate inputs\n",
    "X, y = gen_inputs(m, s0, s1, scaling)\n",
    "\n",
    "# bins for istograms\n",
    "vbins = np.linspace(-3,3,30)\n",
    "\n",
    "# plot\n",
    "plt.figure(figsize=[6,7])\n",
    "plt.subplot(411)\n",
    "i = 0\n",
    "plt.hist(X[:s0,i], histtype='step', bins=vbins, color='r')\n",
    "plt.hist(X[s0:,i], histtype='step', bins=vbins, color='b')\n",
    "plt.axis(xmin=-3, xmax=3)\n",
    "plt.legend(['class 0', 'class 1'], fontsize=10)\n",
    "plt.title('input {}'.format(i), loc='left')\n",
    "plt.subplot(412)\n",
    "i = int((m-1)*0.33)\n",
    "plt.hist(X[:s0,i], histtype='step', bins=vbins, color='r')\n",
    "plt.hist(X[s0:,i], histtype='step', bins=vbins, color='b')\n",
    "plt.axis(xmin=-3, xmax=3)\n",
    "plt.title('input {}'.format(i), loc='left')\n",
    "plt.subplot(413)\n",
    "i = int((m-1)*0.66)\n",
    "plt.hist(X[:s0,i], histtype='step', bins=vbins, color='r')\n",
    "plt.hist(X[s0:,i], histtype='step', bins=vbins, color='b')\n",
    "plt.axis(xmin=-3, xmax=3)\n",
    "plt.title('input {}'.format(i), loc='left')\n",
    "plt.subplot(414)\n",
    "i = m-1\n",
    "plt.hist(X[:s0,i], histtype='step', bins=vbins, color='r')\n",
    "plt.hist(X[s0:,i], histtype='step', bins=vbins, color='b')\n",
    "plt.axis(xmin=-3, xmax=3)\n",
    "plt.title('input {}'.format(i), loc='left')\n",
    "plt.xlabel('X values')\n",
    "plt.tight_layout()\n",
    "plt.savefig('ex_contrast_X')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b9a9c63",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create matrix rgb\n",
    "mat_rgb = np.zeros([m, vbins.size-1, 3])\n",
    "for i in range(m):\n",
    "    mat_rgb[i,:,0] = np.histogram(X[:s0,i], bins=vbins)[0]\n",
    "    mat_rgb[i,:,2] = np.histogram(X[s0:,i], bins=vbins)[0]\n",
    "mat_rgb /= mat_rgb.max() / 2.0\n",
    "\n",
    "plt.figure(figsize=[6,7])\n",
    "plt.imshow(mat_rgb)\n",
    "plt.xlabel('values')\n",
    "plt.ylabel('input index')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46880c12",
   "metadata": {},
   "source": [
    "Now let's see how to separate the 2 classes using a classifier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d25b05e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Classifiers and learning parameters\n",
    "clf = make_pipeline(StandardScaler(), \n",
    "                    LogisticRegression(C=10000.0, penalty='l2', solver='lbfgs', max_iter=500) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "845c6ded",
   "metadata": {},
   "source": [
    "What is the reference as \"chance\" level: $50\\%$ for $2$ classes?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "832441a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = pd.DataFrame(columns=['score'])\n",
    "\n",
    "# repetitions\n",
    "n_rep = 20\n",
    "for i_rep in range(n_rep):\n",
    "    # generate data\n",
    "    X, y = gen_inputs(m, s0, s1, scaling)\n",
    "    \n",
    "    # Train and test classifiers with subject labels\n",
    "    clf.fit(X, y)\n",
    "    # accuracy on train set\n",
    "    d = {'score': [clf.score(X, y)]}\n",
    "    acc = acc.append(pd.DataFrame(data=d))\n",
    "\n",
    "# plot\n",
    "sb.violinplot(data=acc, y='score', scale='width', palette=['brown']) # cut=0\n",
    "plt.text(0, 1.05, str(acc['score'].mean())[:4], horizontalalignment='center')\n",
    "plt.yticks([0,1])\n",
    "plt.axis(ymax=1.02)\n",
    "plt.xlabel('classifier')\n",
    "plt.ylabel('accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28622c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# theoretical chance level\n",
    "chance_level = 1.0 / 2\n",
    "\n",
    "sb.violinplot(data=acc, y='score', scale='width', palette=['brown']) # cut=0\n",
    "plt.plot([-1,1], [chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.axis(ymax=1.02)\n",
    "plt.xlabel('classifier')\n",
    "plt.ylabel('accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae55e427",
   "metadata": {},
   "source": [
    "We can play with the contrast between the 2 classes, for example a more difficult classification with lower contrast / separability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5452ae09",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = pd.DataFrame(columns=['score'])\n",
    "\n",
    "# change contrast\n",
    "scaling = 0.5 # try 0.2, 0.1, 0.0\n",
    "\n",
    "# repetitions\n",
    "n_rep = 20\n",
    "for i_rep in range(n_rep):\n",
    "    # generate data\n",
    "    X, y = gen_inputs(m, s0, s1, scaling)\n",
    "    \n",
    "    # Train and test classifiers with subject labels\n",
    "    clf.fit(X, y)\n",
    "    # accuracy on train set\n",
    "    d = {'score': [clf.score(X, y)]}\n",
    "    acc = acc.append(pd.DataFrame(data=d))\n",
    "\n",
    "# plot\n",
    "sb.violinplot(data=acc, y='score', scale='width', palette=['brown']) # cut=0\n",
    "plt.text(0, 1.05, str(acc['score'].mean())[:4], horizontalalignment='center')\n",
    "plt.plot([-1,1], [chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.axis(ymax=1.02)\n",
    "plt.xlabel('classifier')\n",
    "plt.ylabel('accuracy')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9d4a9a2",
   "metadata": {},
   "source": [
    "Even for `scaling=0`, the classification accuracy is above the expected chance level $0.5$...\n",
    "\n",
    "Let's try with similar (i.e. with same scaling), but new data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65df5e9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new, y_new = gen_inputs(m, s0, s1, scaling) # also change the scaling to play with the code\n",
    "\n",
    "print(clf.score(X_new, y_new))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cddebf1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = pd.DataFrame(columns=['type','score'])\n",
    "\n",
    "# input dimensionality (number of features)\n",
    "m = 10 # try 5, 20\n",
    "\n",
    "# class contrast\n",
    "scaling = 1.0 # try 0.5, 0.0\n",
    "\n",
    "# loop with training on a dataset and testing on a new dataset\n",
    "for i_rep in range(n_rep):\n",
    "    # generate data\n",
    "    X, y = gen_inputs(m, s0, s1, scaling)\n",
    "\n",
    "    # train and calcluate accuracy\n",
    "    clf.fit(X, y)\n",
    "    d = {'type': ['training'],\n",
    "         'score': [clf.score(X, y)]}\n",
    "    acc = acc.append(pd.DataFrame(data=d))\n",
    "\n",
    "    # generate new data\n",
    "    X_new, y_new = gen_inputs(m, s0, s1, scaling)\n",
    "    \n",
    "    # only test classifier that was trained on other data\n",
    "    d = {'type': ['new'],\n",
    "         'score': [clf.score(X_new, y_new)]}\n",
    "    acc = acc.append(pd.DataFrame(data=d))\n",
    "    \n",
    "sb.violinplot(data=acc, x='type', y='score', scale='width', palette=['brown','orange'])\n",
    "plt.text(0, 1.05, str(acc[acc['type']=='training']['score'].mean())[:4], horizontalalignment='center')\n",
    "plt.text(1, 1.05, str(acc[acc['type']=='new']['score'].mean())[:4], horizontalalignment='center')\n",
    "plt.plot([-1,2], [chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.axis(ymax=1.02)\n",
    "plt.xlabel('classifier')\n",
    "plt.ylabel('accuracy')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "154c6878",
   "metadata": {},
   "source": [
    "The classifier tends to extract specific \"information\" from the data it is trained with, which corresponds to the notion of overfitting.\n",
    "\n",
    "The \"real\" accuracy that should be taken into account is the accuracy for the new data, which quantifies the generalization capability of the classifier to new data from the same class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3783c050",
   "metadata": {},
   "source": [
    "## Cross-validation scheme\n",
    "\n",
    "The idea is to generalize the previous observation by splitting the data into a training set and a testing set, ofr a number of repetitions. The relevant result to report is the test accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5895003c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# number of repetitions and storage of results\n",
    "n_rep = 10\n",
    "\n",
    "# Cross-validation scheme\n",
    "cvs0 = ShuffleSplit(n_splits=n_rep, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d69f322",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate n_rep splits\n",
    "ind_split = np.zeros([n_rep,s0+s1])\n",
    "i_rep = 0\n",
    "for train_ind, test_ind in cvs0.split(X, y):\n",
    "    ind_split[i_rep, test_ind] = 1\n",
    "    i_rep += 1\n",
    "\n",
    "# calculate the size of the test set for each split\n",
    "test_size = np.vstack((ind_split[:,:s0].sum(axis=1),\n",
    "                       ind_split[:,s0:].sum(axis=1)))\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.imshow(ind_split, cmap='binary', interpolation='nearest', aspect=40)\n",
    "plt.xlabel('sample index')\n",
    "plt.ylabel('split index')\n",
    "plt.subplot(122)\n",
    "plt.plot(test_size[0,::-1], np.arange(n_rep), 'b')\n",
    "plt.plot(test_size[1,::-1], np.arange(n_rep), 'r')\n",
    "plt.xlabel('test size per class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bd9e7e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# wrapper to test cross-validation scheme (cvs) \n",
    "def test_cvs(cvs, cvs_lbl):\n",
    "    \n",
    "    acc = pd.DataFrame(columns=['type', 'cv', 'score'])\n",
    "    \n",
    "    for train_ind, test_ind in cvs.split(X, y):\n",
    "    \n",
    "        # train and test classifier\n",
    "        clf.fit(X[train_ind,:], y[train_ind])\n",
    "        # accuracy on train set\n",
    "        d = {'type': ['train'],\n",
    "             'cv': [cvs_lbl], \n",
    "             'score': [clf.score(X[train_ind,:], y[train_ind])]}\n",
    "        acc = pd.concat((acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "        # accuracy on test set\n",
    "        d = {'type': ['test'],\n",
    "             'cv': [cvs_lbl],\n",
    "             'score': [clf.score(X[test_ind,:], y[test_ind])]}\n",
    "        acc = pd.concat((acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "\n",
    "    return acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e594d7e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a dataset of features and labels\n",
    "m = 10 # input dimensionality\n",
    "s0 = 100 # number of samples for class 0\n",
    "s1 = 100 # number of samples for class 1\n",
    "scaling = 1.0 # class contrast\n",
    "X, y = gen_inputs(m, s0, s1, scaling)\n",
    "\n",
    "# evaluate classification accuracy on train and test sets\n",
    "acc = test_cvs(cvs0, 'no strat')\n",
    "\n",
    "# theoretical chance level\n",
    "chance_level = 0.5\n",
    "\n",
    "sb.violinplot(data=acc, x='cv', y='score', hue='type', split=True, scale='width', palette=['brown','orange'])\n",
    "plt.plot([-1,2], [chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy')\n",
    "plt.axis(ymax=1.02)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ea00c4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# matrix of predicted class versus true class in test set\n",
    "cm = np.zeros([2,2])\n",
    "\n",
    "# repeat classification\n",
    "for train_ind, test_ind in cvs0.split(X, y):\n",
    "    clf.fit(X[train_ind,:], y[train_ind])\n",
    "    cm += confusion_matrix(y_true=y[test_ind], \n",
    "                           y_pred=clf.predict(X[test_ind,:]))        \n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(cm, cmap='jet', vmin=0)\n",
    "plt.colorbar(label='sample count')\n",
    "plt.xticks([0,1], ['class 0','class 1'])\n",
    "plt.yticks([0,1], ['class 0','class 1'])\n",
    "plt.xlabel('predicted label')\n",
    "plt.ylabel('true label')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04873ea0",
   "metadata": {},
   "source": [
    "# Baseline accuracy as reference \n",
    "\n",
    "How to interpret the test accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c49c4763",
   "metadata": {},
   "outputs": [],
   "source": [
    "# wrapper to test cross-validation scheme (cvs) \n",
    "def test_cvs(cvs, cvs_y):\n",
    "    \n",
    "    acc = pd.DataFrame(columns=['type', 'cv', 'score'])\n",
    "    \n",
    "    for train_ind, test_ind in cvs.split(X, y):\n",
    "    \n",
    "        # train and test classifiers\n",
    "        clf.fit(X[train_ind,:], y[train_ind])\n",
    "        # accuracy on test set\n",
    "        d = {'type': ['test'],\n",
    "             'cv': [cvs_y],\n",
    "             'score': [clf.score(X[test_ind,:], y[test_ind])]}\n",
    "        acc = pd.concat((acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "\n",
    "        # shuffling\n",
    "        train_ind_rand = np.random.permutation(train_ind)\n",
    "        clf.fit(X[train_ind,:], y[train_ind_rand])\n",
    "        # accuracy on test set\n",
    "        d = {'type': ['shuf'],\n",
    "             'cv': [cvs_y],\n",
    "             'score': [clf.score(X[test_ind,:], y[test_ind])]}\n",
    "        acc = pd.concat((acc, pd.DataFrame(data=d)), ignore_index=True)\n",
    "\n",
    "    return acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e0f1248",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = test_cvs(cvs0, 'no strat')\n",
    "\n",
    "# theoretical chance level\n",
    "chance_level = 0.5\n",
    "\n",
    "sb.violinplot(data=acc, x='cv', y='score', hue='type', split=True, scale='width', palette=['orange','gray'])\n",
    "plt.plot([-1,2], [chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy')\n",
    "plt.axis(ymax=1.02)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "082f8d49",
   "metadata": {},
   "source": [
    "We find distributed values of accuracies for the shuffling surrogates around the expected value. It looks fine, but is all the variability coming from the classifier?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21edac79",
   "metadata": {},
   "source": [
    "# Unbalanced data\n",
    "\n",
    "We now look into a different version of the dataset, a bit less artificial. We consider the case of unbalanced classes, with different number of samples in the two classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "902d7b27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate inputs with unbalanced classes\n",
    "m = 10 # input dimensionality\n",
    "s0 = 100 # number of samples for class 0\n",
    "s1 = 300 # number of samples for class 1\n",
    "scaling = 1.0 # class contrast\n",
    "\n",
    "X, y = gen_inputs(m, s0, s1, scaling) # also change the scaling to play with the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1629a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc = test_cvs(cvs0, 'no strat')\n",
    "\n",
    "# theoretical chance level\n",
    "naive_chance_level = 1.0 / 2 # equal probability for each category\n",
    "greedy_chance_level = max(s0,s1) / (s0+s1) # always predict the larger class\n",
    "\n",
    "sb.violinplot(data=acc, x='cv', y='score', hue='type', split=True, scale='width', palette=['orange','gray'])\n",
    "plt.plot([-1,2], [naive_chance_level]*2, ':k')\n",
    "plt.plot([-1,2], [greedy_chance_level]*2, '--k')\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy')\n",
    "plt.axis(ymax=1.02)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f9a0bd6",
   "metadata": {},
   "source": [
    "The baseline accuracy has changed in this case... A \"stupid\" classifier predicting always the larger class will perform above $50\\%$ accuracy for unbalanced classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b662b586",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate n_rep splits\n",
    "ind_split = np.zeros([n_rep,s0+s1])\n",
    "i_rep = 0\n",
    "for train_ind, test_ind in cvs0.split(X, y):\n",
    "    ind_split[i_rep, test_ind] = 1\n",
    "    i_rep += 1\n",
    "\n",
    "# calculate the size of the test set for each split\n",
    "test_size = np.vstack((ind_split[:,:s0].sum(axis=1),\n",
    "                       ind_split[:,s0:].sum(axis=1)))\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.imshow(ind_split, cmap='binary', interpolation='nearest', aspect=40)\n",
    "plt.xlabel('sample index')\n",
    "plt.ylabel('split index')\n",
    "plt.subplot(122)\n",
    "plt.plot(test_size[0,::-1], np.arange(n_rep), 'b')\n",
    "plt.plot(test_size[1,::-1], np.arange(n_rep), 'r')\n",
    "plt.xlabel('test size per class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a26b002",
   "metadata": {},
   "source": [
    "Does the variability of the class ratio in the test set matter?\n",
    "\n",
    "Let's compare to other possibility for splitting the data in train-test sets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c462266b",
   "metadata": {},
   "source": [
    "## Stratification\n",
    "\n",
    "Let's consider a different splitting scheme that preserves the ratio of classes from the original dataset in each train-test sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19850420",
   "metadata": {},
   "outputs": [],
   "source": [
    "# stratified shuffle split: preserving ratio of classes in train-test sets\n",
    "cvs1 = StratifiedShuffleSplit(n_splits=n_rep, test_size=0.2)\n",
    "\n",
    "# generate n_rep splits\n",
    "ind_split = np.zeros([n_rep,s0+s1])\n",
    "i_rep = 0\n",
    "for train_ind, test_ind in cvs1.split(X, y):\n",
    "    ind_split[i_rep, test_ind] = 1\n",
    "    i_rep += 1\n",
    "\n",
    "# calculate the size of the test set for each split\n",
    "test_size = np.vstack((ind_split[:,:s0].sum(axis=1),\n",
    "                       ind_split[:,s0:].sum(axis=1)))\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(121)\n",
    "plt.imshow(ind_split, cmap='binary', interpolation='nearest', aspect=60)\n",
    "plt.xlabel('sample index')\n",
    "plt.ylabel('split index')\n",
    "plt.subplot(122)\n",
    "plt.plot(test_size[0,::-1], np.arange(n_rep), 'b')\n",
    "plt.plot(test_size[1,::-1], np.arange(n_rep), 'r')\n",
    "plt.xlabel('test size per class')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e9b65a9",
   "metadata": {},
   "source": [
    "## Leave-one-out scheme\n",
    "\n",
    "Another popular scheme consists in leaving out one sample for testing, especially in the case of small datasets. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "005d13cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# leave-two-out scheme\n",
    "cvs2 = LeaveOneOut()\n",
    "\n",
    "# generate n_rep splits\n",
    "ind_split = np.zeros([s0+s1,s0+s1])\n",
    "i_rep = 0\n",
    "for train_ind, test_ind in cvs2.split(X, y):\n",
    "    ind_split[i_rep, test_ind] = 1\n",
    "    i_rep += 1\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(ind_split, cmap='binary', interpolation='nearest', aspect=1)\n",
    "plt.xlabel('sample index')\n",
    "plt.ylabel('split index')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f2ae9a7",
   "metadata": {},
   "source": [
    "## Comparison of cross-validation schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4790171d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# labels for cvs\n",
    "cvs_lbl = ['no strat', 'strat', 'loo']\n",
    "\n",
    "# accuracy\n",
    "acc = pd.DataFrame(columns=['type', 'cv', 'score'])\n",
    "\n",
    "# loop over cross-validation schemes\n",
    "for i, cvs in enumerate([cvs0, cvs1, cvs2]):\n",
    "    acc_tmp = test_cvs(cvs, cvs_lbl[i])\n",
    "    acc = pd.concat((acc, acc_tmp))\n",
    "    \n",
    "# theoretical chance level\n",
    "naive_chance_level = 1.0 / 2 # equal probability for each category\n",
    "greedy_chance_level = max(s0,s1) / (s0+s1) # always predict the larger class\n",
    "\n",
    "# calculate mean of accuracy for leave-two-out scheme\n",
    "sel_df = np.logical_and(acc['cv']=='loo', acc['type']=='shuf')\n",
    "mean_small_test = acc[sel_df]['score'].mean()\n",
    "\n",
    "# plot\n",
    "sb.violinplot(data=acc, x='cv', y='score', hue='type', split=True, scale='width', palette=['orange','gray'], cut=0)\n",
    "plt.plot([-1,3], [naive_chance_level]*2, ':k')\n",
    "plt.plot([-1,3], [greedy_chance_level]*2, '--k')\n",
    "plt.plot([2.5], [mean_small_test], marker='x', markersize=16, color='k', ls='')\n",
    "plt.yticks([0,1])\n",
    "plt.ylabel('accuracy')\n",
    "plt.legend(loc='lower left')\n",
    "plt.tight_layout()\n",
    "plt.savefig('acc_cv')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a5cbd1",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "We see that the leave-one-out (loo) only makes sense when averaging the results (see the cross).\n",
    "\n",
    "The non-stratified splitting inflates the distribution of chance-level accuracy.\n",
    "\n",
    "A good default (for sufficient data) is $80\\%$-$20\\%$ split for the train-test sets, repeated $100$ times.\n",
    "\n",
    "Importantly, this approach is *conservative*: the risk to find \"information\" by chance where there is none (akin to a \"false positive\") can controlled by comparing the classifier with teh chance-level accuracy, in particular its spread. Instead, mistakes tend to yield a large variability in the accuracy distributions, with an absence of positive outcome. In other words, when you find some effect, you can trust it, but it is also likely that you have to try different classifiers to get the best accuracy. This is in fact not a problem, but can inform on the structure of the data (where there is information related to the predicted label).\n",
    "\n",
    "*Refs:*\n",
    "- cross-validation: https://scikit-learn.org/stable/modules/cross_validation.html\n",
    "- chance-level evaluation: https://scikit-learn.org/stable/modules/generated/sklearn.dummy.DummyClassifier.html\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}