Synthetic_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import randint,rv_discrete\n",
    "import matplotlib.pyplot as plt\n",
    "from estimator import *\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def true_delta(supports, myrv1_pmf, myrv2_pmf, epsilon):\n",
    "    \"\"\" computes delta given the P and Q\n",
    "    \"\"\"\n",
    "    delta = 0\n",
    "    for i in range(len(supports)):\n",
    "        p1 = myrv1_pmf[i]\n",
    "        p2 = myrv2_pmf[i]\n",
    "        delta += max(p2-np.exp(epsilon)*p1,0)\n",
    "    return delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "S = 100\n",
    "\n",
    "supports = np.arange(1,S+1,1)\n",
    "\n",
    "\n",
    "pmfs_uniform = S*[1/S]\n",
    "myuniform = rv_discrete(name='MyUniform', values=(supports, pmfs_uniform))\n",
    "\n",
    "\n",
    "alpha = -0.6\n",
    "pmfs_zipf = 1/supports**alpha /np.sum(1/supports**alpha)\n",
    "myzipf = rv_discrete(name='MyZipf', values=(supports, pmfs_zipf))\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "epsilon = 0.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_delta12 = true_delta(supports, pmfs_uniform, pmfs_zipf, epsilon = epsilon)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_once(n_mean):\n",
    "    np.random.seed() \n",
    "\n",
    "    n1 = np.random.poisson(n_mean)\n",
    "    n2 = np.random.poisson(n_mean)\n",
    "    n12 = np.random.poisson(n_mean)\n",
    "    n22 = np.random.poisson(n_mean)\n",
    "\n",
    "\n",
    "    rvs1 = myuniform.rvs(size = n1)\n",
    "    count_1 = {}\n",
    "    for support in supports:\n",
    "        count_1[support] = np.count_nonzero(rvs1 == support)\n",
    "\n",
    "\n",
    "    rvs2 = myzipf.rvs(size = n2)\n",
    "    count_2 = {}\n",
    "    for support in supports:\n",
    "        count_2[support] = np.count_nonzero(rvs2 == support)\n",
    "\n",
    "\n",
    "    rvs12 = myuniform.rvs(size = n12)\n",
    "    count_12 = {}\n",
    "    for support in supports:\n",
    "        count_12[support] = np.count_nonzero(rvs12 == support)\n",
    "\n",
    "\n",
    "    rvs22 = myzipf.rvs(size = n22)\n",
    "    count_22 = {}\n",
    "    for support in supports:\n",
    "        count_22[support] = np.count_nonzero(rvs22 == support)\n",
    "    \n",
    "    for event in supports:\n",
    "        if event not in count_1.keys():\n",
    "            count_1[event] = 0\n",
    "        if event not in count_2.keys():\n",
    "            count_2[event] = 0\n",
    "        if event not in count_12.keys():\n",
    "            count_12[event] = 0\n",
    "        if event not in count_22.keys():\n",
    "            count_22[event] = 0\n",
    "        \n",
    "    delta_mle12 = mle_estimator(count_1, count_2,n1, n2, epsilon = epsilon)\n",
    "    delta_poly12 = opt_estimator(count_1,count_2,count_12,count_22,n1,n2,n12,n22, n_mean, epsilon = epsilon, c_1 = 4, c_2 = 0.1, c_3 = 1.5)\n",
    "    return [delta_mle12,delta_poly12]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [03:36<00:00,  2.18s/it]\n",
      "100%|██████████| 100/100 [08:10<00:00,  4.38s/it]\n",
      "100%|██████████| 100/100 [09:02<00:00,  4.80s/it]\n",
      "100%|██████████| 100/100 [01:03<00:00,  1.58it/s]\n",
      "100%|██████████| 100/100 [01:09<00:00,  1.54it/s]\n",
      "100%|██████████| 100/100 [01:02<00:00,  1.60it/s]\n",
      "100%|██████████| 100/100 [01:14<00:00,  1.40it/s]\n"
     ]
    }
   ],
   "source": [
    "n_list = [100, 316, 1000,3162, 10000,31622, 100000]\n",
    "r_runs = 100\n",
    "\n",
    "final_results = []\n",
    "for n in n_list:\n",
    "    results = []\n",
    "    for i in tqdm(range(r_runs)):\n",
    "        results.append(test_once(n))\n",
    "    results = np.array(results,dtype = np.float32)\n",
    "    MSE_mle12 = np.mean((results[:,0]-true_delta12)**2)\n",
    "    MSE_poly12 = np.mean((results[:,1]-true_delta12)**2)\n",
    "    std_error_mle12 = np.std((results[:,0]-true_delta12)**2)/np.sqrt(r_runs)\n",
    "    std_error_poly12 = np.std((results[:,1]-true_delta12)**2)/np.sqrt(r_runs)\n",
    "    \n",
    "    output = [MSE_mle12,MSE_poly12,std_error_mle12,std_error_poly12]\n",
    "    final_results.append(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_results = np.array(final_results)\n",
    "\n",
    "np.save(\"synthetic_results.npy\", final_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a1624bc88>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEeCAYAAACdYvI/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VNXWwOHfSiOQQCiB0DvSkRKkQ4J0QQULIKJ0UcGCeG3c79q9NlQuFpCqIqgISJMiGHoXkC6olFAFpYSeZH9/nIkMMSGZZCZnZrLe55nHnDOnrLCdWdnl7C3GGJRSSqnMCrA7AKWUUr5FE4dSSimXaOJQSinlEk0cSimlXKKJQymllEs0cSillHKJJg7lc0TkExH5t9P2wyJyXEQSRKSInbHlFBHZISIxmTy2qohsFpFzIvKYY187EZnl0SCv3X+GiHTIiXupnCH6HIfKaSJigCrGmH1O+14EKhtj7nfxWsHAWaCxMWarWwO1iYg8DzyfejeQD3jQGPOZi9cbD5w1xjzptG8jMMQYszadcwT4LzDAsWs88IxJ5wtDRIoCHwCdAAPMN8b0crx3C/CxMaaBK3Er76U1DuXrooBQYIerJ4rF6z4DxpjXjTHhzi/gPWAn8G0WLlkOp38fEWkIRKSXNBwGAXcCNwN1gM7AQzc4fgZwzHGvYsA7Tr/PeqCAiERnIXblhbzuQ6OUiMSISLyIPCUiJ0TkqIj0dXp/koi8KiI3AXscu0+LyFLH+01FZIOInHH8t6nTuXEi8pqIrAIuABUd+14VkdWO5q45IlJERKaIyFnHNcqnE+sCERmSat9WEenmSEzvOX6HMyLys4jUysK/RyfgMeBuY8x5x779ItLG8fOLIjJdRL5yNEf9JCI3O95bCsQCox2/201AR2BZBrd9EHjXGBNvjDkMvAv0SSe+dkAZ4GljzBljzFVjzOZUh8UBt7n6uyvvpIlDeaviQARQCugPfCgihZwPMMb8AtR0bBY0xrQWkcLAPGAUUAQYCcxL1ffRG+sv6vzAAce+Ho79pYBKwBpgIlAY2AX8J504vwR6pmyISA2sv7rnAe2AlsBNQEGgO3DKlX8ER8L6HBhkjNl1g0PvAL5xxPslMEtEgo0xrYEVWM1S4Y5/s9pcS7jpqQk4N/1t5dq/dWqNHdebLCKnHIm2VapjdmHVXpQf0MShvNVV4GXHX6/zgQSgaibOuw3Ya4z53BiTaIyZCuwGujgdM8kYs8Px/lXHvonGmF+NMWeA74FfjTE/GGMSsb6Q66Vzv5lAXREp59juBcwwxlx2/A75gWpY/Ym7jDFHM/sPICJ5HPeeYoyZlsHhm4wx0x2/z0is5rvG6RxbEDiXwfXCgTNO22eAcEffR2qlsZLkj1gJ/13gOxGJdDrmnOO+yg9o4lB2SAKCU+0LxvqiTXHK8aWd4gLWl1lGSnKtFpHiAFZNIsWhNM477vTzxTS207y3MeYcVu2ih2NXD2CK472lwGjgQ+C4iIwVkQKZ+B1SfAAkAk9l4ti/fydjTDIQj/VvkZa/sBIaYHXGO5qxEkTkE8fuBMA51gJAQjqd4xeB/caY8Y5EP80RTzOnY/IDpzPxeygfoIlD2eEgUD7Vvgr88ws/K45gNRU5Kwscdtp291DCqUBPEWkC5MX6y9u6kTGjHKOJamI1WT2dmQuKSG/gLuBep1rRjZRxOjcAqxZwJJ1jf3bEkhKjc2f8YMfuHVzftHQz6Q9A+JmM/02rc33Tl/JhmjiUHb4CRohIaREJcHTydgGmu+Ha84GbROQ+EQkSke5ADWCuG659o3uWA14GvnL8xY+INBSRRo4hw+eBS1i1rRtydKB/BPQyxqRVO0pLA0eHfBDwBHAZSG/U1HwgdR9Eap8Bw0SklIiUxKr1TErn2JlAIRF5UEQCReRurBreKqdjWmE1ASo/oIlD2eFlYDWwEqvZ5C2sL8nt2b2wMeYU1tDRp7A6ov8FdDbGnMzutW9wz8tYw1HbYHVMpygAfIr1Ox5wxPMO/N08lN4X6TAgDJjh1ISU8kr9fEeK77A63//C6uTvll5NxRjzE3BGRBrd4NcaA8wBtgHbsZrjxqS86YilheN6fwK3A8Ox+kKeBe5I+Td3DP897xiWq/yAPgColI/LysOTjiG0jxhj7vRYYNfu9S0w3jHIQfmBILsDUErlPGPMImBRDt3rrpy4j8o52lSllFLKJdpUpZRSyiVa41BKKeUSv0ocItLF8ZBVl4yPVkoplRV+2VQVGRlpypcvn6Vzz58/T1hYmHsDUtmm5eJ9tEy8T3bLZNOmTSeNMUUzOs4vR1WVL1+ejRs3ZuncuLg4YmJi3BuQyjYtF++jZeJ9slsmIpKp2Rv8qqlKKaWU52niUEop5RK/ShwpneNnzpzJ+GCllFJZ4ld9HMaYOcCc6OjogXbHopTyXVevXiU+Pp5Lly7ZHYpLIiIi2LXrRut9WUJDQyldujTBwalXN8gcv0ocSinlDvHx8eTPn5/y5cuT9tpV3uncuXPkz5//hscYYzh16hTx8fFUqFAhS/fRxAFEv7qYkwlXru1YMA+AyPAQNo5oa1NUSim7XLp0yeeSRmaJCEWKFOGPP/7I8jX8qo8jq65LGpnYr5Tyf64mje5j1tB9zBoPReNe2U2ImjiUUsoLiQi9e/f+ezsxMZGiRYvSuXNnACZNmsSQIUP+cV758uWpXbs2devWpW7dujz22GNuj02bqpRSKhtSN3WXf9Y9Td1hYWFs376dixcvkjdvXhYvXkypUqUyde6PP/5IZGRklu+dEb+qcXhiOO59n65l04G/3HY9pZR/8WRTd8eOHZk3z0pEU6dOpWfPntm+pjv4VY3DE8Nxfzl+jrs+Xk3rasV4qt1N1CwZ4a5LK6V8wEtzdrDzyNksnZten0eNkgX4T5eaGZ7fo0cPXn75ZTp37szPP/9Mv379WLFiRYbnxcbGEhgYCMCDDz7Ik08+6VrgGfCrxJFVkeEhaf51EBkewvJ/xTJp9X7GLPuN20at5LbaJXiybRUqF7vxkDellMquOnXqsH//fqZOnUqnTp0yfZ6nm6o0ccB17ZBpTRL2SExlejUqx/iVvzN+xW98v/0od9YrxRO33kTZIvlyOFqlVE7KqGaQ0qeRlq8eapLt+99+++0MHz6cuLg4Tp06le3ruYMmjkyKyBvMsLY30adpeT5Z9iuTV+9n9pYjdG9YhqGtq1A8ItTuEJVSfqhfv35ERERQu3Zt4uLi7A4H8LPO8ZxQOCyE5ztVZ/m/Yul5S1m+3niIlm//yCtzd3Iq4bLd4SmlclhkeIhL+11VunRpHn/88TTfmzRpEqVLl/77dfjwYcDq40gZjvvAAw+4JQ5nWuPIoqgCobxyZy0GtazIqCV7mbjqd6auP0i/ZhUY2LIiEXmzNgeMUsq3pDR1p3SEu6N5CiAhIeEf+2JiYv5uSu/Tpw99+vS57v1z586xf/9+t9z/RjRxZFOZwvl4+56bGRxTifcW/8LoH/fx2Zr9PNSqEn2alicsj/4TK5UbuCth+AJtqnKTSkXDGX1ffeY/1oJbKhTm7YV7aPnWj4xb8RuXribZHZ5SSrmN1ycOEakoIuNFZLrdsWRGjZIFGPdgQ2Y+0pTqJQrw6rxdxLwdx5R1B7ialGx3eEoplW0eTRwiMkFETojI9lT7O4jIHhHZJyLP3ugaxpjfjDH9PRmnJ9QrW4gvBjRi6sDGlCqUlxdmbufWd5cx46d4kpKN3eEppVSWebrGMQno4LxDRAKBD4GOQA2gp4jUEJHaIjI31auYh+PzuCaVijB9cBMm9mlI/tAghn29lfbvL2f+tqMkawJRSvkgMcazX14iUh6Ya4yp5dhuArxojGnv2H4OwBjzRgbXmW6MufsG7w8CBgFERUU1mDZtWpbiTUhIIDw8PEvnZiTZGDYdT2Lm3iscOW8oVyCAblWCqRMZ6Jfz/ruTJ8tFZY0/l0lERASVK1d26Zy8X1lfTxe729eqnpSU9PdUIxnZt28fqef1i42N3WSMic7oXDuG/JQCDjltxwON0jtYRIoArwH1ROS59BKMMWYsMBYgOjrapH76O7PSenLcnVoDw5IN3205zPs/7OW9TRdoUK4QT7W7iaaVPDdFgK/zdLko1/lzmezatSvDlfT+IdD6OnX5vPQuFxhI7dq1SUxMpHr16kyePJl8+dKfqSI8PJyjR49m+v6hoaHUq1cvS7HZ0Tme1p/W6VZ7jDGnjDGDjTGVMlErcfvsuJ4QGCB0q1+aJU+14rWutTj810Xu+3QdvcatZfNBnYlXKZ+UeAWO/Qznjrvlcnnz5mXLli1s376dkJAQPvnkE7dc1x3sSBzxQBmn7dLAEXdc2BgzxxgzKCLCN2awDQ4MoFejcsQ9HcO/O9dg99FzdP1oNQMmb8jybJxKKZucOQiXz8KyN91+6RYtWrBv3z4ARo4cSa1atahVqxbvv//+P47t3bs333333d/bvXr1Yvbs2W6Nx46mqg1AFRGpABwGegD3uePCItIF6OJq26TdQoMD6d+8Aj0alnHMxPsrnUat4LY6JRjW9iYqFfXPdmSlfML3z8KxbTc+JvEKJByzft400To+8AZTjhSvDR3/m6nbJyYm8v3339OhQwc2bdrExIkTWbduHcYYGjVqRKtWra5rchowYADvvfced9xxB2fOnGH16tVMnjw5U/fKLE8Px50KrAGqiki8iPQ3xiQCQ4CFwC7ga2PMDnfcz9dqHKmF5Qni0djKrHimNUNiK/Pj7hO0HbmM4d9s5dCfF+wOTymVnjMHr/1sDJw+mP6xmXTx4kXq1q1LdHQ0ZcuWpX///qxcuZKuXbsSFhZGeHg43bp1+8f6HK1atWLfvn2cOHGCqVOnctdddxEU5N46gkdrHMaYNJerMsbMB+Z78t6+LCJvMMPbV6VPs/J8Evcrn609wHdbDtOjYVmGtK5MVAGdiVepHJNRzeDcMfjgZqcdBi6dhrsnQP6oLN82pY/DWWZHwfbu3ZspU6Ywbdo0JkyYkOUY0uP1T467wlc6xzMrMjwPIzrXYPnTsdwbXYap6w/S8q0feW3eTv48n/1lKZVSbrDsLTCpZoUwyR7p62jZsiWzZs3iwoULnD9/npkzZ9KiRYt/HNenT5+/+z9q1sx4pUFX+VXi8PWmqvQUjwjlta61WfpUDJ3rlGT8yt9p8eZSRi7aw9lLV+0OT6ncLX49JKX6Qy7pirXfzerXr0+fPn245ZZbaNSoEQMGDEhzSG1UVBTVq1enb9++bo8B/Gx2XF/tHM+sskXy8e69N/NwTEXeW7yXUUv3MXnNAR5qVZE+TcuTL8SvilMp3zB4pfXfibdZ/+2b/oqArkhrWnWAYcOGMWzYsDSPP3fuHAAXLlxg79699OyZZm9BtvnVN40xZg4wJzo6eqDdsXhS5WL5+bBXfR4+fIaRi3/hrQV7mLDydx6JqcxHcfvSXT/deYlcpZSbuSlhZNcPP/xAv379GDZsGJ5qffGrxJFt545Rd/Pz0GBmtjq1ckqtUhFM6NOQTQf+5J2Fv/Dy3J3pHptWMlFK+Z82bdpw8GD2R3XdiF/1cWS7c/yT5kSc2eGRTi1PalCuMFMHNebLAenO3KKUUm7jV4kjW53j547B+T+s+VA2f+62aQNyUtPKOteVUu7i6Qlg7ZTd382vEke2LHuLv6fRSroC0z0zGsFOJxMu2x2CUj4hNDSUU6dO+WXyMMZw6tQpQkOz/jyY9nGAVdvYMoXr5lo8sAoWPAftXwc/mfI89u04Hm9ThQealCckSP9mUCo9pUuXJj4+nj/++MPuUFxy6dKlTCWE0NBQSpcuneX7+FXiyPJw3LQe4JEAWPsRXDoDXT6AwGC3xelJkeEhaXaEF8oXTJ3SBXl13i6+XH+Qf3euQWxVn18nSymPCA4OpkKFCnaH4bK4uLgsT5XuCr9KHFkejpvWAzwmGcKirJrIuaNw72eQxz3z7HtSRkNuf9x9glfm7qTvxA3EVi3KvzvXoKJOoqiUcoFfJY4sS3mAhzQWp9n8Bcx+DCZ2hF7TIX/xnI/PjWKrFaNZ5Ugmr97PqCV7affecvo2K8/QW6tQINQ3alVKKXtpQ3dG6t0P930Np36DcW3hjz12R5RtIUEBDGxZkaXDY7i7QWnGrfyd2LfjmLb+IEm6DrpSKgOaODKjShvoOx8SL8H4dnBgjd0RuUXR/Hn47111mDOkORWLhvHsjG3c8eFKNuz/0+7QlFJezK8Sh0dnxy1ZFwb8AGFF4bM7YMcs99/DJrVKRfD1Q00Y1bMepxKucM8naxg6dTOHT1+0OzSllBfyq8Th8dlxC5WD/ougZD34pg+s+cgz97GBiHD7zSVZ+lQMj99ahUU7jnHru3G8/8MvXLySZHd4Sikv4leJI0fkKwwPzILqnWHhc7DwBUhOzvg8H5E3JJAn297EkqdacWv1KN7/YS9tRi5j7s9H/PJhKKWU6zRxZEVwXrhnMjQaDGtGw7f94Oolu6Nyq9KF8vHhffX5alBjIvIGM+TLzXQfs5bth/1jkSylVNZp4siqgEDo8F9o9yrsmAlfdIOLf9kdlds1qliEOUOb83rX2uz7I4Euo1fy3IyfdfoSpXIxTRzZIQJNh8Jd4yF+A4xvD6cP2R2V2wUGCPc1KsuPw2Po16wC32yMJ/btOMat+I0rif7TTKeUyhxNHO5Q+264f4Y159W4NnD0Z7sj8oiIvMH8u3MNFjzRkvrlCvHqvF10+GA5P+45YXdoSqkc5FeJw6PDcTNSoQX0XwgBQTCxE/y6NOdjyCGVi4Uzud8tTOzTEAz0nbiBvhPX89sfaS91qZTyL36VODw+HDcjxarDgMXWsN0p98DWafbEkUNiqxVjwRMteaFTdTbu/4t27y3ntXk7OXvpqt2hKaU8yK8Sh1coUNJ6yrxcM5j5ECx/B/x4GKtOX6JU7qOJwxNCI6wJEet0h6WvwNwnISnR7qg8KmX6ktmPNqdCpE5fopQ/08ThKUEh0HUMNB8GmybCV/fDlfN2R+VxtUtH8M1gnb5EKX+micOTRKDNf+C2d2HvQpjcBRJ8a0WxrNDpS5Tyb5o4ckLDAdB9ChzfCePbwqlf7Y4oR+j0JUr5J00cOaVaJ3hwDlw+ayWP+I12R5RjdPoSpfyLJo6cVKYh9F8MeQrApM6we77dEeWo9KYvOaXTlyjlU3wicYjInSLyqYh8JyLt7I4nW4pUspJHserwVS/YMN7uiHJUWtOXxLyj05co5Us8njhEZIKInBCR7an2dxCRPSKyT0SevdE1jDGzjDEDgT5Adw+GmzPCi0KfuVClHcwbBkte9utnPdJy3fQlZXX6EqV8SU7UOCYBHZx3iEgg8CHQEagB9BSRGiJSW0TmpnoVczp1hOM83xcSZnWYN+gDK96FmYMh8YrdUeU4nb5EKd8jOTG6RUTKA3ONMbUc202AF40x7R3bzwEYY95I53wB/gssNsb8kM4xg4BBAFFRUQ2mTcvadB8JCQmEh4dn6dwsMYayB7+h4u9T+LPQzeyo+SxJQfly7v5eJDHZsPhAIrN/vcKVJAgUuJJG61WBEBjVOiznA1TXyfHPispQdsskNjZ2kzEmOqPjgrJ8h+wpBTjPPx4PNLrB8UOBNkCEiFQ2xnyS+gBjzFhgLEB0dLSJiYnJUmBxcXFk9dysi4UtzSk8eygt9r4Gvb6xpi7JhdoAw89d5t1Fe5i2Ie0p6s9ewYYyUqnZ81lRN5JTZWJX57iksS/dqo8xZpQxpoExZnBaSePvi9o5O2521b0P7vsa/toP49rCiV12R2SblOlLlFLeya7EEQ+UcdouDRzJ7kVtnx03uyrfCn2/h+REmNAe9q+0OyKllPoHuxLHBqCKiFQQkRCgBzA7uxf16RpHihJ1rKnZw4vD511h+7d2R+SVluw6bncISuVaOTEcdyqwBqgqIvEi0t8YkwgMARYCu4CvjTE7snsvn69xpChYFvotgFLRML0frB6d64brZqT/5I08/MUmjp25ZHcoSuU6Hu8cN8b0TGf/fCB3PTrtinyFofdMa02PRS/AmXho/xoEBNodWY6JDA/hZMI/hyhHhofQt1kFRi3Zy4q9Jxne7iZ6NylPYEBaXWdKKXeza1SVR4hIF6BL5cqV7Q7FPYJD4e6JsKgUrP0Qzh6GbmMhOK/dkeWIjSPa/v1zWqNFOtcpwb+/28GLc3YyY/NhXu9am1qlfLy2qZQP8IkpRzLLb5qqnAUEQIfXof3rsGsOfHYnXNDFkQDKFQljct+GjOpZjyOnL3H76JW8PGcnCZf9e9EspezmV4nDLzrH09PkUbhnIhzZbI24+uuA3RF5hZS1P5Y81Yqet5RlwqrfaTtyGYt2HLM7NKX8ll8lDr+scTir2RUemAUJJ6yp2Y9utTsirxGRN5jXutbm24ebEpE3mEGfb2LgZxs5oisPKuV2fpU4coVyTaH/IggMgYmdYF+aM7DkWg3KFWLO0OY817EaK/b+QZuRyxi34jcSk3TmXaXcxa8Sh183VTkrWtWamr1wBZhyL2yeYndEXiU4MICHWlVi8ZOtaFShMK/O28UdH65i66HTdoemlF/wq8Th901VzgqUgD7zoUJL+O4RWPaWPuuRSpnC+ZjQpyEf9arPH+cuc+dHq/jPd9s5d+mq3aEp5dP8KnHkOqEFrAkRb74PfnwN5jwGSTqiyJmI0Kl2CX54qhUPNC7HZ2sP0GbkMuZvO6rrniuVRZo4fF1gMNz5EbR8Gn76DKb1hMsJcO4YTOwI53RqDoACocG8dEctZj7SjCJheXhkyk/0n7yRQ39esDs0pXyOXyWOXNPHkZoItB4Bnd+3Ossnd4ZPWsCB1bDsTbuj8yp1yxRk9pBmjLitOmt/O0W795YzZtmvXNXOc6Uyza8SR67q40hLdF/oMRVO7IbzjiVYt0zRWkcqQYEBDGhRkcXDWtGsciRvfL+bLv9byU8H/7I7NKV8gl8lDgVU7QCVr03VQeIlGNMCfnwD9q+CxMv2xeZlShXMy7gHoxnTuwGnL1zlro9XM2LWNs5c1M5zpW7Er+aqUlh9G/sWXb8v4YTVZLXsvxCUF8o2skZjlW8JJetBYO7+36B9zeI0qxzJyEW/MGn17yzccZz/61yDznVKYK1arJRylru/MfzRsrfApGqvDwyGOt2haif4fTnsXwFLXrbeC8kP5ZpYiaRCS4iqbc2PlcuE5wni/7rUoGu9Ujw/cxtDp27mm03xvHpHLcoWyZ1rwCuVHr9KHH43O25WxK+HpFRTkSddgaNb4I7RUK2Tte/8SSuB/L7CSiZ7HbWU0IJQvvm1RFK0mtX5nkvULh3BrEeb8dma/byzcA9t31vGY7dWYWCLioQE5b6EqlRa/CpxGGPmAHOio6MH2h2LbQZncrnZsEhr7quaXa3ts0cdiWSZlUx2z3UcVxTKt7iWSApX9PtEEhgg9G1WgQ61ivPS7J28vXAP3205zGtda9OwfGG7w1PKdn6VOFQ2FCgBde61XgB/7bcSyH5HjWTHDMdxpRz9I45kUrBMupf0dSUi8vJJ7wb8sPM4/5m9g3s+WUPPW8rwTIdqFMwXYnd4StlGE4dKW6Hy1qt+b2sqk1O/OmojjmatrVOvHVehJVRoZSWT/FE2Bu0ZbWpE0aRSET5YspfxK39n0Y7j/LtzDe6oW1I7z1WupIlDZUwEIitbr4b9ITkZ/thlJZHfV8CO76yn1gEiq0KFFtdqJfn8o2knLE8Qz3eqzh11S/L8zO088dUWpm+K55U7a1EhMszu8JTKUZo4lOsCAiCqpvVq/DAkJ8Gxnx2JZDlsmQobxlnHRtW+1j9SrgmE+vbDmTVLRjDj4aZ8ue4Aby3YQ/v3lzMktjIPtapInqDcsx68yt00cajsCwi0ngcpWQ+aPQ5JV+HwT46hv8utJLL2Q5AA65iU2kjZxhDie3+tBwYIvZuUp13N4rw8dycjF//yd+d544pF7A5PKY/zq8Shw3G9RGCw9ZBh2UbQ6mm4eskaJpwy9Hf1/2DlexAQDKWjr9VISjeEoDz/vN65Y9Td/Dw0mOlVfShRBUL58L763N3gBP+etZ0eY9dyT4PSPNepOoXDtPNc+S+/Shw6HNdLBYdeSw68YM3ee2jttT6S5W9bT7YHhUKZRo4+klaOp9qDYUwrIhKOWcd0Hmn3b/MPsVWLsfjJVnywZC/jVvzGD7uO88JtNbirfintPFd+ya8Sh/IRecKhchvrBXDxtDWTb8rQ36WvAq9CSDiUrA8JxxGwJmxs9YxX1TpS5A0J5NmO1bizXkmen7GN4d9sZfqmQ7zWtTaViobbHZ5SbqWJQ9kvb0HrifbrnmpfaSWRbdMBx4JLyUleW+tIUa14AaYPbsrUDQd58/vddHx/BUGBwoUrSf84NjI8hI0j2qZxFaW8m86hoLxPWCTUvBNa/QuSnGbzTb5qDfs9d8y+2DIhIEDo1agcS56KoWPt4mkmDYCTCVfS3K+Ut9PEobxXWhM2Jl+FCR3g0ll7YnJB0fx5+KBHPbvDUMrtNHEo75XWhI0Af/0OY1vB0a05H5NSSvs4lBdzmrAxLi6OmJgYa+PAGpjeD8a1hQ5vQHQ/n5148fjZS0QVCLU7DKVcojUO5XvKNbGSSoUWMG+YlUR8oOkqLW1GLuOrDQcxxtgdilKZdsPEISL3O/3cLNV7QzwVVKr7VBeRT0Rkuog8nBP3VD4grAjc9w20eRF2fufVTVeR4Wk/DFgoXzA1ShTgmW+30WvcOg6eupDDkSmVNRk1VQ0DvnD8/D+gvtN7/YDRNzpZRCYAnYETxphaTvs7AB8AgcA4Y8x/07uGMWYXMFhEAoBPM4hX5SYBAdD8SSjT2Knp6nWI7u9VTVc3GnKbnGyYuuEgb8zfTfv3lzO8fVX6NC1PYID3xK9Uahk1VUk6P6e1nZZJQIfrThIJBD4EOgI1gJ4iUkNEaovI3FSvYo5zbgdWAksycU+V2/zddNUS5j0F0/v6TNNVytDdxcNa0qRSEV6Zu5O7Pl7NL8fP2R2aUumSG7WtishPxpj6qX9Oa/sG1ygPzE2pcYhIE+BFY0x7x/ZzAMaYNzJxrXnGmNvSeW8QMAggKiqqwbRp0zIrgMivAAAf7ElEQVS6XJoSEhIID9cnfb1NpsrFJFPm0Ewq/vYFl0KLsaPmMyTkr5gzAbqBMYa1R5OYsusyFxPh9krB3FYxmCAvrX3oZ8X7ZLdMYmNjNxljojM6LqOmqmoi8jNW7aKS42cc21n9RJYCDjltxwON0jtYRGKAbkAeYH56xxljxgJjAaKjo83fI3BcdN3oHeU1Ml8ureHAfeSd3o/oLc84Rl15V9PVjcQCgxIu8+KcnczceoRdCaG8eVcdbi5T0O7Q/kE/K94np8oko8RR3QP3TOsTnG61xxgTB8Rl6sI6O66Ca01XMx+ymq72r4QuoyC0gN2RZUqR8Dz8r2c9br+5JCNmbaPrR6sY2KIiT7S5ibwhuuaHst8N+ziMMQecX0ACVgd5pGM7K+IB54WqSwNHsnit6xhj5hhjBkVE+PZiQcoNworAfV87Rl3NhjEtvXbUVXra1ohi8bBWdG9YhjHLf6PjB8tZ+9spu8NSKsPhuHNFJKVvogSwHWs01eci8kQW77kBqCIiFUQkBOgBzM7itVLH20VExp45c8Ydl1O+LmXUVd/5kHgZxrWB9Z9aa6j7iAKhwbzRrQ5fDmhEsoEeY9fywsxtnLt01e7QVC6W0aiqCsaY7Y6f+wKLjTFdsPok+mV0cRGZCqwBqopIvIj0N8YkAkOAhcAu4GtjzI4s/wZOtMah0lS2sdV0VTEG5g+Hb/rAJd/646Jp5UgWPtGSAc0rMHX9Qdq9t5ylu4/bHZbKpTJKHM5/1tyKo3PaGHMOSE7zDCfGmJ7GmBLGmGBjTGljzHjH/vnGmJuMMZWMMa9lNfjUtMah0hVWBHp+BW1egl1zYEwrOLLF7qhckjckkBGda/Dtw03JHxpEv0kbeWLaZv48r7PsqpyVUeI4JCJDRaQrVt/GAgARyQsEezo4V2mNQ91QQAA0f8Jqukq6AuPb+lzTFUC9soWYO7QFj99ahXnbjtJm5DJmbz2i05aoHJNR4ugP1AT6AN2NMacd+xsDEz0Yl1KeU7YxPLTCqenqQZ9rugoJCuDJtjcxZ2hzyhTKy2NTNzPws00cO3PJ7tBULpDRqKoTxpjBxpg7jDGLnPb/aIx5x/PhuUabqlSmXdd0NdcadeVjTVdgrTg445FmvNCpOiv3/UHbkcuYul4nTVSeldGoqtk3euVUkJmlTVXKJdc1XV312aarwABhYMuKLHi8JTVLFeC5Gdu479N1HDh13u7QlJ/K6AHAJlhPeU8F1pG5+amU8i0po65mPmQ1Xe1fAbf/D0J96w+Q8pFhfDmgMdM2HOL1+busSRPbVaVvswo6aaJyq4z6OIoDzwO1sGazbQucNMYsM8Ys83RwrtKmKpVl+QpbTVdtX3Zqutpsd1QuCwgQ7mtUlsXDWtK0UiSvzttFt49Xs+eYTpqo3CejPo4kY8wCY8yDWB3i+4A4ERmaI9G5SJuqVLYEBECzx6Hv946mq3awbqzPNV0BlIjIy/gHo/mgR10O/XmBzv9bwfs//MKVxAxH0SuVoQxXABSRPCLSDWtdjkeBUcAMTwemlG3KNnI8MBgL3z8NXz/gc6OuAESEO+qWYvGTLelYqwTv/7CX20evZOuh0xmfrNQNZNQ5PhlYjfUMx0vGmIbGmFeMMYdzJDql7JKvMPScZjVd7Z7ns01XYE2aOKpnPcY9EM3pC1fp+tEqXp+/i4tXkuwOTfmojGocvYGbgMeB1SJy1vE6JyJet1KO9nEot7qu6SrRp5uuANrUiGLRsJZ0b1iWsY5JE9f8qpMmKtdl1McRYIzJ73gVcHrlN8Z43RzV2sehPKJsIxi8wuebriBl0sTafDmwEQbo+elanp+5jbM6aaJyQYZ9HEopnJquXrnWdHX4J7ujyrKmlSJZ8HhLBraowLT1B2k3UidNVJmniUOpzAoIgGaPpWq6GuOzTVd5QwJ54bYazHikGRF5g+k3aSOPT9vMqYTLdoemvJwmDqVcldJ0Vak1fP8v+Lo3XPTdkUp1yxRkztDmPNGmCvO3HaXte8t10kR1Q36VOLRzXOWY65qu5vt801VIUABPtLmJuUNbUKZwPsekiRt10kSVJr9KHNo5rnKUc9NVcpLPN10BVC2enxkPN2XEbdVZue+kTpqo0uRXiUMpW/hZ01VggDCgRUUWPtGSWqUidNJE9Q8ZTXKolMqMlKarNaNhyUtW09U9E6FUA7sjy7JyRcL4cmAja9LEedakiYEinHd+cHDBPAAiw0PYOKKtTZGqnKY1DqXc5R9NV+1h7Sc+3XQlIvS8pSyLh7WieeXI65OGk5MJunxtbqKJQyl3K3OL1XRV+VZY8Ax8db/VdHXuGEzsCOd873mJ4hGhfPpAtN1hKC+hiUMpT0hpumr3KvyywGq6+qgpHFgNy960O7osEdE1PZTFrxKHDsdVXkUEmg61mq4Sr8BFx7xQW6b4ZK0jI7/9kWB3CCqH+FXi0OG4yiuVucVqtkqReNlqwvIz7d9fzpsLdnP+cqLdoSgP86vEoZRXOncMtk932mFgx0yY8wRc9q2V+SLDQ9LcXzgshDvqluLjuF+59d1lzNEnz/2aDsdVytOWvQUm1cp7EgCbJsLeRdDxTajW2Wra8nLOQ27j4uKIiYm57v2et5TlP7O3M3TqZqasO8BLt9eiavH8ORyl8jStcSjlafHrISnVcFWTDIUrQWhBa9TV1B7w1wF74nOjBuUK8d2jzXmtay12HztHp1EreGXuTp223c9ojUMpTxu8Mv33kq7Cuk/gxzfgw0bQ6l/QZAgEpd0k5AsCA4RejcrRqVYJ3l60hwmrfue7LUd4rmM1utYrRUCA99es1I1pjUMpOwUGWyOvhqy3OtCXvARjWljDdn1cobAQXu9am9mPNqd0obw89c1W7hmzhu2HddSjr9PEoZQ3iCgNPaZYz35cuWA9KDjrUTjv+0u71i4dwYyHm/L23XXYf/I8XUavZMSsbZy+oE+b+ypNHEp5k6od4dG10OwJ+HkajG4AP30OyckZn+vFAgKEe6LLsHR4DH2almfq+kPEvhPHl+sOkpSso698jU8kDhEJE5FNItLZ7liU8riQMGj7Ejy0AopWg9lDYFInOL7T7siyLSJvMP/pUpO5Q5tTJSo/z8/cRtePVrH54F92h6Zc4NHEISITROSEiGxPtb+DiOwRkX0i8mwmLvUM8LVnolTKS0XVgD7z4fbR8Mduq+9j8f/BFd+f3rx6iQJ8NagxH/Soy/Gzl+j60Wr+NX0rJ3XZWp/g6RrHJKCD8w4RCQQ+BDoCNYCeIlJDRGqLyNxUr2Ii0gbYCfjfHA1KZSQgAOr3hiGb4OYesOoDa/TVnu/tjizbRIQ76pZiyVMxPNSqIjN+OkzsO3FMWvU7iUm+3TTn78TTT3eKSHlgrjGmlmO7CfCiMaa9Y/s5AGPMG+mc/xoQhpVkLgJdjUn9NBWIyCBgEEBUVFSDadOmZSnehIQEwsPDs3Su8hwtF0vE6R3c9MsnhF04yB+RjdhXeSCXQ4vaEou7y+RIQjJTdl1mx6lkSocLvWvkoWrhQLddPzfIbpnExsZuMsZkOA2yHYnjbqCDMWaAY7s30MgYMySD6/QBThpj5mZ0z+joaLNx48YsxZvW07DKflouThKvwNoPIe5N6wn02Oeg0WBraG8O8kSZGGNYuOM4r8zdyeHTF7mzbkme61SdqAKhbr2Pv8pumYhIphKHHZ3jaT39k2H2MsZMyihp6Oy4KlcICoHmT8Kj66BCC1g0AsbGwKH1dkeWbSJCh1rF+WFYKx5rXZn524/R+p04xi7/lSuJ2nzlLexIHPFAGaft0sARd1xYZ8dVuUqhctZzH92nwMW/YHxbmP0YXPjT7siyLW9IIMPaVWXxky1pUqkIr8/fTccPlrNy70m7Q1PYkzg2AFVEpIKIhAA9gNnuuLDWOFSuIwLVO8Oj662pSjZ/AaMbwpapPr1kbYpyRcIY92BDJvSJJjHZcP/4dTz8xSYOn75od2i5mqeH404F1gBVRSReRPobYxKBIcBCYBfwtTFmhzvupzUOlWvlCYf2r8FDy6BwBZg1GCZ3gT/22B2ZW7SuFsXCJ1rydPuq/LjnBLe+G8fopXu5dDXtNdCVZ3k0cRhjehpjShhjgo0xpY0x4x375xtjbjLGVDLGvOau+2mNQ+V6xWtDv0XQ5QM4tg0+bgZLXoGrvv8XemhwII/GVmbJUzG0rlaMdxb9Qvv3l7N0t47Uz2k+8eR4ZmmNQymsZz8a9IEhG6H23bDiHevZj72L7Y7MLUoVzMtHvRrwRf9GBAUI/SZtpP+kDRw45fsPRvoKv0ocSikn4UWh6yfw4BwIDIEpd8PXD8BZt4xFsV3zKpF8/3hLXuhUnbW/naLte8t5d9EeLl7R5itP86vEoU1VSqWhQkt4eBW0HgG/LLQ6z9d+DEm+vzZ4SFAAA1tWZOnwGG6rXYL/Ld1Hm5HLWLD9qC5d60F+lTi0qUqpdATlgZZPwyNroGxjWPAsfBoL8ZvsjswtogqE8l73unz9UBPyhwYx+IufeGDCevadSLA7NL/kV4lDKZWBwhWh13S4ZxKc/wPG3QrznoKLp+2OzC1uqVCYuUOb89LtNdl66DQd3l/OG/N3kXDZ92tX3sSvEoc2VSmVCSJQs6v17EejwbBxgtV8tW26Xzz7ERQYwINNy7N0eAx31S/NmOW/0fqdOL7bclibr9zErxKHNlUp5YLQAtDxvzDwR2sFwm/7w+d3wqlf7Y7MLSLD8/Dm3XWY+UhTikeE8vi0LXQfu5ZdR8/aHZrPC7I7AKWUzUrWhQE/WDWPJS/DR42h+TBrPqxg359csF7ZQsx6pBlfbTzEWwt20/l/K+nduBxzth7h1Pl/Ll8bGR7CxhFtbYjUd/hVjUMplUUBgXDLQOvZjxp3wLL/wsdN4NeldkfmFgEBQs9byvLj8Bjuu6Usn63Zn2bSADiZoGuhZ8SvEof2cSiVTfmj4K5x0HsWIPB5V5jeH875x9PZBfOF8MqdtZg9pLndofg0v0oc2sehlJtUioWHV0PMc7BrNoyOhvWfQrJ/PFxXq5R+R2SHXyUOpZQbBYdCzLPwyFooVR/mD7eG7x7ZAp+2pfmKnn5TE0lNR1/dmCYOpdSNFalkNV3dNR7OHLYeHDyxncCkC7DsTbuj84h7x6xh9zEdfZUeTRxKqYyJWBMmDtkAdXrC1QvWUp4/TYbjO+2OLksiw0PS3B+eJ4h9JxK4bdRKXp27Ux8eTINfDccVkS5Al8qVK9sdilL+KW/B64foJidaU7c3eACaDIVI3/ns3WjI7V/nr/DWwj2MX/U7c34+wojbatC5TglE0lr5OvfxqxqHdo4r5WHnjsGWKdfvE4HNX1od6NN6wcF19sTmRoXCQnijW21mPtKMovnzMHTqZnqPX8+vf+jcV+BniUMp5WHL3gKTfP2+gECofQ+0HA77V8KEdjC+HeyaC8nJaV/HR9QtU5DvHm3OK3fUZGu8NffVWwt25/qp2zVxKKUyL349JKV6QC7pChz/2Zq2fdhO6PgWnDsKX/WCDxvCxolw9ZI98bpBYIDQu0l5lj4VQ5ebS/JR3K+0GbmMhTuO5drRV37Vx6GU8rDBK//+MS4ujpiYmOvfDwmDRg9BdH/Y9R2sGgVzn4AfX7u2P1/hnI3ZTYrmz8PIe+vSo2FZ/j1rOw99vonW1YrxYpealC2Sz+7wcpTWOJRS7hcYBLXugkFx1gqEJerC0lfhvVrw/TPw1wG7I8yyWyoUZu5jzRlxW3XW/XaKtu8t44Mf9nLpau5pvtLEoZTyHBFrBcL7p8PDa6x5sDaMh1H1YHo/62FCHxQcGMCAFhVZ8lQMbWtE8d4Pv9D+/eXE7Tlhd2g5wq8Sh85VpZQXi6oBXT+Gx7dCk0fgl0UwthVM7gJ7f/DJtUCKR4Qy+r76fNG/EYEBQp+JGxj8+SaOnL5od2ge5VeJQ4fjKuUDIkpBu1dh2A5o+zKc3AtT7rKeB9kyFRJ9b3ba5lUi+f7xFjzdvipxv5zg1neX8XHcr1xJ9O1RZenxq8ShlPIhoRHQ7HF4/Ge482PAwKzB8MHNVqf6Jd+a8iNPUCCPxlZm8ZOtaF4lkjcX7KbTqBWs/vWk3aG5nSYOpZS9gkKg7n3WbLy9pltzYy3+N7xXExb9G84esTtCl5QpnI9PH4hmQp9oLicmcd+n63h82mZOnPXdIcmpaeJQSnkHEajSFvrMtZazrdwG1oyG9+vAzId9bk6s1tWiWPxkKx67tQrfbzvGre8uY8LK30lM8v3mK00cSinvU6o+3DMRHtsM0f1g5yxrRcIv7obfl/tMR3pocCDD2t7EwidbUq9cIV6eu5Muo1ex6cCfdoeWLZo4lFLeq1B56PQWPLkDYkfA0S3WKKxPY2H7DEjyjZlrK0SGMblvQz7uVZ/TF65w18drePqbrZxKuGx3aFmiiUMp5f3yFYZWT8MT26Dz+3D5HEzvC/+rD+vGwpXzdkeYIRGhY+0S/DCsFQ+1qsjMzYdp/e4ypqw7QFKyb9SgUmjiUEr5juC8EN0XHt0A3adAeBR8/7TVkb70NUj4w+4IMxSWJ4jnOlbn+8dbUL1Efl6YuZ1uH63i5/jTdoeWaV6fOEQkRkRWiMgnIhJjdzxKKS8QEADVO8OAxdBvIZRtCsvfhvdrwZwn4NSvdkeYoSpR+Zk6sDEf9KjL4dOXuOPDVYyYtY0zF67aHVqGPJo4RGSCiJwQke2p9ncQkT0isk9Ens3gMgZIAEKBeE/FqpTyUWUbQ88vHasTdoctX8L/Glhrgxxab3d0NyQi3FG3FEuHt+LBJuX5ct1BWr8bxzcbD5Hsxc1Xnq5xTAI6OO8QkUDgQ6AjUAPoKSI1RKS2iMxN9SoGrDDGdASeAV7ycLxKKV8VWQVuHwVPbocWT1lrg4xvC+Pbw+55Xr02SIHQYF68vSZzhjanXJF8PD39Z+4ds4ZdR73zIUjx9HzyIlIemGuMqeXYbgK8aIxp79h+DsAY80YG1wkBvjTG3J3O+4OAQQBRUVENpk2blqV4ExISCA8Pz9K5ynO0XLyPt5dJYOJFih/7gTKHZhN6+QQX8pbiUJk7OR4VQ3Jg2uuNe4NkY1h5OJGv91zhQiK0KRtE1yoh5A3KeNna7JZJbGzsJmNMdEbH2ZE47gY6GGMGOLZ7A42MMUPSOb8b0B4oCHxsjInL6J7R0dFm48aNWYo3zTUGlO20XLyPz5RJUqL1HMjqUXB0K4QVg0aDrq0N8mlbOLkbhmyE/FF2R/u3v85f4e1Fe5i6/iBFw/MwonMNumSw7nl2y0REMpU47OgcT+u3Tjd7GWNmGGMeMsZ0zyhp6Oy4Sql/CAyC2nfDoGXwwGwoUef6tUFO/w6Xz8KyN+2O9DqFwkJ4vau17nmxAnl4bOpmeo1bx74T9q97bkfiiAfKOG2XBtwyGY3OjquUSpcIVGwF938Lg1dB9S6w/lM47xjC+9NnXjkay3nd8+2Hz9Dxg+W8uWA3F67Y9/CjHYljA1BFRCo4+i16ALPdcWGtcSilMqV4Leg2xqqJpEi+CqMbwle9YcdMuHLBvvhS+Xvd8+Ex3H5zKT6O+5W2I5ezYLs96557dM1xEZkKxACRIhIP/McYM15EhgALgUBggjFmhzvuZ4yZA8yJjo4e6I7rKaX82LljsPO7f+7fvwp2zYbgMKjWyVoCt1JrCMqT8zGmEhmeh3fvvZnuDcvwf99tZ/AXmwgOFK4mOSWPBfMcx4awcURbj8Th0cRhjOmZzv75wHx3309EugBdKleu7O5LK6X8zbK3wKQaohsQCDXvtJa43f6tlUC2fQN5IqymrVrdoEIrq9/ERrdUKMycoc2ZvHo/r87bleYxJxM8tyCW1z857grt41BKZVr8ekhK9eWadMXaX7GV9UzI8L3WGiHVOllJ5Itu8O5NMPdJ6zmR5CR7Yufauud2sDdtKqWUXQavzPiYwGBrjZAqbeHqJdj3g1UT2ToNNk6A/CWgxp1Wc1bpaKsDPhfwq8ShTVVKKY8JDrXmx6re2ZqN95cF1tTuGyfAuo+hYFmo2c1qzipex6+TiF8lDu0cV0rliJAwq5ZR6y64dAZ2z7dqImtGw6r3oUhl672a3aBYNbujdTu/ShxKKZXjQiOgbk/rdeFPqy9k+7dW5/uyN6FYTasWUqsbFHZ/n0RkeEiaHeGR4Z6bVsWvEoc2VSmlbJWvMDToY71ShvtunwFLX7FeJes5aiJdIaK0W27pPOQ2p6aB0VFVSinlCfmLQ6OHoP9CeGI7tH3FWit90Qhr4anx7a3VC88dtztSl/lVjUMppbxSwTLQ7DHrdepX2DHDqol8/zQseAbKt7CasqrfbtVavJxf1TiUUsrrFakELZ+GR9bAI2uhxXA4Ew9zHod3qsCUe2DLVLjknWtxgJ/VOLSPQynlU4pVh9YvQOzz1pTvKTWRvYMhMI/1/Eitu+Cm9tZILi/hV4lDh+MqpXySCJSsa73avATxG6wEsmMm7J4LwfmgakcriVRuY/u8WX6VOJRSyueJQJlbrFf71+DAaqsmsmOWNcw3TwRUu81KIhVbWU+3A5w7Rt3Nz0ODmR5fkEoTh1JKeauAQKjQwnp1fAt+X2bVRHbNha1fQt7C1oSMtbrBtwOJSDhmPTvSeaRnw/Lo1ZVSSrlHYLDVTHXnR/D0Xugx1Zru/eevYXIXSDhmLa+6+QuPD/H1q8ShCzkppXKFoDzWjL13j4en90HF2GvvmWSPL4PrV4lDHwBUSuU6l8/CwTXXtpOvwpYpHq11+FXiUEqpXCetBak8XOvQxKGUUr7sRgtSeYiOqlJKKV/mtCCVTnKolFLKK2niUEop5RK/Shw6HFcppTzPrxKHDsdVSinP86vEoZRSyvM0cSillHKJJg6llFIu0cShlFLKJX73AKBjFcCTInIg1VsRQOrhVmntiwROeii8jKQVT05cJ7PHZ3Tcjd5P7z1vLxe7yiSz52TnGF8tE3BPuXiqTDJznKc+K9ktk3KZOsoY41cvYGxm96ezb6O3xe7p62T2+IyOu9H7vloudpVJZs/JzjG+WibuKhdPlUlmjvPUZyWnysQfm6rmuLA/vWPt4q54XL1OZo/P6Lgbve+r5WJXmWT2nOwc46tlAu6Jx1NlkpnjfPqzIo4spRxEZKMxJtruONT1tFy8j5aJ98mpMvHHGkd2jbU7AJUmLRfvo2XifXKkTLTGoZRSyiVa41BKKeUSTRxKKaVcoolDKaWUSzRxZEBE7hSRT0XkOxFpZ3c8CkSkuoh8IiLTReRhu+NR14hImIhsEpHOdseiQERiRGSF4/MS467r5srEISITROSEiGxPtb+DiOwRkX0i8iyAMWaWMWYg0AfobkO4uYKLZbLLGDMYuBfQ4aAe5Eq5ODwDfJ2zUeYuLpaJARKAUCDeXTHkysQBTAI6OO8QkUDgQ6AjUAPoKSI1nA4Z4XhfecYkXCgTEbkdWAksydkwc51JZLJcRKQNsBM4ntNB5jKTyPxnZYUxpiNWQn/JXQHkysRhjFkO/Jlq9y3APmPMb8aYK8A04A6xvAl8b4z5KadjzS1cKRPH8bONMU2BXjkbae7iYrnEAo2B+4CBIpIrv188zZUyMcYkO97/C8jjrhj8bpLDbCgFHHLajgcaAUOBNkCEiFQ2xnxiR3C5VJpl4mir7Yb1QZhvQ1y5XZrlYowZAiAifYCTTl9ayvPS+6x0A9oDBYHR7rqZJo5rJI19xhgzChiV08EoIP0yiQPicjYU5STNcvn7B2Mm5VwoyiG9z8oMYIa7b6ZVyWvigTJO26WBIzbFoixaJt5Jy8X75GiZaOK4ZgNQRUQqiEgI0AOYbXNMuZ2WiXfScvE+OVomuTJxiMhUYA1QVUTiRaS/MSYRGAIsBHYBXxtjdtgZZ26iZeKdtFy8jzeUiU5yqJRSyiW5ssahlFIq6zRxKKWUcokmDqWUUi7RxKGUUsolmjiUUkq5RBOHUkopl2jiUEop5RJNHEoppVyiiUMpDxCRBDdcY7U7YlHK3fTJcaU8QEQSjDHhdsehlCdojUPlCo61sOeJyFYR2S4i3R37ZznWyN4hIoMc+8qLyG4RGec4doqItBGRVSKyV0RuSXXcZBH52bEGer407n2/iKwXkS0iMsaxWltmYktw/Hew49wtIvK7iPyYmes6XX+miLzqWHv6mGOlPqWyTBOHyi06AEeMMTcbY2oBCxz7+xljGmCtXf6YiBRx7K8MfADUAaphrWrXHBgOPO903arAWGNMHeAs8IjzTUWkOtZa9c2MMXWBJP65amF6sQFgjPnEcW5DrOmzR2byuilqAaeNMS0c8emqiSpbNHGo3GIb0EZE3hSRFsaYM479j4nIVmAt1noGVRz7fzfGbHOsYrcDWGKsdt1tQHmn6x4yxqxy/PwFVnJxdivQANggIlsc2xUzGVtqHwBLjTFzMnldHDWgCOA9x64g4HQ611cqU3QFQJUrGGN+EZEGQCfgDRFZBCzHWha4iTHmgojEAaGOUy47nZ7stJ3M9Z+b1J2EqbcFmGyMec6V2IwxL193EWs51nJYU2dn6roONYFNxpgkx3YdYHsG5yh1Q1rjULmCiJQELhhjvgDeAepj/SX+lyNpVAMaZ+HSZUWkiePnnsDKVO8vAe4WkWKOOAqLSLlMxOb8fgOsJrL7ndbxzvC6DrWALU7bdYCfXfwdlbqO1jhUblEbeFtEkoGrwMNYTUSDReRnYA9Wc5WrdgEPisgYYC/wsfObxpidIjICWCQiAY57PwocyCA2Z0OAwsCPIgKw0RgzIBPXTbn2OqftWmiNQ2WTDsdVKotEpDww19GhrVSuoU1VSimlXKI1DqWUUi7RGodSSimXaOJQSinlEk0cSimlXKKJQymllEs0cSillHKJJg6llFIu0cShlFLKJZo4lFJKueT/AWqihUtQzQXyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.errorbar(n_list, final_results[:,0],final_results[:,2],marker='s',label=\"MLE\")\n",
    "plt.errorbar(n_list, final_results[:,1],final_results[:,3],marker = '^' ,label=\"Poly\")\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.ylabel(\"MSE\")\n",
    "plt.xlabel(\"sample size $n$\")\n",
    "plt.title(\"Uniform vs. Zipf(-0.6)\")\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": false,
   "bibliofile": "ref.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}