RyotaBannai · January 12, 2019 08:56
diff --git a/dffits.ipynb b/dffits.ipynb
 {
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.linalg import inv\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dataset_ (N=200, idx_outlier=0, ydistance=5):\n",
    "    rng = np.random.RandomState(4)\n",
    "    data  = np.dot(rng.rand(2, 2), rng.randn(2, N)).T\n",
    "    data[idx_outlier:idx_outlier+1,1] = ydistance\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=50\n",
    "inx=10 #an influential point\n",
    "i = dataset_(N, inx)\n",
    "p = len(i.T)-1\n",
    "Z = i[:,0].reshape(len(i[:,0]),1)\n",
    "y = i[:,1].reshape(len(i[:,1]),1)\n",
    "f_H = lambda i: np.dot(i, inv(np.dot(i.T, i))).dot(i.T)\n",
    "b = lambda Z,y: inv(np.dot(Z.T,Z)).dot(Z.T).dot(y)\n",
    "ypred = lambda Z,b: Z.dot(b)\n",
    "H = f_H(i)\n",
    "original_ypred = ypred(Z, b(Z,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "outliers,Dfs = [],[]\n",
    "t = 2*np.sqrt((p+1)/(N-p-1)) #cutoff value\n",
    "\n",
    "errs = y - original_ypred\n",
    "var_= sum(errs**2)/N       #for the internally Studentized residuals\n",
    "\n",
    "for n in range(N):\n",
    "    rstudent = errs[n]/(var_ * np.sqrt(1-H[n,n])) #the internally Studentized residuals\n",
    "        \n",
    "    Df = rstudent*(np.sqrt(H[n,n]/ 1-H[n,n]**2))\n",
    "    Dfs.append(abs(Df))\n",
    "    if abs(Df) > t: \n",
    "        outliers.append(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.408248290463863"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[10]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "outliers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([10.38312882]),\n",
       " array([0.13462557]),\n",
       " array([0.12596528]),\n",
       " array([0.11952226]),\n",
       " array([0.10717736])]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dfs.sort(reverse=True);Dfs[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x11e07cf28>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFpCAYAAABJdYvCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3WlsXNd99/HfXebeGe7LiMtwkUTttC0pqu20QY1mUdEgCIwiL1zHSICmL4LCRdKicOs0aewXrgu2tZEAhYO8Sdoi3ZAUcYMCAYryCQI/bYIniR0lrmRL1i6K63C4c7a7PC8o06aphZoZce5wvh8jgHh5585/Tij9eM499xwjDMNQAAAgEsxqFwAAAN5BMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARYlfrjcfHx6v11pGTTCaVTqerXUbNod1KQ7uVhnYrDe22JpVKbflceswAAEQIwQwAQIQQzAAARAjBDABAhBDMAABECMEMAECEEMwAAEQIwQwAQIQQzAAAREhFVv76gz/4A8XjcZmmKcuyNDIyUonLAtipPE/mtWsyPU/ewIAUj1e7IiAyKrYk57PPPquWlpZKXQ7ADmVkMnL/+78VhqFkWbL/939VPH5c/t691S4NiASGsgFsK+enP1UYj0uJhOQ4ChsbFfvFL6RCodqlAZFQsR7z888/L0n6zd/8TZ08ebJSlwWwk2SzMlZXFTY2bjxuGDKnphQMDFSnLiBCjDAMw3Ivkslk1NHRoYWFBf3FX/yFPvOZz2h4eHjDOaOjoxodHZUkjYyMqMBvx+ts25bnedUuo+bQbqWparsVizK++13pvcG8sqLwN35D6umpTl1bwM9baWi3NY7jbPncigTzu337299WPB7Xo48+etvz2PbxHWyLVhrarTTVbrfYj38sM5ORYrG1A0EghaHyH/2oZBhVq+tOqt1utYp2W7Ot2z7mcjlls9n1P//yl7/U4OBguZcFsEMVH3pIQTIp4+1h7Xhchd/4jUiHMrCdyr7HvLCwoBdeeEGS5Pu+fv3Xf13Hjx8vuzAAO5Rtq/j+96t4o6csy6p2RUCklB3M3d3d+pu/+ZtK1AKgnpg8FALcDH8zAACIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiJCKBXMQBPrTP/1TjYyMVOqSAADUnYoF8/e//3319fVV6nIAANSligTz7OysXnvtNX3kIx+pxOUAAKhbFQnmv//7v9enPvUpGYZRicsBAFC37HIv8Oqrr6q1tVVDQ0M6ffr0Lc8bHR3V6OioJGlkZETJZLLct94xbNumPUpAu5WGdisN7VYa2u3uGWEYhuVc4J//+Z/1yiuvyLIsFQoFZbNZPfzww/r85z9/29eNj4+X87Y7SjKZVDqdrnYZNYd2Kw3tVhrarTS025pUKrXlc8vuMT/xxBN64oknJEmnT5/Wf/zHf9wxlAEAwM3xHDMAABFSdo/53e677z7dd999lbwkAAB1hR4zAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEWKXe4FCoaBnn31WnufJ93396q/+qh577LFK1AYAQN0pO5hjsZieffZZxeNxeZ6nZ555RsePH9fBgwcrUR8AAHWl7KFswzAUj8clSb7vy/d9GYZRdmEAANSjsnvMkhQEgZ5++mlNTk7qt37rt3TgwIFKXBYAgLpjhGEYVupiKysreuGFF/SZz3xGg4ODG743Ojqq0dFRSdLIyIgKhUKl3rbm2bYtz/OqXUbNod1KQ7uVhnYrDe22xnGcLZ9b0WCWpO985ztyXVePPvrobc8bHx+v5NvWtGQyqXQ6Xe0yag7tVhrarTS0W2lotzWpVGrL55Z9j3lxcVErKyuS1mZov/766+rr6yv3sgAA1KWy7zHPzc3ppZdeUhAECsNQv/Zrv6Zf+ZVfqURtAADUnbKDeffu3frrv/7rStQCAEDdY+UvAAAihGAGACBCKvIcMwAAtS4IAk3PzmhhaVlNDQ3q2dUly7K2vQ6CGQBQ93zf16k3z8jzPMVsW/OLCxqfntL7hu+XbW9vVDKUDQCoe2NTE/JvhLK0tg+EDENXxse2vRaCGQBQ95aWlzf1jC3T1NKNdTq2E8EMAKh7MTum9y6EGYbheg96OxHMAIC6N5jqU76Q33AsVyhoMLX9K1kSzACAupeIx/XAwSNrm274vkzL0n37D6q5sWnba2FWNgAAkpqbmnT00JFql0GPGQCAKCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBAelwIAlG1ucUFjE+Mqep4SiYSG+gfkOm61y6pJ9JgBAGXJzM/rjfPnVPQ8SdLKyopOvXFanu9XubLaRDADAMpybeK64m58/WvTNGXI0PjUVBWrql0EMwCgLAWvuOmYbdtaza1WoZraRzADAEoSBIHGp6eUzmQ0Oz+nIAjWv1coFtXW0nLT13m+r7HJCZ2/clmLy0vbVW7NYPIXAOCueb6vU2+clu/7am1u0YWrlzQ7P699g7sV+L5cx1FXR3LT61ZWV/XLs2/INE3ZlqWp2bS6Ozu1f/feKnyKaKLHDAC4a9fGrysMAsVsW3HX1aGh/WqMJ7SwvKj+nl4dOzws09wcMW9dvSwnFlPMtmUYhhKuq+nZtJZXVqrwKaKJHjMA4K4trSzLsqz1r23LVl9Pj2zLUl9P701fE4ahVrOrir/nMSon5ig9P6emxsZ7WnOtoMcMALhrsVhMYRhuOBaGoWKx2C1fYxiGbGtzf9APfMVdnnl+G8EMALhrAz0p5Qv5DcdyhbwGe/tu+7qe5C4Viu/M4g7DUIZhqKuj857UWYsIZgDAXWtqbNTwgUMyLUue78u0LN1/4JAaGxpu+7qB3pRSXd3yA1+FYkFuzNEDh47c9H50veIeMwCgJG3NLTp+ePiuXmMYhgZTfRpM3b5nXc/4FQUAgAghmAEAiJCyh7LT6bReeuklzc/PyzAMnTx5Uh/72McqURsAAHWn7GC2LEuf/vSnNTQ0pGw2qy984Qs6evSo+vv7K1EfAAB1peyh7Pb2dg0NDUmSEomE+vr6lMlkyi4MAIB6VNF7zNPT07p06ZL2799fycsCAFA3jPC9S7eUKJfL6dlnn9UnPvEJvf/979/0/dHRUY2OjkqSRkZGVCgUKvG2O4Jt2/JubDCOraPdSkO7lYZ2Kw3ttsZxnC2fW5Fg9jxPf/VXf6Vjx47p4x//+JZeMz4+Xu7b7hjJZFLpdLraZdQc2q00tFtpaLfS0G5rUqnUls8teyg7DEN9/etfV19f35ZDGQAA3FzZs7LPnj2rV155RYODg/qTP/kTSdInP/lJnThxouziAACoN2UH8+HDh/Xtb3+7ErUAAFD3WPkLAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAaCOrWZXlZ6fU6FYrHYpuKHsx6UAALUnCAKdPn9OSyvLMmQoVKjUrm7t6R+odml1jx4zANShS9evKZvNKu64ch1HccfV9ekpLSwtVru0ukcwA0Adml9YlG1vHDSNO44mZqarVBHeRjADQB0yjVsdv8U3sG0IZgCoQ8mOzk0TvnKFvFJdPVWqCG9j8hcAVIjn+8oX8nIdV7ZlVbuc2+rv6VWhWNB0ZlZhEMi2be0b2KOmxsZql1b3CGYAqIDLY9c0kZ5WGAQyDFO9u7oiPcPZMAztG9yj3X0D8jxPTiwm02QQNQoIZgAo03QmrYmZKbmOu35sYmZKDQ0JdXUkq1jZndmWFfnefb3h1yMAKNNUOr0hlCXJdVxNpdNVqgi1jGAGACBCCGYAKNOujk7lC4UNx/KFgnZ1dFSpItQy7jEDQJm6O5NaWl7WzNyswjCUYRja1d6p7s5d1S4NNYhgBoAyGYahA3v2andfv7K5rMJQml2Y0+Xr19TX3SMn5lS7RNQQghkAKsSJxTSdTuvK+DW5jqtQ0sT0lA7u3adkO8Pa2BruMQNAhRS9oq5NjisRT8g0TVmmqbgb16VrVxWGYbXLQ40gmAGgQhaWl296vOB5myaHAbdCMANAhcRjjoKb9IxNw9i0kxNwKwQzANylMAy1mssql89vON7U2KjGeFx+EKwf83xPnW1trK6FLeNXOAB1byYzq8n0jMIwVLK9Q727umTcYvvD5ZUVvXHxvArFgiRDDfG4hvcfWF/56/5DR/TWlUtaXFqSaRrqbO/Q3r7orpmN6CGYAdS1q+PXdX1qUq7j3Ph6TEsryzq0d9+mc4Mg0Jnz52TbthJufP3YGxfO6/iR+yStrT19ZGj/9n0A7DgMZQOoW77va3x6aj2UJcmJOZqdy2wappakxeVleb6/4ZhhGFrNrt7oQQPlI5gB1K1CsbjhfvA7DK1mszc5Hko3RrjDMJTneQrDUKEM8TQUKoWhbAB1y4nFZN10D+JQjQ0Nm462NDXLNC3NZGaVWVhQEASyTFNdyaScWExXxseUzswpCH21NjVraHAPk75w1yoSzF/72tf02muvqbW1VS+++GIlLgkA95xlWUp1dWtsakLxG5O38oWCDJl648I5eX6glsYmDQ3uliSZpqld7W1669JFmeZa19m0TBmGodfPndVqdkVBECoIAs2HS3r97Bt63/D9Vft8qE0VCeYPfvCD+uhHP6qXXnqpEpcDgG2xms2q6HtyHVdBGMp1HBmScoWCgiCUaRhaWF7SL8+eUXdXlyRpJZvT8P4DyhfykmEo7rhrk8LeOivTNFX0ipIMWaapzvYOLS4vqaWpuaqfE7WlIveYh4eH1dTUVIlLAcC2mJiZ1s/P/K/m5ucVBoHy+bw629pV8P0Nk8Es01Q+X9DsXEbS2oQxwzAUd+PrvewgDDQ1m5a0NnnMicVkWZYmZ6a1mrvZvWrg1pj8BaDuBEGgK9fHlIjHZZqmTNNUIh7X5bFr8r3ipvMty9LSyookqb21VUXP2/D9XD6vhOtuevbZD3yZBv/M4u5s2+Sv0dFRjY6OSpJGRkaUTCa3660jz7Zt2qMEtFtpaDdpaXlZDQ0NakgkNhy3s1lJoRoSGyd+reZy6tnVpUQ8rvaODp06/boWl5dkWbZ839Oe5G4ZpqmZzKxisZgMGSoUC0r19Kinu0cdbW3b+OmihZ+3u7dtwXzy5EmdPHly/et0Or1dbx15yWSS9igB7VYa2m1tgtfy8pK84sbecS6f09DAbl28dlVOLCbTNJUrrA1xJ+Lx9Xbb3dunpeVlrWRX1NLUrIZEgzKZjGzT0tzignzfV1dHpxw7Jq9QqOv25udtTSqV2vK5PC4FoO64jqPW5hatrK6uby7heZ5am1vUs6tLrS0tGp+aVKHoabCvX52tm3u8zU1Nan7X3Jrh/Qd0+q1zamtuVihDTiymI0P7Zd70cSzg1ioSzF/96ld15swZLS0t6fd///f12GOP6cMf/nAlLg0AZcnl81rNZtXU2CAn9s6kriP7DujitSuaW1yQJHW0tWloYO2xqIQb177BPXf1Pk7M0fuG71cun1cQBkq48Vuutw3cTkWC+Y/+6I8qcRkAqJggCHT20gVlFhbeXqxLXZ2d2r97r6S1Z5Lf/nMlxV234tdEfWEoG8CONDY5qcXlZSXeFZQzmYzamluU7OisYmXA7XHzA8COlFmYU8ze2PdwHUdTTERCxBHMAHYk8yb3d8MwlGFy3xfRxlA2gMjzfX9tkpZhqKOldX2mcxAEujY5rszcnCSps71DA70pGYahrmRSl8auyX3XhK9cobA+wQuIKoIZQKTNLS7o7MXz69sqGqah4aEDamlu1rnLl7S4tLj+yNPEzLRy+ZwO7t2nnmSXcrm8JtMz8oNAtmVqd2+f2ltbq/hpgDsjmAFEVhAEOnfpomJ2TIViUbZlybIsnb10QYeG9un1s2cUhKFMw1RHW5vaW1qVnp/TnmJRTiymPf0DGuhNqVAsynUcnilGTSCYAUTW8uqqZucyWlheVhD4MmSoualRrU0t+n+/+LnyhYLcGxtJTKVnFAaBEvGEijeCWVpb5zrBnsioIfz6CCCyCsWCpmbTsi1LTsxRLBbTSnZVVybG5DiOTPOdwI3ZMWUWFmRZFs8So6YRzADuuTAMtbSyrIWlRQVBsOXXzS8uqCGRUKhw/ZhpmMrmcmpKNKizrU2FYlHhjf+y+az6untk0UNGDWMoG8A9tZrN6vT5cyoUi5JCWZatw3v3qa2l5Y6vDYJQg719Gp+eVD5fkAwp7rjav3uv8oWCdnV0qjHRoNmFecmQ+ntSGujd+mYBQBQRzADuqTcvviXTMDaswHX20gU99MCxO07G6t3VpZm5jPb2D8r3fUlSKKmjtVWFoqeFpUU1JBJyHUee7+uBg4fv5UcBtgXBDOCeyeXzyuXzirvxDcd939fi8vIde83NTU0a6O3V9clJheHagHZLU5OGBnbLMAxlFhaUzszKdR2lunrWJ3wBtYxgBnDP3Gp3pTAMtdWNlwZ7+5Ta1a3l1RW5jqtE/J2Q72xrU2fb5i0ZgVpGMAO4Z1zHUUOiQb7vr4d0GIaKxWJqaWre8nVs21ZbCwuDoD4wKxvAPTW8/6BMw9D84qIWV5alMNR9+w+yVzFwCwQzgHtqbnFe+WJhbaJXKDmuy3PGwG0QzAC2JAgCjU9P6c2L5zU2ObE+S/p28oWCzl++LCfmqKWpSS1NTcrlcrpw9co2VAzUJu4xA1gXBIFWc1nZlr2hV+v7vk69eUae5ylm21pYWtLE9JSOD9+nmH3rmdDj05NyHWfDMcuy1naKAnBTBDMASdLkzIwuX78q3w8kw1BLU5OG9x2QZVlrPeQboSxJMduWHwS6cn1M+3fvrXLlwM7CUDYA5fJ5Xbh2RU7MUSIeV8J1lX3XkPPSyvL61opvs0xTy6urt71uqqtHhWJhwzHf99XODGvglghmAJqYnpL7nsU5bMvS/NLakHPMjikMww3fD8Pwjgt6uI6jfbv3qFAsKpvLKZvPKx6Pa9/g7sp+AGAHYSgbwG35vi8Z0vmrl9QQb9Cujk45sZhyhYIO7h264+u7O3dpV3vnTe9dA9iMHjMA9XZ1K18sbjjm+b5aGpt16s0zWlhcVF93r4rFos5fvqh8oaD79h9Uc2PTlq5vmqaaGhoJZWAL6DEDUNx1tX9wt85evKCF5SXF3bh6d+1Sc2OjFpYXFbNjsm1be/oHJL29Etedd4cCcPcIZgAKw1DzS4syTFOtN5bKLHqe5leWNz0OVSwWdWXiuoIgkOs4GuztU2NDQzXKBnYkhrIBaHZhXrPzc2pMJNTU2KimxkZ5nqfFxQUV3zXE7XmeLoxdle/5a888Z7P6xdkzWrnD7GwAW0cwA9DMbFpxZ+P9X8uyZFsxmaYpPwjWzpuble/56u5MSlrbPcqNOboyPrbtNQM7FcEMQJZlK7gRvu9mxyy9b/h+dbS2yjRNhZKGBgc3DF0bhqFCobDptQBKwz1mAOrv7tHPM7NKvGvWdKFQ0O6+ftm2rX2DeyRJbS0tmp6d3fDaMAzlxhPbWS6wo1Wkx3zq1Cn94R/+oT73uc/p3//93ytxSQDbqCGR0KG9QwrCULlCXp7vKdXdo96u7g3n9Xf3StL6BhZBEKhQLGpPf/+21wzsVGX3mIMg0De+8Q39+Z//uTo7O/Vnf/ZnevDBB9XPX1SgpiTbO5Rs75Dn+zINY22bxvewbVsnhu/XtYlxrWRXFXddDab65MScm1wRQCnKDubz58+rp6dH3d1rv1l/4AMf0E9/+lOCGahRtmXd/vu2rb0Dg9tUDVB/yh7KzmQy6uzsXP+6s7NTmUym3MsCAFCXyu4xv3dhe2ltluZ7jY6OanR0VJI0MjKiZDJZ7lvvGLZt0x4loN1KQ7uVhnYrDe1298oO5s7OTs2+a5bm7Oys2tvbN5138uRJnTx5cv3rdDpd7lvvGMlkkvYowU5rN9/3tbSyItu21NTQeM/eZ6e123ah3UpDu61JpVJbPrfsYN63b58mJiY0PT2tjo4O/ehHP9LnP//5ci8L1JXp2VlduHpZQRhIMhR3XT1w8PAdt1UEsPOUHcyWZen3fu/39PzzzysIAn3oQx/SwMBAJWoD6kLRK+r81UsbVt4KgkDnLl7Q/YcOV7EyANVQkQVGTpw4oRMnTlTiUkDdmclkZJkbZ0Kbpqml1RUFQXDTx5ZuJwgCZRbmFYahOlrbZN1hljWAaGHlL6DKTEPSzSZRavMkyjtZXFrSmQtvSQrXL3lo7z51tLVtOC8zP6drE+Mqep7i8biG+gfVkGD1LiAKWCsbqLJkR1LBe4LZDwK1trTcVW85DEOdvXRBTiwmJ+bIddb+d+7yxQ3rYM9kZvXmxfPyb/TG8/m8fvHmGRW94m2uDmC7EMxAldmWpSP7Dsj3feXyOeXyOcVdVwf27N1w3sLSkq6Mj60PU79XNp9TwfM2HfcDX8urK+tfXx67prgbX//aMAzZlqXrU5MV/FQASsVQNhABbS0tevCBY8rl87Isa8Ns7CAIdObCW1pYWpITi2l8akpx19XRQ0c23D9eu0+9ObCljctrejcJb8uylM3lKvmRAJSIHjMQEYZhKBGPb3pEanp2RssrK0q4rizTlOs4KnqeLr9nD2TXcdTc0Lhh2DoMQ8VdZ8Nz0YlEYlOPu1Asqr2l9R58KgB3i2AGIm52bn49rH3f1/j0lK5NXNfrb76hhaXFDece2X9AccdVrpBXrpBXzLZ1/4FDG845tHdIRc+Td2OHqEKxqLjrqquT1ZmAKGAoG4g4y7YUFkKFYaiL165KhmQaphQGOv3WOR3ed0AdrWu93Zgd0/2HDq+FbhjKtjf/FU/EE3rw/qO6PjWpbC6nVHePujo67/qxLAD3BsEMRFx/d69+efaMVrJZhWEgy7RV9IpKdfUo7rq6Oj62Hsxv28oOUbv72AEOiCJ+RQYirqmxUQf37lM2l1MQhgrDQN0dSbU0NUlaG4oGsHPQYwZqQLK9Qw8fe58uXbsq13E2fC/uurd4FYBaRI8ZqBFdHZ2Ku2szsqW1x6hyhbz29rM2PbCT0GMGaoRpmjp6aFjXp6a0uLwoJxbTQG+fEvH4nV8MoGYQzEANsSxLg6mUpK3v7QqgtjCUDQBAhBDMAABECMEMAECEEMxAGQrFolZz2Q3rUwNAOZj8BZQgCAKdvXRRc4trWzDatq2h/kHt6uisdmkAahw9ZqAEl8auaWl5SXHHVcKNK2bZOnf5ovKFQrVLA1DjCGagBLPzc5s2iHBijsanJ6tUEYCdgmAGShJuOmJICsLNxwHgbhDMqAuFYkGT6WnNLazdEy5Xa3OL/PdM+MoXC0rt6i772gDqG5O/sONdnbiusckJWYYpPwjkxGI6eviInJhz5xffwr7BPTrz1lktr64o1Nr+yHv7BlkeE0DZCGbsaNlcTtcmxpVw1wIzJikMQ527eFH3Hzpc8nVty9LRw8NazWVVKBTV1Nh4xz2QAWArCGbsaFOz6U09Y8MwtLS6UpHrN8QTaognKnItAJAIZuxwMctSEASyzI3TKUzz5tMrrk9Nano2LT8I1NLYpKHB3fSEAWwrJn+h5oVhqMmZaZ29eF7j01PyfX/9e927ujatyuX5npLtHZuuc3X8uq5OXFcYhjINQwvLS3r97JmKTBYDgK2ix4yatrC0pP/zo/+rxdVluTFHcddVqqtb7xt+QE4sJtuydN/+g3rr6mXl83mZhqnO9nbt7R/YcJ0wDDUxMy33XcPelmkqm89rfnFR7a2t2/3RANQpghk1y/M8/ejnP9NKLquV1azmCguSpKXlFbU0NevIvgOSpJbmZv3KfQ8oXyjI9325jrNpKNsPAvmBr9h7/krYpqXV7CrBDGDbEMyoWRMz01rJZpWZn5dhaP1ecHpuVheuXV0PZkkam5zQ9alJ+YEv0zDVs2uX9vS902u2Leumj08VfY9QBrCtyrrH/OMf/1h//Md/rN/5nd/RhQsXKlUTsCWFYlGe78n3fRky1o8bpqmVlaX1e8tzCwu6On5dMdtW3HHlxGKanJ7WVHpmw/WG+geUzefWX5fP59XV0amGRMP2fSgAda+sYB4YGNBTTz2lI0eOVKoeYMv6ursVj8U2LIP59sSttua29eMTM1OKu+6G1zqOo6nZ9IZjHW3tOjH8gNpbW5VIJHRwaL8O7Bm69x8EAN6lrKHs/v7+StUB3LW4G9fRQ0c0PTurbD4nz/Nk2zEdGdqvZGfn+tD2rfZKvtls60Q8rn2De+5l2QBwW9t2j3l0dFSjo6OSpJGRESWTye1668izbZv2uInp2VmNTVyX5/vqbO/Qnv6BDc8j27ath088qJaWVp0+f1ZNiUYlEnEFQajj9z2gtuZmSdKh4IDevHh+ffUvaW0YfHeqvy7bnZ+30tBupaHd7t4dg/m5557T/Pz8puOPP/64HnrooS2/0cmTJ3Xy5Mn1r9Pp9G3Ori/JZJL2eI/xqSldHh9T3FmbkDWTTuvK1Ss6enh4/Zy32y3Z1q6H7juqyXRaTiymnuQuefm80vm8JMk2TMXtmKZmpmUapsIwUHtLmxKOU5ftzs9baWi30tBua1Kp1JbPvWMwf/nLXy6rGOBuhWGosamJ9VCWpJhta3l1VQtLS2q90RN+t4ZEg4YGBm96PcMwdGDPkAZT/VpaXVFjPMFmEwAii8elUBVBECg9N6vM/ILicVeprh45sZiktWeKPd/btBSmbVlaXlm+aTBvhes4cp3Sd5QCgO1QVjD/5Cc/0Te/+U0tLi5qZGREe/bs0Ze+9KVK1YYdKgxN22O8AAAMoUlEQVRDnT53VsvZVbmOo6WVZY1PT+vYoSNqbGiQZZqK2Zt/NIu+r5YSQxkAakVZwfzwww/r4YcfrlQtqBOzcxnNzGeUzWVlmbY62lrlxmK6NHZV9x88LMMwNJjq0/krlxV3XBmGoUKxqLaWFjU3NlW7fAC4pxjKxrY7ff6cJqanFIvFFIahMgtz6u9JyX/XY03dnbvU4CY0NjUhz/eV6upSd7KrilUDwPYgmLGt8oWCFpeXZJmWDBkyDEOmbWpyZkoH9m5czKO5qUlHmg7c4koAsDOx7SO21dzigjrbOhSG4YYFPrK5nLo7dlWxMgCIBoIZd201u6rM/LyKXvGuX5u4sTTm3oFBuY6jIAxkGoZSXT3q2UUwAwBD2dgy3/d15vxbWlpZkmRIhtTf3avBVN+Wr9HS1CzXdRUGgQZ61x64L3qeku3tm7ZiBIB6xL+E2LJLY1eVzecUd+OKu67ijqtrkxNaWl7e8jUMw9DRQ4fV2NC4tjNU4Ks7mdTQwO57WDkA1A56zNiy+aXFTYt+xB1HEzPTam7a+mNMMTumI/v2V7o8ANgRCGZs2Wo2q5nMrIpFT7ZtKdneqaaGBpmmcecXS1pZXdXVievKFwpqiMe1u69fruPe+YUAUEcYysaWrGazmltYUNEryrIshaF0fWpS8wuLSnX13PH1K6ur+sWbZ7SazSoIAi0uL+vnb5xWoXj3E8gAYCcjmLEl1ybG1dfdo6ZEo4peUYViQZZpyo07akgk7vj6K+Njch1HhrHWuzZNU5Zp6drk+L0uHQBqCkPZ2JKiV5Rpmkp196jb9+UHvmJ2TNZ77jnfSqFQWA/lt1mmqVwudy/KBYCaRY8ZW9LW0qrijWFny7LkxBwFYajmxsYtvd514xsWFJEkz/fV2NBQ8VoBoJYRzNiSVFe3XNdVvlCQtNaDNgxDu1P9W3r9nv5+FYpFBTfWw/Z9X1Ko/u7ee1UyANQkhrKxSRAEunz9mjILC5KkztY27e7r19FDRzQ7P6f5xQU1NzWrq6Nzy4uCJNy4Ttz3gK6OX1e+kFdDS4sGe/tk32R7RwCoZ/yriE3evHheSysr63siz8xllMvndWT/Ae3q6NSujs6Srus6jg7s2VvJUgFgx2EoGxvkCwXNLS6sh7Ik2ZalucV5FYqFKlYGAPWBHjM2KHpFKdx8PAylQrG4NukrCOT5vmK2vWmm9c3k8nldnRhXPp9TQ0ODBntTitmxe1A9ANQ+ghkbJNy4LHPzI1CWZSrhxnV14rompqfl3wjmwb4+dXfeeleobC6rU2+cVsyOyTRNZebmNDuX0YnhB7i/DAA3wVA2NrAsSwOplLL53Pqeydl8TgO9fUrPZXR9clIx21bcdWVZls5fuaKllVtvYnH5+picmLM+SWztuWdDY1MT2/SJAKC2EMzYpK+7R8cODaulqVktTc06dmhYfd09mkrPyHWcDefGHUdjk5O3vFYun9803G1bllZWV+9J7QBQ6xhLxE01NTbqQOPGGdT+jWeQ380wDAWBf8vruI6jbC63IZx931ciHq9csQCwg9Bjxpa1t7bK87wNx3KFvHZ1dNzyNbtT/coXCuurfgVBID8MNNCbuqe1AkCtIpixZYO9fUq4cWXzeRU9T7l8Xp1t7drVkbzlaxobGnTs8LAaEg0yTUstTc06MXw/s7IB4BYYyq4Tnu8rnUkrXyyqu3OX4u7d74NsmqYeOHxEi8tLWs2uqqWpWQ2JO6913djQoCP79pdSNgDUHYK5DiyvrOj1c2/KNAyZpqmxiQntTvWrv7e0darfnhQGAKg8hrLrwLkrF+U6jmKxtW0aE/G4rk5eZyUvAIgggnmH827cC34vy7SUmZ+rQkUAgNshmHc40zRlGpv/b/YDX7HY3d9nBgDcWwTzDmeapjra2uX57zzmFIahYpat9paWKlYGALiZsiZ/fetb39Krr74q27bV3d2tJ598Uo2NjZWqDRWyf3C3Lo0ZmslkFISBmhoadXDP0Jb3UgYAbJ+ygvno0aN64oknZFmW/vEf/1Evv/yyPvWpT1WqNlSIaZraN7hHQwO7JWlLO0IBAKqjrC7TsWPHbmxKIB08eFCZTKYiReHeMAyDUAaAiKvYWOYPfvADHT9+vFKXAwCgLt1xKPu5557T/Pz8puOPP/64HnroIUnSd7/7XVmWpUceeeSW1xkdHdXo6KgkaWRkRMnkrZdxrDe2bdMeJaDdSkO7lYZ2Kw3tdveM8O3dBUr0wx/+UP/1X/+lZ555Ru5dLPM4Pj5eztvuKMlkUul0utpl1BzarTS0W2lot9LQbmtSqa1v3FPWUPapU6f0ve99T08//fRdhTIAALi5smZlf+Mb35DneXruueckSQcOHNBnP/vZihQGAEA9KiuY//Zv/7ZSdQAAALHyFwAAkUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEWJXu4ByhGGomblZzc7NyYnF1Nfdq7jrVrssAABKVtPB/Mb5t7SwvCTXcbSyuqqp2bSG9x9UW3NLtUsDAKAkNTuUvbC0pPmlRbmOI0kyTVNxx9Wla1erXBkAAKWr2WDOzK8NX79XNp+rQjUAAFRGzQZzY2ODPN/bdNyxN4c1AAC1omaDOdnWIduyFQTB+rF8saBUV08VqwIAoDw1G8ymaerYkWG1NDUpCEMZpqmhgd1KdXdXuzQAAEpW07OyY3ZMB/fuq3YZAABUTM32mAEA2IkIZgAAIqSsoex//dd/1c9+9jMZhqHW1lY9+eST6ujoqFRtAADUnbKC+dFHH9Xjjz8uSfr+97+vf/u3f9NnP/vZihQGAEA9Kmsou6GhYf3P+XxehmGUXRAAAPWs7FnZ//Iv/6JXXnlFDQ0NevbZZytREwAAdcsIwzC83QnPPfec5ufnNx1//PHH9dBDD61//fLLL6tYLOqxxx676XVGR0c1OjoqSRoZGVGhUCin7h3Ftm153uZVzHB7tFtpaLfS0G6lod3WODf2ddiKOwbzVs3MzGhkZEQvvvjils4fHx+vxNvuCMlkUul0utpl1BzarTS0W2lot9LQbmtSqdSWzy3rHvPExMT6n3/2s5/d1RsDAIDNyrrH/E//9E+amJiQYRhKJpPMyAYAoExlBfNTTz1VqToAAIAqeI8ZAACUjyU5I+ALX/hCtUuoSbRbaWi30tBupaHd7h7BDABAhBDMAABECMEcASdPnqx2CTWJdisN7VYa2q00tNvdY/IXAAARQo8ZAIAIKXsTC1TGt771Lb366quybVvd3d168skn1djYWO2yIu/HP/6xvvOd7+j69ev6y7/8S+3bt6/aJUXaqVOn9Hd/93cKgkAf+chH9Nu//dvVLinyvva1r+m1115Ta2vrlpcchpROp/XSSy9pfn5ehmHo5MmT+tjHPlbtsmoCPeaIOHr0qF588UW98MIL6u3t1csvv1ztkmrCwMCAnnrqKR05cqTapUReEAT6xje+oS9+8Yv6yle+ov/5n//R2NhYtcuKvA9+8IP64he/WO0yao5lWfr0pz+tr3zlK3r++ef1n//5n/y8bRHBHBHHjh2TZVmSpIMHDyqTyVS5otrQ39/PGu1bdP78efX09Ki7u1u2besDH/iAfvrTn1a7rMgbHh5WU1NTtcuoOe3t7RoaGpIkJRIJ9fX18e/aFhHMEfSDH/xAx48fr3YZ2GEymYw6OzvXv+7s7OQfSmyL6elpXbp0Sfv37692KTWBe8zbaCt7W3/3u9+VZVl65JFHtru8yNrqnuC4vZs9gGEYRhUqQT3J5XJ68cUX9bu/+7tqaGiodjk1gWDeRl/+8pdv+/0f/vCHevXVV/XMM8/wD+a73KndsDWdnZ2anZ1d/3p2dlbt7e1VrAg7ned5evHFF/XII4/o/e9/f7XLqRkMZUfEqVOn9L3vfU9PP/20XNetdjnYgfbt26eJiQlNT0/L8zz96Ec/0oMPPljtsrBDhWGor3/96+rr69PHP/7xapdTU1hgJCI+97nPyfO89UkmBw4cYH/rLfjJT36ib37zm1pcXFRjY6P27NmjL33pS9UuK7Jee+01/cM//IOCINCHPvQhfeITn6h2SZH31a9+VWfOnNHS0pJaW1v12GOP6cMf/nC1y4q8N998U88884wGBwfXRwA/+clP6sSJE1WuLPoIZgAAIoShbAAAIoRgBgAgQghmAAAihGAGACBCCGYAACKEYAYAIEIIZgAAIoRgBgAgQv4//l/4g58O13oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "c=['#415952', '#f35134', '#243AB5']\n",
    "carray = np.array([c[0]]*N)\n",
    "carray[outliers]='red'\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.scatter(i[:,0], i[:,1], label=r'$df=%i$' % 2, c=carray, alpha=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.5.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
 }
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 16,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"import numpy as np\n",
	"from numpy.linalg import inv\n",
	"import matplotlib\n",
	"import matplotlib.pyplot as plt\n",
	"plt.style.use('ggplot')"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": [
	"def dataset_ (N=200, idx_outlier=0, ydistance=5):\n",
	" rng = np.random.RandomState(4)\n",
	" data = np.dot(rng.rand(2, 2), rng.randn(2, N)).T\n",
	" data[idx_outlier:idx_outlier+1,1] = ydistance\n",
	" return data"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 7,
	"metadata": {},
	"outputs": [],
	"source": [
	"N=50\n",
	"inx=10 #an influential point\n",
	"i = dataset_(N, inx)\n",
	"p = len(i.T)-1\n",
	"Z = i[:,0].reshape(len(i[:,0]),1)\n",
	"y = i[:,1].reshape(len(i[:,1]),1)\n",
	"f_H = lambda i: np.dot(i, inv(np.dot(i.T, i))).dot(i.T)\n",
	"b = lambda Z,y: inv(np.dot(Z.T,Z)).dot(Z.T).dot(y)\n",
	"ypred = lambda Z,b: Z.dot(b)\n",
	"H = f_H(i)\n",
	"original_ypred = ypred(Z, b(Z,y))"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 9,
	"metadata": {},
	"outputs": [],
	"source": [
	"outliers,Dfs = [],[]\n",
	"t = 2*np.sqrt((p+1)/(N-p-1)) #cutoff value\n",
	"\n",
	"errs = y - original_ypred\n",
	"var_= sum(errs**2)/N #for the internally Studentized residuals\n",
	"\n",
	"for n in range(N):\n",
	" rstudent = errs[n]/(var_ * np.sqrt(1-H[n,n])) #the internally Studentized residuals\n",
	" \n",
	" Df = rstudent(np.sqrt(H[n,n]/ 1-H[n,n]*2))\n",
	" Dfs.append(abs(Df))\n",
	" if abs(Df) > t: \n",
	" outliers.append(n)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 10,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"0.408248290463863"
	]
	},
	"execution_count": 10,
	"metadata": {},
	"output_type": "execute_result"
	}
	],
	"source": [
	"t"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 11,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"[10]"
	]
	},
	"execution_count": 11,
	"metadata": {},
	"output_type": "execute_result"
	}
	],
	"source": [
	"outliers"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 13,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"[array([10.38312882]),\n",
	" array([0.13462557]),\n",
	" array([0.12596528]),\n",
	" array([0.11952226]),\n",
	" array([0.10717736])]"
	]
	},
	"execution_count": 13,
	"metadata": {},
	"output_type": "execute_result"
	}
	],
	"source": [
	"Dfs.sort(reverse=True);Dfs[:5]"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 17,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"<matplotlib.collections.PathCollection at 0x11e07cf28>"
	]
	},
	"execution_count": 17,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAFpCAYAAABJdYvCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3WlsXNd99/HfXebeGe7LiMtwkUTttC0pqu20QY1mUdEgCIwiL1zHSICmL4LCRdKicOs0aewXrgu2tZEAhYO8Sdoi3ZAUcYMCAYryCQI/bYIniR0lrmRL1i6K63C4c7a7PC8o06aphZoZce5wvh8jgHh5585/Tij9eM499xwjDMNQAAAgEsxqFwAAAN5BMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARYlfrjcfHx6v11pGTTCaVTqerXUbNod1KQ7uVhnYrDe22JpVKbflceswAAEQIwQwAQIQQzAAARAjBDABAhBDMAABECMEMAECEEMwAAEQIwQwAQIQQzAAAREhFVv76gz/4A8XjcZmmKcuyNDIyUonLAtipPE/mtWsyPU/ewIAUj1e7IiAyKrYk57PPPquWlpZKXQ7ADmVkMnL/+78VhqFkWbL/939VPH5c/t691S4NiASGsgFsK+enP1UYj0uJhOQ4ChsbFfvFL6RCodqlAZFQsR7z888/L0n6zd/8TZ08ebJSlwWwk2SzMlZXFTY2bjxuGDKnphQMDFSnLiBCjDAMw3Ivkslk1NHRoYWFBf3FX/yFPvOZz2h4eHjDOaOjoxodHZUkjYyMqMBvx+ts25bnedUuo+bQbqWparsVizK++13pvcG8sqLwN35D6umpTl1bwM9baWi3NY7jbPncigTzu337299WPB7Xo48+etvz2PbxHWyLVhrarTTVbrfYj38sM5ORYrG1A0EghaHyH/2oZBhVq+tOqt1utYp2W7Ot2z7mcjlls9n1P//yl7/U4OBguZcFsEMVH3pIQTIp4+1h7Xhchd/4jUiHMrCdyr7HvLCwoBdeeEGS5Pu+fv3Xf13Hjx8vuzAAO5Rtq/j+96t4o6csy6p2RUCklB3M3d3d+pu/+ZtK1AKgnpg8FALcDH8zAACIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiJCKBXMQBPrTP/1TjYyMVOqSAADUnYoF8/e//3319fVV6nIAANSligTz7OysXnvtNX3kIx+pxOUAAKhbFQnmv//7v9enPvUpGYZRicsBAFC37HIv8Oqrr6q1tVVDQ0M6ffr0Lc8bHR3V6OioJGlkZETJZLLct94xbNumPUpAu5WGdisN7VYa2u3uGWEYhuVc4J//+Z/1yiuvyLIsFQoFZbNZPfzww/r85z9/29eNj4+X87Y7SjKZVDqdrnYZNYd2Kw3tVhrarTS025pUKrXlc8vuMT/xxBN64oknJEmnT5/Wf/zHf9wxlAEAwM3xHDMAABFSdo/53e677z7dd999lbwkAAB1hR4zAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEWKXe4FCoaBnn31WnufJ93396q/+qh577LFK1AYAQN0pO5hjsZieffZZxeNxeZ6nZ555RsePH9fBgwcrUR8AAHWl7KFswzAUj8clSb7vy/d9GYZRdmEAANSjsnvMkhQEgZ5++mlNTk7qt37rt3TgwIFKXBYAgLpjhGEYVupiKysreuGFF/SZz3xGg4ODG743Ojqq0dFRSdLIyIgKhUKl3rbm2bYtz/OqXUbNod1KQ7uVhnYrDe22xnGcLZ9b0WCWpO985ztyXVePPvrobc8bHx+v5NvWtGQyqXQ6Xe0yag7tVhrarTS0W2lotzWpVGrL55Z9j3lxcVErKyuS1mZov/766+rr6yv3sgAA1KWy7zHPzc3ppZdeUhAECsNQv/Zrv6Zf+ZVfqURtAADUnbKDeffu3frrv/7rStQCAEDdY+UvAAAihGAGACBCKvIcMwAAtS4IAk3PzmhhaVlNDQ3q2dUly7K2vQ6CGQBQ93zf16k3z8jzPMVsW/OLCxqfntL7hu+XbW9vVDKUDQCoe2NTE/JvhLK0tg+EDENXxse2vRaCGQBQ95aWlzf1jC3T1NKNdTq2E8EMAKh7MTum9y6EGYbheg96OxHMAIC6N5jqU76Q33AsVyhoMLX9K1kSzACAupeIx/XAwSNrm274vkzL0n37D6q5sWnba2FWNgAAkpqbmnT00JFql0GPGQCAKCGYAQCIEIIZAIAIIZgBAIgQghkAgAghmAEAiBAelwIAlG1ucUFjE+Mqep4SiYSG+gfkOm61y6pJ9JgBAGXJzM/rjfPnVPQ8SdLKyopOvXFanu9XubLaRDADAMpybeK64m58/WvTNGXI0PjUVBWrql0EMwCgLAWvuOmYbdtaza1WoZraRzADAEoSBIHGp6eUzmQ0Oz+nIAjWv1coFtXW0nLT13m+r7HJCZ2/clmLy0vbVW7NYPIXAOCueb6vU2+clu/7am1u0YWrlzQ7P699g7sV+L5cx1FXR3LT61ZWV/XLs2/INE3ZlqWp2bS6Ozu1f/feKnyKaKLHDAC4a9fGrysMAsVsW3HX1aGh/WqMJ7SwvKj+nl4dOzws09wcMW9dvSwnFlPMtmUYhhKuq+nZtJZXVqrwKaKJHjMA4K4trSzLsqz1r23LVl9Pj2zLUl9P701fE4ahVrOrir/nMSon5ig9P6emxsZ7WnOtoMcMALhrsVhMYRhuOBaGoWKx2C1fYxiGbGtzf9APfMVdnnl+G8EMALhrAz0p5Qv5DcdyhbwGe/tu+7qe5C4Viu/M4g7DUIZhqKuj857UWYsIZgDAXWtqbNTwgUMyLUue78u0LN1/4JAaGxpu+7qB3pRSXd3yA1+FYkFuzNEDh47c9H50veIeMwCgJG3NLTp+ePiuXmMYhgZTfRpM3b5nXc/4FQUAgAghmAEAiJCyh7LT6bReeuklzc/PyzAMnTx5Uh/72McqURsAAHWn7GC2LEuf/vSnNTQ0pGw2qy984Qs6evSo+vv7K1EfAAB1peyh7Pb2dg0NDUmSEomE+vr6lMlkyi4MAIB6VNF7zNPT07p06ZL2799fycsCAFA3jPC9S7eUKJfL6dlnn9UnPvEJvf/979/0/dHRUY2OjkqSRkZGVCgUKvG2O4Jt2/JubDCOraPdSkO7lYZ2Kw3ttsZxnC2fW5Fg9jxPf/VXf6Vjx47p4x//+JZeMz4+Xu7b7hjJZFLpdLraZdQc2q00tFtpaLfS0G5rUqnUls8teyg7DEN9/etfV19f35ZDGQAA3FzZs7LPnj2rV155RYODg/qTP/kTSdInP/lJnThxouziAACoN2UH8+HDh/Xtb3+7ErUAAFD3WPkLAIAIIZgBAIgQghkAgAghmAEAiBCCGQCACCGYAaCOrWZXlZ6fU6FYrHYpuKHsx6UAALUnCAKdPn9OSyvLMmQoVKjUrm7t6R+odml1jx4zANShS9evKZvNKu64ch1HccfV9ekpLSwtVru0ukcwA0Adml9YlG1vHDSNO44mZqarVBHeRjADQB0yjVsdv8U3sG0IZgCoQ8mOzk0TvnKFvFJdPVWqCG9j8hcAVIjn+8oX8nIdV7ZlVbuc2+rv6VWhWNB0ZlZhEMi2be0b2KOmxsZql1b3CGYAqIDLY9c0kZ5WGAQyDFO9u7oiPcPZMAztG9yj3X0D8jxPTiwm02QQNQoIZgAo03QmrYmZKbmOu35sYmZKDQ0JdXUkq1jZndmWFfnefb3h1yMAKNNUOr0hlCXJdVxNpdNVqgi1jGAGACBCCGYAKNOujk7lC4UNx/KFgnZ1dFSpItQy7jEDQJm6O5NaWl7WzNyswjCUYRja1d6p7s5d1S4NNYhgBoAyGYahA3v2andfv7K5rMJQml2Y0+Xr19TX3SMn5lS7RNQQghkAKsSJxTSdTuvK+DW5jqtQ0sT0lA7u3adkO8Pa2BruMQNAhRS9oq5NjisRT8g0TVmmqbgb16VrVxWGYbXLQ40gmAGgQhaWl296vOB5myaHAbdCMANAhcRjjoKb9IxNw9i0kxNwKwQzANylMAy1mssql89vON7U2KjGeFx+EKwf83xPnW1trK6FLeNXOAB1byYzq8n0jMIwVLK9Q727umTcYvvD5ZUVvXHxvArFgiRDDfG4hvcfWF/56/5DR/TWlUtaXFqSaRrqbO/Q3r7orpmN6CGYAdS1q+PXdX1qUq7j3Ph6TEsryzq0d9+mc4Mg0Jnz52TbthJufP3YGxfO6/iR+yStrT19ZGj/9n0A7DgMZQOoW77va3x6aj2UJcmJOZqdy2wappakxeVleb6/4ZhhGFrNrt7oQQPlI5gB1K1CsbjhfvA7DK1mszc5Hko3RrjDMJTneQrDUKEM8TQUKoWhbAB1y4nFZN10D+JQjQ0Nm462NDXLNC3NZGaVWVhQEASyTFNdyaScWExXxseUzswpCH21NjVraHAPk75w1yoSzF/72tf02muvqbW1VS+++GIlLgkA95xlWUp1dWtsakLxG5O38oWCDJl648I5eX6glsYmDQ3uliSZpqld7W1669JFmeZa19m0TBmGodfPndVqdkVBECoIAs2HS3r97Bt63/D9Vft8qE0VCeYPfvCD+uhHP6qXXnqpEpcDgG2xms2q6HtyHVdBGMp1HBmScoWCgiCUaRhaWF7SL8+eUXdXlyRpJZvT8P4DyhfykmEo7rhrk8LeOivTNFX0ipIMWaapzvYOLS4vqaWpuaqfE7WlIveYh4eH1dTUVIlLAcC2mJiZ1s/P/K/m5ucVBoHy+bw629pV8P0Nk8Es01Q+X9DsXEbS2oQxwzAUd+PrvewgDDQ1m5a0NnnMicVkWZYmZ6a1mrvZvWrg1pj8BaDuBEGgK9fHlIjHZZqmTNNUIh7X5bFr8r3ipvMty9LSyookqb21VUXP2/D9XD6vhOtuevbZD3yZBv/M4u5s2+Sv0dFRjY6OSpJGRkaUTCa3660jz7Zt2qMEtFtpaDdpaXlZDQ0NakgkNhy3s1lJoRoSGyd+reZy6tnVpUQ8rvaODp06/boWl5dkWbZ839Oe5G4ZpqmZzKxisZgMGSoUC0r19Kinu0cdbW3b+OmihZ+3u7dtwXzy5EmdPHly/et0Or1dbx15yWSS9igB7VYa2m1tgtfy8pK84sbecS6f09DAbl28dlVOLCbTNJUrrA1xJ+Lx9Xbb3dunpeVlrWRX1NLUrIZEgzKZjGzT0tzignzfV1dHpxw7Jq9QqOv25udtTSqV2vK5PC4FoO64jqPW5hatrK6uby7heZ5am1vUs6tLrS0tGp+aVKHoabCvX52tm3u8zU1Nan7X3Jrh/Qd0+q1zamtuVihDTiymI0P7Zd70cSzg1ioSzF/96ld15swZLS0t6fd///f12GOP6cMf/nAlLg0AZcnl81rNZtXU2CAn9s6kriP7DujitSuaW1yQJHW0tWloYO2xqIQb177BPXf1Pk7M0fuG71cun1cQBkq48Vuutw3cTkWC+Y/+6I8qcRkAqJggCHT20gVlFhbeXqxLXZ2d2r97r6S1Z5Lf/nMlxV234tdEfWEoG8CONDY5qcXlZSXeFZQzmYzamluU7OisYmXA7XHzA8COlFmYU8ze2PdwHUdTTERCxBHMAHYk8yb3d8MwlGFy3xfRxlA2gMjzfX9tkpZhqKOldX2mcxAEujY5rszcnCSps71DA70pGYahrmRSl8auyX3XhK9cobA+wQuIKoIZQKTNLS7o7MXz69sqGqah4aEDamlu1rnLl7S4tLj+yNPEzLRy+ZwO7t2nnmSXcrm8JtMz8oNAtmVqd2+f2ltbq/hpgDsjmAFEVhAEOnfpomJ2TIViUbZlybIsnb10QYeG9un1s2cUhKFMw1RHW5vaW1qVnp/TnmJRTiymPf0DGuhNqVAsynUcnilGTSCYAUTW8uqqZucyWlheVhD4MmSoualRrU0t+n+/+LnyhYLcGxtJTKVnFAaBEvGEijeCWVpb5zrBnsioIfz6CCCyCsWCpmbTsi1LTsxRLBbTSnZVVybG5DiOTPOdwI3ZMWUWFmRZFs8So6YRzADuuTAMtbSyrIWlRQVBsOXXzS8uqCGRUKhw/ZhpmMrmcmpKNKizrU2FYlHhjf+y+az6untk0UNGDWMoG8A9tZrN6vT5cyoUi5JCWZatw3v3qa2l5Y6vDYJQg719Gp+eVD5fkAwp7rjav3uv8oWCdnV0qjHRoNmFecmQ+ntSGujd+mYBQBQRzADuqTcvviXTMDaswHX20gU99MCxO07G6t3VpZm5jPb2D8r3fUlSKKmjtVWFoqeFpUU1JBJyHUee7+uBg4fv5UcBtgXBDOCeyeXzyuXzirvxDcd939fi8vIde83NTU0a6O3V9clJheHagHZLU5OGBnbLMAxlFhaUzszKdR2lunrWJ3wBtYxgBnDP3Gp3pTAMtdWNlwZ7+5Ta1a3l1RW5jqtE/J2Q72xrU2fb5i0ZgVpGMAO4Z1zHUUOiQb7vr4d0GIaKxWJqaWre8nVs21ZbCwuDoD4wKxvAPTW8/6BMw9D84qIWV5alMNR9+w+yVzFwCwQzgHtqbnFe+WJhbaJXKDmuy3PGwG0QzAC2JAgCjU9P6c2L5zU2ObE+S/p28oWCzl++LCfmqKWpSS1NTcrlcrpw9co2VAzUJu4xA1gXBIFWc1nZlr2hV+v7vk69eUae5ylm21pYWtLE9JSOD9+nmH3rmdDj05NyHWfDMcuy1naKAnBTBDMASdLkzIwuX78q3w8kw1BLU5OG9x2QZVlrPeQboSxJMduWHwS6cn1M+3fvrXLlwM7CUDYA5fJ5Xbh2RU7MUSIeV8J1lX3XkPPSyvL61opvs0xTy6urt71uqqtHhWJhwzHf99XODGvglghmAJqYnpL7nsU5bMvS/NLakHPMjikMww3fD8Pwjgt6uI6jfbv3qFAsKpvLKZvPKx6Pa9/g7sp+AGAHYSgbwG35vi8Z0vmrl9QQb9Cujk45sZhyhYIO7h264+u7O3dpV3vnTe9dA9iMHjMA9XZ1K18sbjjm+b5aGpt16s0zWlhcVF93r4rFos5fvqh8oaD79h9Uc2PTlq5vmqaaGhoJZWAL6DEDUNx1tX9wt85evKCF5SXF3bh6d+1Sc2OjFpYXFbNjsm1be/oHJL29Etedd4cCcPcIZgAKw1DzS4syTFOtN5bKLHqe5leWNz0OVSwWdWXiuoIgkOs4GuztU2NDQzXKBnYkhrIBaHZhXrPzc2pMJNTU2KimxkZ5nqfFxQUV3zXE7XmeLoxdle/5a888Z7P6xdkzWrnD7GwAW0cwA9DMbFpxZ+P9X8uyZFsxmaYpPwjWzpuble/56u5MSlrbPcqNOboyPrbtNQM7FcEMQJZlK7gRvu9mxyy9b/h+dbS2yjRNhZKGBgc3DF0bhqFCobDptQBKwz1mAOrv7tHPM7NKvGvWdKFQ0O6+ftm2rX2DeyRJbS0tmp6d3fDaMAzlxhPbWS6wo1Wkx3zq1Cn94R/+oT73uc/p3//93ytxSQDbqCGR0KG9QwrCULlCXp7vKdXdo96u7g3n9Xf3StL6BhZBEKhQLGpPf/+21wzsVGX3mIMg0De+8Q39+Z//uTo7O/Vnf/ZnevDBB9XPX1SgpiTbO5Rs75Dn+zINY22bxvewbVsnhu/XtYlxrWRXFXddDab65MScm1wRQCnKDubz58+rp6dH3d1rv1l/4AMf0E9/+lOCGahRtmXd/vu2rb0Dg9tUDVB/yh7KzmQy6uzsXP+6s7NTmUym3MsCAFCXyu4xv3dhe2ltluZ7jY6OanR0VJI0MjKiZDJZ7lvvGLZt0x4loN1KQ7uVhnYrDe1298oO5s7OTs2+a5bm7Oys2tvbN5138uRJnTx5cv3rdDpd7lvvGMlkkvYowU5rN9/3tbSyItu21NTQeM/eZ6e123ah3UpDu61JpVJbPrfsYN63b58mJiY0PT2tjo4O/ehHP9LnP//5ci8L1JXp2VlduHpZQRhIMhR3XT1w8PAdt1UEsPOUHcyWZen3fu/39PzzzysIAn3oQx/SwMBAJWoD6kLRK+r81UsbVt4KgkDnLl7Q/YcOV7EyANVQkQVGTpw4oRMnTlTiUkDdmclkZJkbZ0Kbpqml1RUFQXDTx5ZuJwgCZRbmFYahOlrbZN1hljWAaGHlL6DKTEPSzSZRavMkyjtZXFrSmQtvSQrXL3lo7z51tLVtOC8zP6drE+Mqep7i8biG+gfVkGD1LiAKWCsbqLJkR1LBe4LZDwK1trTcVW85DEOdvXRBTiwmJ+bIddb+d+7yxQ3rYM9kZvXmxfPyb/TG8/m8fvHmGRW94m2uDmC7EMxAldmWpSP7Dsj3feXyOeXyOcVdVwf27N1w3sLSkq6Mj60PU79XNp9TwfM2HfcDX8urK+tfXx67prgbX//aMAzZlqXrU5MV/FQASsVQNhABbS0tevCBY8rl87Isa8Ns7CAIdObCW1pYWpITi2l8akpx19XRQ0c23D9eu0+9ObCljctrejcJb8uylM3lKvmRAJSIHjMQEYZhKBGPb3pEanp2RssrK0q4rizTlOs4KnqeLr9nD2TXcdTc0Lhh2DoMQ8VdZ8Nz0YlEYlOPu1Asqr2l9R58KgB3i2AGIm52bn49rH3f1/j0lK5NXNfrb76hhaXFDece2X9AccdVrpBXrpBXzLZ1/4FDG845tHdIRc+Td2OHqEKxqLjrqquT1ZmAKGAoG4g4y7YUFkKFYaiL165KhmQaphQGOv3WOR3ed0AdrWu93Zgd0/2HDq+FbhjKtjf/FU/EE3rw/qO6PjWpbC6nVHePujo67/qxLAD3BsEMRFx/d69+efaMVrJZhWEgy7RV9IpKdfUo7rq6Oj62Hsxv28oOUbv72AEOiCJ+RQYirqmxUQf37lM2l1MQhgrDQN0dSbU0NUlaG4oGsHPQYwZqQLK9Qw8fe58uXbsq13E2fC/uurd4FYBaRI8ZqBFdHZ2Ku2szsqW1x6hyhbz29rM2PbCT0GMGaoRpmjp6aFjXp6a0uLwoJxbTQG+fEvH4nV8MoGYQzEANsSxLg6mUpK3v7QqgtjCUDQBAhBDMAABECMEMAECEEMxAGQrFolZz2Q3rUwNAOZj8BZQgCAKdvXRRc4trWzDatq2h/kHt6uisdmkAahw9ZqAEl8auaWl5SXHHVcKNK2bZOnf5ovKFQrVLA1DjCGagBLPzc5s2iHBijsanJ6tUEYCdgmAGShJuOmJICsLNxwHgbhDMqAuFYkGT6WnNLazdEy5Xa3OL/PdM+MoXC0rt6i772gDqG5O/sONdnbiusckJWYYpPwjkxGI6eviInJhz5xffwr7BPTrz1lktr64o1Nr+yHv7BlkeE0DZCGbsaNlcTtcmxpVw1wIzJikMQ527eFH3Hzpc8nVty9LRw8NazWVVKBTV1Nh4xz2QAWArCGbsaFOz6U09Y8MwtLS6UpHrN8QTaognKnItAJAIZuxwMctSEASyzI3TKUzz5tMrrk9Nano2LT8I1NLYpKHB3fSEAWwrJn+h5oVhqMmZaZ29eF7j01PyfX/9e927ujatyuX5npLtHZuuc3X8uq5OXFcYhjINQwvLS3r97JmKTBYDgK2ix4yatrC0pP/zo/+rxdVluTFHcddVqqtb7xt+QE4sJtuydN/+g3rr6mXl83mZhqnO9nbt7R/YcJ0wDDUxMy33XcPelmkqm89rfnFR7a2t2/3RANQpghk1y/M8/ejnP9NKLquV1azmCguSpKXlFbU0NevIvgOSpJbmZv3KfQ8oXyjI9325jrNpKNsPAvmBr9h7/krYpqXV7CrBDGDbEMyoWRMz01rJZpWZn5dhaP1ecHpuVheuXV0PZkkam5zQ9alJ+YEv0zDVs2uX9vS902u2Leumj08VfY9QBrCtyrrH/OMf/1h//Md/rN/5nd/RhQsXKlUTsCWFYlGe78n3fRky1o8bpqmVlaX1e8tzCwu6On5dMdtW3HHlxGKanJ7WVHpmw/WG+geUzefWX5fP59XV0amGRMP2fSgAda+sYB4YGNBTTz2lI0eOVKoeYMv6ursVj8U2LIP59sSttua29eMTM1OKu+6G1zqOo6nZ9IZjHW3tOjH8gNpbW5VIJHRwaL8O7Bm69x8EAN6lrKHs/v7+StUB3LW4G9fRQ0c0PTurbD4nz/Nk2zEdGdqvZGfn+tD2rfZKvtls60Q8rn2De+5l2QBwW9t2j3l0dFSjo6OSpJGRESWTye1668izbZv2uInp2VmNTVyX5/vqbO/Qnv6BDc8j27ath088qJaWVp0+f1ZNiUYlEnEFQajj9z2gtuZmSdKh4IDevHh+ffUvaW0YfHeqvy7bnZ+30tBupaHd7t4dg/m5557T/Pz8puOPP/64HnrooS2/0cmTJ3Xy5Mn1r9Pp9G3Ori/JZJL2eI/xqSldHh9T3FmbkDWTTuvK1Ss6enh4/Zy32y3Z1q6H7juqyXRaTiymnuQuefm80vm8JMk2TMXtmKZmpmUapsIwUHtLmxKOU5ftzs9baWi30tBua1Kp1JbPvWMwf/nLXy6rGOBuhWGosamJ9VCWpJhta3l1VQtLS2q90RN+t4ZEg4YGBm96PcMwdGDPkAZT/VpaXVFjPMFmEwAii8elUBVBECg9N6vM/ILicVeprh45sZiktWeKPd/btBSmbVlaXlm+aTBvhes4cp3Sd5QCgO1QVjD/5Cc/0Te/+U0tLi5qZGREe/bs0Ze+9KVK1YYdKgxN22O8AAAMoUlEQVRDnT53VsvZVbmOo6WVZY1PT+vYoSNqbGiQZZqK2Zt/NIu+r5YSQxkAakVZwfzwww/r4YcfrlQtqBOzcxnNzGeUzWVlmbY62lrlxmK6NHZV9x88LMMwNJjq0/krlxV3XBmGoUKxqLaWFjU3NlW7fAC4pxjKxrY7ff6cJqanFIvFFIahMgtz6u9JyX/XY03dnbvU4CY0NjUhz/eV6upSd7KrilUDwPYgmLGt8oWCFpeXZJmWDBkyDEOmbWpyZkoH9m5czKO5qUlHmg7c4koAsDOx7SO21dzigjrbOhSG4YYFPrK5nLo7dlWxMgCIBoIZd201u6rM/LyKXvGuX5u4sTTm3oFBuY6jIAxkGoZSXT3q2UUwAwBD2dgy3/d15vxbWlpZkmRIhtTf3avBVN+Wr9HS1CzXdRUGgQZ61x64L3qeku3tm7ZiBIB6xL+E2LJLY1eVzecUd+OKu67ijqtrkxNaWl7e8jUMw9DRQ4fV2NC4tjNU4Ks7mdTQwO57WDkA1A56zNiy+aXFTYt+xB1HEzPTam7a+mNMMTumI/v2V7o8ANgRCGZs2Wo2q5nMrIpFT7ZtKdneqaaGBpmmcecXS1pZXdXVievKFwpqiMe1u69fruPe+YUAUEcYysaWrGazmltYUNEryrIshaF0fWpS8wuLSnX13PH1K6ur+sWbZ7SazSoIAi0uL+vnb5xWoXj3E8gAYCcjmLEl1ybG1dfdo6ZEo4peUYViQZZpyo07akgk7vj6K+Njch1HhrHWuzZNU5Zp6drk+L0uHQBqCkPZ2JKiV5Rpmkp196jb9+UHvmJ2TNZ77jnfSqFQWA/lt1mmqVwudy/KBYCaRY8ZW9LW0qrijWFny7LkxBwFYajmxsYtvd514xsWFJEkz/fV2NBQ8VoBoJYRzNiSVFe3XNdVvlCQtNaDNgxDu1P9W3r9nv5+FYpFBTfWw/Z9X1Ko/u7ee1UyANQkhrKxSRAEunz9mjILC5KkztY27e7r19FDRzQ7P6f5xQU1NzWrq6Nzy4uCJNy4Ttz3gK6OX1e+kFdDS4sGe/tk32R7RwCoZ/yriE3evHheSysr63siz8xllMvndWT/Ae3q6NSujs6Srus6jg7s2VvJUgFgx2EoGxvkCwXNLS6sh7Ik2ZalucV5FYqFKlYGAPWBHjM2KHpFKdx8PAylQrG4NukrCOT5vmK2vWmm9c3k8nldnRhXPp9TQ0ODBntTitmxe1A9ANQ+ghkbJNy4LHPzI1CWZSrhxnV14rompqfl3wjmwb4+dXfeeleobC6rU2+cVsyOyTRNZebmNDuX0YnhB7i/DAA3wVA2NrAsSwOplLL53Pqeydl8TgO9fUrPZXR9clIx21bcdWVZls5fuaKllVtvYnH5+picmLM+SWztuWdDY1MT2/SJAKC2EMzYpK+7R8cODaulqVktTc06dmhYfd09mkrPyHWcDefGHUdjk5O3vFYun9803G1bllZWV+9J7QBQ6xhLxE01NTbqQOPGGdT+jWeQ380wDAWBf8vruI6jbC63IZx931ciHq9csQCwg9Bjxpa1t7bK87wNx3KFvHZ1dNzyNbtT/coXCuurfgVBID8MNNCbuqe1AkCtIpixZYO9fUq4cWXzeRU9T7l8Xp1t7drVkbzlaxobGnTs8LAaEg0yTUstTc06MXw/s7IB4BYYyq4Tnu8rnUkrXyyqu3OX4u7d74NsmqYeOHxEi8tLWs2uqqWpWQ2JO6913djQoCP79pdSNgDUHYK5DiyvrOj1c2/KNAyZpqmxiQntTvWrv7e0darfnhQGAKg8hrLrwLkrF+U6jmKxtW0aE/G4rk5eZyUvAIgggnmH827cC34vy7SUmZ+rQkUAgNshmHc40zRlGpv/b/YDX7HY3d9nBgDcWwTzDmeapjra2uX57zzmFIahYpat9paWKlYGALiZsiZ/fetb39Krr74q27bV3d2tJ598Uo2NjZWqDRWyf3C3Lo0ZmslkFISBmhoadXDP0Jb3UgYAbJ+ygvno0aN64oknZFmW/vEf/1Evv/yyPvWpT1WqNlSIaZraN7hHQwO7JWlLO0IBAKqjrC7TsWPHbmxKIB08eFCZTKYiReHeMAyDUAaAiKvYWOYPfvADHT9+vFKXAwCgLt1xKPu5557T/Pz8puOPP/64HnroIUnSd7/7XVmWpUceeeSW1xkdHdXo6KgkaWRkRMnkrZdxrDe2bdMeJaDdSkO7lYZ2Kw3tdveM8O3dBUr0wx/+UP/1X/+lZ555Ru5dLPM4Pj5eztvuKMlkUul0utpl1BzarTS0W2lot9LQbmtSqa1v3FPWUPapU6f0ve99T08//fRdhTIAALi5smZlf+Mb35DneXruueckSQcOHNBnP/vZihQGAEA9KiuY//Zv/7ZSdQAAALHyFwAAkUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEUIwAwAQIQQzAAARQjADABAhBDMAABFCMAMAECEEMwAAEWJXu4ByhGGomblZzc7NyYnF1Nfdq7jrVrssAABKVtPB/Mb5t7SwvCTXcbSyuqqp2bSG9x9UW3NLtUsDAKAkNTuUvbC0pPmlRbmOI0kyTVNxx9Wla1erXBkAAKWr2WDOzK8NX79XNp+rQjUAAFRGzQZzY2ODPN/bdNyxN4c1AAC1omaDOdnWIduyFQTB+rF8saBUV08VqwIAoDw1G8ymaerYkWG1NDUpCEMZpqmhgd1KdXdXuzQAAEpW07OyY3ZMB/fuq3YZAABUTM32mAEA2IkIZgAAIqSsoex//dd/1c9+9jMZhqHW1lY9+eST6ujoqFRtAADUnbKC+dFHH9Xjjz8uSfr+97+vf/u3f9NnP/vZihQGAEA9Kmsou6GhYf3P+XxehmGUXRAAAPWs7FnZ//Iv/6JXXnlFDQ0NevbZZytREwAAdcsIwzC83QnPPfec5ufnNx1//PHH9dBDD61//fLLL6tYLOqxxx676XVGR0c1OjoqSRoZGVGhUCin7h3Ftm153uZVzHB7tFtpaLfS0G6lod3WODf2ddiKOwbzVs3MzGhkZEQvvvjils4fHx+vxNvuCMlkUul0utpl1BzarTS0W2lot9LQbmtSqdSWzy3rHvPExMT6n3/2s5/d1RsDAIDNyrrH/E//9E+amJiQYRhKJpPMyAYAoExlBfNTTz1VqToAAIAqeI8ZAACUjyU5I+ALX/hCtUuoSbRbaWi30tBupaHd7h7BDABAhBDMAABECMEcASdPnqx2CTWJdisN7VYa2q00tNvdY/IXAAARQo8ZAIAIKXsTC1TGt771Lb366quybVvd3d168skn1djYWO2yIu/HP/6xvvOd7+j69ev6y7/8S+3bt6/aJUXaqVOn9Hd/93cKgkAf+chH9Nu//dvVLinyvva1r+m1115Ta2vrlpcchpROp/XSSy9pfn5ehmHo5MmT+tjHPlbtsmoCPeaIOHr0qF588UW98MIL6u3t1csvv1ztkmrCwMCAnnrqKR05cqTapUReEAT6xje+oS9+8Yv6yle+ov/5n//R2NhYtcuKvA9+8IP64he/WO0yao5lWfr0pz+tr3zlK3r++ef1n//5n/y8bRHBHBHHjh2TZVmSpIMHDyqTyVS5otrQ39/PGu1bdP78efX09Ki7u1u2besDH/iAfvrTn1a7rMgbHh5WU1NTtcuoOe3t7RoaGpIkJRIJ9fX18e/aFhHMEfSDH/xAx48fr3YZ2GEymYw6OzvXv+7s7OQfSmyL6elpXbp0Sfv37692KTWBe8zbaCt7W3/3u9+VZVl65JFHtru8yNrqnuC4vZs9gGEYRhUqQT3J5XJ68cUX9bu/+7tqaGiodjk1gWDeRl/+8pdv+/0f/vCHevXVV/XMM8/wD+a73KndsDWdnZ2anZ1d/3p2dlbt7e1VrAg7ned5evHFF/XII4/o/e9/f7XLqRkMZUfEqVOn9L3vfU9PP/20XNetdjnYgfbt26eJiQlNT0/L8zz96Ec/0oMPPljtsrBDhWGor3/96+rr69PHP/7xapdTU1hgJCI+97nPyfO89UkmBw4cYH/rLfjJT36ib37zm1pcXFRjY6P27NmjL33pS9UuK7Jee+01/cM//IOCINCHPvQhfeITn6h2SZH31a9+VWfOnNHS0pJaW1v12GOP6cMf/nC1y4q8N998U88884wGBwfXRwA/+clP6sSJE1WuLPoIZgAAIoShbAAAIoRgBgAgQghmAAAihGAGACBCCGYAACKEYAYAIEIIZgAAIoRgBgAgQv4//l/4g58O13oAAAAASUVORK5CYII=\n",
	"text/plain": [
	"<Figure size 576x432 with 1 Axes>"
	]
	},
	"metadata": {},
	"output_type": "display_data"
	}
	],
	"source": [
	"c=['#415952', '#f35134', '#243AB5']\n",
	"carray = np.array([c[0]]*N)\n",
	"carray[outliers]='red'\n",
	"fig, ax = plt.subplots(figsize=(8,6))\n",
	"ax.scatter(i[:,0], i[:,1], label=r'$df=%i$' % 2, c=carray, alpha=0.3)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {
	"collapsed": true
	},
	"outputs": [],
	"source": []
	}
	],
	"metadata": {
	"hide_input": false,
	"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.5.6"
	},
	"toc": {
	"base_numbering": 1,
	"nav_menu": {},
	"number_sections": true,
	"sideBar": true,
	"skip_h1_title": false,
	"title_cell": "Table of Contents",
	"title_sidebar": "Contents",
	"toc_cell": false,
	"toc_position": {},
	"toc_section_display": true,
	"toc_window_display": false
	},
	"varInspector": {
	"cols": {
	"lenName": 16,
	"lenType": 16,
	"lenVar": 40
	},
	"kernels_config": {
	"python": {
	"delete_cmd_postfix": "",
	"delete_cmd_prefix": "del ",
	"library": "var_list.py",
	"varRefreshCmd": "print(var_dic_list())"
	},
	"r": {
	"delete_cmd_postfix": ") ",
	"delete_cmd_prefix": "rm(",
	"library": "var_list.r",
	"varRefreshCmd": "cat(var_dic_list()) "
	}
	},
	"types_to_exclude": [
	"module",
	"function",
	"builtin_function_or_method",
	"instance",
	"_Feature"
	],
	"window_display": false
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}
No results found