detect_prediction_bias_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predicting utilization of medical services: bias detection\n",
    "\n",
    "In this tutorial we illustrate that empirical quantiles can be used to detect prediction bias. Follow the instructions in `cqr/get_meps_data/` folder to download the data set. This experiments replicates Figure 1 in the equalized coverage paper.\n",
    "\n",
    "The Medical Expenditure Panel Survey (MEPS) 2016 data set, provided by the Agency for Healthcare Research and Quality, contains information on individuals and their utilization of medical services. The features used for modeling include age, marital status, race, poverty status, functional limitations, health status, health insurance type, and more. The goal is to predict the health care system utilization of each individual; a score that reflects the number of visits to a doctor’s office, hospital visits, etc. After removing observations with missing entries, there are $n=15656$ observations on $p=139$ features. We set the sensitive attribute $A$ to race, with $A=0$ for non-white and $A=1$ for white individuals, resulting in $n_0=9640$ samples for the first group and $n_1=6016$ for the second.\n",
    "\n",
    "The code starts by randomly splitting the data into train (80%) and test (20%) sets and standardized the features to have zero mean and unit variance; the means and variances are computed using the training examples. We also transform the response variable by $Y = \\log(1 + (\\text{utilization score}))$ as the raw score is highly skewed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data set name: meps_21\n",
      "Data set dimensions: (n=15656, p=139)\n",
      "A=0 : non-white\n",
      "A=1 : white\n",
      "Samples per group: A=0 (n=9640) ; A=1 (n=6016)\n",
      "Train set (n=12524, p=139) ; Test set (n=3132, p=139)\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "import torch\n",
    "import random\n",
    "import numpy as np\n",
    "np.warnings.filterwarnings('ignore')\n",
    "\n",
    "from cqr import helper\n",
    "from datasets import datasets\n",
    "from matplotlib import pyplot\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# set some arbitrary seed for reproducability \n",
    "seed = 1234\n",
    "\n",
    "random.seed(seed)\n",
    "np.random.seed(seed)\n",
    "torch.manual_seed(seed)\n",
    "if torch.cuda.is_available():\n",
    "    torch.cuda.manual_seed_all(seed)\n",
    "\n",
    "save_figures = 0\n",
    "dataset_name = \"meps_21\"\n",
    "\n",
    "random_state_train_test = seed\n",
    "random.seed(seed)\n",
    "np.random.seed(seed)\n",
    "torch.manual_seed(seed)\n",
    "if torch.cuda.is_available():\n",
    "    torch.cuda.manual_seed_all(seed)\n",
    "\n",
    "# path to meps data set, download by following the instructions in ./cqr/get_meps_data/\n",
    "dataset_base_path = '/Users/romano/mydata/regression_data/'\n",
    "\n",
    "# used to determine the size of test set\n",
    "test_ratio = 0.2\n",
    "\n",
    "# load the dataset\n",
    "X, y = datasets.GetDataset(dataset_name, dataset_base_path)\n",
    "\n",
    "# display basic information\n",
    "print(\"Data set name: %s\" % (dataset_name))\n",
    "print(\"Data set dimensions: (n=%d, p=%d)\" % \n",
    "      (X.shape[0], X.shape[1]))\n",
    "\n",
    "print(\"A=0 : non-white\")\n",
    "print(\"A=1 : white\")\n",
    "print(\"Samples per group: A=0 (n=%d) ; A=1 (n=%d)\" % \n",
    "      (X[X[:,-1]==0].shape[0], X[X[:,-1]==1].shape[0]))\n",
    "\n",
    "dataset_name_group_0 = dataset_name + \"_non_white\"\n",
    "dataset_name_group_1 = dataset_name + \"_white\"\n",
    "\n",
    "# divide the dataset into test and train based on the test_ratio parameter\n",
    "x_train_orig, x_test_orig, y_train_orig, y_test_orig = train_test_split(X,\n",
    "                                                                        y,\n",
    "                                                                        test_size=test_ratio,\n",
    "                                                                        random_state=random_state_train_test)\n",
    "\n",
    "\n",
    "# zero mean and unit variance scaling \n",
    "scalerX = StandardScaler()\n",
    "scalerX = scalerX.fit(x_train_orig)\n",
    "\n",
    "# scale\n",
    "x_train = scalerX.transform(x_train_orig)\n",
    "x_test = scalerX.transform(x_test_orig)\n",
    "\n",
    "y_train = np.log(1.0 + y_train_orig)\n",
    "y_test = np.log(1.0 + y_test_orig)\n",
    "\n",
    "# reshape the data\n",
    "y_train = np.squeeze(np.asarray(y_train))\n",
    "y_test = np.squeeze(np.asarray(y_test))\n",
    "\n",
    "# compute input dimensions\n",
    "n_train = x_train.shape[0]\n",
    "in_shape = x_train.shape[1]\n",
    "\n",
    "print(\"Train set (n=%d, p=%d) ; Test set (n=%d, p=%d)\" % \n",
    "      (x_train.shape[0], x_train.shape[1], x_test.shape[0], x_test.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we fit a neural network regression function on the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test MSE = 0.77782035\n"
     ]
    }
   ],
   "source": [
    "# pytorch's optimizer object\n",
    "nn_learn_func = torch.optim.Adam\n",
    "\n",
    "# number of epochs\n",
    "epochs = 1000\n",
    "\n",
    "# learning rate\n",
    "lr = 0.0005\n",
    "\n",
    "# mini-batch size\n",
    "batch_size = 64\n",
    "\n",
    "# hidden dimension of the network\n",
    "hidden_size = 64\n",
    "\n",
    "# dropout regularization rate\n",
    "dropout = 0.1\n",
    "\n",
    "# weight decay regularization\n",
    "wd = 1e-6\n",
    "\n",
    "# seed for splitting the data in cross-validation.\n",
    "cv_test_ratio = 0.1\n",
    "\n",
    "# ratio of held-out data, used in cross-validation\n",
    "cv_random_state = 1\n",
    "\n",
    "model = helper.MSENet_RegressorAdapter(model=None,\n",
    "                                       fit_params=None,\n",
    "                                       in_shape = in_shape,\n",
    "                                       hidden_size = hidden_size,\n",
    "                                       learn_func = nn_learn_func,\n",
    "                                       epochs = epochs,\n",
    "                                       batch_size=batch_size,\n",
    "                                       dropout=dropout,\n",
    "                                       lr=lr,\n",
    "                                       wd=wd,\n",
    "                                       test_ratio=cv_test_ratio,\n",
    "                                       random_state=cv_random_state)\n",
    "    \n",
    "model.fit(x_train,y_train)\n",
    "predictions = model.predict(x_test)\n",
    "MSE = np.mean((predictions - y_test)**2)\n",
    "\n",
    "print('Test MSE = ' + str(MSE))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we compute the signed residuals of the test samples,\n",
    "$$R_i = Y_i - \\hat{Y}_i,$$\n",
    "where $\\hat{Y}_i $ is the predicted response, and then compute the resulting empirical cumulative distribution functions $\\mathbb{P}\\{R \\leq r | A=0\\}$ and $\\mathbb{P}\\{R \\leq r | A=1\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A=0: median = 0.5127405096203849; upper quantile = 1.5930867; lower quantile = -1.0478824\n",
      "A=1: median = 0.45492142266335817; upper quantile = 1.83082; lower quantile = -1.2507485\n"
     ]
    }
   ],
   "source": [
    "# function that extracts the group identifier\n",
    "def condition(x, y=None):\n",
    "    return int(x[0][-1]>0)\n",
    "\n",
    "# extract groups\n",
    "category_map = np.array([condition((x_test[i, :], None)) for i in range(x_test.shape[0])])\n",
    "categories = np.unique(category_map)\n",
    "\n",
    "# compute quantiles\n",
    "signed_residual = (y_test - predictions)\n",
    "signed_residual_0 = signed_residual[category_map == 0]\n",
    "signed_residual_1 = signed_residual[category_map == 1]\n",
    "\n",
    "# compute alpha/2 and 1-alpha/2 quantiles\n",
    "def compute_lo_hi_q(residuals, alpha):\n",
    "    err_high = np.reshape(residuals, (residuals.shape[0],1))\n",
    "    err_low = np.reshape(-residuals, (residuals.shape[0],1))\n",
    "    residuals = np.concatenate((err_low,err_high),1)\n",
    "    residuals = np.sort(residuals,0)\n",
    "    index = int(np.ceil((1 - alpha / 2) * (residuals.shape[0]))) - 1\n",
    "    return -residuals[index,0], residuals[index,1]\n",
    "\n",
    "# compute 0.05th and 0.95th quantiles per group\n",
    "lo_0, up_0 = compute_lo_hi_q(signed_residual_0, alpha=0.1)\n",
    "lo_1, up_1 = compute_lo_hi_q(signed_residual_1, alpha=0.1)\n",
    "med_0 = sum(signed_residual_0<=0)/len(signed_residual_0)\n",
    "med_1 = sum(signed_residual_1<=0)/len(signed_residual_1)\n",
    "\n",
    "print(\"A=0: median = \" + str(med_0) + \"; upper quantile = \" + str(up_0) + \"; lower quantile = \" + str(lo_0))\n",
    "print(\"A=1: median = \" + str(med_1) + \"; upper quantile = \" + str(up_1) + \"; lower quantile = \" + str(lo_1))\n",
    "\n",
    "# compute CDFs\n",
    "nbins = 100\n",
    "bins = np.linspace(min(signed_residual),max(signed_residual),nbins)\n",
    "\n",
    "counts_0, bin_edges_0 = np.histogram((signed_residual_0), bins=bins, normed=True)\n",
    "cdf_0 = np.cumsum(counts_0)\n",
    "\n",
    "counts_1, bin_edges_1 = np.histogram((signed_residual_1), bins=bins, normed=True)\n",
    "cdf_1 = np.cumsum(counts_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell display the two CDF functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAElCAYAAAAcHW5vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXt0XNWZ4PvbVaWn9bZkW35bxjYmYEAyCZAE0yDnAabTDz9ImGR65oLFcNPpGUjjkJs12FlZAdNt39WZTia2M3clZBIeVuhJY0iCFGITAgQsAcY8bCPZlm3JT6lkvVWPff84p6SqUj1OSfXW91vrrKqz9z7nfKf2rvOdvb+9v09prREEQRCEWLClWgBBEAQh8xDlIQiCIMSMKA9BEAQhZkR5CIIgCDEjykMQBEGIGVEegiAIQsyI8hAEQZgESqntSqnNqZYjVShZ5yEIghAbSqmHgUagBijTWjemWKSkI8pDEAQhRpRSZVprZ/D36YQoD0EQBCFmxOaR5iilapRSe5VSbUopbX7u8svfbqb3+G0tocZi/cqG2mqmKOd6PxnXx3hsm1Jq+1SuP1mUUg8rpZqSdK2I9xlPWcK1gSmes95si9psSw+bn3tjrfMYr1ujlNplXrdJKVUblF9m3m/c7zkdrp+uOFItgBAe8w+5F9hibt3AamBtUNFWrXWd33G1wF6l1AattX/ZsuCycZKzBtgD1Gmt25VSZTGeYgvQHk+Z0pSMvk+tdbPZtsq01lv888yXFhIx9q+1bgcalFL1GO23NSjfqZRqDpYpW66frkjPI00xH8B7gQat9RNa63attVNrHbWRaq1btdZLgdVJehOqBw6afzJiHf/VWjcG/yGzkSy5z7VAqN5RN3BDgq+9CwjVo344SQ/uVF8/rRDlkb48ArRrrXdP4RxbgJQMBwlZSz3Q7J9gvpFXYDxcE8luoMy8nu/am7XWTyT4uuly/bRClEf6Uo8xFXAqNGM09qnaM8rMMd8ec9z+Yb+87RgPjXozP2xPx+8c2hwfrjHTA8aKzTHmJrNskzmm3uO7rpm22a9Mi/89mvLu9ZO3NujcLb7xa2CphftvMm0625VSPRauEct9hpUlRPkAm4j5G7T5XSdgLN7Kbx8L5vmdvt6TKf9moAFzyDLWc8aC2aNtNq/nG9Z9NpHXTKfrpxti80hfaoC2qZzAtD/4zuX7Y9cqpYKn2D0Rpdv9e6BZa13uG05TSi3VWjdorbcopdqAYPtKAOYfbbXWutzcr8UY6gjFXmCX1nqtqZxqfceZVGD0qOrMc/weo5fV4Hd8uylvvZnvO74JaAVu97u3cHL4X28PxoPCZy8KeY0Y73MysvjTjfHQdpoPcf/7HCNGmSJRD7SrceP4DRj2jw3RDlR+kzyisCXKsOdeYJcpQ3sKpsim+vppgyiP9KUdC2/FkfB7u/R/I4zJYG7+ScYMpOaDagPQo5SK9kcPpsZ80B6MMvZfy/gb3TOEGGcGnvW96SqlnsGcRGA+GFf7FJlp5G32G2qoMe1BBB8bhYNaa98bZ9RrRLtPM3+ysmBet9Hv+26zd1Eb5ppWf/tIrMVQ6r7rNpq9r4ejDd34frs48CxGT7ciRQvzUn39tEGUR/pyEONNbyrUwthskcni32vBPJ/Tr0dj6UGktW40ldkujAdZM0ZvJZTyaQU2YowxTxhjN2kJc6nVGEN1wb02nyKd7IPT30gc9hrmQ9zKfVr+7cJh9gK3m+eqCFcuxt8+EvWM9+58lDHFl5wYqQFDWcZykPkStClKsW4LSm5S189GRHmkL1sw3u7XT+ENZ7t5nqnQjvGwHCNMjyYq5tvpE+ZD7/cYPYpQb6zdGFMjt2Mo0ajDIkHHhuxd+d72Y5HZD/8HbdhrgOX7bJ+ELDN9X8xz92AogUYzLeyK3xh++5CYv50zxIvIakLPvgo+Pl7DVuFeJiJi/kbx6ClM6vrZiBjM0xTzD7QBw77wsPmn91+wFHYthVKqVinVgjEmO6WZIOafrtt3TfO6u4Ddsby5KmOBmX9PKtKYew1wH4Z9IaY3ZFPeMhVo1F9vDuc0m/u7zM8aJr5JT+kaVu/ToiztmDYWM99/IV6FXxnfm3VIIslktqewhnY/1mIo8mDKrMhg2sesbNHqei1T7LFNkVRfP20Q5ZHGmA+ppRiGyePmm2UT0Bb0J6tVfivMMYy7z0QyYMdIHcbD6ri5tU5yDHuLKd9xjLfYcIrNiTEs1YbR+9IqthXodcANfr/HWsZ7SXUY6196MHpmjUxu4V6ka1i9z2iyPANsNIfHtmAobSeMDUU2Ai1m/g3mseEevuFkWo/RXkJivjA8jKkk1MTZdFvM32EzCXyomop5O8abf02QMkw4qb5+OiK+rYS0wnwIBczcMt+MWzBmFslbXxwxf9tHrMyYEgR/pOchpCMVauI6BCfh36qFybOaxC/uE7IQ6XkIaYffwjOfAjkIbPfZCQRBSD2iPARBEISYkWErQRAEIWaydp1HZWWlXrx4carFyD7Odxqfs+aGzD5zwQ3AvKrMalrd3cbs1YqKsGvtphWW6zFKexAyj5aWlota66po5TLrHx4Dixcv5uDBUNPSBUEQhHAopU5aKSfDVoIgCELMiPIQBEEQYkaUhxAb937B2MJw2wMd3PZARxIFig/btm1j27ZtqRYjbbBcj1Hag5C9iPIQBEEQYibp6zz8HLzNtBL31/Sr047pCM6qK+TVq1drMZgLgiDEhlKqRWu9Olq5pPY8TGditRjO/sJ6hfUrvx3DM2yjqTSWRvLcKQiCICSHpE7V9XNDfQMWlAewOah30oThxXNaR/ASBGFqaK0ZcXsZcXsJNfqiNbg8XvqHRhgcGsY1MoTHNYx2j+L1ukKeU3m92N392F0D2N2DwPh5vV6Nx6vxWhzp8WqNVxPTMTavC7t3FJt3FLt3lLzqlVx9852Wjp0MabvOI0yMgW6mHl1PmAo/eNT4/EZo4/K3f3QegO8/MCtZEsWFp556CoAvf/nLKZYkPbBcj1HaQzLQWnN5yM2lgRF6Bl1cHnLRO+Sib9iFy2M8tN1ejUaD1uR4BvG6hrF5RrF7Xdj1+Geeu48Cl5MCt5MczxA27YlJlthKR8dGbMND2pTBA+jCxP4H01Z5YNg4ggPpRPSqajrU2wywcOHCBIk1zTn054jZbxweTpIg8eXo0aOpFiGtsFyPUdrDZHB7vDiHXHQPjDIw4mbY5WXE7Rn7HHF5GfF48Xi8uL2aYZcHl0ejtGfswZ/vclLi7kfhBcCmveS7e8l39WLXoXsO2YbNO5rQ86ez8gg1rNUNRoCaUBHHTLvIbjAM5okVb5ry9a0Rs793f2Vy5Igzd999d6pFSCss12OU9hCOEbeH7oFRLvWP0j0wSv+Im4ERN/0jbi4PucMO1Tg8wxS4neS6+8nRLvK9o8z09FM0coFC16WYewrZjC3M8Fq8SGfl4WQ81KYPcTyUaq67MWL2zasKkyRIfFmxYkWqRUgrLNdjlPbg42L/CCcuDnDu8gjn+4ZxDoZ4sGlNnqePElcPhS6j95DjHcLuHcWu3eS5+3F449eztdsUDptCqdD5NqXIcdix5+TjyM1DOXJR9jyw50DIYxReRyGenBl4c4rQNrvfuWzYbcY1w1wu8ExKYbcp7MqQL5yMgQLnou05aHse2p5LfkGRhYMmTzorj24m9j7KYCy+tyAIaYpv6OnkpUE+OnuZ85dHJpRxeIYoGj1P0cgF43P0Ag7vxHJWcNgUhbl2CnLs5OfYycuxkeewjykHmxpXEtqehz2/GFtOHthzjc1hfs+dAYUzjS2/1Eiz9OSefqSt8tBatyqlgpVEBSABgVLJgReNzzV3hMze92o/AOs+k9i3nnjT0tICQF1dXYolSQ8s1+OBF/FqzaW6ejqdQ3Q6hzjfN4Jz0DVh6ElpDyXDnVQMnaBs+BR57v6wp81z2CjItZPvsOGw2bDbjV6Cw6aw22w4bAqbqRjsSuGwKxQK8opgxixTAZSDLce8uIK8YiM9t0gUQhxIK+VhLiCs1Vr7puI+q5Ra77e/FgmZmVp+/gPjM4zy2PlLY45DpimPffv2AaI8fESrx2GXh47uQZb//AfYgP/tWjaxkNbku3spHe6kdOQMpcNnsJtGXAXk59gpyLVTmGsnz2Ej12Ej1258Omxh5hjZHIZSKCgHRz7Y84xeQ9FsKJ4D+SVxuHvBCklVHub0200YK8x9iwCf0Vq3mkXWYyiIRgCtdYNS6mFzYWAN0OanSIRUcMsXI2bf+ekZSRIkvtTWhpoZPn0JV48X+kZ455STj7ou4/ZqRmpuDsh3eIYpGz5FyXAnpcNnyPMYvQsFlBbkUFZYSFGeg8I8B/ZQb/+OXEMRFFaODx05cg0lkZMPeaUQTrEISSVrw9CKexJBmBper8Y55OJ83zDnL4/Q1TtEpzO0wXrGyHnm9L9P5WAbShvTY31DT6UFOVQW5ZFrD3ro2+wwowpK5kJxtfFZOFOGlFKMVfckaTVsJQhC6vB6NccvDXDsXB8X+kfpGRjF4w3xcqm9lA91UDx6zlhX4eqhyNtHcb6D4vJ8ivMdFObaQw895RbCzGVQuRzKFxkzl4SMxL5169ZUy5AQtm7durW/v59bb711LO2pp57iueeeo7q6mspKYx57S0sLe/bsoa+vb2y6Zl9fH48//jgtLS3cfPN4t3z37t3s27eP5cuXU1xcDMD+/fv52c9+BhjRCwE6OzvZuXMnR48eDRhD37FjB01NTdTV1ZGXlwfA888/z9NPP01xcTFz5xqhPI8cOcIPf/hDurq6uOaaa8aO37ZtGwcOHEjtPX34AXXLaiC/MOQ9XXS6eeGF52nc+0zm3NPRoyxfvpzR0VHy8vKyo55iaHv23AI6Rwt46YOztB76gBN/fA7nhfN4ihbhMJ/tl197ipFTh5kzs5jll5qp7n+f9s7LHD03wnWlbmoXFVNVnM+xHhtPvudlYBRWzLSBzUFf3lwef/kiLRfzuHnjf4Wq5VBYwe6f/C/5P6XhPe3fv79r69atUR3QyuChEBunj8M37wmbvfHbnTS/OZhEgeLDzp072blzZ6rFSBp9wy7OXzaGoF5ru8Qb7ZfoG3aP5V+8YKfpt/nGjt/Q9oqLLzHH3seK2UWUFBiapeDgx6jgoaai2XDdl+Ez/w2u/isjTdnEXpFFiM1DiI2HvmJ87vhlyOwNj5wBYO9j85IlUVzYsWMHAA899FCKJUksvUMuDp7o5v3Oy6GHpExe+o3x1vuf6z/iiqH3qMDJjFwHRfkO8h3ji9947nXj829vhvIlMH81VNSI3SKDsWrzEOUhCNOAwVE3b7Rf4vCZyEojx65YMaeEq8rcVJ1pJre/M3RBpWDWVYai8C2qc+QmSHohmYjBXBAE3B4vb59y8ubxbkbd3pBllIL55YWsmF3M8tkzyDt/CD5+GTxhfCNVLoMla6CoKoGSC+mOKA9ByFJOXBzg5Y/O0zsUWgmUFORQt6ic5bOLKMx1wPBl+OBX0N0+sbBSULUCFt5kLMYTpj2iPITY+O7Xjc///q8hsxse6wJg1yPVyZIoLuzebUwu2bx5c4olmToDI24OHL3AkbN9IfNLC3L45JIKVlaXYLcpwyB+9jAcewnchm+phqYvALBr7W9h5lK4oh4KQ/gljdIehOxFlIcQGx0fR8w+diozYyV0dXWlWoQp0zfs4p1TTg6d7g05RJWfY+dTNRVcO7/MUBoAg93Q9jJcPBZQ9phzpvFlxReh+trwBvAo7UHIXsRgLsTGCfMhsziELyPgaIfhu2j5wswynnZ2GoZh3zz6TGLY5eHVYxf5oCu0MVwpuHZBGTfVzCQ/x5wp5R6Bk6/B6bfAOzEGxlHXSlj8GZYvjxLXI0p7EDIPMZgLiSHKQyLTlIaPTFQaAGd7h3nhvS4uh7FrzCrJo37lbGaX5I8ndh+HIy8aNo5gbA6oWcPy+TdYm24rSmPaMu2UR29vLxcvXmR0NLEhGoXkkpubS2VlJaWlpakWJSlorXn3dC+vHL0QsrdRmGtn9eIKrl9Qhs03ROVxQdsf4ExL6JOWzoMVd8CMzIwGKSSXaaU8hoeHOXfuHPPnz6egoGDiqlghOs5LxmfZzJDZF53GKuXKsuQ1La01Q0NDnD59mry8PPLz86MfFMT+/fsBAtw6pCsjbg/NH5zn6LmJBvFScwbVVXNLyPF3RDh8GQ49AwMXJ54wrwhq/gJmf2Kst/HTfUYonb9bFyoatB+//rnx+aWvTupehMxlWimPCxcuUFVVRWFhZoZKTQuiKI9Lvcb4eTKVh1KKwsJCKisruXDhAgsWLIj5HAcOHADSX3lc6BvhhUOd9IQI43r9wjI+u6xq3BjuY/gyvPNLGOoJTFcK5t8Aiz87YYHfky8aQ1pRlcfzvzA+RXlMO6aV8hgeHmbOHJmjPiXCKA0fM0vtEfMTSXFxMZcuXZrUsWvWrImzNPHn2Lk+fnv4LO6gYapch43PXTWbZbOLJx4UTnEUlMGV66AstKL92h0WgyrdFd7PmZDdTCvl4Xa7cTim1S3HnyjKI5k9jmAcDgdutzt6wRCke4/jxMUBXnzv7ITQrpXFeay7ppryGSEmKgw54d2nJyqOymWw8i8juhOJ2uPwIT2Oacu0e5KKnSN7yda67eodYt+hzgmK4+p5pdy6oirQtuGj75xh4xgdCEyvWg5X/ZURiEkQpsC0Ux7CFBkxI8nlhTZKD48Yi9Py8zLL9Xa6rvO41D/C/3m7E5cnUHHUr5zNNfPDzCzrOQGHfwXuoBmFMSgOy+t1ZJ3HtEWUhxAbXR3G5+LlIbNPnjUMuSsW5SVLoriwZ88eAB599NEUSzJO98Aoz7WeYdgVuIjv1hVV4RVHdzu81zhx4d+slbDyLss9jvsfPwvAyz9aGLng9/7e+PzJby2dV8geRHkIsZEbeRpsfm5mDh1VV6eXL67ugVF+1XKa/pFAG86naiq4fmF56IOGeuCDX09UHPNvgCtujynGxrIFFsPDLrzC8jmF7EKUhxAbcyO/iS6qNoY56urqaGhoCOtocOnSpaxfv57t27fHXcTJkE4OES/1j/Cr1tMMjAQqgVXzS7mpJsyEBY8LDj8HruHA9KW3wcJPxSyDZceW4hBx2iLKQ0gJ27dvp6amJtVipB29Qy6eaz0zQXFcPa+U266cFXpSgNZw5DfQfz4wvebWSSkOQbDCtFYe/2/T0VSLEMB/WxvajpCNrF+/PmB/9+7d7N27l6amphRJlHqGXR7+z9tnJgxVXTOvlNtXRlAcJ16Fc+8HpletgIU3JlBaYbozrZWHMAlOmYGCFoTuNbSdzkyfYamOYe7yePn3dzrpHgj8/VbNj9LjaHsZTr0ZmD6jEq68c0pxxC3Hoo8S017IXjJrPqUwxtq1a9m9ezdr166lvLycuro62tvHI8A5nU4aGhooLy9n6dKlPPHEEzEdH8yGDRtoaGgAjxs8bp544gnKywMNt0op3n23Fbc5rdTpdLJhw4aQ56+rqxsLwOQ7d3Nz85i8/vfhO8fSpUtpbW2d/I8Wgf7+fvr7+xNy7mhorfnd+2c54xwKSF8+uziy4jj6u4mKw5ELn/gbcExtttulXs+Yq5mI9HYbmzDtEOWRoXR3d7NlyxZ27drF8ePHAQKMz7fffjtlZWX09PTQ0tJCU1OT8fC3eHwwmzZtorm5GebXwPwampqaWL169djDvLm5mbKyMu76/KeomWcYzR977DG2b98e9fx79+5l165d1NfX09PTQ1tb21jehg0bqKiooKenh127dnH77bdP5ueKyoMPPsiDDz6YkHNH483j3Rw7F6i45pcX8PlPzI5g43gROt8OTHfkwapNMCOyFwArPPv9uTz7fQtrXv75F8YmTDum9bBVptsYNm7cOGZ03rRp05i9oLGxEafTOfawLisrY+/evZSXl7N9+3bKysoiHh+K9evXs2HDBto7OqipqaG9vZ2GhgaeeeYZamtraWpqYuPGjeQ4xh92sZw/FK2trRw8eHDsuPr6eurr62lubqa+vj6mc0WjuDiEX6gkcPLSAK+3B/rjqizK5a5r5+IItXJca/i4GboOBabnFMC1d8ctvrhlNzNR3NUI2Yv0PDKYurq6kOnt7e0TZjL5FEbw0FG485aXl49tzc3NgPHwbmxspLGxcexB3tjYCBgKa8OGDZbks8rBgwdxOp0sXbp0bGttbY04vJZJXB528ZvDZ/H3OlKQa+dL188bj/gXzPEDcDooQmZeEVz/H+KmOATBCinpeSilHgbagQoArfVuC+Wd5m6Z1vqJSOWnOzU1NQSH4PU9cK1Mj21pCR0saMOGDez95S+oXXUNa9eupba2lu7ubtrb22lvb6e+vp6zlybnmDAUFRUV1NbWhpUnnjz//PMA3HXXXQm/FsCo28sLh7oYGh23KygFX7x6DiX5YRbodfwZTr4emJY7A667Bwor4irfjl8YvaGH7onSs3jyX4zPr/1DXK8vpD9J73kopbYD7VrrRlNpLFVKrY9Q/mGt9RNa691m+WZTmQhhWL9+PRUVFTQ0NOB0OseM55s3bx7rgUyGjRs3cvDtd2j+w/6xqbYbN25ky5YtY/u9/R56+y0YWoOoqKgYU3C+ns769etxOp0Bxv7GxsaEGM1bW1sTZowPZsRtTMk92xu4oO+mmpksmjkj9EHd7dD+h8C0nHxjqCrOigPghT8N8MKfBqIXfOU3xiZMO1IxbLVZa93ot98ENIQrDGzy39FatwI3JEKwbKKlpYXu7m6WLFnCkiVLqK2tZdeuXVM6Z1lZGTU1S6ioqhpLW7t2LY2NjWzaZFTT7AoHsyti79D6bBjl5eUBcra0tPDWW2+NDaE1NTUlZHHhunXrWLduXdzPG8ywy8NzrWcmzKyqqZrBJ5eEUQI+tyP+41uOXMM4XjQrIXI++JUKHvyKBaX01W8YmzDtUFpPjH+csIspVQv8XmtdHpTWorUOOSldKdUEOLXWG8z9zUB3kAKawOrVq3Xw0M2HH37IypUrp3gXQjqTznU87PLQ2HKaC30jAemzSvL429r5oe0cHhe0Phm4elwpWLURKmSFvhB/lFItWuvV0colu+dRAQRPCneGKuhHA1CrlOoxh6vCKg6l1Gal1EGl1MELFy7EQVxBiB+vHL0wQXFUl+aHVxwAR38b2u2IKA4hxSRbeYQacO8GUEqFHIzXWrcDu4CDwHYiDFmZdpHVWuvVVX5DK0IcGew3tjD0D3roH4zd5pFqjhw5wpEjRxJ2/hMXB3i/83JA2vzyAv66NsLMqnPvw9nDgWlVK2BB4v1VvXZokNcODUYv+M4bxiZMO5KtPJyYM6z8iDiwqpTaBTRrrdcCa4HNSqm9CZJPiMb5TmMLw5kLbs5ciN+Mq2Tx9NNP8/TTTyfk3CNuD80fngtIm1mUy19dP488RxjFMdxrrCD3Jw5uR6zynR9f5Ds/vhi94L9uNTZh2pHsqbrdTOx9lAForScMX5n2EKdpJEdr3ayUWgIcT7SgQhgKwswGMplRkJlLh5YvT9yC0T99fJG+4XGFqhSsvWp26PCxAF4vfLgP3H5DXDaHEQVwim5HrHLj1ZHjtoyxSrz2TleSqjy01q1KqWAlUQE0hzmkAghYfqu1diqlwpUXEs3syI7y5s+yGEQozfjyl7+ckPOe7hnk3VO9AWm1C8upLi2IcNCb4OwITKu5FYqSNxT7/QcszuL6xrbECiKkLal4TXw2aF3HWgybBgBKqRpfvta62czHL78MY4GhIKQ9r7cFuh4pL8zhpqURFt45T0H7gcC0iiUwP+rkF0FIKklfYa61blBKPWwqiBqgLWj21HoMheFLazAXFrb5nWNL0gQWhEnS1TvE6Z7A9Rz1kYarRvrg/X8D7R1Py8lPmp1DEGIhJe5JIrkXMfOe8NtvB0RZpAsnzABai0PbCI6cNMbpVyxKzth8vNi2zRh+efTRR+N2zoMnegL2F1YUMr+8MHRhr8dQHKNBq7pX3Al5yXfaeNsDxrDZyz+KHHaYe79gfP7ktwmWSEg3MtO6KQhpTvfAKG0XAqc0r15cHqY08PHvofdMYNqim6Aqsz0/C9nLtHbJLkyCMD0OH5nW4/ARzx4HQMvJngBvIrNK8lhYEabX0X0czgQ5f6yogcW3xFWmWIja4/AhPY5pi/Q8BCHO9I+4+bArcEHg6kUVoQM7eb1GfA5/Csrgqr8Em/w9hfRFWuc0wT/sayiWLl3Kli1iWooHLSd78HjHux2lBTksm1UUunDX2zAQtBhv5V8awZ0EIY2RYSsBMELEWvJWe84clw+z3uP0eReQees9nnrqKWDq6z16B128eypwKVPdonJsthC9DtcQHP9jYNqcq6E08lqaZPDtHxn+tKKu9/iBOdwn6z2mHdNbefzhsVRLEMhfPJKyS/vicfjYvXs3e/funRg6dihyjIeBIW/E/HTl6NGjcTnPqx9fDOh1FOc7uGpuSejCJ141FIgPe46xGDANeOPwcPRCAIf+nFhBhLRleisPIXZmzY2YPa8qM5vU3XffPeVzdDqHOHquLyDt5qWVodd19J+HM0HBpxbdnJJpuaH43v2V1gp+fWtC5RDSF7F5ZCgbNmygoWE8htYTTzxBeXngVFClVEB0PKfTyYYNGygvL6eurm5CPHOfTcR37ubmZsrLy1m6dOn4OUbdbPiP/2ksPTj6XlGhnaLCMM7+0pgVK1awYsWKSR+vteaVo4FhAGaV5LGyOoQy8Hrgo32BiwHzS2H+Jyd9/Xhz86pCbl4VZnaYP9fdaGzCtEOUR4ayadOmsXCtAE1NTaxevXrsYd7c3ExZWRm1tbVjZR577DG2b9/O8eOGX8nt27eHPPfevXvZtWsX9fX19PT00NY2trifDRs2UFFRQU9PD7t27eL2229PxO1lHEfP9dMVFFb2lmVVoWdYnfwT9AV62eWK28Gemb02YXoyvVtrCm0MU2X9+vVs2LCB9vZ2ampqaG9vp6GhgWeeeYba2lqamprYuHFjwDEbN24cM4pv2rRpoj0jCq2trRw8+BZNzxmeY+rr66mvr6e5uXksjKyzz4jlUVacWb2PlhZjnUVdXV3Mx2qt+fPxQB9WS2cVsSBZ5y0UAAAgAElEQVTUuo7LXXDy9cC0WSuNOB1pxL5XjQWO6z4TZpaYjwMvGp9r7kiwREK6Mb2VR4ZTX19PY2MjNTU1Yw/yDRs2sH37dhobGyfELJ/Mg9GfgwcP4nT2svTqa8AxPpvKf/jrXLfhejzTlMe+ffuAyf1Gp7qHuNQ/OrZvU4rPXhHCZuBxTxyuyp0Byz4X8zUTzc5fGgE/oyqPn//A+BTlMe0Q5ZHBbNiwgb1791JbW8vatWupra2lu7ub9vZ22tvbx3oD8aKiooLaVato+f1LUDk7ZJnSosxSGj78h/di5Z3TgVNzr5hVRPmM3IkFO16fuKZjxR2Qa8G2kGTu/HTkuC1j3PLFxAoipC1i88hgNm7cyMGDB2lubh6bartx40a2bNkyYeptrFRUVIz1KHy2lfXr1+Ps7+eJ/+9nY+UaGxsDjOZzZjqYMzPz3knuuusu7rrrrpiP6x1y0R7kw+raBaUTCw71QEdQuNbqVVB5RczXTAYP3TOTh+6J4Drex9f+wdiEaUfMykMpdZ1S6m+UUt9USt1rfr8uEcIJkSkrK6OmpoaKivFIvmvXrqWxsZFNmzZN6dy+Xkt5eXnA8FdLSwtvvfUW5eXllJeX09TUZG1xYZZy6LQzwIdVVXEe88pCrA7/+Pfg9QvPm1sIS2WygZC5KO3f8sMVMkK/bgGWYARiaseIRw5GGNmlmLE5gO1a6xOJEDYWVq9erQ8ePBiQ9uGHH7Jy5coUSZQluM0HoCN078LlNtpTjiM18ScmW8d9fcb6jOJi6+ssXB4vP/njcYZdnrG0tVfN5up5QT2PS21w6NnAtCvvNHoeacpFp1HPlWVRepFOc6JAmYVeipARKKVatNZRo49FHV9QSt0HlAJbtNa9UcqWApuVUj1a659YllbIHE6bxvEw3nXbzxiG40zzrrtz504gNu+6R872BSiO/Bw7K+YEKR+Pe6Ljw9J5MOeaScuaDDZ+uxOw4F33m/cYn+Jdd9oRUXmYiqNZa33cyslM5fJPSqklSql7RYFkIVHWIjjsmRnxrqgoyqyiIFweLwdPdAekXT2vZOJq8lN/hkG/ckoZs6vSPDLgzFKLEx9KK6KXEbKSiE8CrfWeyZzUVDaiOLKRBZHtG0vnh5hllAE89NBDMZV/ve0SPYOusX2lYNX8ssBC/ReMBYH+VF8HxXMmK2bS2PuYReeMO36ZWEGEtEVmWwlCjJxxDtHaERhi9uq5pZQW+HkS9nrhyIuGKxIfuYWwJHUBngQhnsRNeSil0sxFbWisTBAQMpNk1K3L4+Wl988GzLAqznfw2eVBiwLPHITLnYFpyz6Xlms6BGEyWJqQr5QqASowZlTVMD7DqsJvvwZIa38fOTk5DA0NUVgof+BJ09lhfM4NbUg92WUYzBdVJ3/4amhoiJycycUR8TmF3Lx5c8Ryf/r4Ik6/4SqAz101hzyHn41gsBvaDwQeWLUcqq6clGypoOGxLgB2PVIdueB3v258/vd/TbBEQroRzWC+BGjBmG0FoDCm6LYD3RjK4xmgF5jaqrQkMGvWLM6cOcO8efMoKCgI7bROiMxo5DgPw6PJ79lprRkaGuLMmTPMnh165Xs0urq6opZxDo7y7qnACYfXLihl4Uy/lxGt4ehvA9d0OPIywkjuz7FTruiFADo+TqwgQtoSzWB+XClVAzwONGmtf+Wfr5T6R631P5nf20KdI50oKTGC8nR2duJyWfxzCIG4zd/tww9DZvvWeXz4YXIflDk5OcyePXusjmPlvvvui1rmjfZLeP3Gq0oKcvjMFVWBhc6+Bz0nA9OuqE+bOB1W+fG3LBr1v/M/EiuIkLZEHbbSWjuB+5VS1yulvgns1lpf9mX7lft9gmSMKyUlJZN+wAjZy9y5kYNcXeof4aOzgYGePn3FTHIdfmbD0QFoC/obVCxJ+zUdoVi+0OKw4+JliRVESFssOyHSWr8NvK2Uuk8ppWUNhzCdeKO9O8BIXlmUy4rZQb2Jj5vB5TesZ3fA8s9n1HCVIFglZg92Wus9SqlSpdQ/Ytg8hOnEr39ufH7pqyGzf7rP8Frzd+vKQuanK/v37wfg1ltvnZB3oW9kQnjZG2tmBtrMLrXBuQ8CD1z8WSgIjO6YKViuxyjtQcheJjVVV2vda9o6njUdJJYopf4mzrIJ6cjzvzC2MDz54mWefPFy2Px05cCBAxw4cCBk3uvtgYGeqorzuGKW34p0raHt5cCDimalVVjZWLFcj1Hag5C9TMl3tt9Q1j9iGNUt+TRQSj2MMWOrwjzP7ijlyzCmAfuM8ge11q0RDhESxV33RMz+2h2ZaU9as2ZNyPTzl4dpOx/ocv2mpUG9jgsfBcbpUApWfBFsmbsG13I9RmkPQvZiyauupRMp9azWeqOFctuBt7TWjaH2Q5QvA/Zqrdea+5uBtVrrDZGuE8qrriDEyr+/2xmgPOaU5nP3DQvGlYfW8NZPApXHnKthZeyxQQQhHbDqVTeer0ZbLJbbHKQomoCGCOX3AP7xVJ+N4VqCMGnO903sdUywdYTqdSz6dJIkFITUETflYcXzrlIqVKzPbiBSvNT1QKtSqkYpVau1dmqt2yOUFxLJiWPGFoajHaMc7RgNm5+udHZ20tkZ6E7kzeOBXnNnl+SzOHhB4IlXA080+xNQmPnzSCzXY5T2IGQv0VaYT2VG1SWt9T8HpVVgKAt/nITBT9nUAq1m2l7gPnP9SXD5zcBmgIULo8QhECbH9/7e+AwTv+H+x88CFuJApBl79hgOpH3xPC70jXDsXHCvo2La9Dos12OU9iBkL9EM5hEN2ZMg1Ly/bjBsGyEUgs//d7uvt6GUegZjKGuCzcM0vO8Gw+YRL6EFPxZGjrm9bMHkfEulmurqQB9OoXodSypnjCe4Ryf6r8qSXgfEUI9R2oOQvcTNYG7pYkrVYxi/y/3SfOFry4OVh9nzaNFaq0hpoRCDuTBZegZG+dnrJwIWBf7ldXNZWuU3PffoS3CmZXxfKfjk5qxRHsL0JWEGc3OF+f+cnFh0M7H3UQZjblCCCZtmzsIShLjzfuflAMVRVZxHjX+vo/t4oOIAmFsrikOYVkzWYD4pg7W5NiNYIVQAzSGKYw5VOc3eiY8ywBlG2QjClPB6NR90BXrOvX5h2bitwzVsBHnyp6Acam5NinyCkC5MdoX5P03hms8qpfzdt6/FbyquOavKP/8xAmdjbTLThFTw0FeMLQwbHjnDhkfOJFGg+LBjxw527NjB8UsDDIyMR//LddhYNsvPh9XHzTDst/JaKbjyTnBkZvjdcFiuxyjtQcheJrPC/Fml1O+Ax7XWf4j1YK11g1LqYVNB1ABtQes+1mMolEaz/BNm+YfN/Eta6ycmIbcQD3qDJ8sFcqnXEzE/XenvN2ZWHT4T2OtYMbt43HNu1yHD5bo/82+AsgXJEDGpWK7HKO1ByF5iNpgrpX6MMfRUC9QBBzEW+jVqrU/EW8DJIgbzBOE0/TyVzQyZfdFpBEGqLJuS55uk09fXx8CIm2feuRgQs+PLn1zInNJ86L8ArT8Fj1+QpxmVUPefDO+5WYbleozSHoTMw6rBfDKtvgkjMNRl80K3Mz709PlJnE/IJKI8JDJNafgoLi7mo0vdAYqjsiiX2SV5xrTc9/8tUHHYHHDVl7JScUAM9ShKY9oyGZfsvzIDQ2mt9TtmEKiMCAQlCOHQWvN+0JDVJ+aVGobyY7+DwUDPuiz/nOE5VxCmKZN6bTK96QrTkSf/xfj82j+EzN7xC+Mh+9A9mfVG+syv/g+dFwcoWGq4UbfbFCvnlEDfOTh7OLDwnGug+toUSJk8LNdjlPYgZC+Z6zNaSA2v/MbYwvDCnwZ44U8DSRQoPhx5/xCuc21j+0uriijItcOZILvZjEojOmCWY7keo7QHIXvJzgFbIXF89RsRsx/8SuYtlHN7vBQt+yRuz7i946q5JUZM8uDogEvWgD0zXbDEguV6jNIehOxFlIcQG2vuiJi97jNFEfPTkROXBrBVLcW3UqMw187CikLoeA28fkby/FKYOT18OVmuxyjtQcheZNhKmPZ82BUYn3z5nGLseKEzKFjlvLqMjg4oCPFE/glCbLzzhrGF4bVDg7x2aDCJAk2NYZeH4xcHcHWfwdVtrKheOacELhyBET+X7HZH1hvJ/bFcj1Hag5C9WFIe0Rwh+ucrpcR1SDbzr1uNLQzf+fFFvvPji2Hz042j5/rweDVDH73C0EevUF6YY6ztOP1WYME5qyAnPzVCpgDL9RilPQjZi9WeR7tS6rpQGWa6v6PE8lDlks3RjlFue6AjIO3bPzrPbQ90BLxR7Xu1n9se6BibmgjG6trbHuiY4Nun4bEubnugIyDC2k/3ObntgQ5+um/cT6Pv2g2PdQUcv+GRM9z2QMfY6l0wpkTe9kAH+14df8t97dAgtz3Qwbd/dD7g+Nse6Ej5PR0qvBZWfSrsPd14dT6zyu0Zc0//9fuXeeUPuTjK5+Ion8uV1SVs/NYpbvtff8HFoYKxsjtevyZj7ikeba+s2MaNV48ry3D39LpaRffi8cXI6XxP2VhPibgnq1hVHruBR8LkNRAYY1yCMGUxP6z+r/CNbWHzv//ALG64KjPe0PuHx/+ghSvXULhyDStnF4NnJLBgRQ048pIsXWq5clEu338g+iLI/yfn63z0l99OgkRCumHZt5VS6nHgx/7+q5RSS4AGrfW3/NL+p9b6v8Rb0FgR31ZCNF77+CJ/9osYOLcsn01LRuG9vYEFr/8PWen8UBBCkYhgUI8BW4LSHga+H4tggpAOuDxeDgW5I1k5pxiO7w8sOPMKURyCEALL6zy01r3KoERrfVkpVYrRc7kc9WAhe7j3C8bnT34bMts3ZvryjxYmS6JJceRsH0Oj427HL7/2FP/2Gqy6JSgux5JbkixZemC5HqO0ByF7iXWq7nbGbR/fAh6PrziCkHi01rzd0RO94OyroHh24gUShAwkphXmWuvjSqmlSqkSYGk6xe8QkkSUN8x073EAnOoe4mL/+GwYm1I89Hd3UdT+u/FCygaLP5sC6dIDy/UoPY5py2QWCT4OtAA/DpOvJi+OICSet08F9jquqJpB0fmg1eTV10Jh5vnpEoRkEbPy0Fq3Ak9orV8OU6RpaiIJQuLoGRil/UKgt9jVpb1GpEAfSsHCTyEIQngm5Z5Ea70nQt6vJi+OkPb84FFjC8O3f3R+wsKldCK41zGnNJ9Zznd46rCLpw67jMSqFVCQFmtdU4bleozSHoTsRbzqCrFx6M8Rs984PJwkQWJncNTN+2cCJweurnShTp7gaLffeqf5n0yyZOmH5XqM0h6E7CV7HSOePDY+jdDHDx410vwduR140UjzRUQDcF4y0h76SuDx3/26kX7i2Hjar39upP365+NpJ8xrf/frgcc/9BUj3ekX0vTJfzHSDrw4nvbOG0Za8BvdvV9I/T1VVsPXt4a9p+/dX8n3VryelvdU+MA6Vh96YSxp8UAny77/f8NvWrj7Ew7u/oQDSufDd/8x8+tpim3ve/Ne4nv3V0a/J4A7v5wR95SN9ZSQe7KI9DyE2CicAdfdGDb75lWF8M7ZJApkDY/W2IPSVs4cT1kx03yPWvApQCLj3Vx8ClYVWiu8ZEVihRHSEsvuSTINcU8i+PPOKSd/+Gh8DD8vx8a9c0+Se+bN8UKFFfDJzYbBXBCmKYlwTyIIRnfXv/scxL5X+wO8f6YDXq+m9WSgofzaOYXknj80tt/S5aGlv0oUh4nleozSHoTsxdKwlVLqb4AbMBwjngyR/7dAiywanAb8/AfGZ5jwozt/aTgaTKdwtB9f6Kd3yDW2b7cprs87Da5xo/C+Yx449h51t/9NKkRMOyzXY5T2IGQvUZWHqRj2AN3AZqVUbbAC0Vr/Sil1L/CTxIgppA23fDFi9p2fnpEkQazzzilnwP6Vs4soPB+4TKn2imooqU6mWGmN5XqM0h6E7CWqzUMp9YzWepP5vR64T2u9SSn1TWAp0IYRDOoGrXW4mB9JR2weAkD3wCg/e+1EQNp/XOGmou3X4wk2O9z4XyCvOLnCCUIaYtXmYWXY6rjvi9a6WSm1xAw7OxMoAzYBpUx01y4IKee9ILfr88oLqOj+Y2ChqitFcQhCjFhRHpeC9psxnCJOOuCTUuphjN5KBYDWencMx+7SWjdM9trCFPHNPy+bGTLbFz6zsiz1s8DdHi8fdgUuCryufAQ6gsx2Cz5JX18fAMXFokQghnqM0h6E7MXKbKsA73Ba6+PAW5O9oFJqO9CutW40lcZSpdT6GI6tmey1hTjwzXuMLQwbv93Jxm93JlGg8LRdGAiI2ZGfY6emP8gBYtkCKJ7Dzp072blzZ5IlTF8s12OU9iBkL1ZeD7copWoxHB42aa3fZWpxyjdrrf2HuJowhrwaIx1kyiCkmtLInmZnlgYvxUsdwUNW1xf34uhuCyxkuiIpKkqf2WHpgOV6jNIehOzFisH8WYyexlqgHujBGHLaBTT7pucqpe7VWkecbWUqgN9rrcuD0lq01hEn2CulNpvX3aK1XhvlvsRgPs3pGRjlp/6Gcq25b8YBikb9RmFL58H1X5W1HYLgRzwXCe7SWv+T1vpzWmsbhhJ5FtgItCulLimlfgfUWThXBcaUX3+coQr6Yw5rPWuh3Gal1EGl1MELFy5EKy5kMYc7A3sdV9k7AhUHwNLbRHEIwiSJqjy01r8P2m8NUiabgLctXq8sRFo3gFIqVJ4v3am1jqpktNa7tdartdarq6qqLIokZBser+aDznFDuc3r5lrXu4GFZl1pOEEUBGFSTNk9ida6WWv9LYwhpWg4CTLAh9gPZqPWunlSwgnx57tfn+gJ1I+Gx7poeKwriQJNpP1CP4N+hvKFQ+9TlTsedhabHWpuDThm9+7d7N5tedJf1mO5HqO0ByF7ied8Siv/vG4m9j7KAEL1LEx7iCiOdKLj44jZx065IuYnA/8hK7tnmFXeD7GrnPEC8+omBHvq6kqtwks3LNdjlPYgZC9xUx5a614LZVqVUsFKooLwCqICqFfj49I3ADXmOpFGrbWV3o4QT77zPyJm//hbc5IkSGh6h1ycvDQ4tj+v713mlPjZNXLyYdGnJxx33333JUO8jMFyPUZpD0L2ElF5KKXuw5ieeyKWkyqllgC3h5l99axSar3W2jc1dy3GzC3fsTVArbkOpBk/xWLOuKrRWj8RizxCHFm8LGL28oW5SRIkNB90XsY3gTDHPcAVox8xI9dvGu7CmwwFEsTcuXOTJGFmYLkeo7QHIXuJaPMwY5WvVUo9ppRaHO1kSqkSpdTjwPXhpu2aq8NrlFLrzR5Em58iAVgPTFhBbioO37HbwxnYhemL16t532/Iav7lVuYU+b0f5RUZQ1aCIEyZqMNWWus9Zk/ifqXU9RiGcSeGQ8Qyxn1c+ZwkPmGuQo90zrA9BzNvQr65Gl0smqnGF0rzS18Nmf3Tfcao5N+tS75uP9k9SN+w4VYjz9VL9eARZlaWjhdY9Gmw54Q8dv/+/QDceuutCZYyM7Bcj1Hag5C9WLJ5mMrgWzA2JFVjbr0YjhPbtdZWp+sKmczzvzA+wzwsnnzRmCKbCuVx2G9F+YLeFmbOyMFhMzvXBeVQfW3YYw8cOACI8vBhuR6jtAche4nZYG4qkuPA76OVFbKQuyL7MfraHSVJEiSQy8Mu2i8MAFAw2k3lYBuz5/o5OVzyWWOKbhjWrFmTaBEzCsv1GKU9CNmL5RjmSqkSAK315Whl0wFxTzK9+NPHF3nzuOG8YPnFJua7O1g1vxSFghmVcMO9sppcECwQ1xjmSqlnMHxa9SilfjtV4QQhnrg93jEniIWjF6kYPM6cknxDcQAsWSOKQxDiTFTlYU7X3QAoc/ucUuqhRAsmpCknjhlbGI52jHK0YzRsfiI4cq5vzPX6gt4WHDZFZVGekVk8ByqjTyft7OykszM9XMmnA5brMUp7ELIXKz2PBqBBa20zfVltBO5OrFhC2vK9vze2MNz/+Fnuf/xs0sTRWo/FKC8aOUf50EmqivOw23y9jlss9Tr27NnDnj17EilqRmG5HqO0ByF7sTrbao/f90ZzzYUwHVl4RcTsZQtCT4VNFF29w5y/PALA/N4WFDCnxFwEWDoPKqzFDquurk6QhJmJ5XqM0h6E7MVKPI+XtNafC0r7sdb6/oRKNkXEYD49ePG9Lo6c7aNw9BKrzv6K8sIcrpxjzhS69m6oWJJaAQUhw4inwbwnRNoEjWMa1QUhafQNuzh2rh+A6r73AJhT6tfrKF+cIskEIfuxMmxVo5T6C0AFpd0bXC5+YglCdA6d7sWrNTnuASoHP6Ygx05pgTncsuBGmWElCAnEivKoI7TX2/o4yyJkAg99xfjc8cuQ2RseOQPA3sfmJVQMl8fLodPG9Nzq/sMo7WVO6Qxjem5hBcyMbSx+x44dADz0kEwkhBjqMUp7ELIXK8rDCVgxHojHuelAb3AU4UAu9Xoi5seLj7r6GHZ5sHtHmd3/IQ6boso3PXf+DWCLLc5Zf39/AqTMXCzXY5T2IGQvVpRHs9Z6Y7RCSqmoMcaFLOCffxEx+9nvJ961udaat08ZpriqgSOGAikrMKbn5hTAnGtiPueDDz4YbzEzGsv1GKU9CNmLFeWxK3qRmMoJmUzZzIjZlWXxDE4Zmo7uQS71j5Ln7mNe7zsoYHaJ2euYVxvWc24kiouLoxeaRliuxyjtQcheovbttdaWHCBaLScIU+XtDic2r4sVF14ixzvEzKJc8hx2sDkkXocgJInYBoYF4cl/MbYw7PjFJXb84lLCLn+hb4TjF/q54tIfKHQZ1xmbnrvoJsidManzPv/88zz//PPxEjPjsVyPUdqDkL2I8hBi45XfGFsYXvjTAC/8aSBhl3/rRDfze1uoGDoBQHGeg+K8HKhaHjI2uVVaW1tpbW2Nk5SZj+V6jNIehOwl8QPUQnbx1W9EzH7wKxUJu3TPwChnOtq4+vJ43LG5ZQVQNAuuvGtK6zrWrVsXDxGzBsv1GKU9CNmL5XgemYa4J8k+mg53Ymv92dhwVWGunVVL5qLq/iMUSEh7QYgHcY3nIQip5vKwC+fRP40pDoB5ZQWoletEcQhCChDlIcTGO28YWxheOzTIa4cG437Zd48eZ66zZWw/P8dGRc11MHNpXM5/5MgRjhw5EpdzZQOW6zFKexCyF7F5CLHxr1uNz5+EDij5nR9fBODlHy2M2yWHRtwMHX6RYu0eS6ueWY5t2dq4XePpp58G4NFHH43bOTMZy/UYpT0I2YsoDyE2Vn0qYvaNV+fH/ZIfHztM8dDpsf08h42q67446Wm5oVi+fHnczpUNWK7HKO1ByF7EYC6kNVprXm38Hzgud4ylzVm4nCX194nXXEFIAGIwF7KCro62AMVhUzDn+s+L4hCEFCPKQ0hrzh7eH7CfX7WEgspFqRFGEIQxRHkIsXHvF4wtDLc90MFtD3SEzY+FwZ4uXOc+Ckir+sRfxOXcwWzbto1t27Yl5NyZiOV6jNIehOwlJQZzpdTDQDtQAaC13h2hbBmw2dy9AWiKVF7IHjoP/QGvn0lOF82mevGVqRNIEIQxkm4wV0ptB97SWjeG2g9VXmu9xW+/DdgeTYGIwTyz0f0XePu5HYy4xqfnlqzeyCeu/WQKpRKE7CedDeabgxRFE9AQqqDZ6wiOjb4L2BKiuJAteL30vP3rAMUxkltGzZW1KRRKEAR/kqo8lFKh/v3dhI+HXgHUK6X8FYiTiQpFyCbOtNDTdTwgyV5zCwV5sixJENKFZP8bKzCUhT/OcIW11u1AeVDyWqA5VHml1GZM+8jChfFb4Sz48QNzBfY3QhuXv/2j8wB8/4FZkzv/YDeetj/QPTA6ltRTsIirlie21/HUU08B8OUvfzmh18kULNdjlPYgZC/JVh6hPNh1gzFEpbUOq0h8ZTB6KbeHyjftILvBsHlMTVQhJIf+HDH7jcPDkz+31nD0t/T2D+E2LeUeWy5nZ6/h8zPjt5o8FEePHk3o+TMNy/UYpT0I2UuylYcTc4aVH7EEgNgD3K61lqg9qeLrWyNmf+/+ysmfu+td6DnJxf6RsaQTZTexZN4c7LbELgq8++67E3r+TMNyPUZpD0L2kmzl0c3E3kcZgIVex8PALlEcKea6GyNm37yqcHLnHemHtpdxe730DBpDVr3587kwYzm3zSme3DljYMWKFQm/RiZhuR6jtAche0mqwdx88AcriQrC2DB8KKXWA61a62ZzP5yBXchUPm4G9wg9Ay68GrzKQXvFZykpzKW6NP7OFgVBmBqpmKr7rKkMfKzFmH4LgFKqxj/fVBQVwEEzr8Y8RkgFB140tjDse7Wffa/2x3bOix/D+Q+NrwPGkNWp0tWMOIq5ck4xKgl+rFpaWmhpaYlecJpguR6jtAche0n63EetdYNS6mFTQdQAbUHrPtZjKIdG00DeZKbv8isTckGhkAR+/gPjc80dIbN3/tKYTLfuM0XWzucehWO/A2DU46V30MVAbiVdxVcDsCIJQ1YA+/btA6Curi4p10t3LNdjlPYgZC8pmTivtX4iSt4T5ncnIO5T04lbvhgx+85Pxzgr6vgBGL4MwMW+ETSK9opbQNmoLM6jsihvspLGRG2tLED0x3I9RmkPQvYi8TyE1OE8BW//bwC8WvPOKSfHC66mo9wwwn5mWSU3LI5lMp4gCFMlnd2TCAJ4XHBkfKy8e2CUy8zgdKkxbJRjV1w9tzRV0gmCEAVRHkJsOC8ZWxguOt1cdLrD5o9x4o8waIyrazRdvcO0VdyC15YDwFVzSyjItcdFZCv09fXR19eXtOulO5brMUp7ELIXcRYkxMY37zE+f/LbkNkbv90JwMs/iuAexnkKTr05tts37KYtZxmX8+eNpV23INgrTWLZuXMnAI8++mhSr5uuWKpHiNoehOxFlIcQG6WRbQqKjRkAAA2LSURBVBAzS6P0FgYuwuFGwxWJyakBBx1lnxrbr6maQcWM3CmJGStFRRZnh00TotajjyjtQchexGAuJI/hy/D2z8dmVwEMuzw8PXIjPfnjb7jr6+azoGKSK9UFQZgSYjAX0gvXEBx6JkBxaDRv6pUBiqOqOI/55QWpkFAQhBgQ5SEkHo8L3ms0hqz8aHMso0VdHZBWt6g8KSvKBUGYGqI8hNj47teNLQwNj3XR8FjXeILXCx/8GnpPB5S7XLSEF0dWgZ+iWFBRyJVJWlEezO7du9m9O2Jk42nFhHoMR5T2IGQvYjAXYqPj44jZx065xnfM+BxcPBZQxlMyn38fWY1He8fS8nJsfO4Ts1PW6+jqsvCgnEYE1GMkorQHIXsRg7kQGydMRbB4Wcjsox2GO/XlczF8Vp09HJCvZ1TyUu5aPrgwGpB+xzXVSfNjFYrOTmNq6ty5c1MmQzoxVo8Lo8x6i9IehMzDqsFceh5CbER5SCxfmAt956Dl1zAYuHhM5xXzx4Lb+OBsoOJYWV2cUsUBojSCiao0fIjSmLaI8hDiy7kP4KMXwBu0Ojknn9ayz9NyxhOQXF6Yw60rJhnvXBCElCEGcyE2fv1zYwvFhSP8dO9xfvreysD0GZV8OOcveeVM4BBpUZ6Dv75+Pvk5yXNDEo79+/ezf//+VIuRNvx0n5Of7osY3NMgUnsQshpRHkJsPP8LYwvmUht88Gue/OAanvxg1Xh69bWcXLyel054A4rn5dj4q+vnUVqYk2CBrXHgwAEOHDiQajHShidfvMyTL16OXjBcexCyHhm2EmLjrnsmpvWchMPPgdfD1646ZKQpBSvu4ELRCvYdPIXXb2KGw6b40nXzqCpOTqwOK6xZsybVIqQVX7ujxFrBUO1BmBbIbCthapx9D478BryBtgyuvJP+iqt4+s0O+obH7R9KwbpVc7lilviSEoR0RGZbCYlFazj+Cpx8bWLess/RX3EVv2o5HaA4AD67rEoUhyBkAaI8hNg4cQxcgzB0FC5NXCB2dMbnGBxeyTsHT9EzGLjQ7LoFZdQuLEuWpDEh6zwCkXUeQjREeQix8b2/Nz7vCbIR2Byw4ovc/90S4CJ3/XWg4qipmsGa5VVp67dqz549gMTz8HH/42cBC/E8fO1B4nlMO0R5CNZwj8Cxl6AixJBTbiEDy/+KNy7mU1o2NCG7pmoGd15Tjc2WnooDoLq6OtUipBXLFlicBbfwisQKIqQt9q1bt6ZahoSwdevWrf39/dx6661jaU899RTPPfcc1dXVVFZWAtDS0sKePXvo6+tjxYoVgBGS9PHHH6elpYWbb7557Pjdu3ezb98+li9fTnGxsSJ6//79/OxnPwNg8eLFgDEEsnPnTo4ePUpdXd3Y8Tt27KCpqYm6ujry8oyZRs8//zxPP/00xcXFY0MmR44c4Yc//CFdXV1cc801Y8dv27aNAwcOJP+e+s7S+cf/zc7ftXO0Yg51t8wfv6c/e2g6NsSHajFn+z0sWuJhlvd1hj56BZVTwLIlC7jzmmraPj6WXvcUVE+bN28eq6uMrSfi1/bqlvbxyH+5Kfo9XRqm+u7/KyPuKRvrKRH3tH///q6tW7dG9RIqPQ8hPNoLHX+G4wdgOMhRns1O/5wbGeYNYAiXV09YNDSrJI91q+ZiT+MehyAIk0Om6gqh6W6Hj38/IQYHwEBOOR+UruHNizmMur0T8ksLcrhp6UyunFOctjYOQRBCI1N1hdgZ7oULR+HCRwHxNzSaoVEPzkEXM5vfQ6scXv3SbUCg4sjPsfO7F3Nx2BT/+XGLi8zShB07dgDw0EMPpViS9GDDI2cA2PvYvMgFH/qK8bnjlwmWSEg3RHlMV7SG/nPQewb6OuFy1wQvuG6vl7O9w1zoG2HAa6e9/LN8beRN8phoFF9cWcjaq+aw99lOIPN6s/39/akWIa241OuJXgigtzuxgghpiwxbTRe8Hug/D5c7wXkSnB1GXPHgYloz4vZysX+Es73DuL1wfsZyTpXdgMteyIyhXgAGCkoByHXY+MwVlayaX4pSiotOY1FgZVlmvZf09fUBjBk5pzuW69FpvnCUzUywREKykGErAQa74eJRI5JfX9eYCxGv1gy5PAyNehh2eRh2exlxeRhxexl1e8f6DZfz5nCy7CYG8qrGTulTGvPLC1hZXcIVs4oCvOJmmtLwIUojEMv1KEpj2pKSf7pS6mGgHagA0FpHnBYWa/lsxOvVuL0az8ggnsFutHsI5R4Fzyju0WFco8O4R4ZhtB+bawDbSC8M9+L2ePF4NaNub4CSCN3fVPTlzaG7YBHdBYsZySkdy7HbFAsqClhSWURN1QxK8tPDG64gCKkh6cpDKbUdeEtr3ejbV0qt9+1PtfxUGO7pQgcN5WgNbq2Nh7f5IPZo8HijD/dpbTzwR9weRt068BjtweYZRXld4HEx6jHKuDxePOb1PB43erjfUAijlylw9ZLjnTjUNFk8tlx68+bSnzeb/twqBnIr8dgC3VHk59ipXVjGtQvKjB7Gk/9iZHztH0Kec8cvjGGMh+7JrDfS559/HoC77rorxZKkB5brMUp7ELKXVPQ8Nmutt/jtNwFbgHDKINbyk+bwq/+O6/yxeJ/WMjYgkU7KR+xFDORW0p83i968uQzkVoIKXJ2hlBGkqbQghyWVM7hmfil5Dr9gTa/8xvgM87B44U8DQOYpj9bWVkCUhw/L9RilPQjZS1KVh1KqNkRyN1Afj/JCIBpFX341lwqW0FOwiFHHuGuRkoIclszIpXxGLmUFOZSaW3G+A4c9Qoywr34j4jUf/EpFvMRPKuvWrUu1CGmF5XqM0h6E7CXZPY8KjIe/P5FiXcZUXim1GdgMsHBhFIduGYhdgbbl4Motxe0owGvLxWPLRdnzsOfkYs/NQ+cU4nbMwO0oQueVkpOXz6IcOyty7GMKorQgh1zHJINIrrkjYva6z2Smu3V/FxFCDPUYpT0I2UuylUcof9zdAEqpMq11sGKIqbxpSN8NxlTdWIXzFlQyPGMgIE0Z18KuwGazmZ8Km1JYWTttsyly7DYcNhXgpkMrhdeWi7bngj2HHJuNXIcNh11hVzbsNrDbbDjyi8kpLCF3RimqoBzyy8Am0YMFQUgtyVYeTswZU35E6h/HWn5KfLL+bxN16uzhnTeMz+tuDJn92qFBAG5eVZgsieLCkSNHAMacxE13LNdjlPYgZC/JVh7dTOxNlAGE6HVMpryQaP51q/EZJn7Dd35s+MKKGgcizXj66acBiefhw3I9RmkPQvaSVOWhtW5VSgU/9CuA5niUF5LAqk9FzL7x6vwkCRJfli9fnmoR0grL9RilPQjZS9LdkyildgFN/us2CFzHUQPU+u1HLB8OcU8iCIIQO1bdkyTd8qq1bgBqlFLrzZXjbUGKYD3QEEN5QRAEIcmIY0RBEARhjLTteQgZzr1fMLYw3PZAB7c90JFEgeLDtm3b2LZtW6rFSBss12OU9iBkL6I8BEEQhJjJ2mErpdQF4GSq5UgAlcDE2LBCPJDfNjHI75oYEvW7LtJaV0UrlLXKI1tRSh20Mh4pxI78tolBftfEkOrfVYatBEEQhJgR5SEIgiDEjCiPzGPaRVFMIvLbJgb5XRNDSn9XsXkIgiAIMSM9D0EQBCFmRHkIgiAIMSPKQxAEQYiZZMfzEIS0w3S42Y4ZaMyMSClMAaVUGWZIaOAGDM/Y8rvGGaXULtN5bNIR5ZGhyJ8zPoQICbBdKbVePDdPmUe01lt8O0qpNqWUKOY4YrbdmlRdX5RH5iJ/zviw2f93BJqALYAoj0livtgEP9R2Yfyu0j7jgFKqNtUyiM0jA4ny5xQsEuYP2A3UJ1uWLKMCqDcDu/lwksK35CxkNcaLTsoQ5ZGZyJ8zPlRgKAt/gsMeCzGitW7XWpdrrdv9ktci4aPjglJqPfBsquWQYasMxPxTlgcly58zdspCpHWD0bvTWosiiQNmT7keuD3VsmQ65m/p1Fo7lVIplUWURxYgf85J48ScYeVH8L4wdfYAt2utW1MtSBawMV3smjJslR3In3NydDOx91EGIL2O+GBOg94lbXPqmDa6tBldkJ5HmmCOY26KUqw7eE63/Dknj9a6VSkVrCQqSKM/aCZjtulWrXWzuV/v+y5MCp+t07d/A1BjPgMag2xMCUccI2Yw5p/TKX/OyaOU2oWxRmZsnQd+6z6EyaGUqseYwPEs40OBDUHTooUpoJTajPGb1qXk+qI8MhP5c8YPvxXmNRjKOC3GlDMV0wbXEyKrUWu9IdnyZCM+xYHRZncDjyV7qFWURwYif05BEFKNKA9BEAQhZmS2lSAIghAzojwEQRCEmBHlIQj/f3t3bIMwDIRh9LIBrJIdWIGdYBVGYAV2gQ1CgZFofxLZAr3XJZW7T06kOyAmHgDExAOAmHgAEBMPAGLiAUBMPACIiQcAMSPZoZM2BflQr13zx/d7gxj5RW4e0MHHIp9HVV1aMHZVNWScNqzl5gGdtL3Tc1Wd2vN58JHga6bqQkfTNN2XZdmPPges5bMVdNIWeN1GnwO2IB7Qz6GqrqMPAVsQD+hnriq70fkL/nkAEHPzACAmHgDExAOAmHgAEBMPAGLiAUBMPACIiQcAsSf3OM6flQhcAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# display results\n",
    "\n",
    "SMALL_SIZE = 16\n",
    "MEDIUM_SIZE = SMALL_SIZE\n",
    "BIGGER_SIZE = SMALL_SIZE\n",
    "\n",
    "pyplot.rc('font', size=SMALL_SIZE)          # controls default text sizes\n",
    "pyplot.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n",
    "pyplot.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n",
    "pyplot.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "pyplot.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n",
    "pyplot.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n",
    "pyplot.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title\n",
    "\n",
    "pyplot.rc('text', usetex=True)\n",
    "pyplot.plot(bin_edges_0[1:], cdf_0/cdf_0[-1], alpha=0.5, linewidth=4, label='non-white')\n",
    "pyplot.plot(bin_edges_1[1:], cdf_1/cdf_1[-1], alpha=0.5, linewidth=4, label='white')\n",
    "pyplot.axvline(0, color='gray',linestyle=':',linewidth=2)\n",
    "pyplot.axhline(med_0, color='royalblue',linestyle=':',linewidth=2)\n",
    "pyplot.axhline(med_1, color='tomato',linestyle=':',linewidth=2)\n",
    "pyplot.axvline(up_0, color='royalblue',linestyle=':',linewidth=2)\n",
    "pyplot.axvline(lo_0, color='royalblue',linestyle=':',linewidth=2)\n",
    "pyplot.axvline(up_1, color='tomato',linestyle=':',linewidth=2)\n",
    "pyplot.axvline(lo_1, color='tomato',linestyle=':',linewidth=2)\n",
    "pyplot.axhline(0.05, color='gray',linestyle=':',linewidth=2)\n",
    "pyplot.axhline(0.95, color='gray',linestyle=':',linewidth=2)\n",
    "\n",
    "pyplot.title('CDF of signed residuals, $R=Y-\\hat{Y}$')\n",
    "pyplot.ylabel(r'$\\mathbf{P}(R \\leq r | A)$')\n",
    "pyplot.xlabel(r'$r$')\n",
    "pyplot.legend(loc='center right', bbox_to_anchor=(0.4, 0.75))\n",
    "\n",
    "if save_figures:\n",
    "    filename = \"CDF_combined_q.png\"\n",
    "    pyplot.savefig(filename, bbox_inches='tight', dpi=300)\n",
    "    \n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that $\\mathbb{P}\\{R \\leq r | A=0\\} \\neq \\mathbb{P}\\{R \\leq r | A=1\\}$. In particular, when comparing the two functions at $r=0$, we see that the fitted regression function overestimates the response of the non-white group and underestimates the response of the white group, as\n",
    "$$\n",
    "\\mathbb{P}\\{Y \\leq \\hat{Y} | A=0\\} = 0.51 > 0.45 = \\mathbb{P}\\{Y \\leq \\hat{Y} | A=1\\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}