Analysis_25JUN2019_modular.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Data-generation-modules\" data-toc-modified-id=\"Data-generation-modules-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Data generation modules</a></span></li><li><span><a href=\"#Decider-modules\" data-toc-modified-id=\"Decider-modules-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Decider modules</a></span></li><li><span><a href=\"#Evaluator-modules\" data-toc-modified-id=\"Evaluator-modules-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Evaluator modules</a></span><ul class=\"toc-item\"><li><span><a href=\"#Convenience-functions\" data-toc-modified-id=\"Convenience-functions-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Convenience functions</a></span></li><li><span><a href=\"#Contraction-algorithm\" data-toc-modified-id=\"Contraction-algorithm-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Contraction algorithm</a></span></li><li><span><a href=\"#Evaluators\" data-toc-modified-id=\"Evaluators-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Evaluators</a></span></li></ul></li><li><span><a href=\"#Performance-comparison\" data-toc-modified-id=\"Performance-comparison-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Performance comparison</a></span><ul class=\"toc-item\"><li><span><a href=\"#Without-unobservables-in-the-data\" data-toc-modified-id=\"Without-unobservables-in-the-data-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Without unobservables in the data</a></span></li><li><span><a href=\"#With-unobservables-in-the-data\" data-toc-modified-id=\"With-unobservables-in-the-data-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>With unobservables in the data</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Refer to the `notes.tex` file for explanations about the modular framework."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scs\n",
    "import scipy.integrate as si\n",
    "import seaborn as sns\n",
    "import numpy.random as npr\n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Settings\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams.update({'font.size': 16})\n",
    "plt.rcParams.update({'figure.figsize': (10, 6)})\n",
    "\n",
    "# Suppress deprecation warnings.\n",
    "\n",
    "import warnings\n",
    "\n",
    "\n",
    "def fxn():\n",
    "    warnings.warn(\"deprecated\", DeprecationWarning)\n",
    "\n",
    "\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    fxn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data generation modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    '''Return value of sigmoid function (inverse of logit) at x.'''\n",
    "\n",
    "    return 1 / (1 + np.exp(-1 * x))\n",
    "\n",
    "\n",
    "def coinFlipDGWithoutUnobservables(N_total=50000):\n",
    "\n",
    "    df = pd.DataFrame()\n",
    "\n",
    "    # Sample feature X from standard Gaussian distribution, N(0, 1).\n",
    "    df = df.assign(X=npr.normal(size=N_total))\n",
    "\n",
    "    # Calculate P(Y=0|X=x) = 1 / (1 + exp(-X)) = sigmoid(X)\n",
    "    df = df.assign(probabilities_Y=sigmoid(df.X))\n",
    "\n",
    "    # Draw Y ~ Bernoulli(1 - sigmoid(X))\n",
    "    # Note: P(Y=1|X=x) = 1 - P(Y=0|X=x) = 1 - sigmoid(X)\n",
    "    results = npr.binomial(n=1, p=1 - df.probabilities_Y, size=N_total)\n",
    "\n",
    "    df = df.assign(result_Y=results)\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "def thresholdDGWithUnobservables(N_total=50000):\n",
    "\n",
    "    df = pd.DataFrame()\n",
    "\n",
    "    # Sample the variables from standard Gaussian distributions.\n",
    "    df = df.assign(X=npr.normal(size=N_total))\n",
    "    df = df.assign(Z=npr.normal(size=N_total))\n",
    "    df = df.assign(W=npr.normal(size=N_total))\n",
    "\n",
    "    # Calculate P(Y=0|X, Z, W)\n",
    "    probabilities_Y = sigmoid(beta_X * df.X + beta_Z * df.Z + beta_W * df.W)\n",
    "\n",
    "    df = df.assign(probabilities_Y=probabilities_Y)\n",
    "\n",
    "    # Result is 0 if P(Y = 0| X = x; Z = z; W = w) >= 0.5 , 1 otherwise\n",
    "    df = df.assign(result_Y=np.where(df.probabilities_Y >= 0.5, 0, 1))\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "def coinFlipDGWithUnobservables(N_total=50000,\n",
    "                                beta_X=1.0,\n",
    "                                beta_Z=1.0,\n",
    "                                beta_W=0.2):\n",
    "\n",
    "    df = pd.DataFrame()\n",
    "\n",
    "    # Sample feature X, Z and W from standard Gaussian distribution, N(0, 1).\n",
    "    df = df.assign(X=npr.normal(size=N_total))\n",
    "    df = df.assign(Z=npr.normal(size=N_total))\n",
    "    df = df.assign(W=npr.normal(size=N_total))\n",
    "\n",
    "    # Calculate P(Y=0|X=x) = 1 / (1 + exp(-X)) = sigmoid(X)\n",
    "    probabilities_Y = sigmoid(beta_X * df.X + beta_Z * df.Z + beta_W * df.W)\n",
    "\n",
    "    df = df.assign(probabilities_Y=probabilities_Y)\n",
    "\n",
    "    # Draw Y from Bernoulli distribution\n",
    "    results = npr.binomial(n=1,\n",
    "                           p=1 - df.probabilities_Y,\n",
    "                           size=N_total)\n",
    "\n",
    "    df = df.assign(result_Y=results)\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decider modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def humanDeciderLakkaraju(df,\n",
    "                          result_Y,\n",
    "                          featureX_col,\n",
    "                          featureZ_col,\n",
    "                          nJudges_M=100,\n",
    "                          beta_X=1,\n",
    "                          beta_Z=1,\n",
    "                          hide_unobserved=True):\n",
    "\n",
    "    # Assert that every judge will have the same number of subjects.\n",
    "    assert df.shape[0] % nJudges_M == 0, \"Can't assign subjets evenly!\"\n",
    "\n",
    "    # Compute the number of subjects allocated for each judge.\n",
    "    nSubjects_N = int(df.shape[0] / nJudges_M)\n",
    "\n",
    "    # Assign judge IDs as running numbering from 0 to nJudges_M - 1\n",
    "    df = df.assign(judgeID_J=np.repeat(range(0, nJudges_M), nSubjects_N))\n",
    "\n",
    "    # Sample acceptance rates uniformly from a closed interval\n",
    "    # from 0.1 to 0.9 and round to tenth decimal place.\n",
    "    acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 10)\n",
    "\n",
    "    # Replicate the rates so they can be attached to the corresponding judge ID.\n",
    "    df = df.assign(acceptanceRate_R=np.repeat(acceptance_rates, nSubjects_N))\n",
    "\n",
    "    probabilities_T = sigmoid(beta_X * df[featureX_col] + beta_Z * df[featureZ_col])\n",
    "    probabilities_T += np.sqrt(0.1) * npr.normal(size=nJudges_M * nSubjects_N)\n",
    "\n",
    "    df = df.assign(probabilities_T=probabilities_T)\n",
    "\n",
    "    # Sort by judges then probabilities in decreasing order\n",
    "    # Most dangerous for each judge are at the top.\n",
    "    df.sort_values(by=[\"judgeID_J\", \"probabilities_T\"],\n",
    "                   ascending=False,\n",
    "                   inplace=True)\n",
    "\n",
    "    # Iterate over the data. Subject will be given a negative decision\n",
    "    # if they are in the top (1-r)*100% of the individuals the judge will judge.\n",
    "    # I.e. if their within-judge-index is under 1 - acceptance threshold times\n",
    "    # the number of subjects assigned to each judge they will receive a\n",
    "    # negative decision.\n",
    "    df.reset_index(drop=True, inplace=True)\n",
    "\n",
    "    df['decision_T'] = np.where((df.index.values % nSubjects_N) <\n",
    "                                ((1 - df['acceptanceRate_R']) * nSubjects_N),\n",
    "                                0, 1)\n",
    "    \n",
    "    if hide_unobserved:\n",
    "        df.loc[df.decision_T == 0, result_Y] = np.nan\n",
    "\n",
    "    return df\n",
    "\n",
    "\n",
    "def coinFlipDecider(df,\n",
    "                    featureX_col,\n",
    "                    featureZ_col,\n",
    "                    nJudges_M=100,\n",
    "                    beta_X=1,\n",
    "                    beta_Z=1,\n",
    "                    hide_unobserved=True):\n",
    "\n",
    "    # Assert that every judge will have the same number of subjects.\n",
    "    assert df.shape[0] % nJudges_M == 0, \"Can't assign subjets evenly!\"\n",
    "\n",
    "    # Compute the number of subjects allocated for each judge.\n",
    "    nSubjects_N = int(df.shape[0] / nJudges_M)\n",
    "\n",
    "    # Assign judge IDs as running numbering from 0 to nJudges_M - 1\n",
    "    df = df.assign(judgeID_J=np.repeat(range(0, nJudges_M), nSubjects_N))\n",
    "\n",
    "    # Sample acceptance rates uniformly from a closed interval\n",
    "    # from 0.1 to 0.9 and round to tenth decimal place.\n",
    "    #acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 10)\n",
    "    \n",
    "    # No real leniency here???\n",
    "    acceptance_rates = np.ones(nJudges_M)*0.5\n",
    "    \n",
    "    # Replicate the rates so they can be attached to the corresponding judge ID.\n",
    "    df = df.assign(acceptanceRate_R=np.repeat(acceptance_rates, nSubjects_N))\n",
    "\n",
    "    probabilities_T = sigmoid(beta_X * df[featureX_col] + beta_Z * df[featureZ_col])\n",
    "    #probabilities_T += np.sqrt(0.1) * npr.normal(size=nJudges_M * nSubjects_N)\n",
    "\n",
    "    df = df.assign(probabilities_T=probabilities_T)\n",
    "\n",
    "    # Draw T from Bernoulli distribution\n",
    "    decisions = npr.binomial(n=1, p=1 - df.probabilities_T, size=df.shape[0])\n",
    "\n",
    "    df = df.assign(decision_T=decisions)\n",
    "    \n",
    "    if hide_unobserved:\n",
    "        df.loc[df.decision_T == 0, 'result_Y'] = np.nan\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluator modules\n",
    "\n",
    "### Convenience functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fitPredictiveModel(x_train, y_train, x_test, class_value, model_type=None):\n",
    "    '''\n",
    "    Fit a predictive model (default logistic regression) with given training \n",
    "    instances and return probabilities for test instances to obtain a given \n",
    "    class label.\n",
    "    \n",
    "    Arguments:\n",
    "    ----------\n",
    "    \n",
    "    x_train -- x values of training instances\n",
    "    y_train -- y values of training instances\n",
    "    x_test -- x values of test instances\n",
    "    class_value -- class label for which the probabilities are counted for.\n",
    "    model_type -- type of model to be fitted.\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    (1) Trained predictive model\n",
    "    (2) Probabilities for given test inputs for given class.\n",
    "    '''\n",
    "\n",
    "    if model_type is None or model_type in [\"logistic_regression\", \"lr\"]:\n",
    "        # Instantiate the model (using the default parameters)\n",
    "        logreg = LogisticRegression(solver='lbfgs')\n",
    "\n",
    "        # Check shape and fit the model.\n",
    "        if x_train.ndim == 1:\n",
    "            logreg = logreg.fit(x_train.values.reshape(-1, 1), y_train)\n",
    "        else:\n",
    "            logreg = logreg.fit(x_train, y_train)\n",
    "\n",
    "        label_probs_logreg = getProbabilityForClass(x_test, logreg,\n",
    "                                                    class_value)\n",
    "\n",
    "        return logreg, label_probs_logreg\n",
    "\n",
    "    elif model_type in [\"random_forest\", \"rf\"]:\n",
    "        # Instantiate the model\n",
    "        forest = RandomForestClassifier(n_estimators=100, max_depth=3)\n",
    "\n",
    "        # Check shape and fit the model.\n",
    "        if x_train.ndim == 1:\n",
    "            forest = forest.fit(x_train.values.reshape(-1, 1), y_train)\n",
    "        else:\n",
    "            forest = forest.fit(x_train, y_train)\n",
    "\n",
    "        label_probs_forest = getProbabilityForClass(x_test, forest,\n",
    "                                                    class_value)\n",
    "\n",
    "        return forest, label_probs_forest\n",
    "\n",
    "    elif model_type == \"fully_random\":\n",
    "\n",
    "        label_probs = np.ones_like(x_test) / 2\n",
    "\n",
    "        model_object = lambda x: 0.5\n",
    "\n",
    "        return model_object, label_probs\n",
    "    else:\n",
    "        raise ValueError(\"Invalid model_type!\", model_type)\n",
    "\n",
    "\n",
    "def getProbabilityForClass(x, model, class_value):\n",
    "    '''\n",
    "    Function (wrapper) for obtaining the probability of a class given x and a \n",
    "    predictive model.\n",
    "\n",
    "    Arguments:\n",
    "    -----------\n",
    "    x -- individual features, an array of shape (observations, features)\n",
    "    model -- a trained sklearn model. Predicts probabilities for given x. \n",
    "        Should accept input of shape (observations, features)\n",
    "    class_value -- the resulting class to predict (usually 0 or 1).\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    (1) The probabilities of given class label for each x.\n",
    "    '''\n",
    "    if x.ndim == 1:\n",
    "        # if x is vector, transform to column matrix.\n",
    "        f_values = model.predict_proba(np.array(x).reshape(-1, 1))\n",
    "    else:\n",
    "        f_values = model.predict_proba(x)\n",
    "\n",
    "    # Get correct column of predicted class, remove extra dimensions and return.\n",
    "    return f_values[:, model.classes_ == class_value].flatten()\n",
    "\n",
    "\n",
    "def cdf(x_0, model, class_value):\n",
    "    '''\n",
    "    Cumulative distribution function as described above. Integral is \n",
    "    approximated using Simpson's rule for efficiency.\n",
    "    \n",
    "    Arguments:\n",
    "    ----------\n",
    "    \n",
    "    x_0 -- private features of an instance for which the value of cdf is to be\n",
    "        calculated.\n",
    "    model -- a trained sklearn model. Predicts probabilities for given x. \n",
    "        Should accept input of shape (observations, features)\n",
    "    class_value -- the resulting class to predict (usually 0 or 1).\n",
    "\n",
    "    '''\n",
    "\n",
    "    def prediction(x):\n",
    "        return getProbabilityForClass(\n",
    "            np.array([x]).reshape(-1, 1), model, class_value)\n",
    "\n",
    "    prediction_x_0 = prediction(x_0)\n",
    "\n",
    "    x_values = np.linspace(-15, 15, 40000)\n",
    "\n",
    "    x_preds = prediction(x_values)\n",
    "\n",
    "    y_values = scs.norm.pdf(x_values)\n",
    "\n",
    "    results = np.zeros(x_0.shape[0])\n",
    "    print(\"en loop\")\n",
    "    for i in range(x_0.shape[0]):\n",
    "\n",
    "        y_copy = y_values.copy()\n",
    "\n",
    "        y_copy[x_preds > prediction_x_0[i]] = 0\n",
    "        \n",
    "        results[i] = si.simps(y_copy, x=x_values)\n",
    "    print(\"jlk loop\")\n",
    "    return results\n",
    "\n",
    "\n",
    "def bailIndicator(r, y_model, x_train, x_test):\n",
    "    '''\n",
    "    Indicator function for whether a judge will bail or jail a suspect.\n",
    "    Rationale explained above.\n",
    "\n",
    "    Algorithm:\n",
    "    ----------\n",
    "\n",
    "    (1) Calculate recidivism probabilities from training set with a trained \n",
    "        model and assign them to predictions_train.\n",
    "\n",
    "    (2) Calculate recidivism probabilities from test set with the trained \n",
    "        model and assign them to predictions_test.\n",
    "\n",
    "    (3) Construct a quantile function of the probabilities in\n",
    "        in predictions_train.\n",
    "\n",
    "    (4)\n",
    "    For pred in predictions_test:\n",
    "\n",
    "        if pred belongs to a percentile (computed from step (3)) lower than r\n",
    "            return True\n",
    "        else\n",
    "            return False\n",
    "\n",
    "    Arguments:\n",
    "    ----------\n",
    "\n",
    "    r -- float, acceptance rate, between 0 and 1\n",
    "    y_model -- a trained sklearn predictive model to predict the outcome\n",
    "    x_train -- private features of the training instances\n",
    "    x_test -- private features of the test instances\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    (1) Boolean list indicating a bail decision (bail = True) for each \n",
    "        instance in x_test.\n",
    "    '''\n",
    "\n",
    "    predictions_train = getProbabilityForClass(x_train, y_model, 0)\n",
    "\n",
    "    predictions_test = getProbabilityForClass(x_test, y_model, 0)\n",
    "\n",
    "    return [\n",
    "        scs.percentileofscore(predictions_train, pred, kind='weak') < r\n",
    "        for pred in predictions_test\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contraction algorithm\n",
    "\n",
    "Below is an implementation of Lakkaraju's team's algorithm presented in [their paper](https://helka.finna.fi/PrimoRecord/pci.acm3098066). Relevant parameters to be passed to the function are presented in the description."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contraction(df, judgeIDJ_col, decisionT_col, resultY_col, modelProbS_col,\n",
    "                accRateR_col, r):\n",
    "    '''\n",
    "    This is an implementation of the algorithm presented by Lakkaraju\n",
    "    et al. in their paper \"The Selective Labels Problem: Evaluating \n",
    "    Algorithmic Predictions in the Presence of Unobservables\" (2017).\n",
    "\n",
    "    Arguments:\n",
    "    -----------\n",
    "    df -- The (Pandas) data frame containing the data, judge decisions,\n",
    "        judge IDs, results and probability scores.\n",
    "    judgeIDJ_col -- String, the name of the column containing the judges' IDs\n",
    "        in df.\n",
    "    decisionT_col -- String, the name of the column containing the judges' decisions\n",
    "    resultY_col -- String, the name of the column containing the realization\n",
    "    modelProbS_col -- String, the name of the column containing the probability\n",
    "        scores from the black-box model B.\n",
    "    accRateR_col -- String, the name of the column containing the judges' \n",
    "        acceptance rates\n",
    "    r -- Float between 0 and 1, the given acceptance rate.\n",
    "\n",
    "    Returns:\n",
    "    --------\n",
    "    (1) The estimated failure rate at acceptance rate r.\n",
    "    '''\n",
    "    # Get ID of the most lenient judge.\n",
    "    most_lenient_ID_q = df[judgeIDJ_col].loc[df[accRateR_col].idxmax()]\n",
    "\n",
    "    # Subset. \"D_q is the set of all observations judged by q.\"\n",
    "    D_q = df[df[judgeIDJ_col] == most_lenient_ID_q].copy()\n",
    "\n",
    "    # All observations of R_q have observed outcome labels.\n",
    "    # \"R_q is the set of observations in D_q with observed outcome labels.\"\n",
    "    R_q = D_q[D_q[decisionT_col] == 1].copy()\n",
    "\n",
    "    # Sort observations in R_q in descending order of confidence scores S and\n",
    "    # assign to R_sort_q.\n",
    "    # \"Observations deemed as high risk by B are at the top of this list\"\n",
    "    R_sort_q = R_q.sort_values(by=modelProbS_col, ascending=False)\n",
    "\n",
    "    number_to_remove = int(\n",
    "        round((1.0 - r) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])))\n",
    "\n",
    "    # \"R_B is the list of observations assigned to t = 1 by B\"\n",
    "    R_B = R_sort_q[number_to_remove:R_sort_q.shape[0]]\n",
    "\n",
    "    return np.sum(R_B[resultY_col] == 0) / D_q.shape[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def contractionEvaluator(df, featureX_col, judgeIDJ_col, decisionT_col,\n",
    "                         resultY_col, accRateR_col, r):\n",
    "\n",
    "    train, test = train_test_split(df, test_size=0.5)\n",
    "\n",
    "    B_model, predictions = fitPredictiveModel(\n",
    "        train.loc[train[decisionT_col] == 1, featureX_col],\n",
    "        train.loc[train[decisionT_col] == 1, resultY_col], test[featureX_col],\n",
    "        0)\n",
    "\n",
    "    test = test.assign(B_prob_0_model=predictions)\n",
    "\n",
    "    # Invoke the original contraction.\n",
    "    FR = contraction(test,\n",
    "                     judgeIDJ_col=judgeIDJ_col,\n",
    "                     decisionT_col=decisionT_col,\n",
    "                     resultY_col=resultY_col,\n",
    "                     modelProbS_col=\"B_prob_0_model\",\n",
    "                     accRateR_col=accRateR_col,\n",
    "                     r=r)\n",
    "\n",
    "    return FR\n",
    "\n",
    "\n",
    "def trueEvaluationEvaluator(df, featureX_col, decisionT_col, resultY_col, r):\n",
    "\n",
    "    train, test = train_test_split(df, test_size=0.5)\n",
    "\n",
    "    B_model, predictions = fitPredictiveModel(train[featureX_col],\n",
    "                                              train[resultY_col],\n",
    "                                              test[featureX_col], 0)\n",
    "\n",
    "    test = test.assign(B_prob_0_model=predictions)\n",
    "\n",
    "    test.sort_values(by='B_prob_0_model', inplace=True, ascending=True)\n",
    "\n",
    "    to_release = int(round(test.shape[0] * r / 10))\n",
    "\n",
    "    return np.sum(test[resultY_col][0:to_release] == 0) / test.shape[0]\n",
    "\n",
    "\n",
    "def labeledOutcomesEvaluator(df, featureX_col, decisionT_col, resultY_col, r):\n",
    "\n",
    "    train, test = train_test_split(df, test_size=0.5)\n",
    "\n",
    "    B_model, predictions = fitPredictiveModel(\n",
    "        train.loc[train[decisionT_col] == 1, featureX_col],\n",
    "        train.loc[train[decisionT_col] == 1, resultY_col], test[featureX_col],\n",
    "        0)\n",
    "\n",
    "    test = test.assign(B_prob_0_model=predictions)\n",
    "\n",
    "    test_observed = test.loc[test[decisionT_col] == 1, :]\n",
    "\n",
    "    test_observed = test_observed.sort_values(by='B_prob_0_model',\n",
    "                                              inplace=False,\n",
    "                                              ascending=True)\n",
    "\n",
    "    to_release = int(round(test_observed.shape[0] * r / 10))\n",
    "\n",
    "    return np.sum(\n",
    "        test_observed[resultY_col][0:to_release] == 0) / test.shape[0]\n",
    "\n",
    "\n",
    "def humanEvaluationEvaluator(df, judgeIDJ_col, decisionT_col, resultY_col,\n",
    "                             accRateR_col, r):\n",
    "\n",
    "    # Get judges with correct leniency as list\n",
    "    is_correct_leniency = df[accRateR_col].round(1) == r / 10\n",
    "\n",
    "    correct_leniency_list = df.loc[is_correct_leniency, judgeIDJ_col]\n",
    "\n",
    "    # Released are the people they judged and released, T = 1\n",
    "    released = df[df[judgeIDJ_col].isin(correct_leniency_list)\n",
    "                  & (df.decision_T == 1)]\n",
    "\n",
    "    # Get their failure rate, aka ratio of reoffenders to number of people judged in total\n",
    "    return np.sum(released[resultY_col] == 0) / correct_leniency_list.shape[0]\n",
    "\n",
    "\n",
    "def causalEvaluator(df, featureX_col, decisionT_col, resultY_col, r):\n",
    "\n",
    "    train, test = train_test_split(df, test_size=0.5)\n",
    "\n",
    "    B_model, predictions = fitPredictiveModel(\n",
    "        train.loc[train[decisionT_col] == 1, featureX_col],\n",
    "        train.loc[train[decisionT_col] == 1, resultY_col], test[featureX_col],\n",
    "        0)\n",
    "\n",
    "    test = test.assign(B_prob_0_model=predictions)\n",
    "\n",
    "    released = cdf(test[featureX_col], B_model, 0) < r / 10\n",
    "\n",
    "    return np.mean(test.B_prob_0_model * released)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance comparison\n",
    "\n",
    "Below we try to replicate the results obtained by Lakkaraju and compare their model's performance to the one of ours."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Without unobservables in the data\n",
    "\n",
    "The underlying figure is attached to the preliminary paper. When conducting finalization, last analysis should be conducted with a preset random seed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "hidden": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# f_rates = np.zeros((8, 5))\n",
    "# f_sems = np.zeros((8, 5))\n",
    "\n",
    "# nIter = 15\n",
    "\n",
    "# #npr.seed(0)\n",
    "\n",
    "# for r in np.arange(1, 9):\n",
    "\n",
    "#     print(\"[\", r, \"]\", sep='', end=\" \")\n",
    "\n",
    "#     s_f_rate_true = np.zeros(nIter)\n",
    "#     s_f_rate_labeled = np.zeros(nIter)\n",
    "#     s_f_rate_human = np.zeros(nIter)\n",
    "#     s_f_rate_cont = np.zeros(nIter)\n",
    "#     s_f_rate_caus = np.zeros(nIter)\n",
    "\n",
    "#     for i in range(nIter):\n",
    "\n",
    "#         print(i, end=\" \")\n",
    "\n",
    "#         s_train_labeled, s_train, s_test_labeled, s_test, s_df = dataWithoutUnobservables(sigma=2)\n",
    "\n",
    "#         s_logreg, predictions = fitPredictiveModel(\n",
    "#             s_train_labeled.dropna().X,\n",
    "#             s_train_labeled.dropna().result_Y, s_test.X, 0)\n",
    "#         s_test = s_test.assign(B_prob_0_model=predictions)\n",
    "\n",
    "#         s_logreg, predictions_labeled = fitPredictiveModel(\n",
    "#             s_train_labeled.dropna().X,\n",
    "#             s_train_labeled.dropna().result_Y, s_test_labeled.X, 0)\n",
    "#         s_test_labeled = s_test_labeled.assign(\n",
    "#             B_prob_0_model=predictions_labeled)\n",
    "\n",
    "#         #### True evaluation\n",
    "#         # Sort by actual failure probabilities, subjects with the smallest risk are first.\n",
    "#         s_sorted = s_test.sort_values(by='B_prob_0_model',\n",
    "#                                       inplace=False,\n",
    "#                                       ascending=True)\n",
    "\n",
    "#         to_release = int(round(s_sorted.shape[0] * r / 10))\n",
    "\n",
    "#         # Calculate failure rate as the ratio of failures to successes among those\n",
    "#         # who were given a positive decision, i.e. those whose probability of negative\n",
    "#         # outcome was low enough.\n",
    "#         s_f_rate_true[i] = np.sum(\n",
    "#             s_sorted.result_Y[0:to_release] == 0) / s_sorted.shape[0]\n",
    "\n",
    "#         #### Labeled outcomes\n",
    "#         # Sort by estimated failure probabilities, subjects with the smallest risk are first.\n",
    "#         s_sorted = s_test_labeled.sort_values(by='B_prob_0_model',\n",
    "#                                               inplace=False,\n",
    "#                                               ascending=True)\n",
    "\n",
    "#         to_release = int(round(s_test_labeled.dropna().shape[0] * r / 10))\n",
    "\n",
    "#         # Calculate failure rate as the ratio of failures to successes among those\n",
    "#         # who were given a positive decision, i.e. those whose probability of negative\n",
    "#         # outcome was low enough.\n",
    "#         s_f_rate_labeled[i] = np.sum(\n",
    "#             s_sorted.result_Y[0:to_release] == 0) / s_sorted.shape[0]\n",
    "\n",
    "#         #### Human error rate\n",
    "#         # Get judges with correct leniency as list\n",
    "#         correct_leniency_list = s_test_labeled.judgeID_J[\n",
    "#             s_test_labeled['acceptanceRate_R'].round(1) == r / 10].values\n",
    "\n",
    "#         # Released are the people they judged and released, T = 1\n",
    "#         released = s_test_labeled[\n",
    "#             s_test_labeled.judgeID_J.isin(correct_leniency_list)\n",
    "#             & (s_test_labeled.decision_T == 1)]\n",
    "\n",
    "#         # Get their failure rate, aka ratio of reoffenders to number of people judged in total\n",
    "#         s_f_rate_human[i] = np.sum(\n",
    "#             released.result_Y == 0) / correct_leniency_list.shape[0]\n",
    "\n",
    "#         #### Contraction\n",
    "#         s_f_rate_cont[i] = contraction(s_test_labeled, 'judgeID_J',\n",
    "#                                        'decision_T', 'result_Y',\n",
    "#                                        'B_prob_0_model', 'acceptanceRate_R',\n",
    "#                                        r / 10)\n",
    "#         #### Causal model\n",
    "\n",
    "#         #released = bailIndicator(r * 10, s_logreg, s_train.X, s_test.X)\n",
    "#         released=0\n",
    "#         #released = cdf(s_test.X, s_logreg, 0) < r / 10\n",
    "\n",
    "#         s_f_rate_caus[i] = np.mean(s_test.B_prob_0_model * released)\n",
    "\n",
    "#         ########################\n",
    "#         #percentiles = estimatePercentiles(s_train_labeled.X, s_logreg)\n",
    "\n",
    "#         #def releaseProbability(x):\n",
    "#         #    return calcReleaseProbabilities(r * 10,\n",
    "#         #                                     s_train_labeled.X,\n",
    "#         #                                     x,\n",
    "#         #                                     s_logreg,\n",
    "#         #                                     percentileMatrix=percentiles)\n",
    "\n",
    "#         #def integrand(x):\n",
    "#         #    p_y0 = s_logreg.predict_proba(x.reshape(-1, 1))[:, 0]\n",
    "\n",
    "#         #    p_t1 = releaseProbability(x)\n",
    "\n",
    "#         #    p_x = scs.norm.pdf(x)\n",
    "\n",
    "#         #    return p_y0 * p_t1 * p_x\n",
    "\n",
    "#         #s_f_rate_caus[i] = si.quad(lambda x: integrand(np.ones((1, 1)) * x),\n",
    "#         #                           -10, 10)[0]\n",
    "\n",
    "#     f_rates[r - 1, 0] = np.mean(s_f_rate_true)\n",
    "#     f_rates[r - 1, 1] = np.mean(s_f_rate_labeled)\n",
    "#     f_rates[r - 1, 2] = np.mean(s_f_rate_human)\n",
    "#     f_rates[r - 1, 3] = np.mean(s_f_rate_cont)\n",
    "#     f_rates[r - 1, 4] = np.mean(s_f_rate_caus)\n",
    "\n",
    "#     f_sems[r - 1, 0] = scs.sem(s_f_rate_true)\n",
    "#     f_sems[r - 1, 1] = scs.sem(s_f_rate_labeled)\n",
    "#     f_sems[r - 1, 2] = scs.sem(s_f_rate_human)\n",
    "#     f_sems[r - 1, 3] = scs.sem(s_f_rate_cont)\n",
    "#     f_sems[r - 1, 4] = scs.sem(s_f_rate_caus)\n",
    "\n",
    "# x_ax = np.arange(0.1, 0.9, 0.1)\n",
    "\n",
    "# plt.errorbar(x_ax,\n",
    "#              f_rates[:, 0],\n",
    "#              label='True Evaluation',\n",
    "#              c='green',\n",
    "#              yerr=f_sems[:, 0])\n",
    "# plt.errorbar(x_ax,\n",
    "#              f_rates[:, 1],\n",
    "#              label='Labeled outcomes',\n",
    "#              c='magenta',\n",
    "#              yerr=f_sems[:, 1])\n",
    "# plt.errorbar(x_ax,\n",
    "#              f_rates[:, 2],\n",
    "#              label='Human evaluation',\n",
    "#              c='red',\n",
    "#              yerr=f_sems[:, 2])\n",
    "# plt.errorbar(x_ax,\n",
    "#              f_rates[:, 3],\n",
    "#              label='Contraction, log.',\n",
    "#              c='blue',\n",
    "#              yerr=f_sems[:, 3])\n",
    "# # plt.errorbar(x_ax,\n",
    "# #              f_rates[:, 4],\n",
    "# #              label='Causal model, ep',\n",
    "# #              c='black',\n",
    "# #              yerr=f_sems[:, 4])\n",
    "\n",
    "# plt.title('Failure rate vs. Acceptance rate without unobservables')\n",
    "# plt.xlabel('Acceptance rate')\n",
    "# plt.ylabel('Failure rate')\n",
    "# plt.legend()\n",
    "# plt.grid()\n",
    "# plt.show()\n",
    "\n",
    "# print(f_rates)\n",
    "# print(\"\\nMean absolute errors:\")\n",
    "# for i in range(1, f_rates.shape[1]):\n",
    "#     print(np.mean(np.abs(f_rates[:, 0] - f_rates[:, i])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With unobservables in the data\n",
    "\n",
    "Lakkaraju says that they used logistic regression. We train the predictive models using only *observed observations*, i.e. observations for which labels are available. We then predict the probability of negative outcome for all observations in the test data and attach it to our data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0 en loop\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:78: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jlk loop\n",
      "1 en loop\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:78: RuntimeWarning: invalid value encountered in long_scalars\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-97-03cd8a3c6103>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     65\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     66\u001b[0m         f_rate_caus[i] = causalEvaluator(df_labeled, 'X', 'decision_T',\n\u001b[0;32m---> 67\u001b[0;31m                                          'result_Y', r / 10)\n\u001b[0m\u001b[1;32m     68\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m     \u001b[0mfailure_rates\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mr\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf_rate_true\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-96-4ef6c6b281a2>\u001b[0m in \u001b[0;36mcausalEvaluator\u001b[0;34m(df, featureX_col, decisionT_col, resultY_col, r)\u001b[0m\n\u001b[1;32m     90\u001b[0m     \u001b[0mtest\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massign\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mB_prob_0_model\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpredictions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     91\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 92\u001b[0;31m     \u001b[0mreleased\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtest\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mfeatureX_col\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mB_model\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mr\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     93\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     94\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtest\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mB_prob_0_model\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mreleased\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-94-f0303c92af6c>\u001b[0m in \u001b[0;36mcdf\u001b[0;34m(x_0, model, class_value)\u001b[0m\n\u001b[1;32m    123\u001b[0m         \u001b[0my_copy\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mx_preds\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mprediction_x_0\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    124\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 125\u001b[0;31m         \u001b[0mresults\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msi\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msimps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_copy\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mx_values\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    126\u001b[0m     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"jlk loop\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    127\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mresults\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadrature.py\u001b[0m in \u001b[0;36msimps\u001b[0;34m(y, x, dx, axis, even)\u001b[0m\n\u001b[1;32m    477\u001b[0m                 \u001b[0mfirst_dx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mslice2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mslice1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    478\u001b[0m             \u001b[0mval\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mfirst_dx\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mslice2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mslice1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 479\u001b[0;31m             \u001b[0mresult\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0m_basic_simps\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    480\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0meven\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'avg'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    481\u001b[0m             \u001b[0mval\u001b[0m \u001b[0;34m/=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.7/site-packages/scipy/integrate/quadrature.py\u001b[0m in \u001b[0;36m_basic_simps\u001b[0;34m(y, start, stop, x, dx, axis)\u001b[0m\n\u001b[1;32m    358\u001b[0m         \u001b[0mh0divh1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh0\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mh1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    359\u001b[0m         tmp = hsum/6.0 * (y[slice0]*(2-1.0/h0divh1) +\n\u001b[0;32m--> 360\u001b[0;31m                           \u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mslice1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mhsum\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mhsum\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mhprod\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    361\u001b[0m                           y[slice2]*(2-h0divh1))\n\u001b[1;32m    362\u001b[0m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "failure_rates = np.zeros((8, 5))\n",
    "failure_sems = np.zeros((8, 5))\n",
    "\n",
    "nIter = 8\n",
    "\n",
    "for r in np.arange(1, 9):\n",
    "\n",
    "    print(\"[\", r, \"]\", sep='', end=\" \")\n",
    "\n",
    "    f_rate_true = np.zeros(nIter)\n",
    "    f_rate_label = np.zeros(nIter)\n",
    "    f_rate_human = np.zeros(nIter)\n",
    "    f_rate_cont = np.zeros(nIter)\n",
    "    f_rate_caus = np.zeros(nIter)\n",
    "\n",
    "    for i in range(nIter):\n",
    "\n",
    "        print(i, end=\" \")\n",
    "\n",
    "        # Create data\n",
    "        df = coinFlipDGWithUnobservables()\n",
    "\n",
    "        # Decider\n",
    "        df_labeled = coinFlipDecider(df,\n",
    "                                     featureX_col=\"X\",\n",
    "                                     featureZ_col=\"Z\",\n",
    "                                     nJudges_M=100,\n",
    "                                     beta_X=1,\n",
    "                                     beta_Z=1,\n",
    "                                     hide_unobserved=True)\n",
    "\n",
    "        df_unlabeled = coinFlipDecider(df,\n",
    "                                       featureX_col=\"X\",\n",
    "                                       featureZ_col=\"Z\",\n",
    "                                       nJudges_M=100,\n",
    "                                       beta_X=1,\n",
    "                                       beta_Z=1,\n",
    "                                       hide_unobserved=False)\n",
    "\n",
    "        # True evaluation\n",
    "\n",
    "        f_rate_true[i] = trueEvaluationEvaluator(df_unlabeled, 'X',\n",
    "                                                 'decision_T', 'result_Y',\n",
    "                                                 r / 10)\n",
    "\n",
    "        # Labeled outcomes only\n",
    "\n",
    "        f_rate_label[i] = labeledOutcomesEvaluator(df_labeled, 'X',\n",
    "                                                   'decision_T', 'result_Y',\n",
    "                                                   r / 10)\n",
    "\n",
    "        # Human evaluation\n",
    "\n",
    "        f_rate_human[i] = humanEvaluationEvaluator(df_labeled, 'judgeID_J',\n",
    "                                                   'decision_T', 'result_Y',\n",
    "                                                   'acceptanceRate_R', r / 10)\n",
    "\n",
    "        # Contraction\n",
    "\n",
    "        f_rate_cont[i] = contractionEvaluator(df_labeled, 'X', 'judgeID_J',\n",
    "                                              'decision_T', 'result_Y',\n",
    "                                              'acceptanceRate_R', r / 10)\n",
    "\n",
    "        # Causal model - empirical performance\n",
    "\n",
    "        f_rate_caus[i] = causalEvaluator(df_labeled, 'X', 'decision_T',\n",
    "                                         'result_Y', r / 10)\n",
    "\n",
    "    failure_rates[r - 1, 0] = np.mean(f_rate_true)\n",
    "    failure_rates[r - 1, 1] = np.mean(f_rate_label)\n",
    "    failure_rates[r - 1, 2] = np.mean(f_rate_human)\n",
    "    failure_rates[r - 1, 3] = np.mean(f_rate_cont)\n",
    "    failure_rates[r - 1, 4] = np.mean(f_rate_caus)\n",
    "\n",
    "    failure_sems[r - 1, 0] = scs.sem(f_rate_true)\n",
    "    failure_sems[r - 1, 1] = scs.sem(f_rate_label)\n",
    "    failure_sems[r - 1, 2] = scs.sem(f_rate_human)\n",
    "    failure_sems[r - 1, 3] = scs.sem(f_rate_cont)\n",
    "    failure_sems[r - 1, 4] = scs.sem(f_rate_caus)\n",
    "\n",
    "x_ax = np.arange(0.1, 0.9, 0.1)\n",
    "\n",
    "plt.errorbar(x_ax,\n",
    "             failure_rates[:, 0],\n",
    "             label='True Evaluation',\n",
    "             c='green',\n",
    "             yerr=failure_sems[:, 0])\n",
    "plt.errorbar(x_ax,\n",
    "             failure_rates[:, 1],\n",
    "             label='Labeled outcomes',\n",
    "             c='magenta',\n",
    "             yerr=failure_sems[:, 1])\n",
    "plt.errorbar(x_ax,\n",
    "             failure_rates[:, 2],\n",
    "             label='Human evaluation',\n",
    "             c='red',\n",
    "             yerr=failure_sems[:, 2])\n",
    "plt.errorbar(x_ax,\n",
    "             failure_rates[:, 3],\n",
    "             label='Contraction, log.',\n",
    "             c='blue',\n",
    "             yerr=failure_sems[:, 3])\n",
    "plt.errorbar(x_ax,\n",
    "             failure_rates[:, 4],\n",
    "             label='Causal model, ep',\n",
    "             c='black',\n",
    "             yerr=failure_sems[:, 4])\n",
    "\n",
    "plt.title('Failure rate vs. Acceptance rate with unobservables')\n",
    "plt.xlabel('Acceptance rate')\n",
    "plt.ylabel('Failure rate')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n",
    "print(failure_rates)\n",
    "print(\"\\nMean absolute errors:\")\n",
    "for i in range(1, failure_rates.shape[1]):\n",
    "    print(np.mean(np.abs(failure_rates[:, 0] - failure_rates[:, i])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,