Bayes_Error_Rate_Estimation.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "mount_file_id": "17cckAyOPGvw2iAdsINZfMzHMQt_9qNQ0",
      "authorship_tag": "ABX9TyPEvXpPzsDqTCwgvUMTzGPV",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/rflameiro/projects/blob/main/Bayes_Error_Rate_Estimation.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introduction\n",
        "\n",
        "## Bayes error rate\n",
        "\n",
        "Analogous to the irreducible error, the Bayes error rate (BER) is the lowest possible error rate (which is 1 - accuracy) any classifier can achieve on a fixed probability distribution. If the BER is non-zero, as it usually is, then the two classes overlap, and the best possible model will still make a few wrong predictions.\n",
        "\n",
        "Note that the bit about the \"fixed probability distribution\" is particularly relevant for cheminformatics, because as novel molecules are developed, the \"distribution\" that generated the samples changes, and the BER will likely change as well (it will probably increase).\n",
        "\n",
        "## Estimating the Bayes error rate\n",
        "\n",
        "Calculating an exact value for the BER is not possible for real problems. However, estimation methods applied to artificial datasets with known BERs have been shown to present good approximations to the real value.\n",
        "\n"
      ],
      "metadata": {
        "id": "wWvUbsDpV0tJ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Some methods with code available\n",
        "\n",
        "## Learning to Benchmark\n",
        "\n",
        "[Learning to Benchmark: Determining Best Achievable Misclassification Error from Training Data](https://arxiv.org/abs/1909.07192)\n",
        "\n",
        "This is an ensemble method that computes the Bayes error rate directly from training data based on an ensemble of ε-ball estimators and Chebyshev approximation. Yeah, this paper is more \"theory heavy\" than what I'm used to reading and I do not claim to fully understand its limitations. It seems to require >1000 samples to achieve a stable prediction of the error rate. The paper discusses classification problems with simulated datasets of up to 100 dimensions and states that\n",
        "\n",
        "> In larger dimensions, solving the optimization problem can be computationally difficult\n",
        "\n",
        "However, they estimate the BER for the MNIST dataset, with 70k samples described by 784 features. Therefore, for large chemical datasets descbribed as 1024-bit fingerprints, it might be a good idea to project the datapoints using PCA before using this method, but keep in mind that this might increase the estimated error rate.\n",
        "\n",
        "The full code implementation is open-source and available on [GitHub](https://github.com/mrtnoshad/Bayes_Error_Estimator/blob/master/Example.ipynb).\n",
        "\n",
        "## FeeBee\n",
        "\n",
        "In [Evaluating Bayes Error Estimators on Real-World Datasets with FeeBee](https://arxiv.org/pdf/2108.13034.pdf), the authors propose FeeBee, a framework for analyzing and comparing BER estimators on real-world datasets with unknown probability distributions.\n",
        "\n",
        "The idea of this work is to evaluate BER estimators by injecting \"label noise\" to existing datasets and estimating the corresponding BER on different noise levels to measure to what degree a BER estimator under/overestimates.\n",
        "\n",
        "Some highlights:\n",
        "\n",
        "- The BER estimators evauated calculate a lower and an upper estimate for the BER of the datasets (not a single value). LD(m) and UD(m) are the scores for the estimates of the lower and upper limits for each estimator.\n",
        "\n",
        "- ne claim of this work is that any method that provides a lower-bound estimate for BER that is lower than the known state-of-the-art for a dataset is \"clearly wrong\". This might not be the case, since \"overfitting the test set\" is a known issue, especially for models trained on widely used benchmarks.\n",
        "\n",
        "- The authors also show examples in which they apply feature transformation techniques to some of the datasets, although admitting that this might increase the BER.\n",
        "\n",
        "- Using the complete FeeBee framework to compare the performance of the estimators was a bit challenging to me, especially because the authors only considered image and text datasets, the code is formatted to be run using `slurm` and the notebook to analyze the results seems a bit convoluted. It also seems to require the value of the state-of-the-art for the dataset. However, the authors do comment that:\n",
        "\n",
        "> For LD(m) we observe that 1NN, kNN, kNN-LOO and GHP are consistently outperforming all the other methods, with either kNN or kNN-LOO being the best choice on each dataset.\n",
        "\n",
        "> For UD(m), we see that there is no method that consistently outperforms the others. Out of the well-performing methods on LD(m), 1NN and GHP are tangibly inferior to kNN and kNN-LOO\n",
        "\n",
        "> kNN and kNN-LOO seem to perform reasonably for determining both the upper and lower limits of BER.\n",
        "\n",
        "> BER estimators are often performing better on a certain type of feature transformations, e.g., 1NN and kNN perform better on pre-trained embeddings, whereas kNN-Extrapolate performs best under transformations that reduce the dimension, such as a low-dimensional PCA.\n",
        "\n",
        "> GHP and 1NN are consistently outperforming the other lower bound estimators whilst enabling easy hyperparameter selection on a single point estimate. We also observe that currently there does not exist an estimator that is able to simultaneously outperform other estimators on both lower and upper bounds.\n",
        "\n",
        "From this discussion, the kNN and GHP estimators seem to be reasonable choices to estimate the upper and lower bounds of the BER. Therefore, I have adapted the code for these methods to run on Colab (I also added the possibility to use Tanimoto as the distance metric in kNN).\n",
        "\n",
        "The full code implementation is open-source and available on [GitHub](https://github.com/DS3Lab/feebee/tree/main).\n",
        "\n",
        "## Bonus: Regression models\n",
        "\n",
        "In [Getting Real with Molecular Property Prediction](https://practicalcheminformatics.blogspot.com/2023/06/getting-real-with-molecular-property.html) and [How Good Could (Should) My Models Be?](http://practicalcheminformatics.blogspot.com/2019/07/how-good-could-should-my-models-be.html) Pat Walters uses a straightforward method to estimate the maximum achievable performance for a regression model trained on a dataset with experimental values and the known (or approximate) standard deviation of the measurements.\n",
        "\n",
        "The method consists in adding random error to every measurement, based on the expected standard deviation, then calculating the correlation coefficient on the \"Simulated x Experimental\" values. Larger experimental variance will result in more noisy plots and smaller correlation coefficients, and the resulting models will also reflect this uncertainty.\n",
        "\n",
        "Errors in measurements can range from around 0.1 log units to in-house measurements (precise assay, measurements made in the same lab on the same conditions) up to 0.5-1.0 log units with data that span multiple assays.\n",
        "\n",
        "As stated in [Getting Real with Molecular Property Prediction](https://practicalcheminformatics.blogspot.com/2023/06/getting-real-with-molecular-property.html):\n",
        "> \"experimental error can dramatically impact our ability to model the data.  If the experimental error is 0.6 logs, we can't expect a Pearson's r greater than 0.77.\"\n",
        "\n",
        "The code for implementing this analysis can be found in [this Google Colab Notebook](https://colab.research.google.com/github/PatWalters/maximum_correlation/blob/master/maximum_correlation.ipynb) and will not be discussed here."
      ],
      "metadata": {
        "id": "6PVNGkQtWd5J"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Code\n",
        "\n",
        "Let's see how to estimate the BER for the MoleculeNet BBBP dataset, available on https://moleculenet.org/datasets-1."
      ],
      "metadata": {
        "id": "xml8W_9XYyF-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Import dataset\n",
        "\n",
        "Dataset: BBBP from MoleculeNet after removing compounds that failed standardization. Descriptor = RDKit implementation of Morgan fingerprints (1024 bits, radius=3)."
      ],
      "metadata": {
        "id": "wFxIxr2E5xCg"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "from sklearn.model_selection import train_test_split"
      ],
      "metadata": {
        "id": "H8ag2QZq6HP8"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# import dataset\n",
        "url = 'https://raw.githubusercontent.com/rflameiro/Python_e_Quiminformatica/main/datasets/BBBP_morganFP_1024_radius3.csv'\n",
        "df_fp = pd.read_csv(url, sep=\";\", index_col=False)"
      ],
      "metadata": {
        "id": "76Ey5b5Q4glp"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "X = df_fp.iloc[:, :-1]\n",
        "y = df_fp.iloc[:, -1]\n",
        "\n",
        "print(X.shape, y.shape)"
      ],
      "metadata": {
        "id": "BeoN9XWJEB2w",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "f23beb93-8eff-4e3c-bc91-ae2761c7829e"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(1934, 1024) (1934,)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "X_train, X_test, y_train, y_test = train_test_split(\n",
        "    X, y, test_size=0.2, random_state=42)"
      ],
      "metadata": {
        "id": "s-gMtqY0EQjF"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "X_train = X_train.to_numpy()\n",
        "X_test = X_test.to_numpy()\n",
        "y_train = y_train.to_numpy()\n",
        "y_test = y_test.to_numpy()"
      ],
      "metadata": {
        "id": "_AaRvFL0po70"
      },
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Learning to Benchmark\n",
        "\n",
        "Adapted from [here](https://github.com/mrtnoshad/Bayes_Error_Estimator/blob/master/Example.ipynb)"
      ],
      "metadata": {
        "id": "g7Kkm3g_5QHQ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "!git clone https://github.com/mrtnoshad/Bayes_Error_Estimator/"
      ],
      "metadata": {
        "id": "v2soaayY5OKj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Important: I found two problems using this code on Colab:\n",
        "\n",
        "First, I had to delete the line `matplotlib.use('Agg')` on BER_estimator_KDtree.py to prevent `ImportError: Cannot load backend 'TkAgg' which requires the 'tk' interactive framework, as 'headless' is currently running`\n",
        "\n",
        "Also, on the same file, I had to to replace `asscalar`, a deprecated function, to `ndarray.item`.\n",
        "\n",
        "Then I deleted the file on the Bayes_Error_Estimator Colab folder and re-uploaded the fixed version. I already created a pull request and an issue on GitHub."
      ],
      "metadata": {
        "id": "UNo5qnBULXk6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import sys\n",
        "\n",
        "sys.path.append(\"/content/Bayes_Error_Estimator\")"
      ],
      "metadata": {
        "id": "_XKfLd8H5mUp"
      },
      "execution_count": 7,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "from BER_estimator_KDtree import ensemble_bg_estimator as BER"
      ],
      "metadata": {
        "id": "KOadLakvLR6K"
      },
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Estimate Bayes Error Rate\n",
        "est_BER = BER(X_train, y_train)\n",
        "\n",
        "print(est_BER)"
      ],
      "metadata": {
        "id": "ghimFx965J-m",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c443a743-83a2-4475-b887-988323fb55b5"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.02134041599271922\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "This indicates that for the BBBP dataset the Bayes error rate is around 2%, which seems a bit low for a cheminformatics dataset. Let's use another method and compare."
      ],
      "metadata": {
        "id": "790ndtwYRqrL"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## FeeBee\n",
        "\n",
        "These functions are used to estimate the lower and upper bounds of the BER for binary classifiers instead of just one value. They were adapted from https://github.com/DS3Lab/feebee/tree/main.\n",
        "\n",
        "`knn_eval_from_matrices_split()` was adapted from https://github.com/DS3Lab/feebee/blob/main/methods/knn.py#L145. It uses the k-NN method and the dataset needs to be split into train and test sets. I added the option to use Tanimoto distances, as this is widely used in cheminformatics.\n",
        "\n",
        "`ghp_eval_from_matrix()` was adapted from https://github.com/DS3Lab/feebee/blob/main/methods/ghp.py#L15. It uses the GHP method and only the train set is used."
      ],
      "metadata": {
        "id": "bFz5zb-TnBfx"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ncKTQT17y6kC",
        "outputId": "73546aa1-5cfa-4635-e3fc-1f88ff79e692"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Cloning into 'feebee'...\n",
            "remote: Enumerating objects: 311, done.\u001b[K\n",
            "remote: Total 311 (delta 0), reused 0 (delta 0), pack-reused 311\u001b[K\n",
            "Receiving objects: 100% (311/311), 812.50 KiB | 8.46 MiB/s, done.\n",
            "Resolving deltas: 100% (165/165), done.\n"
          ]
        }
      ],
      "source": [
        "!git clone https://github.com/DS3Lab/feebee"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import sys\n",
        "\n",
        "sys.path.append(\"/content/feebee\")"
      ],
      "metadata": {
        "id": "POuaUoelzC7m"
      },
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import math\n",
        "import numpy as np\n",
        "import os.path as path\n",
        "import pandas as pd\n",
        "import random\n",
        "from scipy.sparse.csgraph import minimum_spanning_tree\n",
        "from sklearn.model_selection import train_test_split\n",
        "import time\n",
        "\n",
        "from methods.utils import compute_distance_matrix_loo, knn_errorrate, knn_errorrate_loo, load_data, load_embedding_fn\n",
        "# use this instead of:\n",
        "# from methods.utils import *\n",
        "# because I reimplemented compute_distance_matrix\n",
        "\n",
        "import transformations.label_noise as label_noise"
      ],
      "metadata": {
        "id": "vOFKCTMKzj0m"
      },
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Function to calculate distances for kNN"
      ],
      "metadata": {
        "id": "C8Fu55DaJJsr"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# from sklearn.metrics import jaccard_score\n",
        "\n",
        "def compute_distance_matrix(x_train, x_test, measure=\"tanimoto\"):\n",
        "    \"\"\"Calculates the distance matrix between test and train.\n",
        "\n",
        "    Args:\n",
        "    x_train: Matrix (NxD) where each row represents a training sample\n",
        "    x_test: Matrix (MxD) where each row represents a test sample\n",
        "    measure: Distance measure (not necessarly metric) to use\n",
        "\n",
        "    Raises:\n",
        "    NotImplementedError: When the measure is not implemented\n",
        "\n",
        "    Returns:\n",
        "    Matrix (MxN) where element i,j is the distance between\n",
        "    x_test_i and x_train_j.\n",
        "\n",
        "    ------\n",
        "\n",
        "    Removed the option to use tensorflow\n",
        "    Added the option to use Tanimoto distances: measure=\"tanimoto\"\n",
        "\n",
        "    \"\"\"\n",
        "\n",
        "    # Type check\n",
        "    if x_train.dtype != np.float32:\n",
        "        x_train = np.float32(x_train)\n",
        "    if x_test.dtype != np.float32:\n",
        "        x_test = np.float32(x_test)\n",
        "\n",
        "    # First version, took 18 min to run on the BBBP dataset. I remade it using matrix operations (below)\n",
        "    # if measure == \"tanimoto\":\n",
        "        # M = x_test.shape[0]\n",
        "        # N = x_train.shape[0]\n",
        "\n",
        "        # # Initialize an empty MxN matrix for Jaccard distances\n",
        "        # x_xt = np.zeros((M, N))\n",
        "\n",
        "        # # Iterate through each pair of rows from x_test and x_train\n",
        "        # for i in range(M):\n",
        "        #     for j in range(N):\n",
        "        #         # Calculate the Jaccard distance between the i-th row of x_test and the j-th row of x_train\n",
        "        #         jaccard_distance = 1.0 - jaccard_score(x_test[i], x_train[j])\n",
        "\n",
        "        #         # Store the Jaccard distance in the matrix\n",
        "        #         x_xt[i, j] = jaccard_distance\n",
        "\n",
        "    if measure == \"tanimoto\":\n",
        "        # Assumes x_test and x_train are matrices in which each row is a binary vector.\n",
        "        intersection = np.dot(x_test, x_train.T)\n",
        "        union = np.sum(x_test, axis=1, keepdims=True) + np.sum(x_train, axis=1, keepdims=True).T - intersection\n",
        "        # Compute Jaccard distance matrix. np.where is used to account for division by zero errors\n",
        "        x_xt = np.where(union == 0, 0.0, 1.0 - (intersection / union))\n",
        "\n",
        "    elif measure == \"squared_l2\":\n",
        "        x_xt = np.matmul(x_test, np.transpose(x_train))\n",
        "\n",
        "        x_train_2 = np.sum(np.square(x_train), axis=1)\n",
        "        x_test_2 = np.sum(np.square(x_test), axis=1)\n",
        "\n",
        "        for i in range(np.shape(x_xt)[0]):\n",
        "            x_xt[i, :] = np.multiply(x_xt[i, :], -2)\n",
        "            x_xt[i, :] = np.add(x_xt[i, :], x_test_2[i])\n",
        "            x_xt[i, :] = np.add(x_xt[i, :], x_train_2)\n",
        "\n",
        "    elif measure == \"cosine\":\n",
        "        x_xt = np.matmul(x_test, np.transpose(x_train))\n",
        "\n",
        "        x_train_2 = np.linalg.norm(x_train, axis=1)\n",
        "        x_test_2 = np.linalg.norm(x_test, axis=1)\n",
        "\n",
        "        outer = np.outer(x_test_2, x_train_2)\n",
        "        x_xt = np.ones(np.shape(x_xt)) - np.divide(x_xt, outer)\n",
        "\n",
        "    else:\n",
        "        raise NotImplementedError(\"Method '{}' is not implemented\".format(measure))\n",
        "\n",
        "    return x_xt"
      ],
      "metadata": {
        "id": "4ptzv-0p0Hat"
      },
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### BER estimator: kNN"
      ],
      "metadata": {
        "id": "9SCRa2tkJFOL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# knn.py\n",
        "\n",
        "KEY_PATTERN = \"measure={0}, k={1}\"\n",
        "\n",
        "def knn_get_lowerbound(value, k, classes):\n",
        "\n",
        "    if value <= 1e-10:\n",
        "        return 0.0\n",
        "\n",
        "    if classes > 2 or k == 1:\n",
        "        return ((classes - 1.0)/float(classes)) * (1.0 - math.sqrt(max(0.0, 1 - ((float(classes) / (classes - 1.0)) * value))))\n",
        "\n",
        "    if k > 2:\n",
        "        return value / float(1 + (1.0/math.sqrt(k)))\n",
        "\n",
        "    return value / float(1 + math.sqrt(2.0/k))\n",
        "\n",
        "\n",
        "def knn_eval_from_matrices_split(train_features, test_features, train_labels,\n",
        "                                 test_labels, knn_measure='tanimoto',\n",
        "                                 knn_k=[1,3,5,7,9], knn_subtest=None, knn_subtrain=None):\n",
        "\n",
        "    total_classes = np.unique(np.concatenate((train_labels, test_labels))).size\n",
        "    vals = {}\n",
        "    ks = sorted(set(knn_k), reverse=True)\n",
        "\n",
        "    if knn_subtest is not None:\n",
        "        test_dividable = (test_features.shape[0] % knn_subtest) == 0\n",
        "        test_iterations = (test_features.shape[0] + knn_subtest - 1 ) // knn_subtest\n",
        "    else:\n",
        "        test_dividable = True\n",
        "        test_iterations = 1\n",
        "\n",
        "    for i in range(test_iterations):\n",
        "        if test_iterations == 1:\n",
        "            current_pos = 0\n",
        "            current_samples = test_features.shape[0]\n",
        "        else:\n",
        "            current_pos = i*knn_subtest\n",
        "            current_samples = min(test_features.shape[0],(i+1)*knn_subtest) - current_pos\n",
        "\n",
        "        start = time.time()\n",
        "        d = compute_distance_matrix(train_features, test_features[current_pos:(current_pos + current_samples),:], knn_measure)\n",
        "        end = time.time()\n",
        "\n",
        "        start = time.time()\n",
        "        err = knn_errorrate(d,\n",
        "                            train_labels,\n",
        "                            test_labels[current_pos:(current_pos + current_samples)],\n",
        "                            k=ks)\n",
        "        end = time.time()\n",
        "\n",
        "\n",
        "        if not test_dividable:\n",
        "            err = [e * (current_samples/test_features.shape[0]) for e in err]\n",
        "        else:\n",
        "            err = [e * (1.0/test_iterations) for e in err]\n",
        "\n",
        "        for idx, k in enumerate(ks):\n",
        "            if k not in vals:\n",
        "                vals[k] = err[idx]\n",
        "            else:\n",
        "                vals[k] += err[idx]\n",
        "\n",
        "    results = {}\n",
        "    for k, v in vals.items():\n",
        "        results[KEY_PATTERN.format(knn_measure, k)] = [v, knn_get_lowerbound(v, k, total_classes)]\n",
        "\n",
        "    return results\n"
      ],
      "metadata": {
        "id": "XaiXvlWB0bPg"
      },
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### BER estimator: GHP"
      ],
      "metadata": {
        "id": "xadiocDrJCps"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# ghp.py\n",
        "\n",
        "def ghp_eval_from_matrix(train_features, train_labels, ghp_approx=None):\n",
        "\n",
        "    start = time.time()\n",
        "    d = compute_distance_matrix_loo(train_features, \"cosine\")\n",
        "\n",
        "    if ghp_approx:\n",
        "        indices = np.argpartition(d, ghp_approx, axis=1)[:,:ghp_approx]\n",
        "        for row_i in range(d.shape[0]):\n",
        "            mask = np.ones(d.shape[1],dtype=bool)\n",
        "            mask[indices[row_i,:]] = False\n",
        "            d[row_i, mask] = 0.0\n",
        "\n",
        "    end = time.time()\n",
        "\n",
        "    start = time.time()\n",
        "    d = minimum_spanning_tree(d).tocoo()\n",
        "    end = time.time()\n",
        "\n",
        "    start = time.time()\n",
        "\n",
        "    classes, classes_counts = np.unique(train_labels, return_counts=True)\n",
        "\n",
        "    num_train_samples = train_labels.size\n",
        "    num_classes = len(classes)\n",
        "\n",
        "    mapping = {}\n",
        "    idx = 0\n",
        "    for c in classes:\n",
        "        mapping[c] = idx\n",
        "        idx += 1\n",
        "\n",
        "    deltas = []\n",
        "    for i in range(num_classes-1):\n",
        "        deltas.append([0.0]*(num_classes-i-1))\n",
        "\n",
        "    # Calculate number of dichotomous edges\n",
        "    for i in range(num_train_samples - 1):\n",
        "        label_1 = mapping[train_labels[d.row[i]]]\n",
        "        label_2 = mapping[train_labels[d.col[i]]]\n",
        "        if label_1 == label_2:\n",
        "            continue\n",
        "        if label_1 > label_2:\n",
        "            tmp = label_1\n",
        "            label_1 = label_2\n",
        "            label_2 = tmp\n",
        "        deltas[label_1][label_2 - label_1 - 1] += 1\n",
        "\n",
        "    # Divide the number of dichotomous edges by 2 * num_train_samples to get estimator of deltas\n",
        "    deltas = [[item / (2.0 * num_train_samples) for item in sublist] for sublist in deltas]\n",
        "\n",
        "    # Sum up all the deltas\n",
        "    delta_sum = sum([sum(sublist) for sublist in deltas])\n",
        "\n",
        "    end = time.time()\n",
        "    # logging.log(logging.DEBUG, \"Estimators computed in {} seconds\".format(end - start))\n",
        "\n",
        "    upper = 2.0 * delta_sum\n",
        "\n",
        "    lower = ((num_classes - 1.0) / float(num_classes)) * (1.0 - math.sqrt(max(0.0, 1.0 - ((2.0 * num_classes)/(num_classes - 1.0) * delta_sum))))\n",
        "\n",
        "    return {\"default\": [upper, lower]}"
      ],
      "metadata": {
        "id": "jWkK1ndH44Go"
      },
      "execution_count": 15,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Estimate lower and upper bounds of BER on MoleculeNet BBBP dataset"
      ],
      "metadata": {
        "id": "xbeL6uvoI-LD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# kNN cosine\n",
        "df = knn_eval_from_matrices_split(X_train, X_test, y_train, y_test, knn_measure='cosine')\n",
        "df"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YNgzjqM1IvMG",
        "outputId": "a26a1fc2-abce-4d49-fe48-49aa767a5e1e"
      },
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'measure=cosine, k=9': [0.12144702842377261, 0.09108527131782947],\n",
              " 'measure=cosine, k=7': [0.12144702842377261, 0.08813509404821888],\n",
              " 'measure=cosine, k=5': [0.10852713178294573, 0.07499040371124602],\n",
              " 'measure=cosine, k=3': [0.11627906976744186, 0.07371797630413504],\n",
              " 'measure=cosine, k=1': [0.12661498708010335, 0.0679207173909534]}"
            ]
          },
          "metadata": {},
          "execution_count": 16
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# kNN Tanimoto\n",
        "df = knn_eval_from_matrices_split(X_train, X_test, y_train, y_test, knn_measure='tanimoto')\n",
        "df"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "aUo-wu7lyw2G",
        "outputId": "9a87e466-2483-4290-fc40-60fa2d62b551"
      },
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'measure=tanimoto, k=9': [0.11627906976744186, 0.0872093023255814],\n",
              " 'measure=tanimoto, k=7': [0.12403100775193798, 0.09001030881520225],\n",
              " 'measure=tanimoto, k=5': [0.10852713178294573, 0.07499040371124602],\n",
              " 'measure=tanimoto, k=3': [0.11627906976744186, 0.07371797630413504],\n",
              " 'measure=tanimoto, k=1': [0.13178294573643412, 0.07092130426717413]}"
            ]
          },
          "metadata": {},
          "execution_count": 17
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# GHP - uses cosine distance\n",
        "df = ghp_eval_from_matrix(X_train, y_train)\n",
        "df"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "INXgdUVV4nZZ",
        "outputId": "f211b320-e733-427f-acfe-978a5338123b"
      },
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "WARNING:tensorflow:From /content/feebee/methods/utils.py:104: is_gpu_available (from tensorflow.python.framework.test_util) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use `tf.config.list_physical_devices('GPU')` instead.\n"
          ]
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "{'default': [0.16095669036845509, 0.08826992481023144]}"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "It seems that the Bayes error rate for the BBBP dataset is between 0.07 and 0.12, which translates to an accuracy range of 88-93%. I believe this is a more reasonable estimate than the one from \"Learning to Benchmark\".\n",
        "\n",
        "Let's see what we can get with a Random Forest model:"
      ],
      "metadata": {
        "id": "W6ve9PCD3oZN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Accuracy for a non-optimized Random Forest model"
      ],
      "metadata": {
        "id": "IRsmTNvreEYR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from sklearn.ensemble import RandomForestClassifier\n",
        "from sklearn.model_selection import cross_val_score\n",
        "from sklearn.metrics import accuracy_score\n",
        "\n",
        "rf_classifier = RandomForestClassifier(random_state=42)\n",
        "\n",
        "# cross-validation on the training set\n",
        "cv_scores = cross_val_score(rf_classifier, X_train, y_train, cv=5)\n",
        "\n",
        "# fit, predict test set\n",
        "rf_classifier.fit(X_train, y_train)\n",
        "y_pred = rf_classifier.predict(X_test)\n",
        "\n",
        "# Evaluate accuracy on the test set\n",
        "test_accuracy = accuracy_score(y_test, y_pred)\n",
        "\n",
        "# Print results\n",
        "print(\"Cross-validation scores:\", cv_scores)\n",
        "print(\"Mean CV Accuracy:\", cv_scores.mean())\n",
        "print(\"Test Set Accuracy:\", test_accuracy)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GGZATg38p8Hw",
        "outputId": "4300f59e-3864-46fb-b1f7-ccc7ce7c3b89"
      },
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Cross-validation scores: [0.88064516 0.87741935 0.86731392 0.89967638 0.87055016]\n",
            "Mean CV Accuracy: 0.8791209938406931\n",
            "Test Set Accuracy: 0.875968992248062\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We got an accuracy estimate close to the upper bound of our BER estimate, so there might still be room for improvement on our model."
      ],
      "metadata": {
        "id": "Qj_DJFGL4lad"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Other methods, no code available\n",
        "\n",
        "## Ensemble method\n",
        "\n",
        "In the conference paper [Estimating the Bayes error rate through classifier combining](https://ieeexplore.ieee.org/document/546912), later published as [Bayes Error Rate Estimation Using Classifier Ensembles](https://www.tandfonline.com/doi/abs/10.1080/10255810305042), an approach to estimating the Bayes error rate based on classifier combining techniques (linear combining theory) was proposed.\n",
        "\n",
        "The method, which was shown to outperform classical methods for estimating the BER, uses a majority vote classifier, which outputs the most likely class according to the votes of multiple individual classifiers. However, it relies on classifiers that can reasonably approximate *a posteriori* class probabilities, which requires using ensembles of feed-forward neural networks (universal approximators that can generate class probabilities).\n",
        "\n",
        "The intuition behind the method is that while each individual classifier can make errors, given by the sum of the Bayes error and a random added error:\n",
        "\n",
        "$$ E_{total}= E_{Bayes} + E_{Added}$$\n",
        "\n",
        "by combining their outputs, we can average the added error contribution of each classifier:\n",
        "\n",
        "$$ E_{total}^{avg} = E_{Bayes} + E_{Added}^{avg}$$\n",
        "\n",
        "Then, we can finally apply an equation that calculates an estimate of the BER as a function of the individual classifier error, the combined classifier error, the number of classifiers combined (N) and the correlation among them (ρ).\n",
        "\n",
        "$$ E_{Bayes} = {N E_{total}^{avg} - ((N-1){\\rho}+1)E_{total} \\over (N-1)(1-{\\rho})} $$\n",
        "\n",
        "Therefore, by determining the amount of improvement that can be obtained from a combiner, the BER can be isolated and evaluated. This approach can be used to estimate the BER accurately, even for small training sets.\n",
        "\n",
        "Comment: While the paper seems really interesting and following it was not too difficult, I was not able to fully grasp the methodology and convert it to code. Also, these neural network ensembles seem to require a lot of computational power.\n",
        "\n",
        "## Soft Labels\n",
        "\n",
        "This method was recently published as a conference paper at ICLR 2023 as [Is the Performance of My Deep Network Too Good to Be True? A Direct Approach to Estimating the Bayes Error in Binary Classification](https://openreview.net/pdf?id=FZdJQgy05rz), and is also discussed [here](https://towardsdatascience.com/what-is-bayes-error-4bfadcc9c0ad).\n",
        "\n",
        "Binary classification problems are usually trained with samples and hard labels, that is, samples belong to one of two classes (0/1).\n",
        "\n",
        "Soft labels are values in the interval [0,1] that can indicate our certainty in a prediction. For instance, if a classifier predicts the labels of two samples as 0.01 and 0.3, we might classify both as belonging to class 0, but claim that the model is more certain about the first prediction (closer to 0).\n",
        "\n",
        "The authors work with datasets that are presented with soft labels, and show how they can be used to better estimate the Bayes error. Even though they also claim that it is possible to use their method with hard labels, I could not decipher the details of the implementation. However, they cite two works that involve datasets with hard labels, which are those that I have discussed (Learning to Benckmark and FeeBee)."
      ],
      "metadata": {
        "id": "TosuDPomWZXQ"
      }
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "eslAHG3bUcEM"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}