{"cells":[{"metadata":{"_uuid":"6733556261bf55129c54a6fb610849bea0542082"},"cell_type":"markdown","source":"# Flavours of Physics: τ→ 3μ Kaggle Competition Kernel"},{"metadata":{"_uuid":"61338ce091dde757456f982ba09fcadcb91ebf00"},"cell_type":"markdown","source":"This is a kernel for the above competition. We start by copying the evaluation.py file and defining some useful functions (I would like to eventually put these in a separate file and import them). We then select training variables from the dataset, before training several models (with varying hyperparameters) and selecting the best one (using cross-validation evaluation) for our predictions submission."},{"metadata":{"trusted":false,"_uuid":"6a3308cc5767ad8f322ec47b36c90d703483ac51"},"cell_type":"code","source":"version = '8.0'\nfolder = '../input/'","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ab256de5145a3b960c2ef1d53cf621a6780d7499"},"cell_type":"markdown","source":"## SET UP ------------------------------------------------------------------------------------"},{"metadata":{"trusted":false,"_uuid":"6e9fe1a8ae050fc821702a93f96475be4fc7e37b"},"cell_type":"code","source":"### THIS CELL IS JUST THE EVALUATION PYTHON FILE \n\nimport numpy\nfrom sklearn.metrics import roc_curve, auc\n\n\ndef __rolling_window(data, window_size):\n    \"\"\"\n    Rolling window: take window with definite size through the array\n\n    :param data: array-like\n    :param window_size: size\n    :return: the sequence of windows\n\n    Example: data = array(1, 2, 3, 4, 5, 6), window_size = 4\n        Then this function return array(array(1, 2, 3, 4), array(2, 3, 4, 5), array(3, 4, 5, 6))\n    \"\"\"\n    shape = data.shape[:-1] + (data.shape[-1] - window_size + 1, window_size)\n    strides = data.strides + (data.strides[-1],)\n    return numpy.lib.stride_tricks.as_strided(data, shape=shape, strides=strides)\n\n\ndef __cvm(subindices, total_events):\n    \"\"\"\n    Compute Cramer-von Mises metric.\n    Compared two distributions, where first is subset of second one.\n    Assuming that second is ordered by ascending\n\n    :param subindices: indices of events which will be associated with the first distribution\n    :param total_events: count of events in the second distribution\n    :return: cvm metric\n    \"\"\"\n    target_distribution = numpy.arange(1, total_events + 1, dtype='float') / total_events\n    subarray_distribution = numpy.cumsum(numpy.bincount(subindices, minlength=total_events), dtype='float')\n    subarray_distribution /= 1.0 * subarray_distribution[-1]\n    return numpy.mean((target_distribution - subarray_distribution) ** 2)\n\n\ndef compute_cvm(predictions, masses, n_neighbours=200, step=50):\n    \"\"\"\n    Computing Cramer-von Mises (cvm) metric on background events: take average of cvms calculated for each mass bin.\n    In each mass bin global prediction's cdf is compared to prediction's cdf in mass bin.\n\n    :param predictions: array-like, predictions\n    :param masses: array-like, in case of Kaggle tau23mu this is reconstructed mass\n    :param n_neighbours: count of neighbours for event to define mass bin\n    :param step: step through sorted mass-array to define next center of bin\n    :return: average cvm value\n    \"\"\"\n    predictions = numpy.array(predictions)\n    masses = numpy.array(masses)\n    assert len(predictions) == len(masses)\n\n    # First, reorder by masses\n    predictions = predictions[numpy.argsort(masses)]\n\n    # Second, replace probabilities with order of probability among other events\n    predictions = numpy.argsort(numpy.argsort(predictions, kind='mergesort'), kind='mergesort')\n\n    # Now, each window forms a group, and we can compute contribution of each group to CvM\n    cvms = []\n    for window in __rolling_window(predictions, window_size=n_neighbours)[::step]:\n        cvms.append(__cvm(subindices=window, total_events=len(predictions)))\n    return numpy.mean(cvms)\n\n\ndef __roc_curve_splitted(data_zero, data_one, sample_weights_zero, sample_weights_one):\n    \"\"\"\n    Compute roc curve\n\n    :param data_zero: 0-labeled data\n    :param data_one:  1-labeled data\n    :param sample_weights_zero: weights for 0-labeled data\n    :param sample_weights_one:  weights for 1-labeled data\n    :return: roc curve\n    \"\"\"\n    labels = [0] * len(data_zero) + [1] * len(data_one)\n    weights = numpy.concatenate([sample_weights_zero, sample_weights_one])\n    data_all = numpy.concatenate([data_zero, data_one])\n    fpr, tpr, _ = roc_curve(labels, data_all, sample_weight=weights)\n    return fpr, tpr\n\n\ndef compute_ks(data_prediction, mc_prediction, weights_data, weights_mc):\n    \"\"\"\n    Compute Kolmogorov-Smirnov (ks) distance between real data predictions cdf and Monte Carlo one.\n\n    :param data_prediction: array-like, real data predictions\n    :param mc_prediction: array-like, Monte Carlo data predictions\n    :param weights_data: array-like, real data weights\n    :param weights_mc: array-like, Monte Carlo weights\n    :return: ks value\n    \"\"\"\n    assert len(data_prediction) == len(weights_data), 'Data length and weight one must be the same'\n    assert len(mc_prediction) == len(weights_mc), 'Data length and weight one must be the same'\n\n    data_prediction, mc_prediction = numpy.array(data_prediction), numpy.array(mc_prediction)\n    weights_data, weights_mc = numpy.array(weights_data), numpy.array(weights_mc)\n\n    assert numpy.all(data_prediction >= 0.) and numpy.all(data_prediction <= 1.), 'Data predictions are out of range [0, 1]'\n    assert numpy.all(mc_prediction >= 0.) and numpy.all(mc_prediction <= 1.), 'MC predictions are out of range [0, 1]'\n\n    weights_data /= numpy.sum(weights_data)\n    weights_mc /= numpy.sum(weights_mc)\n\n    fpr, tpr = __roc_curve_splitted(data_prediction, mc_prediction, weights_data, weights_mc)\n\n    Dnm = numpy.max(numpy.abs(fpr - tpr))\n    return Dnm\n\n\ndef roc_auc_truncated(labels, predictions, tpr_thresholds=(0.2, 0.4, 0.6, 0.8),\n                      roc_weights=(4, 3, 2, 1, 0)):\n    \"\"\"\n    Compute weighted area under ROC curve.\n\n    :param labels: array-like, true labels\n    :param predictions: array-like, predictions\n    :param tpr_thresholds: array-like, true positive rate thresholds delimiting the ROC segments\n    :param roc_weights: array-like, weights for true positive rate segments\n    :return: weighted AUC\n    \"\"\"\n    assert numpy.all(predictions >= 0.) and numpy.all(predictions <= 1.), 'Data predictions are out of range [0, 1]'\n    assert len(tpr_thresholds) + 1 == len(roc_weights), 'Incompatible lengths of thresholds and weights'\n    fpr, tpr, _ = roc_curve(labels, predictions)\n    area = 0.\n    tpr_thresholds = [0.] + list(tpr_thresholds) + [1.]\n    for index in range(1, len(tpr_thresholds)):\n        tpr_cut = numpy.minimum(tpr, tpr_thresholds[index])\n        tpr_previous = numpy.minimum(tpr, tpr_thresholds[index - 1])\n        area += roc_weights[index - 1] * (auc(fpr, tpr_cut, reorder=True) - auc(fpr, tpr_previous, reorder=True))\n    tpr_thresholds = numpy.array(tpr_thresholds)\n    # roc auc normalization to be 1 for an ideal classifier\n    area /= numpy.sum((tpr_thresholds[1:] - tpr_thresholds[:-1]) * numpy.array(roc_weights))\n    return area","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8ff1f58debbe3b7a08316f5d0f3875223bed1099"},"cell_type":"markdown","source":"### Check agreement test"},{"metadata":{"trusted":false,"_uuid":"4aa31720a4b14f7744bd04237eed8487559a51f2"},"cell_type":"code","source":"def check_ag_test(model,var):\n    check_agreement = pd.read_csv(folder + 'check_agreement.csv', index_col='id')\n    agreement_probs = model.predict_proba(check_agreement[var])[:, 1]\n    \n    ks = compute_ks(\n        agreement_probs[check_agreement['signal'].values == 0],\n        agreement_probs[check_agreement['signal'].values == 1],\n        check_agreement[check_agreement['signal'] == 0]['weight'].values,\n        check_agreement[check_agreement['signal'] == 1]['weight'].values)\n    print('KS metric', ks, ks < 0.09)\n    return ks<0.09\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e64e77b1a687bea3ea4d2602047f3231cd4b1a72"},"cell_type":"markdown","source":"### Check correlation test"},{"metadata":{"trusted":false,"_uuid":"ec6e6e2c15a9c22ee76106cce94d2f3303835085"},"cell_type":"code","source":"def check_corr_test(model,var):\n    \n\n    check_correlation = pd.read_csv(folder + 'check_correlation.csv', index_col='id')\n    correlation_probs = model.predict_proba(check_correlation[var])[:, 1]\n    cvm = compute_cvm(correlation_probs, check_correlation['mass'])\n    print('CvM metric', cvm, cvm < 0.002)\n    return cvm<0.002\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"35e832f5159c144f5370b2b0e4674a550ca8a37a"},"cell_type":"markdown","source":"### Compute weighted AUC on the training data with min_ANNmuon > 0.4 "},{"metadata":{"trusted":false,"_uuid":"2c40baab8b27c911a785fe781484c5c341d7d9d3"},"cell_type":"code","source":"from sklearn.model_selection import cross_val_score\n\ndef comp_auc(model,var,data):\n    train_eval = data[data['min_ANNmuon'] > 0.4]\n    train_probs = model.predict_proba(train_eval[var])[:, 1]\n    AUC = roc_auc_truncated(train_eval['signal'], train_probs)\n    print('AUC', AUC)\n    return AUC\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b4b1044dcd91daa9aaf6390ec09c8d4404fc71a"},"cell_type":"markdown","source":"### Predict test and create file"},{"metadata":{"trusted":false,"_uuid":"f447fc967559ea85aadba8bd65e566a03d302d6d"},"cell_type":"code","source":"def pred_file(model,var):\n\n    test = pd.read_csv(folder + 'test.csv', index_col='id')\n    \n    result = pd.DataFrame({'id': test.index})\n    result['prediction'] = model.predict_proba(test[var])[:, 1]\n    result.to_csv('prediction %s .csv' % version, index=False, sep=',')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"48f78668b20db097c477c74a22deff0f6b173dc8"},"cell_type":"markdown","source":"## END OF SETUP--------------------------------------------------------------------------------------"},{"metadata":{"_uuid":"65380b224f78f50167ae46690909ba64b87d06e3"},"cell_type":"markdown","source":" "},{"metadata":{"_uuid":"ef75eb5603ee3a93a0f12fb12e2c12bd37d33df7"},"cell_type":"markdown","source":"# Imports"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"ffa76c653b77ef166f641c6198e528f927ea598f"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport sklearn\nimport seaborn as sns","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"dc57dc7db6cb06edf80ad9b292061fc05c5fae91"},"cell_type":"markdown","source":"# Read training data"},{"metadata":{"trusted":false,"_uuid":"5002dc3f80c52bf5807b9f6d72fa734f2541f72f"},"cell_type":"code","source":"train = pd.read_csv(folder + 'training.csv', index_col='id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"7e433a3092bafcb1cd1adc1c98fca8f5e2deffd3"},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"31a9ccc1804154396355bef057b36d920104e1bb"},"cell_type":"code","source":"train.info() # to check for missing values","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6de3a421b89221bd9357e6298191f6ad0fb26733"},"cell_type":"markdown","source":"# Define training features"},{"metadata":{"_uuid":"d80f66b6f134f01b7d32110632b8d48181457a18"},"cell_type":"markdown","source":"\nLets plot a seaborn heatmap to see the correlation of each variable to the signal:"},{"metadata":{"trusted":false,"_uuid":"246df8a70fbb59b5bb79b2c5f6a063eb4ad4de65"},"cell_type":"code","source":"plt.figure(figsize=(5,20))\nsns.heatmap(train.corr()[\"signal\"].to_frame().sort_values(by=\"signal\", ascending=False), annot=True, center=0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"158edafd3b7623c2dc8a4ee5e8937c838d1b9425"},"cell_type":"code","source":"# these are the variables that we want to include \n#(mainly based on common sense and trial and error)\n\nvariables = train.drop([\"production\", \"min_ANNmuon\",\"signal\",\"mass\", # these are not to be included\n                        \"SPDhits\", # including this makes agreement test fail\n                        \"FlightDistanceError\" # this seems to worsen score - perhaps not relevant (noise)\n                       ],axis=1).columns\nvariables","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9e0003f8746d6daab0226e5ed27f6467b19f5fa3"},"cell_type":"markdown","source":"# Let's train and compare some models!"},{"metadata":{"trusted":false,"_uuid":"d9e14e3df24dec3a4370f4e8db8ff7cfef20b5aa"},"cell_type":"code","source":"candidate_models = {}   # we'll store candidate models here\n\n#  split train dataset into 4 subsets for cross-validation\nfrom sklearn.utils import shuffle\ns_train = shuffle(train)\nl = len(s_train)\n\nl_1 = int(l/4)\nl_2 = int(l/2)\nl_3 = int(3*l/4)\n\nind_1=s_train.index[[i for i in range(l_1)]]\nind_2=s_train.index[[i for i in range(l_1,l_2)]]\nind_3 = s_train.index[[i for i in range(l_2,l_3)]]\nind_4 = s_train.index[[i for i in range(l_3,len(s_train))]]\n\nsig_1 = s_train['signal'].drop(ind_1)\nsig_2 = s_train['signal'].drop(ind_2)\nsig_3 = s_train['signal'].drop(ind_3)\nsig_4 = s_train['signal'].drop(ind_4)\n\nvar_1 = s_train[variables].drop(ind_1)\nvar_2 = s_train[variables].drop(ind_2)\nvar_3 = s_train[variables].drop(ind_3)\nvar_4 = s_train[variables].drop(ind_4)\n\ntrain_1 = s_train[:l_1]\ntrain_2 = s_train[l_1:l_2]\ntrain_3 = s_train[l_2:l_3]\ntrain_4 = s_train[l_3:]\n\n\ndef test_model(model):\n    #if the model passes the tests...\n    model.fit(train[variables], train['signal'])\n    if(check_corr_test(model,variables) and check_ag_test(model,variables)):\n       \n        # evaluate the model on the 4 subsets\n        model.fit(var_1,sig_1)\n        val_1 = comp_auc(model,variables,train_1) \n        model.fit(var_2,sig_2)\n        val_2 = comp_auc(model,variables,train_2) \n        model.fit(var_3,sig_3)\n        val_3 = comp_auc(model,variables,train_3) \n        model.fit(var_4,sig_4)\n        val_4 = comp_auc(model,variables,train_4)\n        \n        val =(val_1+val_2+val_3+val_4)/4\n        \n        print(\"Average AUC is: \" + str(val))\n                        \n            \n        model.fit(train[variables], train['signal'])\n        #...add the model trained on all the data to the candidates\n        candidate_models[model] = val\n        print('passed')\n    else:\n        print('failed')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af3c18a0eb85d2e905194fb56e72b0e5e37d3f52"},"cell_type":"markdown","source":"### Gradient Boosting"},{"metadata":{"_uuid":"6a5d9b3f805e0cf36820907957484127b23cb80e"},"cell_type":"markdown","source":"This seemed to be the best ML model, so we try it with many different hyperparameters:"},{"metadata":{"trusted":false,"_uuid":"ff797fc4bfb85958d0beb533823abd2955561677"},"cell_type":"code","source":"from sklearn.ensemble import GradientBoostingClassifier\n\ntest_model(GradientBoostingClassifier())\n\nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.2, n_estimators=200, \n                                max_depth=5))\n\nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.1, n_estimators=300, \n                                max_depth=10,max_features = 10))\n           \nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.2, n_estimators=200, \n                                max_depth=15))\n\nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.1, n_estimators=200, \n                                max_depth=6,max_features = 6))\n\nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.05, n_estimators=300, \n                                max_depth=4))\n\nprint('-----')\ntest_model(GradientBoostingClassifier(learning_rate=0.3, n_estimators=200, \n                                max_depth=6,                                  \n                                 warm_start=True))\n\n","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b2bbfaab200d766691026b7484dbacf5b8fc5f1"},"cell_type":"markdown","source":"### Logistic Regression"},{"metadata":{"trusted":false,"_uuid":"7b37ab529260c6554f90b67fd20d11434e3ab2fd"},"cell_type":"code","source":"from sklearn.linear_model import LogisticRegression\n\nfor c in range(0,5):\n    print('-----')\n    test_model(LogisticRegression(C=0.8 + c*0.2))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"44b8cdb0d60a75f0e1ec951e467c94874b1d1dc5"},"cell_type":"markdown","source":"### Gaussian Naive Bayes"},{"metadata":{"trusted":false,"_uuid":"be8bc54a3a14dd75c13daeb99deeb749f31262c9"},"cell_type":"code","source":"from sklearn.naive_bayes import GaussianNB\n\ntest_model(GaussianNB())","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6e1a6dd3c806f254678299291205c74026fc195d"},"cell_type":"markdown","source":"### K Neighbors"},{"metadata":{"trusted":false,"_uuid":"50d9ac6c68bd054d1d2f78ddb9f7a80f8eec33e0"},"cell_type":"code","source":"from sklearn.neighbors import KNeighborsClassifier\n\n\ntest_model(KNeighborsClassifier())\n           \nprint('-----')\ntest_model(KNeighborsClassifier(n_neighbors=8,  leaf_size=20))           ","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"441af1fc7f3e073a5636dfc48fe5fb58ea0a8bdb"},"cell_type":"markdown","source":"### Decision Tree"},{"metadata":{"trusted":false,"_uuid":"789ec329e05003df2e10d980a6ea650c5b620817"},"cell_type":"code","source":"from sklearn.tree import DecisionTreeClassifier\n\nprint('-----')\ntest_model(DecisionTreeClassifier())\n\nprint('-----')\ntest_model(DecisionTreeClassifier(max_depth = 10,max_features=5))\n\nprint('-----')\ntest_model(DecisionTreeClassifier(max_depth = 12,max_features=8))\n\nprint('-----')\ntest_model(DecisionTreeClassifier(max_depth = 8,max_features=10))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"074fc9f9fa062d9c5e963e4c447884a1d3b2d674"},"cell_type":"markdown","source":"### Random Forest"},{"metadata":{"trusted":false,"_uuid":"5bb6178f4b6984176498ca9af4f810e7c13e2e61"},"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\n\ntest_model(RandomForestClassifier())\n\nprint('-----')\ntest_model(RandomForestClassifier(max_depth = 10,max_features=5))\n\nprint('-----')\ntest_model(RandomForestClassifier(max_depth = 12,max_features=7))\n\nprint('-----')\ntest_model(RandomForestClassifier(max_depth = 8,max_features=10))\n\nprint('-----')\ntest_model(RandomForestClassifier(max_depth = 10,max_features=5,n_estimators=20))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"07229286fbbde3e5c1d01e705edc3067f7891169"},"cell_type":"markdown","source":"### Neural Nets"},{"metadata":{"trusted":false,"_uuid":"01453a531948f6cd5f7fcf13efa639902947fac3"},"cell_type":"code","source":"from sklearn.neural_network import MLPClassifier\n\ntest_model(MLPClassifier())\n\nprint('-----')\ntest_model(MLPClassifier(hidden_layer_sizes=(150,) max_iter=200))","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0e325137b0e7cbaa4bb2e16f1bc974f74a71ac86"},"cell_type":"code","source":"candidate_models","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"d489ec851e14df941172c2959aff6c88615b421f"},"cell_type":"code","source":"best_model = max(candidate_models, key=candidate_models.get)\ntype(best_model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0338d3409f5ec96dd218246941d58eb148e78cef"},"cell_type":"code","source":"candidate_models[best_model]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"d86ddc3e7442394e57d1fc6e91dfaf86fe4a6167"},"cell_type":"code","source":"pred_file(best_model,variables)","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","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.6.5"}},"nbformat":4,"nbformat_minor":1}