{"cells":[{"metadata":{},"cell_type":"markdown","source":"In this notebook I'll make an attempt to anser the following question\n\n# How can we improve the CV estimation of a metric when the train and the test data distributions are not the same?","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"In many cases a metric can be tought as an expected value over the joint distribution of $\\large X, Y$, indicated by $\\large P(X, Y)$\n\n$\\large M = \\int m(x, y) dP(x, y)$\n\nIf we have a finite sample from $\\large P(X, Y)$ indicated as $\\large D = \\{x_i, y_i\\}_{i=1}^N$, we can approximate $\\large M$ using the **Monte Carlo method**\n\n$\\large \\hat{M}^{D}_{MC} = \\frac{1}{N} \\sum_{i=1}^N m(x_i, y_i)$","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# The nail: train and test data distributions are not the same","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"i.e. $\\large P_{train}(X, Y) \\neq P_{test}(X, Y)$\n\nSadly enough the train and the test data distributions are not always the same (thanks Kaggle), so using held out training data can produce a biased estimation of the metric on the test distribution.\n\nCan we do something about it? **Yes!**","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# The hammer 🔨: importance sampling","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We can use a techinque called **importance sampling** to adjust the Monte Carlo estimation with a correction factor given by the density ratio between the two distributions\n\n$\\large \\hat{M}_{IS}^{test} = \\frac{1}{N} \\sum_{i=1}^N m(x_i, y_i) \\frac{p_{test}(x_i, y_i)}{p_{train}(x_i, y_i)}$\n\nwhere $\\large (x_i, y_i)$ are sampled from $\\large P_{train}(X, Y)$\n \nIf $\\large P_{train}(X, Y) > 0$ almost everywhere $\\large P_{test}(X, Y) > 0$ this estimator should be unbiased. In simpler terms, we need to have training examples from all over the test distribution to make a correct estimation.\n\nBut estimating a joint distribution is really hard, especially for non-tabular data... and it is impossible for the test set since we don't have the target variable 😅","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Hypotesys: the conditional distributions of the target variable are the same","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We suppose that, given a point $\\large x$, the conditional distribution of $\\large y$ is the same for training and test data, i.e.\n\n$\\large P_{train}(Y|X) = P_{test}(Y|X)$\n\nThen, using the **chain rule**, we can rewrite the Monte Carlo importance sampling formula as\n\n$\\large \\hat{M}_{IS}^{test} = \\frac{1}{N} \\sum_{i=1}^N m(x_i, y_i) \\frac{p_{test}(x_i)}{p_{train}(x_i)}$\n\nThen all we have to do is to estimate the ratio of densities for $\\large x$, much easier 😉\n\nWe can ironically call this correction factor '*testiness*', since it indicates how much an observation is useful to estimate the real value of the metric on the test distribution and how much more likely it is to be sampled from the test distribution rather then the train one. ","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Estimating the density ratio","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"We can train a learner $\\large L$ to predict wheter a point $\\large x$ is from the train or the test distribution\n\nThen,  appling the **Bayes theorem** on $\\mathbb{P}(Test|x)$, we can prove that\n\n$\\Large \\frac{\\mathbb{P}(x|Test)}{\\mathbb{P}(x|Train)} = \\frac{\\mathbb{P}(Test|x)}{\\mathbb{P}(Train|x)} \\frac{\\mathbb{P}(Train)}{\\mathbb{P}(Test)}$\n\nUsing the trained learner $\\large L$ to estimate $\\large \\mathbb{P}(Test|x) \\approx L(x)$, we obtain\n\n$\\Large \\frac{\\mathbb{P}(x|Test)}{\\mathbb{P}(x|Train)} \\approx \\frac{L(x)}{1 - L(x)} \\frac{ \\mathbb{P}(Train)}{\\mathbb{P}(Test)}$","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Final result","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Replacing the previous expression on the importance sampling formula we obtain\n\n$\\large \\hat{M}_{L}^{test} = \\frac{1}{N} \\sum_{i=1}^N m(x_i, y_i) \\frac{L(x_i) }{1 - L(x_i)} \\frac{\\mathbb{P}(Train)}{\\mathbb{P}(Test)}$\n\nwhere $\\large (x_i, y_i)$ are sampled from the train distribution","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Application to this competition\n\n\nI'll try to use the described technique to improve the CV estimation done in this nice [notebook](https://www.kaggle.com/cdeotte/triple-stratified-kfold-with-tfrecords). The CV AUC is equal to about 0.904, but the LB AUC is equal to 0.9454. This gap is quite big and confusing.\n\nThe AUC can be represented as a double integral multiplied by a constant\n\n$AUC = C \\int \\int \\mathbb{1}[pred(x_0) < pred(x_1)] \\ dp(x_0) \\ dp(x_1)$\n\nso the previous resoning can be similarly applied.\n\nIs important to note that, as pointed [here](https://www.kaggle.com/c/siim-isic-melanoma-classification/discussion/168028), there are regions of the space that contain no train images but test ones, so the importance sampling estimate is still biased, but should be more accurate than the CV one.\n\nTo avoid overfitting the '*testiness*' will be estimated on the hold out data with repeated k-fold splits, using calibrated XGBoost classifiers.\n\nThen, we can simply pass the '*testiness*' as '*sample_weight*' of the chosen sklearn metric (in this case AUC)","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Data reading","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom matplotlib import colors\nfrom sklearn import metrics, datasets, model_selection, linear_model, calibration\nimport xgboost, lightgbm\nfrom tqdm import tqdm\n\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tensorflow.keras import backend as K","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"AUTO = tf.data.experimental.AUTOTUNE\n\ndef read_labeled_tfrecord(example):\n    tfrec_format = {\n        'image'                        : tf.io.FixedLenFeature([], tf.string),\n        'image_name'                   : tf.io.FixedLenFeature([], tf.string),\n        'patient_id'                   : tf.io.FixedLenFeature([], tf.int64),\n        'sex'                          : tf.io.FixedLenFeature([], tf.int64),\n        'age_approx'                   : tf.io.FixedLenFeature([], tf.int64),\n        'anatom_site_general_challenge': tf.io.FixedLenFeature([], tf.int64),\n        'diagnosis'                    : tf.io.FixedLenFeature([], tf.int64),\n        'target'                       : tf.io.FixedLenFeature([], tf.int64)\n    }           \n    example = tf.io.parse_single_example(example, tfrec_format)\n    return example['image'], example['target']\n\n\ndef read_unlabeled_tfrecord(example):\n    tfrec_format = {\n        'image'                        : tf.io.FixedLenFeature([], tf.string),\n        'sex'                          : tf.io.FixedLenFeature([], tf.int64),\n        'age_approx'                   : tf.io.FixedLenFeature([], tf.int64),\n        'anatom_site_general_challenge': tf.io.FixedLenFeature([], tf.int64),\n        'image_name'                   : tf.io.FixedLenFeature([], tf.string),\n    }\n    example = tf.io.parse_single_example(example, tfrec_format)\n    meta = tf.cast(tf.stack([example['sex'], example['age_approx']], 0), tf.float32)\n    return example['image'], meta, example['image_name']\n\n \ndef prepare_image(img, dim=256):    \n    img = tf.image.decode_jpeg(img, channels=3)\n    img = tf.cast(img, tf.float32) / 255.0                      \n    img = tf.reshape(img, [dim,dim, 3])\n            \n    return img\n\n\ndef count_data_items(filenames):\n    n = [int(re.compile(r\"-([0-9]*)\\.\").search(filename).group(1)) \n         for filename in filenames]\n    return np.sum(n)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files = tf.io.gfile.glob(\"../input/melanoma-256x256/train*.tfrec\")\ntest_files = tf.io.gfile.glob(\"../input/melanoma-256x256/test*.tfrec\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"oof = pd.read_csv(\"../input/triple-stratified-kfold-with-tfrecords/oof.csv\") # the kernel out of fold prediction\nbase_train_data = pd.read_csv(\"../input/melanoma-256x256/train.csv\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Since we use XGBoost to estimate the '*testiness*', we need to preprocess the images.\n\nWe will use metadata (age and sex) and a few statistics of the images\n* average and standard deviation\n* channel-wise average and standard deviation\n\nThe target of our classification task will simply be $0$ for images from the train data and $1$ for the ones from the test set.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"ds = tf.data.TFRecordDataset(train_files).map(read_unlabeled_tfrecord, num_parallel_calls=AUTO)\nds = ds.map(lambda i, m, n: (prepare_image(i), m, n), num_parallel_calls=AUTO)\nds = ds.batch(128)\nds = ds.map(lambda i, m, n: (K.concatenate([m, K.reshape(K.mean(i, [1, 2, 3]), [-1, 1]), K.mean(i, [1, 2]), K.std(i, [1, 2])], axis=1), n), num_parallel_calls=AUTO)\nds = ds.prefetch(AUTO)\ndata_train = np.array([i for i, n in tqdm(ds.unbatch())])\nname_train = np.array([n.numpy().decode('utf-8') for i, n in tqdm(ds.unbatch())])\n\nds = tf.data.TFRecordDataset(test_files).map(read_unlabeled_tfrecord, num_parallel_calls=AUTO)\nds = ds.map(lambda i, m, n: (prepare_image(i), m, n), num_parallel_calls=AUTO)\nds = ds.batch(128)\nds = ds.map(lambda i, m, n: (K.concatenate([m, K.reshape(K.mean(i, [1, 2, 3]), [-1, 1]), K.mean(i, [1, 2]), K.std(i, [1, 2])]), n), num_parallel_calls=AUTO)\nds = ds.prefetch(AUTO)\ndata_test = np.array([i for i, n in tqdm(ds.unbatch())])\nname_test = np.array([n.numpy().decode('utf-8') for i, n in tqdm(ds.unbatch())])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"y = np.concatenate([np.zeros(data_train.shape[0]), np.ones(data_test.shape[0])], 0)\nx = np.concatenate([data_train, data_test])\nname = np.concatenate([name_train, name_test])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Minimal EDA","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# train-test data distribution based on features (0: sex, 1: age, 2: img-avg, 3: img-std, 4..6: cw img-avg, 7..9: cw img-std)\ni, j = 2, 3\nbins = 50\nplt.figure(figsize=(20, 5))\nplt.subplot(1, 3, 1)\nc0, a0, b0, _ = plt.hist2d(x[y == 0, i], x[y == 0, j], bins=bins, density=True)\nplt.pcolormesh(a0, b0, np.log(0 + c0).T);\nplt.title(\"Training - Log density\"); plt.xlabel(\"Color avg\"); plt.ylabel(\"Color std\")\nplt.subplot(1, 3, 2)\nplt.title(\"Test - Log density\"); plt.xlabel(\"Color avg\"); plt.ylabel(\"Color std\")\nc1, a1, b1, _ = plt.hist2d(x[y == 1, i], x[y == 1, j], bins=bins, density=True)\nplt.pcolormesh(a1, b1, np.log(0 + c1).T);\nplt.subplot(1, 3, 3)\nplt.title(\"Abs diff\"); plt.xlabel(\"Color avg\"); plt.ylabel(\"Color std\")\nplt.pcolormesh(a0, b0, np.abs(c0 - c1).T);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We see that is sufficient to consider the average and the standard deviation of the images to see a clear difference between the train and the test distributions","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Example of training","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"idxT, idxV = list(model_selection.KFold(5, shuffle=True).split(x, y))[0]\nmodel = xgboost.XGBClassifier()\nmodel.fit(x[idxT], y[idxT], eval_set=[(x[idxT], y[idxT]), (x[idxV], y[idxV])], eval_metric='auc', verbose=0)\nplt.plot(model.evals_result_['validation_0']['auc'], color='tab:blue')\nplt.plot(model.evals_result_['validation_1']['auc'], color='tab:blue', linestyle=\":\")\nplt.xlabel('n_estimators')\nplt.ylabel('AUC');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we see that it is difficult to predict wether an image is from training or test set, but it is clearly better than random guessing","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Calibration example","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"calibrated_model = calibration.CalibratedClassifierCV(model, cv='prefit')\nxcT, xcV, ycT, ycV = model_selection.train_test_split(x[idxV], y[idxV])\ncalibrated_model.fit(xcT, ycT)\nplt.plot(*calibration.calibration_curve(ycV, model.predict_proba(xcV)[:, 1], n_bins=10)[::-1], color='tab:blue', label='XGB')\nplt.plot(*calibration.calibration_curve(ycV, calibrated_model.predict_proba(xcV)[:, 1], n_bins=10)[::-1], color='tab:red', label='calibrated')\nplt.plot(plt.xlim(), plt.ylim(), color='gray', linestyle='--', label='perfect')\nplt.xlabel('Estimated probability'); plt.ylabel('True probability');\nplt.legend();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In this case the calibration seems not to be very useful, so we will skip it for this time","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### Training","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"NREPS = 10\nNFOLDS = 10\ndf_all = []\nmodels = NFOLDS * [xgboost.XGBClassifier()]\nfor rep in range(NREPS):\n    pid = base_train_data[['patient_id']].drop_duplicates().reset_index(drop=True)\n    folds = list(model_selection.KFold(NFOLDS, shuffle=True).split(pid))\n    df_rep = []\n    for fold, (pidIdxT, pidIdxV) in enumerate(folds):\n        pidT, pidV = pid.iloc[pidIdxT]['patient_id'], pid.iloc[pidIdxV]['patient_id']\n        inT = base_train_data[base_train_data['patient_id'].isin(pidT)]['image_name']\n        inV = base_train_data[base_train_data['patient_id'].isin(pidV)]['image_name']\n        idxT = (pd.Series(name).isin(inT) | y == 1)\n        idxV = (pd.Series(name).isin(inV))\n        model = models[fold]\n        model.fit(x[idxT], y[idxT])\n        # calibration should go here\n        pred = model.predict_proba(x[idxV])[:, 1]\n        testiness = pred / (1 - pred) * np.mean(y[idxT] == 0) / np.mean(y[idxT] == 1)\n        df = pd.DataFrame({'rep': rep, 'image_name': name[idxV], 'is_test': y[idxV], 'is_test_pred': pred, 'testiness': testiness})\n        df_rep.append(df)\n        pass\n    df_rep = pd.concat(df_rep, axis=0)\n    df = pd.merge(df_rep, oof, on='image_name')\n    score = metrics.roc_auc_score(df.target, df.pred)\n    corrected_score = metrics.roc_auc_score(df.target, df.pred, sample_weight=df.testiness)\n    print(\"Rep%2d \\t\\tAUC  = %.6f \\tcAUC = %.6f\" % (rep, score, corrected_score))\n    df_all.append(df_rep)\n    pass\ndf_all = pd.concat(df_all, axis=0)\ndf = pd.merge(df_all, oof, on='image_name').groupby('image_name').mean()\nscore = metrics.roc_auc_score(df.target, df.pred)\ncorrected_score = metrics.roc_auc_score(df.target, df.pred, sample_weight=df.testiness)\nprint(\"Overall \\tAUC  = %.6f \\tcAUC = %.6f\" % (score, corrected_score))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"XX = pd.DataFrame({'y': y, 'x2': x[:, 2], 'x3': x[:, 3], 'image_name': name}).merge(df, on='image_name')\nplt.figure(figsize=(15, 5))\nplt.subplot(1, 3, 1)\nc0, a0, b0, _ = plt.hist2d(XX.x2[XX.y == 0], XX.x3[XX.y == 0], 50, density=True);\nplt.subplot(1, 3, 2)\nc1, a1, b1, _ = plt.hist2d(XX.x2[XX.y == 0], XX.x3[XX.y == 0], 50, density=True, weights=XX.testiness[XX.y == 0]);\nplt.subplot(1, 3, 3)\nctest, atest, btest, _ = plt.hist2d(x[y == 1, 2], x[y == 1, 3], 50, density=True);\nkl_pre = keras.losses.kld((c0 / c0.sum()).reshape(-1), (ctest / ctest.sum()).reshape(-1)).numpy()\nkl_post = keras.losses.kld((c1 / c1.sum()).reshape(-1), (ctest / ctest.sum()).reshape(-1)).numpy()\nprint(\"KL divergence pre  = %.4f\" % kl_pre)\nprint(\"KL divergence post = %.4f\" % kl_post)\nprint(\"Difference = %.4f%%\" % (100 * (kl_post - kl_pre) / kl_pre))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_all[['image_name', 'testiness']].groupby('image_name').mean().to_csv(\"../working/testiness.csv\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"# Next steps\n\n* try a better model to estimate the '*testiness*' 😋\n* is it possible to further improve the CV estimation in the case of different support for the two distributions?\n* find a strategy to make calibration work\n* try to use a similar correction directly inside the loss function of a model","execution_count":null}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}