From 823361089a2bc051fe6214ccb605d020a34105a5 Mon Sep 17 00:00:00 2001
From: Riku-Laine <28960190+Riku-Laine@users.noreply.github.com>
Date: Thu, 14 Mar 2019 10:42:35 +0200
Subject: [PATCH] Lakkaraju v1 implemented, has bugs

---
 Bachelors_thesis_analyses.ipynb | 259 ++++++++++++++++++++------------
 1 file changed, 166 insertions(+), 93 deletions(-)

diff --git a/Bachelors_thesis_analyses.ipynb b/Bachelors_thesis_analyses.ipynb
index 8e4b54e..7635956 100644
--- a/Bachelors_thesis_analyses.ipynb
+++ b/Bachelors_thesis_analyses.ipynb
@@ -4,7 +4,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Bachelors thesis: analyses\n",
+    "# Bachelors thesis' analyses\n",
     "\n",
     "*This Jupyter notebook is for the analyses and model building for Riku Laine's bachelors thesis*\n",
     "\n",
@@ -13,7 +13,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
@@ -22,7 +22,7 @@
        "(7214, 53)"
       ]
      },
-     "execution_count": 19,
+     "execution_count": 1,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -319,27 +319,9 @@
     "                                                                    grid = True)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Implement Lakkaraju\n",
-    "\n",
-    "*Below is an implementation of Lakkaraju's algorithm presented in XYZ (link TBA)*\n",
-    "\n",
-    "* M = number of judges\n",
-    "* subj = number of subjects assigned to each judge\n",
-    "* betas $\\beta$ are coefficients\n",
-    "* R = acceptance rates\n",
-    "* X = invidual's features observable to all (models and judges)\n",
-    "* Z = information observable for judges only\n",
-    "* W = unobservable / inaccessible information\n",
-    "* T = decisions"
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -505,10 +487,10 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.axes._subplots.AxesSubplot at 0x1d7fd04c7b8>"
+       "<matplotlib.axes._subplots.AxesSubplot at 0x2d738932400>"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -531,9 +513,27 @@
     "sns.kdeplot(np.array(compas_raw.age))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generate synthetic data set\n",
+    "\n",
+    "In the chunk below, we generate the synthetic data as described by Lakkaraju et al.\n",
+    "\n",
+    "* M = number of judges\n",
+    "* subj = number of subjects assigned to each judge\n",
+    "* betas $\\beta$ are coefficients\n",
+    "* R = acceptance rates\n",
+    "* X = invidual's features observable to all (models and judges)\n",
+    "* Z = information observable for judges only\n",
+    "* W = unobservable / inaccessible information\n",
+    "* T = decisions"
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 63,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -541,48 +541,48 @@
     "\n",
     "npr.seed(0)\n",
     "\n",
-    "M = 100\n",
-    "subj = 500\n",
+    "nJudges_M = 100\n",
+    "nSubjects_N = 500\n",
     "\n",
     "beta_X = 1.0\n",
     "beta_Z = 1.0\n",
     "beta_W = 0.2\n",
     "\n",
-    "judge_IDs = np.repeat(np.arange(0,M), subj)\n",
+    "judgeID_J = np.repeat(np.arange(0, nJudges_M, dtype = np.int32), nSubjects_N)\n",
     "\n",
-    "acceptance_rates = np.round(npr.uniform(.1, .9, M), 1)\n",
+    "acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 1)\n",
     "\n",
-    "R = np.repeat(acceptance_rates, subj)\n",
+    "acceptanceRate_R = np.repeat(acceptance_rates, nSubjects_N)\n",
     "\n",
-    "X = npr.normal(size = M * subj)\n",
-    "Z = npr.normal(size = M * subj)\n",
-    "W = npr.normal(size = M * subj)\n",
+    "X = npr.normal(size = nJudges_M * nSubjects_N)\n",
+    "Z = npr.normal(size = nJudges_M * nSubjects_N)\n",
+    "W = npr.normal(size = nJudges_M * nSubjects_N)\n",
     "\n",
     "probabilities_Y = 1 / (1 + np.exp(-(beta_X * X + beta_Z * Z + beta_W * W)))\n",
     "\n",
-    "Y = np.round(1 - probabilities_Y).astype(int)\n",
+    "result_Y = np.round(1 - probabilities_Y).astype(int)\n",
     "\n",
     "probabilities_T = 1 / (1 + np.exp(-(beta_X * X + beta_Z * Z)))\n",
-    "probabilities_T += npr.normal(.0, .1, M * subj)\n",
+    "probabilities_T += npr.normal(.0, .1, nJudges_M * nSubjects_N)\n",
     "\n",
-    "T = np.ones(M * subj)*(-1)\n",
+    "decision_T = np.zeros(nJudges_M * nSubjects_N) - 1\n",
     "\n",
-    "tmp = pd.DataFrame(np.column_stack((judge_IDs, R, X, Z, W, Y, probabilities_T, T)),\n",
-    "                   columns = [\"judge_IDs\", \"R\", \"X\", \"Z\", \"W\", \"Y\", \"probabilities_T\", \"T\"])\n",
+    "tmp = pd.DataFrame(np.column_stack((judgeID_J, acceptanceRate_R, X, Z, W, result_Y, probabilities_T, decision_T)),\n",
+    "                   columns = [\"judgeID_J\", \"acceptanceRate_R\", \"X\", \"Z\", \"W\", \"result_Y\", \"probabilities_T\", \"decision_T\"])\n",
     "\n",
     "# Sort by judges then probabilities\n",
-    "df = tmp.sort_values(by = [\"judge_IDs\", \"probabilities_T\"])\n",
+    "df = tmp.sort_values(by = [\"judgeID_J\", \"probabilities_T\"])\n",
     "\n",
     "# Iterate over the data. Subject is in the top (1-r)*100% if\n",
     "# his within-judge-index is over acceptance threshold times\n",
     "# the number of subjects assigned to each judge. If subject\n",
     "# is over the limit they are assigned a zero, else one.\n",
-    "for i in range(M * subj):\n",
-    "    index = i % subj\n",
-    "    if index >= df['R'][i] * subj:\n",
-    "        df['T'][i] = 0\n",
+    "for i in range(nJudges_M * nSubjects_N):\n",
+    "    index = i % nSubjects_N\n",
+    "    if index >= df.acceptanceRate_R[i] * nSubjects_N:\n",
+    "        df.decision_T[i] = 0\n",
     "    else:\n",
-    "        df['T'][i] = 1  # TARKISTA!!!!!!"
+    "        df.decision_T[i] = 1  # TARKISTA!!!!!!"
    ]
   },
   {
@@ -594,7 +594,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 64,
    "metadata": {
     "scrolled": true
    },
@@ -605,7 +605,7 @@
      "text": [
       "0.0    25900\n",
       "1.0    24100\n",
-      "Name: T, dtype: int64\n"
+      "Name: decision_T, dtype: int64\n"
      ]
     },
     {
@@ -628,12 +628,12 @@
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
-       "      <th>T</th>\n",
+       "      <th>decision_T</th>\n",
        "      <th>0.0</th>\n",
        "      <th>1.0</th>\n",
        "    </tr>\n",
        "    <tr>\n",
-       "      <th>Y</th>\n",
+       "      <th>result_Y</th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "    </tr>\n",
@@ -654,21 +654,21 @@
        "</div>"
       ],
       "text/plain": [
-       "T      0.0    1.0\n",
-       "Y                \n",
-       "0.0  13053  12122\n",
-       "1.0  12847  11978"
+       "decision_T    0.0    1.0\n",
+       "result_Y                \n",
+       "0.0         13053  12122\n",
+       "1.0         12847  11978"
       ]
      },
-     "execution_count": 12,
+     "execution_count": 64,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "print(df['T'].value_counts())\n",
+    "print(df.decision_T.value_counts())\n",
     "\n",
-    "tab = df.groupby(['Y', 'T']).size()\n",
+    "tab = df.groupby(['result_Y', 'decision_T']).size()\n",
     "tab.unstack()"
    ]
   },
@@ -681,7 +681,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 65,
    "metadata": {},
    "outputs": [
     {
@@ -689,7 +689,8 @@
      "output_type": "stream",
      "text": [
       "(25000, 8)\n",
-      "(25000, 8)\n"
+      "(25000, 8)\n",
+      "(12134, 8)\n"
      ]
     }
    ],
@@ -698,22 +699,18 @@
     "train, test = np.split(df.sample(frac = 1, random_state = 0), 2)\n",
     "\n",
     "print(train.shape)\n",
-    "print(test.shape)"
+    "print(test.shape)\n",
+    "\n",
+    "train_labeled = train[train.decision_T == 1]\n",
+    "\n",
+    "print(train_labeled.shape)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 66,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "['judge_IDs' 'R' 'X' 'Z' 'W' 'Y' 'probabilities_T' 'T' 'B_prob_1']\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# import the class\n",
     "from sklearn.linear_model import LogisticRegression\n",
@@ -722,41 +719,117 @@
     "logreg = LogisticRegression(random_state=0, solver='lbfgs')\n",
     "\n",
     "# fit, reshape X to be of shape (n_samples, n_features)\n",
-    "logreg.fit(train.X.values.reshape(-1,1), train.Y)\n",
+    "logreg.fit(train_labeled.X.values.reshape(-1,1), train_labeled.result_Y)\n",
     "\n",
     "# predict probabilities and attach to data \n",
     "label_probabilities = logreg.predict_proba(test.X.values.reshape(-1,1))\n",
     "\n",
+    "test['B_prob_0'] = label_probabilities[:, 0]\n",
     "test['B_prob_1'] = label_probabilities[:, 1]\n",
-    "# kaks columnia, ekassa nollan tn tokassa ykkösen\n",
-    "\n"
+    "# kaks columnia, ekassa nollan tn, tokassa ykkösen"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "def lakkaraju(x, j, t, y, s, r):\n",
-    "    return 0\n",
-    "#df = df.sort_values(by = [\"R\", \"judge_IDs\"], ascending = False) # ekana isoin acc rate ja pienin tuomari id\n",
-    "#\n",
-    "#D_q = df.iloc[0:500,] # ekat 500 kuuluu sille\n",
-    "#\n",
-    "#R_q = D_q.iloc[(D_q['T'] == 1).values] # valitaan vapaalle päässeet\n",
-    "#\n",
-    "#R_sort_q = R_q.sort_values(by = [\"probabilities_T\"], ascending = False)\n",
-    "#\n",
-    "#u = np.zeros(10)\n",
-    "## Breikkaa\n",
-    "#for judges_approval in range(10):\n",
-    "#    number_to_remove = np.round((1 - judges_approval/10) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])).astype(int)\n",
-    "#    R_B = R_sort_q.head(number_to_remove)\n",
-    "#\n",
-    "#    u[judges_approval] = np.sum(R_B['Y'] == 0)/D_q.shape[0]\n",
-    "#    \n",
-    "#plt.plot(np.arange(0,1,.1), u);plt.show()"
+    "## Implement Lakkaraju (tba)\n",
+    "\n",
+    "*Below is an implementation of Lakkaraju's algorithm presented in [their paper](https://helka.finna.fi/PrimoRecord/pci.acm3098066).*"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "224\n",
+      "199\n",
+      "174\n",
+      "150\n",
+      "125\n",
+      "100\n",
+      "75\n",
+      "50\n",
+      "26\n",
+      "1\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "[<matplotlib.lines.Line2D at 0x2d73bf4f470>]"
+      ]
+     },
+     "execution_count": 89,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOXZ//HPlclGCHvCvi8ubEIJCG64VVErYEUFN3ADF7St2qpP/T22WJ9aW3ex4oJbRaS2z1O0uCtaUJYgm4BAAJUISFhEFhMIuX5/zKBDDGYCSc5M8n2/XnlxlvtkvoxynTP3uefc5u6IiEjtkBR0ABERqT4q+iIitYiKvohILaKiLyJSi6joi4jUIir6IiK1iIq+iEgtoqIvIlKLqOiLiNQiyUEHKC0rK8vbt28fdAwRkYQyb968Te6eXV67uCv67du3Jzc3N+gYIiIJxcw+j6WdundERGoRFX0RkVokpqJvZoPMbLmZ5ZnZrT/SbpiZuZnlRG27LXLccjM7vTJCi4jIwSm3T9/MQsB44KdAPjDXzKa6+9JS7eoBNwCzo7Z1BYYD3YCWwNtmdpi77628v4KIiMQqliv9fkCeu692993AZGBIGe3uBO4BCqO2DQEmu3uRu68B8iK/T0REAhBL0W8FrI1az49s+46Z9QbauPurFT1WRESqTyxF38rY9t10W2aWBNwP3FTRY6N+x2gzyzWz3IKCghgiiYjIwYhlnH4+0CZqvTWwLmq9HtAdmG5mAM2BqWY2OIZjAXD3x4HHAXJycg5q/sY9e0u4980VNK+fRvMGdWjeIJ0WDdLJykwjlFTWuUdEpPaJpejPBbqYWQfgS8I3Zi/ct9PdtwFZ+9bNbDpws7vnmtm3wCQzu4/wjdwuwJzKi/+9LTt389SM1ezZu/85I5RkZGem0bxBOs3rp4f/jF6O/JmeEqqKWCIicaXcou/uxWY2FngDCAET3X2JmY0Dct196o8cu8TMpgBLgWLguqoaudOsfjrL7zyDLbt2s2FbIRu2FbL+m0K+2lbIhm/C63kFO5iRt4kdRcU/OL5hRsoPTgSlTxIN6qQQ+TQjIpKQzP2gelOqTE5Ojlf1Yxh2FBV/d2IInxC+/e7EEP6ziE07in5wXHpK0n4nhGYN0mnx3YmhDs3rp5NdT91JIlL9zGyeu+eU1y7unr1THTLTkuncNJPOTTMP2GZ3cQkbt0efCPZfzv18K199U3jA7qT9TwjptG+SQeem9WjXJIOUkL4ILSLBqJVFPxapyUm0bpRB60YZB2xTUuL7dSdFf1r46ptwd9LMvE1sj+pOSk4y2mfVpXN25ncnns5NM+mUnUmdVN1XEJGqpaJ/CJKSjKzMNLIy0+jeqsEB220v3MOaTTtZ+dUO8gp2kLdxByu+2s5by75ib0n4k4IZtGpYJ3wSKHVCaJiRWl1/JRGp4VT0q0G99BR6tm5Iz9YN99teVLyXzzfvIm/jjv1OCB+t2kxRccl37bIyU78/CWRn0rlpPTo3zaRZ/TTdWBaRClHRD1BacojDmtXjsGb1oMf32/eWOF9u/Za8gu3kbQyfCFZu3MG/Fqxje+H3XUX10pLp2DSTLvudEDJp0zhDN5NFpEy1cvROonJ3CrYXhU8EkU8F+04IBdu/H22UmpxEx6y6dIo+ITTNpENWXdKSdd9ApCbS6J0ayMxoWj+dpvXTOaZz1n77tn27h7yNO1gVdUJYnL+NaYvXs++8nmTQtnFG+MZx00y6NK3HyUc0pXFd3TMQqS1U9GuIBnVS6NOuEX3aNdpve+GevayKnASiTwjvryhgz14nMy2Z0Sd05IrjOlA3Tf87iNR0+ldew6WnhOjWsgHdWu4/uqh4bwmfbtjOw++u5L63VvDcR5/zi1M6M7xfW32PQKQGU5++8PEXW7n7tU+Zs2YL7ZpkcNNph/OzHi1I0s1gkYQRa5++LumEn7RtxEuj+/P0qL7USQlxw4vzOfuRGXywooB4uygQkUOjoi9A+CbxSUc0ZdoNx3P/BUex7ds9XDpxDhc9OZuFa78OOp6IVBIVfdlPUpJxTu/WvHPTQO44uyufbtjOkPEzue6Fj1ldsCPoeCJyiNSnLz9qR1ExT3ywmif+s5qi4hIu6NuGX57Shab104OOJiJRYu3TV9GXmBRsL+KRd1cyac4XhJKMy4/twJiBnWhQJyXoaCKCir5UkS827+Let5bzrwXraJiRwrUnduLSAe0185hIwCp19I6ZDTKz5WaWZ2a3lrH/ajNbbGYLzGyGmXWNbG9vZt9Gti8ws8cq/leReNK2SQYPDu/Nv284jqNaN+R/pn3KyX+ZzpTctd89MVRE4le5V/pmFgJWAD8lPNH5XGCEuy+NalPf3b+JLA8GrnX3QWbWHnjV3bvHGkhX+onlo1Wbufv1T1m49mu6NM3k16cfzk+7NtPTP0WqWWVe6fcD8tx9tbvvBiYDQ6Ib7Cv4EXUBXfLVEgM6NeH/rj2Gxy7+CXvdGf38PIY99hFz1mwJOpqIlCGWot8KWBu1nh/Zth8zu87MVgH3ADdE7epgZvPN7H0zO/6Q0kpcMjMGdW/Bm788gT/+vAf5W3dx/oSPuPyZuXy64Zvyf4GIVJtYin5Zn9N/cCXv7uPdvRNwC3B7ZPN6oK279wZuBCaZWf0fvIDZaDPLNbPcgoKC2NNLXEkOJTGiX1um33wStww6gtzPtnDGg//hxpcWsHbLrqDjiQixFf18oE3Uemtg3Y+0nwwMBXD3InffHFmeB6wCDit9gLs/7u457p6TnZ0da3aJU3VSQ1xzYic++M1JjD6+I/9evJ5T7n2fca8sZcvO3UHHE6nVYin6c4EuZtbBzFKB4cDU6AZm1iVq9SxgZWR7duRGMGbWEegCrK6M4BL/GmakctuZRzL91ydyTu9WPPPhGk645z0eemclO6MmixeR6lNu0Xf3YmAs8AawDJji7kvMbFxkpA7AWDNbYmYLCHfjjIxsPwFYZGYLgZeBq91dd/hqmRYN6vCnYT1581cncGznJtz31goG/nk6z330Gbuj5gIWkaqnL2dJtfv4i6386bVPmb1mC20bZ3DTaYdxds+WepSzyCHQo5Ulbv2kbSMmj+7P05f1JSM1xC8mL+DsR2bwvh7lLFLlVPQlEGbGSYfv/yjnkXqUs0iVU9GXQJV+lPPyyKOcb/+/xervF6kCKvoSF9KSQ1x2bAfe/81JXHlcB/426wsufGIWG7cXBh1NpEZR0Ze4kpmWzO0/68rDI3qzZN03nP3wDBaou0ek0qjoS1w6+6iW/OOaY0gJJXH+Yx8xJXdt+QeJSLlU9CVudW1Zn1fGHkffDo34zcuLuONfn7Bnr/r5RQ6Fir7EtUZ1U3n2sn5cdXwHnv3ocy56cjabdhQFHUskYanoS9xLDiXx27O68uDwXixc+zWDH57B4vxtQccSSUgq+pIwhvRqxT+uOQYzY9hjH/LPj/ODjiSScFT0JaF0b9WAqWOPpXfbhtw4ZSHjXllKsfr5RWKmoi8Jp0lmGs9fcTSXHdueiTPXcOnEOXpks0iMVPQlIaWEkrjj7G785byjyP18K2c/PIMl69TPL1IeFX1JaMP6tObvYwZQ4s65f/2Qfy34MuhIInFNRV8S3lFtGjJ17HH0aNWAX0xewB+nLWNviZ7WKVIWFX2pEbLrpfHClf25pH87JnywmlFPz+HrXernFylNRV9qjNTkJO4c2p27f96D2au3MPiRmXy64ZugY4nElZiKvpkNMrPlZpZnZreWsf9qM1tsZgvMbIaZdY3ad1vkuOVmdnplhhcpy/B+bZk8pj+Fe/ZyzvgPmbZ4fdCRROJGuUU/MrH5eOAMoCswIrqoR0xy9x7u3gu4B7gvcmxXwhOpdwMGAY/umyhdpCr9pG0jXr3+OI5sUY9rX/iYe17/VP38IsR2pd8PyHP31e6+G5gMDIlu4O7Rn6HrAvv+dQ0BJrt7kbuvAfIiv0+kyjWtn86Lo/szol8bHp2+iiuencu2b/cEHUskULEU/VZA9HNt8yPb9mNm15nZKsJX+jdU8NjRZpZrZrkFBQWxZhcpV1pyiD/+vCd3ndOdmXmbGDp+Jiu/2h50LJHAxFL0rYxtP/ic7O7j3b0TcAtwewWPfdzdc9w9Jzs7O4ZIIhVz0dHtmHRVf7YXFjN0/Exe/2RD0JFEAhFL0c8H2kSttwbW/Uj7ycDQgzxWpMr0bd+YV64/ls7N6nH13+Zx31srKFE/v9QysRT9uUAXM+tgZqmEb8xOjW5gZl2iVs8CVkaWpwLDzSzNzDoAXYA5hx5b5OC0aFCHl0b357w+rXnonZWMfj6XbwrVzy+1R7lF392LgbHAG8AyYIq7LzGzcWY2ONJsrJktMbMFwI3AyMixS4ApwFLgdeA6d99bBX8PkZilp4S4Z1hPxg3pxvTlBQwdP5NVBTuCjiVSLcw9vj7e5uTkeG5ubtAxpJaYtXoz173wMbuLS7j/gl6c2rVZ0JFEDoqZzXP3nPLa6Ru5Uqv179iEqdcfR7usDK58LpeH3lmpfn6p0VT0pdZr1bAOL199DOf0bsV9b63gmhfmsaOoOOhYIlVCRV+EcD//fecfxf/7WVfeXraRc8bPZM2mnUHHEql0KvoiEWbGFcd14LnL+7FpRxGDH5nBe8s3Bh1LpFKp6IuUcmznLKaOPY7WjTK4/Jm5PDo9j3gb8CBysFT0RcrQpnEG/7zmGH7WsyX3vL6csZPms1P9/FIDqOiLHECd1BAPDe/FbWccwWufrOfcv37Iao3nlwSnoi/yI8yMMQM78cxl/Vi/rZBBD/yHP7/xKbt266pfEpOKvkgMTjgsmzd/dQJn9WzB+PdWcfJf3mfqwnXq65eEo6IvEqNm9dO5/4JevHz1AJpkpnLDi/MZ/vgslq3XlIySOFT0RSoop31jpo49jrvO6c7yr7Zz1kP/4Y5/fcK2XXpwm8Q/FX2RgxBKMi46uh3Tbz6Ri45ux/OzPueke6fz4pwvNC2jxDUVfZFD0DAjlTuHdufV64+nc3Ymt/1zMUPHz2Te51uDjiZSJhV9kUrQtWV9XhrTnweH92Lj9kLO/euH3DhlARu3FwYdTWQ/KvoilcTMGNKrFe/edCLXnNiJVxau4+S/vM8TH6xmd3FJ0PFEABV9kUpXNy2ZWwYdwZu/Gkjf9o24a9oyznjwAz5YURB0NJHYir6ZDTKz5WaWZ2a3lrH/RjNbamaLzOwdM2sXtW+vmS2I/EwtfaxITdUhqy5PX9aPp0bmUFziXDpxDqOfy2Xtll1BR5NarNyZs8wsBKwAfkp4ovO5wAh3XxrV5iRgtrvvMrNrgBPd/YLIvh3unhlrIM2cJTVR4Z69PDVjDY+8m0eJO2MGduKagZ2okxoKOprUEJU5c1Y/IM/dV7v7bmAyMCS6gbu/5+77Ll9mAa0rGlikJktPCXHdSZ159+aBnNatOQ+9s5JT73ufaYvX61u9Uq1iKfqtgLVR6/mRbQdyBfBa1Hq6meWa2SwzG3oQGUVqjBYN6vDwiN5MHt2feunJXPvCx1z05GxWfLU96GhSS8RS9K2MbWVempjZxUAO8OeozW0jHzkuBB4ws05lHDc6cmLILSjQzS6p+fp3bMKr1x/H7wd345Mvt3HGg/9h3CtL2fatvtUrVSuWop8PtIlabw2sK93IzE4FfgsMdveifdvdfV3kz9XAdKB36WPd/XF3z3H3nOzs7Ar9BUQSVXIoiZHHtGf6r0/i/Jw2PP3hGk65dzpT5q7V5OxSZWIp+nOBLmbWwcxSgeHAfqNwzKw3MIFwwd8Ytb2RmaVFlrOAY4GliMh3GtdN5Y8/78HU646jbeMMfvOPRZzz1w9ZsPbroKNJDVRu0Xf3YmAs8AawDJji7kvMbJyZDY40+zOQCfy91NDMI4FcM1sIvAfcHT3qR0S+16N1A16++hjuPe8ovtz6LUPHz+Q3Ly9k046i8g8WiVG5Qzarm4ZsisD2wj08/G4eE2esoU5qiF+dehiXDGhHSkjfp5SyVeaQTRGpZvXSU/ivM4/k9V+eQK82DRn36lLOeug/fJi3KehokuBU9EXiWOemmTx3eT8mXNKHXbv3cuGTs7n2hXl8+fW3QUeTBKWiLxLnzIzTuzXn7RsH8qtTD+OdZRs55d7pPPTOSgr37A06niQYFX2RBJGeEuIXp3bhnZsGcvIRTbnvrRX89P73eWPJBn2rV2Kmoi+SYFo3yuDRi/rwwpVHk54cYszz8xj74nxd9UtMVPRFEtSxnbOY9ovj+fXphzNt8XqGPz6Lgu0a3ik/TkVfJIGlhJK47qTO/PWiPny64RvOeXQmK/UcH/kRKvoiNcCg7s15afQAiopL+PmjHzJjpYZ2StlU9EVqiKPaNOT/rjuWlg3rMOrpOUye80XQkSQOqeiL1CCtGtbh5WsGcEznLG7952Lufu1TPbxN9qOiL1LD1EtPYeLIHC48ui2Pvb+KsS9+rJE98h0VfZEaKDmUxF1Du/PbM4/ktU82aGSPfEdFX6SGMjOuOqEjj12skT3yPRV9kRru9G7NmTJGI3skTEVfpBbo2VojeyRMRV+kltDIHoEYi76ZDTKz5WaWZ2a3lrH/RjNbamaLzOwdM2sXtW+kma2M/IyszPAiUjH7RvZcpJE9tVa5Rd/MQsB44AygKzDCzLqWajYfyHH3nsDLwD2RYxsDdwBHA/2AO8ysUeXFF5GKSg4l8Yeh3bn9rPDIngs0sqdWieVKvx+Q5+6r3X03MBkYEt3A3d9z912R1VlA68jy6cBb7r7F3bcCbwGDKie6iBwsM+PK48Mje5Zv+Iah4zWyp7aIpei3AtZGredHth3IFcBrB3msiFSjfSN7du/VyJ7aIpaib2VsK/Puj5ldDOQAf67IsWY22sxyzSy3oKAghkgiUln2jexp1Ugje2qDWIp+PtAmar01sK50IzM7FfgtMNjdiypyrLs/7u457p6TnZ0da3YRqSStGtbh71cP4FiN7KnxYin6c4EuZtbBzFKB4cDU6AZm1huYQLjgb4za9QZwmpk1itzAPS2yTUTiTL30FJ7SyJ4aL7m8Bu5ebGZjCRfrEDDR3ZeY2Tgg192nEu7OyQT+bmYAX7j7YHffYmZ3Ej5xAIxz9y1V8jcRkUO2b2RPh6y63DVtGV9+PYsnL80hu15a0NGkkli8Taick5Pjubm5QccQqfXeWLKBX05eQOO6qTx9WV8Oa1Yv6EjyI8xsnrvnlNdO38gVkTKd3q05L43pz+69JZyrkT01hoq+iByQRvbUPCr6IvKjSo/s+eNryzSyJ4Gp6ItIufaN7Lm4f1smvL+a6yZpZE+iUtEXkZgkh5K4c0j4mT2vL9EzexKVir6IxCz6mT0rNmxn6PiZrNAzexKKir6IVFj0M3s0siexqOiLyEHp0bqBRvYkIBV9ETloGtmTeFT0ReSQlDWy59vdGtkTr1T0ReSQlR7ZM/wJjeyJVyr6IlIp9o3smaCRPXFNRV9EKtVpkZE9e/aWcP6Ej1icvy3oSBJFRV9EKl2P1g14+epjyExL5sInZzHv861BR5IIFX0RqRJtm2QwZcwAmtRN5ZKnZvPRqs1BRxJU9EWkCrVsWIcpYwbQqmF4LP8HKzQHdtBiKvpmNsjMlptZnpndWsb+E8zsYzMrNrNhpfbtNbMFkZ+ppY8VkZqtaf10Jo/uT8fsTK58Npe3l34VdKRardyib2YhYDxwBtAVGGFmXUs1+wIYBUwq41d86+69Ij+DDzGviCSgJplpTL6qP0e2qMfVf5vHvxetDzpSrRXLlX4/IM/dV7v7bmAyMCS6gbt/5u6LgJIqyCgiNUCDjBT+duXR9G7bkOtf/Jj/nZ8fdKRaKZai3wpYG7WeH9kWq3QzyzWzWWY2tELpRKRGqZeewrOX96N/xybcOGUhL+p5PdUulqJvZWyryMM12kYm670QeMDMOv3gBcxGR04MuQUFutEjUpNlpCYzcVRfBh6WzW3/XMwzM9cEHalWiaXo5wNtotZbA+tifQF3Xxf5czUwHehdRpvH3T3H3XOys7Nj/dUikqDSU0JMuKQPp3Vtxu9eWcpj768KOlKtEUvRnwt0MbMOZpYKDAdiGoVjZo3MLC2ynAUcCyw92LAiUnOkJYcYf9FPOPuoltz92qc88PYK3PWEzqqWXF4Ddy82s7HAG0AImOjuS8xsHJDr7lPNrC/wv0Aj4Gwz+727dwOOBCaYWQnhE8zd7q6iLyIApISSeOCCXqQlJ/HA2ysp3FPCLYMOx6ysXmWpDOUWfQB3nwZMK7Xtv6OW5xLu9il93IdAj0PMKCI1WCjJuOfcnqSnJPHY+6so3LOXO87uqsJfRWIq+iIiVSkpybhzSHfSkkM8NWMNRcUl3DW0O0lJKvyVTUVfROKCmXH7WUeSnpLE+PdWUbRnL/cM60lySE+LqUwq+iISN8yMX59+BOnJIe59awVFxSU8MLwXKSr8lUZFX0TizvWndCE9JcRd05ZRVFzC+It6k5YcCjpWjaDTp4jEpatO6MidQ7rx9rKvuPLZXM27W0lU9EUkbl0yoD33nNuTGXmbuOyZOewsKg46UsJT0ReRuHZ+3zY8cEEv5n62lUuems03hXuCjpTQVPRFJO4N6dWKR0b0ZvGX27joidls3bk76EgJS0VfRBLCGT1aMOGSPiz/ajsjnpjFph1FQUdKSCr6IpIwTj6iGRNH9uWzzTu5YMJHbNhWGHSkhKOiLyIJ5bguWTx3+dFs2FbIBY9/RP7WXUFHSigq+iKScPp1aMzfrjyarTt3c8GEWXy2aWfQkRKGir6IJKTebRsx6ar+7NpdzPkTPiJv4/agIyUEFX0RSVjdWzVg8ugBlDhcMGEWy9Z/E3SkuKeiLyIJ7fDm9Zgypj8poSRGPDGLRflfBx0prqnoi0jC65idyZQxA8hMS+aiJ2Yz7/MtQUeKWyr6IlIjtG2SwZQxA8iql8YlT83ho1Wbg44Ul2Iq+mY2yMyWm1memd1axv4TzOxjMys2s2Gl9o00s5WRn5GVFVxEpLSWDevw0uj+tGpYh1FPz+H9FQVBR4o75RZ9MwsB44EzgK7ACDPrWqrZF8AoYFKpYxsDdwBHA/2AO8ys0aHHFhEpW9P66Uwe3Z9O2Zlc9Wwuby7ZEHSkuBLLlX4/IM/dV7v7bmAyMCS6gbt/5u6LgJJSx54OvOXuW9x9K/AWMKgScouIHFCTzDRevKo/R7asz7UvfMyri9YFHSluxFL0WwFro9bzI9tiEdOxZjbazHLNLLegQB/HROTQNchI4W9X9KN324bc8OJ8/jEvP+hIcSGWol/WzMQe4++P6Vh3f9zdc9w9Jzs7O8ZfLSLy4+qlp/Ds5f3o37EJN7+8kEmzvwg6UuBiKfr5QJuo9dZArJ+VDuVYEZFDlpGazMRRfTnxsGz+638X8/TMNUFHClQsRX8u0MXMOphZKjAcmBrj738DOM3MGkVu4J4W2SYiUm3SU0I8dkkfTu/WjN+/spQnPlgddKTAlFv03b0YGEu4WC8Dprj7EjMbZ2aDAcysr5nlA+cBE8xsSeTYLcCdhE8cc4FxkW0iItUqLTnEIxf+hLN6tOCuacv496L1QUcKhLnH2j1fPXJycjw3NzfoGCJSQxXu2ctFT87mky+38dKYAfRq0zDoSJXCzOa5e0557fSNXBGpVdJTQky4pA/Z9dK48tlcvvz626AjVSsVfRGpdbIy03h6VF+K9uzlimfmsr0WTbauoi8itVKXZvV49OKfsHLjDm54cT7Fe0t/t7RmUtEXkVrr+C7Z/H5wN95bXsAf/r0s6DjVIjnoACIiQbq4fztWF+xk4sw1dMyuy6UD2gcdqUqp6ItIrffbs47k8807+d3UJbRtnMGJhzcNOlKVUfeOiNR6oSTjoRG9Obx5fcZOms/yDTV3vl0VfRERoG5aMk+NzCEjNcTlz8ylYHtR0JGqhIq+iEhEy4Z1eHJkDpt3FjH6+VwK9+wNOlKlU9EXEYnSs3VDHrigF/O/+Jqb/76QkpL4emrBoVLRFxEpZVD3Ftwy6AheXbSeB95eEXScSqXROyIiZbh6YEfWbNrBQ+/m0SG7Luf0bh10pEqhK30RkTKYGX8Y2oP+HRtzy8uLmftZzXhAsIq+iMgBpCYn8djFfWjVqA5jnp/H55t3Bh3pkKnoi4j8iIYZqUwc1ZcSdy5/Zi7bvk3sh7Op6IuIlKNDVl0eu7gPX2zZxbUvzGNPAj+cLaaib2aDzGy5meWZ2a1l7E8zs5ci+2ebWfvI9vZm9q2ZLYj8PFa58UVEqkf/jk34n3N6MDNvM//9ryXE2wRUsSp39I6ZhYDxwE8JT3Q+18ymuvvSqGZXAFvdvbOZDQf+BFwQ2bfK3XtVcm4RkWp3Xk4b1mzayaPTV9Epuy5XHt8x6EgVFsuVfj8gz91Xu/tuYDIwpFSbIcCzkeWXgVPMzCovpohIfLj5tMM5s0dz7pq2jDeXbAg6ToXFUvRbAWuj1vMj28psE5lIfRvQJLKvg5nNN7P3zez4Q8wrIhKopCTj3vN60bNVA34xeQGffLkt6EgVEkvRL+uKvXRn1oHarAfauntv4EZgkpnV/8ELmI02s1wzyy0oKIghkohIcOqkhnji0hwaZaRw5bO5bNhWGHSkmMVS9POBNlHrrYF1B2pjZslAA2CLuxe5+2YAd58HrAIOK/0C7v64u+e4e052dnbF/xYiItWsaf10nhrVl+2Fe7jyubns2l0cdKSYxFL05wJdzKyDmaUCw4GppdpMBUZGlocB77q7m1l25EYwZtYR6AKsrpzoIiLBOrJFfR6+sDdL133DLycvSIiHs5Vb9CN99GOBN4BlwBR3X2Jm48xscKTZU0ATM8sj3I2zb1jnCcAiM1tI+Abv1e5eM77LLCICnHxEM24/qytvLv2KP73+adBxyhXTA9fcfRowrdS2/45aLgTOK+O4fwD/OMSMIiJx7bJj27Nm004mfLCaDll1Gd6vbdCRDkhP2RQROURmxh1nd+XzLbu4/f8+oW3jDI7pnBV0rDLpMQwiIpUgOZTEIxf2pkNWXa7+2zzyNu4IOlKZVPRFRCpJ/fQUJo7qS0ooiSuencuWnbuDjvQDKvoiIpWoTeMMHr+pvC9RAAAHRElEQVQ0h/XbCrn6+XkUFcfXPLsq+iIilaxPu0b85byjmPPZFm775+K4ejibbuSKiFSBwUe1ZE3BTu5/ewUds+oy9uQuQUcCVPRFRKrMDad0Zs2mHfzlzRW0z6rLz3q2DDqSundERKqKmXH3uT3JadeIm6YsZP4XW4OOpKIvIlKV0lNCTLikD83qp3PVc/PI37or0Dwq+iIiVaxJZhoTR+VQVLyXK57JZXthcPPsquiLiFSDzk3r8deL+pBXsIPrX5xPcUDz7Kroi4hUk+O6ZHHnkO5MX17AH/69LJAMGr0jIlKNLjy6LasLdvDkjDV0zK7LpQPaV+vrq+iLiFSz2848ks827+J3U5fQtnEGJx7etNpeW907IiLVLJRkPDi8F0c0r8/YSfNZvmF7tb22ir6ISADqpiXz1KgcMlJDXP7MXAq2F1XL68ZU9M1skJktN7M8M7u1jP1pZvZSZP9sM2sfte+2yPblZnZ65UUXEUlsLRrU4amRfdmyczdXPZdL4Z6qfzhbuUU/MsfteOAMoCswwsy6lmp2BbDV3TsD9wN/ihzblfCcut2AQcCj++bMFRER6NG6Afdf0IuF+V9z898XVvk8u7Fc6fcD8tx9tbvvBiYDQ0q1GQI8G1l+GTjFzCyyfbK7F7n7GiAv8vtERCRiUPfm3DroCDpmZ2JWta8Vy+idVsDaqPV84OgDtXH3YjPbBjSJbJ9V6thWB51WRKSGGjOwU7W8TixX+mWdd0p//jhQm1iOxcxGm1mumeUWFBTEEElERA5GLEU/H2gTtd4aWHegNmaWDDQAtsR4LO7+uLvnuHtOdnZ27OlFRKRCYin6c4EuZtbBzFIJ35idWqrNVGBkZHkY8K6Hp4qZCgyPjO7pAHQB5lROdBERqahy+/QjffRjgTeAEDDR3ZeY2Tgg192nAk8Bz5tZHuEr/OGRY5eY2RRgKVAMXOfu8TVhpIhILWLxNHcjQE5Ojufm5gYdQ0QkoZjZPHfPKa+dvpErIlKLqOiLiNQiKvoiIrVI3PXpm1kB8Pkh/IosYFMlxalMylUxylUxylUxNTFXO3cvd8x73BX9Q2VmubHczKhuylUxylUxylUxtTmXundERGoRFX0RkVqkJhb9x4MOcADKVTHKVTHKVTG1NleN69MXEZEDq4lX+iIicgAJWfRjmL7xBDP72MyKzWxYHOW60cyWmtkiM3vHzNrFUbarzWyxmS0wsxllzI4WSK6odsPMzM2sWkZcxPB+jTKzgsj7tcDMroyHXJE250f+P1tiZpPiIZeZ3R/1Xq0ws6/jJFdbM3vPzOZH/l2eGSe52kVqxCIzm25mrSvtxd09oX4IP/RtFdARSAUWAl1LtWkP9ASeA4bFUa6TgIzI8jXAS3GUrX7U8mDg9XjIFWlXD/iA8IQ8OfGQCxgFPFId//0qmKsLMB9oFFlvGg+5SrW/nvCDGwPPRbgP/ZrIclfgszjJ9XdgZGT5ZOD5ynr9RLzSL3f6Rnf/zN0XASVxlus9d98VWZ1FeH6BeMn2TdRqXcqY7CaIXBF3AvcAhdWQqSK5qlssua4Cxrv7VgB33xgnuaKNAF6Mk1wO1I8sN6CM+T4CytUVeCey/F4Z+w9aIhb9sqZvjIcpGCua6wrgtSpN9L2YspnZdWa2inCBvSEecplZb6CNu79aDXlizhVxbuTj98tm1qaM/UHkOgw4zMxmmtksMxsUJ7mAcLcF0AF4N05y/Q642MzygWmEP4XEQ66FwLmR5XOAembWpDJePBGLfkxTMAYg5lxmdjGQA/y5ShNFvWQZ236Qzd3Hu3sn4Bbg9ipPVU4uM0sC7gduqoYs0WJ5v14B2rt7T+Bt4NkqTxVbrmTCXTwnEr6iftLMGsZBrn2GAy979cyrEUuuEcAz7t4aOJPwvCBVXRdjyXUzMNDM5gMDgS8Jz0lyyBKx6Mc0BWMAYsplZqcCvwUGu3tRPGWLMhkYWqWJwsrLVQ/oDkw3s8+A/sDUariZW+775e6bo/77PQH0qeJMMeWKtPmXu+9x9zXAcsIngaBz7TOc6unagdhyXQFMAXD3j4B0ws+/CTSXu69z95+7e2/C9QJ331Ypr17VNy2q4CZIMrCa8EfEfTdBuh2g7TNU343ccnMBvQnfwOkSb+9ZdCbgbMKzogWeq1T76VTPjdxY3q8WUcvnALPiJNcg4NnIchbhboQmQeeKtDsc+IzI94Pi5P16DRgVWT6ScPGt0nwx5soCkiLLdwHjKu31q+PNr4I37UxgRaSA/jaybRzhq2eAvoTPpjuBzcCSOMn1NvAVsCDyMzWO3rMHgSWRXO/9WPGtzlyl2lZL0Y/x/fpj5P1aGHm/joiTXAbcR3iK0sXA8HjIFVn/HXB3deSpwPvVFZgZ+e+4ADgtTnINA1ZG2jwJpFXWa+sbuSIitUgi9umLiMhBUtEXEalFVPRFRGoRFX0RkVpERV9EpBZR0RcRqUVU9EVEahEVfRGRWuT/AxHzVXtBLQrTAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "def contraction(df, j_name, t_name, y_name, s_name, r_name, 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",
+    "    Parameters:\n",
+    "    df = The (Pandas) data frame containing the data, judge decisions,\n",
+    "    judge IDs, results and probability scores.\n",
+    "    j_name = String, the name of the column containing the judges' IDs\n",
+    "    in df.\n",
+    "    t_name = String, the name of the column containing the judges' decisions\n",
+    "    s_name = String, the name of the column containing the probability\n",
+    "    scores from the black-box model B.\n",
+    "    r_name = 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",
+    "    u = The estimated failure rate at acceptance rate r.\n",
+    "    '''\n",
+    "    # Sort first by acceptance rate and judge ID.\n",
+    "    sorted_df = df.sort_values(by = [r_name, j_name], ascending = False)\n",
+    "\n",
+    "    most_lenient_ID = int(sorted_df[j_name].head(1)) # \"hot mess\"\n",
+    "\n",
+    "    D_q = sorted_df[sorted_df[j_name] == most_lenient_ID]\n",
+    "\n",
+    "    R_q = D_q[D_q[t_name] == 1]\n",
+    "    \n",
+    "    R_sort_q = R_q.sort_values(by = s_name, ascending = False)\n",
+    "    \n",
+    "    number_to_remove = int(np.round((1 - r) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])))\n",
+    "    print(number_to_remove)\n",
+    "    R_B = R_sort_q.head(number_to_remove)\n",
+    "    \n",
+    "    return 1 / D_q.shape[0] * np.sum(R_B[y_name] == 0) \n",
+    "\n",
+    "failure_rates = np.zeros(10)\n",
+    "\n",
+    "for r in range(10):\n",
+    "    failure_rates[r] = contraction(test, 'judgeID_J', 'decision_T',\n",
+    "                                   'result_Y', 'B_prob_0', 'acceptanceRate_R',  r / 10)\n",
+    "    \n",
+    "plt.plot(np.arange(0.1,1,.1), failure_rates[1:])\n"
    ]
   }
  ],
-- 
GitLab