{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we'll generate a random matrix" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#Number of columns (features)\n", "K = 5\n", "\n", "#Number of records\n", "N = 1000\n", "\n", "#Generate an NxK matrix of uniform random variables\n", "X = np.random.random([N, K])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's peak at our data to confirm it looks as we expect it" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.9786862 , 0.50933013, 0.47327265, 0.57154566, 0.17725616])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X[100, :]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(1000, 5)" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This exercise is about designing a scoring function for a logistic regression. As we are not concerned with fitting a model to data, we can just make up a logistic regression.

\n", "\n", "For quick intro, the Logistic Regression takes the form of $\\hat{Y} = f(x * \\beta^T)$, where $x$ is the $1xK$ vector of features and $\\beta$ is the $1xK$ vector of weights. The function $f$, called a 'link' function, is the inverse logit:

\n", "\n", "
$f(a)=\\frac{1}{1+e^{-a}}$


\n", "\n", "In this notebook we'll write a function that, given inputs of $X$ and $\\beta$, returns a value for $\\hat{Y}$.\n", "

\n", "First let's generate a random set of weights to represent $\\beta$.\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-0.21854798, 0.15267024, -0.07139126, 0.19223845, 0.99029314])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beta = np.random.uniform(low= -1.0, high= 1.0, size=(K,))\n", "beta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how we applied a neat NumPy trick here. The numpy.random.random() function returns an array, yet we applied what appears to be a scalar operation on the vector. This is an example of what NumPy calls vectorization and broadcasting, which offers us both a very fast way to do run vector computations as well as a clean and concise method of coding. \n", "\n", "\n", "\n", "Question: we designed the above $beta$ vector such that $E[\\beta_i]=0$. How can we confirm that we did this correctly?" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.20905252029504284" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#start by taking the mean of the beta we already calculated\n", "beta.mean()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "#It is likely the above is not equal to zero. Let's simulate this 100k times and see what the average mean is\n", "\n", "sims = 100000\n", "means = []\n", "for i in range(sims):\n", " means.append(2 * (np.random.random(K) - 0.5).mean())\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's use matplotlibs hist function to plot the histogram of means here. " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAASpElEQVR4nO3df6zddX3H8edrVAn+gPGjYG27lc3OCGTiaLpuLhkL26gYByaSXbJIF0mqBBJN3JKyJZuJaQLLlIRskNRBKMSJxB+jmTBFcDFOBC8EKKUiV6lS20EVAuwPWIrv/XE+1317e3rvub/O6ZXnIzk53/P+fj/n+z7fc9pXvz/OaaoKSZJ+ZdQNSJKODgaCJAkwECRJjYEgSQIMBElSs2zUDczVKaecUmvWrBl1G5K0pDz44IM/rarl/eYt2UBYs2YN4+Pjo25DkpaUJD860jwPGUmSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJKAJfxNZelotWbLV0a27j1Xv3dk69bS5x6CJAkwECRJjYEgSQIGCIQkq5N8I8nuJLuSfLTVP5HkJ0kebrcLOmOuSjKR5Ikk53fq5yTZ2eZdlyStfmySz7f6/UnWLPxLlSRNZ5A9hIPAx6vqHcAG4IokZ7R511bV2e12J0CbNwacCWwErk9yTFv+BmAzsLbdNrb6ZcDzVfU24Frgmvm/NEnSbMwYCFW1v6oeatMvAbuBldMMuRC4rapeqaqngAlgfZIVwPFVdV9VFXALcFFnzPY2/QXgvMm9B0nScMzqHEI7lPMu4P5WujLJo0luSnJiq60Enu4M29tqK9v01PohY6rqIPACcHKf9W9OMp5k/MCBA7NpXZI0g4G/h5DkTcAXgY9V1YtJbgA+CVS7/xTwIaDfv+xrmjozzPv/QtU2YBvAunXrDpsvdY3y+wDSUjTQHkKS19ELg89W1ZcAquqZqnq1qn4OfAZY3xbfC6zuDF8F7Gv1VX3qh4xJsgw4AXhuLi9IkjQ3g1xlFOBGYHdVfbpTX9FZ7P3AY216BzDWrhw6nd7J4weqaj/wUpIN7TkvBe7ojNnUpj8A3NvOM0iShmSQQ0bvBj4I7EzycKv9DXBJkrPpHdrZA3wYoKp2JbkdeJzeFUpXVNWrbdzlwM3AccBd7Qa9wLk1yQS9PYOx+b0sSdJszRgIVfUt+h/jv3OaMVuBrX3q48BZfeovAxfP1IskafH4TWVJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJwACBkGR1km8k2Z1kV5KPtvpJSe5O8mS7P7Ez5qokE0meSHJ+p35Okp1t3nVJ0urHJvl8q9+fZM3Cv1RJ0nQG2UM4CHy8qt4BbACuSHIGsAW4p6rWAve0x7R5Y8CZwEbg+iTHtOe6AdgMrG23ja1+GfB8Vb0NuBa4ZgFemyRpFmYMhKraX1UPtemXgN3ASuBCYHtbbDtwUZu+ELitql6pqqeACWB9khXA8VV1X1UVcMuUMZPP9QXgvMm9B0nScMzqHEI7lPMu4H7gtKraD73QAE5ti60Enu4M29tqK9v01PohY6rqIPACcHKf9W9OMp5k/MCBA7NpXZI0g4EDIcmbgC8CH6uqF6dbtE+tpqlPN+bQQtW2qlpXVeuWL18+U8uSpFkYKBCSvI5eGHy2qr7Uys+0w0C0+2dbfS+wujN8FbCv1Vf1qR8yJsky4ATgudm+GEnS3A1ylVGAG4HdVfXpzqwdwKY2vQm4o1Mfa1cOnU7v5PED7bDSS0k2tOe8dMqYyef6AHBvO88gSRqSZQMs827gg8DOJA+32t8AVwO3J7kM+DFwMUBV7UpyO/A4vSuUrqiqV9u4y4GbgeOAu9oNeoFza5IJensGY/N8XZKkWZoxEKrqW/Q/xg9w3hHGbAW29qmPA2f1qb9MCxRJ0mj4TWVJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpGaQXzuVtESs2fKVkax3z9XvHcl6tbDcQ5AkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpmTEQktyU5Nkkj3Vqn0jykyQPt9sFnXlXJZlI8kSS8zv1c5LsbPOuS5JWPzbJ51v9/iRrFvYlSpIGMcj/qXwz8E/ALVPq11bVP3YLSc4AxoAzgbcCX0/yW1X1KnADsBn4DnAnsBG4C7gMeL6q3pZkDLgG+PM5vyIddUb1//xKmp0Z9xCq6pvAcwM+34XAbVX1SlU9BUwA65OsAI6vqvuqquiFy0WdMdvb9BeA8yb3HiRJwzOfcwhXJnm0HVI6sdVWAk93ltnbaivb9NT6IWOq6iDwAnByvxUm2ZxkPMn4gQMH5tG6JGmquQbCDcBvAmcD+4FPtXq/f9nXNPXpxhxerNpWVeuqat3y5ctn17EkaVpzCoSqeqaqXq2qnwOfAda3WXuB1Z1FVwH7Wn1Vn/ohY5IsA05g8ENUkqQFMqdAaOcEJr0fmLwCaQcw1q4cOh1YCzxQVfuBl5JsaOcHLgXu6IzZ1KY/ANzbzjNIkoZoxquMknwOOBc4Jcle4O+Bc5OcTe/Qzh7gwwBVtSvJ7cDjwEHginaFEcDl9K5YOo7e1UV3tfqNwK1JJujtGYwtxAuTJM3OjIFQVZf0Kd84zfJbga196uPAWX3qLwMXz9SHJGlx+U1lSRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkScAAgZDkpiTPJnmsUzspyd1Jnmz3J3bmXZVkIskTSc7v1M9JsrPNuy5JWv3YJJ9v9fuTrFnYlyhJGsQgewg3Axun1LYA91TVWuCe9pgkZwBjwJltzPVJjmljbgA2A2vbbfI5LwOer6q3AdcC18z1xUiS5m7GQKiqbwLPTSlfCGxv09uBizr126rqlap6CpgA1idZARxfVfdVVQG3TBkz+VxfAM6b3HuQJA3PXM8hnFZV+wHa/amtvhJ4urPc3lZb2aan1g8ZU1UHgReAk/utNMnmJONJxg8cODDH1iVJ/Sz0SeV+/7KvaerTjTm8WLWtqtZV1brly5fPsUVJUj9zDYRn2mEg2v2zrb4XWN1ZbhWwr9VX9akfMibJMuAEDj9EJUlaZHMNhB3Apja9CbijUx9rVw6dTu/k8QPtsNJLSTa08wOXThkz+VwfAO5t5xkkSUO0bKYFknwOOBc4Jcle4O+Bq4Hbk1wG/Bi4GKCqdiW5HXgcOAhcUVWvtqe6nN4VS8cBd7UbwI3ArUkm6O0ZjC3IK5MkzcqMgVBVlxxh1nlHWH4rsLVPfRw4q0/9ZVqgSJJGx28qS5IAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCRjg104laSZrtnxlZOvec/V7R7buXzbuIUiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQL8/xBeM0b5e/WSlgb3ECRJwDwDIcmeJDuTPJxkvNVOSnJ3kifb/Ymd5a9KMpHkiSTnd+rntOeZSHJdksynL0nS7C3EHsIfVdXZVbWuPd4C3FNVa4F72mOSnAGMAWcCG4HrkxzTxtwAbAbWttvGBehLkjQLi3HI6EJge5veDlzUqd9WVa9U1VPABLA+yQrg+Kq6r6oKuKUzRpI0JPMNhAK+luTBJJtb7bSq2g/Q7k9t9ZXA052xe1ttZZueWj9Mks1JxpOMHzhwYJ6tS5K65nuV0bural+SU4G7k3xvmmX7nReoaeqHF6u2AdsA1q1b13cZSdLczGsPoar2tftngS8D64Fn2mEg2v2zbfG9wOrO8FXAvlZf1acuSRqiOQdCkjcmefPkNPCnwGPADmBTW2wTcEeb3gGMJTk2yen0Th4/0A4rvZRkQ7u66NLOGEnSkMznkNFpwJfbFaLLgH+tqv9I8l3g9iSXAT8GLgaoql1JbgceBw4CV1TVq+25LgduBo4D7mo3SdIQzTkQquqHwDv71H8GnHeEMVuBrX3q48BZc+1FkjR/flNZkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqVk26gYkaT7WbPnKSNa75+r3jmS9i8k9BEkSYCBIkhoDQZIEGAiSpMaTykM2qhNgkjQT9xAkSYCBIElqDARJEmAgSJIaA0GSBBxFgZBkY5Inkkwk2TLqfiTpteaouOw0yTHAPwN/AuwFvptkR1U9PtrOJKm/UV5Cvli/o3RUBAKwHpioqh8CJLkNuBBYlEDwuwCSdLijJRBWAk93Hu8FfnfqQkk2A5vbw/9J8kRn9inATxetw4VjnwtnKfQI9rmQlkKPsMh95pp5Df/1I804WgIhfWp1WKFqG7Ct7xMk41W1bqEbW2j2uXCWQo9gnwtpKfQIS6fPqY6Wk8p7gdWdx6uAfSPqRZJek46WQPgusDbJ6UleD4wBO0bckyS9phwVh4yq6mCSK4GvAscAN1XVrlk+Td9DSUch+1w4S6FHsM+FtBR6hKXT5yFSddiheknSa9DRcshIkjRiBoIkCVhigZDk4iS7kvw8yREv6TrSz2AkOSnJ3UmebPcnLlKfM64nyduTPNy5vZjkY23eJ5L8pDPvglH12Zbbk2Rn62V8tuMXu8ckq5N8I8nu9vn4aGfeom7LmX5yJT3XtfmPJvmdQccOsce/aL09muTbSd7Zmdf3vR9Rn+cmeaHzXv7doGOH3Odfd3p8LMmrSU5q84a2PeekqpbMDXgH8HbgP4F1R1jmGOAHwG8ArwceAc5o8/4B2NKmtwDXLFKfs1pP6/m/gV9vjz8B/NUQtudAfQJ7gFPm+zoXq0dgBfA7bfrNwPc77/mibcvpPmudZS4A7qL3XZsNwP2Djh1ij78PnNim3zPZ43Tv/Yj6PBf497mMHWafU5Z/H3DvsLfnXG9Lag+hqnZX1RMzLPaLn8Goqv8FJn8Gg3a/vU1vBy5anE5nvZ7zgB9U1Y8WqZ8jme/2GMb2nHEdVbW/qh5q0y8Bu+l9+32xTfdZm3QhcEv1fAf41SQrBhw7lB6r6ttV9Xx7+B163wMatvlsj2Fty7ms6xLgc4vUy4JbUoEwoH4/gzH5l8NpVbUfen+JAKcuUg+zXc8Yh39ormy78Dct1qEtBu+zgK8leTC9nw+Z7fhh9AhAkjXAu4D7O+XF2pbTfdZmWmaQscPqsesyens0k4703i+0Qfv8vSSPJLkryZmzHLsQBl5XkjcAG4EvdsrD2p5zclR8D6ErydeBt/SZ9bdVdccgT9GntuDX1k7X5yyf5/XAnwFXdco3AJ+k1/cngU8BHxphn++uqn1JTgXuTvK9qvrmXPrpZwG35Zvo/eH7WFW92MoLti37rbJPbepn7UjLDOVzOpv1JPkjeoHwB53yor73s+zzIXqHVf+nnQv6N2DtgGMXymzW9T7gv6rquU5tWNtzTo66QKiqP57nU0z3MxjPJFlRVfvbbvuzc13JdH0mmc163gM8VFXPdJ77F9NJPgP8+yj7rKp97f7ZJF+mt9v8TRZoey5Ej0leRy8MPltVX+o894Jtyz4G+cmVIy3z+gHGDqtHkvw28C/Ae6rqZ5P1ad77offZCXmq6s4k1yc5ZZCxw+yz47A9/yFuzzn5ZTxkNN3PYOwANrXpTcAgexxzMZv1HHaMsf3FN+n9wGML2t3/m7HPJG9M8ubJaeBPO/0MY3sO0mOAG4HdVfXpKfMWc1sO8pMrO4BL29VGG4AX2qGvYf1cy4zrSfJrwJeAD1bV9zv16d77UfT5lvZek2Q9vb+/fjbI2GH22fo7AfhDOp/XIW/PuRn1We3Z3Oj9gd4LvAI8A3y11d8K3NlZ7gJ6V5r8gN6hpsn6ycA9wJPt/qRF6rPvevr0+QZ6H+gTpoy/FdgJPErvw7ZiVH3Su5rikXbbNeztOWCPf0Bvt/1R4OF2u2AY27LfZw34CPCRNh16//nTD1of66Ybu0jv80w9/gvwfGfbjc/03o+ozytbH4/QO/n9+8PeloP02R7/JXDblHFD3Z5zufnTFZIk4JfzkJEkaQ4MBEkSYCBIkhoDQZIEGAiSpMZAkCQBBoIkqfk/ViamNUmM2yEAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plt.hist(means)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We should expect the distribution to be centered around zero. Is it? As fun technical side, let's dive a little deeper into what this distribution should look like. The histogram shows a distribution of the average of a sample of 5 uniformly distributed random variables taken over N different samples. Can we compare this to a theoretical distribution?
\n", "\n", "Yes we can! We sampled each $\\beta_i$ from a uniform distribution over the interval $[-1, 1]$. The variance of a sample of uniformly distributed variables is given by $(1/12) * (b - a)^2$, where $b$ and $a$ are the min/max of the support interval. The standard error (or the standard deviation of the mean) of a sample of size K with with $Var(X) = \\sigma^2$ is $\\sigma / \\sqrt(K)$.
\n", "\n", "Given the above knowledge, we should expect our distribution of averages to be normally distributed with mean = 0 and var = $(12 * 5)^{-1} * (1 - (-1))^2 = 0.66667$. Let's compare this normal distribution to our sample above." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwU9Z3/8denaw65zxG5BA9u8YABOVRQo6BcxisqnpGHS9Rsstlk40+TbDauWXNrVtQlxiRGo/HmZvBGQYRBuYdL5BhRGW4FhJmqz++PauJkGJiemer59vF5Ph796KO+XfWmuudD9beqviWqijHGmPQXcx3AGGNMNKygG2NMhrCCbowxGcIKujHGZAgr6MYYkyGsoBtjTIbIqamBiDwOjAa2qeppR2kzHHgAyAW2q+qwmubbtm1b7dq1a63CGmNMtlu8ePF2VS2oblqNBR34M/AQ8ER1E0WkJfAwMFJVN4vI8YmE6tq1K8XFxYk0NcYYEycim442rcYuF1WdC+w8RpPrgBdVdXO8/bZaJzTGGFNvUfShdwdaicibIrJYRG48WkMRuU1EikWkuKysLIJFG2OMOSyKgp4D9AdGASOAH4tI9+oaqupkVS1U1cKCgmq7gIwxxtRRIn3oNSkl3BG6D9gnInOBM4C1EczbGGNMgqLYQp8CnCsiOSLSGDgbKIlgvsYYY2ohkcMWnwaGA21FpBT4T8LDE1HVR1W1RERmA8uAAHhMVVckL7Ixxpjq1FjQVfXaBNr8CvhVJImMMcbUiZ0pakxVFYfArhNg0lAUO0WNSW8HP4cPnoLV02HrEjj0OXj50K43dLsYzroBWnZ2ndKYGllBN9krCKD4j/DavXBwD7TrC2deC02Ohy93w9YP4K1fwju/g0HfgvP+A/Kbuk5tzFFZQTfZ6cAueO5m2PAmnDwcLvgJdOp/ZLvdm+HN+2Heg7C2CK5+Agp6NGxWYxJkfegm++zaBI99DTbNhzEPwg0vV1/MAVqeCJc9DDdOgX3b4bGLYMvChs1rTIKsoJvssncrPDEW9pXBjVOh/80gUvP7Th4Ot70JTdrAE5fB5gVJjWlMXVhBN9nj4Ofw18th3w64/iXoMrh272/ZGW6ZBc1OgKevgR0fJienMXVkBd1khyCAlybC9rVwzZNH72KpSbMTYPxzIDF46ir4cm+0OY2pByvoJjvMeyA8LPHie8Puk/pocwp840nY9RFM/zc7Zt2kDCvoJvN9shTeuA/6fB0G3R7NPLsMgeF3w4rnYclT0czTmHqygm4yW8XBsKulcVsY9dvEdoAm6tzvQddzYdZdsOfj6OZrTB1ZQTeZ7e3fwLZVMPZ/oXHraOcd88L5BhUw6z+inbcxdWAF3WSunR/BOw9A36ug+8XJWUbrk2D4XWH/fMm05CzDmARZQTeZq+geiOXART9L7nIG3wntToOiu8MuHmMcsYJuMtP6V2HNDBj2A2jeIbnL8nLg4v8Ohwl47/+SuyxjjsEKusk8QQBzfgKtToruqJaanHJ+ODLj3F+HJy4Z44AVdJN5Vr4I21bCBT+CnPyGW+5FPwuH3p37y4ZbpjGVWEE3mcWvgDd+Dsf3gT6XN+yyj+8FZ46H4sfDMWOMaWA1FnQReVxEtonIMa8TKiIDRMQXkSuji2dMLS39G+z8EC64B2IOtlfO+z5oEA63a0wDS2Q89D8DDwFPHK2BiHjAL4CiaGIZUwd+RdiH3aEf9Li0Vm/teteMyGL8Iuccxi34I+e+1ZcyWlXbZuP9oyJbnjGH1bgJo6pzgZ01NPs28AKwLYpQxtTJqpdh9yY499+jPSO0lib548jBZ2LOdGcZTHaq929SEekIfB14NIG2t4lIsYgUl5WV1XfRxnxFNRyAq233Wm+dR22ztmNKMJRrvddpwRdOs5jsEkUn4wPAD1XVr6mhqk5W1UJVLSwoKIhg0cbEffg6fLochvyrm77zKiZXjKKxHGS895rrKCaLRPHNLwSeEZGNwJXAwyJyWQTzNSZx8x6AZu3h9KtdJwFgjZ7IXL8vN+UUkUuF6zgmS9S7oKvqSaraVVW7As8Dt6vqy/VOZkyiPlkKH82FQd9q2OPOa/CYfyntZDejY++6jmKyRCKHLT4NvAv0EJFSEblVRCaKyMTkxzMmAQsnQ25j6HeT6yT/ZG5wOmuDjkzImQnYRTBM8tV42KKqXpvozFT15nqlMSYu0cMIW/AF7+X/nRf9c7j7v+YlOVVtCY/5l/LL3D8wOLaKd4M+rgOZDOd+75Ex9XCV9xbHSTl/9ZM0PG49TfGHskub2s5R0yCsoJu0JQRc773KwqAHJdrFdZxqHSSPF/xzGRFbRFv2uI5jMpwVdJO2hsWW0TX2GX+tuMh1lGP6m38hueJzlfeW6ygmw1lBN2nrBu8VyrQFs4OBrqMc0wbtwLt+b67zXkMIXMcxGcwKuklL7dnB+bElPOOfT3lCQxK59ZR/IZ1jZZwXW+46islgVtBNWrrce5uYKM/6w1xHSUhRMIDt2pzx3quuo5gMZgXdpCHlKu8t5vu92aLtXIdJSDk5POcP48LY+xzPLtdxTIaygm7SzkBZTdfYZzyXJlvnhz3nD8MT5TLvHddRTIaygm7SztU5b/G5NmJWiu8MrWqDdmBx0I0rvLfD0SGNiZgVdJNWmnCAS2PvMc0fzJekzrgtiXrBP48esVLY+oHrKCYDWUE3aWWUt4DGcjDtulsOm+4P4qDmwtKnXUcxGcgKukkrV3tvsS7oyAd6qusodbKXJswJ+sPy56DioOs4JsNYQTdp40T5jMLYWp73zwPcXWKuvp73h8GBXbDWLsFromUF3aSNcbF5BCpM8Ye4jlIvbwd9oekJ1u1iImcF3aQJ5TJvHgu1J5/SxnWYegmIhVdWWjcH9m13HcdkECvoJi30kU2cEvsk7bfO/+H0b0BQAavs4l4mOlbQTVoY683jkHrM9M92HSUa7fpAQU9Y/oLrJCaDWEE3KS9GwFjvXd4KzmQPTV3HiYYInHYlbJ4Pe0pdpzEZIpFrij4uIttEZMVRpo8XkWXx23wROSP6mCabDYytpr3szJzulsNOuzy8X/Gi2xwmYySyhf5nYOQxpn8EDFPV04F7gckR5DLmH8bG5rFP83k16Oc6SrTanAIdzoIVz7tOYjJEjQVdVecCO48xfb6qHh4+bgHQKaJsxpBHOZd6CykKBqTlqf41Ou1K+GQpbF/vOonJAFH3od8KzDraRBG5TUSKRaS4rKws4kWbTDQstpSWso+pmdbdcthplwMCK2znqKm/yAq6iJxPWNB/eLQ2qjpZVQtVtbCgoCCqRZsMNsZ7lx3ajHeC01xHSY7mHaDL0LDbxUZgNPUUSUEXkdOBx4Bxqrojinkak88hLoh9QJFfSEUaXGauzvpeAdvXwqd2eTpTP/Uu6CJyIvAicIOqrq1/JGNCQ2IraSpfUpRm457XWq9xIJ51u5h6S+SwxaeBd4EeIlIqIreKyEQRmRhv8hOgDfCwiCwRkeIk5jVZZERsEXu1EfODPq6jJFeTNnDSeVAy1bpdTL3U+DtWVa+tYfoEYEJkiYwBPHwu8hbzenAW5Znc3XJY73Ew/bvw2Qo4oa/rNCZN2ZmiJiUNiK2hjXzObD/Du1sO6zkaJAarprpOYtKYFXSTkkbEFvGl5vJWcLrrKA2jaUF4tMuqKa6TmDRmBd2kHlVGeIuYG5zOAY5znabh9B4H29fAttWuk5g0ZQXdpJ6t79NBdjLbH+A6ScPqNQYQ20o3dWYF3aSekmmUq8drmTZ2S02anQAnDraCburMCrpJLapQMo0FQa/MGSq3NnqPhW0rYfs610lMGrKCblJL2RrYsZ6iIMu6Ww7rNSa8t610UwdW0E1qKZkGCEV+oeskbrToBJ0GWEE3dWIF3aSWkqnQaQBltHKdxJ3e4+DTZbBzg+skJs1YQTepY9emsJAd7nbIVr3Ghvd2kpGpJSvoJnWsnh7e9xrtNodrrbqEVzIqsYJuascKukkdJdOgXV9ofbLrJO71HA0fL4a9W10nMWnECrpJDV9sg80LbOv8sJ7x9bB6htscJq1YQTepYfUMQK3//LCCHtDm1K+6oYxJgBV0kxpKpoVdLcf3dp0kNYiEW+kb34EDu2pubwxW0E0qOLAbPpobH0JWXKdJHT1HQ1ABa+e4TmLShBV04966ORCUf3W4ngl17A9NT4DV01wnMWnCCrpxr2QqNGsfFjDzlVgMeo6C9a9B+QHXaUwaSOSaoo+LyDYRWXGU6SIivxeR9SKyTESybIg8Uy+H9sO6V8PuhZhtXxyh5ygo3w8fvuE6iUkDifwF/RkYeYzplwDd4rfbgEfqH8tkjQ9fh4oDdrji0XQ9F/Jb2OGLJiE1FnRVnQvsPEaTccATGloAtBSR9lEFNBmuZBo0ahVefs0cKScPul8Ma2aCX+E6jUlxUfzG7QhsqfS8NP7aEUTkNhEpFpHisrKyCBZt0ppfDmtnQY9Lwct1nSZ19RwNB3bC5nddJzEpLoqCXt1xZlpdQ1WdrKqFqlpYUFAQwaJNWtv4Nny556uzIk31Tv0aePnW7WJqFEVBLwU6V3reCbABKEzNSqZBbhM45XzXSVJbftNwHa2eHl7RyZijiKKgTwVujB/tMgjYo6qfRDBfk8kCH0qmQ7eLILeR6zSpr+do2LMlHF7YmKPIqamBiDwNDAfaikgp8J9ALoCqPgrMBC4F1gP7gVuSFdZkkNJFsG9b1o7d0vWu2nWftCaHRfnCQ5N+x+8qrqr18jbeP6rW7zHpp8aCrqrX1jBdgTsiS2SyQ8k08PKg28Wuk6SFnTSnWHswIlbM76h9QTfZwc7kMA1PNSzoJw+H45q7TpM2ivwB9IxtoYt86jqKSVFW0E3D+3Q57N6Utd0tdTUnCC+cfXGs2HESk6qsoJuGt3o6SCw8/twkrFQLWBl0YYRnBd1Uzwq6aXgl0+DEIdCkreskaafIH0A/WUcBu11HMSnICrppWDs+hG2rrLuljoqCQmKifM1b7DqKSUFW0E3DKomP7W2DcdXJGu3MpuB4Rlg/uqmGFXTTsEqmQYezoEUn10nSlFAUDGBIbAXN2O86jEkxVtBNw9nzMXxcbN0t9TTbH0Ce+JwfW+I6ikkxVtBNwzk8uJRdaq5ePtBT2aYtGeEtdB3FpBgr6KbhrJ4GBT2hbTfXSdKaEmOO35/hsaXkc8h1HJNCrKCbhrFvB2ycZ0PlRqQoGEATOcg5seWuo5gUYgXdNIy1s0B96z+PyIKgN3u1MSNji1xHMSnECrppGCXTocWJ0P4M10kyQjk5vBr042ve+3j4ruOYFGEF3STfwc/Di0H3Gg1S3QWuTF0U+YW0ki8YGFvtOopJEVbQTfKtewX8g9bdErG5wekc0DxGWLeLibOCbpJv9XRoUgCdz3adJKMc4DjmBqczwitGCFzHMSmgxgtcGHNYba+yA5DPIRbnz2CqP4S7756dhFTZrcgvZIRXzOmygaV6qus4xrGEttBFZKSIrBGR9SJyVzXTW4jINBFZKiIrRcQuQ2cAGBpbQVP5kqJggOsoGem1oB/l6tmQugZIoKCLiAdMAi4BegPXikjvKs3uAFap6hmE1x/9jYjkRZzVpKGRsUXs1cbMD/q4jpKR9tCUBUGveD+6uo5jHEtkC30gsF5VN6jqIeAZYFyVNgo0ExEBmgI7gYpIk5q0k0MFF3mLeTXoR7n17iVNUTCAU2KfcKp87DqKcSyRgt4R2FLpeWn8tcoeAnoBW4HlwHdU1fbSZLmzYyW0ki+Y7Vt3SzLN8cNL09mQuiaRgl7dgcNVf9uNAJYAHYAzgYdE5Iir/4rIbSJSLCLFZWVltQ5r0sslsYXs03zeCuxkomTaRiveD05lpA3WlfUSKeilQOdKzzsRbolXdgvwoobWAx8BPavOSFUnq2qhqhYWFBTUNbNJAzECRnjFvBGcyUFsd0qyzfYH0De2kY7YhlI2S6SgLwK6ichJ8R2d1wBTq7TZDFwIICLtgB7AhiiDmvTST9ZSIHsosu6WBnH4KCI72iW71VjQVbUCuBMoAkqAZ1V1pYhMFJGJ8Wb3AkNEZDnwGvBDVd2erNAm9V3iLeKg5vJ6cJbrKFlhk55ASdCZEZ6dNZrNEjr0QFVnAjOrvPZopcdbgYujjWbSlzLCW8TcoC/7aOQ6TNaYEwzgTu8l2rCHHbRwHcc4YKf+m8idLhvoJNvtZKIGVuQX4onyNe9911GMI1bQTeRGeosoV49X/P6uo2SVVdqFzUGBDdaVxaygm4gpI2MLeTfozR6aug6TZYSiYABDYytoxn7XYYwDVtBNpHrIFk6OfcrsYKDrKFlppn82+VLB12KLXUcxDlhBN5G6xFtIoGLdLY58oKfysbZhlLfAdRTjgBV0E6kRsUUs0h6U0dJ1lCwlzPTP5tzYcpqzz3UY08CsoJvInCxb6RXbYmO3ODbDH2TdLlnKCrqJzOjYAgINtxCNO0v0FEq1LaO891xHMQ3MCrqJzChvAYu0B5/R2nWULCfM8M/m3NgymvOF6zCmAVlBN5HoJqX0iJUy3R/kOoohPNolT3wu9qzbJZtYQTeRGO29i6/CbN8OV0wFSw93u8TsaJdsYgXdREAZHVvAe0EvO7olZQjT/UGcE1th3S5ZxAq6qbfesolTYp8wPRjsOoqpZIY/iFzxbUjdLGIF3dTbaG8BFRpjlh2umFKW60lsDgoYFbOjXbKFFXRTT8qo2ALmB33YxRFXHTROCTODQQyNrYD9O12HMQ3ACrqpl77yEV1i25ge2NEtqWi6fza54sPq6a6jmAZgBd3Uy2jvXcrVs0vNpagVehIbg3aw/HnXUUwDsIJu6kwIGOW9x9tBXxsqN2UJU4Ih8NFc2PuJ6zAmyRIq6CIyUkTWiMh6EbnrKG2Gi8gSEVkpIm9FG9Okov6ylk6ynWm+Hd2Syqb4QwGFFS+4jmKSrMaCLiIeMAm4BOgNXCsivau0aQk8DIxV1T7AVUnIalLM17157Nd8u9RcitugHaDDWbD8WddRTJIlsoU+EFivqhtU9RDwDDCuSpvrgBdVdTOAqm6LNqZJNXmUM8pbwJygP/s5znUcU5O+V8MnS6FsjeskJokSKegdgS2VnpfGX6usO9BKRN4UkcUicmN1MxKR20SkWESKy8rK6pbYpIThsSW0lH287J/jOopJxGlXgMRgmW2lZ7JECrpU85pWeZ4D9AdGASOAH4tI9yPepDpZVQtVtbCgoKDWYU3quMybx3ZtzttBX9dRTCKatYOTh8Py50Cr/vmaTJFIQS8FOld63gnYWk2b2aq6T1W3A3OBM6KJaFJNc/ZxYewDpvmD8fFcxzGJ6ns17N4EWxa6TmKSJJGCvgjoJiIniUgecA0wtUqbKcC5IpIjIo2Bs4GSaKOaVDHSW0i+lPOSdbekl16jIacRLPu76yQmSWos6KpaAdwJFBEW6WdVdaWITBSRifE2JcBsYBmwEHhMVVckL7Zx6euxeXwYtGeZnuw6iqmN/GbQ81JY+RL45a7TmCTISaSRqs4EZlZ57dEqz38F/Cq6aCYVdWA7g71V/Kb8SqrfvWJSWt+rw+PR178GPUa6TmMiZmeKmloZ680H4OVgqOMkpk5OvRAat4GlT7tOYpLACrqpBeVKby7FQXe2aDvXYUxdeLnhVvqamTYCYwaygm4S1k/WcWpsK8/6w1xHMfVx1njwD4WHMJqMYgXdJOwq7y32aT4z7ELQ6e2EvtD+DPjgr66TmIhZQTeJOfgFY7x3meEPYh+NXKcx9XXWDfDp8nA4AJMxrKCbxKyaQlP50rpbMsVpV4CXDx885TqJiZAVdJOYD57kw6A9xdrDdRIThcatoeeocATGioOu05iIWEE3Ndu+HjbP5zl/GHbseQY563o4sCs84sVkBCvopmZLngTxeME/13USE6WTh0PzjvDBk66TmIhYQTfH5lfAkqeh20WU0cp1GhOlmAdnjg/PGt21yXUaEwEr6ObY1syELz6FftUOcW/SXf+bQAQW/9l1EhMBK+jm2Ir/CM07QXcb9yMjtYh/th/8FSoOuU5j6skKujm67ethw5vQ/+bw57nJTIW3wr4yKKk6KrZJN1bQzdEVPw6xHOtuyXSnXACtuoaft0lrVtBN9coPwJKnoNeY8PJlJnPFYtD/Ftg0D7bZdWnSmRV0U70VL8KXu2HABNdJTEM463rw8mwrPc1ZQTfVW/QYFPSELjbueVZo0hZ6XwZLn4GDX7hOY+rICro5UmkxbH0/3FkmdmZo1hgwAQ7utYtfpLGECrqIjBSRNSKyXkTuOka7ASLii8iV0UU0De7dhyC/BZx5neskpiF1Hggd+sGCRyAIXKcxdVDjNUVFxAMmARcBpcAiEZmqqquqafcLwotJm3S1ezOsmgqD74D8pq7TmIh0vWtGQu3GxIbwv3kP8c0f3cfrQb86L2/j/aPq/F5Td4lsoQ8E1qvqBlU9BDwDjKum3beBF4BtEeYzDe29/wvvz/4XtzmME7OCgWzV1kzwbMCudJRIQe8IbKn0vDT+2j+ISEfg68Cjx5qRiNwmIsUiUlxWVlbbrCbZDn4O7z8BfS4LzyA0WaeCHP5cMYIh3ip6y0bXcUwtJVLQq9srplWePwD8UFX9Y81IVSeraqGqFhYUFCSa0TSUD54Md4oNusN1EuPQM/4F7NN8bs2Z5TqKqaVECnop0LnS807A1iptCoFnRGQjcCXwsIhcFklC0zACP9wZ1nkQdOrvOo1xaC9NeNYfzpjYfI5nl+s4phYSKeiLgG4icpKI5AHXAP806IOqnqSqXVW1K/A8cLuqvhx5WpM8K1+C3ZtgyJ2uk5gU8Cd/JB4Bt+ZYX3o6qbGgq2oFcCfh0SslwLOqulJEJorIxGQHNA0gCODt34QnEvWwoxMMbNZ2TA2GcL33Ki353HUck6CEjkNX1Zmq2l1VT1HV++KvPaqqR+wEVdWbVfX5qIOaJFo7C7atgnO+F47rYQzwcMU4mshBbsmZ7TqKSZD99WY7VZj763C0vdOucJ3GpJB12olZ/gBu8Ypoxn7XcUwCrKBnuw1vhKf5D/0ueDWeZ2ayzKSKcTSX/dzgveI6ikmAFfRsdnjrvFkHO83fVGuFnswb/hncmjOTRnzpOo6pgRX0bLbhzXAM7HO+Czn5rtOYFPW/FV+njXxuW+lpwAp6tlKF134GLTqHl5gz5ije1+686Z/Bt3KmWV96irOCnq1Wzwj7zof90LbOTY1+VXE1reQLJuQkNsiXccMKejYKfHjjPmhzKpxxres0Jg2s1JOY7p/NBG8mbdjjOo45CjusIY0lOiRqVeNi7/Bg3iruOPSvzLjHRjs2ifltxVWMzFvEHTlT+FmFXTg8FdkWepbJ5xA/yH2WFUFXZgYDXccxaWSDduB5/zzGe6/SSWyU7FRkBT3LTPBm0km2898V16P28Zta+l3Flfh43J3zN9dRTDXsLzqLHM8ubs+Zwix/AAuC3q7jmDT0Ga15uGIsl3oLGRRbVfMbTIOygp5FfpDzd3Lw+Z8KO4nI1N0f/FGUalv+M+cJYti1R1OJFfQs0Vc2cFXOXP7kj2SztnMdx6Sxg+TxP+XX0Su2mWu8N1zHMZVYQc8CHj4/z32MbdqSSRV23RFTfzOCs3kv6MkPcv5Oa/a6jmPirKBngZu8OfSNbeSn5TfyOY1dxzEZQfhR+TdpwgHuyX3SdRgTZwU9w3VgO/+e8yyv+2cyMzjbdRyTQdZpJx71x3CF9w7nxJa7jmOwgp7hlP/K/QsC/KTiFqq/3rcxdTep4jI+DNpzX84fOY6DruNkPSvoGWxcbB4XeYv5bcWVlGqB6zgmAx0kj7vLJ9Alto3v5diFylxLqKCLyEgRWSMi60XkrmqmjxeRZfHbfBE5I/qopjbas4N7c//MwqAHj/uXuI5jMth72osnKy5kgjeTwbGVruNktRoLuoh4wCTgEqA3cK2IVD0r5SNgmKqeDtwLTI46qEmcEPDr3EeJEfDv5RMJ7IeYSbL7KsbzkZ7Ar3MfpTn7XMfJWon8pQ8E1qvqBlU9BDwDjKvcQFXnq+qu+NMFQKdoY5rauMUrYqi3kp9V3MAWO+bcNIADHMe/ld/O8ezmZ7l/ch0nayVS0DsCWyo9L42/djS3ArOqmyAit4lIsYgUl5WVJZ7SJOwMWc9dOX/jFb8/z/rDXccxWWSZnsKDFZdzmTcflthYLy4kUtCrOzRCq20ocj5hQf9hddNVdbKqFqpqYUGB7aSLWiv2Minv93ymrfl++b9gR7WYhvaIP5b5fm+Y/m/wyTLXcbJOIgW9FOhc6XknYGvVRiJyOvAYME5Vd0QTzyQqRsCDuZMoYDffKv8Oe2jqOpLJQj4e3y7/NjRqDX+/Hg7sqvlNJjKJFPRFQDcROUlE8oBrgKmVG4jIicCLwA2qujb6mKYm38t5jvO85fy04iZW6Mmu45gstoMWcPUTsHcrvHhbeIUs0yBqLOiqWgHcCRQBJcCzqrpSRCaKyMR4s58AbYCHRWSJiBQnLbE5wlXem9yZM4WnK87naf8C13GMgc4D4JL7Yd0cKLrHdZqskdAl6FR1JjCzymuPVno8AZgQbTSTiKGx5fw854/M9fvyYzsb1KSSARNgxwZYMAlangiDb3edKOPZNUXTWE/ZzCO5D/ChduCO8u9QYR+nSTUX/zfs2QJFd0PzDtDHRvtMJjvjJF1tW81Tefexj0bceuj7NoqiSU2xGFw+GToPhBcmwFq7KHkyWUFPR9vXwxNj8fG47tA9fIwdAmpSWG4juO5ZaNcnPPJl3auuE2UsK+jpZttq+MsYCHyuO3Q3H2l714mMqVmjlnDDS1DQA565Dta94jpRRrKCnk62LITHR4D6cNNU1quNsGDSSOPWcONUKOgOT18DS552nSjjWEFPF2tmw1/Ghn8Ut84Jf74ak24at4abZ0KXofDyRHjnd6DVnnhu6sAKeqpThbm/DrdoCnrAN+dAq66uUxlTd8c1h/HPw2lXwqs/hZdvh/IDrlNlBDvOLZUd/Bxe/haUTIO+V8GY30OeHc1iMkBOHlz+B2jbDd68Hz5bDt940jZW6skKeqravCA8bXpPKYz4OQy6HcROGjLpoRsCdL8AAAnsSURBVOtdMxJs2Zfhse/z4CeT4IHB/Lj8FqYGQ6jNCXIb7x9Vp4yZyLpcUk35l/Dqf8Gf4lcZunkGDL7DirnJWG8GZzHm0H2s1478Pm8Sk3IfpBV7XcdKS1bQU8n6V+GRwfDOb+HM8fCtedBlsOtUxiTdZm3HVYf+k1+UX8NFscW8lv99rvFeJ0bgOlpasYKeCsrWwDPj4ckrQGJw/Ysw7iHIb+Y6mTENJiDGI/5Yxhy6j3XaiftzH+PlvB8zUEpcR0sb1ofu0s4N8OYvYPmzkNsYLvgxDPk25OS7TmaMM2v0RL5x6MeMjb3L/8v9G8/m38s7fh9+V3Eli7WH63gpTdTRMaCFhYVaXJxZo+wmtiNIGSiruSVnNhfHijlELn/xL+b/Kkazi+ZJz2hMOsnnEOO91/hWzhQKZC/vBT35U8VIXgn64+MB2bdTVEQWq2phddNsC72BtGUPY7z5XOnNpU9sE7u1CZP90TxeMZIyWrmOZ0xKOkgej/uX8LR/Ptd5r3GzN4dH8x7gY23Dc/4wXvaHuo6YUmwLPUJVt9DbsZPh3lIuiS3knNhyciRgRdCVp/wLeck/hy+xrhVjaiNGwIWx97nRm8PQ2EpiotChH5x2BXQfAW1OzfgjwmwLvYE05wv6xdYzILaa4bGl9IltAqBU2/KoP4aX/aE2/oox9RAQ45WgkFeCQk5gB6O9BfwoWA5z7glvLbtAt4vgpPOg00Bonl2D19kWel3t2w7bVsG2EvhsBWxZBGXh3vgKjbFYu/OGfyavB2exVjthVxIyJjk23j8Kdm2C9a+EQ/N+9BaU7w8ntjgxvBxe+zPg+N5wfC9o3jGtt+LrvYUuIiOBBwEPeExV768yXeLTLwX2Azer6vv1Su2SXxFerXz/9vBCt7s3h1dd2b0lvN+xHvaVfdW+UWvo2J9fbz2NxdqdpcEp7Oc4d/mNyTatuoSXvBswASoOwafLwtFJt7wXnnW94oWv2uY3D7tmWp74z7dm7aFJATRuEw5NkIZqLOgi4gGTgIuAUmCRiExV1VWVml0CdIvfzgYeid9Hzy8PxzgJKsA/FD4PKuL35eF95ceH2x3aD+X74NDh2xdfPT74BezfEb9thwO7gSq/XMSDFh3D//G7j4j/bx+/NT0eRHgo4dOdjTFJk5MHnQrD2+HrmO7fCWWrv/pVveND+HQ5rJkF/sEj53FcC2jcFpq0DR/nNQ3PC8lvFn/cNLzPbQReXnjLyf/n+388zoVYTlhDYjkQ88L35TaK/p+eQJuBwHpV3QAgIs8A44DKBX0c8ISG/TcLRKSliLRX1U8iT1wyDZ6/pd6zOaB57OM49ms++ziO3dqMHbRmp3ZlF83Yoc3Yqc35TFvxsbblM1rhH/Dg08pz2Q+kcbeRMRkg8XFjANrHbxcAIAS0ZQ+dpYwC2U0b+Zw27KFNxV7a7NtLm7K9NJPPaMKXNJUvacoBGks1/wHU1tDvwEU/q/98qkikoHcEtlR6XsqRW9/VtekI/FNBF5HbgNviT78QkTW1SvuVtsD2Or43mVI1F6RuNstVO5ardmrMtREXm2X3toV767q+uhxtQiIFvbq9B1X3pCbSBlWdDExOYJnHDiRSfLSdAi6lai5I3WyWq3YsV+1kW65ExnIpBTpXet4J2FqHNsYYY5IokYK+COgmIieJSB5wDTC1SpupwI0SGgTsSUr/uTHGmKOqsctFVStE5E6giPCwxcdVdaWITIxPfxSYSXjI4nrCPYX132t5bPXutkmSVM0FqZvNctWO5aqdrMrl7MQiY4wx0bLx0I0xJkNYQTfGmAyRsgVdRK4SkZUiEojIUQ/vEZGRIrJGRNaLyF2VXm8tIq+IyLr4fSRj1CYyXxHpISJLKt32ish349N+KiIfV5p2aUPlirfbKCLL48suru37k5FLRDqLyBsiUhL/zL9TaVqk6+to35dK00VEfh+fvkxE+iX63iTnGh/Ps0xE5ovIGZWmVfuZNlCu4SKyp9Ln85NE35vkXD+olGmFiPgi0jo+LZnr63ER2SYiK44yPbnfL1VNyRvQC+gBvAkUHqWNB3wInAzkAUuB3vFpvwTuij++C/hFRLlqNd94xk+BLvHnPwW+n4T1lVAuwvMo2tb33xVlLsJT9/rFHzcD1lb6HCNbX8f6vlRqcykwi/DcikHAe4m+N8m5hgCt4o8vOZzrWJ9pA+UaDkyvy3uTmatK+zHA68leX/F5nwf0A1YcZXpSv18pu4WuqiWqWtOZpP8YlkBVDwGHhyUgfv+X+OO/AJdFFK22870Q+FBVN0W0/KOp77/X2fpS1U80Ppibqn4OlBCeaRy1Y31fKud9QkMLgJYi0j7B9yYtl6rOV9Vd8acLCM/1SLb6/Judrq8qrgWejmjZx6Sqc4Gdx2iS1O9Xyhb0BB1tyAGAdho/Fj5+f3xEy6ztfK/hyC/TnfGfW49H1bVRi1wKzBGRxRIOxVDb9ycrFwAi0hU4C3iv0stRra9jfV9qapPIe5OZq7JbCbfyDjvaZ9pQuQaLyFIRmSUifWr53mTmQkQaAyOBSsMtJm19JSKp3y+nF7gQkVeBE6qZdI+qTklkFtW8Vu/jMI+Vq5bzyQPGAv+v0suPAPcS5rwX+A3wzQbMNVRVt4rI8cArIrI6vlVRZxGur6aEf3jfVdW98ZfrvL6qW0Q1ryU6jEVSvms1LPPIhiLnExb0cyq9HPlnWotc7xN2J34R37/xMuGoqymxvgi7W+apauWt5mStr0Qk9fvltKCr6tfqOYtjDTnwmcRHfIz/pNkWRS4Rqc18LwHeV9XPKs37H49F5A/A9IbMpapb4/fbROQlwp96c3G8vkQkl7CYP6WqL1aad53XVzXqM4xFXgLvTWYuROR04DHgElXdcfj1Y3ymSc9V6T9eVHWmiDwsIm0TeW8yc1VyxC/kJK6vRCT1+5XuXS7HGpZgKnBT/PFNQCJb/ImozXyP6LuLF7XDvg5Uuzc8GblEpImINDv8GLi40vKdrS8REeCPQImq/rbKtCjXV32GsUjkvUnLJSInAi8CN6jq2kqvH+szbYhcJ8Q/P0RkIGFN2ZHIe5OZK56nBTCMSt+5JK+vRCT3+5WMPb1R3Aj/eEuBg8BnQFH89Q7AzErtLiU8KuJDwq6aw6+3AV4D1sXvW0eUq9r5VpOrMeEXu0WV9/8VWA4si39g7RsqF+Ee9KXx28pUWV+E3QcaXydL4rdLk7G+qvu+ABOBifHHQnhBlw/jyy081nsj/L7XlOsxYFel9VNc02faQLnujC93KeHO2iGpsL7iz28GnqnyvmSvr6cJhw0vJ6xftzbk98tO/TfGmAyR7l0uxhhj4qygG2NMhrCCbowxGcIKujHGZAgr6MYYkyGsoBtjTIawgm6MMRni/wNxKeqHomPwXgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Compute a vector from the normal distribution specified above\n", "from scipy.stats import norm\n", "mu = 0\n", "sig = np.sqrt(4 / 60.0) \n", "xs = np.linspace(-1, 1, 1000)\n", "ys = norm.pdf(xs, mu, sig) \n", "\n", "plt.hist(means, density = True)\n", "plt.plot(xs, ys)\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's write our scoring function. Let's try to use as much of Numpy's inner optimization as possible (hint, this can be done in two lines and without writing any loops). The key is that numpy functions that would normally take in a scalar can also take in an array, and the function applies the operations element wise to the array and returns an array. i.e.:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 1])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ex_array = np.array([-1, 1])\n", "np.abs(ex_array)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's use this feature to write a fast and clean scoring function" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def score_logistic_regression(X, beta):\n", " '''\n", " This function takes in an NxK matrix X and 1xK vector beta.\n", " The function should apply the logistic scoring function to each record of X.\n", " The output should be an Nx1 vector of scores\n", " '''\n", " \n", " #First let's calculate X*beta - make sure to use numpy's 'dot' method\n", " xbeta = X.dot(beta)\n", " \n", " #Now let's input this into the link function\n", " prob_score = 1 / (1 + np.exp(-1 * xbeta))\n", " \n", " return prob_score" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So how much faster is it by using Numpy? We can test this be writing the same function that uses no Numpy and executes via loops." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def score_logistic_regression_NoNumpy(X, beta):\n", " '''\n", " This function takes in an NxK matrix X and 1xK vector beta.\n", " The function should apply the logistic scoring function to each record of X.\n", " The output should be an Nx1 vector of scores\n", " '''\n", " #Let's calculate xbeta using loops\n", " xbeta = []\n", " for row in X:\n", " \n", " xb = 0\n", " for i, el in enumerate(row):\n", " xb += el * beta[i]\n", " \n", " xbeta.append(xb)\n", " \n", " #Now let's apply the link function to each xbeta\n", " prob_score = []\n", " for xb in xbeta:\n", " prob_score.append(1 / (1 + np.exp(-1 * xb)))\n", " \n", " return prob_score" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before doing any analysis, let's test the output of each to make sure they equal" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mean Absolute Deviation = 0.0\n" ] } ], "source": [ "diff = np.abs(score_logistic_regression_NoNumpy(X, beta) - score_logistic_regression(X, beta))\n", "\n", "print(f'Mean Absolute Deviation = {np.round(diff.sum(), 8)}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If they equal then we can proceed with timing analysis" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6.17 ms ± 55.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "%timeit score_logistic_regression_NoNumpy(X, beta)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "19.2 µs ± 476 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit score_logistic_regression(X, beta)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.8.3" } }, "nbformat": 4, "nbformat_minor": 1 }