{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"31ae9e5f-9955-d971-c2c5-141a7c290be5"},"source":"# Matthews correlation coefficient\n\nSee https://en.wikipedia.org/wiki/Matthews_correlation_coefficient\n\nAuthor: [CPMP](https://www.kaggle.com/cpmpml)\n\nA fast implementation of [Anokas mcc optimization code](https://www.kaggle.com/c/bosch-production-line-performance/forums/t/22917/optimising-probabilities-binary-prediction-script).\n\nThis code takes as input probabilities, and selects the threshold that yields the best MCC score.  It is efficient enough to be used as a custom evaluation function in xgboost"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0a633cf2-8deb-f33a-a30e-a435661c7072"},"outputs":[],"source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"3d3eb384-ce43-a4c3-c4ca-a37e8e4c1346"},"outputs":[],"source":"%matplotlib inline\n\nfrom sklearn.metrics import matthews_corrcoef\n\nimport matplotlib.pyplot as plt\nimport numpy as np"},{"cell_type":"markdown","metadata":{"_cell_guid":"10575e78-0a7f-d203-b973-5c18a7950f08"},"source":"We compile the code with Numba.  If you don't have it installed, then just comment out the import and the @jit lines below."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d62f429d-342e-bc25-d97a-6a1ff2ef8bdc"},"outputs":[],"source":"from numba import jit\n\n@jit\ndef mcc(tp, tn, fp, fn):\n    sup = tp * tn - fp * fn\n    inf = (tp + fp) * (tp + fn) * (tn + fp) * (tn + fn)\n    if inf==0:\n        return 0\n    else:\n        return sup / np.sqrt(inf)\n\n@jit\ndef eval_mcc(y_true, y_prob, show=False):\n    idx = np.argsort(y_prob)\n    y_true_sort = y_true[idx]\n    n = y_true.shape[0]\n    nump = 1.0 * np.sum(y_true) # number of positive\n    numn = n - nump # number of negative\n    tp = nump\n    tn = 0.0\n    fp = numn\n    fn = 0.0\n    best_mcc = 0.0\n    best_id = -1\n    prev_proba = -1\n    best_proba = -1\n    mccs = np.zeros(n)\n    for i in range(n):\n        # all items with idx < i are predicted negative while others are predicted positive\n        # only evaluate mcc when probability changes\n        proba = y_prob[idx[i]]\n        if proba != prev_proba:\n            prev_proba = proba\n            new_mcc = mcc(tp, tn, fp, fn)\n            if new_mcc >= best_mcc:\n                best_mcc = new_mcc\n                best_id = i\n                best_proba = proba\n        mccs[i] = new_mcc\n        if y_true_sort[i] == 1:\n            tp -= 1.0\n            fn += 1.0\n        else:\n            fp -= 1.0\n            tn += 1.0\n    if show:\n        y_pred = (y_prob >= best_proba).astype(int)\n        score = matthews_corrcoef(y_true, y_pred)\n        print(score, best_mcc)\n        plt.plot(mccs)\n        return best_proba, best_mcc, y_pred\n    else:\n        return best_mcc"},{"cell_type":"markdown","metadata":{"_cell_guid":"e39bdb23-7b6d-ae71-5e80-0e74e6ee4e25"},"source":"Let's see what this gives on an example.  We also compute the MCC using scikit-learn function to make sure we compute the same."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d6164f6f-5b4b-50f5-40c6-8d68a1cfb07c"},"outputs":[],"source":"y_prob0 = np.random.rand(1000000)\ny_prob = y_prob0 + 0.4 * np.random.rand(1000000) - 0.02\n\n\ny_true = (y_prob0 > 0.6).astype(int)\nbest_proba, best_mcc, y_pred = eval_mcc(y_true, y_prob, True)"},{"cell_type":"markdown","metadata":{"_cell_guid":"82d2f2c0-defb-6935-72e7-9a766aabf8a2"},"source":"It is fast.  For a one million items input:"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"5c8300d0-9d30-e851-82d5-c3d39b64d04b"},"outputs":[],"source":"%timeit eval_mcc(y_true, y_prob)"},{"cell_type":"markdown","metadata":{"_cell_guid":"6d6d8048-9282-6be2-6fca-226b7de9ef78"},"source":"Probabilities can be identical for several values as pointed out by @Commander.  Let's test it."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"3062db82-a103-392d-0320-5ef349e9f82c"},"outputs":[],"source":"def roundn(yprob, scale):\n    return np.around(y_prob * scale) / scale\n\nbest_proba, best_mcc, y_pred = eval_mcc(y_true, roundn(y_prob, 100), True)\n            "},{"cell_type":"markdown","metadata":{"_cell_guid":"fa0e4867-e79c-1fee-3bd6-3368167b0ec7"},"source":"We see that the best mcc is lower than when probabilities were all different."},{"cell_type":"markdown","metadata":{"_cell_guid":"2f8eded0-bb73-7124-a68e-065bed3edc0f"},"source":"For use with xgboost, we wrap it to get the right input and output. "},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"ba7f396c-145f-e5fa-20df-039aca9bae1f"},"outputs":[],"source":"def mcc_eval(y_prob, dtrain):\n    y_true = dtrain.get_label()\n    best_mcc = eval_mcc(y_true, y_prob)\n    return 'MCC', best_mcc"},{"cell_type":"markdown","metadata":{"_cell_guid":"e713d9b8-acd4-b1f1-c1ae-eff0eb18bc36"},"source":"We can then use it with xgboost, for instance passing it as feval parameter as follow:"},{"cell_type":"markdown","metadata":{"_cell_guid":"5d15ef64-561a-01ce-72f4-554d0d5af43e"},"source":"    bst = xgb.train(params, Xtrain, \n                        num_boost_round=num_round, \n                        evals=watchlist,\n                        early_stopping_rounds=early_stopping_rounds, \n                        evals_result=evals_result, \n                        verbose_eval=verbose_eval,\n                        feval=mcc_eval, \n                        maximize=True,)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"bb466f27-cf2c-a760-cc4a-eda7899b1f8a"},"outputs":[],"source":""}],"metadata":{"_change_revision":0,"_is_fork":false,"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.5.2"}},"nbformat":4,"nbformat_minor":0}