Illustration.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from robust_mean import M_estimator\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist, euclidean\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "from neurtu import timeit, delayed  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.11174249649047852\n"
     ]
    }
   ],
   "source": [
    "a = time.time()\n",
    "X = np.random.normal(size=[1000,1000])\n",
    "result = M_estimator(X, beta=1)\n",
    "print(time.time()-a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0155168291969072\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.norm(result)) # distance to the mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D Illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 60 # Sample size\n",
    "name = 'Huber'\n",
    "epsilon = 0.1\n",
    "center = [20, -1]\n",
    "\n",
    "np.random.seed(2020)\n",
    "\n",
    "# Sample dataset from a skewed distribution.\n",
    "X = np.vstack([np.random.normal(size= [3*N//4,2]), \n",
    "              3+np.random.normal(size= [(N//4),2])])\n",
    "\n",
    "# Contaminate with outliers\n",
    "for f in np.arange(int(len(X)*epsilon)):\n",
    "    X[f] = np.random.normal(size=2)+np.array(center)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAE/CAYAAACdNPbaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABFGUlEQVR4nO3deZwcdZ3/8denJycwQzImHJMAARxFBAQ3wqKogCMC4kJUvFdls2Z1PR+LuIrsCuiq+3NXdF1Fo0FABRdFRPGCkUtAwIDcIMNpmAAJTCYzIQkkM5/fH99qpqbT50x3V3XN+zmPfvRUdR3fqq6uqk99L3N3REREREREWkku6QSIiIiIiIjUSoGMiIiIiIi0HAUyIiIiIiLSchTIiIiIiIhIy1EgIyIiIiIiLUeBjIiIiIiItJwpHciYmZvZC5uwHjOz75vZOjO7uYrpF0Vpm9botNWbmb3BzH4eG95gZntNYDll94GZnW5mP5xEUvPL2dnM7jWzmZNdVpXrO9zMHmvCek41s+81ej0TUa/vLglmdr2ZHRT9n9g+LjyOzOxmM3tpEmlplkadr83sbjM7vE7LerWZ/aUey0obM3vEzHqSTkcpZrZ7dL1pSzothep1Ta90PZ3sd2Rm55rZFyY6fxqZ2W/M7H1NWlfLXttaWcVAxswOM7MbzGy9mQ1EF/JXTGalZvZ+M7uuYFxqf0DF0lujw4DXAwvd/eA6JQto6o1xtev5D+DL+QF338HdH2pcyibH3Z8ErgKWJZ2WiSr23bj7F939Hxu0vlTf0DSKmb0JGHb3P0Nj9/EE/BdwZi0zmNkrzexKMxuOzu+/NLN9G5S+1HL3l7r71dVMW+nYd/c/uPuL65Y4Kanwu3D3v0bXm5EGrCsVN6jx6+lk75nqcF/TEtz9GHc/r5ppG3lti67To1EwusHM+s3sjBrmT+09ctLKBjJm1gFcBnwD6AQWAGcAzzY+abVJee7FHsAj7v5M0glppCjA3dHdb0w6LTX6EfBPSSdCUu+DwA8aseA6PEX+BXCEme1S5foOBS4HLgW6gD2B24HrJ5KDWquUn69FRBphdRSM7kB4wL3UzE5IOE2tz91LvoDFwGCFaT4A3AsMA/cAL4/Gfxp4MDZ+STT+JcBmYATYAAwSnoZvAZ6Lxv0ymrYLuBhYCzwMfCy23tOBnwI/BIaAfyyStnOBbwNXROm4Btgj9rkDL4z+3xE4P1rXo8BphEBvm/SW2A9dhJuJAeAB4APR+KUF859RZN42whPVp4CHgA9HaZsWfX5SbB8/BPxTNH57YBMwGi17Q5SOg4E/Rvv2ceB/gRnRPAacBayJ9tudwH7RZzOjdPwVeDLad7NLrafIdvw78L2CcfF9fC7wTeBX0bbcBOxdYn8uiuZdBqyOtuOTBd//D6P/DwceK5j/EaAn+j/H2PH4NHAR0BmbdhqwMX5sFCzrjcCfo/21Cji9SDrfF+23p4DPxj6fHW33OsLv4JTCtBasax/C8ToA/AV4W+yzY6NlDAP9wCfLHAPx/ZNP40lR+tcRbspfAdwRHSf/G1vP3sCV0b56ihDozYk++0G0rk3Ruj4Vjf9b4IZoWbcDh8eWtyfhtzccbdv/5tNWZPsPBx4DPkU4Rh8HToi2/f5ov5wam77Sd/sT4AlgPXAt8NKC80O1x+OMaJsXljgGyx4HJc5NZwO/Bp4Beih/vqt4HEX79n3lztexaf8AfKvI+N8A5xd8FyfHvouTYtMWPV+UWN/7gesJ556ngS9Umj/axscJv/9/IHYuKbL8q4EvATcTfqeXFhwHfwfcTTg+rwZeUuJccXp0DJ0fHRN3A4vLHfvFjt/Y8L8SfqvDhN/z60qkv+j1J7bvrov21bro2DimYN4V0b7qj/ZtW5F1dEVpj++XgwjH6nTK/O6L7KdzgS+U2e6Sx3KRdJU8DoB5hIepg4Tf/h8Iv/ltvgvGfoP56+bV0b64IZrml8ALou0aAv4ELIql4+uE8+MQcAvw6mj80YR7ky3Rcm6vtN+pcE0v2P6TiO55ouE+4Cex4VXAgdH/DryQ0vdMjxCuC3cQznn/B8wqss6i9zVUOCdS5vpUYru2uW8p972W+81Ex8nXCOeD1dH/M2PLPB64Lfr+HgSOjh0H/xj9n/i1rWDcRYy/nhXdv2W+76L32VPtVf5D6Ii+8POAY4C5BZ+fGB1wryDcIL+Q6GYw+qyLcNJ5O+FivWv02fuB6wqWdS7jT4w5wsnk3wk3EXtFP4Y3RJ+fHn2xJ0TTbnMBjZY5DLwm+hF8Pb5ext9kn0+4+LUTToj3A0tLpbfIuq4FvgXMAg4knMCPrGZ+wk3lfcBuhJyvqxh/Qn4j4QdowGsJN935gLHYj+NvCD++adG23At8IvrsDdF+nRMt7yWx7+UsQjDWGe2HXwJfKrWeItvxE+CUgnHxfXwu4Xg6OErbj4Afl1jWomjeCwk36/tH+zR+w1FtIPNx4EZgYXQcfAe4sGD6O4C/K5GWw6P154ADCBfbEwrS+V3CzebLCDmWL4k+/zLhJN0Zfb93ldqP0XauIlwApjF2k7Fv9PnjjF1c51Y4BuL7J5/GbxOOz6MIF7GfAzsRclrXAK+Npn8hoSjkTGA+4dj+WrF9Gw0viL7XY6N99PpoeH70+R+Br0bLew3hN1nuZL+V8LufTnhQsha4gHBMvpRwodmzmu+WcAPczthF8LaC80O1x+NLgWeq2MdFj4MS56b1wKuifbYd5c93FY8j4H+Ar5b7jUbTbUe4gTmiyGcnAY8XfBdnRt/FsYRzz9xK54siy31/tKyPRvt6drn5CTePTwL7EX4XF1A5kOmPTX9x7Lt5EeH68/poOz5FeNiUf7jzCOPPK5ujbW0jBEc3ljr2Sxy/j0X/v5jwe+6KHSOlAuVK158thN9CG/Ahwo2cRZ9fQjjutyf8nm8mdtNYsJ4riR6yRcNfAb5d6++eMoEMFa7dRdJU7jj4EuG8NT16vTq23eO+C4oHMg8Qrp07Em707ic8NJgW7fPvx+Z/DyHQmUYI3p8gCgKI/dZj05fc71S4phcsZy/CTXKOcM/0aGxf7kUIXvM3+YXX0y8ULOuRKB1d0XrvBT5Y5jdZ7D6s6DmRCtenIssvd99S9HulzG+GcB66MdrX8wnBxeejzw4mnE9fH+3HBcA+seMgH8gkfW2LB/vdhHPWkdXs3xLfd8n77Kn0qjxBuNE9l/BkbivhhLNz9NnvgI9XtaIQKR9f4QcUPzEeAvy1YJrPEJ14CCeWayus81xiNybADoQL+G7RsEcHdhsh0t03Nu0/AVeXSm/BenaLltseG/cl4Nwq57+S2MmGcKNZ9KQXff7z/H4v/HGUmP4TwCXR/0cSTuZ/S3RyjMZb9COIP305FHi4hvVcQcFJk21PvN+LfXYscF+JZS2K5t0nNu7/ASti33+1gcy9xJ6EArsSbgymxcZdD7y3ymP5a8BZBemMP6m/GXhH9P9DRE+GouFlpfYj4UT0h4Jx3wE+F/3/1+i47CiYptj2x/dPPo0LYp8/Dbw9NnwxUbBbJF0nAH8utm+j4X8FflAwz+8IuRO7E84b28c+u4DyJ/tNjD3ZbI/SfkhsmlsYCyQrfrexz+ZEy9pxAsfjq4AnqtjHRY+DIss7lyjnIxqudL6reBwR6qedU8Xxu5CC31bss6OBLQXfRfx3soZw7ih7viiy3PfHt6/S/MA5wJdjn72IyoFMfPp9Cef0NuDfgItin+UINxCHFx7P0XfaW7CcTaWO/RLHb/4m9IXR/uoBppeZp5rrzwOxz7aL9sUuwM6EgDmek/VO4KoS6/pH4MrYd7AKeE2JaU+gxO+e8oFM2WO5YHyl4+BMQoC3zfde+F1QPJCJ547/N/Cb2PCbiD3YKLL8dcDLYsfFD2Ofld3v1H5NXwW8HHgHsJxw7tiHcFP7i9h0hdfTYoHMe2LD/48oUC3xmyx2H1b0nEiF61OlF+PvW4p+r5T5zRByHo6NDb+BUGQ/n46zSqz3aoqU2Kl0jEfD9b62jRKC1qHou/wZYw9UKl3/t/m+i6zjNqL77Kn0qljZ393vdff3u/tCwtOuLsKNHIQb+AeLzWdm7zWz28xs0MwGo3nnVVpfzB5AV37+aBmnEk4geauqWM7z07j7BkKWXVfBNPMITwUejY17lBCNV6MLGHD34UnMH9+WeDows2PM7MaosYVBwsml5L40sxeZ2WVm9oSZDQFfzE/v7lcSsj+/Cawxs+VRXaj5RE+FY/v7t9H4aq0j3HiW80Ts/42E4LKcwv1S+N1VYw/gkth23UsIPOPHUjvhBLMNMzvEzK4ys7Vmtp7wtK1w/5farrLfbZF0HlJwzL+bcMMC8BbCd/+omV0T1XOoxZOx/zcVGd4Bnm/J7cdRZcQhQvHNcr/dPYATC9J9GCGo6ALW+fj6YeX2AcDTPlZhd1OJtOf3b8nv1szazOzLZvZgtB2PRPPEt6Xa47GaY7uW5cH446LS+a6a46jkMVxgHeGCumuRz3YlPAXMe9rdt8aG89s0kfNFPP2V5q/ld1Ns+Y8SzunzGHvKDYC7j0bTljo/F36HsyZSp8fdHyA8RDqdcK79sZkVO39Vc/15Pk3uvjH6dwfCcTMdeDy2H79DeGpdzMXAoWa2K+EJ8ighp28iv/tSqrl251U6Dr5CyFW53MweMrNP15iWqs55AGb2SQstWK6P0rEjpbe/0n6v9fi9hnCj+5ro/6sJuRivjYZrUes1ttr5K12fxqlw31L0e63wmxn3O2b8/UDJe9GCNCV9bVvt7nPcvYPwYG0TocRTfl1V799oeyZ7n50JNTW/7O73EaLC/aJRqwhZh+OY2R6EIhYfAV7g7nMIRSEsv6hiiy8YXkV4KjMn9mp392PLzFPMbrF07UDIbl1dMM1ThKe4e8TG7U54alfNelYDnWYWv9GJz1/J4/F0RvPm0zyTcPH5L0JO2BxCufpy+/JsQrZ2d/SDOTU2Pe7+P+7+N4SnjS8ilEV/ivCjemlsf+/ooVJaqfUUuiNaXj0V7pfC7w7CE73t8gNRxen4DdUqQpny+LE0y937o+mnEZ4E3V4iDRcQciJ3c/cdCVniVmLaQiW/2yJWAdcUpHMHd/8QgLv/yd2PJ1wsf04oXwvVfTe1+GK0zP2j4+c9jN/eYr/VHxSke3t3/zJh++ea2fax6cvtg1qV+27fRSg33UO4KVkUzVPtdxf3AKEl9WofTlQjvh8rne+qOY5eQuljeGyl4cL7R0KxhEJvA35fRdornS+KrrqG+Wv53eQVTr8lWs9qYud2M7No2mrPz3E1/dbc/QJ3PyxavwP/WWSyStefclYRcgbmxfZjh7sXbYrb3dcRGnl4O+H38WN3z29Tpd993LhzLuNvtqq5dueVPQ7cfdjdT3b3vQj1nP7FzF6X35xyO6YWZvZqQpHDtxGKTs4hFFUqdZ2ttN9rPX7zgcyro/+voXIgM9ntr3X+stenuEr3LeW+1zK/mXG/Y8bfDxS9Fy0iNdc2d19PuLd4U2xd5fbvuLRVcZ89ZVRqtWwfMzvZzBZGw7sRsk/zrVJ9D/ikmf2NBS+Mdu72hJ2+NprvJMaCHwhPRRaa2YyCcfHWcm4Ghs3sX81sdvR0dT+rvennYy00IT0D+DyhvPO4nJzo6e9FwH+YWXu0Df9CiNZLpTc+/ypCec0vmdksMzuAUMm/2uYaLwI+ZmYLzWwuoQJX3gxC+cu1wFYzO4aQTZ33JPACM9sxNq6dkHW5wcz2IZSpBkLLYhZyGKYTLkabgdHoKeV3gbPMbKdo2gVm9oYy6yn0a8KJt57+zcy2s9BHxkmEyouF7ic8NX1jtF2nEfZZ3rcJ3+0eAGY238yOj31+MCGLutTTlHZCjttmMzuYcANQrYuAz5jZ3Oh39NEy014GvMjM/t7MpkevV5jZS8xshpm928x2dPcthO93NJqvmu+mFu2ECoXroxv3Uwo+L/yt/hB4k4U+hNqi38DhZrYw2qcrgTOibTiMsRN3PZT7btsJNxtPE266vjjRlbj7c0Av9T++8yqd78oeR2Y2i1A37ooq1/dp4H1m9rHonDfXQtOehxJapiyrivPFZOe/CHi/me1rZtsBn6tise+JTX8m8NPYuf2NZva66PxwMuG4uKGatBYoPPZLMrMXm9mR0U3dZsYa5RiniutPSe7+OCEw+W8z6zCznJntbWbljtMLgPcCb43+z6v0u4+7jXBt7bTQUt4nYp9Vfe2udByY2XHRfYURAosRxp/36tXCXjuhmNBaYJqZ/TuhjnDek8AiM8tF6a6038td04u5BjiCUFTtMUIu2dGEOjt/LjHPZLe/7H1NESWvT0WmLXvfUup7rfCbuRA4LTrHzyPUwcr/RlYAJ0W/8Vx0DO1TJF2pubZZeLD+DkKDIlB5/xamrdJ99pRRKUdmmFDe9SYze4YQwNxFuBDg7j8hlMu+IJr254QWUe4hlEf9I2Hn70+og5B3JeHLe8LM8sUYVgD7Wsgi+3l0cj+OUHH+YcKTm+8RnqzW4gLCRXCAcKF/T4npPkq4sX+I0ELMBYRy2qXSW+idhCe+qwmVAD/n7r1VpvG7hHKXtwO3EspNAuHJBfAxwolxHeEm+hexz+8j/MAfivZdF6HVkncRvpPvMv7mvyMat46QDfo0IZsXQnnQB4AbLWS79hIq35VazzjufivhBHFIldtdjWuiNP0e+C93v7zIetcD/0w4PvoJ32O8X5WvE/bZ5WY2TDiO42l8N+GGuJR/Bs6M5v13xnJCqnEGYT8/TLjwlWy+N/qujyKc3FYTsvj/k7Gg7O+BR6Lv5oNRuqv6bmp0BqG89npC6zU/K/j8S4QLyqCZfTIK5I8n5PytJTxZOoWx88u7CPt7gPBbPH+S6Ysr992eT9j3/YSKvpNtFvw7hO+g7qo431U6jt5EqFOxGsZ1Dlj0CaG7X0coY/5mwpPFRwmVSw9z974qk13yfDHZ+d39N4QizFdG01xZxfJ+QCgx8AShUYuPRcv6C+G8/w3Cfn0T8KYoOK3VuGO/wrQzCY00PBWlaSdCXZFiyl1/Knkv4cbxHsJ5/acULzaY9wtCReMn3D2eg1fpdx/3A8L16hHC8fj8NWYC1+5yx1F3NLyBcD/xLXe/Kvqslu+ikt8RirTdT/gtbGZ80bCfRO9Pm9mt0f/l9nvJa3ox7n4/YRv/EA0PEY6F6710vzjj7pmq2srxqrmviaex0vWpcNqS9y2U/l7L/Wa+QAgc7iC0tnprNA53v5nwoPMswvF7DeNzb/KSvrZ1ReflDYTjrJOx63il/Vt4j1z2PttC57wbKqQnE/Ktf2SSmZ1LqIB4WtJpmSrM7Cjgn939hKTTUo3oKeA1wEHuvjnp9Ei6mdn1wEc86hQzLczsJkIrV3clnZYkmNnVhEq230s6LSIi0jzqlEzqKsox2SbXJK3cfQ2hboFIRe7+qqTTUIy71zMXVEREpCXUVNlfREREREQkqjv0ZzO7LLE0ZLlomYiIiIiI1J+Z/QuwmNC/3XFJpEE5MiIiIiIiUjULLWi+kdCYR2IUyIiIiIiISC2+Ruh/aZtm5ZspM5X9582b54sWLUo6GSIiIiKScbfccstT7j6/8pTJ2XenXXzDc7W3NL9q/bq7Cc2Q5y139+X5ATM7Dljj7reY2eGTTedkZCaQWbRoEStXrkw6GSIiIiKScWZWqhPt1Hjmuef49Kt7ap7vw5f9ZLO7Ly4zyauAvzOzYwn9dnWY2Q/dvVRfjQ2jomUiIiIiIhmUM6v5VYm7f8bdF7r7IkInnlcmEcRAhnJkRERERERkTOWwpLUpkBERERERyRij8UWv3P1q4OoGr6YkBTIiIiIiIhmkHBkREREREWk5VVR5aWkKZEREREREMijrrXopkBERERERyRhDOTIiUkZfH/T2Qn8/LFgAPT3Q3Z10qkRERESyL+s5TiIN09cHK1bA8DAsXBjeV6wI40VERESSlpvAq5W0WnpFUqO3Fzo7Yc4cyOXCe2dnGC8iIiKSPJvAq3UokBGZoP5+6OgYP66jA1avTiY9IiIiInn5fmSUIyMi21iwAIaGxo8bGoKurmTSIyIiIvI8C5X9a321EgUyIhPU0wMDAzA4CKOj4X1gIIwXERERSZpyZESkqO5uWLoU2ttDMbP29jCsVstEREQkDcys5lcrUfPLIpPQ3a3ARURERNKn9aru106BjIiIiIhIBuUyHskokBERERERyaCMxzEKZEREREREssYwchkPZVIbyJjZI8AwMAJsdffFyaZIRERERKR1ZDuMSXEgEznC3Z9KOhEiIiIiIq1GgYyIiIiIiLScrFf2T3M/Mg5cbma3mNmypBMjIiIiItIqbIKvVpLmHJnD3L3fzHYCrjCz+9z92vgEUYCzDGD33XdPIo0iIiIiIuljypFJjLv3R+9rgEuAg4tMs9zdF7v74vnz5zc7iSIiIiIiqZX1HJlUBjJmtr2Ztef/B44C7ko2VSIiIiIikhZpLVq2M3CJmUFI4wXu/ttkkyQiIiIi0jpaLYelVqkMZNz9IeBlSadDRERERKQVGahDTBERERERaT2W7ThGgYyIiIiISBY1Io4xs1nAtcBMQizxU3f/XANWVZECGRERERGRDGpQjsyzwJHuvsHMpgPXmdlv3P3GhqytDAUyIiIiIiIZ06g6Mu7uwIZocHr08rqvqAqpbH5ZREREREQmp1H9yJhZm5ndBqwBrnD3m+qb8uookBERERERyaAJBjLzzGxl7LWscLnuPuLuBwILgYPNbL+Gb0wRKlomIiIiIpJBEyxY9pS7L65mQncfNLOrgKNJoPN65ciIiIiIiGSMATmr/VVxuWbzzWxO9P9s4PXAfY3cllKUIyMiIiIiItXaFTjPzNoImSIXuftlSSREgYyIiIiISAY1ovVld78DOKgBi66ZAhkRERERkQxqTDcy6aFARkREREQkixrUI2ZaKJAREREREcmgbIcxCmRERERERDKnlg4uW5WaXxYRERERkZajHBkRERERkayxzFeRUSAjIiIiIpJFlvHCZQpkREREREQyKNthjOrIiIiIiIhIC1KOjIiIiIhIxhjZz7HI+vaJiIiIiEgGKUdGRERERCSDsl5HRoGMtLy+Pujthf5+WLAAenqguzvpVImIiIgkK+vNL6tombS0vj5YsQKGh2HhwvC+YkUYLyIiIiLZpRwZaWm9vdDZCXPmhOH8e2+vcmVERERkKjNyGS9cphwZaWn9/dDRMX5cRwesXp1MekRERETSwCb4aiWpzpExszZgJdDv7sclnR5JnwULYGhoLCcGwnBXV2JJUp0dERERSYdWi0xqlPYcmY8D9yadCEmvnh4YGIDBQRgdDe8DA2F8ElRnR0RERNIi6zkyqQ1kzGwh8Ebge0mnRdKruxuWLoX29pAD0t4ehpPKAYnX2cnlwntnZxgvIiIi0kw2gb9WkuaiZV8DPgW0J5wOSbnu7vQU3ervDzkxcR0dYbyIiIhIM7VWWFK7VObImNlxwBp3v6XCdMvMbKWZrVy7dm2TUidSWr7OTlzSdXZERERkalLRsmS8Cvg7M3sE+DFwpJn9sHAid1/u7ovdffH8+fObnUaRbaStzo6IiIhIVqUykHH3z7j7QndfBLwDuNLd35NwskQqSludHREREZGsSnMdGZGWlKY6OyIiIjI1haJirVZYrDapD2Tc/Wrg6oSTISIiIiLSUizbcUw6i5aJiIiIiMjkNKKyv5ntZmZXmdk9Zna3mX28IYmvQupzZEREREREpHYNypDZCpzs7reaWTtwi5ld4e73NGZ1pSmQERERERHJoEYULXP3x4HHo/+HzexeYAHQ9EBGRctERERERKRmZrYIOAi4KYn1K0dGRERERCRjJtHB5TwzWxkbXu7uy7dZvtkOwMXAJ9x9qPDzZlAgIyIiIiKSQRNsfvkpd19cdrlm0wlBzI/c/WcTWUk9KJAREREREckaa0wdGTMzYAVwr7t/tf5rqJ7qyIiIiIiIZJDlrOZXFV4F/D1wpJndFr2ObeyWFKccGRERERGRDGpE88vufl2DFl0zBTIiIiIiIlmUinCjcRTIiIiIiIhkjNGYOjJpokBGRERERCSDMh7HKJARybK+Pujthf5+WLAAenqguzvpVImIiEjDWYOaLUsRtVomklF9fbBiBQwPw8KF4X3FijBeREREpNUpkBHJqN5e6OyEOXMglwvvnZ1hvIiIiGRfPlOmllcrUSAjklH9/dDRMX5cRwesXp1MekRERKS5bAJ/rUSBjEhGLVgAQ0Pjxw0NQVdXMukRERGR5lKOjIi0pJ4eGBiAwUEYHQ3vAwNhvIiIiGRbvvllBTIi0nK6u2HpUmhvD8XM2tvDsFotExERmRqyHsio+WWRDOvuVuAiIiIyJU2BHjEVyIiIiIiIZFDG4xgFMiIiIiIiWdRqrZDVSoGMiIiIiEgWZTuOUSAjUqivL3Qa2d8fmjDu6VE9ExEREWktBuQyXrZMrZaJxPT1wYoVMDwMCxeG9xUrwngRERERSQ/lyMiUVSznpbcXOjthzpwwTf69t1e5MiIiItJa0p4hY2bzgX8F9gVm5ce7+5HVzK8cGZmSSuW83HEHdHSMn7ajA1avTiadIiIiIhMygT5kEgh8fgTcC+wJnAE8Avyp2plTG8iY2Swzu9nMbjezu83sjKTTJNkRz3nJ5cJ7ZycMDsLQ0Phph4agqyuBRIqIiIhMQgsEMi9w9xXAFne/xt3/AagqNwbSXbTsWeBId99gZtOB68zsN+5+Y9IJk9bX3x9yYuI6OkIwMzAwNjw0FIaXLGl+GkVEREQmLpnIpEZbovfHzeyNwGqgs9qZUxvIuLsDG6LB6dHLk0uRZMmCBSFIydeBgTC8335jdWX6+0NOzJIlqh8jIiIirSf9cQxfMLMdgZOBbwAdwCeqnTm1gQyAmbUBtwAvBL7p7jclnCTJiJ6eUCcGts156e5W4CIiIiKtzWiJDjHXuft6YD1wBICZvaramVNbRwbA3Ufc/UBgIXCwme0X/9zMlpnZSjNbuXbt2kTSKK2puxuWLoX29pDz0t4ehhXAiIiIiDTNN6ocV1Sqc2Ty3H3QzK4Cjgbuio1fDiwHWLx4sYqdNVAWO4lUzouIiIhkWVqLlpnZocArgflm9i+xjzqAtmqXk9ocGTObb2Zzov9nA68H7ks0UVOUOokUERERaTHpbn55BrADIVOlPfYaAt5a7ULSnCOzK3BeVE8mB1zk7pclnKYpaSp2EpnFHCgRERGRyTKzc4DjgDXuvl+l6Ytx92uAa8zsXHd/1My2c/eNtS4ntTky7n6Hux/k7ge4+37ufmbSaZqq+vunVieRyoESERGRLDCzml9VOJdQ3aMeuszsHqJSV2b2MjP7VrUzpzlHRlKiVFPFWe0kslE5UMrlERERkWYxGlNUzN2vNbNFdVrc14A3AL+Iln27mb2m2plTmyMj6dHTE5omHhyE0dHwPjAQxmdRI3KglMsjIiIizdagHJm6cvdVBaNGqp1XgYxUNNWaKs7nQMVNNgcqnsuTy4X3zs4wXkRERKQhbAIvmJfv3iR6LWtgCleZ2SsBN7PpZvZJ4N5qZ1bRMqnKVGqquFxnmRPV3x9yYuI6OsJ4ERERkbqbeCtkT7n74jqnppQPAl8HFgD9wOXAh6udWYGMZEa96qDkc6Dyy+rqCkHMZAK5qVbPSERERFIgrR3JRNz9KeDdE51fgYxkQr4OSmdnyPkYGgrDEy0CV+8cqEbk8oiIiIiU1piOYczsQuBwQhG0x4DPufuKGpcxC3g7sA74JXAK8BrgQeDzUYBTkQIZyYS093XTiFweERERkXKM+gcy7v7OOizmfGALsD1wMnAX8L/AYYTmnY+rZiEKZKSpGtUEcSvUQZlK9YxEREQkBdLbrNe+7r6fmU0DHnP310bjf2tmt1e7kPRunmROI5sgbkRLYyIiIiItK9+RTK2v5ngOwN23AoUdXFTd/LJyZKRpGln8q9XqoKhzTBEREWm49Fb2X2hm/0MIt/L/Ew0vqHYhCmSkaRpZ/KuV6qDUu2ECERERkaJSG8dwSuz/lQWfFQ6XpEBGmqbRTRC3Sh2UtDdMICIiItlgKc2Rcffz6rEc1ZGRpunpCcW9BgdhdDS8DwyE8VNJf3/IiYrr6IDVhSVERURERKQkBTLSNPniX+3t4Wa+vX1qFqdSwwQiIiLSeBOo6J/SHJxSVLRMmqpVin81Uqs1TCAiIiItyEhtHRkz+093/1czO9HdfzLR5ShHRqTJlDMlIiIijRZaX7aaX01yrIWVfWYyC1GOjEgClDMlIiIiDZfeomK/BdYBO5jZECHu8vy7u3eUmzlPOTIiIiIiIlmU0joy7n6Ku88BfuXuHe7eHn+vdjnKkRERERERyaLUZsgE7n68me0MvCIadZO7r612fuXIiIiIiIhkjU3w1cwkmp0I3AycCLwNuNnM3lrt/MqRERERERHJoLR2iBlzGvAKd18DYGbzgV7gp9XMrEBGRERERCRzWqJfmFw+iIk8TQ0lxhTIiIiIiIhkUforkfzWzH4HXBgNvx34dbUzK5AREREREckgS3ltf3c/xczeDBwWjVru7pdUO78CGRERERGRrDEgl+5ABsDdfwb8bCLzKpARSaG+Pujthf5+WLAAenrUgaaIiIhInAIZqRvdfNdHXx+sWAGdnbBwIQwNheGlS7U/RUREpHot0GrZpKSyCpCZ7WZmV5nZPWZ2t5l9POk0SXn5m+/h4XDzPTwchvv6kk5Z6+ntDUHMnDmQy4X3zs4wXkRERKRqZrW/mp5Em21mL57IvKkMZICtwMnuvi/wt8CHzWzfhNMkZejmu376+6GjY/y4jg5YvTqZ9IiIiEgrmkAQ0+RAxszeBNwG/DYaPtDMflHt/KkMZNz9cXe/Nfp/GLgXWJBsqqQc3XzXz4IFoThZ3NAQdHUlkx4RERGRBjkdOBgYBHD324A9q5059XVkzGwRcBBwU8JJkTLyN99z5oyN0833tqqpR9TTE4rlQQgGh4ZgYACWLGl+ekVERKRFGSnNshhni7uvL6jL49XOnOrNM7MdgIuBT7j7UJHPl5nZSjNbuXbt2uYnUJ7X0xNutgcHYXQ0vA8MhPESVFuPqLs7VOxvbw8BT3u7KvqLiIjIRNgEXk11t5m9C2gzs24z+wZwQ7UzpzZHxsymE4KYH0XtS2/D3ZcDywEWL15cdfQm9Ze/+c7nNnR1hRwE3XyPidcjgrH33t5t91N3t/adiIiITI41qB8ZMzsa+DrQBnzP3b88wUV9FPgs8CxwIfA74PPVzpzKQMZC/tIK4F53/2rS6ZHq6Oa7vDvvDDlVg4Mwdy7ssw/Mnx8CPxEREZFWYGZtwDeB1wOPAX8ys1+4+z21LsvdNxICmc9OJC2pDGSAVwF/D9xpZrdF4051918nlySRievrgwcfhLa2EMRs3gzXXw/77w977ZV06kRERCRzDGhMjszBwAPu/hCAmf0YOB6oOZAxs6soUifG3Y+sZv5UBjLufh0JFNITaZTeXjjggJArs2kTzJ4d3u+4A5YtSzp1IiIikjWhxktDbqcXAKtiw48Bh0xwWZ+M/T8LeAuhG5aqpDKQEcma/n7Ye+/QCtl994XiZTvuGHJnVBxPREREGmJizXrNM7OVseHlUb30unP3WwpGXW9mN1c7vwIZkSbIN0+9887hBSGYaW9PNFkiIiKSWRPu4PIpd19c5vN+YLfY8MJoXM3MrDM2mAP+Btix2vkVyIg0gfqGERERkeZrSNGyPwHdZrYnIYB5B/CuCS7rFkIdGSMUKXsYWFrtzApkRJpAzVOLiIhI000sR6Ysd99qZh8hNJXcBpzj7ndPcFl7TiYtCmQEqK7HeZkcNU8tIiIiTWM0JJABiFoSnnBrwmb25grLL9qHZCEFMvJ8j/OdnaHH+aGhMKze5EVERERaWIMCmTp4U5nPHFAgI9Wppcd5EREREWkRKY1j3P2keixHgYzQ3x9yYuI6OtTjvIiIiEjrMlIbycSY2RuBlxL6kQHA3c+sZl4FMvJ808D5nBgIw11diSVJRERERCYrl+5Axsy+DWwHHAF8D3grUHU/MhPrJkcypacnNAU8OAijo+F9YCCMFxEREZEWlK/sX+uruV7p7u8F1rn7GcChwIuqnVmBjDzfNHB7eyhO1t6uiv4iIiIi0nCboveNZtYFbAF2rXZmFS0TQE0Di4iIiGRPuouWAZeZ2RzgK8CthBbLvlvtzApkRERERESyKKXNL5vZr4ELgLPcfQNwsZldBsxy9/XVLkdFy0REREREMmcC9WOaF/h8B3gj8JCZXWRmSwCvJYgBBTIiIiIiItmU0kDG3S9193cCi4CLgfcCfzWz75vZ66tdjgIZEREREZGssQm+msjdN7r7/7n7EuAo4EDgt9XOrzoyIiIiIiKZlM46MnlmtjPwNuAdhNbKLgLeX+38CmRERERERLIopR1imtkHgHcCLyYULTvF3W+odTkKZEREREREsiilrZYROr78EvB7dx+d6EIUyEhD9fVBb2/oaHPBAujpUX81IiIiIo2XQKWXKrn7P9RjOarsLw3T1wcrVsDwMCxcGN5XrAjjRURERKTBUtpqWb0okJGG6e2Fzk6YMwdyufDe2RnGi4iIiEgDGakNZMxsz3osR4GMNEx/P3R0jB/X0QGrVyeTHhEREZEpJaWBDPDTkDz7/WQWojoy0jALFsDQUMiJyRsagq6uxJIkIiIiMnWkt6RYzsxOBV5kZv9S+KG7f7WqhdQ9WSKRnh4YGIDBQRgdDe8DA2G8iIiIiDRaanvEfAcwQshUaS/yqopyZKRhurth6dKxVsu6umDJErVaJiIiItJ4ltp+ZNz9L8B/mtkd7v6biS5HgYw0VHe3AhcRERGRpstX9k+3G8zsq8BrouFrgDPdfX01M6e2aJmZnWNma8zsrqTTIiIiIiIidXcOMAy8LXoNAd+vdubUBjLAucDRSSdCRERERKQ1pbaOTN7e7v45d38oep0B7FXtzKktWubu15rZoqTTIdvq6xur97JgQai8r+JjIiIiIimT0joyMZvM7DB3vw7AzF4FbKp25tQGMtUws2XAMoDdd9894dRkT7GABWDFitCx5cKFoTnlFStCpX4FMyIiIiJp0txAxsxOBE4HXgIc7O4rK8zyQeB8M9sxGl4HvK/a9bV0IOPuy4HlAIsXL/aEk5MpfX3FA5bZs8O4fN8w+ffeXgUyIiIiIunR1A4u8+4C3gx8p5qJ3f124GVm1hEND9WyspYOZKRxenuLByxXXQXHHz9+2o6OkGvTDCrWJiIiIlKlJscx7n4vgNUYQNUawOSlubK/JKi/PwQocfnhoYJDbWgo9BHTaPlcouHhkEs0PByG+/oav24RERGRlmKA5Wp/tZDUptbMLgT+CLzYzB4zs6VJp2kqWbCgeMByyCEwMACDgzA6Gt4HBsbqzzRSPJcolwvvnZ1hvIiIiIgUmFijZfPMbGXstWzcIs16zeyuIq+CMjuNl9qiZe7+zqTTMJX19ITcDgg5MUNDIWBZGoWT+eJdXV2wZElzinf194ecmLhmFmsTERERaSkTqyPzlLsvLvWhu9ft8bWZvbnI6PXAne6+ptL8qQ1kJFnd3SFoKRWw1Dtw8ZGtkGt7vkylb3yGkXtuxzdtom3vF5PrWvh8LlG+vg40r1ibiIhIVqn+aVYZKS58lbcUOBS4Kho+HLgF2NPMznT3H5SbWYGMlNTd3fgT2eiqR9hy6QWMPnQ/zNqOaa/uwfZ+Cc8tPwt/9lnwUTBj+uFH0fO6d7PiHGPGliFetOF6pj+9ir+yN6/96KHAdo1NqIiISAaVaqVU3SpkRJPjGDNbAnwDmA/8ysxuc/c3lJllGvASd38ymn9n4HzgEOBaQIHMVJfWJy2ja5/k2W99OVS2mTsPtm5ly+8uZXToB9gL5mOdnQD46Chbrr6cvfZ9GR84vpOtXzsDNg4zbcY0Dmq/mtkXXsrowtPJze3cZh2+cSNb77sLgGkvfim2/fZN3UYREZE0K9VKqbpVyIrmRjLufglwSQ2z7JYPYiJronEDZral0swKZDIuzU9atv7xKnzLFnJzogBk+nR85ix83cPkdhmrDGO5HJ7LsXXlDXQ9/TRbX7CJ3N7zn/989Om1PPernzHrPf84bvlbbruFTWd/LRRbw7G2acz64MeZcdArmrF5IiIiqaf6pxk2Vnk/za42s8uAn0TDb43GbQ8MVpo59QXnZHLS1NLX6JrVbPnxcp794ik89+3/ZOSe27DpM8ZNY5YDDLY8N35mg9GtW9l6753YjnPGf9Qxh6233Dh+XUPr2fitr8KMGeTmdpKb+wKYOZNNZ3+N0cF1VaV36yMPs+GC8zh/j27+bfoczpw5l9u+clatmy0iIpJapVopVf3TLLBWaH75w8D3gQOj13nAh939GXc/otLMCmQyrlR/MKtXNzcdo2seZ8s3v8jIHSvx0RFGV/8Vf/h+fPDpcdP5zJmQy0W5KNG40VEYGWH64kOxadNhZGT8wkdGsJmzxo3aesefYesWbNbYeJs5C7ZsCZ9VsOnKy3noox/gjPd8lBv/+iRPbR3liedGWP6pz3FWp87uIiKSDT09yXWrIE0wseaXm8bdHbgOuBL4PXBtNK4qCmQyLi1PWkau+Q2+5Vlsx05s+gxs+3Zsp13x4fX44NP4yFZ88yYYHmT6MSeA5fCBp/Cn1sLgOqYd9jra9juIaa98LT64jvwx7u74hiGmH37U+BVu2QLFfgbueGFuT4HRDcM8c/45/Pz3t/JMtIwonwiAB9ZtoP/3105md4iIiKRCvpXS9vbw8LO9PR3Fz6VOUp4jY2ZvA24mFCl7G3CTmb212vlVRybjSvUHs2RJc9Mx+nAfzN5h3LjcDu2wSxe57n0ZfeQBbPt2ph19AtNeeSRs2sjInbfimzfTtveLsIV7YGbMfPO78LVPsvW+u0POzegI0w4+jOk9x4xbdtu++4EZvnUrNi0c5r51K7S1MW3f/cumdcv994GPsnpzyBWKP5wwQnx0wZITOWXoyWKzi4iItJRmtFIqCZlYPzLN9FngFfk+Y8xsPtAL/LSamRXIZFyl/mCaxebthD90P8yY+fw437qV3OzZzDjpo1hsPADb78C0v33NtsuZPZtZH/8Mo/1/xZ9aS26XLnK7bJu91Lbzrsx889t59mf/93zuDWbMPOFE2nZdUD6tM2dRKW91dHS07OciIiIiiTJaIZDJFXR8+TQ1lBhTIDMFJPGkxUdHGL3jT4z+6VrY8hy2y0LY+hy+eRM2aza+dQsMryd3+DHbBjEVmBltC/eAhXuUnW7WcW9m+n4HsuXWPwHO9Je/grZFe1dc/vQX70Ouo4N509tY/dwIzlhYky+tdvw3/7umNIuIiIjINn5rZr8DLoyG3w78utqZFcg0QVr7cWmkkV/9mJGVf8Bmzg7lLZ98DJszF9zw9QPYtOnkjjyOaa87rqHpaFu0F22L9qppHps2nY6TT+Wdg4N840e/4zmHfP6LATtPz7HP+95b97SKiIiI1I+lPkfG3U8xs7cAr4pGLY/6oqmKApkGS3M/Lo3gqx9l5KpfMnLD76F9R9i+HZs+E585C4bWMe3EpeT2eBHMmhVaIEupaYv2ZK/vX8ip/3wjP3jDW3hiwyZywMGHHcxbrk2g7WoRERGRmqU7kAFw94uBiycyrwKZBpsqPeaOPjPE6JWX4tf9Ft86AiNbYXgQ3zAEu+2FzQj1TnzVQ9hL/ybp5FbFpk1jp0MP42RV6hcREZFWlEtnIGNmwxRvX9YIrTJ3FPlsGwpkGizrPeb6yAgjl/8fftUvYMOGMDLXFjq0nL59aJT+6TWw6+6Aw5wXTGg9U7F4noiIiMjkpDOQcff2eixH/cg0WFr6cWkEd2fkyx+Cy38CmzeGoMVHYXQEcjnYtBFybfjGDfgzwzBzFm371Z4bky+eNzwcgsLh4TDc19eAjRIRERHJCrPaXy1EgUyDZbnH3JGvfBTWD8LzzRsDDuaOzZwBbbmQMwPY3Bcw7b0fx3bYseb1xIvn5XLhvbMzjBcRERGRIiYSxLRYIKOiZQ2Wln5c6m3kgTtgw9DzQYyZ4W052DoSghkc2242vuN8cse+g7ZXvBab4I8j68XzRERERBqixQKTWimQaYJ4Py75uh7nnde8uh75dd55J6xbF3I0Djhg/LrL1UEp9tlef7g0fGgWsklGouJkbcDIaAhwZm1P2zFvJzeJIAbGiuflG0qA7BTPExEREWmYjAcyKlrWREnU9civ86GH4MEHYf16ePjhMJxfd7l0lfrs2Rmzoa0tvHI5yOUwM2xaG8yYBgccQtup36DtkCMmFcRAtovniYiIiDSMipZJvSTRFHN+nbffDjvsALNnw8aNsHo1vOxlY/VMSqWr1GfnDb6PD7R9BmbMCAf9yEiIMgCbvxPT/um0um1DVovniYiIiDRUiwUmtVIg00RJ1PXIr3PdOpg7N4ybPTvkauTX7V46XaU/m48ddBB+z59DjkwUxLD9Dkz71Dfrvh3x4nkiIiIiUkEL5rDUSoFMEyVR1yO/zrlzYfPmEMRs2hTSEF93uXSV+qztzR9m5NVr4LLlofnlfQ6i7YgTG7cxIiIiIlKDbAcyqiPTREnU9civs6sr9Fc5MBCKlnV1ja27XLoqpblt/k60nXQabR/6ooIYERERkTRRHRkpZiI9zSdR1yO+zo0bQxAydy7stdf4NJdLV7nPJrIfRERERKQJcq0VmNRKgcwE5Fvy6uwM9UeGhsLw0qXVBTPNvtGvtM5KwUip+SezH0RERESk0RTISIEkWh9rlHLBCJQPcLK0H0RERESypfWKitUqtYGMmR0NfJ3QxeL33P3LCSfpea3W03y5HJdSwcgFF4RGAcrltqRxP6iom4iIiMjUkMrK/mbWBnwTOAbYF3inme2bbKrG5FsCi0trT/OVOuHs7w/BR1xHB9x001iAk8uF987Osb5lIH37IYkOR0VERERSych8Zf9UBjLAwcAD7v6Quz8H/Bg4PuE0PS9NPc339cHZZ8Npp4X3wpv2eI5LsYCkVDDiXjzAWb16bDhN+wEqb6uIiIjIlNLkQMbMvmJm95nZHWZ2iZnNqc+GFJfWQGYBsCo2/Fg0LhXyLYG1t4ccjfb2ZCq4V5MDUSrHJR+QlApGDjmkcm5LWvZDXqVtFREREZlSLFf7a3KuAPZz9wOA+4HPTHobykhtHZlqmNkyYBnA7rvv3tR1p6Gn+UqV7fv64MEH4YYbYNddYZ99YOedxwckpZqEhhAUQQgGhoZCgJP/LC8N+yEviQ5HRURERCRw98tjgzcCb23k+tIayPQDu8WGF0bjxnH35cBygMWLF3tzkpYe5Srb53NrFi6EdevC67rr4IADoK1tfEBSKhhpdp83k9XTU13wJSIiIpJ9FsraJ+cfgP9r5ArSGsj8Ceg2sz0JAcw7gHclm6T0KZcDEc+t6eiA++6DJ56Axx6DM8+sLiBJU25LNZLocFREREQklYyJdiMzz8xWxoaXR5kHYbFmvcAuReb7rLtfGk3zWWAr8KMJpaBKqQxk3H2rmX0E+B2h+eVz3P3uhJOVOuVyIM47byy3Zuedw2t0NNzgZ/nGvtWCLxEREZHGmVCOzFPuvrjUh+5etlknM3s/cBzwOndvaImpVAYyAO7+a+DXSacjzcrlQKi+iIiIiMgU1+TWlKN+ID8FvNbdNzZ6fakNZKQ6pXIg0l5fRB1XioiIiDSS1aMVslr9LzATuMJCU843uvsHG7WytDa/LJOUtqaR49RxpYiIiEgTNLkfGXd/obvv5u4HRq+GBTGgHJlMS2t9kUrNRouIiIhIHUwyMEk7BTJ1lpUiU43cjnLNRouIiIhIvWQ7kFHRsjrKSpGpRm9HviGCODVEICIiIlJHRtOLljWbApk6iheZyuXCe2dnGN9KGr0dPT2h4YHBwdAk9OBgGO4p25ifiIiIiFTPMKv91UoUyNRRf38oIhXX0QGrVyeTnolq9HakuSECERERkczIeI6M6sjUUVb6bqlmOyZbhyatDRGIiIiIZEe28yyyvXVNlpUiU5W2Iyt1gURERESkdSlHpo7yRabyORVdXaEDysnkPCTRClql7Uhr88lZaTFOREREpC4ynmWhQKbO6llkKp/z0dkZcj6GhsLwZOqTVHuzX2470th8ciP2lYiIiEjrMrIeyWR761pcvVsPq1eRsDQ2n5yVFuNERERE6sIm+GohypFJsWI5H5s2wdVXT6z4VL2KhPX0hAAIQk7M0FCoQ7NkSfXLqLc05hKJiIiIJMqynWeR7a1rcYU5H08+CVddBTNnTixHpV7NKqex+eQ05hKJiIiISOMoRybFCnM+br01/P/yl48Vn4Lqc1Tq2Tx02ppPTmMukYiIiEiilCMjSSnM+di8GY44AnbeeWyaWnJUstI8dDFpzCUSERERSY5BbgKvFqIcmZSL53ycfXYoThZXS45KI5qHTpO05RKJiIiIJKu1ApNaKZBpIfUoPqWbfREREZEpwADLdiCjomUtRMWnRERERKRqZrW/WohyZFqMclREREREpDqtFZjUSoFMhvT1jdV/qbWPGREREZFWoPudarVeDkutVLQsI/r6Qv2Z4eGJ9TEjIiIikna636mR5Wp/tZDWSq2U1NsLnZ2hj5h8HzOdnWG8iIiISBbofqc2Ga8io0AmK/r7Q0tmcbX0MSMiIiKSdrrfqYGhHBlpDQsWhOaY42rpY0ZEREQk7XS/UyObwKuFKJDJiJ6e0KfM4CCMjob3gYEwfirq6wsdiJ52WnhX2VkREZHWp/sdiVMgkxHqY2aMKgKKiIhkk+53amGEW/1aX60jdc0vm9mJwOnAS4CD3X1lsilqHepjJohXBISx995e7R8REZFWp/udGjS5qJiZfR44HhgF1gDvd/eG1WBKXSAD3AW8GfhO0gmR1tTfH3Ji4jo6wngRERFpHPXxkjLNr7z/FXf/NwAz+xjw78AHG7Wy1OUfufu97v6XpNMhrUsVAUVERJpPRbtTqMntL7t7/A5se8AntcAKUhfIiEyWKgKKiIg0n/p4SZmJBDF16EjGzP7DzFYB7ybkyDRMIoGMmfWa2V1FXsfXuJxlZrbSzFauXbu2UcmVFqOKgCIiIs2nPl7SaELtL8/L319Hr2XjlljhPt7dP+vuuwE/Aj7SyK1LpI6Mu9fl2bi7LweWAyxevLihWVfSWlQRUEREpLnyRbvzjeyAina3qKfcfXGpD2u4j/8R8Gvgc3VJVREqWiYiIiIik6ai3SnU5KJlZhZ/jHw8cN+kFlhB6gIZM1tiZo8BhwK/MrPfJZ0mERERESlPRbtTqPl1ZL4cFTO7AzgK+PjkN6K01DW/7O6XAJcknY4sSluTiGlLj4iIiEyOinZPbe7+lmauL3WBjDTGFVfAf/0XbN0KO+0Ezz4bmkRs9JOSUsFKvonGzs7QROPQUHPSIyIiIjJl1KEVsjRLXdEyqb++vhDEtLXBLruEIObOO2FkpLFNIpZrT15NNIqIiGRDXx+cfTacdlp4V78x0iwKZKaA3l7YsgXmzg2B+XbbhddjjzW2ScRywYqaaBQREWl96gQzzSbS9HJr5eCoaFkGFRbnuvPOkBOzeTPMnh2mmT0bnngCXvOaxqWjvz+c1OI6OsbSpSYaRUREWlv8oSWMvff2qqh4KqhombSSYk9GHnww5MA88wxs2gTusG4dTJvW2CYR88FKXD5YURONIiIire+OO+DPf4ZLL4VrroEnn1QJi9TIfoaMApmsKVac64ADYNUq2H9/mDkz5MSMjMAnP9nYpyXlghU10SgiItLa+vrgoYdg/fpwv7FpE1x/fXiAqhIWaZHtSEZFyzKmWHGuvfcOuTF77QWzZsGrX92cpo7zwUq+mFtXFyxZMrZeNdEoIiLSunp7w8PSO+8cK76+eXPIpVm2LOnUyVSgQCZjStU92X9/+NCHmp8eBSsiIiLZ1N8fHpZ2dMB994WSF3PmhMaFdO1Pi9bKYamVApmM6ekJdWQgnFiGhkJxriVLql+GOqoUERGRSvIPT3feObwgBDPt7YkmS6YQ1ZFpsGa3rT7ZuidqRlFERESqoYZ7WoBZ7a8WohyZBkqq9/rJFOdSM4oiIiJSjUp1YSVprReY1EqBTAO1YlBQru8XERERkTjVhZUkqWhZA7Vi7/Xl+n4RERERkdZhE/hrJcqRaaBW7L2+Ho0FiIiISOOpcR6pKONFy5Qj00CtWAlOHVWKiIiknxrnkeqoQ0yZoFatBKfyriIiIunWivVwpclaLy6pmQKZBlNQICIiIvWmxnmkKipaJiIiIiJposZ5pLKJFCtrrcBHgYyIiIhIi2nFergi9aZARkRERKTFqHEeqUq2M2RUR0ZERESkFakerlTWYpFJjRTIiIiIiIhkkSr7i4iIiIiIpItyZEREREREMkk5MiIiIiIiIqmiQEZEREREJHMs1JGp9VWPNZudbGZuZvPqssASVLRMRERERCRrjEQq+5vZbsBRwF8bvS7lyIiIiIiISL2cBXwK8EavSDkyIiIiIiKZ1NwcGTM7Huh399utCblB5t7wYKkpzGwt8GjS6Wgx84Cnkk6EpIaOBymkY0LidDxI3FQ/HvZw9/lJJ6IcM/st4Xuq1Sxgc2x4ubsvjy23F9ilyHyfBU4FjnL39Wb2CLDY3Rt2nGQmkJHamdlKd1+cdDokHXQ8SCEdExKn40HidDxIITPbH/g9sDEatRBYDRzs7k80Yp0qWiYiIiIiIpPi7ncCO+WHm5Ejo8r+IiIiIiLScpQjM7UtrzyJTCE6HqSQjgmJ0/EgcToepCx3X9TodaiOjIiIiIiItBwVLRMRERERkZajQGYKMrOjzewvZvaAmX066fRI8szsETO708xuM7OVSadHmsvMzjGzNWZ2V2xcp5ldYWZ90fvcJNMozVXimDjdzPqj88RtZnZskmmU5jGz3czsKjO7x8zuNrOPR+N1npBEKZCZYsysDfgmcAywL/BOM9s32VRJShzh7geqOc0p6Vzg6IJxnwZ+7+7dhOY09dBjajmXbY8JgLOi88SB7v7rJqdJkrMVONnd9wX+FvhwdO+g84QkSoHM1HMw8IC7P+TuzwE/Bo5POE0ikiB3vxYYKBh9PHBe9P95wAnNTJMkq8QxIVOUuz/u7rdG/w8D9wIL0HlCEqZAZupZAKyKDT8WjZOpzYHLzewWM1uWdGIkFXZ298ej/58Adk4yMZIaHzGzO6KiZypGNAWZ2SLgIOAmdJ6QhCmQERGAw9z95YQihx82s9cknSBJDw/NW6qJSzkb2Bs4EHgc+O9EUyNNZ2Y7ABcDn3D3ofhnOk9IEhTITD39wG6x4YXROJnC3L0/el8DXEIogihT25NmtitA9L4m4fRIwtz9SXcfcfdR4LvoPDGlmNl0QhDzI3f/WTRa5wlJlAKZqedPQLeZ7WlmM4B3AL9IOE2SIDPb3sza8/8DRwF3lZ9LpoBfAO+L/n8fcGmCaZEUyN+wRpag88SUYWYGrADudfevxj7SeUISpQ4xp6CoycyvAW3AOe7+H8mmSJJkZnsRcmEApgEX6JiYWszsQuBwYB7wJPA54OfARcDuwKPA29xdlb+niBLHxOGEYmUOPAL8U6x+hGSYmR0G/AG4ExiNRp9KqCej84QkRoGMiIiIiIi0HBUtExERERGRlqNARkREREREWo4CGRERERERaTkKZEREREREpOUokBERERERkZajQEZERERERFqOAhkREREREWk5CmRERERERKTl/H+kCGUnmebFNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(15,5))\n",
    "\n",
    "# define the estimators\n",
    "estimated_values = []\n",
    "for beta in np.logspace(-4,4):\n",
    "    estimated_values += [M_estimator(X, beta, name, p=4)]\n",
    "\n",
    "plt.scatter(X[:,0],X[:,1],alpha = 0.4, c='b')    \n",
    "im = plt.scatter(np.array(estimated_values)[:,0], \n",
    "            np.array(estimated_values)[:,1], \n",
    "            cmap = plt.get_cmap('Reds'), \n",
    "            c = np.linspace(-4,4), alpha = 0.7)\n",
    "cbar = fig.colorbar(im, ax=ax)\n",
    "cbar.set_label('log of Value of Beta')\n",
    "\n",
    "\n",
    "ax.set_title('Scatter plot of dataset (in blue) and estimated mean (in red).'\n",
    "             ' One red point is one value estimated with the associated Beta.')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Heavy tail illustration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "def geometric_median(X, eps=1e-5):\n",
    "    y = np.mean(X, 0)\n",
    "\n",
    "    while True:\n",
    "        D = cdist(X, [y])\n",
    "        nonzeros = (D != 0)[:, 0]\n",
    "\n",
    "        Dinv = 1 / D[nonzeros]\n",
    "        Dinvs = np.sum(Dinv)\n",
    "        W = Dinv / Dinvs\n",
    "        T = np.sum(W * X[nonzeros], 0)\n",
    "\n",
    "        num_zeros = len(X) - np.sum(nonzeros)\n",
    "        if num_zeros == 0:\n",
    "            y1 = T\n",
    "        elif num_zeros == len(X):\n",
    "            return y\n",
    "        else:\n",
    "            R = (T - y) * Dinvs\n",
    "            r = np.linalg.norm(R)\n",
    "            rinv = 0 if r == 0 else num_zeros/r\n",
    "            y1 = max(0, 1-rinv)*T + min(1, rinv)*y\n",
    "\n",
    "        if euclidean(y, y1) < eps:\n",
    "            return y1\n",
    "\n",
    "        y = y1\n",
    "        \n",
    "def pareto(df, n, d):\n",
    "    return (np.random.pareto(df, size=[n,d])+1)\n",
    "\n",
    "def gmom(X, K):\n",
    "    perm = np.random.permutation(len(X))\n",
    "    blocks = np.array_split(perm, K)\n",
    "    return geometric_median(np.array([np.mean(X[b], axis=0) for b in blocks]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19"
     ]
    }
   ],
   "source": [
    "res = []\n",
    "alpha = 2.1\n",
    "power = 2\n",
    "n = 1000\n",
    "d = 100\n",
    "mean = alpha/(alpha -1)*np.ones(d)\n",
    "M = 20\n",
    "\n",
    "for iteration in range(M):\n",
    "    print(\"\\r\"+str(iteration), end=\"\" )\n",
    "    X = pareto(alpha, n, d)\n",
    "    X = np.vstack([X, 300*np.ones([2,d])])\n",
    "    res += [ [np.linalg.norm(M_estimator(X)-mean),\n",
    "              np.linalg.norm(M_estimator(X, name= \"Catoni\")-mean),\n",
    "              np.linalg.norm(M_estimator(X, name= \"Polynomial\", p=5)-mean),\n",
    "            np.linalg.norm(gmom(X, 9) - mean), \n",
    "            np.linalg.norm(geometric_median(X)-mean)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Error')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAW40lEQVR4nO3df5DcdX3H8efLS8AWjdyZa3UIENpiPeZELSu25apcq53otESrFm5aKXI2tZUYtTSFnpVIG1TUcaYRayNJiZ12gTpqo6WDikfxKJRcBEOSK5iCyqGWk1xFUyMhvPvHfo8sl729vdx997t7n9djZuf2+/1+dve93937vvb7+f5SRGBmZul6RtEFmJlZsRwEZmaJcxCYmSXOQWBmljgHgZlZ4pYUXcBcLV++PFauXFl0GWZmbWXnzp3fj4juWtPaLghWrlzJ6Oho0WWYmbUVSd+aaZq7hszMEucgMDNLnIPAzCxxuQWBpK2SHpG0e4bpvydpl6R7Jf2HpBfnVYuZmc0szzWC64BVdaY/CLwyIl4E/BWwOcdazMxsBrkFQUTcBuyvM/0/ImIyG7wTWJFXLWbWPOVymd7eXjo6Oujt7aVcLhddks2iVXYfHQT+baaJktYAawBOOeWUZtVkZnNULpcZGhpiy5Yt9PX1MTIywuDgIAADAwMFV2czUZ6noZa0EvhCRPTWadMPfBzoi4hHZ3vOUqkUPo7ArDX19vayadMm+vv7nxo3PDzM2rVr2b275uZCaxJJOyOiVGtaoWsEks4ErgVe00gImFlrGxsbo6+v72nj+vr6GBsbK6gia0Rhu49KOgX4DPDmiLi/qDrMbOH09PQwMjLytHEjIyP09PQUVJE1Is/dR8vAHcAvShqXNCjpbZLeljV5L/Bc4OOS7pHk/h6zNjc0NMTg4CDDw8McOnSI4eFhBgcHGRoaKro0qyO3rqGIqLtlKCLeCrw1r9c3s+ab2iC8du1axsbG6OnpYePGjd5Q3OJy3VicB28sNjObu3obi32KCTOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwSl+fF67dKekTS7hmmv1DSHZJ+IunSvOowM7P68lwjuA5YVWf6fuAdwIdzrMHMzGaRWxBExG1UFvYzTX8kInYAh/KqwczMZtcW2wgkrZE0Kml0YmKi6HLMzBaVtgiCiNgcEaWIKHV3dxddjpnZotIWQWBmZvlxEJiZJW5JXk8sqQycCyyXNA5cASwFiIhPSHoeMAosA56U9E7gjIh4LK+azMzsaLkFQUQMzDL9e8CKvF7fzMwa464hM7PEOQjMzBLnIDAzS5yDwMwscbltLDazxUvSgjxPRCzI89j8OAjMbM4aWYBL8oK+TbhryMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PEOQjMzBLnIDAzS5yDwMwsJ+Vymd7eXjo6Oujt7aVcLhddUk25BYGkrZIekbR7humS9DeS9knaJemX8qrFzKzZyuUyQ0NDbNq0iYMHD7Jp0yaGhoZaMgzyXCO4DlhVZ/prgNOz2xrgb3OsxcysqTZu3MiWLVvo7+9n6dKl9Pf3s2XLFjZu3Fh0aUfJLQgi4jZgf50mq4FPRcWdwImSnp9XPWZmzTQ2NkZfX9/TxvX19TE2NlZQRTMrchvBScBDVcPj2bijSFojaVTS6MTERFOKMzObj56eHkZGRp42bmRkhJ6enoIqmllbbCyOiM0RUYqIUnd3d9HlmJnNamhoiMHBQYaHhzl06BDDw8MMDg4yNDRUdGlHKfJSlQ8DJ1cNr8jGmZm1vYGBAQDWrl3L2NgYPT09bNy48anxraTIINgOXCLpeuDlwA8i4rsF1mNmtqAGBgZacsE/XW5BIKkMnAsslzQOXAEsBYiITwA3Aa8F9gH/B7wlr1rMzGxmuQVBRNSNwYgI4O15vb6ZmTWmLTYWm1lzdXV1IWleN2Bej+/q6ip4LqSjyG0EZtaiJicnqay0F2cqTCx/XiMwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0ucg8DMLHHefdTMjhJXLIMNzym+BmsKB4GZHUXve6wljiOIDYWW0JCFOt6hyPntIDAzm4dGFuCSCg/WeryNwMwscQ4CM7PEOQjMzBLnbQRW12LYEGZm9TkIrK7FsCHMzOpz15CZWeIcBGZmics1CCStknSfpH2SLqsx/VRJt0jaJelWSSvyrMfMzI6WWxBI6gCuAV4DnAEMSDpjWrMPA5+KiDOBK4H351WPmZnVlucawdnAvoh4ICIeB64HVk9rcwbwlez+cI3pZmaFSuH6zXkGwUnAQ1XD49m4al8Hfie7/3rg2ZKem2NNZmZzMnX95iJvk5OTub7HojcWXwq8UtLdwCuBh4HD0xtJWiNpVNLoxMREs2s0M1vU8gyCh4GTq4ZXZOOeEhHfiYjfiYiXAkPZuP+d/kQRsTkiShFR6u7uzrFkM7P05BkEO4DTJZ0m6TjgAmB7dQNJyyVN1XA5sDXHeszMrIbcgiAingAuAW4GxoAbI2KPpCslnZc1Oxe4T9L9wM8CG/Oqx8zMalO7nRqgVCrF6Oho0WVYFZ9iYvFphc+0FWoACr9S21M2/GBeD5e0MyJKtab5XENmZnWkcLW2ovcaMjOzgs0aBJKeIelXm1GMmZk136xdQxHxpKRrgJc2oR5rsq6urgU5WGU+1y3o7Oxk//79867BFtZCXYviWHV2dhb6+ilpdBvBLZLeAHwmiu4sswU1ddRkkYpe4NjRFuI70TIbexdA0d/RvEOx0SD4I+DdwGFJPwYEREQsy60yM7MWkEIoNhQEEfHsvAsxM7NiNLz7aHYQ2CuywVsj4gv5lGRmZs3U0O6jkj4ArAP2Zrd1knztADOzRaDRNYLXAi+JiCcBJG0D7qZyfiAzM2tjczmg7MSq+y1yzLWZmc1Xo2sEVwF3SxqmssfQK4CjrkFsZmbtZ9YgyE4T/STwy8DLstF/HhHfy7MwMzNrjkaPLF4fETcy7XoCZmbW/hrtGvqypEuBG4ADUyMjwucFMEtQo0faztaulQ+ySkmjQXB+9vftVeMC+LmFLcfM2oEX4ItLo9sILouIG5pQj5lZW1kMa0ez7j6aHTvwZ02oxcys7UTEgtyK1OhxBF+WdKmkkyV1Td1yrczMzJqi0SA4n8r2gduAndlt1gsHS1ol6T5J+yQdddyBpFMkDUu6W9IuSa+dS/FmZjZ/jZ599LS5PrGkDuAa4NXAOLBD0vaI2FvV7D3AjRHxt5LOAG4CVs71tczM7NjVXSOQtL7q/pumTbtqluc+G9gXEQ9ExOPA9cDqaW0CmLqmwXOA7zRStJmZLZzZ1gguAK7O7l8O/HPVtFXAX9R57EnAQ1XD48DLp7XZAHxR0lrgBOBVtZ5I0hpgDcApp5wyS8k2F3HFMthQ7Kmj4gpf38isSLMFgWa4X2v4WAwA10XERyT9CvAPknqnznI6JSI2A5sBSqWSd2BeQHrfY4XvsSCJ2FBoCWZJm21jccxwv9bwdA8DJ1cNr8jGVRsEbgSIiDuAZwLLZ3leMzNbQLMFwYslPSbph8CZ2f2p4RfN8tgdwOmSTpN0HJVupunnKvo28BsAknqoBMHEnN+FmZkds7pdQxHRcaxPHBFPSLoEuBnoALZGxB5JVwKjEbEd+FPgk5LeRWUN46Ioup/CzCwxDV+z+FhExE1UdgmtHvfeqvt7gXPyrMHMzOqbyxXKzMxsEXIQmJklzkFgZpa4XLcRWHto9DS6eens7Cz09c1S5yBI3ELspCWp8IPSzOzYuWvIzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8T5FBM1LNS5dxbDaRcanReztVsM88JssXIQ1NDIQiuV8+uk8B7NUueuITOzxDkIzMwSl2sQSFol6T5J+yRdVmP6RyXdk93ul/S/edZjZmZHy20bgaQO4Brg1cA4sEPS9uyC9QBExLuq2q8FXppXPWZmVlueawRnA/si4oGIeBy4Hlhdp/0AUM6xHjMzqyHPIDgJeKhqeDwbdxRJpwKnAV+ZYfoaSaOSRicmJha8UDOzlLXKxuILgE9HxOFaEyNic0SUIqLU3d3d5NLMbC7K5TK9vb10dHTQ29tLuewV/VaXZxA8DJxcNbwiG1fLBbhbyKztlctl1q1bx4EDBwA4cOAA69atcxi0uDyDYAdwuqTTJB1HZWG/fXojSS8EOoE7cqzFzJpg/fr1LFmyhK1bt3Lw4EG2bt3KkiVLWL9+fdGlWR25BUFEPAFcAtwMjAE3RsQeSVdKOq+q6QXA9eFDWM3a3vj4ONu2baO/v5+lS5fS39/Ptm3bGB8fL7o0qyPXU0xExE3ATdPGvXfa8IY8azAzs/paZWOxmS0CK1as4MILL2R4eJhDhw4xPDzMhRdeyIoVK4ouzepIMgi6urqQNK8bMK/Hd3V1FTwXzBbe1VdfzeHDh7n44os5/vjjufjiizl8+DBXX3110aVZHUmefXRycrLws2ou1KmuzVrJwMAAABs3bkQSJ5xwAlddddVT4601qegF4lyVSqUYHR2d13O0wimkW6EGM0uHpJ0RUao1LcmuITMzO8JBYGaWOAeBmVniHARmZolzEJiZJS7J3UfjimWw4TnF12Bm1gKSDAK977HCd92UhE+uYWatwF1DZmaJcxCYmSXOQWBmljgHgZlZ4hwEZmaJS3KvISj+7J+dnZ2Fvr6Z2ZQkg2Ahdh312UPNbLFw15CZWeJyDQJJqyTdJ2mfpMtmaPO7kvZK2iPpn/Ksx2w+5ntVu+qr25m1kty6hiR1ANcArwbGgR2StkfE3qo2pwOXA+dExKSkn8mrHrP5aqQr0F2G1o7yXCM4G9gXEQ9ExOPA9cDqaW3+ELgmIiYBIuKRHOsxM7Ma8gyCk4CHqobHs3HVXgC8QNLtku6UtKrWE0laI2lU0ujExERO5ZqZpanojcVLgNOBc4EB4JOSTpzeKCI2R0QpIkrd3d3NrdDMbJHLMwgeBk6uGl6Rjas2DmyPiEMR8SBwP5VgMDOzJskzCHYAp0s6TdJxwAXA9mltPkdlbQBJy6l0FT2QY01mZjZNbkEQEU8AlwA3A2PAjRGxR9KVks7Lmt0MPCppLzAM/FlEPJpXTWZmdjS1265upVIpRkdHiy7DuwkuQl1dXUxOThZaQ2dnJ/v37y+0BlucJO2MiFKtaUmeYsKslsnJycLD3QecWRGK3mvIzMwK5iAwM0ucg8DMLHEOAjOzxDkIzMwS572Gamh0z43Z2hW9B4rNTVyxDDY8p/gazJrMQVCDF+Bp0vseK/yzl0RsKLQES5C7hszMEucgMDNLnIPAzCxxDgIzs8R5Y7FZlaLP9dPZ2Vno61uaHARmmYXYY8hnpbV25K4hM7PEOQjMzBLnIDAzS5yDwMwscQ4CM7PE5RoEklZJuk/SPkmX1Zh+kaQJSfdkt7fmWY+ZmR0tt91HJXUA1wCvBsaBHZK2R8TeaU1viIhL8qrDzMzqy3ON4GxgX0Q8EBGPA9cDq3N8PTMzOwZ5BsFJwENVw+PZuOneIGmXpE9LOrnWE0laI2lU0ujExEQetZqZJavojcWfB1ZGxJnAl4BttRpFxOaIKEVEqbu7u6kFmpktdnkGwcNA9S/8Fdm4p0TEoxHxk2zwWuCsHOsxM7Ma8gyCHcDpkk6TdBxwAbC9uoGk51cNngeM5ViPmZnVkNteQxHxhKRLgJuBDmBrROyRdCUwGhHbgXdIOg94AtgPXJRXPWZmVpva7UyJpVIpRkdHiy7DErRQp6hut/85Wxwk7YyIUq1pPg21WYO8ALfFqui9hszMrGAOAjOzxDkIzMwS5yAwM0ucg8DMLHEOAjOzxDkIzMwS5yAwM0tc2x1ZLGkC+FbRdQDLge8XXUSL8Lw4wvPiCM+LI1phXpwaETVP39x2QdAqJI3OdLh2ajwvjvC8OMLz4ohWnxfuGjIzS5yDwMwscQ6CY7e56AJaiOfFEZ4XR3heHNHS88LbCMzMEuc1AjOzxDkIzMwSl3QQSPrRtOGLJH1slsecK+kL+VbWfJKeJ+l6Sf8taaekmyS9YIa2J0r6k3m+3rWSzpjPczSTpMOS7pG0W9I/S/rpOm1n/R6ZzaSI70/SQVAESR1F1zCdKtdg/Cxwa0T8fEScBVwO/OwMDzkRmFcQRMRbI2LvfJ6jyX4cES+JiF7gceBtRRdktlAcBDOQdJ2kN1YNV689LJP0r5Luk/QJSc/I2vympDskfS371fisbPw3JX1Q0teANzX3nTSkHzgUEZ+YGhERXwfulnRL9n7ulbQ6m/wB4OezX8gfUsWHsl/L90o6H55ae7pV0qcl/Zekf8xCh2x8yx5gM4uvAr8gqUvS5yTtknSnpDOrG0l6tqQHJS3NhpdVD7cDSX+Zfc9HJJUlXZp9dh+VNCppTNLLJH1G0jck/XXVY9+dfSd2S3pnNm5l9l24TtL92XfiVZJuzx5/dmFvtkHznCe/L+mu7H/n76Z+GEp6SzY/7gLOafqbiohkb8Bh4J6q27eBj2XTrgPeWNX2R9nfc4GDwM8BHcCXgDdSOYT8NuCErN2fA+/N7n8TWF/0+60zH94BfLTG+CXAsuz+cmAfIGAlsLuq3Ruy+dBBZS3i28Dzs3n1A2AFlR8ddwB92WNuBUpFv/c5zKOpz38J8C/AHwObgCuy8b8O3JPdv6jqe/T3wOuy+2uAjxT9Xubwnl+W/V88E3g28A3g0uyz+2DWZh3wnezzPh4YB54LnAXcC5wAPAvYA7w0++48Abwo+07sBLZm36vVwOeKft85zpMe4PPA0qzdx4ELs3bfBrqB44Dbp74/zbqlfvH6H0fES6YGJF0ENPIr9a6IeCB7TBnooxIOZwC3Zz96j6Oy4Jtyw8KU3FQCrpL0CuBJ4CRqdxf1AeWIOAz8j6R/p/IP8xiVeTUOIOkeKguCkfxLX3A/ldUPlTWCLcB/UglBIuIrkp4radm0x10LrAc+B7wF+MOmVLswzgH+JSIOAgclfb5q2vbs773Anoj4LoCkB4CTqXwnPhsRB7LxnwF+LXvcgxFxbzZ+D3BLRISke6l8P1rZfOfJWcCObBnxU8AjwMupdMtOZO1vAGpun8tL6kFQzxNkXWdZ189xVdOmH3wRVBaaX4qIgRme78CCV7hw9lBZq5nu96j8SjkrIg5J+iaVX0Jz8ZOq+4dp3+/c0340AGT/zHVFxO1Zd8i5QEdE7M6luuab+lyf5Omf8ZPM/hlPb1/9XO36/YDZ54mAbRFxefWDJL2uKdXV4W0EM/smlfQGOA+o7tc9W9JpWUCcT+UX7p3AOZJ+AUDSCZphr5sW9BXgeElrpkZk/d2nAo9kIdCfDQP8kMpq8ZSvAudL6pDUDbwCuKs5pRfqq1TCkmxB//2IeKxGu08B/0Slm6id3A78tqRnZtu7fmsOj/0q8DpJPy3pBOD12bh2N595cgvwRkk/A5BtYzqVyprlK7M1yqUUsB3RQTCzT1L5cL4O/ApP/0W/A/gYMAY8SGUVeIJK33BZ0i4q3UIvbGrFxygqHZavB16lyu6je4D3AzcBpWyV/ULgv7L2j1LpAtst6UNU9jjaBXydSqisj4jvFfBWmm0DcFb2eX8A+IMZ2v0j0AmUm1TXgoiIHVS6O3YB/0aly+MHDT72a1S2s91FZUF3bUTcnU+lzTPPebIXeA/wxew78yXg+VkX0gYqy4zbqSxXmsqnmDDLmSp7n62OiDcXXctcSXpWRPxIleMmbgPWZAv5ZC3GedLO/XFmLU/SJuA1wGuLruUYbVblwL9nUunfbusF3gJZdPPEawRmZonzNgIzs8Q5CMzMEucgMDNLnIPAzCxxDgIzs8T9P1i4R5U9TjleAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.boxplot(np.array(res), labels= ['Huber', 'Catoni', 'Poly', \"gmom\", \"gmed\"])\n",
    "plt.ylabel('Error')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}