{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "31ae9e5f-9955-d971-c2c5-141a7c290be5"
      },
      "source": [
        "# Matthews correlation coefficient\n",
        "\n",
        "See https://en.wikipedia.org/wiki/Matthews_correlation_coefficient\n",
        "\n",
        "Author: [CPMP](https://www.kaggle.com/cpmpml)\n",
        "\n",
        "A 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",
        "\n",
        "This 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",
        "\n",
        "import numpy as np # linear algebra\n",
        "import 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",
        "\n",
        "from sklearn.metrics import matthews_corrcoef\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import 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\n",
        "def 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\n",
        "def 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)\n",
        "y_prob = y_prob0 + 0.4 * np.random.rand(1000000) - 0.02\n",
        "\n",
        "\n",
        "y_true = (y_prob0 > 0.6).astype(int)\n",
        "best_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",
        "\n",
        "best_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
}