{"cells":[{"metadata":{"_uuid":"c908cd93441353464f793cf862e99a8af2a6b534"},"cell_type":"markdown","source":"# Flavours of Physics: τ→ 3μ Kaggle Competition Kernel"},{"metadata":{"_uuid":"df52bf91d784041433258eaa66d6c053a6d27045"},"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 certain training variables from the dataset, before training several models and selecting the best one for our predictions submission."},{"metadata":{"_uuid":"ef208fa93d91472caecd6926b0e60ed6c707325d"},"cell_type":"markdown","source":"## SET UP ------------------------------------------------------------------------------------"},{"metadata":{"trusted":false,"_uuid":"0cf01d37531813affa74f8df40927f212f0ce978"},"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":"a8f36fe27a7c2f89daf57a615b8d3ceaedb222a1"},"cell_type":"markdown","source":"### Check agreement test"},{"metadata":{"trusted":false,"_uuid":"a46340ada5d711ce655b4d86f3252ee3557139d5"},"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":"02fe6f5e37b43adef154e5ed8f7f7c4067234e01"},"cell_type":"markdown","source":"### Check correlation test"},{"metadata":{"trusted":false,"_uuid":"fcada28bcf3dca558dfc72e59de4f7a515707a23"},"cell_type":"code","source":"def check_corr_test(model,var):\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":"9b47ad16574e898e85700e11c48c2ff4bd0ca7cc"},"cell_type":"markdown","source":"### Compute weighted AUC on the training data with min_ANNmuon > 0.4"},{"metadata":{"trusted":false,"_uuid":"64b994545c072aacd920a63d4cb854f27e6af99f"},"cell_type":"code","source":"def comp_auc(model,var):\n    train_eval = train[train['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":{"trusted":false,"_uuid":"86257833a4553b6131438c1cbac5e7463d3a6dc1"},"cell_type":"code","source":"## combine tests into one function\n\ndef eval(model,var):\n    check_ag_test(model,var)\n    check_corr_test(model,var)\n    comp_auc(model,var)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d7b0e9c7888912649f8f6f45e2dff1d1132481b6"},"cell_type":"markdown","source":"### Predict test and create file"},{"metadata":{"trusted":false,"_uuid":"311866d1de9cef330586e812e12f7cc881d18aee"},"cell_type":"code","source":"def pred_file(model,var):\n\n    test = pd.read_csv(folder + 'test.csv', index_col='id')\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":"d46da5bb2e0bc66bdbb58a03e3089fc35ddf9585"},"cell_type":"markdown","source":"## END OF SETUP--------------------------------------------------------------------------------------"},{"metadata":{"_uuid":"b8031b80ebadce49ace29d102fa8d0902164e4e8"},"cell_type":"markdown","source":" "},{"metadata":{"_uuid":"564de0d8d18eab946102d7741f79bc7581d68fa1"},"cell_type":"markdown","source":"# Imports"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"5b146e57a9fb41b60a56201b99f810317239831e"},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport sklearn\nimport seaborn as sns\n\n#import evaluation\n#import evalfunctions\n\nversion = \"1.0\"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e9c492eb5d16b1f32cb7a82ae272da8d664d32d7"},"cell_type":"markdown","source":"# Read training data"},{"metadata":{"trusted":false,"_uuid":"1accca5255af10e269c7a6b3fef6e4cac8d3091a"},"cell_type":"code","source":"folder = '../input/'\ntrain = pd.read_csv(folder + 'training.csv', index_col='id')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0da061fb0233c52b711eeaadc67f41717f462221"},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"af1046385d07af189b6a5fb9dc674784259e42e7"},"cell_type":"markdown","source":"# Define training features"},{"metadata":{"_uuid":"c1bb2611553201785328a669aa46ac905a59561b"},"cell_type":"markdown","source":"Lets plot a seaborn heatmap to see the correlation of each variable to the signal:"},{"metadata":{"trusted":false,"_uuid":"60214f10d5be04b1cd30a8ee483190e4ff015255"},"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":"bf99434e09a61541128aaef8e6ec1bee45ecba49"},"cell_type":"code","source":"# this function returns a list of variables with a magnitude correlation with signal greater than n\n# excluding some variables not to be included in the training\n\ndef significantFactors(n):\n    x =[]\n    corr = train.corr()[\"signal\"]\n    for i in range(len(corr)):\n        if abs(corr[i])>n:\n            if(corr.index[i] not in [\"production\", \"min_ANNmuon\",\"signal\",\"mass\"]):\n                x.append(corr.index[i])\n    return x","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"4ff441ece75d89e355f543415b0d702489a2708c"},"cell_type":"code","source":"# these are some variables that we would definitely like 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":"bf1a129c18d32897180e4b31ff2c6a27ba6d9e88"},"cell_type":"markdown","source":"## Feature Scaling"},{"metadata":{"trusted":false,"_uuid":"fc41efbcf87b80351332975309296baefa203092"},"cell_type":"code","source":"# ↓↓↓ I'm not sure this improves efficiency\n\n#from sklearn.preprocessing import MinMaxScaler\n#scaler = MinMaxScaler()\n#scaled_data=train\n#scaled_data[variables] = scaler.fit_transform(scaled_data[variables])\n#train=scaled_data","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"633618f1a6dec6b0eeda83ca0f26cff87d558ed1"},"cell_type":"markdown","source":"# Let's train and compare some models!"},{"metadata":{"trusted":false,"_uuid":"04a4c66a57b8aae6d6a083613bb8b83c5ca29f81"},"cell_type":"code","source":"candidate_models = {}   # we'll store candidate models here\n\ndef test_model(model):\n    #if the model passes the tests...\n    if(check_corr_test(model,variables) and check_ag_test(model,variables)):\n        #...add it to the candidates\n        candidate_models[svc] = comp_auc(model,variables)\n        print('passed')\n    else:\n        print('failed')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cc5504f63675338d8e617472e56ce8026c23f52d"},"cell_type":"markdown","source":"### Gradient Boosting Classifier"},{"metadata":{"trusted":false,"_uuid":"996de949fe25a815296a29c07e6844b79e0632be"},"cell_type":"code","source":"from sklearn.ensemble import GradientBoostingClassifier\n\ngbc = GradientBoostingClassifier(n_estimators=40, learning_rate=0.01, subsample=0.7,\n                                      min_samples_leaf=10, max_depth=7, random_state=11)\ngbc.fit(train[variables], train['signal'])\n\n#eval(gbc,variables)\n\nif(check_corr_test(gbc,variables) and check_ag_test(gbc,variables)):\n    candidate_models[gbc] = comp_auc(gbc,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a9dc5c6b92a45d8eee06d3b48b5837809119463f"},"cell_type":"markdown","source":"### Logistic Regression"},{"metadata":{"trusted":false,"_uuid":"73b0154cb4ab8c1a04387c76ef7527686a7d1a44"},"cell_type":"code","source":"from sklearn.linear_model import LogisticRegression\n\nlr = LogisticRegression()\nlr.fit(train[variables], train['signal'])\n\n#eval(lr,variables)\n\nif(check_corr_test(lr,variables) and check_ag_test(lr,variables)):\n    candidate_models[lr] = comp_auc(lr,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"85bdab7e86fcf3e7cd242a11d79acb30ba6db08d"},"cell_type":"markdown","source":"### Gaussian Naive Bayes"},{"metadata":{"trusted":false,"_uuid":"430647bf69331ea03702c8bde083a672f516affe"},"cell_type":"code","source":"from sklearn.naive_bayes import GaussianNB\n\ngnb = GaussianNB()\ngnb.fit(train[variables], train['signal'])\n\n#eval(gnb,variables)\n\nif(check_corr_test(gnb,variables) and check_ag_test(gnb,variables)):\n    candidate_models[gnb] = comp_auc(gnb,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"476a587fffae79fa055a63d5bd1febf0717c42c7"},"cell_type":"markdown","source":"### K Neighbors Classifier"},{"metadata":{"trusted":false,"_uuid":"c8197818a764bb1b8cbf09763908a6a4a03e208f"},"cell_type":"code","source":"from sklearn.neighbors import KNeighborsClassifier\n\nknc = KNeighborsClassifier()\n\nknc.fit(train[variables], train['signal'])\n\n#eval(gbc,variables)\n\nif(check_corr_test(knc,variables) and check_ag_test(knc,variables)):\n    candidate_models[knc] = comp_auc(knc,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d4dc1fb89fa8a656828e050c5c797e3418b7d883"},"cell_type":"markdown","source":"### Decision Tree"},{"metadata":{"trusted":false,"_uuid":"af5096c6db1c4657497be9314ce503daf0ed1b18"},"cell_type":"code","source":"from sklearn.tree import DecisionTreeClassifier\n\ndtc = DecisionTreeClassifier(max_depth = 10,max_features=5 )\ndtc.fit(train[variables], train['signal'])\n\n#eval(dtc,variables)\n\nif(check_corr_test(dtc,variables) and check_ag_test(dtc,variables)):\n    candidate_models[dtc] = comp_auc(dtc,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"6d1297845242400622cdd75988bb57b12c036ae6"},"cell_type":"markdown","source":"### Random Forest"},{"metadata":{"trusted":false,"_uuid":"7319b1fc267f507152ef2013dbe6ccf74f4e0538"},"cell_type":"code","source":"from sklearn.ensemble import RandomForestClassifier\n\nrfc = RandomForestClassifier(max_depth = 15,max_features=5) # <<<< this was just from trial-and-error tweaking\n\nrfc.fit(train[variables], train['signal'])\n\n#eval(rfc,variables)\n\nif(check_corr_test(rfc,variables) and check_ag_test(rfc,variables)):\n    candidate_models[rfc] = comp_auc(rfc,variables)\n    print('passed')\nelse:\n    print(\"failed\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"317e5e615e5e33a5515013b1f78bd3e6d1ef1e8e"},"cell_type":"code","source":"candidate_models","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"a38f00cdc84c3107ff08da6b8e56a1a776e1c310"},"cell_type":"code","source":"best_model = max(candidate_models, key=candidate_models.get)\ntype(best_model)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"772f33d21f7d4c8763c9d5eec271e9b4390f782a"},"cell_type":"code","source":"candidate_models[best_model]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"b68fc39edd98f7f25ea3635e5cd13cdabf0b1356"},"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":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}