{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Bias Mitigation with Disparate Impact Remover\n", "\n", "This notebook is an example of bias mitigation with an pre-processing algorithm called Disparate Impact Remover. \n", "\n", "The library we are using to implement this algorithm is [AI Fairness 360](http://aif360.mybluemix.net/). \n", "\n", "The algorithm was introduced in the paper [\"Certifying and removing disparate impact\" by M. Feldman, S. A. Friedler, J. Moeller, C. Scheidegger, and S. Venkatasubramanian](https://arxiv.org/abs/1412.3756)\n", "\n", "### Data\n", "Import required libraries" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from matplotlib import pyplot as plt\n", "\n", "import numpy as np\n", "from numpy.random import choice as np_choice\n", "import pandas as pd\n", "from tqdm import tqdm\n", "import seaborn as sns\n", "from IPython.display import Markdown, display\n", "\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.linear_model import LogisticRegression\n", "from sklearn.svm import SVC as SVM\n", "from sklearn.preprocessing import MinMaxScaler\n", "from sklearn.decomposition import PCA\n", "\n", "from aif360.algorithms.preprocessing import DisparateImpactRemover\n", "from aif360.datasets import BinaryLabelDataset\n", "from aif360.metrics import BinaryLabelDatasetMetric" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are using a random seed so that the comments remain correct when rerunning this notebook" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "seed = 123\n", "np.random.seed(seed)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "We are going to create two unequal datasets\n", "- The first is an unprivileged group where their mean is 5.5 and standard deviation 0.5. The `output` variable indicates if the individual received a favourable outcome (1) or unfavourable (0), the probability for a favourable outcome for this group is 0.5.\n", "- The second a priviledged group, their mean is higher, 6, and the standard deviation 0.4. The probability for a favourable outcome for this group is 0.7." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [], "source": [ "unpriv = pd.DataFrame({'group':[0 for i in range(0, 400)],\n", " 'value':np.random.normal(5.5, 0.6, 400), \n", " 'output':np.random.choice([0, 1], size=(400), p=[0.5, 0.5])})" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": true }, "outputs": [], "source": [ "priv = pd.DataFrame({'group':[1 for i in range(0, 1000)],\n", " 'value':np.random.normal(6, 0.4, 1000),\n", " 'output':np.random.choice([0, 1], size=(1000), p=[0.3, 0.7])})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can plot the distributions on a graph to show the disparity between the values for the two groups" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEKCAYAAADpfBXhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd8XNWZ+P/PmVEvo96LJcu994JpoUMIJhsCmCUJJAGSQHaz2Wx+5Le7SZYN2WyyCZtCCkswCQQIkAbElAA2YIyL3Lsty+pW7xpJ0873jysTocjWSJo7dzTzvF8vvaSZuXPPg7AeHT333OcorTVCCCHCi83qAIQQQgSeJHchhAhDktyFECIMSXIXQogwJMldCCHCkCR3IYQIQ5LchRAiDElyF0KIMCTJXQghwlCUVQNnZmbqkpISq4YXQogpaffu3a1a66yxjrMsuZeUlFBeXm7V8EIIMSUppar9OU7KMkIIEYYkuQshRBiS5C6EEGFIkrsQQoQhSe5CCBGGJLkLIUQYkuQuhBBhSJK7EEKEIUnuQggRhiy7Q1UIEQDlG8/92oo7gxeHCDljztyVUo8ppZqVUofO8bpSSv1IKVWhlDqglFoW+DCFEEKMhz9lmceBa87z+rXAzKGPu4GfTT4sIYQQkzFmctdavw20n+eQ9cCvtWE7kKqUygtUgEIIIcYvEBdUC4DaYY/rhp77G0qpu5VS5Uqp8paWlgAMLYQQYjSBSO5qlOf0aAdqrR/RWq/QWq/IyhqzHbEQQogJCkRyrwOKhj0uBBoCcF4hhBATFIjk/gLwyaFVM2uALq31mQCcVwghxASNuc5dKfU0cCmQqZSqA74BRANorX8ObAKuAyoAJyCLa4UQwmJjJnet9YYxXtfAvQGLSIhIIzciCRNI+wEhhAhDktyFECIMSXIXQogwJMldCCHCkCR3IYQIQ5LchRAiDElyFyIUDPZARxW4B6yORIQJ2axDCCtVvA7vfB+6hvXeS8qF/CVQsHz092gN/R3Gh/aBkjma+FuS3IWwysnX4ZkNEJcCs68zknpvI7SegBOvwolXYO8TULgSouLB64LW49ByAjz9xjliHZC70Hh/TKK1/z0ipEhyFyIYRt6F2lYBO35uJPQ1X4CYhKEXFsHMq4xZeeNBI6HXlYPPA8oOmTOMu1ZTiiA6HnY/DrXboavOOE9UbLD/y0SIkuQuRLD5vHDgWYhPgzWfH5bYh4lPg9KL/Ws/0HjA+OVR/hisugts8mMt5IKqEMFXtwv6mmHuRwJTSsldBItuMUo2VVsnfz4RFiS5CxFMXrdRS0+dBjkLA3fe4jWQXgaVW4wSjoh4ktyFCKbqd2GgE+ZcD2q0TcwmYcYVxrnrygN7XjElSXIXIli0huptkFYKmTMDf/6sOeAohFNvGEskRUST5C5EsHTWGLX2otXmnF8pY/be1wJNh80ZQ0wZktyFCJa6XWCLhrzF5o2RuxCiE+DMfvPGEFOCJHchgsHngYY9kLvAWJ9uFpsdcuZD8xHj4q2IWJLchQiG5iPgdhp3m5otd6ExVvU288cSIUuSuxDBUL8HYpIgc7b5Y2XNMco/x14yfywRsiS5C2E2rxtajkHOAqNsYjZ7DGTNhmN/NlboiIgk9ykLYbbaHeAZgOy543/vyJ40/spdBPufgoa9ULBsYucQU5rM3IUw28nXhpp+BaEkc1bOPOPzqTeDN6YIKZLchTDbidcgfTpExwVvzJgkyJ4nF1UjmCR3IczUWQMtR/86kw6maRcYJSGv9JqJRJLchTDTydeMz9kWJXdXr9ESWEQcSe5CmKniDUgthsTs4I9dfIHxWUozEUlWywhhFp8Xqt6F+TcGvgOkP068AgmZsP/pv+0b788mIGJKk5m7EGY5sx8Gu4wdlaySUQbtldIlMgLJzF2IQBm5Jv3sMsSuOmMTbCuklxkXVXsawZFvTQzCEjJzF8IsrScgKdu6xA7GzB2M2buIKH4ld6XUNUqp40qpCqXU/aO8XqyU2qyU2quUOqCUui7woQoxhfi8RkLNMGFTjvGITzfWvHdWWxuHCLoxk7tSyg48DFwLzAM2KKVGruv6N+BZrfVS4Fbgp4EOVIgppasGvC7rk7tSxn6tnTXWxiGCzp+Z+yqgQmtdqbV2Ac8A60ccowHH0NcpQEPgQhRiCmqtMD5nzrA2DjCWYvY2g7vf6khEEPmT3AuA2mGP64aeG+6bwO1KqTpgE/DFgEQnxFTVVgHJeUZJxGpp0wANXbVjHirChz/JfbQFuiP7iG4AHtdaFwLXAU8opf7m3Eqpu5VS5Uqp8paWlvFHK8RU4PNCx2ljpUooSCkyPktpJqL4k9zrgKJhjwv527LLZ4BnAbTW7wFxQObIE2mtH9Far9Bar8jKyppYxEKEuu76oXp7iCT3mERIzIIOuagaSfxJ7ruAmUqpUqVUDMYF0xdGHFMDXA6glJqLkdxlai4iU/sp43P6dGvjGC61WGbuEWbM5K619gD3Aa8CRzFWxRxWSj2glLph6LB/Bu5SSu0Hngbu0Fq2gBERqu2Ucdu/levbR0otNu6W7e+0OhIRJH7doaq13oRxoXT4c18f9vURYF1gQxNiCtI+aD8NOfOtjuSDUqcZnztrID7V2lhEUMgdqkIEUm8zuPtCp95+lqMAlE1uZoogktyFCKRQrLcD2KMhOd/ocyMigiR3IQKpvRJiHUbNPdSkFhlr3eVyWESQ5C5EILWdMmbtVvRvH0tKEbid4GyzOhIRBJLchQgUZzsMdIZeSeas1KHbVeRO1YggyV2IQDnbVjdUk3tyHtjs0CnJPRJIchciUNpPQVRc6G6KYYuC5AKZuUcISe5CBEp7JaSVGksOQ9XZi6o+2XYv3IXwv0IhppC+NuhtgowQLcmclVIEngHZmSkCSHIXIhBq3jM+h0onyHNJLTY+N+y1Ng5hOknuQgRCzXtGTTul2OpIzi8pB2zRUL/b6kiEySS5CxEI1duMWbHdr3ZN1rHZjbq7JPewJ8ldiMka7IEz+0O/JHNWarERr9dtdSTCRJLchZis2h2gvZARAvul+iO1GLyD0HTY6kiEiSS5CzFZVe8a9fa0Eqsj8c/Z9r/15dbGIUwlyV2IyareBvlLISrW6kj8E59uNDar32N1JMJEktyFmAyX07g4Oe0CqyPxn1JQsBzqZOYeziS5CzEZdbvA54ZpF1odyfgUroDWEzDQZXUkwiSS3IWYjOptRruB4tVWRzI+BcsADQ37rI5EmESSuxCTUf0u5C4Mrc2w/ZG/zPhct8vaOIRpJLkLMVHuAajdCSUXWR3J+CWkQ+YsYxmnCEuS3IWYqLqdxnrxqZjcAYrXQs0O8HmtjkSYQJK7EBN1+h2j3j5trdWRTEzxWhjsguYjVkciTCDJXYiJqnoH8pZMvXr7WWd/KdVstzYOYQpJ7kJMhMtprBMvnaIlGTDuVE3ON1b8iLAjyV2IiajdbqxvL7nY6kgmTilj9l7zHmhtdTQiwCS5CzERp982+skUr7E6kskpXgs9Z6CjyupIRICFePNpIYLvqR01Yx5z1cE30I4F/GVfO9AOQFlN+/uvry5NNyu8wDrbNqHmPUgvtTYWEVAycxdinKLdPaR3HaIpY6XVoUxe1lyIT4OqrVZHIgJMkrsQ45TdXo4NH40ZU3QJ5HA2G0y/FCrekLp7mJHkLsQ45bbtwGOLozV1sdWhBMaMK6C3UTbvCDOS3IUYp5y2HbSkLcVnj7E6lMAou8z4fOoNa+MQAeVXcldKXaOUOq6UqlBK3X+OY25WSh1RSh1WSj0V2DCFCA1xg62k9lbQmDHFukCejyMfsudDxetWRyICaMzkrpSyAw8D1wLzgA1KqXkjjpkJfA1Yp7WeD3zJhFiFsFxOm9Foqylzii+BHGnG5VD9Hgz2Wh2JCBB/Zu6rgAqtdaXW2gU8A6wfccxdwMNa6w4ArXVzYMMUIjTktu3AFZVMh2OO1aEE1ozLjZuyZNVM2PAnuRcAtcMe1w09N9wsYJZS6l2l1Hal1DWjnUgpdbdSqlwpVd7S0jKxiIWwitbktO2gKWMVWtmtjiawitdCdAKcfM3qSESA+HMTkxrluZFrpqKAmcClQCHwjlJqgda68wNv0voR4BGAFStWyLorEXIG3F4qW3qpbO2jpWeQ7n433QNu3F7NHKq4TTXwlPsSKna8yczEfmYkDpBg91kd9uRFxRqrZo69BDnzjW6X57LizuDFJSbMn+ReBxQNe1wINIxyzHattRs4rZQ6jpHsZZsXEfKcLg+vHm7k93vq2XG6HZfHhwLSE2NIiY9mWkYi0XYbV7ftAzds9izkvYYsAKKUj+UpvVya2cViRx/20aZCU8X8G+HoC9B+GjLKrI5GTJI/yX0XMFMpVQrUA7cCt4045o/ABuBxpVQmRpmmMpCBChFoFc09PPrOaV7c30Cfy0thWjyfXDMNl8fHtIxE4mM+WHq5+t19uHQyX5rVz13eE5xyxrGnK4mt7Q52dDpIi3ZzY24bl2dO0U2nZ14NUXFwZp8k9zAwZnLXWnuUUvcBrwJ24DGt9WGl1ANAudb6haHXrlJKHQG8wL9ordvMDFyIidpd3c7Ptpzi9aPNxEXbWL+4gI8tL2TFtDRsNjV6bxntw9FXRWfSDFCKxCgfixxOFjmc3F7QzN7uJP7clM7G2lxeaMzgX2z9fLxkYGrN5GOTjNLM6Xdg/kfPX5oRIc+vxmFa603AphHPfX3Y1xr48tCHECHp6JluvvvKMTYfbyEtIZovXTGTT64tIT1x7JuRUntOEu110p30t821omywMrWXFSm9HOxJ4LcNWdy/28Ezp+P5r2XdzE2dQtvYzbvRqLt3VEH6dKujEZMgXSFF2GvsGuC/XznGH/fVkxwbxf3XzuFTa0v+puxyPrltxm5F3Ynn7pyoFCxyOFmYXE1jVB7/uT+Zj7yRzt2znPzjvD5ip8ICm9nXGK2MG/ZJcp/iJLmLsOX2+vjVtioe+ssJ3D7NPReX8flLykhJiB73uXLadtAfk4Er2jHmsUrBjcWDXJLj4tsHk/jp8UTeborhx6u7KU0O8Vl8bLJxt2rDHpi3HmxT4TeSGI0kdxF2ntpRQ31HP8/vqaWpe5DZOclcvyiPjKRY/nzwzLjPp3xustt30+6YO673pcVqvreih6vyB/mXcgfXv5HGg8t6uLF4cNwxBFXRSmjcDy1HIWeB1dGICZIrJiKsuL0+3jjWxM/eqqDf5eX21dP45NppZCTFTvicGV2HjXr7eUoy53NlvotNV7QzL8XDl3am8I29SXhCeWl81lyISYLanVZHIiZBZu4ibNS2O7nvqT3sr+tiSVEqH1mUP666+rlktxu3a3QnTpvwOfITfDx9SSffOZjEoycTqOiJ4uE1XaTGhOC9fDY7FCw3WhG4+iAm0eqIxATIzF2Eha0nW/nIT7ZyurWPDauKuXlFUUASO0BO2046kmfiiZpckouywb8t7uW7K7rZ1RrNjW+mUdkTojXtolWgvUbtXUxJMnMXIcefPUzP0lqztaKVVw41kpUcyyfWTa4EM5LN6yKrYx+nij42rvftON1+ztemAf82M57/OVXA+jdSuH9GHTMSBz5wjOV7sDoKjI/aHVBykbWxiAmRmbuYsnxa8+KBM7x8qJF5+Q4+f0lZQBM7QEbXQaJ8AzRmrAroeWcn9fPA7GribT4eOFHMvq4QLH0UrYauOuiutzoSMQGS3MWU5NOaP+ypZ3tlGxfOyOS2VcXERk++xFFW89wHPuZWbkQDic66yQc9Ql6cm/+cU01erIvvVhTybntywMeYlILlRv29ZofVkYgJkOQuphyvT/PbXbXsrungsjnZXLsgF6XMuc/f0VeFMy4Prz3elPOnRnv5xuwaZif18+PT+bzVNvY6+qCJSYTcRVBfDl6P1dGIcZKau5hStNb8YW8dB+u7uGZ+LhfPyjJtLOXzkNRfR1P6StPGAEiw+7h/Zi3fqyjkZ1V5eLRi9cRWXfqvfKN/xxWtgYa90HQA8peZG5MIKJm5iynljWPN7Knp5PK52aYmdoCk/jps2kt3Yomp4wDE2jRfnVHHYkcfj1Tn8VRlnOlj+iVzJsSnGRdWxZQiyV1MGeVV7bx5rJkV09K4bHa26eM5+qrRQE9CseljAcTYNF8pq2dZSi//uieZP1QH9uLwhCgbFK6ElhPQ3zn28SJkSHIXU0JFcy9/3FfPzOwk1i8pMK3GPlyyswZnXC5ee/Bm0dE2zT9Nr2dNlpuvlDt4pX7sjpWmK1wJaKP2LqYMSe4i5HUPuPlteS2ZSbHctqoYu838xK58XpKctUGbtQ8XY9M8uq6LxWkevrg9hbcaLU7wiVmQVgp15aBD8I5aMSpJ7iKk+bTm2V21uDxeNgRouaM/EgcasGvPpFoOTGr8KM3GCzuZ4fDwhe0ODnVYvPahcCX0NkJXrbVxCL9JchchbfPxZipb+7hhcT45juCVR5L7qoHg1dtHkxKjefxCo//Mne+mUNtn4Y9r/hKjz3udbIs8VUhyFyHrdGsfbx5tZmlRKsuK04I6tsNZjTM2a9L9ZCYrJ97HxnWdDHgVd25Npctl0b590QlG+9+GvbLmfYqQ5C5Cksvj43d76khLjOGGJflBuYD6Pu0j2VlLT4I1JZmRZqV4eWRtFzV9dj73Xop17YILloGrF06/ZVEAYjzkJiYRkl4/2kR7n4vPXlhKbFRwOycmDDRh97noSbSuJDOy8ZgN+GxxPz+tyucftkbxqaLm877flMZjWXMhKg4O/R5mXB7484uAkpm7CDl1HU7erWhlVUk607OSgj5+stPoSmllvX00l2R0c212O5ua03nbijYF9mjIXQhHXwRPiO8mJSS5i9Di8vj4/Z56kuOiuGZBriUxOJw1DESn+rVfarDdXtjMvCQnj1TnUum04Can/GUw2AUVrwd/bDEuktxFSHnk7VM0dg+wfkkBcUFa9vgBWpPkrKE3oSj4Y/shSsE/Ta8nJdrLD04V4PQG+Uc4cxYkZMDB54M7rhg3Se4iZDR09vOTzRXMz3cwN8+aWXOsq50YTx/dIVaSGc4R7eUfS+tpc0XzaHVOcO8rstlh7g1w4lVw9wdxYDFektxFyPivl4+hNVy3MM+yGBxD9fbeEE7uALOSBvh4fivvdqTwdnuQfxHOWw/uPqh4I7jjinGR5C5Cwo7KNl7c38DnLikjLcG62+2TnLW47fH0x2ZaFoO/bsxtY15SH7+syaVhIDp4A5dcCPHpcORPwRtTjJskd2E5r0/zzRePkJ8Sx+cuKbM0FsfZensw19VPkE3BfaVniFaaH5/Oxxus8szeJ43a+9E/wY7/M3rDn/0QIUOSu7DcM7tqOHqmm3/98DziYyy4iDokbrCVOFd7yC2BPJ+MGA+fLW6k0hnPpqYgbqqdt8hYDtl6PHhjinGR5C4s1Tvo4QevnWBVaTrXLbRm6eNZWR17gdBb3z6WNWk9rEjp4dmGTBqDVZ7JnGXc0HRmf3DGE+MmyV1Y6hdvnaKtz8W/Xjc3uC0GRpHVsQefiqIvzroLuhOhFHymuAm70jxSkxuc1TO2KOOGpqaD4PMGYUAxXpLchWWaugf4v3cq+cjifBYXpVodDlkde+mNL0DbrCsNTVR6jIfbC5s53JPIm20pwRk0d7GxHLLtZHDGE+PiV3JXSl2jlDqulKpQSt1/nuNuUkpppdSKwIUowtUPXjuB16f56tWzrQ4Fu8dJWvcxekL05iV/XJbZxbwkJ7+py6Z9MAh/BWXNAnuMlGZC1JjJXSllBx4GrgXmARuUUvNGOS4Z+AdAdtIVYzre2MNzu2v55NoSitITrA6HzK6D2LR3ytXbh7Mp+HRxI/1eGz84HIRWxfYYyJ4HjQdBW9WqUpyLPzP3VUCF1rpSa+0CngHWj3LcfwLfBQYCGJ8IU995+ShJsVF88bIZVocCQFb7HjSK3oRCq0OZlKJ4F1dldfBUZTxHO4NQXspbbLQBbj9t/lhiXPxJ7gXA8L216oaee59SailQpLV+KYCxiTD13qk2Nh9v4QsfmkGqhTcsDZfVsZfO5FlB3QzbLB/Pb8URo3lgf7L5F1ez54EtWkozIcif5D5a8e79fzJKKRvwEPDPY55IqbuVUuVKqfKWlhb/oxRhQ2vNd145Rl5KHHdcUGJ1OAAon4fMzv20pC21OpSASIry8eV5fbzXEsOrDSb/8oyKhazZ0LhfSjMhxp/kXgcMv8pUCDQMe5wMLAC2KKWqgDXAC6NdVNVaP6K1XqG1XpGVlTXxqMWU9fKhRvbXdvJPV86ypuvjKFJ7jhPtdYZNcge4bXo/sxweHjyQjMvsnJu/FAa6oKPa5IHEePizE9MuYKZSqhSoB24Fbjv7ota6C3i/EYdSagvwFa11eWBDFVPFUztqRn3e69P87+snyE6OxeXxnfO4YMtu3wNAc9oy8lvesTiawNhd3c7Hsgf5r4oivrPLx9XZneN6/7h2csqZb6x7b9g7ziiFmcacuWutPcB9wKvAUeBZrfVhpdQDSqkbzA5QhI9dVe209bm4Zn4uthDq3XJ2fXt/vLV3yAbaYkcfc5Kc/KExg0Gfid/vqDij9n5mH/ikNBMq/FrnrrXepLWepbUu01o/OPTc17XWL4xy7KUyaxcjDbq9vHGsmZKMBGbnJlsdzl9pTVbHblrSllkdScApBbfmt9DhjubV5jRzB8tbAoPdULvd3HGE3+QOVREUb59soW/Qw7UL8ixvMzBcsrOaeFc7zenhl9wB5ib3s8TRy58aM8zdtSlnvrFq5vAfzBtDjIskd2G6rn43WytaWVSYEhI3LA2XNVRvD8eZ+1m35LfQ67XzZzO7RkbFGqWZw38Er8e8cYTfJLkL0/3lSCM+DVfPC72adnbHbgZi0ulOLLU6FNNMTxxkdWo3LzWl0esx8Ue+YBn0NUPlFvPGEH6T5C5M1dDZz96aTi4oyyAtMTRuWBouq32PsQQyhEpFZrgpv5UBn52Xm02cvWfPh7hUOPCMeWMIv0lyF6bRWrPp0Bniou1cOivb6nD+RvxAM8n9dWG1vv1ciuNdrEzt4eXmNPNq7/YoWPB3cPQlGOwxZwzhN0nuwjSHGrqpbOnjirnZlu6wdC7Z7bsAaEpfaXEkwfHR3Db6vHb+0mJie+VFt4KnH46+aN4Ywi+S3IUpBj1eNh08Q15KHKtKM6wOZ1Q57btwRSXT6bC+5XAwlCUOsNjRy0tN6eatey9aBWmlsF9KM1aT5C5MseV4C139bm5YnI/dFpr17Oy2cprTl6NV6P1VYZaP5rbR7YniDbNm70rB4g1w+m3oqDJnDOEXSe4i4Fp6Btl6spWlRalMywhCX/EJiB9owuGspjk9svaVmZvcz9wkJy82peM2a/a+9HYjye/+lTnnF36R5C4CSmvNSwcaiLIrrlkQeksfz8puN26ibkpfZXEkwXdjbhvt7mi2tjvMGSClAGZeDXufBK/bnDHEmCS5i4B6fncdJ5t7uXJeDslx0VaHc05/rbfPsjqUoFvs6KM4foCXmtLxmdXvfcWdxpr345tMGkCMRZK7CJjGrgEeeOkIJRkJrJkemhdRz8pp2xVx9fazlIIbctqpG4hlb5dJZbMZV4CjEMo3mnN+MSZJ7iIgtNZ87fcHcHt9fGxZYUh1fRwpvr+RZGdNxCyBHM3a9G4yot282GTSL2GbHZZ/Cio3Q8sJc8YQ5yXJXQTE7/bUs/l4C1+9eg4ZSbFWh3Nei078GIAYdxdlNc994CNSRCn4cE47R3sTONln0taCy+8wNtHe+Qtzzi/OS5K7mLTadif/8eJhVpWkh8zWeeeT0leJKyoRZ2yO1aFY6rLMLhLsXl5sNKklQVI2LLgJ9j0F/R3mjCHOyZ+dmIQ4pwG3ly/8Zg8K+J+PL8YWomva36d9pPRW0pVUFvb9ZMYSb/dxVVYHf2rMoHEgmtw4E1a2rPkc7H8K9jwBsefp47/izsCPHeFk5i4m5Vt/PsLB+i6+f/MSijNCq53vaNK6jf1Su5KmWx1KSLgmuwO70mwyq6FY3mKYtg52PgI+rzljiFFJchcT9qd99Ty5vYZ7Lp7OlfOmRokjt+09ALoSJbkDpEV7WZfezZa2FHrMage85vPQVQuNB8w5vxiVJHcxIYfqu/ja7w+ysiSNr1w9dXqz5LVuwxmbjTs6hLb6s9j1Oe0M+mz8pcWkrfhmX2f0m6ncDNqshfViJEnuYtzqO/v59OO7SEuI4eHblhFtnxr/jOzefrLa90hJZoTieBeLHb282pxmTksCmx3W3gudNdBxOvDnF6OaGj+VImR09bu5c+NO+t1eNt65kmyHScvoTJDTtgu7dktyH8X1Oe10eqLMa0mw5DaITpBdmoJIkrvw26DHy+ee2M3p1j5+8YnlzMqZWqWNguYtuO3xdCeUWB1KyFmY7Hy/JYEplZOYROPCauNB6Gs1YQAxkiR34ZdBj5fPP7mH9yrb+N5Ni7mgLNPqkMZHawqa3+JM5jq0TVYAj6SUMXuvG4hlf7dJLQlKLjQGqnrHnPOLD5DkLsbk8vi49zd7efNYM9/+6EJuXFpgdUjjltZ9hITBZuqzL7U6lJC1Lq2btGg3LzaZtCwyLgXyl0LtdnAPmDOGeJ8kd3Febq+PLz69h9ePNvGf6+dz2+piq0OakMLmt9AoGrIusjqUkBVlg+uyOzjUk8ihDpP+uim9BDyDULvDnPOL90lyF+c04PZyzxO7efVwE9/4yDw+sbbE6pAmrKB5Cy1pSxiMNWlWGiauyOok3ublFydMuiEttdhYFln1DmifOWMIQJK7OIfeQQ93bNzJ5uPNfOvGBdy5rtTqkCYsvr+R9O6j1GddYnUoIS/B7uOKrE7+XBtLbZ9J6aH0InC2Qssxc84vAOktI0bR6XTxqY27OFTfxUM3L5mSNfbhipteB6Au5zKLI5kars3u4OXmdH55MoFvLukd35v96d+euwhikqB6G2TPm1iQYkyS3MUHNHUPcMNPttLa62LDymKcLi9P7aixOqxJmXbmZdqT59CTNHX/+gimjBgPNxQP8NvT8fzj3D7SYgO8NtIWBcVroOIN6O+EeJM2645wUpYR76tu6+Omn2+jw+nmU2tLmJdg6RzBAAAVrUlEQVRv0g0tQZTorCez8wDVeddYHcqUcs8sJ/1exa9OxZszQPFa43PNe+acX0hyF4Zjjd3c9PP36B3w8NkLS5mRnWR1SAFR3PgqADV5V1scydQyO8XLFXmDbDyZQK/bhJYECRmQNQdqtku3SJNIchfsr+3kll9sx6bg2XvWUpgW+q17/TXtzCu0piyiL6HQ6lCmnPvm9tHltvFkpUmz92kXwGCXXFg1iV/JXSl1jVLquFKqQil1/yivf1kpdUQpdUAp9YZSalrgQxVm2FHZxt8/ugNHfBTP3XMBM6dYS4HzSe6rIr37qJRkJmhJuocLs108eiKeATMm19nzjAurtTtNOLkYM7krpezAw8C1wDxgg1Jq5CXuvcAKrfUi4Hngu4EOVATe2yda+NTGneQ4YnnungumxGYb4zG97k/4lF1KMpNw75w+Wgft/Pa0CbN3mx0KlkPTIXC2B/78Ec6fmfsqoEJrXam1dgHPAOuHH6C13qy1dg493A7I38AhbsvxZj7763JKMhL57T1ryU2ZOt0d/aF8HqbX/5GGrIvoj8u2Opwpa02WmxUZLn5xPAGXGfccFa0G7YWDz5tw8sjmT3IvAGqHPa4beu5cPgO8PNoLSqm7lVLlSqnylpYW/6MUAbX5WDN3/3o3M7KSePquNWQmxVodUsAtO/Id4gdb6YvNoazmuQ98CP8pBffOcdLQb+d3VSZMABz54CiEfb8J/LkjnD/JfbRL5aMufFVK3Q6sAL432uta60e01iu01iuysrL8j1IEzOZjzdzzxG5m5Sbx1F2rSUuMsTokU2R37sUVlURn8kyrQ5nyLs11sTjNzY+OJppTey9aBWf2QdMRE04eufxJ7nVA0bDHhUDDyIOUUlcA/wrcoLUeDEx4IpDePtHCPU/uZnZuMr/5zBpSE8IzsccPNJPac5LW1MWgZEHYZCkFX13Qy5l+O78xY+VM/lJQdjj4bODPHcH8+Ze/C5iplCpVSsUAtwIvDD9AKbUU+AVGYm8OfJhist471cZdvy6nLCuJJz6zipSEaKtDMs2M2udQaJpTl1odSthYl+NmXbaLnx5LDPy699hkKLvMqLv7pJlYoIzZfkBr7VFK3Qe8CtiBx7TWh5VSDwDlWusXMMowScBzSimAGq31DSbGLc5jZLuA6rY+Nr5bRWpCNB9dWsCmg40WRWY+m3eQmTXP0pE0UzpABthX5vfy0c3pbKyI54tznWO/YTwW3QK//6xxx2rJusCeO0L51VtGa70J2DTiua8P+/qKAMclAqS+o5/Ht1XhiI/iMxeWkhQb3u2ESs5sIs7Vzum8D1sdSthZmuHhyvxBHjmewCfK+kmNCWDPmTnXQXSiUZqR5B4QUpAMY03dA2zcdpqEGDufuXA6yXHhW4oBQGtmVz1BR/IsuhNLrI4mLH1lfi+9HsWPjwZ4K76YRJjzYTj8R2MzDzFpktzDVFvvII9tPU2UTfHpdaWkxId5Ygdy2raT1nOS4yW3G1cBRcDNTvFya+kAv6qIp6LbHtiTL7oZBjrh5F8Ce94IJck9DHU6Xfxy62m8WvPpdaVkhOE69tHMP/UoztgsqvKuszqUsPaVBb0kRGm+uS8ZHchuwNMvhYRMOCj3IgSCJPcw09wzwC+3nmbA4+XT60rJdoTXnafnktmxl9z2nRwtvROfPTJ+mVklI1bz5fl9bG2O4dWGAC6ntUfD/I/CiVdgoDtw541Q4X11LcJ09Ln4xKM76Rnw8Ol1JeSnmtTNLwSMvNN0dvVvcNsT8Cm73IUaBLdP7+fpyni+tT+ZS3PbiAtUhWbRzbDr/+DYn2HJhgCdNDLJzD1MdA+4+eRjOznd1scn1k6jOCPAF7xCWKKzntTeU5zJWIvPFv7XFkJBlA2+uaSHOqedhwN5cbVwJaROkxuaAkCSexjoHfTwqcd2cqyxm1/cvpyyrPDYaMNfRc2bcdvjaUpfYXUoEWVttpu/m9bPT48nsL89QEUApWDhx6FyC/Q0BeacEUqS+xTndHn49MZdHKzr4ie3LeNDcyKrA6Kjt5KUvkoaMi+SWrsFvrG4l+w4H1/e5Qhc35lFN4P2wSHpFDkZUnOfwvpdXj77q3LKq9v50YalXD0/1+qQgktriprfZDDaIbP2ANtx2v/+6ncW9vPtk8X887tRfLLI6D6yunQSdwdnzTb6zex/GtbeO/HzRDhJ7lOU0+Xh04/vYufpdn5w8xKuX5RvdUhBl959lKT+Bk7l34C2yT9lqyx2OLkyq4NNzWmsSO1hXnL/+E9SvvGDj9Omw+HfwZsPGm2BV9wZmGAjiJRlpqC+QQ93bDQS+0O3LOHGpedrrx+elM9NUdPrOGOzaU1dZHU4Ee/2gmZyYt38sLKAdlcAftEWLDM6RdbJFnwTJcl9iukecHPHxp3sru7gh7cuZf2SyEvsAHlt24lzd1Kdd7W09Q0BcXbNV8rqGPAp/udUweTr7zGJkDMf6neDz4wm8uFPfiqmkLbeQTY8sp19tZ386NalfGRx5JViAOL7G8lv2Uq7Yy7diaVWhyOGFMW7uK/0DKec8Xxtt2Pyd68WroTBHmg+GpD4Io0UKkPMyHa9Z3U6XTz2bhVd/S5uWzWNrn73OY8Nd8uPfheFpiZHmpGGmpWpvdyc38KzNVnMcHi4d84kWgNnz4PYFKjeGrgAI4gk9ymgqXuAx7dVMeD2cscFpZRmRs4NSiMVNL1JcdNfqM2+jMGYNKvDEaP4u9w2BuzJfO9QEjE2zV2zJnCBFcBmh2lrjXYE7ZWQPj2wgYY5KcuEuFMtvfzi7VP4fJq7Lpoe0Yk92t3DysMP0pE8izOZa60OR5yDUvD9ld18uHCABw8k88jxSbTBKF5rXFPZ9cvABRghJLmHsL01HTz+bhWOuGg+d2lZWPeK8ceyY98jbrCVHQv+A60C3G5WBFS0DX64qpvrCwf49sFkfnw0YWI1+LgUyF0Ee58E9wT/AohQktxDkNenefnQGZ7bXce0jATuubiMtDDdzNpfRWdeo6zuDxwp+wztqQusDkf4IcoG/7uqmxuLB/j+4STu2+GY2P6r09YZfd73PRX4IMOYJPcQ4xz08KttVbxzspXVpencsa6E+JjInqUm9Dey+tA3aU1ZyMEZn7c6HDEOUTZ4aGU3X1vYy8t1sdz4Ztr4N/nImGGsnHnnB+BxmRNoGJLkHkL21Xby8JYKTrf18XdLC1i/pIAoW2T/L7J5B7lw7z+jtIdti7+Dlq6PU45ScM9sJ09e3En7oI0Pv57Oj48mMOjv8nWl4NL7obsO9j1paqzhJLIzR4jw+TQ/3VLBTT/bhtZw90XTWVEyid4c4UJrVh1+gMyuA7y36Nv0JhZbHZGYhAuy3Wy6sp0r8gf5/uEkrvlLOpvPxPhXiy+7XGbv4yTJ3WJ1HU5u/+UOvvvKca5ekMsXL5tJUXqC1WGFhLmnH2d6/QscnPF56nJlTXs4yI338fCabn51YScauPPdVNa/mcYr9bH4zpfkz87eu2ph5yPBCndKk+RuEa9P8/i7p7nqobfZV9vJdz+2iJ9sWBrx9fWzymp/x9LjP6A692oOzvic1eGIALsk18WrV7bz7WXddLkUn3svhctfTednxxJo7j9HWiq7HGZeDZsfhI6qoMY7FUlyt8Ch+i4+/vNtfPPFI6wsSee1f7qYm1cWodQEVhKEoWkNm1h16D9oyFzHe4v/S3rHhKlYO9w2fYA3rm7nR6u7yIz18d+Hkli7KYM7t6bwu+o4ulzDfiaUgut/YDQUe/FLBHZ37vAjd6gGUXP3AN979TjP76kjLSGGh25ZzI1LCiSpDzOj5llWHv4WzWnLeGfZQ7JtXgSIssENRYPcUDRIZY+d56ri+FNNHJsbY7ErzaLkPkpO7WNunoOk2ChmzvgHVh55kPJnv8OJkr8/53lvWx3Z12gkuQdBR5+LR7dWsvHdKtxeH3ddNJ17PzSDlHhJXO9vZq01BS1bKGx5h46kmdTkXoXXHtk3bUWi6cle/r+FffzLgj6eOuxke0cyOzoc7N1bzx/31jMtI4EFeZeRmbGVZUe/S198HvU5l1kddkiS5B5gw5t59Qy42Xaqjfcq23B7fMwvSOHqeTlkJMXy5wNnLIwytNi9A5TV/YG03pM0py6hKv/DcgdqhLMpmJk4wMzEAW4vaGFryvUcOdPNkYZuXjrUxOvcwe/iz7B271fZtOB/cRZeaHXIIUeSuwkaOvvZdqqV/XVd+HyaBQUpXDYnmxxHnNWhhRxHbyXTG14k2t3D6bxraU5bYdRWGTarFxFNKchPjSc/NZ4r5ubQ1jvI4YZu7q////kf579x/YF7eeTQTdSnreWCjB4yYzzGG+3pEb2DkyT3AOnqd/Pi/gZ+/tYp6jr6ibHbWFmSxtrpmWQly8bNI8UNtrLoxI+ZUfd7+mMyOFp6B70JhVaHJQJkPHuwjldGUiwXz8qCWVm81f0Unl338QXXs2xtOcy3znwCnZDPBendlObZiKzt4j9Ikvsk9A162Hy8mU0Hz/DG0WYGPT5yHLFctzCP5cVpsqxxFLGDbcyqfpo5Vb/G5nPTkLGWuuxL5c5TMSFJjjR6Z6ynqr2AVc1b2GT/Gpu9K/lJ3YdZUzuDNUe3c93CPK6an0N2cmT95SzJfZzqOpxsOd7CluMtbK1oYcDtIzMplltWFnHT8kIO1nXJ6peRtI+sjj1Mr/sTJWc2Yfe5qMm5kn2zv0RO2w6roxNTnVI0ZaykNWUB+a1buaRjD5fH7qQ5pohnW67moT8u59//lMKKaWlcOS+Hy+ZkU5aVFPY/p0r7sVZUKXUN8EPADjyqtf7OiNdjgV8Dy4E24BatddX5zrlixQpdXl4+wbCDQ2tNXUc/5dXt7KhsZ+fpdipb+wAoTIvn8jnZXLswj5Ul6dhtxj+USN0daSS7x0l2x27mVm4krec4se5uvLYYWlMW0JixhoHYTKtDFGHK5h0ks+sApc5D0FWLVnZqUlbyu8FVPNExjw4cFKbFc9HMLNZMT2d1aQa5KVNnVq+U2q21XjHmcWMld6WUHTgBXAnUAbuADVrrI8OO+QKwSGv9OaXUrcBHtda3nO+8oZbcu/rdVLb0UtHcS0VLL0caujlY30Wn0w1AclwUq0rSWVuWwaWzsynLShz1N39EJnetSRhoIr37MJkd+8ns3EdG50Hs2oNXRdGVVEa7Yy4djjn4bJHdulgEz+rSdChaDYeeh4PPQ2c1WtloSV3COyzlufYyygcL8RBFQWo8CwtSWFDgYG6eg+lZSRSlxRNlD70b6PxN7v6UZVYBFVrryqETPwOsB44MO2Y98M2hr58HfqKUUtqfPwsmSWuNx6fxDn24PD7cXh+DHh8Dbi9Ol/HRM+CmZ8BDV7+b9j4Xrb2DtPYOUt85QF2Hk54Bz/vnjLHbmJGdxDXzc1lYmMKSolTm5Dren51HBK1R2oPd58LuHSDa4yTa00Ocq4M4VxsJ/Y0k9jfg6DuNo7eSOHcnAF4VRUfKPI6XfILGjDUk9VVJPV1YJ2ce5HwdLvt3OLMfdfRFsk++xscaf8nHFPiS4mlJnMVJijhQm82ho0ls0al0kIxTJZKWmkp6ioPslCSyUuJIT4ghLSGGlIRokmKjSIyNIjHGTmyUndhoGzF2G1F2RbTdRpRNYbcpy8o//iT3AqB22OM6YPW5jtFae5RSXUAG0BqIIId79J1KvvPyMXxan7/R0HnYbYr0xBgyk2LJT4ljZUkaBanxlGYmMiM7ieL0hJD8jW2mrPY9fGjX3SjtQ+HDpsfux9ofk05PYgl1uZfTkTybDsdc2h1z8dn/ujqorL/ezLCF8I9SkL/E+Lj836GnCWq2YavZQU7jAXKatnKhqxNG/mHpHPo4Ax5tw4uNHhJYMfjzcQ1vtykUYFMKFPzHDfPZsMrcO2j9Se6j/doZmVb9OQal1N3A3UMPe5VSx/0Y3xSVVg38V5mY8MsvuLqBKmBLIE8aBt8XU8j3ZXRjfF8+bcKQ7cD1kzrDbQ/CbRN/+zR/DvInudcBRcMeFwIN5zimTikVBaRgfAc+QGv9CCD9OgGlVLk/dbNII9+X0cn3ZXTyfTk3f2oPu4CZSqlSpVQMcCvwwohjXgA+NfT1TcCbwai3CyGEGN2YM/ehGvp9wKsYSyEf01ofVko9AJRrrV8Afgk8oZSqwJix32pm0EIIIc7Pr5uYtNabgE0jnvv6sK8HgI8HNrSwJ+Wp0cn3ZXTyfRmdfF/Owa+bmIQQQkwtkbXeTwghIoQkdwsopexKqb1KqZesjiWUKKWqlFIHlVL7lFKhc/uyxZRSqUqp55VSx5RSR5VSa62OyWpKqdlD/07OfnQrpb5kdVyhRBqHWeMfgaOAw+pAQtCHtNaynvuDfgi8orW+aWjFWoLVAVlNa30cWALvt0ipB/5gaVAhRmbuQaaUKgQ+DDxqdSwi9CmlHMDFGCvS0Fq7tNad1kYVci4HTmmtq60OJJRIcg++/wW+CvisDiQEaeA1pdTuobuZBUwHWoCNQ6W8R5VSiVYHFWJuBZ62OohQI8k9iJRS1wPNWuvdVscSotZprZcB1wL3KqUutjqgEBAFLAN+prVeCvQB91sbUugYKlPdAMiejCNIcg+udcANSqkq4BngMqXUk9aGFDq01g1Dn5sx6qerrI0oJNQBdVrrs7uaPI+R7IXhWmCP1rrJ6kBCjST3INJaf01rXai1LsH4U/JNrfXtFocVEpRSiUqp5LNfA1cBh6yNynpa60agVik1e+ipy/lgu+1ItwEpyYxKVsuIUJED/GGo93UU8JTW+hVrQwoZXwR+M1SCqATutDiekKCUSsDYROgeq2MJRXKHqhBChCEpywghRBiS5C6EEGFIkrsQQoQhSe5CCBGGJLkLIUQYkuQuxDBKqV6rYxAiECS5CyFEGJLkLsKaUuq/lVJfGPb4m0qpbyil3lBK7RnqH79+lPddOrzfvlLqJ0qpO4a+Xq6UemuowdmrSqm8oPzHCDEOktxFuHsGuGXY45uBjcBHh5qUfQj4vhq6NXYsSqlo4MfATVrr5cBjwIOBDVmIyZP2AyKsaa33KqWylVL5QBbQAZwBHhrqOukDCjDaHzT6ccrZwALgL0O/D+xD5xMipEhyF5HgeeAmIBdjJv/3GIl+udbaPdSlM27Eezx88C/bs68r4LDWOuK3uhOhTcoyIhI8g9GF8yaMRJ+C0VffrZT6EDBtlPdUA/OUUrFKqRSMbowAx4Gss/uYKqWilVLzTf8vEGKcZOYuwp7W+vBQO+F6rfUZpdRvgBeHNuHeBxwb5T21SqlngQPASWDv0PMupdRNwI+Gkn4Uxu5ah4P0nyOEX6QrpBBChCEpywghRBiS5C6EEGFIkrsQQoQhSe5CCBGGJLkLIUQYkuQuhBBhSJK7EEKEIUnuQggRhv4fpb4cx5NlvbkAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.distplot(unpriv[\"value\"], hist=True, rug=False)\n", "sns.distplot(priv[\"value\"], hist=True, rug=False)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "We create the final dataset combining the privileged and unprivileged groups" ] }, { "cell_type": "code", "execution_count": 78, "metadata": { "collapsed": true }, "outputs": [], "source": [ "groups = priv.append([unpriv]).reset_index(drop=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want a model to be able to find some sort of pattern so we decide that if `value` is less than 5.3 we assign the `output` value as 0 with probability 0.8" ] }, { "cell_type": "code", "execution_count": 79, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for i in np.arange(len(groups)):\n", " if (groups.at[i,\"value\"] < 5.3):\n", " groups.at[i,\"output\"] = np.random.choice([0, 1], p=[0.7, 0.3])\n", " elif (groups.at[i,\"value\"] > 6.2):\n", " groups.at[i,\"output\"] = np.random.choice([0, 1], p=[0.2, 0.8])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Disparate Impact\n", "\n", "We wish to evaluate if the unprivileged group is receiving similar treatment to the privileged group. Of course, as I purposely created this data to be unfair, we expect this to be untrue.\n", "\n", "Disparate Impact is a metric to evaluate fairness. It compares the proportion of individuals that receive a positive output for two groups: an unprivileged group and a privileged group.\n", "\n", "```Pr(Y=1|D=unprivileged) / Pr(Y=1|D=privileged)```\n", "\n", "We can visualise the proportion of favourable outcomes for each of the groups" ] }, { "cell_type": "code", "execution_count": 80, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAEN1JREFUeJzt3X+QXWddx/H3pymxAkXALFNs0iZoYIzAyLAGGFEqUg0DJAygpOrITyMzBhgRavBHdeKMIwFhRiczEqRQcNpQ6gCLRjMOP2SQX9mWCiY1uhMp3cZMlx9iW5Wy8PWPvXl6e3OTvU1zcjfd92tmZ8/znOec+91Mdj/3nHPPc1JVSJIEcN64C5AkLR2GgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNeePu4D7a9WqVbV27dpxlyFJ55Qbb7zxa1U1sdi4cy4U1q5dy/T09LjLkKRzSpJbRxnn6SNJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWrOuZvXJD34XXnllRw7doyLLrqIXbt2jbucZcVQkLTkHDt2jNtvv33cZSxLnj6SJDWGgiSpMRQkSY2hIElqDAVJUmMoSJKaTkMhyaYkh5PMJNkxZP07ktzc+/q3JP/VZT2SpFPr7D6FJCuA3cDlwCxwIMlUVR06PqaqfrNv/GuBp3RVj3Qu+OrOJ427hCVh/huPBs5n/hu3+m8CXHLVl8/aa3V5pLARmKmqI1V1D7AX2HKK8VcA13VYjyRpEV2GwsXAbX3t2V7fCZJcCqwDPn6S9duSTCeZnpubO+OFSpIWdBkKGdJXJxm7Fbihqr47bGVV7amqyaqanJiYOGMFSpLuq8tQmAXW9LVXA0dPMnYrnjqSpLHrMhQOAOuTrEuykoU//FODg5I8AXgU8NkOa5EkjaCzUKiqeWA7sB+4Bbi+qg4m2Zlkc9/QK4C9VXWyU0uSpLOk06mzq2ofsG+g76qB9h92WYMkaXQ+T0HSkrPqgu8B873vOpsMBUlLzhuf7OQG4+LcR5KkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpKbTUEiyKcnhJDNJdpxkzC8mOZTkYJJru6xHknRqnT2OM8kKYDdwOTALHEgyVVWH+sasB94M/GRVfTPJY7qqR5K0uC6PFDYCM1V1pKruAfYCWwbG/Bqwu6q+CVBVd3RYjyRpEV2GwsXAbX3t2V5fv8cDj0/yT0k+l2TTsB0l2ZZkOsn03NxcR+VKkroMhQzpq4H2+cB64DLgCuAvkzzyhI2q9lTVZFVNTkxMnPFCJUkLugyFWWBNX3s1cHTImI9U1Xeq6j+AwyyEhCRpDLoMhQPA+iTrkqwEtgJTA2M+DPwMQJJVLJxOOtJhTZKkU+gsFKpqHtgO7AduAa6vqoNJdibZ3Bu2H/h6kkPAJ4A3VdXXu6pJknRqnX0kFaCq9gH7Bvqu6lsu4A29L0nSmHlHsySpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNZ3evKal7corr+TYsWNcdNFF7Nq1a9zlSFoCDIVl7NixY9x+++3jLkPSEuLpI0lSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVLTaSgk2ZTkcJKZJDuGrH95krkkN/e+Xt1lPZKkU+tsmoskK4DdwOXALHAgyVRVHRoY+oGq2t5VHZKk0XV5pLARmKmqI1V1D7AX2NLh60mSHqAuQ+Fi4La+9myvb9CLk3wpyQ1J1gzbUZJtSaaTTM/NzXVRqySJbmdJzZC+Gmh/FLiuqr6d5DXANcCzT9ioag+wB2BycnJwH/fbU9/0vge6iweFC792JyuAr37tTv9NgBvf+qvjLkEauy6PFGaB/nf+q4Gj/QOq6utV9e1e813AUzusR5K0iC5D4QCwPsm6JCuBrcBU/4Akj+1rbgZu6bAeSdIiOjt9VFXzSbYD+4EVwNVVdTDJTmC6qqaA1yXZDMwD3wBe3lU9kqTFdfrktaraB+wb6Luqb/nNwJu7rEGSNDrvaJYkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqel0mgstbd9b+bD7fJckQ2EZu3v9z427BElLzEinj5K8fpQ+SdK5bdRrCi8b0vfyM1iHJGkJOOXpoyRXAL8ErEvS/4CcC4Gvd1mYJOnsW+yawmeA/wRWAX/a138n8KWuipIkjccpQ6GqbgVuBZ5xdsqRJI3TSJ8+SnInUL3mSuAhwN1V9YiuCpMknX0jhUJVXdjfTvJCYGMnFUmSxua07miuqg8Dz15sXJJNSQ4nmUmy4xTjXpKkkkyeTj2SpDNj1NNHL+prngdMcu/ppJNtswLYDVwOzAIHkkxV1aGBcRcCrwM+fz/qliR1YNQ7ml/QtzwPfAXYssg2G4GZqjoCkGRvb5tDA+P+CNgFvHHEWiRJHRn1msIrTmPfFwO39bVngaf1D0jyFGBNVf1NEkNBksZs1GkuHpfko0nmktyR5CNJHrfYZkP62imnJOcB7wB+a4TX35ZkOsn03NzcKCVLkk7DqBearwWuBx4L/BDwQeC6RbaZBdb0tVcDR/vaFwJPBD6Z5CvA04GpYRebq2pPVU1W1eTExMSIJUuS7q9RQyFV9f6qmu99/RWLXGgGDgDrk6xLshLYCrSpMqrqW1W1qqrWVtVa4HPA5qqaPo2fQ5J0BowaCp9IsiPJ2iSXJrkS+Nskj07y6GEbVNU8sB3YD9wCXF9VB5PsTLL5zJQvSTqTRv300Ut73399oP+VLBwxDL2+UFX7gH0DfVedZOxlI9YiSerIqKHwo1X1f/0dSS4Y7JMkndtGPX30mRH7JEnnsMWep3ARC/cbfH/vnoLjHzN9BPDQjmuTJJ1li50++nkWnrC2Gnh7X/+dwO90VJMkaUwWe57CNcA1SV5cVX99lmqSJI3JqBean5jkxwY7q2rnGa5HkjRGo4bCXX3LFwDPZ+HeA0nSg8ioE+L1P5+ZJG+j7+5kSdKDw2k9ZIeFTx4tNiGeJOkcM+pDdr7MvXMdnQc8hoXnIEiSHkRGvabwfOBRwE8BjwT2VdWNnVUlSRqLUU8fbQHeD6wCHgK8J8lrO6tKkjQWox4pvBp4elXdDZDkLcBngT/vqjBJ0tk38vMUgO/2tb/L8CerSZLOYaMeKbwH+HySD/XaLwTe3U1JkqRxGfU+hbcn+STwTBaOEF5RVV/ssjBJ0tk36pECVXUTcFOHtUiSxux0b16TJD0IGQqSpKbTUEiyKcnhJDNJdgxZ/5okX05yc5JPJ9nQZT2SpFPrLBSSrAB2A88FNgBXDPmjf21VPamqfhzYxX0f5CNJOsu6PFLYCMxU1ZGqugfYy8Kd0U1V/Xdf82HcO7+SJGkMRv700Wm4GLitrz0LPG1wUJLfAN4ArASe3WE9kqRFdHmkMOyO5xOOBKpqd1X9MPDbwO8N3VGyLcl0kum5ubkzXKYk6bguQ2EWWNPXXg0cPcX4vSzcKX2CqtpTVZNVNTkxMXEGS5Qk9esyFA4A65OsS7IS2MrA09qSrO9rPg/49w7rkSQtorNrClU1n2Q7sB9YAVxdVQeT7ASmq2oK2J7kOcB3gG8CL+uqHknS4rq80ExV7QP2DfRd1bf8+i5fX5J0/3hHsySpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJajoNhSSbkhxOMpNkx5D1b0hyKMmXknwsyaVd1iNJOrXOQiHJCmA38FxgA3BFkg0Dw74ITFbVk4EbgF1d1SNJWlyXRwobgZmqOlJV9wB7gS39A6rqE1X1P73m54DVHdYjSVpEl6FwMXBbX3u213cyrwL+rsN6JEmLOL/DfWdIXw0dmPwKMAk86yTrtwHbAC655JIzVZ8kaUCXRwqzwJq+9mrg6OCgJM8BfhfYXFXfHrajqtpTVZNVNTkxMdFJsZKkbkPhALA+ybokK4GtwFT/gCRPAd7JQiDc0WEtkqQRdBYKVTUPbAf2A7cA11fVwSQ7k2zuDXsr8HDgg0luTjJ1kt1Jks6CLq8pUFX7gH0DfVf1LT+ny9eXJN0/3tEsSWoMBUlSYyhIkhpDQZLUGAqSpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJKaTkMhyaYkh5PMJNkxZP1PJ7kpyXySl3RZiyRpcZ2FQpIVwG7gucAG4IokGwaGfRV4OXBtV3VIkkZ3fof73gjMVNURgCR7gS3AoeMDquorvXXf67AOSdKIujx9dDFwW197ttcnSVqiugyFDOmr09pRsi3JdJLpubm5B1iWJOlkugyFWWBNX3s1cPR0dlRVe6pqsqomJyYmzkhxkqQTdRkKB4D1SdYlWQlsBaY6fD1J0gPUWShU1TywHdgP3AJcX1UHk+xMshkgyU8kmQV+AXhnkoNd1SNJWlyXnz6iqvYB+wb6rupbPsDCaSVJ0hLgHc2SpMZQkCQ1hoIkqTEUJEmNoSBJagwFSVJjKEiSGkNBktQYCpKkxlCQJDWGgiSpMRQkSY2hIElqDAVJUmMoSJIaQ0GS1BgKkqTGUJAkNYaCJKnpNBSSbEpyOMlMkh1D1n9fkg/01n8+ydou65EknVpnoZBkBbAbeC6wAbgiyYaBYa8CvllVPwK8A3hLV/VIkhbX5ZHCRmCmqo5U1T3AXmDLwJgtwDW95RuAn02SDmuSJJ1Cl6FwMXBbX3u21zd0TFXNA98CfrDDmiRJp3B+h/se9o6/TmMMSbYB23rNu5IcfoC16V6rgK+Nu4ilIG972bhL0H35f/O4PzgjJ1AuHWVQl6EwC6zpa68Gjp5kzGyS84EfAL4xuKOq2gPs6ajOZS3JdFVNjrsOaZD/N8ejy9NHB4D1SdYlWQlsBaYGxkwBx9+evQT4eFWdcKQgSTo7OjtSqKr5JNuB/cAK4OqqOphkJzBdVVPAu4H3J5lh4Qhha1f1SJIWF9+YL29JtvVOz0lLiv83x8NQkCQ1TnMhSWoMhWVqsSlIpHFJcnWSO5L8y7hrWY4MhWVoxClIpHF5L7Bp3EUsV4bC8jTKFCTSWFTVpxhyv5LODkNheRplChJJy5ChsDyNNL2IpOXHUFieRpmCRNIyZCgsT6NMQSJpGTIUlqHeNOXHpyC5Bbi+qg6OtyppQZLrgM8CT0gym+RV465pOfGOZklS45GCJKkxFCRJjaEgSWoMBUlSYyhIkhpDQZLUGArSaUjS2aNspXEyFKQhkvx+kn9N8g9JrkvyxiSfTPLHSf4ReH2SS5N8LMmXet8v6W373iQv6dvXXb3vlyX5VJIPJTmU5C+S+DuoJcV3O9KAJJPAi4GnsPA7chNwY2/1I6vqWb1xHwXeV1XXJHkl8GfACxfZ/UYWnmFxK/D3wIuAG874DyGdJt+lSCd6JvCRqvrfqroT+Gjfug/0LT8DuLa3/P7edov5Qu85Ft8FrhtxG+msMRSkEw2bWvy4u0+x7vicMfP0freSBFg5ZMzJ2tJYGQrSiT4NvCDJBUkeDjzvJOM+w8IMswC/3NsO4CvAU3vLW4CH9G2zsTc77XnAS/u2kZYErylIA6rqQJIp4J9ZOPc/DXxryNDXAVcneRMwB7yi1/8u4CNJvgB8jPseXXwW+BPgScCngA918kNIp8lZUqUhkjy8qu5K8lAW/nhvq6qbHuA+LwPeWFXPPxM1Sl3wSEEabk+SDcAFwDUPNBCkc4VHCpKkxgvNkqTGUJAkNYaCJKkxFCRJjaEgSWoMBUlS8/8uJShVM+coHAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sns.barplot(x=groups[\"group\"], y=groups[\"output\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To calculate disparate impact we first calculate the proportions of the groups receiving favourable outcomes. We calculate a little function to do this" ] }, { "cell_type": "code", "execution_count": 81, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def calc_prop(data, group_col, group, output_col, output_val):\n", " new = data[data[group_col] == group]\n", " return len(new[new[output_col] == output_val])/len(new)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We calculate the proportion of the unprivileged group receiving the favourable outcome" ] }, { "cell_type": "code", "execution_count": 82, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.4975\n" ] } ], "source": [ "pr_unpriv = calc_prop(groups, \"group\", 0, \"output\", 1)\n", "print(pr_unpriv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we calculate the proportion receiving the unfavourable outcome" ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.709\n" ] } ], "source": [ "pr_priv = calc_prop(groups, \"group\", 1, \"output\", 1)\n", "print(pr_priv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, to calculate Disparate Impact, we divide the former by the latter" ] }, { "cell_type": "code", "execution_count": 84, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "0.7016925246826516" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pr_unpriv / pr_priv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The industry standard is a four-fifths rule: if the unprivileged group receives a positive outcome less than 80% of their proportion of the privilege group, this is a disparate impact violation. However, you may decide to increase this for your business.\n", "\n", "In this scenario, we are below the threshold of 0.8 so we deem this to be unfair." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Building an ML Model\n", "\n", "If we build an ML model on biased data, it's predictions will replicate the bias, let's demonstrate that here.\n", "\n", "We split our data into training and testing. The purpose of splitting the data is to be able to assess the quality of a predictive model when it is used on unseen data. When training, you will try to build a model that fits to the data as closely as possible, to be able to most accurately make a prediction. However, without a test set you run the risk of overfitting - the model works very well for the data it has seen but not for new data.\n", "\n", "We will only be randomly splitting our data into test and train, with a 80/20 split." ] }, { "cell_type": "code", "execution_count": 160, "metadata": { "collapsed": false }, "outputs": [], "source": [ "train, test = \\\n", " train_test_split(groups, stratify=groups[\"group\"], test_size = 0.2, random_state = seed)" ] }, { "cell_type": "code", "execution_count": 161, "metadata": { "collapsed": false }, "outputs": [], "source": [ "train.reset_index(drop=True, inplace=True)\n", "test.reset_index(drop=True, inplace=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We build a logistic regression model passing the feature `value` and label `output`" ] }, { "cell_type": "code", "execution_count": 162, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] }, { "data": { "text/plain": [ "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n", " intercept_scaling=1, max_iter=100, multi_class='warn',\n", " n_jobs=None, penalty='l2', random_state=123, solver='warn',\n", " tol=0.0001, verbose=0, warm_start=False)" ] }, "execution_count": 162, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr=LogisticRegression(random_state=seed)\n", "lr.fit(train[[\"value\"]], train[\"output\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can calculate the accruacy of our model on the test data" ] }, { "cell_type": "code", "execution_count": 163, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.7107142857142857" ] }, "execution_count": 163, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr.score(test[[\"value\"]], test[\"output\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can capture the predictions for our test data to evaluate fairness" ] }, { "cell_type": "code", "execution_count": 164, "metadata": { "collapsed": false }, "outputs": [], "source": [ "preds = lr.predict(test[[\"value\"]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We join the predictions with the group so we can filter accordingly" ] }, { "cell_type": "code", "execution_count": 165, "metadata": { "collapsed": false }, "outputs": [], "source": [ "pred_df = pd.DataFrame({\"group\": test[\"group\"], \"preds\": preds})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before we calculate the proportion of the unprivileged group receiving the favourable outcome" ] }, { "cell_type": "code", "execution_count": 168, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.75\n" ] } ], "source": [ "lr_pr_unpriv = calc_prop(pred_df,\"group\",0,\"preds\",1)\n", "print(lr_pr_unpriv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly for the privileged group" ] }, { "cell_type": "code", "execution_count": 169, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.995\n" ] } ], "source": [ "lr_pr_priv = calc_prop(pred_df,\"group\",1,\"preds\",1)\n", "print(lr_pr_priv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then divide the former by the latter to get our Disparate Impact value" ] }, { "cell_type": "code", "execution_count": 170, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.7537688442211056" ] }, "execution_count": 170, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr_pr_unpriv / lr_pr_priv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As with our initial data, we see a substantially greater proportion of the privileged group recieving the favourable output and identify this to be unfair.\n", "\n", "If this weren't a toy example and we were to put this model into production, this biased output could harm many people in the unprivileged group" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Disparate Impact Remover\n", "\n", "As we saw in the initial diagram, the distributions of `value` for the two groups are significantly different. If you were to pick a data point at random, you'd be able to predict with reasonable confidence which group you selected from. This means that even though you're not explicitly including `group` as a feature in your model, it is still present. Disparate Impact Remover removes the ability to distinguish between the two.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To use AIF360's Disparate Impact Remover, we need to create a `BinaryLabelDataset` for the training and testing data" ] }, { "cell_type": "code", "execution_count": 171, "metadata": { "collapsed": true }, "outputs": [], "source": [ "train_BLD = BinaryLabelDataset(favorable_label='1',\n", " unfavorable_label='0',\n", " df=train,\n", " label_names=['output'],\n", " protected_attribute_names=['group'],\n", " unprivileged_protected_attributes=['0'])\n", "test_BLD = BinaryLabelDataset(favorable_label='1',\n", " unfavorable_label='0',\n", " df=test,\n", " label_names=['output'],\n", " protected_attribute_names=['group'],\n", " unprivileged_protected_attributes=['0'])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "We then create a Disparate Impact Remover and specify the `repair_level`, this indicates how much you wish for the distributions of the groups to overlap. We will try with two values: 1.0 and 0.8\n", "\n", "We create our first DisparateImpactRemover and fit and transform the train and test data" ] }, { "cell_type": "code", "execution_count": 172, "metadata": { "collapsed": false }, "outputs": [], "source": [ "di = DisparateImpactRemover(repair_level=1.0)\n", "rp_train = di.fit_transform(train_BLD)\n", "rp_test = di.fit_transform(test_BLD)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly we create the second DisparateImpactRemover and again fit and trainsform the train and test data" ] }, { "cell_type": "code", "execution_count": 173, "metadata": { "collapsed": true }, "outputs": [], "source": [ "di_2 = DisparateImpactRemover(repair_level=0.8)\n", "rp2_train = di_2.fit_transform(train_BLD)\n", "rp2_test = di_2.fit_transform(test_BLD)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create new training and test datasets with the repaired values" ] }, { "cell_type": "code", "execution_count": 174, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "['group', 'value']" ] }, "execution_count": 174, "metadata": {}, "output_type": "execute_result" } ], "source": [ "train_BLD.feature_names" ] }, { "cell_type": "code", "execution_count": 175, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rp_train_pd = pd.DataFrame(np.hstack([rp_train.features,rp_train.labels]),columns=[\"group\",\"value\",\"labels\"])\n", "rp_test_pd = pd.DataFrame(np.hstack([rp_test.features,rp_test.labels]),columns=[\"group\",\"value\",\"labels\"])" ] }, { "cell_type": "code", "execution_count": 176, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rp2_train_pd = pd.DataFrame(np.hstack([rp2_train.features,rp2_train.labels]),columns=[\"group\",\"value\",\"labels\"])\n", "rp2_test_pd = pd.DataFrame(np.hstack([rp2_test.features,rp2_test.labels]),columns=[\"group\",\"value\",\"labels\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We remind ourselves of the distributions of the values for the two groups" ] }, { "cell_type": "code", "execution_count": 177, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 177, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4XNWd//H3mRmN6qj3YhVb7t3CxpheEgzBkAKYkJDAb8OmkL7JJrtZQkjZ9IVNCAkkkGoIhBBMi0MwYONubGzcZEu2mtV7GUnTzu+PK7PCyNZYmpk7M/q+nkfPaEZX93wtSx8dnXvuOUprjRBCiOhiMbsAIYQQgSfhLoQQUUjCXQghopCEuxBCRCEJdyGEiEIS7kIIEYUk3IUQIgpJuAshRBSScBdCiChkM6vhzMxMXVJSYlbzQggRkd544412rXXWeMeZFu4lJSXs3r3brOaFECIiKaVq/TlOhmWEECIKSbgLIUQUknAXQogoJOEuhBBRSMJdCCGikIS7EEJEIQl3IYSIQhLuQggRhSTchRAiCpl2h6oQIkB2P3rmj1XcHro6RFiRnrsQQkQhCXchhIhCEu5CCBGFJNyFECIKSbgLIUQUknAXQogoJOEuhBBRSMJdCCGi0LjhrpR6RCnVqpQ6cIaPK6XU/yqlqpRS+5VSSwNfphBCiHPhT8/9t8DVZ/n4aqB85O1O4MHJlyWEEGIyxg13rfUmoPMsh1wP/F4btgOpSqm8QBUohBDi3AVizL0AqB/1vGHktXdRSt2plNqtlNrd1tYWgKaFEEKMJRDhrsZ4TY91oNb6Ia11hda6IisrKwBNCyGEGEsgwr0BKBr1vBBoDMB5hRBCTFAgwn09cNvIrJnzgR6tdVMAziuEEGKCxl3PXSn1GHApkKmUagC+CcQAaK1/CbwAXANUAU5AFpAWQgiTjRvuWutbxvm4Bj4TsIqEEEJMmtyhKoQQUUjCXQghopCEuxBCRCHZIFuIYDvbBtYgm1iLoJCeuxDhYLgP2o+B1212JSJKSM9dCDP1NMB9C6G71nielAtLb4Pld0KS3MUtJk7CXQiz9NTD9gchPh2uvAcSMuDws7DpR7D9F7DyLiPkEzP+73O0ho4qqNkMddvBFgdDvZC/GOJSzPqXiDAk4S6EGYZ6jGC3xcLHn4P0UuP1pbcZwzMbvw2vfd8I+pILIaUQPENQvwt66oxjE7PB54bBLqjeCOd/Chy55v2bRFiRcBfCDNUvG2G96vP/F+ynZJbDTb+HloPw1l+g8kXoPA4WK+QthAs/D2WXQXqZcfzG78DOh2Dr/8KKT0Fq0bvbE1OOhLsQoTbUA7XboPA8SMo583E584y3K7959vOlFMIFn4NtP4d96+Dir4CSuRJTnYS7EKFW/TJoL8y4yr/jx5tKCZCYCXPWwN7fw8k9UFgxuRpFxJNf70KE0nD///XaEzMDe+78xZCcD0dfBJ8nsOcWEUfCXYhQanrTuAhaemngz60sMOtacHZA/Y7An19EFAl3IUKpca8xlz05SNsMZ8+FlCKo3Rqc84uIIWPuQoTKUI8x62Xme9/5uj9j6v5SCgqWwaG/QX9L4M4rIo703IUIlaY3AQ35S4LbTv4SQBkXVsWUJeEuRKg0vgmO/LNPfwyEuBTImG4MAekx96oXU4CEuxChMNgNXSeC32s/JX8JDLRC8/7QtCfCjoS7EKHQdsR4zJkfmvbyFhmzZw48FZr2RNiRcBciFNorITYldGu/2JMgYwYceyk07YmwI+EuRLBpH7QdhayZxmyWUMmcCa2HoE9mzUxFEu5CBFvPSXAPQOas0LZ7qr3jr4a2XREWJNyFCLb2SuMxK8ThnlJgrBV//JXQtivCgoS7EMHWVmlMgYx1hLZdZYGyS4yeu0yJnHIk3IUIJpcTuo6Hvtd+Stll0Ndk/IIRU4qEuxDBVL8dfF7j4qYZyi41HmXcfcqRcBcimGq3GsMjp++2FCppxcaOTSdeM6d9YRoJdyGCqXYrJBcaG1mbpfgCYzNtn8+8GkTISbgLESzuIWjYbazzYqZpK2GwE9qPmluHCCm/wl0pdbVSqlIpVaWU+toYH5+mlHpFKbVXKbVfKXVN4EsVIsI07gHvMKSHQbgD1G0ztw4RUuOGu1LKCjwArAbmArcopeaedtg3gCe01kuAtcAvAl2oEBGndovxaNZ4+ynpZZCYZQzNiCnDn577cqBKa31ca+0CHgeuP+0YDSSPvJ8CNAauRCEiVO1WY2cke6K5dSgF086XnvsU40+4FwD1o543jLw22j3AR5RSDcALwGcDUp0QkcrrgbodxsXMcDDtAuiuhV7pd00V/oT7WCsdnX672y3Ab7XWhcA1wB+UUu86t1LqTqXUbqXU7ra2tnOvVohI0bzPWE/m1Hi32aadbzxK733K8CfcG4CiUc8Lefewy/8DngDQWm8D4oDM00+ktX5Ia12hta7IysqaWMVCRIK6HcZjuPTccxdCTCLUSrhPFf6E+y6gXClVqpSyY1wwXX/aMXXAFQBKqTkY4S5dczF11W2D1GmQnG92JQarDQqXQf0OsysRITJuuGutPcBdwAbgMMasmINKqXuVUmtGDvsy8Aml1D7gMeDjWstKRWKK0toI0aLzza7knYpWQMtBGO43uxIRAjZ/DtJav4BxoXT0a3ePev8QsCqwpQkRobpOQH/L/41zh4vC5aC9xvz70ovNrkYEmdyhKkSgnRpvD7twrzAe63eaW4cICQl3IQKtbpuxX2rWHLMreaeEdGN1Sgn3KUHCXYhAq98BRcvBEoY/XkXLoWGnbN4xBYThd58QEczZCW1Hwm9I5pTC5TDYBR1VZlcigsyvC6pCCD+dGvIIl3Df/eg7n/c1G4+bf2LMnqm4PfQ1iZCQnrsQgVS3FSwxULDM7ErGlpQNMfHQecLsSkSQSbgLEUh12yF/iRGg4UhZILUYumrMrkQEmYS7EIHiHoSTe8JnSOZM0sugvxlcA2ZXIoJIxtyFmKxT49od1eBzg8v57rHucJI2sr689N6jmvTchQiUzuPGo9mbc4wnrdgYnumScfdoJuEuRKB0HoekXPM35xiP1W5s2i0XVaOahLsQgaB9Rk843Hvtp6SXQncdeFxmVyKCRMJdiEDoawLPkHGxMhKklxnXB5r2mV2JCBIJdyECoaPaeIyUcD91UbVeNs2OVhLuQgRCRzXEp0FChtmV+CcuGRIyjXn5IipJuAsxWVpDZzWkTze7knOTUQa1W8HnM7sSEQQS7kJMVn8LuPohY4bZlZybjBkw2Amth8yuRASBhLsQk3VqhcWMCOu5p4/8Mqp53dw6RFBIuAsxWZ3VxuYcCZlmV3JuEtKNdWZqNptdiQgCCXchJkNr42JqxnRQyuxqzl3pRVC7Rcbdo5CEuxCT0Xkchnsjb0jmlJKLjM07Wg+aXYkIMAl3ISbjxCbjMdIupp5ScqHxKOPuUUfCXYjJOPGaMd6emG12JROTUmjc0HRCxt2jjYS7EBPl8xmhmFkemePtp5RdYvwF4nWbXYkIIFnPXYjTrNtR59dxqb2VXONspzptFe0nOt/xsRWl6cEoLTimXwFv/BYadkHxBWZXIwJEeu5CTFBOxw4AehNLzC1kssouAWWFqn+aXYkIIAl3ISYot2MnvQnFuGJSzC5lcuJSoGiFhHuUkXAXYgKUz0N2525aMlaYXUpgzLjCWP63v83sSkSASLgLMQEZPQeI8Q7QnLHc7FICY8YVxmP1RnPrEAEj4S7EBOS2b0OjaI2WcM9dZCyfIEMzUcOv2TJKqauB+wEr8Gut9ffHOOYm4B5AA/u01h8OYJ1CmGZ63ZPveq208VkG4vIobI6SMLRYYMaVcGwDeD1glYl0kW7cnrtSygo8AKwG5gK3KKXmnnZMOfB1YJXWeh7whSDUKkRYsHqHSHI20JMUIbsu+Wv2NcZSBLI7U1TwZ1hmOVCltT6utXYBjwPXn3bMJ4AHtNZdAFrr1sCWKUToebw+upwujvbHcaAvgaP9cdQ4Y1E99Sg0PUkRup7MmUy/AqyxcOQFsysRAeDP314FQP2o5w3A6VMEZgIopbZgDN3co7X+e0AqFCIEPF4fe+u72VXTyfP7m2joGqRn8NQdmyXvOPbbtleYbY3jGw3nU+pwc0FaL4XxrpDXHHCxScac9yPPwXu/G9l33Qq/wn2s/2E9xnnKgUuBQmCzUmq+1rr7HSdS6k7gToBp06adc7FCBJLL4+OVylY2HGjmlcpWupxGmKcn2inOSCDbEYsjLoaZfduJtWhcPoXLp7i2eS9HVTmd3lj2NiXzVFMmJfFDXJLRwxVZ3cRaTv/xiCCzroFj/zB2Z8qZZ3Y1YhL8CfcGoGjU80KgcYxjtmut3cAJpVQlRtjvGn2Q1voh4CGAioqKCP4JEJHseFs/j++q5697Gmjvd5GaEMNls7K5Yk42q6Zn8uKB5nccP71u4O33Y12dpPs66M1dzg8yauh2W9nWlczmjmR+15DD+pZ0PpTXzrJisEXiXLRZq+G5LxhDMxLuEc2fcN8FlCulSoGTwFrg9JkwfwNuAX6rlMrEGKY5HshChZisQ429/PyVY7x4oBmrUlwxJ5u1503jovJMbFb/kjilvxrg7YupqTFeVmd3sTq7i8N98aw7mcXDdXm83OXhxxW9LM3wBO3fExSOXCg8D448C5d8xexqxCSMG+5aa49S6i5gA8Z4+iNa64NKqXuB3Vrr9SMfe49S6hDgBb6ite4IZuFCnMnpC391Drh44a0mDjX1EmuzcMnMLFaWZeCIi6GpZ4gndjf4fe7UviqGYtIYsr97YbA5jkHunVXH7p4kHmvK48ZX07hrtpPPzhmIrF78nOvgpbuhqwbSSsyuRkyQX5NZtdYvAC+c9trdo97XwJdG3oQICx6vj03H2ni1sg2LxeipX1CWSbzdOqHzKZ+H5IEa2tIWn/Fio1JwXmo/t83v5Jt7Hdx/OJFNLXZ+cX4PeQkRspXd3OuNcD+0HlZ9zuxqxATJnQoiKjV2D/L4rnra+4dZUJDCNQvySImPmdQ5Hc46rNrt1xTI5BjN/yzv5fK8Yb7+hoMbNqbxyIU9zEsNs2Ga3Y+O/XpKEex6WMI9gkm4i6jzRm0nz7zZSILdyscvKGFmjiMg503tr8KnrH4t8btjZH33bODu8k5+UFXIBzam8sWykyxJGTj7JxMG68HnLTKmRHbXQarMbItEkTQSKMRZuTw+vvbUfp7ac5JpGQncdXl5wIIdjHDvSyjGZ7Gf0+cVJwzzndm15MW6+EFVIZs6kgNWU9DkLTIeDz9rbh1iwqTnLsKOvzshjeb2+li3o47Klj4unZnFlXNzsATwJhy7q5v44XZaU5dO6PPT7R6+NauWH1UX8ouaPDRwSUZvwOoLuMQsSC6Ag3+DlZ8xuxoxAdJzFxHP5fHxh+21VLb0cf3ifN4zLzegwQ5Grx2g2zFjwueIs2q+OqOBeQ4nD9bk8Vq49+DzFkPDTujxfzaRCB8S7iKiuTw+frethurWfj64tIAVpRlBaSe1v4qhmFSG7JM7f6xF8+8zGpg/EvCvd4ZxwOctNh4PrTe3DjEhEu4iYvm05sk36qlpH+DGikKWFQfnIqTyeUjuP0G3ozwg663YLUYPfk6Sk1+cyOPNnsQAVBkESVmQuwAOPm12JWICJNxFxPrn4RYONvayekEei4vSgtZOsrMGq3bTnTTxIZnT2S2ar8w4SWH8MD89XsCxgbiAnTug5t4gQzMRSsJdRKS9dV28WtlGRXEaq6YHZyjmlNS+KnzK5tcUyHORYPXxH+X1pNo8/OBYISeHzm0WTkjMe7/xKEMzEUfCXUSchi4nf917ktLMRNYszkcFeWna1P4qehJL0JbJ3QQ15rljvPzHzHosCn5wrJBe98Tung2ajOmQI0MzkUjCXUSUYY+XP++qJynWxq3Lp2GzBPdb2DFQS5yrk56k8qC1kRvr5t+mN9DptvHj6gJcvjBbR33e9cbQTG+T2ZWIcyDhLiLKC2810zng4sZlhSTEBv82jfy2zcDkpkD6Y2bSEJ8paaJyIIFf1uaiw2lB7DlrjMcjz5lbhzgnEu4iYhxu6mVXTScXlWdSlpUUkjbzWzcxaM9k2B68C7anrEzvY21+G1s6U7j/cELQ2/Nb1izInAmHZdw9ksgdqiIi9A25eWpPA3kpcVw5Jyckbdo8TrI7d9Oaviwk7QHckNvBySE79x1KYXaKl6sLhkPW9phOLSyWWgzVL8PWn4F91C/WitvNqUuMS3ruIiI8t78Jl8fHTRVFfm+sMVk5HTtGpkAGb7z9dErBncXNLEpz86WdDo70hMkF1ryFoH3QfMDsSoSfJNxF2DvW0sdbJ3u4ZFYWOcmhmw+e37YZtzWRvoTQropot2geuqCHpBjNv2xJpXM4DC6wJhdCfDo07ze7EuEnCXcR1txeH8/sayQzyc4l5Vmha1hr8ts20ZS5Em0Jfe85J97HQxf00Dpk4a7tKXjM3udDKeNu1fZK8AyZXIzwh4S7CGuvHW2jc8DFmkUFIRuOAUjpP0biUAuNWReFrM3TLU738N2lfWxts/OjA2GwREHuQvB5ofWI2ZUIP0i4i7DV3jfMa0fbWFyUyozs0MyOOSW/1ZgC2ZS1KqTtnu7GkiE+Ot3Jr44m8nxDrKm1kFYCMYnQIuPukUDCXYSt595qxGZRrJ6fG/K289s20+mYzWBcaGbmnM1/Lepnabqbr+xycNTMC6wWK+TMhdZDRg9ehDUJdxGWjrb0cbSln8tnZ+OIC/xt/2cT4+4lq/tNGrPNG5IZzW6BB1f2kBij+eT2FPrdJl5gzVkAbid0HjevBuEXCXcRdrw+zQtvNZGeaGdlWXAXBRtLbvt2LNpr6nj76XLiffxsRS81fVb+/Q2HeXewZs0CSwy0vGVSAcJfchOTCDu7aztp7Rvm1hXTQnoR9ZT8ts24bA46UhaEvO1TTm2wPZoCbi5w81hDNtmqm6uzu8/4+UHbYNsWa9yt2nwA5r4/OG2IgJCeuwgrvUNu/nmohdLMRObmmbBLkfaR3/46TVmr0Jbw6/usyelkaUo/v2/IocqsNeBz5sNgJ/Q1m9O+8IuEuwgrD7xShdPl5doFeUFfyncsab1HiB9uD6shmdEsCj5T0kh6jIf7jufj9JrwI5wz13hsPRj6toXfJNxF2DjZPcijW2pYXJRKfmq8KTWcWgWyKfMCU9r3R5LNx+dKT9LhiuFhM1aQjEuBlCJokXAPZxLuImz85B+VAFw117zph/ltr9ORMo+h2EzTavDHzKQhPpTfztauZF7rSAl9AdlzoasGBtpD37bwi4S7CAsHG3t4eu9J7lhVSmqCOdvN2V09ZHTvpzHzQlPaP1fvz+1gbpKTR+pzaBwK7XRRcuYDGo69FNp2hd8k3EVY+P6LR0iJj+FTl04PedvT655ket2TLKr8KRZ8WHzut1+bXvdkyOvxl0XBXaWNxCjNz07k4wnl8ExKAcQmw9EXQ9ioOBcS7sJ0m462sflYO3ddNoOU+BD3QEdJ7avCbY2nPz7ftBrOVYbdwyeKmznujOeZ5hDeE6AskDMPqjaCxxW6doXf/Ap3pdTVSqlKpVSVUuprZznuQ0oprZSqCFyJIpp5fZr/fvEIhWnxfHRlsXmFaG1shJ003QiuCHJ+Wh8XpPXyVGMmNc4Qrj+TMw9cfVC7JXRtCr+N+12slLICDwCrgbnALUqpuWMc5wA+B+wIdJEiej299ySHm3r56tWzibWZt25K4lATMV4n3UnB3Ss1WO6Y1ozD5uWBmjzcodpgO3Mm2OLg6IbQtCfOiT9dlOVAldb6uNbaBTwOXD/Gcd8GfgjIYs/CL4MuLz/eUMmiolSuW5hnai2pfcfQYPTcI5DD5uPO4ibqBuN4qilEwzNWO5ReYoy7h9WO3gL8C/cCoH7U84aR196mlFoCFGmtZXt04bdHtpyguXeI/7xmjik3LI2W0l/FQHw+HlsYrJs+QctSB7gko5tnmjM43B2iv4JmvteYEtl+NDTtCb/5E+5j/dS9/WtaKWUB/gf48rgnUupOpdRupdTutrY2/6sUUae9f5gHX63mPXNzWB6sdVD8ZPM4SRo8GbFDMqPdVthKks3L1/ck4wtFZ3rm1cZjpcyaCTf+hHsDUDTqeSHQOOq5A5gPvKqUqgHOB9aPdVFVa/2Q1rpCa12RlRXCLdNE2Lnvn0cZcnv52urZZpdCSn81CqIi3JNsPm4rbOXNzhj+dDwEd/mmFBjb78m4e9jxJ9x3AeVKqVKllB1YC6w/9UGtdY/WOlNrXaK1LgG2A2u01ruDUrGIeIcae1m3o46PnF9MWVZod1gaS2p/FW5rAgMRNAXybC5M72VVtosfvpVI62AIZv7MXA3122GgI/htCb+N+z+vtfYAdwEbgMPAE1rrg0qpe5VSa4JdoIguWmvuWX+Q1AQ7X7xyptnlgPaR0l8dkVMgz0Qp+M6SPoZ9im/tC8Evz9nXgPbBMem9hxO/1jTVWr8AvHDaa3ef4dhLJ1+WiFbP7m9iZ00n//2BBaQkmHfD0ikZPQcjegrkmbS2t3FDro8nGrJYur+V+Q7nOX3+Oa0Hn7cYkgvgyPOw+MPnWKkIlujoqoiI4HR5+N7zh5lfkMxNFUXjf0II5LdtjugpkGdzXU4nWXYXv6vPxhvMi6tKwaxroOplcJ3bLxERPBLuImQeeKWK5t4h7rluHlaLuVMfT8lv20x/fCEeW4LZpQSc3aL5SGEbdYNxbGxPDW5js68BzyAcfzW47Qi/SbiLkDjc1MuvXjvOB5YUUFFi7tTHU+KG28noOUC3I7qGZEZbkdrH3CQnfz6ZSb8niD/uxRcaC4lVPh+8NsQ5kXAXQefx+vj3p/aTEh/DN973rpUrTJPXvhWA7qRykysJHqXgY0Ut9Hut/KUpiGvU2+xQ/h5jvrvXE7x2hN8k3EXQ/eb1E+xv6OFb188jPdGctdrHkt+6icHYTJxxuWaXElQlCcNckdnNhta04K77Puc6cHZA3dbgtSH8JuEugupE+wA/fekoV83N4doF5q4fM5ryechr32rslWry0gehcGN+OzEWzeMng3jzYPlVYIuHg38LXhvCbxLuImg8Xh9f/cs+7DYL37lhvunrx4yW2b0Pu6cvbDfCDrTUGC/X5XSyozuZYwNxwWnEnggz3wOHnwWfNzhtCL/5Nc9diHOxbkcdAC8damZXTRc3Livk5cOtJlf1Tvltm/ApG80Z5zOt6e9mlxMS78vp5KW2VNY1ZHP3zLrg/MEy9wY49AzUbYOSyNiuMFpJuIugONbSx6uVbSwrTmPJtDSzy3mXgtZNtKYvwx3jMLuUkIm3+vhgXjuP1OfyZm8iS1IGJn/S3Y++87lnGCwx8Or3Yf4HoeL2ybchJkSGZUTA9Q66eWJ3PVmOWK5bGH7rtSQ6G0jtr+Jk1sVmlxJyV2R2kxPr4k8NWcFZNdIWC9lzoGmfsSSBMI2Euwgol8fH47vqcXl9fHj5NOy28PsWK2jbBEBj9iUmVxJ6NguszW+jfiiOLZ3JwWkkfwkM90JHVXDOL/wSfj95ImL5fJqv/mUfNR0DfGBJIdnJQbpwN0n5rZvoTSyhL9HEPVtNdH5aH8XxQ/ylKRNPMHrvOfOMHvxJWRjWTBLuImB+9I9K/vZmI++Zm8OioiDf7j5BNo+TnI6dU3JI5hSLgpvz22ketrOpIyXwDVjtkLfIGJpxDwb+/MIvEu4iIP6wvZYHX63m1hXTuGRm+G7EktuxDat2c3IKDsmMtjSlnxkJgzzVlBmcDbULKoyLq7JDk2kk3MWkPbazjrufOcAVs7P51pp5YTWf/XQFLa/isjloS1tidimmUgpuLmij3RXDxvYg9N4zZkBcCux/IvDnFn6RqZBiUh55/QT3PneIS2dl8cCtS7FZw7O/ML3uSdA+pjX/g56kMsoa5C7KBQ4nc5Kc/LUpk0sze4i1BHAAXlkgfxlUvQT9bZAUvn/NRavw/EkUYU9rzQOvVHHvc4d477wcfvXRZcTFWM0u66ySnA3EeJ10OWaZXUpYUApuym+j22PjpbYgXCMpWg4+D+x/PPDnFuOScBfnzO318R9PH+BHGyq5fnE+D3x4KbG28A52gLS+SnzKQk+U7bo0GXMdgyxwDPBMcwZD3gAPpzlyoXA57PkD6GDuFiLGIuEuzkm308Vtv9nJYzvr+NSl0/mfmxaH7VDMO2hNWl8lvYmleK2xZlcTVm7Mb6fXY2NDWxDuJF76UWivhPqdgT+3OKsI+KkU4aKyuY8bHtjCG7Vd/PSmRfz71bOxhMmOSuOJc7UT7+qUIZkxzEoaZHFyP882pzPoDXAkzHs/xCTC3t8H9rxiXBLuwi/PvHmSGx7YwoDLy7pPrOADSwvNLumcpPVWAtDlmGlyJeHpxvx2+rw2/t4a4LH3WAfMfz8ceBqG+wJ7bnFWEu7irFweH9985gCff/xN5hck8/xnLwybbfLORXrvEfrj83HHBOmW+wg3I3GIpSn9PNuSgTPQvfdlt4N7QKZFhpiEuzij5p4h1j60jd9tq+VfLixl3SfOD9slBc4m0XmSpKFGOpPDZ4u/cHRjXhsDXivPtwR47L1gGeQuMFaQlAurISPz3MW7rNtRx/G2fh7bVY/b4+OW5dMoy0riyd0NZpc2IUUt/wSgM3mOyZWEt7LEYc5L7eP5lnS+4eok1R6gIFYKKu6A574IDbuMKZIi6KTnLt5Ba83W6nYe2XKChBgrn750OgsKgnAHYwgVNb/EQFwew/bwW1c+3NyU386Qz8JDlQmBPfGCG8HugN2PBPa84oyk5y7e5vb6uGf9QZ7b38ScvGRuWlZIbJjfmDSe+MFmsrr3UZ99udmlRIRp8cOsTOvj0SoHd5Q7yYybZO999GYeeQvhrSchZ76xJR/IZh5BJD13AUCP083HH93Jn3bUccnMLG5dMS3igx2gqOVlQIZkzsWN+e0Me+GXlYmBPXHxKuOO1fodgT2vGJOEu6C1b4ibH9rGzhOd/PjGRbx3Xi6WMF7861wUN71Il2MmQ7EZZpcSMfLjXHygeIjfV8fTPBjAiEjOh/QyqN0iuzSFgIT7FNfQ5eSmX26jtsPJIx8c+SZ6AAATsklEQVQ/jw8ti6z562eT6Kwnq3sfNXnXmF1KxPn83AG0hvsPBbr3fiE4O6DtSGDPK95Fwn0Kq27r50MPbqNzwMUf/2UFF5VH18p9xU1/B6A2f7XJlUSeokQft04f5ImaOKp6Azg8l7fQuLGpZkvgzinG5Fe4K6WuVkpVKqWqlFJfG+PjX1JKHVJK7VdKvayUmpr7l0WQug4nH354Ox6fjz//60qWFUfZTBKtKWl8nta0pTjjw2+T7khw1+wB4qyaHx8MYO/dYoNpK6H1kNGDF0Ez7mwZpZQVeAC4CmgAdiml1mutD406bC9QobV2KqU+BfwQuDkYBYvxrdtRd9aP9wy6eWhTNUNuH5+4qIy9dd3sresOUXWhkdp3lNT+anbO+y+zS4lYmXGaT8x0ct+hJPZ0OFma4QnMiaddAFX/hNqtgTmfGJM/PfflQJXW+rjW2gU8Dlw/+gCt9Staa+fI0+1A9AzcRpn+YQ+/ef0ETpeX21eVkJsSeXec+qOk8Xl8ykZ97lVmlxLR/qV8kMxYH99/KylwN5fGp0L2PKjfbmzFJ4LCn3AvAOpHPW8Yee1M/h8gGyeGIZfHx++21tAz6OJjK0soTAvwjSphQvk8lDY+S2PWhXLj0iQlxWg+N2eAne12NjbZA3fikgvBNQCHngncOcU7+BPuY82JG/N3uFLqI0AF8KMzfPxOpdRupdTutrY2/6sUk+bTmid219PYPcja86ZRkhngWRBhJK99C/HD7VQXvt/sUqLC2rJBypI8fHd/Eq5AzWDMLIfELNj16wCdUJzOn3BvAIpGPS8EGk8/SCl1JfCfwBqt9Zh/a2mtH9JaV2itK7KyomtmRrh78a0mDjX1cu3CPObkRffKiNMb/sqgPYPGrIvMLiUq2C3wX4v6Od5v47fH4gNzUmUxbmqq3wHNbwXmnOId/An3XUC5UqpUKWUH1gLrRx+glFoC/Aoj2FsDX6aYjO3HO9hS3cHK6RlcMD3T7HKCKm64nYLWTZwoWIO2xJhdTtS4LM/FZbnD/O/hRFqHAjSDuvA8sMXJejNBMu7/ktbaA9wFbAAOA09orQ8qpe5VSq0ZOexHQBLwpFLqTaXU+jOcToTYifYBntvfyKwcB9cuyDO7nKAraXwei/ZwvPAGs0uJOv+1qJ9hr+JHBwI0pGdPhPkfNNZ5l408As6vX8Fa6xe01jO11tO11t8dee1urfX6kfev1FrnaK0Xj7ytOfsZRSj0DLpZt7OOtAQ7N1UURc2SAmekfZTXPUFb6mJ6k8rMribqlDm83F7u5MmaePZ1BmjNwYo7wNUP+/8cmPOJt8kdqlHK4/Xxpx21uL0+PnJ+MfH2yF8EbDx57VtxOOs4WnyL2aVErc/OcZIV5+Xrexy4A3FxtWAZ5C6EXY/IRh4BJkv+Rqln9zfS0DXIrSumkROBuyedq+l1TzKzdh0uWxIxrl6m1z1pdklRyRGj+faSfj65LYWHjybw6dnO8T/pbN74LWTNgbf+DC99E9JL3/lxWRJ4wqTnHoX21Haxq6aLS2ZmMS8/sjfa8FfscCep/VW0pi1FW6L/rxQzXV0wzOqCIe47lEh1XwC+1gVLwRZrrBYpAkbCPcocae7lmX0nKctM5Mo5OWaXEzI5nbvQWGhNW2Z2KVPCt5b0E2fVfP0NB77JjqbYYqHgPGjaa4y/i4CQcI8i/cMePv3HPcTZrNx8XhFWS5RfQB1hd/WQ3b2HzpR5uGMcZpczJWTH+fjGon52ttv5Q3UA5r6XrAKfF+p3Tv5cApBwjxpaa7721H5qOga4eXkRjripM8d7Zu06rD43jZkXmF3KlHJj8RCX5g7z3f1JHOqe5OU7R97IRh5bZSOPAJFwjxJ/3F7Lc/ub+Lf3zqIsM8nsckLG6nEys3YdXUnlDMZNnWGocKAU/OS8XlLtPu7akcyAZ5J/KRZfCM52aKsMTIFTnMyWiQJvNfTw7ecOc9msLD558XQe31U//idFiRkNfyXO3S3ryATYjhOdfh/7r9MG+c7RIj692c6nS5oAWFGafu6N5i2EQ8lQsxmyZc/byZKee4TrGXTz6XVvkJlk56c3LcYyRcbZAazeQeYcf4SWtGX0JxSN/wkiKOY7nHwgr4PXOlJ4rWMS6xa9vZHHYRiQhQUnS8I9gmmt+epf9tHUPcTPPryUtMQALskaAWbWPk7CcBv7Z37W7FKmvA/mtTPPMcBDtbkc6pvEBdbiC4zxnprXA1fcFCXhHsEe3nycDQdb+Nrq2dG3Td44Ytx9zD3+GxqzLqQtXaY/ms2q4EtlJ8mOdfPj6sKJz3+PS4G8RcZqkbKRx6RIuEeo7cc7+MHfK7lmQS7/78LS8T8hysw+8Tti3T3sK5dee7hIsvn4+owGbEpz++uptA9NcIiw5GLwDBkBLyZMwj0CtfQOcde6vRRnJPDDDy1CRfuCYKdJGGxizonfUZt3NV0pc80uR4ySHevmqzMaaB2y8LHXU+kcnsD3ZnoppJXC8VfBG6B9W6cgCfcIM+zx8pk/7cHp8vCrjywjKXbqTXhacuQnAOyd9SWTKxFjmZE4xC9X9lDVa2Pta2m0Dk4gZqZfDoOdcOhvgS9wipBwjyBaa77x9AF213bxgw8upDxn6t2Nmd2xi+LmDRwquwNnfPSvTx+pLs118dsLu2kYsHDTa6mcdJ5j1OTMg8Rs2HK/rBY5QRLuEeQ3r5/gyTca+NzlM7huUb7Z5YScxeui4tD3GIjL43DZx80uR4xjZbabP1zcTcewhfdvTGNX+zncNa0sRu+9eT9U/TN4RUYxCfcI8UplK9974TCr5+fyhStnml2OKeYdf5jU/ip2zfsGXmuA9vIUQbUsw8OTl3aRYNPc8loqvzkW739HvLACUovh5XvBJ0sSnKupN2AbgQ6c7OGz6/YyOzeZn9y0aErdqHTK/KM/Z171w7SnLCB+qEXWa48gs1O8rL+iiy/vSubb+xzsao/h3sX9ZMePE9gWG1z2n/D0nXDoaWNLPuE3CfcwtG5H3dvvt/cN86tN1cRYLVy3KJ+/7W00sTJzWLwupp98Bq81jtrc95pdjpiA5BjNr1b28PDRBH5yMJEtLXa+umCAD5cNYj1bX2XBh2DLfbDxuzBnDVinzoJ4kyXDMmGsZ9DNI1tPAHDHqlJS4qfmN/bio/9D4lAzJ/Kvw2NLMLscMUEWBf86y8mGqzpZmO7hv/Y6uGFjGhub7GceqrFY4Yq7obMatv8ipPVGOgn3MNU35ObRLScYdHn5+AWlZDpizS7JFPmtrzG75o80py+nK3mW2eWIACh1ePnjRd3cv7yHrmELd2xJ5fqNafyj0Y53rJCfeTXMuhZe+R50VIe83kglwzJhqHPAxSNbTtA/5OG2lcUUpE3Ni4eOgRou2Pd1uhyzqMu50uxyRAApBddPG+aawmGero3j50cSuXNrKgUJXtaWDlJubSItxku11xiijM//EtdWb6LrT3fy8vLfGLNpzuLDK6aF4p8R1qTnHmaOtfTx0KZqBl1e7riwlLKsqbM2+2gx7h4u2X0XPmVj09L70Bbph0SjGAvcVDrEy+/t4MHzeyhN8vKTg0l8av8MvnuskD11XQy7vQzGZbN39pfJ6dzNvOqHzS47IshPTBjZfKyNzz62F63hExeVkZsSZ3ZJprB6h7h4zxdJHDzJxuW/YSCh0OySRJDFWGB14TCrC4c50Wfl/n2aLZ3J7H+jgb9ZFLNyHczPv5yMvF0sOvZz+hKLqcu72uyyw5qEexjw+TQ/21jFfS8fpTw7iesW5pORNDXH2C1eFxft+QLZnbvZtvB7tKUvNbskEWKlDi+3FHSyNr+dV5KuYV9DDwcbezjY2Muzlpt5IuEEK/b9J04SaM+72Oxyw5aEu8maegb52lNv8drRNj6wpIDvvH/+lJzuCGDzOFn15r+R376F7fO/RU3B+8wuSZhIKSjOSKQ4I5H3LcyjtsPJgZM9fKrxS/zS922u2HsXjxxYS3/WUpam9BNnHXU11poOFbebV3wYkHA3idvr49EtJ7jvn8fw+jTfuWE+t66YNuVWeDwlbqiNS9/4DKm9leyYdzfHiz5gdkkijFiUojQzkdLMRHwL83ih9VGGDv4Hnxh+jA0Nx/hm7UfIT4ljVXovi5P7zS43LEi4h5jXp/n7gWbuf/koR1v6uWJ2NvesmUdR+tSdv53XtoVVb/4bVt8wx6bdBCB3oIozsihFQU4WR7N/SeaeL3N52yYute7nz/2X8cuua/mFNY/3dblYk9rO+WUZWKfgHd0g4R4yPYNu1r95koc3n6Cu00lpZiIP31bBVXNzzC7NNLHDnSw89gDl9U/gjM3iSPGtDMZN3a+HODdaWWnOWkVnynyKWjfy0Z5/cKv1JfZaF/DIySu57dd20pISWD0/l9ULclleko7NOnUmCCpt0nKaFRUVevfu3aa0HSrt/cNsPtbG8/ub2HS0HZfXx+KiVD55SRlXzc09Y49i9PID0Sh2uIPyuieYc+J3WH1DHC2+hd74IrRlat6BKwLD7uohp3MnWd37ifEOMBybyea4S3mgcxl73dNIT4zlyjnZXD47mwvLsyJ2LwSl1Bta64pxj/Mn3JVSVwP3A1bg11rr75/28Vjg98AyoAO4WWtdc7ZzRlu4D7q8HG3p40hzL/saeth5opOqVmPsLy8ljmsX5PG+RfksKkwZd1w9GsPd6nGS376FouaXKGr+J1btpiH7Ut6c9UV6k8pkGEYEjNJelic0QX8rHN0APjd9jhm8FnsxD7Yv4uBQFjFWxbLiNFaWZbJyegaLilKItU1w39cQC1i4K6WswFHgKqAB2AXcorU+NOqYTwMLtdafVEqtBd6vtb75bOeNpHDXWtM/7KGtb5iW3mGaewdp7B6ioWuQus4BajucNHYP4hv5UsbaLBRnJFCamURZZiIFafFYptCFUuVz43DWk9ZbSVrvYbK69pLecwCr9jAUk0pt/mqOTVtLb1LZ258j4S4CaUXpyGwZ58huTvufgLptADhTZ7I3fiXP9s3ibx35DGk7dquF2XkOFhamMDcvhRnZSczITiI90W7yv+TdAhnuK4F7tNbvHXn+dQCt9X+POmbDyDHblFI2oBnI0mc5eTDCXWuNx6fxjrx5vBq3z4fHq3F5fLi8XobcPoY9XgZdPpwuD06Xl75hDwPDHnoH3fQOuekd9NDldBlvA27a+4cZ9rx7edK0hBimZSRSnJ5ASWYic/MczM5N5vWq9ugJc62xaA8Wnxurdwirb4gYTz92dz92dw+xri7ih9tIGGolcfAkiYMncTgbsGhj70uvstGZMg+3LZGepOn0JhSPe+u4EJP1driP1tMAh56ByhehditoL9pqpze5nBpbGQeHstnT66DWlUwnDvp1PLa4JDJSU8hJTSQ7JYGMRDvpiXbSEu044mwkxcaQGGslLsZKrM1CrM2K3WrBZlXYrAqrUlgtKqCz4PwNd38GnQqA+lHPG4AVZzpGa+1RSvUAGUC7f+X679ebj/PDDZWgQaPRGnxav91rngyrRZESH4MjzkZagp1sRxyzcpLJSLKTlRRLpsNOVlIcealx5KXEkWAf+8u3tbpj8sWYIL3nIFdu/xgKDdqH0j4s+LdJwlBMKs74PHqSymnIuZLexBK6kmfRmzQdnyVGeubCfCmFsPIzxttgN9RtR9VuIaVpH4tatrHI2c6HFXD6/YPdxpsXhVdbeM63ks+7P31OTStlzPJRGI/3rJkX9PVv/An3sX7lnB6l/hyDUupO4M6Rp/1KqUo/2o9WmQThl595eoE6YMdkTxRlX5eAkK/J2Mb5utwRpGZfGHmbuFu/B7dO/NOL/TnIn3BvAIpGPS8ETr+F8tQxDSPDMilA5+kn0lo/BDzkT2HRTim1258/raYa+bq8m3xNxiZfl7PzZ/BzF1CulCpVStmBtcD6045ZD3xs5P0PARvPNt4uhBAiuMbtuY+Mod8FbMCYCvmI1vqgUupeYLfWej3wG+APSqkqjB772mAWLYQQ4uz8msWvtX7XIJPW+u5R7w8BNwa2tKgnw1Njk6/Lu8nXZGzydTkL0+5QFUIIETwy4VgIIaKQhLsJlFJWpdRepdRzZtcSLpRSNUqpt5RSbyqlIuPW5RBQSqUqpf6ilDqilDo8clPhlKaUmjXyfXLqrVcp9QWz6wo3kblyTuT7PHAYSDa7kDBzmdZa5nO/0/3A37XWHxqZrTZ114YeobWuBBbD28ujnASeNrWoMCQ99xBTShUC1wK/NrsWEd6UUsnAxRiz0dBau7TW3eZWFXauAKq11rVmFxJuJNxD7z7gq+Dnff1Thwb+oZR6Y+ROZgFlQBvw6Mgw3q+VUolmFxVm1gKPmV1EOJJwDyGl1PuAVq31G2bXEoZWaa2XAquBzyilZOdjY9h0KfCg1noJMAB8zdySwsfIMNUaQBYuGoOEe2itAtYopWqAx4HLlVJ/NLek8KC1bhx5bMUYP11ubkVhoQFo0FqfWrDnLxhhLwyrgT1a6xazCwlHEu4hpLX+uta6UGtdgvHn5Eat9UdMLst0SqlEpZTj1PvAe4AD5lZlPq11M1CvlJo18tIVwKGzfMpUcwsyJHNGMltGhIMc4OmRNa9twDqt9d/NLSlsfBb408gQxHHg9nGOnxKUUgkYGwj9q9m1hCu5Q1UIIaKQDMsIIUQUknAXQogoJOEuhBBRSMJdCCGikIS7EEJEIQl3IYSIQhLuQggRhSTchRAiCv1/14O26u7tTAsAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "target_0_orig = train.loc[train['group'] == 0]\n", "target_1_orig = train.loc[train['group'] == 1]\n", "\n", "sns.distplot(target_0_orig[['value']], hist=True, rug=False)\n", "sns.distplot(target_1_orig[['value']], hist=True, rug=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We next look at the distribution of the repaired data with `repair_level` 1.0" ] }, { "cell_type": "code", "execution_count": 178, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 178, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd81fW9x/HX54zsHUIIJCHMMGQogThxKzigrXvUUVusrdZrW0eX13qH1tpbW0ctrlbrxqpoWQ5wsMOUlTASkhBC9h7n5Jzv/YNoIwZygJP8zjn5PB8PHuac8z0n70Ty5pvf+f6+PzHGoJRSKrTYrA6glFLK/7TclVIqBGm5K6VUCNJyV0qpEKTlrpRSIUjLXSmlQpBP5S4iM0QkX0R2ich93Tz+RxHZ2PmnQETq/B9VKaWUr6Snde4iYgcKgPOBUmAtcI0xZtthxt8BnGiM+Z6fsyqllPKRLzP3acAuY8weY4wLeA2YfYTx1wCv+iOcUkqpY+PwYcwQoKTL7VIgt7uBIjIUGAZ8fJjH5wBzAKKjo6eMGTPmqMIqpVR/t27duipjTEpP43wpd+nmvsMdy7kamGeM8XT3oDFmLjAXICcnx+Tl5fnw6ZVSSn1JRPb6Ms6XwzKlQEaX2+lA2WHGXo0eklFKKcv5Uu5rgVEiMkxEwjhY4PMPHSQi2UAisNK/EZVSSh2tHsvdGNMB3A4sBrYDbxhjtorIgyIyq8vQa4DXjG4zqZRSlvPlmDvGmAXAgkPuu/+Q2w/4L5ZSSqnjoWeoKqVUCNJyV0qpEKTlrpRSIUjLXSmlQpCWu1JKhSCfVsso1a/kvXD0z8m52f85lDoOOnNXSqkQpOWulFIhSMtdKaVCkJa7UkqFIC13pZQKQVruSikVgrTclVIqBOk6d6UOsbqw5qifs9tT/LXb1+Zm+iuOUsdEZ+5KKRWCtNyVUioEabkrpVQI0nJXSqkQpOWulFIhSMtdKaVCkJa7UkqFIF3nrtSRGENsy15iWvcR2V5FZHsVXrHhdsThcsZSFzOChujhVqdU6ht8KncRmQH8CbADzxpjHu5mzJXAA4ABNhljrvVjTqX6lNPdQErtRlLqNhLhrgPA5YilNTwZMYbotv0kNu4grXoVbWFJOL2t7Mq4gg5HlMXJlTqox3IXETvwJHA+UAqsFZH5xphtXcaMAn4BnGaMqRWRgb0VWKle1VxF5v7FpNbmYTMe6qOHUZRyDnsjsvHYIxEgyu4hwm4QbwdJDdtJrc3jpB2PMrL4TZZPfoTa+HFWfxVK+TRznwbsMsbsARCR14DZwLYuY34APGmMqQUwxlT4O6hSvcrdCiseh+V/YpCrmc0RU3lJLmFNSzql1eF4ka+GCoZhUW2Mi21hUlwcE7Im0BSTxamb7uOCldexMfsumPpLsOlbWso6vpT7EKCky+1SIPeQMaMBRGQ5Bw/dPGCMWXToC4nIHGAOQGam7r2hAoAxkL8Qs+g+pG4vG2Kmc2/TJRS0pRPv6GB4VBs58U2khLsxncOrXE62N0WxqCKR9w8kMyKqlTMmj6H2tHmcvOU/mbLj9+x8oYi1438DIj1G6I7uTaOOly/l3t3fTtPN64wCzgLSgc9E5ARjTN3XnmTMXGAuQE5OzqGvoVTfqimEhffCzsWUOjK5x/UrtrdMZlpyBd9NKmJUdNsRu9nlFZbXxPHW/mT+tqKIoclRVJz0e86KeYrxe57HY49g/Zi7j7nglToevpR7KZDR5XY6UNbNmFXGGDdQKCL5HCz7tX5JqZQ/dbTDij/j/eRR2r3Co+7rWOicxW2XjuFvUzPY+PZjPr1MmM1w9oB6zkiq5w1zDou3lvPEst1U5dyCY2g7Y4peosMeyebRd/TyF6TUN/lS7muBUSIyDNgHXA0cuhLmHeAa4G8iMoCDh2n2+DOoUn6xdwVm/p1IdQFLvLk8aruJa2eewse5mUQ47cf0kg4b5GYmMzIlhpdXF/O3lXspG3sTd6e3csLuuTRHprE743I/fyFKHVmP5W6M6RCR24HFHDye/rwxZquIPAjkGWPmdz52gYhsAzzA3caY6t4MrtRRaa2DD+6H9X+nwpbKva57CBtzIS9/6wRS4yL88imSY8L54ZkjeHtDKYu3V9I+ag4PDShnyraHqUqYRH3sKL98HqV84dM6d2PMAmDBIffd3+VjA/y0849SgaVoOfxzDqZxPy/JLJ7wXM79V0/l4glpiJ+Ph4c5bFyZk0G4086ynTX8YfRd/FfjbZy+4ecsOvVVPLoOXvURPUNVBZ+8F3wb5/VAwSLY9SENURnc4PotjUkTePWGHEakxPg10ojiN792+65EeKI+jXkF8UxKvZzr65/m7LW3UjhkFgC7M6/w6+dX6lBa7io0dbTDuuehMp9tg2ZzedFsTh2byUtXTSY2wtnrn94m8KOs/bR7bPzmwHTGpuwkp+4j6mJHUxs3ptc/v1J6loUKPa5mWPUUVBawdOANXFR0FRdMHs5fv5vTJ8X+JYfAncPLGBbVxvdqvkt9WCpD9y/C5nH1WQbVf+nMXQWcV1YXH/HxEcWHv4C1093ISWWvQksVi9N+yK2FZzB78mD+cOVk7La+X28eZjPcNXwf923P4tfu7/G4eYghlZ+wc9h1fZ5F9S86c1chQ7wdjC55HVprWDHsDm4tPINL0tv4wxWTLCn2L6WGu7lt6H7ea53AZ85TSateRXxjgWV5VP+g5a5CxtDyRcS0lrFv9A3csiOHnGQX/zetAYfd+r/m0xKbuHhgDXc03kibLZKpW/8bjNfqWCqE6WEZFRJSajeQWruevUmnc9PO6cQ4vTx1cgNh1vf6V65Nr2BLYxa/67iGB2qfY9qW31Idf0L3g+1J37wv5+beDahCipa7CnpRreVk7V9AXfQw7m66juImO/dnF1NY3kohsNtz5GP4hxrROzFxCPwwaz+/3n4W34taRPqBpdTEjsXYju3MWKWOJIDmNUodA2PIKl+IxxbBC1E3saYhnuvTKxgT02p1sm4Nj2rnkkF13N92HRHuWgbWrrM6kgpRWu4qqCU1bCO2pYT8AefybNlwxsW0MHNgrdWxjujytCoKnGNYZ8YwuPIzXRqpeoWWuwpa4nWTeeBDmiMG8XDDhbR7he8PLQ/4HXbDbIZbsw7wX65rCfM0k1a90upIKgRpuauglVa9knB3PR/HzuLz2kS+NaiGIRHBMQseE9NKQtIAFnqmklq1EkdHs9WRVIjRcldByeluZHDlcqpix/LfB04lLbyd2YOCayPSa4ZU8oTnMuzGrbN35Xe6WkYFpUHVq7CZDl60X8YBVxi/GVVMmK37i3sduqlXoEhwejhxcATvHTiFi6rz2J98Kh26a6TyE525q6Bj87oYWLuBytixvFA1mklxTZwQ12J1rGMyM6WG12wXY/e6GVi1yuo4KoRouaugM6BuMw5vG29xPo0dDq4cXGV1pGPmsMFZGQ7+5c1lYPVaHB3B+Y+UCjxa7iq4GC+DqlfTEDGYJ2qmcFJ8EyOj26xOdVwmxzezJHwGTuMiqXK11XFUiNByV0ElrWo5ka5q3rOdS7PHwZWDK62O5BdnZDhZ4J1GWs0a7B2BeQKWCi5a7iqoZBe9TLs9ht/VncPUhEaGRbVbHckvhka1szLmfCJoJ65yrdVxVAjQcldBI65pD4OrlrMsbDoN3nCuSAveY+3dOTU9jCWeKQypWQ3u4D7UpKyn5a6CRlbZ+3ix8fumC5kY28zQEJm1f2lAWAeb4s4hmlYO7FhhdRwV5LTcVXAwhqyyBeyKmcIudwozUw9/NaZgNiU9ms+8E4ksXgq654w6Dj6Vu4jMEJF8EdklIvd18/hNIlIpIhs7/3zf/1FVfzagbiMxrft403Uqg8JdTI4LzdP1ox1etiacRZxppGS7rpxRx67HchcRO/AkMBMYB1wjIuO6Gfq6MWZy559n/ZxT9XNZZf/CbQvnlYaJzBhYi4VXzet1Ywcnss6MIXrvR+BxWx1HBSlfZu7TgF3GmD3GGBfwGjC7d2Mp9W/idZO5fzGrnLl0OKI5M7ne6ki9KsJuqMyYQZKpY/e2PKvjqCDlS7kPAUq63C7tvO9Ql4nIZhGZJyIZfkmnFJBWtZIIdx0vNk3jpKGJRNlD/9qjZ58wlC8YSVzxBxhPh9VxVBDypdy7+wX40B2a3gOyjDETgQ+Bv3f7QiJzRCRPRPIqK0Pj5BPV+7LK3qfJFssyz0ROGZ5sdZw+Ee4QqjNnkmJq2L51g9VxVBDypdxLga4z8XSgrOsAY0y1MebLdWnPAFO6eyFjzFxjTI4xJiclJeVY8qp+xtHRQnrFUhaaU8gYEM+AmHCrI/WZ08YPZzvDSSxZgvF4rI6jgowvW/6uBUaJyDBgH3A1cG3XASKSZozZ33lzFrDdrylV6Mp74Rt3jSj+9zLHpPptODxtzGvPZWbqTkYU959rjjrtQs3QGYzd+xRfbN3IhFyrE6lg0uPM3RjTAdwOLOZgab9hjNkqIg+KyKzOYT8Rka0isgn4CXBTbwVW/UtiYz4NxLBVRpCb2Gh1nD6XO24UBQwlqWQJ3g5dOaN859M6d2PMAmPMaGPMCGPM/3Ted78xZn7nx78wxow3xkwyxpxtjNnRm6FV/yDGQ3zjTj7wnEhuUjPhh7kYRyhz2IXqoTMZYg6wddEzVsdRQUTPUFUBK7a5GKe3jcWeHM4K8eWPRzJ1XDYFZDFg/WN43XrWqvKNlrsKWImN+bTjZHdYNqOCfM/24+GwC1VDLyLNe4Bti562Oo4KElruKjAZQ1xDPp96JnBychsSwmek+mLauFFss40mZf2fMbpjpPKBlrsKSFFtB4jqqOcj7xSm9+NDMl9y2IXqqT8n1VSyY+FfrI6jgoCWuwpICY35eI2wPyqbBKeu8QY4+fzL2WQby8CNj2Nceq1VdWRa7iogRdXvZL0ZxYT+cUKqT5wOO5U5d5PsraZw0Z+tjqMCnC8nMSnVp8Lc9SS7ynjOezXTEpqsjmOJ1YXf3K9+t6eYjoQpLGcSJ6x/gjeSZtHhjOn2+dfmZvZ2RBXgdOauAk5s4y4ADkRnE9kPNgk7Gg67jeWZtxFPI6nbnrc6jgpgWu4q4NhrCynxpjA8OcrqKAEpdcwpfGCmkVP2D8JcdVbHUQFKy10FFk8HaW27+dxMZHJ8aF5t6Xg57TaWZ9xKhGlj6La/Wh1HBSgtdxVQ2ip3E0E7FdGjcejfzsPKyJ7CfHMGE/e/QUSbbp+tvkl/fFRAKS0qoN04SRwwyOooAS0yzM5nQ27BZjyMzNfZu/omLXcVUKJrtpLHGEbF6dr2nozOnsA871mMLfsnUa37e36C6le03FXAaKyrJs1bzr7I7JC+ALa/xEU6WZp6I8YYsgt0zxn1dVruKmDs3F0AQHjyUIuTBI9xY8byiudcRpe9Q0xzsdVxVADRclcBw165jWIGkRofbXWUoDEwNoIPkq/DZRyM36l7zqh/03JXAaG6xU22ewcVsSf0+x0gj9bY0aP4h+c8svYvJLqlxOo4KkBouauAsGlnERHiZmDmaKujBJ1hydG8F/VtPAhj93zzmrSqf9JyVwHBVb6DVsLIyBhudZSgIyKMGpnNvI7pDC99R9e9K0A3DlMBoLxFGN++mf0xYxju0L+Sx2Jiejwvbv0WV5lljCl6Ec6cAnnHMIvPudn/4ZQldOauLPf5njoybJXEDhlrdZSg5bTbGDxsHO95Tmbk3teh5Zu7Sqr+RctdWa5h3w4AUjKyLU4S3HKHJfFXz2zCvK2wZq7VcZTFfCp3EZkhIvkisktE7jvCuMtFxIhIjv8iqlC2t8nO6LbN1IQNhqgkq+MEtdgIJxFDJrDUeyLeNc+Cp8PqSMpCPR7gFBE78CRwPlAKrBWR+caYbYeMiwV+AqzujaAqNC3aa7jZtoP2QWdZHSXgjSh+s8cx34mK4G8dF3B2y++gfBMMmdIHyVQg8mXmPg3YZYzZY4xxAa8Bs7sZ91/AI4Beml35xBhDWcluwsRD7OAxVscJCaOi2ygKG0WZpELR51bHURbypdyHAF3PjCjtvO8rInIikGGMed+P2VSI21HeSHbbF7hsEZCkSyD9QQTOTWngOdf5UFsI9aVWR1IW8aXcuztf0Hz1oIgN+CPwsx5fSGSOiOSJSF5lpa7F7e/e3bCPs+0bYUA22OxWxwkZZyQ38D7TcRGms/d+zJdyLwUyutxOB8q63I4FTgCWiUgRcDIwv7s3VY0xc40xOcaYnJSUlGNPrYKe12vYsXE5aVJD2CBdAulPUXYvZ2XYecd7GmbfOnC1WB1JWcCXcl8LjBKRYSISBlwNzP/yQWNMvTFmgDEmyxiTBawCZhlj8nolsQoJ64trmdyyHINA6nir44Sca4e38oL7AsTrhlJd49Af9VjuxpgO4HZgMbAdeMMYs1VEHhSRWb0dUIWm+ZvKON++Hm9CFoTHWh0n5ExM7MAeP5htMhJTvBqM6flJKqT4tM7dGLPAGDPaGDPCGPM/nffdb4yZ383Ys3TWro6kvcPDmo2bGC9F2NMmWB0nJInA1cNaedk1HWkqh4Z9VkdSfUw38lB9bumOCqa5VoMTPSTTS1YX1pDmtfGYdxoP8CJV2z6jeNCFPT5vt+fgBT+uzc3s7Yiql+n2A6rPzVu3j4vDNmCSR0FMqtVxQlaMw8uoBPjEO4nkui1gvFZHUn1Iy131qcrGdtbmF5HDVmTMRVbHCXlnJtfzZsd0wjzNxDftsTqO6kNa7qpPvbtxH9PZgN14IFvLvbdNjGtmve0EGoliQP1mq+OoPqTlrvrUW+v3cUXMZogaAOlTrY4T8uwCJye3ML/jFBIb8rF52q2OpPqIlrvqM1vL6tmzv4qTveshe4aeldpHpifX85bnDOzGTVLjDqvjqD6i5a76zLx1pZzn2ExYRxOM/7bVcfqNzEgXtREZlDGA5PqtVsdRfUTLXfWJNreHtzfs45aEdQcPyQw7y+pI/cqZAxp4ryOXuKZC7B7duLU/0HJXfeK9TWV0tNQzuXXVwVm7XU+x6EunJDaw2DMVGx4SGndaHUf1AS131SdeWrWX7yZsPfiG3oQrrI7T7yQ4Pbij06gwCSQ2bLc6juoDWu6q120sqWNzaT3Xx6yB+EzImGZ1pH7p1OQmFnqmEt+0G5vXbXUc1cu03FWve2nlXtLDmhlcvQomXHZw4xPV56YlNPKBdyoO4ya+abfVcVQv03JXvaqm2cV7m8u4J3MHYjx6SMZCMQ4vrph06ky0HprpB/RdLeV3r6wu/urjTwsqcXV4mVL/IXUxI1lQFAtF/358RHGNFRH7rZOTm1lSksPshtWI14PRcw1Cls7cVa/xeA2rCqs5O7GKIY2bKByi2/9bbUp8Ex+ZHMJNG3HNhVbHUb1Iy131mk0lddS1uLkt6mM6bOHsTtcTl6wWYTe0xWbRbCJIaMy3Oo7qRVruqld4jWFZQQUj4zycWLeYvWkzcYUlWB1LAVOTW/nMO4GYhl16haYQpuWuesWWffVUNbm4c0AeTk8rBUOvsTqS6jQprpnPzCRiPPVEtldYHUf1Ei135XfGGJblV5IS7eTMuneoTJhEbfw4q2OpTmE2Q23sSADiGwosTqN6i5a78rsd5Y2UN7TxgyF7iWvZq7P2ADQ20cYm73DC6/UCHqFKy135lTGGpfkVJEY5uah1Pq1hyZQMusDqWOoQk+Ob+cQ7mYGuYhwdLVbHUb1Ay1351b++2E9pbSs3ZVSQXvkpBUOvxWtzWh1LHSLMZjgQnY0NQ1zjLqvjqF7gU7mLyAwRyReRXSJyXzeP/1BEvhCRjSLyuYjoAdZ+qM3t4aEFO0iLC+fqhudoDUtmR9b1VsdShzEwKYEKk4C9Tg/NhKIez1AVETvwJHA+UAqsFZH5xphtXYa9Yox5unP8LOD/gBm9kFcFsOc+L2RfXStzM5aQWrmewrSLyCr7l9Wx1GFMjm/hk9JJXNyyhirjwYierRpKfJm5TwN2GWP2GGNcwGvA7K4DjDENXW5GA7p4tp+paGjjqaW7uDCtldMbFtAWlkRl4olWx1JHEGE37I0YQxStRDeXWB1H+Zkv5T4E6Pp/vrTzvq8RkR+LyG7gEeAn/omngsWjS/Jxebz8d+pSotorKBl4ts4Eg0BEUjouY8dbU2R1FOVnvpR7d/uzfmNmbox50hgzArgX+HW3LyQyR0TyRCSvsrLy6JKqgLVqTzVv5JXy42nxpBS+S1PEYGri9G2XYHBCooc8k01is76pGmp8KfdSIKPL7XSg7AjjXwO+1d0Dxpi5xpgcY0xOSkqK7ylVwGpu7+CeeZvJTIzkjsbHwN3KniGzdM/2IBFp97LTOZZ0bxlOV0PPT1BBw5dyXwuMEpFhIhIGXA3M7zpAREZ1uXkxoBdp7Cd+t2gHJbUtvDhxM/ZdS2DsLFojBlodSx2FjoRhALhq9lqcRPlTj6tljDEdInI7sBiwA88bY7aKyINAnjFmPnC7iJwHuIFa4MbeDK0Cw4pdVby4ci/3nmTIWvcQjDwPss6Aolqro6mjkJUczf6qJCIbdElkKPHpYh3GmAXAgkPuu7/Lx3f6OZcKcPWtbu6et5mTklzcuv9BCIuB2U9B/oKen6wCSozTsNk2ntNceWz1uMAeZnUk5Qd6hqo6asYY7n5zE+0Nlbwc/hC2pnK4+hWITbU6mjpGDXEjiJFWHGVrrY6i/ETLXR21Zz7bw4pthSwc8BiRDUVwzauQmWt1LHUcEgek4TZ2okuWWR1F+YmWuzoqawpreH7RKt6P/z0Dmgrgyr/D8LOsjqWOU2yEk+0yglENq6yOovxEy135rKKhjcde/ifvht/PUG8JcuVLkD3T6ljKT8qjRjPaFNFaXWp1FOUHWu7KJ21uD3Of+wvPuH9FUqQduXkhjLnI6ljKj8KTMwFwFH5scRLlDz6tllH9m/F6+fiZ+/hl3TM0Jo4l+ntvQdxgq2MpP4uITaaCZDKrP7c6ivIDnbmrI3M1s+fpq7ioYi4FKRcQ/6OPtNhDlQj5sblM8WxkX7WerRrsdOauvuGV1cUAhLtqyV35Q7Kat/Nc1M2ET/kP1m+oOuJzRxTX9EVE1Uvq088mbvsCPluxmCGXXmF1HHUcdOauuhXVup+zV9xASvNO7nPei+30O7HZ9K9LqGsZcjpuHLi3L7I6ijpO+tOqviG2uYjzV36XsNYKvu/9JUNPvYJwh27f2x90OGPYFTGBMU2rqWhoszqOOg5a7urrGg9w9ppbcbvauLL9N4yaNoOkaD0dvT+pHHQGY2wlfJa30eoo6jhouat/a2+Ely/H0V7Dd9vuJmt8LiMHxlidSvWx+vSzAajdrPsEBTN9Q1Ud5HHDGzdiDmzlNtfPMGmTuSFsGVJsdTDV1xpjRlAfnkZmzXJqml36m1uQ0pm7Ouij38Luj3jIfiubIqbxnZPS9Xob/ZUIHcPP5VTZwoeb9V/3YKXlrmD3UljxOEtjL+WF1jO4emoGEU59A7U/S5p0MTHSxu51H1gdRR0jLff+rqUG3rmNuuhh3Fb5HX550VjSE6OsTqUsJsPPxC3hDCpfSnVTu9Vx1DHQY+79Sd4LX79tDKx7AdNUwS2u35KbCjeFLWNNkZ6I1O+FRdOWcQbnFa1n4Rf7uf6ULKsTqaOkM/f+rGw9lG/mH87LKLBl8fCURj3Orr4SM2kWGbZKNq9fbnUUdQy03PurjjbY9i5VEUP5z4ZZ/GZSE2lRXqtTqQAio2dgEFL3L6WiUU9oCjZa7v1VwRJob+DHTd9j+iA3VwzVH151iNhU2lJP5DzbOhZ+UW51GnWUtNz7o6YDmMJP+Mx5OlsZwUMn6eEY1b3IEy5hkm0PKzZstjqKOkpa7v2NMbD1bTrEyV2N13HX+GY9HKMOL/tiAJLLllFer7/dBROfyl1EZohIvojsEpH7unn8pyKyTUQ2i8hHIjLU/1GVX1Rsh8odPGUuIykumhtGtFqdSAWylGzccVmcb8vj/c1lVqdRR6HHchcRO/AkMBMYB1wjIuMOGbYByDHGTATmAY/4O6jyA2OgYAH1jgE83nohv53ciFN/d1NHIoJz/CWcbt/KgnW7rE6jjoIv69ynAbuMMXsAROQ1YDaw7csBxpilXcavAq73Z0h1dL682MahcrasZHR9Kf/rnsO0xBZszQdYXdjH4VTwyZ6Jc+UTpFV8Rn75KWQPirU6kfKBL/O2IUBJl9ulnfcdzi3AwuMJpXqB8ZJesYwySeU9cxrXp1dYnUgFi8xT8EYP5FLHKv65odTqNMpHvpR7d+soTLcDRa4HcoDfH+bxOSKSJyJ5lZWVvqdUxy2zfAlR7RU83H4FM1PrSQ7rsDqSChY2O7bx3+Ic+0Y+WL8bj7fbH38VYHwp91Igo8vtdOAb76yIyHnAr4BZxphuN6Mwxsw1xuQYY3JSUlKOJa86BmI8TNj5FEUM5hNbDpek6vYC6iiN/w5hxsUJzStYubva6jTKB76U+1pglIgME5Ew4GpgftcBInIi8FcOFrv+vh9gMvcvIb65kN+5ruDbabVE2XXpozpKGbmY2DS+FbZaD80EiR7L3RjTAdwOLAa2A28YY7aKyIMiMqtz2O+BGOBNEdkoIvMP83KqrxnD2MLn2StDWO+YxPkDaq1OpIKRzYaM/zbTZROfb9lDi0sP6wU6n3aFNMYsABYcct/9XT4+z8+5lJ8Mql5JUsMOHnb/gCsza3Do0kd1rMZ/B8eqpzi9YzWLt07j2yemW51IHYH+qIe4sbufp4IkVkadx6mJDVbHUcEsPQcTn84VkWt5dU1Jz+OVpbTcQ1hi/VbSalbzrPtCpo8bgk33j1HHQwQZ/22meTeRX1jMzgONVidSR6AX6whWh154o4sRxQdXw4womUeTieRz52n8xr2k+0WtSh2NCVdgX/E433au4pU1J/Cfl463OpE6DJ25h6hwVy1JDdt5yXMeMwe36Kxd+UfaJBg0gVuiPuetdaW0uT1WJ1KHoeUeogZWr8FrbHzgOIuchCar46hQcuINZLQXkN6+i/c377c6jToMLfcQZPe0k1y7kfe9uZw92KOzduVfE6/A2MP5QcxyXl691+o06jC03ENQcu1Gwk07C+3n6Kxd+V9kIjL2Ui7GUv+jAAAPxUlEQVQyn7GtuIJtZboKKxBpuYca4yWxKo8872jGpsXrrF31jpO+S3hHAxc71/HSqiKr06huaLmHGHNgGwmeat6Rc8lN1KVqqpdkTYeEofwofgVvrd+nF9AOQLoUMsTU5X9Gi0kmNnUYNtFDMurojCh+8+AH9qSeBw8cx8iChQz2lPH3FUXcfeGY3g2njoqWewgxDftJbMznH+ZKTh+gxa6O3erCnncOdTKGSfIhP4/7iJ9/PpgB0eGEO+0AXJub2dsRVQ/0sEwIqdixnDbjpC1lAg491q56mdsZS9GQS7jQ/SER7nrW7tVN6QKJlnuocLWQULGGxZzCyQN1xz7VN3Zk3YDT284dcZ+yfFeVXsgjgGi5h4j9BWsJx0V75pmE2fQHTPWN+tiRlKWczlXeBbS1NrOptM7qSKqTlnsoMF4cxZ+z3oxm5rgBVqdR/cz2YTcR01HLjdGrWbqjQmfvAULLPQSUFeWT4q2kbOBZxDr1B0v1rQNJ06iJG8sPHAuoaW5jbZFexjEQaLmHgIadyzlgEjl1oi5FUxYQYfuwm0lp38vN8ev5eEcFze36vo/VtNyDXHllFWNcW9gWP52kSP3fqayxN+1CamOz+Qmv0d7exvOfF1odqd/TNghyRVtW4DIOxk6canUU1Z+JjY3Zd5HQXsZPE5fz10/3UNPssjpVv6blHsQqm9qZ0LSCLyKnMighxuo4qp/bP+BUypNzudH9OjZXA3/+aKfVkfo1Lfcgtn7jBqKljbRxp1kdRSkQYePo/yCqo47HMpfz4soiNuvSSMtouQep6jZDdu0yihzDGTxYr0KvAkNNwgnsHXQhZ9e8zgnR9dz31hd0eLxWx+qXtNyD1MebC8mSciJGnmF1FKW+ZsOYnyJi47mkf7Btfz3P6ZurlvCp3EVkhojki8guEbmvm8eni8h6EekQkcv9H1N1VdfiYlD5Mupt8QwaPsHqOEp9TUvkYDjvAVIqlnN/+kb++GEBxdUtVsfqd3osdxGxA08CM4FxwDUiMu6QYcXATcAr/g6ovumdD5Zyhm0zrozTwaYbe6oAlHMLZJ7KTY1zGWSr5963NuuZq33Ml5n7NGCXMWaPMcYFvAbM7jrAGFNkjNkM6MG1XlbX4iJmw9O4cJKSfYrVcZTqns0Gsx7H5mnn5UGvs3JPla6e6WO+lPsQoKTL7dLO+46aiMwRkTwRyausrDyWl+j3XvpwLZeaT2kelAthuvxRBbABI+GcXzPkwMf8aehy/vzxTj7fWWV1qn7Dl3LvbmfwY/r9yhgz1xiTY4zJSUlJOZaX6NcqGtuw5T2LUzwkjpludRylenbK7TD2UmZVPM0Vibu587UNHGjQS/L1BV8O2JYCGV1upwNlvRNHHcncD7fwI1lC27ALiIoZaHUcpf59Wb5Ddb1M39DTkJK1PNT2MBe0/Q8/enk9L38/l4jOqzap3uHLzH0tMEpEholIGHA1ML93Y6lDldS04F7/MknSRNRZd1kdRynfOSIg5xbseHg79vfs2lvCf7y2Ud9g7WU9lrsxpgO4HVgMbAfeMMZsFZEHRWQWgIhMFZFS4ArgryKytTdD90d//mA7N9sW4Bp0ImSebHUcpY5OTApMuZk4VzkfpTzGiq27uf/dLRijBd9bfFpHZ4xZACw45L77u3y8loOHa5QfvLK6+Gu3y+pa8Wx6k6ywcj4dfB+la0oYUax7Zqsgk5INU77HgLznWRj3v1y4+lcMbNzGneN6WAOfc3Pf5AsxeoZqgDPGsGBTCT9xvkN1zGhKB55jdSSljl3qeJhyE4Pde/lXzP/w0jYPf9wajU7g/U/LPcBtLq1nYt1HZMl+to76EUh3i5eUCiKDJiBTbibTW8qSqF+zcEcN/7s5Rgvez7TcA5irw8viLfv4afg71MaOpjT1bKsjKeUfgyYgp95BotPN/Ij/ZOfuAn65PpYOPQ3Sb7TcA9iyggrOdH1Kpinji5G3gej/LhVCEjKR039KeGwyfwt7hIySd5nzeTQNbv3t1B+0LQLUgYY2Vu0s557Id6iNHUVpqh5rVyEoMgE57U7IPIUfOeZze90j/PAjD8VNWk3HS7+DAcjjNcxbV8r3HIsZ4iljY/ZdOmtXocseBhOvgpNuYKKzhL+6f8WzH29h2X6n1cmCmjZGAPqkoIL2unJ+4nibfSnT2Z+ie7arfmDwSTjOvBtnwmAetM2lfc0LPLXJoyc7HSMt9wCzZV89H++o4KH4t3GadtaPvdvqSEr1nahkIk7/Me7s2Zxr38hVxQ/w5yce1f1ojoFuBh5A2twefvbGJqaGF3Ne+wdsH3YjjdFZVsdSqm+JjfWOSUSOSCNl7/vcVfPfvP+HpWw64VcMy8zo+fnAtbmZvRwy8OnMPUAYY7h73mb2VNTyf1F/py0ska0j5lgdSynLtEYMpGT0jazMnMMMWcX3t1xL4ap3ae/wWB0tKGi5B4gnl+7ivU1lvDrqE4a0bCdv3K9wO2OtjqWUpYzYKRx/B4tOeRm3M55f1f6amA/vofSAXg+iJ1ruAWDRlnIeXVLAz0ZXMKX4eXalf4eStAusjqVUwGhIGM/n57zFykHXcZn5gCvyriE/7yNcetbTYekxd4utKazhrtc3ctoQO7fXPoIkj2D92HutjqXUcVld6P+N7bz2cApPvI/aynPI2fAL7q+4i+c+vIx9J97JsNQEv3++YKczdwt9WlDJDc+vJj3ewfOxc5HmSrjsWTocUVZHUypg1aVMY9nZ77JlwEzmmHlcmncDq9csp8XVYXW0gKIzd4ss3lrOHa9sYFRKJP9MfZ7w/I/gksdg8IlQUtzzCyjVTxzuak9tg3LYFhHD8LIFPFJ1B48vuYzqtClMH9AIuT/r45SBR2fufczrNTz9yW5+9PJ6Thgcw9sZrxOe/y6c/1+6b7VSR6kxYQz5o+dQFT2Sn9tf4/zyv/LcDgdflNZbHc1yWu59qLy+jeufW83DC3dw8Zh4Xk99ibAvXoEz74XTfmJ1PKWCUocjhtKsy9k5+FuMsZfxlOdBPvjLT/n5q2vYV9dqdTzL6GGZPuD2eHl9bQl/WJJPm9vLXy6IZMb2nyF7dsBZvzhY7kqpYydCTeJEGmJHkF62mJ82zqNwx+c8sO1GMnO/xa1nDmdgbITVKfuUWHUNw5ycHJOXl2fJ5+4rHq/h/c1l/N8HBeytbuH0zEgeG7meAWsehfAY+M5cGNG522PeC189rzdWGijVn+ROnoR7wT04a3fzsfdE/ui9hpzc05kzfThp8ZFWxzsuIrLOGJPT0zidufeCoqpm3lxXwlvr9lHe0EZOqo3ncjcyYufzyIoqGHUBzHocYgdZHVWp0DTqPJw/XgWrn+asZb/jHPc9/GtNLjevuowR46Zy46lZTM1KREL4ymY6c/cDV4eXjSV1LMuv4JOCSraW1TNMyrkh4QsucaxhQNMOxHgPXiB41AxIGnbE19OZu1LHJ/eKLqtlWmpg5RN4Vz2Nzd3McibxvOs8SpNPZ9ZJmcyaNJiMpOBZfuzrzF3L/Si1uT3sqmii4EAjW/Y1kL+3lNbyArK8JYy2lzEtaj/Zrm1Em+aD452J1MSNoTp+PC2Rgy1Or1T/8LVy/1JLDax5BpP3AtK0n0rbQN52TWWRZyreITmcmZ3KmdkpTEpPwG4L3Bm9X8tdRGYAfwLswLPGmIcPeTwceBGYAlQDVxljio70moFa7sYYGlo72FfXyr66Vsqq66krL6K1qghqi4lqKSVDKsiUCrLkAMnS8O/n2pzIgFFUuMJoihxCU1QGreEpelFrpfpYt+X+JY8b8hfAhn9gdi9FvG5qbIksd2ezxjuGHY6xRKePY3zmQCYMSWDkwBiGJkfhtAfG4kK/HXMXETvwJHA+UAqsFZH5xphtXYbdAtQaY0aKyNXA74Crji360THG4DUH37z0GoPb48XtMXR4vLS7O2h3tdPuctPW1kZbawvtbS20tjTR3tJIe0sj7uZ6PK210FqHva2GiPZqEk09qVLDJKnlXOqxyb//AfQ6bLRFDYKEoUSkngrJwyFpOKSMQZKGgd1J4Zt/6IsvXSl1LOxOGDcbxs1G2uqhYDFJBYu5uGg5lzatAsCzz8be0lR2e9NYbpJ5W5LxRKXiiEkiIjaZ6LhEoqJjiYmNIzYmhsjwCCKjIokMDyfcaSfcYSfMYcNpF+w2wWGzYRP69Bi/L2+oTgN2GWP2AIjIa8BsoGu5zwYe6Px4HvCEiIjphWM+z3y6h0cW78BrwGsMXT/DM85HOdO2iUgMDjn6DYW82GgJS6A9PBlP9BDsCdNoS04nMiULSciE+Axs8RlEOcL8+BUppSwTEQ8Tr4SJV2IzBur2Qmke9sp8hh7YRlrFTuxNuwhzN0A7B/9U9/yyHcbGv7wnc6f79q/dLwI2ER6cPZ7rcof2ypf01efqqX9F5HJghjHm+523vwvkGmNu7zJmS+eY0s7buzvHVB3yWnOALzcpzwby/fWFBKEBQFWPo/of/b50T78v3euP35ehxpiUngb5MnPv7veIQ/9F8GUMxpi5wFwfPmfIE5E8X46b9Tf6femefl+6p9+Xw/PlHYJSoOu1rdKBssONEREHEA/oej6llLKIL+W+FhglIsNEJAy4Gph/yJj5wI2dH18OfNwbx9uVUkr5psfDMsaYDhG5HVjMwaWQzxtjtorIg0CeMWY+8Bzwkojs4uCM/ereDB0i9PBU9/T70j39vnRPvy+HYdlJTEoppXpPYKzKV0op5Vda7kopFYK03C0gInYR2SAi71udJZCISJGIfCEiG0Uk8PamsIiIJIjIPBHZISLbReQUqzNZTUSyO/+efPmnQUT+w+pcgUS3/LXGncB2IM7qIAHo7ENPflP8CVhkjLm8c8Va8Gxh2EuMMfnAZPhqi5R9wNuWhgowOnPvYyKSDlwMPGt1FhX4RCQOmM7BFWkYY1zGmDprUwWcc4Hdxpi9VgcJJFrufe8x4B7g6De/CX0GWCIi6zq3qlAwHKgEXug8lPesiERbHSrAXA28anWIQKPl3odE5BKgwhizzuosAeo0Y8xJwEzgxyIy3epAAcABnAT8xRhzItAM3GdtpMDReZhqFvCm1VkCjZZ73zoNmCUiRcBrwDki8g9rIwUOY0xZ538rOHj8dJq1iQJCKVBqjFndeXseB8teHTQTWG+MOWB1kECj5d6HjDG/MMakG2OyOPir5MfGmOstjhUQRCRaRGK//Bi4ANhibSrrGWPKgRIRye6861y+vt12f3cNekimW7paRgWKVODtzosZOIBXjDGLrI0UMO4AXu48BLEHuNniPAFBRKI4eBGhW63OEoh0+wGllApBelhGKaVCkJa7UkqFIC13pZQKQVruSikVgrTclVIqBGm5K6VUCNJyV0qpEPT/6ulJ9WL5QSAAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "target_0_rep1 = rp_train_pd.loc[rp_train_pd['group'] == 0]\n", "target_1_rep1 = rp_train_pd.loc[rp_train_pd['group'] == 1]\n", "\n", "sns.distplot(target_0_rep1[['value']], hist=True, rug=False)\n", "sns.distplot(target_1_rep1[['value']], hist=True, rug=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, you can no longer select a point without inferring which group it belongs to\n", "\n", "We next look at the distribution of the repaired data with `repair_level` 0.8" ] }, { "cell_type": "code", "execution_count": 179, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n", " return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 179, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4nNWV+PHvmRn13oslWbLcO0bYmN57SyixCUkgJEA2JEtCsgtJfoSQsgkJIexCsrAkEKoBA8EQgwkdDBj3IjfJli3J6r1LU+7vj5FBlmVpZE/T6HyeR49n3rnzvgchHd057y1ijEEppVRosQQ6AKWUUt6nyV0ppUKQJnellApBmtyVUioEaXJXSqkQpMldKaVCkCZ3pZQKQZrclVIqBGlyV0qpEGTzpJGIXAA8AFiBR40xvx30eh7wdyCxv80dxpiVw50zNTXV5OfnH03MSik1bq1fv77BGJM2UrsRk7uIWIGHgHOBSmCtiKwwxmwf0OxnwPPGmL+IyExgJZA/3Hnz8/NZt27dSJdXSik1gIjs96SdJ2WZhUCpMWavMaYPWAZcPqiNAeL7HycAVZ4GqpRSyvs8KctMACoGPK8EFg1qczfwpoh8D4gBzvFKdEoppY6KJz13GeLY4KUklwKPG2NygIuAJ0XksHOLyE0isk5E1tXX148+WqWUUh7xJLlXArkDnudweNnlRuB5AGPMJ0AkkDr4RMaYR4wxRcaYorS0Ee8HKKWUOkqeJPe1wBQRKRCRcGAJsGJQm3LgbAARmYE7uWvXXCmlAmTE5G6McQC3AquAHbhHxRSLyD0icll/s9uBb4vIZuBZ4Hqju4AopVTAeDTOvX/M+spBx+4a8Hg7cLJ3Q1NKKXW0dIaqUkqFIE3uSikVgjwqyyilvGTdY561K7rBt3GokKc9d6WUCkGa3JVSKgRpWUap4XQ1QfHL0FgKTWXuY/FZkJALk86ArPlg0T6SCj6a3JUaSusBWP0AbHwS7F1gi4LkAhArVH4GXY3w9i8gJh3mXgOLb3UnfaWChCZ3pQbb/wksWwq97TDnGlj8XciYBTJgmaXORih9C3a+Bp/+BT57BI77Gpx9F0QlBi52pfppcldqoG0vwsu3QGIefOttSCkcul1MCsz7ivurqczdy1//OJS8CV96GPJ1Tp8KLE3uanwaakhizRZY9zdIngTH3wBlH7i/PJE1D076Hmx8Ch6/GGZeBpPO9G7MSo2C3glSCqCzATY9475Ruug7EB4z+nMk5cNpP4LMubD9FSh73+thKuUpTe5KOR2w4e/ux8dfD9awoz+XLRIWfN2d4Itfhv2rvRKiUqOlyV2pXa9BawXM/ypEpxz7+SxWd4JPnwlbl0PjnmM/p1KjpMldjW8dte66et5iyJzjvfNabLDgGxCdDJufAUeP986tlAc0uavxbcer7jLMtIu8f25bhPvTQFcTFP/D++dXahia3NX41VgKtdug8ByIiPPNNZInQeFZUPEp1G73zTWUGoImdzU+GZd7REtkIkw63bfXmnqheybrjlfc11XKDzS5q/Gpbrv7Juq0i8Aa7ttrWW3u63TUwoH1vr2WUv08Su4icoGI7BKRUhG5Y4jX7xeRTf1fu0WkxfuhKuVFe99z99onHO+f62XNhfgJsPsNcDn9c001ro2Y3EXECjwEXAjMBJaKyMyBbYwxPzDGzDfGzAf+B3jJF8Eq5RVVm9z19oLT3MMW/UEs7t57V6O7/q6Uj3nSc18IlBpj9hpj+oBlwOXDtF8KPOuN4JTyiU//DNYIyDvRv9dNn+mexVrypvbelc95ktwnABUDnlf2HzuMiEwECoB3jvD6TSKyTkTW1dfXjzZWpY5dW5V7cbC8RRAW7d9ri8Dkc6Gn1b2OjVI+5ElylyGOmSO0XQIsN8YM2S0xxjxijCkyxhSlpaV5GqNS3rP2UfeIlfzTAnP99BnuWbD7PgrM9dW44UlyrwRyBzzPAaqO0HYJWpJRwcrldC8ONuU8iEkNTAxigYknQ9MeaKsOTAxqXPAkua8FpohIgYiE407gKwY3EpFpQBLwiXdDVMpL9rwD7dXuWaOBlLvQvTzBfu29K98ZMbkbYxzArcAqYAfwvDGmWETuEZHLBjRdCiwzxhypZKNUYG16GqKSYeoFgY0jPBayj4PKdWDXNWeUb3i0WYcxZiWwctCxuwY9v9t7YSnlZd3NsPOfUPRNsPl40pIn8k+FyrVwYB3knxLoaFQI0hmqanzYuhycfTD/2kBH4paYB3HZ7gSvlA9oclfjw6ZnIGOOezu8YJFTBC37oUOHBSvv0+SuQl9DCVRtgPlLAx3JobIXAOIuzSjlZZrcVegrftn976wvBTaOwaISIXWKezExHYegvEyTuwp9xf+A3BMhPjvQkRxuwvHQ1QAt+wIdiQoxmtxVaGsogbpimHVFoCMZWuY8sIRBpS4FrLxLk7sKbdv7t7ebcdnw7QIlLBIyZ0PVRl1MTHmVJncV2opfgZyFkDDkWnfBIXsB2DuhsSTQkagQosldha7GPVC7NXhLMgelTXdvpl21KdCRqBCiyV2FroOjZGYOt/1AELCGucfg12zR0ozyGo+WH1BqzFj32BePNzwBCXlQ8q9RnWJNWdMxh7GoIHl0b8ie7x7v3rDLvamHUsdIe+4qNPW2Q0s5ZMwKdCSeSZ0OtkgtzSiv0eSuQlNtMWDGTnK32iDzYGnGEehoVAjQ5K5CU20xRCZCfBCPkhksaz44eqB+V6AjUSFAk7sKPU67u3adMdO9b+lYkTYNwqKgWksz6thpclehp7HUvbxv+uxARzI6loOlma3g6A10NGqM0+SuQk9tMVjD3YtyjTVZx7lLM3veCXQkaozzKLmLyAUisktESkXkjiO0uUZEtotIsYg8490wlfKQMe61ZFKnusePjzWpUyEsGra9FOhI1Bg34jh3EbECDwHnApXAWhFZYYzZPqDNFOBO4GRjTLOIpPsqYKWG1VHj3lJv8nmBjuToWKyQORd2rQR7t7sGr9RR8GQS00Kg1BizF0BElgGXA9sHtPk28JAxphnAGFPn7UCV8kj9Tve/6dOP+VQuA/V9YVR2R9DqsCKARQwxVheTontIDvfRkMXs+VDxKZS+DTMu8c01VMjzJLlPACoGPK8EFg1qMxVARFYDVuBuY8wbXolQqdGo3wmx6RCVdFRvr+y08HpdEmtbYintjKLXdeTKZXKYnRlxXZyX1sK0mG7vDcxJmQJRye7lEzS5q6PkSXIf6kd28LYxNmAKcAaQA3woIrONMS2HnEjkJuAmgLy8vFEHq9Sw7N3QuBcmLh7d21zw+oEIHiuJZmOTu06fE9nLWakt5Eb1khvZS3K4A2PAIDTbbezpjKS0K5JNrbGsbkqgMLqbSzKaWJzUjsixL2GwaOZlsOV56O2AiNhjOpcanzxJ7pVA7oDnOUDVEG0+NcbYgTIR2YU72R+ytbsx5hHgEYCioiLdV0x51/6PwWV3r7LogR4nPFEaxd9Ko6nptlIQ6+COOR1kumrJirQf8X3pEXamxXb3n0P4oCmB12uTeKBsAu83dnDzxJpjL9nM/Qqsf9xde597zbGdS41LnoyWWQtMEZECEQkHlgArBrX5B3AmgIik4i7T7PVmoEqNaM877rHiKZOHbeYy8I/yCM5elcJvtsZRGOfksZNbePv8Jm6Z1jVsYh8s0mo4L62F+2aV8c3cGnZ0RPOj7QV82Bh/bP8tuSe6Fz3b8tyxnUeNWyP23I0xDhG5FViFu57+N2NMsYjcA6wzxqzof+08EdkOOIEfG2MafRm4UocpfRuSJ7nHuB+pSZuVH62LZ1NTGLMT7fy+qI2T0j1P5kdiETg/vYW58Z38eV8WD+7LpqI7gqUT6o+uFm+xwNyr4aP7oaPOfR9BqVHwaJy7MWalMWaqMabQGPPr/mN39Sd2jNsPjTEzjTFzjDHLfBm0UodpPQD1O45YknEZeLw0iovfSmZ/h5X7TmhjxdnNXknsA2VF2vnFtHLOTW3mldoUHt6fifNoC5BzrgHjgm0vejVGNT7oeu4qNByc0TlEcm/tE763Jp4PaiM4M7OX3xW1kx7p8lkoFoEb82pJCHOyvDqVdoeV2yZVEWYZZZZPnw5Z89ylmRO/45tgVcjS5QdUaCh7H2LSIS7rkMMVnRa+/G4Sn9aH86vj2vjbya0+TewHicDV2Q18M7eGda1x/HlfFq6j6cHP/Yp78+z63V6PUYU27bmroPPMmvLRvcEYvrT7XWqTT6B2X/PnuyBtbrJx4+pE7C548tQWFqV5twTjifPTW+hxWXjmQDop4Xauy6kf3QlmXwX/ugs2PgHn/co3QaqQpMldjQmF5S8c8bXI3nqiehtwifXzY+sbbVz3QRIpES6Wnd7C5PjA7U16WUYTjX1hvFqbQmq4nQvSW0Z+00FxGTDtQtj4NJz1/9wbaSvlAS3LqDEvvnMfAG0xBQCUtFn55keJZEQ5eems5oAmdnCXaK7PraUooZ3HKzLY1BozuhMUfRO6m2D74BHISh2ZJnc15iV07qM3LIHesEQa+mx8/cNEwq2GJ09t8Ut93RMWge9PqiIvqpcHy7Jo6BvFh+aCMyApH9Y/NlJLpT6nZRk1thlDXOc+WuKm0um08puSXFr64O5p5VTV9R42lTqQIiyG2yYd4M4d+TywN5ufTyvH5skYeIsFjr8e3rob6nZ6ZVE0Ffq0567GtOieWsKc3bRGF/B/5ZlU94Tz48kHyI8Ozp2MsiPt3Dyxht2d0Sw7kOb5G+dfB5Yw7b0rj2lyV2NafGcZAKv65vBJczzXZNczK64rwFEN76Tkds5La+bV2hQ2eFp/j02DmZfDpmegp9W3AaqQoMldjWnxnftoD0vhv6umMSuuk8szj201Rn/5ek4deVE9PLwvi3aHh7+GJ90KvW2wTnvvamSa3NXYZVzEde3nXftswsXFrfnVWLy1prqPhVkM382vpt1h5bHyDM/elH0cFJwOn/5FN9BWI9Lkrsas6J4abK4+/tU3h1vyvbDMrp/lR/dyZXYDq5sT+LQ5zrM3nXKbeyvBLc/7Njg15mlyV2OWta0SAHvcBIoSOwIczdG5PLORSdHdPFqeQYvdOvIbJp0JmXNg9QPgCo5hnio4aXJXY1ZHUzX7TTqX5Pp/WQFvsQl8N7+aHqeFxys8KM+IwMm3QWMJ7HzV9wGqMUuTuxqTtrVGMtm5h7qIfNIixlY5ZrCcqD6+lNXIJ83xbPRk9MysL0HqVHj3N+AK7OxbFbw0uasxx2HgvUoXydJBQkpmoMPxisszGsmJ7OWv5Zn0OEe4K2yxwhl3ujcD37rcPwGqMUeTuxpz3mlIZKJ9DwBdsaGx0brNAt+eWEN9XxgvVKeO/IaZV0DGHHjvv8A5dstSync0uasxxe4S/lGdwtnhxfTZ4ugNSwp0SF4zPbabc1Kb+WdtMtsOjDBRyWKBs34KzWWw6Wn/BKjGFI+Su4hcICK7RKRURO4Y4vXrRaReRDb1f33L+6EqBe82JtBot7HQspO26DyOboPS4HXthHribU5++o9tuEba3WPqBTChCN77Hdi7/ROgGjNGTO4iYgUeAi4EZgJLRWTmEE2fM8bM7/961MtxKvV5r/306P3EONtoj5kY6JC8Lsbm4ms5dWyuaOG5dRXDNxaBc38B7VWw5mH/BKjGDE967guBUmPMXmNMH7AMuNy3YSl1OHevPYwbEjcA0B4dGvX2wU5JbmNhfjK/e2MnzZ19wzfOPwWmnAcf/RG6xsbSC8o/PEnuE4CBXYjK/mODXSkiW0RkuYjkeiU6pfod7LVPjelipqsUhzWS7ohRrKo4hojAPVfMor3Hwb2rdo38hnPuhp42+Oh+X4emxhBP1nMfqqg5uBj4KvCsMaZXRG4B/g6cddiJRG4CbgLIywvNXpfyjfcb42m0h3FzfjVxtRW0R+WGXL19oOmVL3J9YSx/+2w/X4ndxPzkI4zlL7oBMmbBvKXu0syimyEhx7/BqqDkSc+9EhjYE8+BQ/dAMMY0GmMOrmT0f8DxQ53IGPOIMabIGFOUlhaavS7lfS5j+GddMgXRPSyIaiCqr5GO6ND/cHjbzE5SI138fGMcI91b5cyfAAY++IM/QlNjgCfJfS0wRUQKRCQcWAIcspmjiGQNeHoZsMN7IarxrqS2naqeCC5ObyKu210hbB8HyT0uzHDHnE42N4fx4v7I4Rsn5sKCb8DGJ6GpzD8BqqA2YnI3xjiAW4FVuJP288aYYhG5R0Qu62/2fREpFpHNwPeB630VsBp/Vpc2khRmZ3FSG3FdFbjEQkdUdqDD8osv5fUwP9nO77bG0G4foQx16u1gscH79/onOBXUPBrnboxZaYyZaowpNMb8uv/YXcaYFf2P7zTGzDLGzDPGnGmM2enLoNX4UdPaQ2l9B+enNWOzQFx3BV2RWRhLWKBD8wuLwN3z22notfLgjujhG8dnwQnfgi3LoH63fwJUQUtnqKqgtnpPA2FW4dy0FsTlIKa7alyUZAaan+zg6vxu/lYSzd72EZYFPvk2sEXB+7/zT3AqaGlyV0Gro9fB5ooWFuQlEWtzEdNTjcU4x11yB/jx7E4irIZfb4kdvmFsmnsETfHL0LzfP8GpoKTJXQWtz8oacbgMJxW6F9KK6xo/N1MHS490ceuMLt6ujmB17QglqRO/4x4m+ulf/BOcCkqa3FVQcroMa/c1MyU9lrS4CABiuyroCU/CYRuh9xqirp/cRU60k19ticU53NDIhByYfRVseAK6m/0WnwounkxiUsrvdtW009pt59K5/aNsjSGuq4KW2MmBDcxP1pQNvZTAlRl2HiibwO/XOTkztRXK7huy3aKUPLB3wivfhcnnuks1alzRnrsKSp/tayQ+0sa0zHgAIvqaCXN2jYvJS8NZnNTOlJhulh1IG35Tj/hsSJsOZR+Cc2zvVKWOjiZ3FXSaOvsoqe3ghPxkrBZ3Aovrdm+G3R49vqfWi8DXc2ppcdh4pSZl+MaTzoTeNqjZ7J/gVFDR5K6CzmdljYhAUX7y58diuypwWsJDdrGw0Zga28PipDZeq02m2T7M0MjUKRCdAuWf+i84FTQ0uaug0utwsm5/M9Mz40mI+mJUSGx3JR1ROSD6IwuwJLsepxFerBpmSz6xQO6J0FgCjXv8F5wKCvqbooLKG9tq6OpzsmjSF712m6OT6J66cV+SGSgz0s45aS283ZBIVc8wQyNzF7qT/IYn/BecCgqa3FVQeXpNOckx4RSmfTHcMaVlK4KhQ5P7Ib6c1UCYxfBc1TClqsgESJ8Jm57RjbTHGU3uKmjsqe/gs7ImTshPxjJgrfbUFvcNwY4oTe4DJYY5uSSjiU+b4yntHGbVyLzF0FkHu173X3Aq4DS5q6Cx7LNybBZhQV7iIcdTWzbTFZGG0zrCsrfj0CUZTcTbHDx7YJjee/oMiMuGTU/7LzAVcJrcVVDodThZvr6Sc2dmEBc5oIZsXKS2bNZe+xFEW11ckdnItvYYtrdHDd1ILDD7y1D6ts5YHUc0uaugsKq4luYuO0sXHrr9YnznPiLsbVpvH8a5aS0k2hy8MFztffaV4LLDjtf8F5gKKE3uKig8u6acnKQoTpl86NC+g/V2HSlzZOEWwxVZjWzviGZb+xHWfM8+DpIKoPgl/wanAkaTuwq4soZOPtnbyNKFeVgsh06pT23eTG9YPD3hw4znVpyd2kJSmJ0XqlIxQy0qJuIuzex9Hzrq/R6f8j9N7irgln1WjtUiXH384b3z1JbNNCbOdScndUThFsMVmY3sHK73PvtKME7Y8Yp/g1MB4dGqkCJyAfAAYAUeNcb89gjtrgJeAE4wxqzzWpQqZPU6nLywvpKzp6eTHn/oaJgwexsJHXvYn3VBgKIbW85KbWVFTQrLq1P51ryOQ19c9xgYA7GZ8MmfQY6wbIGuHhkyRuy5i4gVeAi4EJgJLBWRmUO0i8O9OfYabwepQtcb22po6uzjuhMnHvZaav/kpYbEeQGIbOwJtxgu6++9r20YYtaqCGTPh6a90NPq/wCVX3lSllkIlBpj9hpj+oBlwOVDtPslcC/Q48X4VIh7Zk05ecnRh91IBXdJxoWFxsQ5AYhsbDoztZU4m4O/7DxCaSZzHmCgttivcSn/8yS5TwAqBjyv7D/2ORE5Dsg1xug4K+Wx0rp21pQ1DXkjFSC1eROtcZNx2GICEN3YFGExXJjezDs1EexsHaL0EpfpXimydpv/g1N+5UlyH+pO1uf340XEAtwP3D7iiURuEpF1IrKuvl7v2I93T68pJ8wqXF00xDBH4yKldSsNifP9H9gYd35aMzE2Fw/vGuKPoghkzIaG3eDo9X9wym88Se6VwMDtb3KAqgHP44DZwHsisg84EVghIkWDT2SMecQYU2SMKUpL03W5x7Meu5MX11dywewsUmMjDns9oWMP4Y4OGpK03j5asTYX1xb0sKIigorOIX7FM2aDywH1u/wfnPIbT5L7WmCKiBSISDiwBFhx8EVjTKsxJtUYk2+MyQc+BS7T0TJqOK9tqaatx8G1g2akHpTavAlAb6YepRundmEBHt09RO09eRKERWlpJsSNmNyNMQ7gVmAVsAN43hhTLCL3iMhlvg5QhR5jDE9+so9JaTGcOGDd9oFSWzbTE5ZEe/TQyV8NLzPKxZcm9vDcviiaewdVVi1WSJsJdcVgXIEJUPmcR+PcjTErgZWDjt11hLZnHHtYKpRtrGhhc2Ur91w+CznC5KS0ls3ukoxOXjoqa8qaOCG6k+edBfx+g+HyzKZDXk+25DOlbz3FxVvoiPniD+geZzkA1y7SP6pjnc5QVX732Op9xEXauHLB0OvFRPQ2Ed+5T0syxygvqpfZcZ2sqkvCOWhJgtbYQlxiIald6+6hSpO78qua1h5e31rNV4pyiYkY+oNjaou73l6fdJw/QwtJF6Q302gPY11L3CHHndZI2qPzSOzQvVVDlSZ35VdPfroPlzF846T8I7ZJa96EU8JoTJjtv8BC1PEJHaSH97GyLumw11pjJxPdW0e4vS0AkSlf0+Su/KbH7uSZNeWcMyOD3OQjzKAE0po30pQwE5f18CGSanQsAuenN7OzI5qyrkO/ny2xkwFI6CgNRGjKxzS5K795ZdMBmrvs3HBywZEb2XtIbi3WkowXnZnSSoTFxeuDeu/dEWn02eK0NBOiNLkrv3C5DH/9qIzpmXFHHP4IQNVGrMZOQ5LOTPWWGJuL01NaWd0UT5tjwJIEIrTETia+Yy9inIELUPmEJnflF29ur2V3bQffOaPwiMMfAaj4FIB6XXbAq85La8ZhLLzfmHDI8Za4ydhcvcR2VQYoMuUrmtyVzxljePDdEvJTorl4TtbwjcvX0BY9kd6IFP8EN07kRvUxLaaLt+sTDtmpqS2mAINo3T0EaXJXPvf+7nq2HWjj386YjM06zI+cMVCxhnotyfjEOWktVPdGUNzxxc1s95DIXK27hyBN7sqnjDH8zzulZCdEcsVxE4Zv3FAC3U3UJy3wT3DjzIlJ7cRYnbxVn3jI8dbYycT01BBm7zjCO9VY5NHyA0qNxjNryj9/vLe+g/X7m7l0XjbL1x9e1y0sf+Hzx2lNG5gERPXUHnJceUe4xXBaSitv1ifRareSEOa+idoSW0hu3TvEd+4NcITKm7TnrnzGGMM7u+qIi7BRNPHwSTSDxXftp88WQ0+41tt95ZzUFpxGeG/AjdWuyEzs1mgSOjS5hxJN7spndtd2sLe+k9OmphE2XK0dwBjiOvfTHj1RFwvzoZyoPqbHdvF2QyKugzdWRWiNKSChcw+H3G1VY5omd+UTTpfh9W3VpMSEs2i4ce39IuwtRDja3Mld+dQ5qS3U9oazfcCN1dbYQsIdnSS27w5gZMqbNLkrn1i3v4m69l7On5WJzTLyj1l85z4A2mI0ufvaoqR2oq1O3mv4ojTTFjsJgMyGTwIVlvIyTe7K63rsTt7aUUd+SjSzsuM9ek9c537s1mi6I3T7RV8LtxhOTmrj0+Y4upzuFNAXFk9XRBpZDR8HODrlLZrcldd9sLuezl4HF83JGn426gDxXfvduy5pvd0vzkhtxW4sfNz0xVLArTGTSG9ej9XZE8DIlLdocldeVVrXwUelDczLSSAn6cgrPw4U3tdChL2Vtph83wanPlcY3UNOZC/vNn4x5r01dhJWVx9pTRsCGJnyFo+Su4hcICK7RKRURO4Y4vVbRGSriGwSkY9EZKb3Q1XBzuky/Hj5ZsKsFi4aaZmBAeK79gNab/cnETgztYXSzigqusMBaI+ZiFPCyGpYHeDolDeMmNxFxAo8BFwIzASWDpG8nzHGzDHGzAfuBf7o9UhV0HtsdRkby1u4dF4WcZFhHr/PXW+Pojsi3YfRqcFOTW7Divl8zLvLEk598gKtu4cIT3ruC4FSY8xeY0wfsAy4fGADY8zArVxiAB0sO86UNXTy+1W7OGdGOvNyEkd+w0HGEN+5T+vtAZAQ5mRBYgcfNibg6P+NrU49icSOUmirDmxw6ph5ktwnABUDnlf2HzuEiHxXRPbg7rl/f6gTichNIrJORNbV19cfTbwqCDmcLv5j+WbCbRZ+/aU5Ht9EBYjoayLS3kJrbKEPI1RHcmZKK60OGxtbYwF3cgdgzzsBjEp5gyfJfajf1MN65saYh4wxhcB/Aj8b6kTGmEeMMUXGmKK0NB3yFiru+9du1u5r5u5LZ5ERHzmq9x5cjbC1f5y18q95CR3E2xyfr/PeEjeV7vAUTe4hwJPkXgnkDnieA1QN034ZcMWxBKXGjjeLa/jLe3tYujCXK4/PGfX7Ezr30hOWRG/4yLNYlffZxF1739AaS7vDAmKhJnUx7H0XXK5Ah6eOgSfJfS0wRUQKRCQcWAKsGNhARKYMeHoxUOK9EFWwKmvo5PbnNzM3J4GfXzpr1O+3uOzEd+7TXnuAnZbSitMIq5vcE86qU0+Crkao2RzgyNSxGDG5G2McwK3AKmAH8LwxplhE7hGRy/qb3SoixSKyCfgh8A2fRayCQkevg+88tR6rVfjzVxcQGWYd+U2DpLRsxurq03p7gOVH95If1fN5aaYmdbH7hdK3AxiVOlYeredujFkJrBxNkpkwAAAaMklEQVR07K4Bj//dy3GpIGZ3uvjOU+spqevg8RtO8Hiy0mBZDR9jEJ28FAROS2nlicoMatp6yIxPhcw5sOddOO1HgQ5NHSWdoapGxRjDf764hQ9LGvivL8/h1ClHf2M8q+FjOqJzcFpHdxNWed8p/WPeN+5vdh8oPMu9WXlP2/BvVEFLk7salfve3M1LGw7wg3Omck1R7shvOIKIvmaSW7fTGqP19mCQEObkuIQONlW04HQZmHwuuBxQ9n6gQ1NHSZO78thfPyrjwXdLWbowl++fPfmYzpVVvxrB0KL19qBxekor7b0OSuvaIe9ECI+DkjcDHZY6SprclUeeW1vOL1/bzoWzM/nl5bNHNVFpKDm179AdkUpn1AibZiu/WZDQQXS4lQ3lLWANg8IzoeQt3Z1pjNLkrka0YnMVd7y0ldOnpvHAkuOwjbRl3giszh6yGj6iMv1MXXIgiNgsMC83ke3VbbR09cGU86C9CmqLAx2aOgoejZZR48sza8o/f7yzuo2n1uxnYnIMZ05LZ/n6ymM+f0bjGsKc3VRmnE1093Dz4ZS/HZ+XxCd7Gnl1cxVfm3WO+2DJm5A5O7CBqVHTnrs6or31HTzzWTlZCVF8ffFEwm3e+XHJrX2bPlsstSkLvXI+5T1ZCZFkxke6/4jHZ7mHRJa+Feiw1FHQ5K6GVNncxROf7ic5JpzrT8o/qklKQxHjZELte1SlnYrL4vmywMo/RIQFE5PYXNnK7tp2d2mm/FPobgl0aGqUNLmrw9S19/D4x/uICbdyw8kFxER4r3qX2ryJSHszlRlnee2cyrvm5yZiswgvrq90J3fj1IXExiBN7uoQDR29/P3jfYgI3zy5gIQo7/auc2vfxilhVKWd6tXzKu+JjbBx5vR0Xtp4AEfWAohKgt1vBDosNUqa3NXneuxOvv3EOjp6HXz9xImkxEZ49fxinOTWvElN6kk4bDFePbfyrquOz6G+vZcP9jTD1Avdyd1pD3RYahQ0uSsAXC7DD57bxKaKFq4pyiU3+ejWixlOeuNaYnpqKZtwidfPrbzrzGnppMaG89zaCph+MfS0wn7dW3Us0eSuAPjT2yW8vq2Gn140g1nZCT65RkHVq/TZ4jiQfoZPzq+8J9xm4coFOby9o466jJPAFgU7/xnosNQoaHJXvLuzjv9+u4Srjs/hxlMKfHINm6OL3Jp/UZ51ni4UNkZcc0IuDpfhxS3N7oXEdv5TZ6uOITqJaZyraOrituc2MSMrnl9dcezLChxUWP7CIc9TWzYT5uym1xp32GsqOBWmxbIwP5nn1pZzyzkXIbv+CdWbIPu4QIemPKA993Gsx+7k357egMsY/ve6o9tww1OpLVvoCUukI/roV5JU/rdkYS77GrtYH7EIxKKlmTFEe+7j2G9f38nWA608+vUiJqb4bvRKeF8r8Z1lHEg7TdeSCXKff6qyuve0vdABPw9L5en3tlCUVAAbn4T4CVB0QwCjVJ7Qnvs49VFJA49/vI/rT8rnnJkZPr1WRvM6QKhPnO/T6yjvi7LBFXk9rKyMoDt9PrTXQHt1oMNSHvAouYvIBSKyS0RKReSOIV7/oYhsF5EtIvK2iEz0fqjKW1q77fx4+WYK02K448LpPr2WxWUnrXkDzXHT6AtP9Om1lG98Jb+HXpfwD/uJgMCBDYEOSXlgxOQuIlbgIeBCYCawVERmDmq2ESgyxswFlgP3ejtQ5T2/WFFMXXsvf7xmvk/r7AApLVsJc3ZTo4uEjVmzkxzMT7bzf+UZmNSpULVBR82MAZ7U3BcCpcaYvQAisgy4HNh+sIEx5t0B7T8FrvNmkGp0Bi7ZO9i2A628tPEAZ01Pp7iqjeIqH+6RaQyZTWvojMykPVo/zI1lXy/s5odr4ymZcAJTG55y995zjg90WGoYnpRlJgAVA55X9h87khuB14d6QURuEpF1IrKuvr7e8yiVV/TYnby6pYrshEjOnJbu8+vFd5YR3VtPTfJCvZE6xl2U00NKhIsHWxaDxQpbdThrsPMkuQ/1WznkZzIRuQ4oAn4/1OvGmEeMMUXGmKK0tDTPo1Re8eb2Gjp6HFxx3ASsFt8n26zGT7Bbo2lM0I0exrpIKywp6Oa1mkS6kmdB8UvgcgY6LDUMT5J7JTBwcHIOcNj2OSJyDvBT4DJjTK93wlPeUtHUxZq9TZxYmEJOkvfXjRkstqucxI49VKcsxlh0xG0o+OqkbgBWsRg6amHfhwGOSA3Hk+S+FpgiIgUiEg4sAVYMbCAixwEP407sdd4PUx0Lp8vwj00HiIu0ce4M3w57PCin7j3s1hhqU07wy/WU72VHuzg3u5ff1S3CRMTDxqcDHZIaxojJ3RjjAG4FVgE7gOeNMcUico+IXNbf7PdALPCCiGwSkRVHOJ0KgE/2NFDd2sMlc7N9PjoGIL3xMxI691GVdjIuS7jPr6f85xuTu6npi2RP5sWw/RXoagp0SOoIPPq8bIxZCawcdOyuAY/P8XJcykvae+y8vbOOaRlxzMqO9/0FjWFuyYP02eKoTSry/fWUXy1OszM9wcEfGk/kf53Pwpbn4MTvBDosNQSdoRri3txei8NpuHhOltcWBRtOTt07pDdv5EDaqVprD0EicPPUTt5oSKM1eS6sf1zHvAcp/e0LYZXNXWzY38wpk1NJjfPurkpDsTq7WbDjXlpiJ1OXtMDn11OBcUluL38odvJ0x/H8W99j8K+fQ/IRlorWNWgCRnvuIcoYw2tbqonu3w/TH2bteZTY7irWzvqpewVBFZLCLPDtqV082Haqe23+8o8DHZIagv4GhqjNla2UN3Vx/swMv9xEje0sZ8bexyjLvpj6ZK21h7pr8ruJCA/nA9tiqNoIve2BDkkNosk9BNmdLlYV15CdGMmCiUm+v6BxsbD4HlyWcDZOu93311MBF22Db0zu4p62SzAup455D0Ka3EPQx6UNtHbbuWh2FhY/3ESduv9ZMhvXsGHGj+iJ1JnH48U3CrupsWSyJXw+7PsIHDp3MZhocg8xjR29vLe7nhmZcUxKi/X59eI79jJ/1/0cSDuVPTlX+vx6KngkRRi+OaWLX7RfDvYuqPgs0CGpATS5h5gH3i7B7nRx/uxMn19LXHZO3PIznNZI1sz+hS4ONg7dNLWLPbbJlFgLoew9MK5Ah6T66VDIELKnvoOn15RzQn4y6XGRPr/ecbvuJ7V1Kx/Ov0/LMSFmTZnnM08vTTf8ofoyHnbeT+nmj2hMmM2igmQfRqc8oT33EPLb13cSFWblbD+sH5NbvYrp+55k18SvUpF1ns+vp4LX+enNrLPOZQ85TKh7DzG6WmQw0J57iFizt5F/ba/lx+dPIzZi9P9bP98YeQR78q4mvmMvJ269i/rEeWycrqNjxrtwi+HK7CZ+VbGUx/g9ac2bgHMDHda4pz33EOByGX6zcgeZ8ZF88+QjzBT0kjB7K6dt+Hec1kg+mv8HXJYwn15PjQ1npLSyK3wGm80UsuvfB2dfoEMa9zS5h4DXtlazubKV28+bSlS47yYsiXFy6sbbiemq5MPj7qc7yvc3bdXYYBW4MbeOX/UtJcLRAWUfBDqkcU+T+xjX63By7xs7mZ4Zx5cX5PjuQsYwsXoVmY1r+Gz23dQn69ox6lCz47uITMriXed8nCVvQU9roEMa1zS5j3FPfLyfyuZufnLRDJ9unZfZtIaM5nVsL7iBspzLfXYdNbZdl1PHva7rcDodmOKXAx3OuKbJfQxr7uzjf94p4dQpqZw21XdDEZPadpFX8yZN8TPYNO02n11HjX2JYU5OmSD8j/0KpHoTlLwV6JDGLU3uY9if3tpNR6+Dn10802fXiOmuorDyJTqjstkz4Qpd7VGN6KzUVtYmnE+ZycLx6g+gryvQIY1LHv2misgFIrJLREpF5I4hXj9NRDaIiENErvJ+mGqw0rp2nlpTztKFeUzLjPPJNcLtbUwtX4bDFs2uvCU6MkZ5xCLw+0Xd3OP6Jra2cpyrfhLokMalEZO7iFiBh4ALgZnAUhEZ3FUsB64HnvF2gGpov/7nDqLDrPzw3Kk+Ob/F1cfU8mVYXX3syluKw+b7dWpU6MiNcXHV8Tn8r+NSrOsfA62/+50nPfeFQKkxZq8xpg9YBhxyR80Ys88YswXQhSX84P3d9by7q55bz5pMSqwPdlgyhsLKfxDdU0tJzpV0R/pnsw8VWi7O6aVywe1scE3G8Y9boaks0CGNK55MZZwAVAx4Xgks8k04aiR2p4tfvbadvORorj853yfXyK7/kOT2nezPOI/WuCmHvObpTFalAH526VxuKruDB9v/nYinryXi26sg0g8btSuPeu5Dja87qh1xReQmEVknIuvq6+uP5hTj3mOryyip6+BnF88gwub9CUsJ7SXk1L9HQ8IcalL0b7g6NpFhVn51/cX81HIb1sZd9D1zLTh09qo/eJLcK4HcAc9zgKqjuZgx5hFjTJExpigtTVcRHK3q1m7+9FYJZ09P59yZ3l8cLKKvicmVL9MVmUlZ9iW6hK/yiryUaL79zZv4metmwss/xPHSLeDSCq6veZLc1wJTRKRARMKBJcAK34alhvLL17bjdBnuvmwW4uXEKy4HUyqWg0BJ7tU6MkZ51dycRM679jZ+51iCbfuLOJZ/S3vwPjZicjfGOIBbgVXADuB5Y0yxiNwjIpcBiMgJIlIJXA08LCLFvgx6PHp/dz0rt9Zw65mTyU2O9vr582r/RUxPDXsmXE5vuB/2XVXjzlnTMyi4/Kf81rHUneCfvgZ6OwIdVsjyaG1YY8xKYOWgY3cNeLwWd7lG+UB3n5Ofv7KNgtQYbjp9ktfPn1v9JplNa6lOOZGWuGleP79SB11zQh6vhP0/7lgex6/LHsXx1/OxLXkKkn27mul4pOu5B6Fn1pQf8vy1LVXsa+zixlMKeHH9Aa9eK6arkkXbfk571AQqMs726rmVGsrl8ycQE/4jbnk2iT/WPUj0w6dj/fLDMO3CQIcWUjS5B7m99R18vKeREyelUOjlDa/F5eCkzXcAQmnOlRjx3XLBSg10zswMYq6/ma88lcsfeu9j5rNLIP80mHEJWMOP/MaiG/wX5BinC4UEsV67kxc3VJISE84Fs7y/dvqc0r+Q1rKZz2bfRV94otfPr9RwFhem8PD3ruSOxPv4m+MC2PcB5sP7oLVi5DerEWlyD2Kvb6uhpcvOVcfnEG7z7v+q9Ma1zNrzf+yZcAXlWRd49dxKeSo3OZrn/u0MNmddw3V9d9Lc2Yv56H7Y/Qa4dC/WY6FlmSC1pbKFz/Y1cerkVCamxHj13OF9LZy0+Q7ao/NYP/NOr55bqdGKCrfyp4VtvLh/Ipds+h0/sfydS3a/QUf5FvbkXEFPROrnbfc4y4c50xeuXZTnq3DHDO25B6Gath5e2nCAvORozp3l5clKxrBo611E9DWxev69OGzeH1ap1GiJwFX5PTx3bi9PxN/Ed/r+HVdPK7P2/B9pTevBHNWk+HFNk3uQae228/Sn+4mwWbh2YR42i3f/F00pf47cunfZNO0HNCf4bh14pY5GboyLZ09vIWvCRC7p+y/WOKcxqfqfTK54AauzO9DhjSma3IOI02W4/flNNHf1ce2iPOKjvDtLNLm1mAU77qUq7RR25V/n1XMr5S1WgfPTW7hzVgv/E/Vv/Np+LQltJUwreZSUli2BDm/M0Jp7kHC5DD95aStv7ajj0nnZPqizt3LKxh/SE5HKJ3N/ozsqKf9Y99hRvzU53MEPJtewqXUet1RM4m7H/3L2J19n9aTbODD1G7r20Qj0NzwIGGP4xavFPLeugu+fPYXFk1K8fAEXJ225g6ieOj467j5dXkCNKfMTOrlhpoXHkn/Ie675nL73Pgrf+y6OrtZAhxbUtOceYMYYfvvGTv7+yX5uOm0SPzhnCs9+5t1xvvN2/zfZ9R/x2cyf0Zg4x6vnVsofbBY4L7uHltRLeX5/Hl/ufoWM965iZdJ1zM2KIcwy6IbrotsDE2gQ0Z57APU6nPx4+RYefn8vXztxIndeON3rqz0Wlr/ArL1/pST3akrzrvHquZXyt8RwFxOnzOH97BuJkx5ubH6AD4r38VFjHC4dUHMI7bkHSF17D7c8uZ4N5S18/+wp3Hb2FK8n9qz6D1lY/EuaYyfTFDeNworlXj2/UoESm5zFvvgbSd+/gp/0/J3XD5zAr2u+xkUTeliQoCtNgib3gPhgdz3/+eIWWrrs/PmrC7hoTpbXr5FVv5pTN/yArsgMSnOu0huoKuQ4bDFUTVqCq+ETzqt7l+NNKT8su4WXoybzs931nDYl1esdprFEk7sf1bX18Mt/7uDVzVUUpMaw/DtFzMpO8Pp1JtS+yykbb6c1tpCy7EtwDbcQk1I+sKasyT8XEqEm7STaYvMprHyZp+S/eMV+Kt/7m5OCnAncetYUzp6ejsUy/pK8Jnc/aOjo5YlP9vPY6jJ67S5uO2cKt5xeSGSYl1dhNIbJ5c9TtOO3NMXP4N0T/kJe9ZvevYZSQagrKptthTeTU/8+lzV8xHnx23mgbQk3P7GYSenx3HByPl8+Loeo8PGz8qkmdx8xxrC5spXn1lbw4oZK7E4X58zI4M4LpzPJy0v3AtgcXSzcdjf51a9zIO1UVs+7F0eY96+jVLAyFhsVGWfTGD+TOX2buaPyIW7OeJf7XV/hpy+38/tVu7hqQQ7XnJDL1Iy4QIfrc5rcvaiz18GG8mbe31XP69tqONDSTbjNwpULcvjWqQVeX48dAGPIrfkXx+36I9Hd1Wya+n22T7pRa+xq3OqKyoKv/R62vUjSW3dzT/vP+c8Jc3kq7Mvc/3EPj35UxrzcRL40P5vzZmWSnRgV6JB9QowHC/KIyAXAA4AVeNQY89tBr0cATwDHA43AV4wx+4Y7Z1FRkVm3bt1Rhh143X1OSus62F7dyo7qdjZWtLDtQCtOl8FqEaakxzI7O4EZWfE++SgoxklW/WpmlD1GRtM6WmIns3bWT6lPLjqkXWH5C16/tlLBbtHV/ePcHX2w+Vn46I/QvA9XTAZb0i7lL43zWVWfBAjzchI4a3oGiwtTmJ+b6PXltb1NRNYbY4pGajdiz11ErMBDwLlAJbBWRFYYY7YPaHYj0GyMmSwiS4DfAV85utCDg93poq69l5rWHg60dFPZ3EVlczfljV3sre+gqrXn87bR4VZmZcdzy+mTWFiQwt66DiK8XU8HxGUntWUzWQ0fk1/1T2K7q+iOSOWzmT9jT+6VGIt+EFPqELZwOP4bMP+rUPImlvWPMb/krzyMoS+jgB1xi1nZNokn357A/W8lEBVmZV5uAnMmJDB7grtzlpcc7f37Y37gSTZYCJQaY/YCiMgy4HJgYHK/HLi7//Fy4EEREePJxwIvMMbgMu6FtxwuFw6XweE09Dlc7i+nkx67i64+J119Djp7nXT2OmjvddDabaelq4+WLjuNnb00dvTR0NFLY2ffYauMJkWHkZcczaJJKRSkxjA5PZaZ/f/zB96NP9A8itXrjAuLy4HV1ev+cnYTbm8n3N5GVG890T21xHZVktheQmJHCTZnNy6xUpdcxMZpP6Qy4yyMxbsLjCkVcqw2mH6R+6utGnatJHzna8zb/xLzHD3cGQE9kWkcCMtnZ1MG2yvieMeZxIvE0UoskTFJJCUlkpiQQFJcLMnxMSRERxEfHU58lI3ocBvR4VaiwqyE2yyEWy2E2SzYLILVItgs4vdhmZ4k9wnAwPnwlcCiI7UxxjhEpBVIARq8EeRAj364l3vf2IXLGAy4/z2GPyEikBAVRkJUGMkx4eQmR3NcXiLpcZFkJkSSGR9JVmIkOUnRxEYce8/4qn+dhNXZg2AQ40QYOfiesCRa4qawJ+fL1CYXUZuyEHtY/DHHotS4FJ8FJ9zo/nL0QtUmqFhDZN0OCuu2U9j8IRdbW91F6IPsQF3/1wBOIzixcErvf1PHyGs2WQSsFuEXl832+YYinmSrof7cDM5InrRBRG4Cbup/2iEiuzy4fkj56hcPU/H4j18bsB94ywcRBZ1RfF/GFf2+DO0I35cf+TmMr42q9Vd/c0guGK2JnjTyJLlXArkDnucAVUdoUykiNiABOGwWgzHmEeARTwILdSKyzpObIuONfl+Gpt+Xoen35cg8uS28FpgiIgUiEg4sAVYMarMC+Eb/46uAd/xVb1dKKXW4EXvu/TX0W4FVuKtQfzPGFIvIPcA6Y8wK4K/AkyJSirvHvsSXQSullBqeR3cIjTErgZWDjt014HEPcLV3Qwt5Wp4amn5fhqbfl6Hp9+UIPJrEpJRSamwJ7qlYSimljoom9wAQEauIbBSR1wIdSzARkX0islVENonI2F2bwstEJFFElovIThHZISKLAx1ToInItP6fk4NfbSJyW6DjCiY6Xz0w/h3YAehMpMOdaYzR8dyHegB4wxhzVf+ItehABxRoxphdwHz4fImUA8DLAQ0qyGjP3c9EJAe4GHg00LGo4Cci8cBpuEekYYzpM8a0BDaqoHM2sMcYsz/QgQQTTe7+9yfgPwBXoAMJQgZ4U0TW989mVjAJqAce6y/lPSoiMYEOKsgsAZ4NdBDBRpO7H4nIJUCdMWZ9oGMJUicbYxYAFwLfFZHTAh1QELABC4C/GGOOAzqBOwIbUvDoL1NdBuja1oNocvevk4HLRGQfsAw4S0SeCmxIwcMYU9X/bx3u+unCwEYUFCqBSmPMmv7ny3Ene+V2IbDBGFMb6ECCjSZ3PzLG3GmMyTHG5OP+KPmOMea6AIcVFEQkRkTiDj4GzgO2BTaqwDPG1AAVIjKt/9DZHLrc9ni3FC3JDElHy6hgkQG83L/mtQ14xhjzRmBDChrfA57uL0HsBW4IcDxBQUSicW8idHOgYwlGOkNVKaVCkJZllFIqBGlyV0qpEKTJXSmlQpAmd6WUCkGa3JVSKgRpcldKqRCkyV0ppUKQJnellApB/x9LwVIz9ZFbxAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "target_0_rep2 = rp2_train_pd.loc[rp2_train_pd['group'] == 0]\n", "target_1_rep2 = rp2_train_pd.loc[rp2_train_pd['group'] == 1]\n", "\n", "sns.distplot(target_0_rep2[['value']], hist=True, rug=False)\n", "sns.distplot(target_1_rep2[['value']], hist=True, rug=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This time, the two distributions don't entirely overlap" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Although the values change to remove bias, the ranking with the groups remain the same. We can test this looking at the top three values for each group and each dataset\n", "\n", "First let's focus on the unprivileged group, viewing the smallest five values" ] }, { "cell_type": "code", "execution_count": 186, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupoutputvalue
701003.872637
412014.215382
175014.233313
1060004.236066
660014.309020
\n", "
" ], "text/plain": [ " group output value\n", "701 0 0 3.872637\n", "412 0 1 4.215382\n", "175 0 1 4.233313\n", "1060 0 0 4.236066\n", "660 0 1 4.309020" ] }, "execution_count": 186, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_0_orig.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So in ascending order, the smallest values are at index 701, 412, 175, 1060 and 660\n", "\n", "We can check if they've changed when repairing our data" ] }, { "cell_type": "code", "execution_count": 187, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupvaluelabels
7010.03.8726370.0
4120.04.2153821.0
1750.04.2333131.0
10600.04.2360660.0
6600.04.3090201.0
\n", "
" ], "text/plain": [ " group value labels\n", "701 0.0 3.872637 0.0\n", "412 0.0 4.215382 1.0\n", "175 0.0 4.233313 1.0\n", "1060 0.0 4.236066 0.0\n", "660 0.0 4.309020 1.0" ] }, "execution_count": 187, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_0_rep1.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Neither the order nor the values have changed. If it hasn't changed with a `repair_level` of 1.0, we don't expect it to change on the second repaired data but will check anyway" ] }, { "cell_type": "code", "execution_count": 188, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupvaluelabels
7010.03.8726370.0
4120.04.2153821.0
1750.04.2333131.0
10600.04.2360660.0
6600.04.3090201.0
\n", "
" ], "text/plain": [ " group value labels\n", "701 0.0 3.872637 0.0\n", "412 0.0 4.215382 1.0\n", "175 0.0 4.233313 1.0\n", "1060 0.0 4.236066 0.0\n", "660 0.0 4.309020 1.0" ] }, "execution_count": 188, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_0_rep2.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we wish to look at the smallest five values for the privileged group" ] }, { "cell_type": "code", "execution_count": 189, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupoutputvalue
401104.861829
546105.015805
294115.020209
394105.064415
326105.091649
\n", "
" ], "text/plain": [ " group output value\n", "401 1 0 4.861829\n", "546 1 0 5.015805\n", "294 1 1 5.020209\n", "394 1 0 5.064415\n", "326 1 0 5.091649" ] }, "execution_count": 189, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_1_orig.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In ascending order we have indexes 401, 546, 294, 394 and 326\n", "\n", "We then compare this with the first repaired dataset" ] }, { "cell_type": "code", "execution_count": 190, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupvaluelabels
5461.03.8726370.0
4011.03.8726370.0
2941.04.2153821.0
3941.04.2153820.0
3261.04.2153820.0
\n", "
" ], "text/plain": [ " group value labels\n", "546 1.0 3.872637 0.0\n", "401 1.0 3.872637 0.0\n", "294 1.0 4.215382 1.0\n", "394 1.0 4.215382 0.0\n", "326 1.0 4.215382 0.0" ] }, "execution_count": 190, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_1_rep1.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have the same five indexes but as 546 and 401 have the same value, the order has slightly change. Also note how these values match the lowest in unprivileged group\n", "\n", "We check with the second repaired dataset" ] }, { "cell_type": "code", "execution_count": 191, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
groupvaluelabels
4011.04.3746180.0
5461.04.4514060.0
2941.04.4771121.0
3941.04.4809530.0
3261.04.5176330.0
\n", "
" ], "text/plain": [ " group value labels\n", "401 1.0 4.374618 0.0\n", "546 1.0 4.451406 0.0\n", "294 1.0 4.477112 1.0\n", "394 1.0 4.480953 0.0\n", "326 1.0 4.517633 0.0" ] }, "execution_count": 191, "metadata": {}, "output_type": "execute_result" } ], "source": [ "target_1_rep2.sort_values(\"value\")[:5]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, the ranking has stayed the same" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Building an Unbiased Model\n", "\n", "Now that we've repaired the values to remove bias, we can build a new model and confirm if the Disparate Impact has been removed\n", "\n", "We build another logistic regression model, this time with the repaired data (`repair_level`=1.0)" ] }, { "cell_type": "code", "execution_count": 114, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] }, { "data": { "text/plain": [ "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n", " intercept_scaling=1, max_iter=100, multi_class='warn',\n", " n_jobs=None, penalty='l2', random_state=None, solver='warn',\n", " tol=0.0001, verbose=0, warm_start=False)" ] }, "execution_count": 114, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr_di=LogisticRegression()\n", "lr_di.fit(rp_train_pd[[\"value\"]], rp_train_pd[\"labels\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can score our test data comparing with the previous accuracy of 71%. If we believe that the labels are biased (which they are as we created them that way), then the accuracy we see at the beginning is untrue. Although it is our ground truth, we are capturing unfairness that we wish to remove. Consequently, we expect our accuracy value to be reduced. Also, in a business context, we may need to consider other performances metrics (precision, recall, F1 score)." ] }, { "cell_type": "code", "execution_count": 115, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.675" ] }, "execution_count": 115, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr_di.score(rp_test_pd[[\"value\"]], rp_test_pd[\"labels\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we create a dataframe with the predictions from our test data and their group assignment" ] }, { "cell_type": "code", "execution_count": 116, "metadata": { "collapsed": false }, "outputs": [], "source": [ "di_preds = lr_di.predict(rp_test_pd[[\"value\"]])" ] }, { "cell_type": "code", "execution_count": 117, "metadata": { "collapsed": false }, "outputs": [], "source": [ "di_pred_df = pd.DataFrame({\"group\": rp_test_pd[\"group\"], \"preds\": di_preds})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We again calculate the proportion of each group receiving the favourable outcome" ] }, { "cell_type": "code", "execution_count": 118, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.925\n" ] } ], "source": [ "di_lr_pr_unpriv = calc_prop(di_pred_df,\"group\",0,\"preds\",1)\n", "print(di_lr_pr_unpriv)" ] }, { "cell_type": "code", "execution_count": 119, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.925\n" ] } ], "source": [ "di_lr_pr_priv = calc_prop(di_pred_df,\"group\",1,\"preds\",1)\n", "print(di_lr_pr_priv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then calculate the Disparate Impact" ] }, { "cell_type": "code", "execution_count": 120, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.0" ] }, "execution_count": 120, "metadata": {}, "output_type": "execute_result" } ], "source": [ "di_lr_pr_unpriv / di_lr_pr_priv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have a perfect ratio, the two groups are being treated equally" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We repeat with our second repaired dataset (`repair_level`=0.8)" ] }, { "cell_type": "code", "execution_count": 121, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/staceyro/anaconda/envs/aif360/lib/python3.5/site-packages/sklearn/linear_model/logistic.py:432: FutureWarning: Default solver will be changed to 'lbfgs' in 0.22. Specify a solver to silence this warning.\n", " FutureWarning)\n" ] }, { "data": { "text/plain": [ "LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,\n", " intercept_scaling=1, max_iter=100, multi_class='warn',\n", " n_jobs=None, penalty='l2', random_state=None, solver='warn',\n", " tol=0.0001, verbose=0, warm_start=False)" ] }, "execution_count": 121, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr_di2=LogisticRegression()\n", "lr_di2.fit(rp2_train_pd[[\"value\"]], rp2_train_pd[\"labels\"])" ] }, { "cell_type": "code", "execution_count": 122, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.6892857142857143" ] }, "execution_count": 122, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lr_di2.score(rp2_test_pd[[\"value\"]], rp2_test_pd[\"labels\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our accuracy is still lower than the original data. However, it is higher than the fully repaired dataset. \n", "\n", "This is a challenge for AI practitioners, when you know you have biased data, you realise that the ground truth you're building a model with doesn't necessarily reflect reality.\n", "\n", "Again we create a dataframe with each data points group and prediction" ] }, { "cell_type": "code", "execution_count": 123, "metadata": { "collapsed": false }, "outputs": [], "source": [ "di2_preds = lr_di2.predict(rp2_test_pd[[\"value\"]])" ] }, { "cell_type": "code", "execution_count": 124, "metadata": { "collapsed": false }, "outputs": [], "source": [ "di2_pred_df = pd.DataFrame({\"group\": rp2_test_pd[\"group\"], \"preds\": di2_preds})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We calculate the proportion of the unprivileged group receiving the favourable outcome" ] }, { "cell_type": "code", "execution_count": 125, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.9\n" ] } ], "source": [ "di2_lr_pr_unpriv = calc_prop(di2_pred_df,\"group\",0,\"preds\",1)\n", "print(di2_lr_pr_unpriv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly for the privileged group" ] }, { "cell_type": "code", "execution_count": 126, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.975\n" ] } ], "source": [ "di2_lr_pr_priv = calc_prop(di2_pred_df,\"group\",1,\"preds\",1)\n", "print(di2_lr_pr_priv)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then our final Disparate Impact" ] }, { "cell_type": "code", "execution_count": 127, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.9230769230769231" ] }, "execution_count": 127, "metadata": {}, "output_type": "execute_result" } ], "source": [ "di2_lr_pr_unpriv / di2_lr_pr_priv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Still fair!\n", "\n", "It is also possible to calculate using AIF360's BinaryLabelDatasetMetric class\n", "\n", "First we create a copy of our repaired BinaryLabelDataset and replace the labels with our logistic regression predictions" ] }, { "cell_type": "code", "execution_count": 198, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rp2_test_pred = rp2_test.copy()\n", "rp2_test_pred.labels = di2_preds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Then we create a new BinaryLabelDatasetMetric with this dataset and specifying the privileged and unprivileged groups" ] }, { "cell_type": "code", "execution_count": 199, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bldm = BinaryLabelDatasetMetric(rp2_test_pred, \n", " privileged_groups=[{\"group\": 1}], \n", " unprivileged_groups=[{\"group\": 0}])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then simply call `disparate_impact`" ] }, { "cell_type": "code", "execution_count": 200, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.9230769230769231" ] }, "execution_count": 200, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bldm.disparate_impact()" ] } ], "metadata": { "kernelspec": { "display_name": "Python (aif360)", "language": "python", "name": "aif360" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.6" } }, "nbformat": 4, "nbformat_minor": 2 }