{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from __future__ import division, print_function\n", "\n", "from collections import defaultdict\n", "import os\n", "\n", "import numpy as np\n", "\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Approximation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "approx = defaultdict(int)\n", "mafs = []\n", "for fname in os.listdir('.'):\n", " if not fname.startswith('tsi-') or not fname.endswith('.frq'):\n", " continue\n", " f = open(fname)\n", " f.readline()\n", " for cnt, l in enumerate(f):\n", " toks = [tok for tok in l.rstrip().split(' ') if tok != '']\n", " maf = float(toks[-2])\n", " mafs.append(maf)\n", " approx[maf] += 1\n", " f.close()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1222126 \n", "417\n" ] } ], "source": [ "print(len(mafs), type(mafs))\n", "print(len(approx))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.2216\n", "CPU times: user 60 ms, sys: 4 ms, total: 64 ms\n", "Wall time: 65.9 ms\n" ] } ], "source": [ "%time print(np.median(mafs))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAKCCAYAAACqK24GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3X+MbOd5H/bvY99IoRuWNK1CpkjqR+pLIHSVSGYsJk0a\nXUcxTRUFybSuyBKxWIRwBTGVAqNoIqaoSNVAEjmxBRmpBAmmJEqtWBKRY8sQQ5OSvagClLq0K8m0\nKUIUYCa8VyaVUqKUQElDVk//2HPD4d6ZvTuzMztzdj8f4PCeeeac9zxnZnbn5bPveU91dwAAAAAY\nr+9bdwIAAAAA7I8CDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDILVTgqao/XlVf\nqKovVdWjVfX3hvgdVXWqqr44LG+e2Oe2qnq8qh6rqqsn4ldW1SPDc++fiL+0qu4Z4g9V1av2c6IA\nAHulrwMAjM1CBZ7u/rdJfqK7X5fkTyf5iar6i0k6yS919+uH5Z8mSVVdkeSGJFckuSbJB6qqhuY+\nmOSW7j6e5HhVXTPEb0nyzBB/X5L3LnaKAADz0dcBAMZm4Uu0uvu7w+pLknx/km8Nj2vK5tclubu7\nn+vuJ5J8LclVVXVxkvO7++Sw3ceTXD+sX5vkrmH9U0netGiuAADz0tcBAMZk4QJPVX1fVX0pydNJ\nfru7/2B46h1V9eWqurOqLhxir0hyamL3U0kumRI/PcQz/PtkknT380m+XVUXLZovAMA89HUAgDHZ\nzwie7w3Dli9N8peq6kS2hyC/JsnrkvxRkl9cRpIAAAdNXwcAGJNj+22gu79dVZ9J8me7e+tMvKp+\nJclvDA9PJ7lsYrdLs/3XrNPD+s74mX1emeTrVXUsyQXd/c2dx6+q3u85AADj0t3TLpNa1bHW1tfR\nzwGAo2mRvs5CBZ6qelmS57v72ao6L8lPJnlPVf1wdz81bPZXkzwyrH86ySer6peyPRz5eJKT3d1V\n9Z2quirJySQ/k+SXJ/a5OclDSX46yedm5XOQnTz2pqru6O471p0HL+Z92Vzem83kfdlMB1H02KS+\njn7OZvL7YTN5XzaT92VzeW8206J9nUVH8Fyc5K6q+r5sX+b1ie7+XFV9vKpel+07TPxhkrclSXc/\nWlX3Jnk0yfNJbu3uMwnfmuRjSc5Lcl933z/E70zyiap6PMkzSW5cMFcAgHnp6wAAo7JQgae7H0ny\nY1Pib91ln7+b5O9Oif9uktdOif+/Sd6ySH4AAPuhrwMAjM3CkyzDOWytOwGm2lp3Asy0te4EmGpr\n3QkAG2tr3Qkw1da6E2CqrXUnwExb606A5akXRg+PU1W1a9MB4Og4St/9R+lcAYBti37/G8EDAAAA\nMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8AD\nAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAj\np8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHIKPAAA\nAAAjp8ADAAAAMHIKPAAAAAAjp8ADAAAAMHLH1p3AslXlw0kuT/LdJDd159k1pwQAAACwUoeuwJPt\n4s4bh/UPJblhjbkAACxdVfW0eHfXQecCAGyGw1jg+e7w78NJ3rbORAAAVmdnjUdtBwCOssM4B89N\nSe5NcrXLswAAAICjYKECT1X98ar6QlV9qaoeraq/N8QvqqoHq+qrVfVAVV04sc9tVfV4VT1WVVdP\nxK+sqkeG594/EX9pVd0zxB+qqlftJbfuPNudGxR3AIBFbXJfBwBgmoUKPN39b5P8RHe/LsmfTvIT\nVfUXk7wryYPdfXmSzw2PU1VXZHsunCuSXJPkA1V1ZhzxB5Pc0t3HkxyvqmuG+C1Jnhni70vy3kVy\nBQCYl74OADA2C1+i1d1n5rp5SZLvT/KtJNcmuWuI35Xk+mH9uiR3d/dz3f1Ekq8luaqqLk5yfnef\nHLb7+MQ+k219KsmbFs0VAGBe+joAwJgsXOCpqu+rqi8leTrJb3f3HyR5eXc/PWzydJKXD+uvSHJq\nYvdTSS6ZEj89xDP8+2SSdPfzSb5dVRctmi8AwDz0dQCAMVn4Llrd/b0kr6uqC5L8ZlX9xI7ne9Yt\nPAEANp2+DgAwJvu+TXp3f7uqPpPkyiRPV9UPd/dTw5DkbwybnU5y2cRul2b7r1mnh/Wd8TP7vDLJ\n16vqWJILuvub03KoqjsmHm5199b+zgoA2BRVdSLJiXUdf919Hf0cADjcltXXqe75//BUVS9L8nx3\nP1tV5yX5zSTvSfJT2Z4s8L1V9a4kF3b3u4aJBz+Z5A3ZHo782SQ/Mvzl6wtJ3pnkZJLPJPnl7r6/\nqm5N8trufntV3Zjk+u6+cUou3d21Mw4AHE4H8d2/KX2dWee6PXJoZx+uok8EAOO3aF9n0RE8Fye5\nq6q+L9vz+Hyiuz9XVV9Mcm9V3ZLkiSRvSZLufrSq7k3yaJLnk9zaL1SWbk3ysSTnJbmvu+8f4ncm\n+URVPZ7kmSRnFXcAAFZEXwcAGJWFRvBsEiN4AOBoOUrf/UbwAMDRs2hfZ+G7aAEAAACwGRR4AAAA\nAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4\nAAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg\n5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcA\nAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZO\ngQcAAABg5BYq8FTVZVX121X1B1X1+1X1ziF+R1WdqqovDsubJ/a5raoer6rHqurqifiVVfXI8Nz7\nJ+Ivrap7hvhDVfWq/ZwoAMBe6esAAGOz6Aie55L8XHf/aJI/l+RvVNWfStJJfqm7Xz8s/zRJquqK\nJDckuSLJNUk+UFU1tPXBJLd09/Ekx6vqmiF+S5Jnhvj7krx3wVwBAOalrwMAjMpCBZ7ufqq7vzSs\n/+skX0lyyfB0TdnluiR3d/dz3f1Ekq8luaqqLk5yfnefHLb7eJLrh/Vrk9w1rH8qyZsWyRUAYF76\nOgDA2Ox7Dp6qenWS1yd5aAi9o6q+XFV3VtWFQ+wVSU5N7HYq252knfHTeaHzdEmSJ5Oku59P8u2q\numi/+QIAzENfBwAYg30VeKrqTyT5x0n+5vDXrQ8meU2S1yX5oyS/uO8MAQDWRF8HABiLY4vuWFV/\nLNvDif+37v61JOnub0w8/ytJfmN4eDrJZRO7X5rtv2adHtZ3xs/s88okX6+qY0ku6O5vzsjljomH\nW929tdhZAQCbpqpOJDmxhuNuRF9HPwcADrdl9XWquxc5eGX7mvFnuvvnJuIXd/cfDes/l+THu/um\nYeLBTyZ5Q7aHI382yY90d1fVF5K8M8nJJJ9J8svdfX9V3Zrktd399qq6Mcn13X3jlFy6u6ddCw8A\nHEIH8d2/KX2dWedaVb093/OLotEnAoDxW7Svs+gInr+Q5K8l+b2q+uIQ+ztJ/puqel22exx/mORt\nSdLdj1bVvUkeTfJ8klv7hcrSrUk+luS8JPd19/1D/M4kn6iqx5M8k+Ss4g4AwIro6wAAo7LQCJ5N\nYgQPABwtR+m73wgeADh6Fu3r7PsuWgAAAACslwIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACM\nnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAA\nAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp\n8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAA\nwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMgp8AAAAACM3EIFnqq6rKp+u6r+oKp+v6re\nOcQvqqoHq+qrVfVAVV04sc9tVfV4VT1WVVdPxK+sqkeG594/EX9pVd0zxB+qqlft50QBAPZKXwcA\nGJtFR/A8l+TnuvtHk/y5JH+jqv5UknclebC7L0/yueFxquqKJDckuSLJNUk+UFU1tPXBJLd09/Ek\nx6vqmiF+S5Jnhvj7krx3wVwBAOalrwMAjMpCBZ7ufqq7vzSs/+skX0lySZJrk9w1bHZXkuuH9euS\n3N3dz3X3E0m+luSqqro4yfndfXLY7uMT+0y29akkb1okVwCAeenrAABjs+85eKrq1Ulen+QLSV7e\n3U8PTz2d5OXD+iuSnJrY7VS2O0k746eHeIZ/n0yS7n4+yber6qL95gsAMA99HQBgDPZV4KmqP5Ht\nvzj9ze7+V5PPdXcn6f20DwCwTvo6AMBYHFt0x6r6Y9nu8Hyiu39tCD9dVT/c3U8NQ5K/McRPJ7ls\nYvdLs/3XrNPD+s74mX1emeTrVXUsyQXd/c0Zudwx8XCru7cWPS8AYLNU1YkkJ9Zw3I3o6+jnAMDh\ntqy+Tm3/8Wnug1e2rxl/prt/biL+C0PsvVX1riQXdve7hokHP5nkDdkejvzZJD/S3V1VX0jyziQn\nk3wmyS939/1VdWuS13b326vqxiTXd/eNU3Lp7q6dcQDgcDqI7/5N6evMOteq6rMHD1X0iQBg/Bbt\n6yxa4PmLSf7PJL+XF3oXt2W743Jvtv8a9USSt3T3s8M+fyfJX0/yfLaHOf/mEL8yyceSnJfkvu4+\ncxvSlyb5RLaveX8myY3DpIU7c1HgAYAj5IAKPBvR11HgAYCj50ALPJtEgQcAjpaj9N2vwAMAR8+i\nfZ1930ULAAAAgPVS4AEAAAAYOQUeAAAAgJFT4AEAAAAYOQUeAAAAgJFT4AEAAAAYOQUeAAAAgJFT\n4AEAAAAYOQUeAAAAgJFT4AEAAAAYOQUeAAAAgJFT4AEAAAAYOQUeAAAAgJFT4AEAAAAYOQUeAAAA\ngJFT4AEAAAAYuWPrTuCgVOXDSS5P8t0kN3Xn2TWnBAAAALAUR6bAk+3izhuH9Q8luWGNuQAALF1V\n9bR4d9dB5wIAHKyjVOD57vDvw0nets5EAABWZ2eNR20HAI6CozQHz01J7k1ytcuzAAAAgMOkuqeO\n5B2NqmrDjgHg6DhK3/2zznX7UqxZI3XOjh+V1wsADoNF+zpHaQQPAAAAwKGkwAMAAAAwcgo8AAAA\nACOnwAMAAAAwcgo8AAAAACOnwAMAAAAwcgo8AAAAACOnwAMAAAAwcgo8AAAAACOnwAMAAAAwcgo8\nAAAAACOnwAMAAAAwcgo8AAAAACO3cIGnqj5SVU9X1SMTsTuq6lRVfXFY3jzx3G1V9XhVPVZVV0/E\nr6yqR4bn3j8Rf2lV3TPEH6qqVy2aKwDAPPRzAICx2c8Ino8muWZHrJP8Une/flj+aZJU1RVJbkhy\nxbDPB6qqhn0+mOSW7j6e5HhVnWnzliTPDPH3JXnvPnIFAJiHfg4AMCoLF3i6+/NJvjXlqZoSuy7J\n3d39XHc/keRrSa6qqouTnN/dJ4ftPp7k+mH92iR3DeufSvKmRXMFAJiHfg4AMDarmIPnHVX15aq6\ns6ouHGKvSHJqYptTSS6ZEj89xDP8+2SSdPfzSb5dVRetIF8AgL3SzwEANtKxJbf3wST/y7D+80l+\nMdtDkFeqqu6YeLjV3VurPiYAcDCq6kSSE2tOI9HPAQBWYFl9naUWeLr7G2fWq+pXkvzG8PB0kssm\nNr0023/ROj2s74yf2eeVSb5eVceSXNDd35xx3DsWzbkqH05yeZLvJrmpO88u2hYAsHxDQWPrzOOq\nun1NeYyunwMAbL5l9XWWeonWcK35GX81yZk7T3w6yY1V9ZKqek2S40lOdvdTSb5TVVcNkxH+TJJf\nn9jn5mH9p5N8bpm5Trg8yRuTvDnJh1Z0DABg5EbazzmTe09bVnlMAOBgLTyCp6ruznZh5GVV9WSS\n25OcqKrXZfsuE3+Y5G1J0t2PVtW9SR5N8nySW7v7TKfi1iQfS3Jekvu6+/4hfmeST1TV40meSXLj\normew3eHfx8+ky8AcLQdon7OhJ31nGnzRQMAY1Uv9D/Gqaq6uxfuoVTlwmyP3Hmby7MAYPPt97t/\nTGad6/bom1kFm73Hj8rrCABjsmhf58gXeA6beecUMgcRAGNzlL77FXgA4OhZtK+z7LtosX5n5hRK\ntkcm3bDk7QEA1mrW/EGzOsPzbg8AY6TAswYrHmUzdU6hXdowBxEAMELzzilkDiIADrel3kXrMKnK\nh6uyVZX7hnl6lmnqnbt2OeY8d/q6Kcm9Sa7eUQia1cas7ada8esCABwS89y5a1l3+XK3MACOMiN4\nZtv3pUsLjJqZdcw9j7IZjjEt16ltzNp+l9zPytG8PwAwTqu/dOnsUTOzCy7LGmGzv3ZczgXAWBnB\nM9syLl2ad9TMrGNO3X7O0TRzjdTZJfdpOa5yRBIAsDI9ZZnPYqNm9nfMZZmd99n5zXueRhMBcNAU\neGbbc1Fll0LGzFEz3blhSqFl6jF32X7PhZJd2phlnmLTuUYk7aVItBYuOQOAZVhNwWb1RZJ5C1zz\nFn42o5AFwNHgNulzHy9beeESpXu7c8O02LDthdkubLxtFZchVeW+bBdPHs7eR+bste095z5r21n5\nzfu6rPKSrlnvHQCby23Sl3eb9OlFh0Ve2lXkskgb+z3m6nNZx2fXZWcA47JoX0eBZ+7jnV20WGWh\n5Ry5rLSAtF/LKuTsUkA7a/sF5gNay3sHwOIUeA6iwLMJxZajU+BZdQFm1uflqPwcAYyNAs+BHe/s\nosWmF1rGYpdCzqyRQGdtP08xaIh77wBGRoFHgWd1xzyoXKbZ/yTTs34uFHgAxmXRvo67aM1p2l2n\ndrlzFfPZbd6faUWYadvPdYcy7x0AcLDmvSxu78UjkzgDHG0mWWaTzDvJ9LTt571D2VQmXwYAxsnE\nzgBHlUu0OBIWmA9oK3Nc6jWjjZVNDg1wlB2l736XaB30MQ8il016zac7Kj9fAJvKJVqwiwUuxdrz\npV67FHKmXham8AMAbIbpxaZpl3rtPr/PlJYViQAOnEu0YLp5LvU6U8h5c7YLObttu9v2U826XMxl\nZADAapx9mVdV9bRlse0BWAUFHphiznl/dpsces/zAe1SsJlVEJqrUAQAsD/zzu+zt+1nFYMUiQDm\no8ADc5hR+FnG5NDJ/COB9jxxtNE+AMAmmGcU0O5xAHYyyTJsiKrcl+3izsOZKP7MmiB6WnzW/D7L\nmDQaYFMcpe9+kywf9DEPIpdNfs0PIpf9H9N8QMBht2hfR4EHNsS8d/qa0cZWphdyZhWPpm4PsMmO\n0ne/As9BH/Mgctnk1/wgclnWMWdxtzBg/Bbt67hECzbELpd0zWPV8wGdZV2Xf7nsDACOstXMBwQw\nZkbwwCEy7yigXS7/2soeL+la1yggo4/g6DpK3/1G8Bz0MQ8il01+zQ8il006/+lm/8ztbdtlcckZ\nHF2L9nWOrSIZYD2GIs2eCx27bH+uW7wn24WhG2Ztu8t8QMua92fPk0wDAJxterFp9p265r1cbBnW\ncUxgrFyiBUwzzyVd894VbFm3d5913D1zmRdsNj+jwHrs/XIut3gHNokRPMBZdhnZc1N2XNK1wCig\npYz4mXe00gzTRiQBm8PPKDACy5ogejVc6gWba9mFXyN4gD2bcyLoWSNsDnzEzy6jAFzmBZvNzyhw\nqKxvtI9JpmFzLe/nU4EHWIlZxaBdikRzjfiZ06wi0VyXeblcBA7cvi/FBNgsO/9HrmfE5+eyMMBd\ntICNsMsdvfZ8Z7BdLvO6L9vFnYezh/9R3KWdrezzzl1LnGQajqyj9N3vLloHfcyDyGWTX/ODyMX5\n7/eYs37/zft7YVPuFgZH2W4/t+6iBYzWrDl15pxrZ9Z8HWfNHbRgO2eNJlrgbmFT21b4AQD2Ylmj\nclZ5tzCFIlgPBR7gMJl6OdcCEzLPuixsWqFoVjFoz0Wic2y/Z/MWidZRVFLIAoBlWMYEzntvY7Gi\n0tntK/zAapmDBzhMljVfx9R2ZswfNO/cQfPcgn6mGfMBTZ1raJe5g+aawHqeOYiWdUwAYFPMmido\nnvmDZs1BtHfzzjW0jrmJzIfEuijwAIfGnHf5WlY7c90tbJe2p24/Z6HkXKODdhZVZt6yfgnFmbmO\nCQAwafciydlFonm3n/+Y+88RVk2BB2AfFrhb2FztZL5Cybyjg+a9Zf3UOYjmvAX9vIWsfXP3s83n\nPQI42pZXJNn79us4Jqyau2gBbLBZdwCb8+5ie9523mPOurPYAsec2s4u2+95Lp9dclzZfEALTL59\npM3//h+d73530TroYx5ELpv8mh9ELs7/4I+511wO62u+11zO1cYs+29nGd9p5jfabLtfojf9s7LI\ne7fwCJ6q+khVPV1Vj0zELqqqB6vqq1X1QFVdOPHcbVX1eFU9VlVXT8SvrKpHhufePxF/aVXdM8Qf\nqqpXLZorwIjNe6nXWRa4dG3fcxAtcMx5L92a53KxeS9dm2raKJMFLmeba56kJc17NNf2yzjmAm1v\n3KV7+jkAvNgyRursbGN5l4vNynGV8wHNatscRLPs/f1f/BDdCy1J/rMkr0/yyETsF5L8rWH9byf5\n+8P6FUm+lOSPJXl1kq/lhdFDJ5O8YVi/L8k1w/qtST4wrN+Q5P+YkUcveg4Wi8Vi2d+S9IVJ35P0\nhatoJ+kPJ72V9H2Tzw2PO+mT5zr2Lm1PbWOXY24N23fS98yKnaPtWfFZ7UyNzzjPPW+7yDGnvS4L\ntDErPtfn6CC++ze9n5OkJ17LYTnTW5wnvjM2K75I2/vNZR3HPIhcNvk1d/6H85hH/Wdu889/nt/z\nu303zpPL7DbmXeZ9Xfa27H6e+99+VW3M91rt/llcpP+y8Aie7v58km/tCF+b5K5h/a4k1w/r1yW5\nu7uf6+4nst3xuaqqLk5yfnefHLb7+MQ+k219KsmbFs0VgNXo1U9sPWskzFmjjGaNDtml7X3PQTQj\ntlvb886TtO95jxYYNTPPiKd521jWiK+V088B4KCsY7TL/HMQzYrPY1rdY5G2l7H92bHFJs1e5es1\nv2VPsvzy7n56WH86ycuH9VckOTWx3akkl0yJnx7iGf59Mkm6+/kk366qi5acLwCbbZ6CwFyXXO1S\nVJinaDPXJXQLFJumxfdc9Fpw+3mKUPO2MSs+Fvo5AKzA3osBy7v8aT0FiL1a7h3N5rHZr8u5HFtV\nw919YNfZVdUdEw+3unvrII4LwMrdlL1P1ryseVymHnNYf9Hkv9Nii5jVzoz4zKLXjFzm2n6Xds56\nXeZtY9HXq6pOJDkx736rpJ8DwPrs/Po5rPMoz3ueY35dtoZlf5Zd4Hm6qn64u58ahiV/Y4ifTnLZ\nxHaXZvsvWqeH9Z3xM/u8MsnXq+pYkgu6+5vTDtrddyzvFADYFHMWBOYpBi3rmOsw73mO/nUZChpb\nZx5X1e3ryCP6OQDASpzIC3/Les/CrSz7Eq1PJ7l5WL85ya9NxG+sqpdU1WuSHE9ysrufSvKdqrqq\nqirJzyT59Slt/XSSzy05VwAOkU2cx2UV5j3Po/K6HBD9HABgYy08gqeq7s72Nf0vq6onk7w7yd9P\ncm9V3ZLkiSRvSZLufrSq7k3yaJLnk9zaw3TU2b6LxMeSnJfkvu6+f4jfmeQTVfV4kmeS3LhorgAA\n89DPAQDGpl7of4xTVXV3j+niOgBgH47Sd/+sc92e/2fWXAPzxKf1A6fFF2l7v7ms45gHkcsmv+YH\nkYvzP/hj7jWXw/qa7zUX53/wx9ykXDbv/Bfp6yz7Ei0AAAAADpgCDwAAAMDIKfAAAAAAjJwCDwAA\nAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwC\nDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAA\njJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAA\nAAAAjJwCDwAwGlX58LpzAADYRAo8AMCYXL7uBAAANpECDwAwJt9ddwIAAJtIgQcAGJOb1p0AAMAm\nUuABAEajO8+uOwcAgE2kwAMAAAAwcisp8FTVE1X1e1X1xao6OcQuqqoHq+qrVfVAVV04sf1tVfV4\nVT1WVVdPxK+sqkeG596/ilwBAOahnwMAbKJVjeDpJCe6+/Xd/YYh9q4kD3b35Uk+NzxOVV2R5IYk\nVyS5JskHqqqGfT6Y5JbuPp7keFVds6J8AQD2Sj8HANg4q7xEq3Y8vjbJXcP6XUmuH9avS3J3dz/X\n3U8k+VqSq6rq4iTnd/fJYbuPT+wDALBO+jkAwEZZ5Qiez1bV71TVzw6xl3f308P600lePqy/Ismp\niX1PJblkSvz0EAcAWCf9HABg4xxbUbt/obv/qKr+oyQPVtVjk092d1dVL+tgVXXHxMOt7t5aVtsA\nwHpV1YkkJ9acxiT9HABgibaGZX9WUuDp7j8a/v2XVfVPkrwhydNV9cPd/dQwLPkbw+ank1w2sful\n2f6L1ulhfTJ+esbx7ljuGQAAm2IoaGydeVxVt68tmejnAADLdiIv/C3rPQu3svRLtKrqB6rq/GH9\nP0hydZJHknw6yc3DZjcn+bVh/dNJbqyql1TVa5IcT3Kyu59K8p2qumqYjPBnJvYBADhw+jkAwKZa\nxQielyf5J8MNIo4l+d+7+4Gq+p0k91bVLUmeSPKWJOnuR6vq3iSPJnk+ya3dfWZY861JPpbkvCT3\ndff9K8iu9HoLAAAgAElEQVQXAGCv9HMAgI1UL/Qxxqmqurt33skCADikjtJ3/6xz3Z7jZ2cf7sxm\n88Sn9QOnxRdpe7+5rOOYB5HLJr/mB5GL8z/4Y+41l8P6mu81F+d/8MfcpFw27/wX6eus8jbpAAAA\nABwABR4AAACAkVvVbdLXpiofTnJ5ku8muak7z645JQAAAICVOowjeC5P8sYkb07yoTXnAgCwIrVj\nAQCOskM3gifbI3eS5OEkb1tnIgAAqzNrYkcA4Cg6jCN4bkpyb5KrXZ4FAAAAHAWHbgTPUNS5Yd15\nAAAAAByUwziCBwAAAOBIUeABAAAAGDkFHgAAAICRU+ABAAAAGDkFHgAAAICRU+ABAAAAGDkFHgAA\nAICRU+ABAAAAGDkFHgAAAICRU+ABAAAAGDkFHgAAAICRU+ABAAAAGDkFHgAAAICRU+ABAAAAGLlj\n607goFTlw0kuT/LdJDd159k1pwQAAACwFEdpBM/lSd6Y5M1JPrTmXAAAVqR2LADAUXBkRvBke+RO\nkjyc5G1ngps2smdaPrNyXFZ8nlyW2T4AsAq947EiDwAcBUdpBM9NSe5NcvWOwsPUkT1V+XBVtqpy\nX1UuPKj4jHxmjT5aSnzOXOaKr+t1PID36SzLaGOZ7cxjlW0vM59NyxNgXDZpZM+0XGblt6z4XvJY\nNBcA2ADdPepl+xT2s3/fl3QnfTLpCyfiW0O8k77nAONn5bNLjsuK7zmXeeNrfB1XFk/6w0P8vn1+\nZlbdzlnxBY6557bP8XM2VztLej9WGl/iue75fVrg99tGnetYc9+092OVn7FVvqfLff/Si+YwtmXW\nuSY587twYkkvFl/Hsqzcl3GuB/06ztpnWfFl5TLr3PfbzqzXbZFc5o2fq+2dOa7ic7jb+e8lx1Xk\ncpA/h+t4zZ3/5h5zk3LZvPNfqN+w7o7Lqjo+e9+/L0z6nuzoUGb1RZVZ8bPy2SXHZcX3nMu88TW+\njiuLZ3kFsVW3c1Z8gWPuue0hvqyi1RiKiMs613nep5UVylZ9riPPfdPej1V+xpb1+i7lNZi1JOlz\nbXNYllnnmoyjc7rZ8XUdc9prtsz4MnPZuSya4842dm6/n1zmjc86/qwcp8Vnva57yWW3819mjnuJ\nTzuXRc9z1cfcy8/yXuLzHnPZue/nPJf9uhxkLus4/03MZWds91wW6jesu+Oyqo7P/ttdeVFlanw9\nr+HqclnX67jKeJZXEFt1O/OMBltWoXNr4pfTfopWYygiLutc53mfZh1zZbksq/2R575p78cqP2Mr\nKzAv0v6sZVXf/Zu4zDrXsXRONzu+rmNOe82WGd+kXHbGz7W+Sbmca33a/zyf65w2aX2vr9Fez3PZ\n7+9uP7PLeN9X0fZecu8VxBfJ/aB/903L4yjmsjO2ey4L9RvW3XFZRccnSxpKbrFMW7Kkgtiq25kW\nn/eY87Q9xJdStFrV+S85vsqRb8squG3MuY489017P1b5GVtZgXmZvwsW7fSMcZl1rmPpnG52fF3H\nnPaaLTO+SbnsjJ9rfZNyWdb6Oo65rNd6Hcc8qPd9k382541vUi7niq/jmJuUy87Y7rks0m+ooZMw\nWlXV3V0vjmUr25P9Jsm93bnhwBODI2iY+PhDSd7Wh/wuaus411nHXHUuy2h/zLnP27afg4N4X8/+\n7j+sZp1rVe3Sgdv5VC0xPu2w0+LLPOaq4us6Zufs12yZ8U3KZWf8XOublMuy1tdxzHnzOqh8D/L9\nPdc5rSOX7IgtK75JuZwrvgm/+9eZy3zfoYv0dQ5rgee+bN/J6eGcfdcsAGDEFHgyFHg2v3O62fF1\nHXNT/idyE4sqm5TLstbXccx58zqofNdRVJl1TuvIJTtiy4pvUi7nim/C7/515jLfd+gifZ3Depv0\nWbdEBwAAADh0jq07gVUYijouywIAAACOhI0fwVNV11TVY1X1eFX97b3tkw9XZasq9w3zAEyNsTpV\ndWLdOXA278vm8t5sJu8Lq7ZIP4dNsbXuBJhqa90JMNXWuhNgpq11J8ASbfQInqr6/iT/KMlfSXI6\nycNV9enu/so5dr08L0yy/KFsj+aZFktVPjw8990kN525pEt8v/G3X1KVH19nPut/DTYxfsHvZfgt\nvhn5LB6fZtNynDN+oio3bVA+o4iv/nN28y9W5V9twrmK7/7zP0b76OecaWHOOMu1leTEmnPgbFvr\nToCpttadADNtxe+yQ2Tdt/88x61B/3yS+ycevyvJu851+9Ck/0W2bzH2bNKvGmJnbtn6naQfyHAr\n1mzfTr2H5Z6JNsT3Fb997fms/zXYxPjP/v5m5bN4fPrvjM3Kcb7Pa+7YpHzGEl/9Md/xxKacq/hk\nPL3b74KxLFmwn3MmPvHaDEt697hlM5ZZ78my4vPs02uI9y7xnc8vEr99D9vvJZdlxac9v4z1na/z\nvPG9fm6WdczbZxzjoD6D6/7cL/v13O/PyWTs9hnxWdvvJz7r+2kV8U3KZef7cK74dhuL9C02egRP\nkkuSPDnx+FSSq/aw379IclmSC5L8QrZH69yU5OtJzk/yk0k+kuS/zPZfBZPk/0vyl6vyqu788w2O\nP5fkJ6vyQJK3bG78e99b8XH/VZIfrMqFvf0X3cn4f1qVf5btv/omyb9McklV/kW2PxsHHb8v25+/\n766w/W8kee0Q/39mx8/7D+fbfi3x/zvJv6nKVpI/meRPTMTflt3t/Hx8LMmrk1yx+fF3nEjyAwu2\ns6zP1MNJ3jYxamLZ7a8iz8nPzZnYf5HkZcuJX/jy1ba/v3iSVOUrSX5kiH/5sMcPmUX7OYN5R/D0\njO02IT6Z80HF99vGHcOy6DHXHZ/2OVl1nOU5rJ+RVcY3KZdlxuEFG32b9Kr6r5Jc090/Ozz+a0mu\n6u53TGyzuScAAKxEH4LbpOvnAACzLNLX2fQRPKezPRLnjMuy/detf+8wdPAAgCNJPwcAWJpNv4vW\n7yQ5XlWvrqqXZPtSq0+vOScAgGXQzwEAlmajR/B09/NV9d8n+c0k35/kzt7znSUAADaXfg4AsEwb\nPQcPAAAAAOe26ZdoAQAAAHAOCjwAAAAAI6fAAwAAADByCjwAAAAAI6fAAwAAADByCjwAAAAAI6fA\nAwAAADByCjwAAAAAI6fAAwAAADByCjwAAAAAI6fAAwAAADByuxZ4quojVfV0VT0yEXtDVZ2sqi9W\n1cNV9eMTz91WVY9X1WNVdfVE/MqqemR47v0T8ZdW1T1D/KGqetXEczdX1VeH5a3LO2UAgG0z+jr/\noKq+UlVfrqpfraoLJp7T1wEANtK5RvB8NMk1O2K/kOR/7u7XJ3n38DhVdUWSG5JcMezzgaqqYZ8P\nJrmlu48nOV5VZ9q8JckzQ/x9Sd47tHXR0PYbhuX2qrpw4bMEAJhuWl/ngSQ/2t1/JslXk9yW6OsA\nAJtt1wJPd38+ybd2hP8oyZm/ZF2Y5PSwfl2Su7v7ue5+IsnXklxVVRcnOb+7Tw7bfTzJ9cP6tUnu\nGtY/leRNw/pPJXmgu5/t7meTPJizO18AAPsyra/T3Q929/eGh19Icumwrq8DAGysYwvs864k/6yq\n/mG2C0R/foi/IslDE9udSnJJkueG9TNOD/EM/z6ZJN39fFV9u6p+aGjr1JS2AAAO0l9Pcvewrq8D\nAGysRQo8dyZ5Z3f/k6r6r5N8JMlPLjetvauqXtexAYD16O4691b7U1X/U5J/192fXPWxdslBPwcA\njqBF+jqLFHje0N1/ZVj/x0l+ZVg/neSyie0uzfZfo07nhaHNk/Ez+7wyyder6liSC7r7mao6neTE\nxD6XJfmtWQkdRCeP+VTVHd19x7rz4MW8L5vLe7OZvC+b6SCKHlX13yb5z/PCJVXJmvo6+jmbye+H\nzeR92Uzel83lvdlMi/Z1FrlN+teq6o3D+l/O9uSDSfLpJDdW1Uuq6jVJjic52d1PJflOVV01TET4\nM0l+fWKfm4f1n07yuWH9gSRXV9WFVfWD2R4h9JsL5AoAMJdhguT/Mcl13f1vJ57S1wEANtauI3iq\n6u4kb0zysqp6Mtt3e/jvkvyvVfXSJP9meJzufrSq7k3yaJLnk9za3WeqTrcm+ViS85Lc1933D/E7\nk3yiqh5P8kySG4e2vllVP5/k4WG79wwTEAIALM2Uvs7t2b5r1kuSPDjcJOv/6u5b9XUAgE1WL/RL\nxqmq2tDlzVNVJ7p7a9158GLel83lvdlM3pfNdJS++4/SuY6N3w+byfuymbwvm8t7s5kW/f5X4AEA\nRuUoffcfpXMFALYt+v2/yBw8AAAAAGwQBR4AAACAkVPgAQAAABg5BR4AAACAkVPgAQAAABg5BR4A\nAACAkVPgAQAAABg5BR4AAACAkVPgAQAAABg5BR4AAACAkVPgAQAAABg5BR4AAACAkVPgAQAAABg5\nBR4AAACAkVPgAQAAABi5Y+tOYBWqqtedw7J1d607BwAAAGAzHcoCz7bDVONR2wEAAABmc4kWAAAA\nwMgp8AAAAACMnAIPAAAAwMgp8AAAAACMnAIPAAAAwMjtWuCpqo9U1dNV9ciO+Duq6itV9ftV9d6J\n+G1V9XhVPVZVV0/Er6yqR4bn3j8Rf2lV3TPEH6qqV008d3NVfXVY3rqc0wUAAAA4fM41guejSa6Z\nDFTVTyS5Nsmf7u7/JMk/HOJXJLkhyRXDPh+oqjP39/5gklu6+3iS41V1ps1bkjwzxN+X5L1DWxcl\neXeSNwzL7VV14X5OFAAAAOCw2rXA092fT/KtHeG3J/l73f3csM2/HOLXJbm7u5/r7ieSfC3JVVV1\ncZLzu/vksN3Hk1w/rF+b5K5h/VNJ3jSs/1SSB7r72e5+NsmD2VFoAgAAAGDbInPwHE/yl4ZLqraq\n6s8O8VckOTWx3akkl0yJnx7iGf59Mkm6+/kk366qH9qlLQAAAAB2OLbgPj/Y3X+uqn48yb1J/uRy\n0wIAAABgrxYp8JxK8qtJ0t0PV9X3qupl2R6Zc9nEdpcO254e1nfGMzz3yiRfr6pjSS7o7meq6nSS\nExP7XJbkt2YlVFV3TDzcmv+UAIBNVVUn8uJ+AQAAO1R3775B1auT/EZ3v3Z4/LYkr+ju26vq8iSf\n7e5XDpMsfzLbkyJfkuSzSX6ku7uqvpDknUlOJvlMkl/u7vur6tYkr+3ut1fVjUmu7+4bh0mWfyfJ\njyWpJL+b5MeG+Xh25tfdXTtjye7nNS6VnecIAEfVtO/+w+oonSsAsG3R7/9dR/BU1d1J3pjkh6rq\nyWzf2eojST4y3Dr93yV5a5J096NVdW+SR5M8n+TWfqF6dGuSjyU5L8l93X3/EL8zySeq6vEkzyS5\ncWjrm1X180keHrZ7z7TiDgAAAAB7GMGz6YzgAYCj5SiNajlK5woAbFv0+3+Ru2gBAAAAsEEWmWQZ\nAIA12h6tfDCMIAKAcVDgAQAYpYOo8ajtAMBYuEQLAAAAYOQUeAAAAABGToEHAAAAYOQUeAAAAABG\nToEHAAAAYOQUeACAI6uqPlJVT1fVIxOxi6rqwar6alU9UFUXTjx3W1U9XlWPVdXVE/Erq+qR4bn3\nT8RfWlX3DPGHqupVE8/dPBzjq1X11oM4XwDg8FLgAQCOso8muWZH7F1JHuzuy5N8bnicqroiyQ1J\nrhj2+UBVnbmP+AeT3NLdx5Mcr6ozbd6S5Jkh/r4k7x3auijJu5O8YVhunywkAQDMS4EHADiyuvvz\nSb61I3xtkruG9buSXD+sX5fk7u5+rrufSPK1JFdV1cVJzu/uk8N2H5/YZ7KtTyV507D+U0ke6O5n\nu/vZJA/m7EITAMCeKfAAALzYy7v76WH96SQvH9ZfkeTUxHanklwyJX56iGf498kk6e7nk3y7qn5o\nl7YAABZybN0JAABsqu7uqup15lBVd0w83OrurTWlAgCsQFWdSHJiv+0o8AAAvNjTVfXD3f3UcPnV\nN4b46SSXTWx3abZH3pwe1nfGz+zzyiRfr6pjSS7o7meq6nRe3JG7LMlvTUumu+/Y3+kAAJts+OPN\n1pnHVXX7Iu24RAsA4MU+neTmYf3mJL82Eb+xql5SVa9JcjzJye5+Ksl3quqqYdLln0ny61Pa+uls\nT9qcJA8kubqqLqyqH0zyk0l+c5UnBQAcbkbwAABHVlXdneSNSV5WVU9m+85Wfz/JvVV1S5Inkrwl\nSbr70aq6N8mjSZ5Pcmt3n7l869YkH0tyXpL7uvv+IX5nkk9U1eNJnkly49DWN6vq55M8PGz3nmGy\nZQCAhdQL/ZJxqqru7toZS8Z9Xi9W2XmOAHBUTfvuP6xmnevB9XX0QQDgoC3a13GJFgAAAMDIKfAA\nAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDIKfAAAAAAjJwCDwAAAMDI7VrgqaqPVNXTVfXIlOf+\nh6r6XlVdNBG7raoer6rHqurqifiVVfXI8Nz7J+Ivrap7hvhDVfWqiedurqqvDstb93+qY1Xb/63q\nw7Cs+cUEAACAQ+lcI3g+muSancGquizJTyb55xOxK5LckOSKYZ8PVFUNT38wyS3dfTzJ8ao60+Yt\nSZ4Z4u9L8t6hrYuSvDvJG4bl9qq6cKEzPBT6kCwAAADAKuxa4Onuzyf51pSnfinJ39oRuy7J3d39\nXHc/keRrSa6qqouTnN/dJ4ftPp7k+mH92iR3DeufSvKmYf2nkjzQ3c9297NJHsyUQhMAAAAAC8zB\nU1XXJTnV3b+346lXJDk18fhUkkumxE8P8Qz/Ppkk3f18km9X1Q/t0hYAAAAAOxybZ+Oq+oEkfyfb\nl2f9+/BSMwIAAABgLnMVeJL8x0leneTLw/Q6lyb53aq6Ktsjcy6b2PbSbI+8OT2s74xneO6VSb5e\nVceSXNDdz1TV6SQnJva5LMlvzUqqqu6YeLg15zkBABusqk7kxf0CAAB2qO7dJ7+tqlcn+Y3ufu2U\n5/4wyZXd/c1hkuVPZntS5EuSfDbJj3R3V9UXkrwzyckkn0nyy919f1XdmuS13f32qroxyfXdfeMw\nyfLvJPmxbI8Q+t0kPzbMx7Mzh+7u2hk7PJP6Vg7Tuex8rwBgXtO++w+rWed6cH0d390AcNAW7evs\nOoKnqu5O8sYkP1RVTyZ5d3d/dGKTf9+z6O5Hq+reJI8meT7Jrf1C9ejWJB9Lcl6S+7r7/iF+Z5JP\nVNXjSZ5JcuPQ1jer6ueTPDxs955pxR0AAAAA9jCCZ9MZwTMm/goIwP4ZwWMEDwAcZov2dea+ixYA\nAAAAm0WBBwAAAGDkFHgAAAAARm7e26TDvmzPGXA4mJMAAACATaHAwwE7LPUdtR0AAAA2h0u0AAAA\nAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4\nAAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZOgQcAAABg\n5BR4AAAAAEZOgQcAAABg5BR4AAAAAEZu1wJPVX2kqp6uqkcmYv+gqr5SVV+uql+tqgsmnrutqh6v\nqseq6uqJ+JVV9cjw3Psn4i+tqnuG+ENV9aqJ526uqq8Oy1uXd8oAAOc29Gv+YOjDfHLot1xUVQ8O\n/ZMHqurCHdsvpR8EADCvc43g+WiSa3bEHkjyo939Z5J8NcltSVJVVyS5IckVwz4fqKoa9vlgklu6\n+3iS41V1ps1bkjwzxN+X5L1DWxcleXeSNwzL7ZMdKACAVaqqVyf52SQ/1t2vTfL9SW5M8q4kD3b3\n5Uk+Nzxeaj8IAGARuxZ4uvvzSb61I/Zgd39vePiFJJcO69clubu7n+vuJ5J8LclVVXVxkvO7++Sw\n3ceTXD+sX5vkrmH9U0neNKz/VJIHuvvZ7n42yYM5u9AE7FNV9WFa1v16AofKd5I8l+QHqupYkh9I\n8vX8/+3df6zd913f8edrcRtSCDGGyfnl/Bh1thoFiYTFhYHIFpZ5EUusUSVGo2TUTFMMpFQM4bAN\nUrE/SCWWJZoSjRKIY5EQ06CQriaNm3JVNEjc0gIBJ4uDZhHfyE7nUJeWodnqe3+c751Pbq7j63Pv\n/Z7v93ufD+nofr+f8/2e9+d9vvfc87nv8/l+z5vHLrs4NaZZznGQJEnSWVvqNXg+AOxtli8GDo/d\ndxi4ZIH22aad5uerAFV1Ejie5Jvf5rEkLbsayE2Slk9VvQH8MvCXjAo7X6qqfcD6qjrabHYUWN8s\nL9c4aN3yZyNJklaDNZPumOTfA/+3qh5dxv5M2pe7x1ZnptQNSZK0ApJcD1zfcsxvBX4KuAI4DvxW\nkh8e36aqWpk9OH+cU1UzKx1TkiS1Z7nGOhMVeJL8a+Am3jyVeBbYMLZ+KaNPrGY5dRrXePvcPpcB\nrzXTny+oqmNJZnlzchuAT5+uP1V197z+LT4ZSZLUaU1BY2ZuPckvtBD2O4E/qKpjTczfBr4LOJLk\nwqo60px+9Xqz/XKNg96Y35H54xxJkjQsyzXWOetTtJoLA/4McEtV/e3YXU8B25K8M8mVwEZgf1Ud\nAb6cZHNzscH3A78zts/tzfL7GF2sEEYXcr4xydok3wT8U+CTZ9tXSZKkCb0EvDfJec345fuBA8DH\nOTV2uR14slleznGQJEnSWXvbGTxJHgO+D/iWJK8Cv8DoW7PeCexrZsr8YVXtqKoDSfYwGvycBHZU\n1dy05R3Aw8B5wN6qerppfwjYneQgcIzRt1NQVW8k+UXgs812H24utixJkrTiqupPkjwCfA74GvB5\n4FeA84E9SbYDh4Bbm+2XbRwkSZI0iZwae/RTkqqqzG8bzkVXg7l0UZj/e9dHQ3utDOGYSDqzhd77\nh+p0ubb393vl/7a29S2Iq+V3RpLUf5OOdSa+yLIkSZK0PFa6xmNtR5I0fBZ4pLM2GiS29YmjFmN4\nx8RPmiVJkiSdDQs80kSGUkcYUg1hKMcEhnVcJEmSJLXhrL9FS5IkSZIkSd1igUeSJEmSJKnnLPBI\nkiRJkiT1nAUeSZIkSZKknrPAI0mSJEmS1HMWeCRJkiRJknrOAo8kSZIkSVLPWeCRJEmSJEnquTXT\n7oAkSZKkxUlSbcSpqrQRR5K0fCzwSJIkSb2y0jUeazuS1EeeoiVJkiRJktRzzuCRJK2Itk4jaIun\nK0iSJKnLLPBIklbQEGo8o7rOUApWFqokSZKGyQKPJElnNIjaDl5XQ5Ikabi8Bo8kSZIkSVLPWeCR\nJEmSJEnqOQs8kiRJkiRJPWeBR5IkSZIkqecs8EiSJEmSJPWcBR5JkiRJkqSee9sCT5JfS3I0yQtj\nbeuS7EvycpJnkqwdu++uJAeTvJTkxrH2a5O80Nx331j7uUkeb9qfS3L52H23NzFeTvIjy5eyJEmr\nV5Lq+23az6EkSVIXnWkGz68DW+a17QT2VdVVwLPNOkk2AbcBm5p9HkiSZp8Hge1VtRHYmGTuMbcD\nx5r2e4F7msdaB/w8cF1z+4XxQpIkSZpUDeAmSZKk+d62wFNVvw/81bzmm4FdzfIuYGuzfAvwWFWd\nqKpDwCvA5iQXAedX1f5mu0fG9hl/rCeAG5rlfwY8U1VfqqovAft4a6FJkgZoVBef9gwJZ1lIkiRJ\n/bJmgn3WV9XRZvkosL5Zvhh4bmy7w8AlwIlmec5s007z81WAqjqZ5HiSb24e6/ACjyVJq8BQaiM5\n8yaSJK2QNj9sqCrf9CRN3SQFnv+vqvyUVpIkSVJHtfGvirUdSd0wSYHnaJILq+pIc/rV6037LLBh\nbLtLGc28mW2W57fP7XMZ8FqSNcAFVXUsySxw/dg+G4BPn65DSe4eW50524QkSVKXzeDbuyRJ0ttL\n1dtXtZNcAXy8qq5u1j/C6MLI9yTZCaytqp3NRZYfZXRR5EuATwHvbmb5PA/cCewHPgHcX1VPJ9kB\nXF1VdyTZBmytqm3NRZY/B1zDqCT+R8A1zfV45vev5k+JHM0qGsrEomAuXTOUPGA4uQwlDzCXLhpK\nHjCcXLJqTodYaJwz197WzISVfq7byWU4vzM+X4s3pNeJpNXldO//Z/K2M3iSPAZ8H/AtSV5l9M1W\nvwTsSbIdOATcClBVB5LsAQ4AJ4Eddap6tAN4GDgP2FtVTzftDwG7kxwEjgHbmsd6I8kvAp9ttvvw\nQsUdSZIkSZIkLWIGT9c5g6dPhpLLUPKA4eQylDzAXLpoKHnAcHJZPZ+WO4Nn2aIM5nfG52vxhvQ6\nkbS6TDqD522/Jl2SJEmSJEndt6Rv0ZIkSZL6oK1vfnUmhyRpWizwSJIkaRXw67IlScPmKVqSJEmS\nJEk9Z4FHkiRJkiSp5zxFS5IkSQsYnW7U1rVrJEnS0ljgkSRJ0ml43RpptWizmOvFyKWVYYFHkiRJ\nkoRFXanfvAaPJEnSApKsTfKxJC8mOZBkc5J1SfYleTnJM0nWjm1/V5KDSV5KcuNY+7VJXmjuu2+s\n/dwkjzftzyW5vO0ctZxOndK2krcpJylJ6jALPJIkSQu7D9hbVe8Bvh14CdgJ7Kuqq4Bnm3WSbAJu\nAzYBW4AHksx9TP0gsL2qNgIbk2xp2rcDx5r2e4F72klLK6dauEmStDALPJIkSfMkuQD43qr6NYCq\nOllVx4GbgV3NZruArc3yLcBjVXWiqg4BrwCbk1wEnF9V+5vtHhnbZ/yxngBuWMGUJEnSwFngkSRJ\neqsrgS8m+fUkn0/y0SRfD6yvqqPNNkeB9c3yxcDhsf0PA5cs0D7btNP8fBVGBSTgeJJ1K5KNJEka\nPAs8kiRJb7UGuAZ4oKquAb5KczrWnKrynBlJWqVW+npbXntLk/BbtCRJkt7qMHC4qj7brH8MuAs4\nkuTCqjrSnH71enP/LLBhbP9Lm8eYbZbnt8/tcxnwWpI1wAVV9cb8jiS5e2x1pqpmlpKYJGm5+K1j\nWh5JrgeuX/LjjD586q8kVVWZ3zacD9SCuXTNUPKA4eQylDzAXLpoKHnAcHIJ89/7VyRK8hngx6rq\n5QFPxAUAABfaSURBVKbI8q7mrmNVdU+SncDaqtrZXGT5UeA6RqdefQp4d1VVkueBO4H9wCeA+6vq\n6SQ7gKur6o4k24CtVbVtXh/eMs6Za1/5Y9nW70sbcYYSo6047fxDudKv4/b+J2jnb9JK8/k6Oz5f\nWkmne/8/E2fwSJIkLewngd9I8k7gL4AfBc4B9iTZDhwCbgWoqgNJ9gAHgJPAjjr1KdoO4GHgPEbf\nyvV00/4QsDvJQeAY8KbijjRdwygiSdJq4gyezhvKp60wnFyGkgcMJ5eh5AHm0kVDyQOGk8vq+TTT\nGTzGmE6cdmI4g6dbfL7Ojs+XVpIzeCRJkiRJUue0ebHo1VwQs8AjSZIkSZJWmBelXmkWeCRJkiS1\nzq9/lqTlZYFHkiRJ0hR4IWdJWk5/Z9odkCRJkiRJ0tJY4JEkSZIkSeo5CzySJEmSJEk9N3GBJ8ld\nSf48yQtJHk1ybpJ1SfYleTnJM0nWztv+YJKXktw41n5t8xgHk9w31n5ukseb9ueSXD55mpIkSZIk\nScM1UYEnyRXAvwGuqaqrgXOAbcBOYF9VXQU826yTZBNwG7AJ2AI8kGTuqmcPAturaiOwMcmWpn07\ncKxpvxe4Z5K+SpIkSZIkDd2kM3i+DJwA3pVkDfAu4DXgZmBXs80uYGuzfAvwWFWdqKpDwCvA5iQX\nAedX1f5mu0fG9hl/rCeAGybsqyRJkiRJ0pIlqZW+Tdq3ib4mvareSPLLwF8C/wf4ZFXtS7K+qo42\nmx0F1jfLFwPPjT3EYeASRkWiw2Pts007zc9Xm3gnkxxPsq6q3pikz5IkSZLUR0v5h0/SSljpl2TO\nvMkCJirwJPlW4KeAK4DjwG8l+eHxbapqSZUnSZIkSeqDdv7v6eY/lJK6Y6ICD/CdwB9U1TGAJL8N\nfBdwJMmFVXWkOf3q9Wb7WWDD2P6XMpq5M9ssz2+f2+cy4LXmNLALTjd7J8ndY6szE+YkSZI6aQbf\n3iV1m8UXSdM3aYHnJeA/JjkP+Fvg+4H9wFeB2xldEPl24Mlm+6eAR5P8Z0anXm0E9jezfL6cZHOz\n//uB+8f2uZ3RqV3vY3TR5gVV1d3j66eu3yxJkvrv+uY258PT6YYkSVKHTXoNnj9J8gjwOeBrwOeB\nXwHOB/Yk2Q4cAm5ttj+QZA9wADgJ7KiquTL3DuBh4Dxgb1U93bQ/BOxOchA4xuhbuiRJkiRJkjRP\nTtVZ+ilJVVXmt638NMm2BHPpmqHkAcPJZSh5gLl00VDygOHkEua/9w/VQuOcufZ2Tglp4/eljThD\nidFWHGN0L85QYoziDOFveHv/c/p8nWWkFX++2noPniSPSb8mXZIkSZIkSR0x6TV4JEmSJEnqHL/N\nWauVBR5JkiRJ0sD4zWZafSzwSJIkSZJa4wwbaWVY4JEkSZIktcjZNV1iwW04LPBIkiRJkrSqWXQb\nAr9FS5IkSZIkqecs8EiSJEmSJPWcp2hJkiRJktRRXiNnsUanga3m58sCjyRJkiRJneX1cRavjdpO\nd58vT9GSJEmSJEnqOQs8kiRJkiRJPWeBR5IkSZIkqecs8EiSJEmSJPWcBR5JkiRJkqSes8AjSZIk\nSZLUcxZ4JEmSJEmSes4CjyRJkiRJUs9Z4JEkSZIkSeo5CzySJEmSJEk9Z4FHkiRJkiSp5yzwSJIk\nSZIk9ZwFHkmSJEmSpJ6zwCNJkiRJktRzExd4kqxN8rEkLyY5kGRzknVJ9iV5OckzSdaObX9XkoNJ\nXkpy41j7tUleaO67b6z93CSPN+3PJbl88jQlSZLOTpJzknwhycebdcc5kiSps5Yyg+c+YG9VvQf4\nduAlYCewr6quAp5t1kmyCbgN2ARsAR5IkuZxHgS2V9VGYGOSLU37duBY034vcM8S+ipJknS2Pggc\nAKpZd5wjSZI6a6ICT5ILgO+tql8DqKqTVXUcuBnY1Wy2C9jaLN8CPFZVJ6rqEPAKsDnJRcD5VbW/\n2e6RsX3GH+sJ4IZJ+ipJknS2klwK3AT8KjBXrHGcI0mSOmvSGTxXAl9M8utJPp/ko0m+HlhfVUeb\nbY4C65vli4HDY/sfBi5ZoH22aaf5+SqMCkjA8STrJuyvJEnS2bgX+Bnga2NtjnMkSVJnTVrgWQNc\nAzxQVdcAX6WZpjynqopTU5olSZJ6IckPAK9X1Rc4NXvnTRznSJKkrlkz4X6HgcNV9dlm/WPAXcCR\nJBdW1ZFmWvLrzf2zwIax/S9tHmO2WZ7fPrfPZcBrSdYAF1TVGwt1JsndY6szE+YkSZI6aYaW396/\nG7g5yU3A1wHfmGQ3cLQL45yqmllKcpIkqWtmWI6xTkYfQE2wY/IZ4Meq6uVm4PGu5q5jVXVPkp3A\n2qra2Vx88FHgOkZTkj8FvLuqKsnzwJ3AfuATwP1V9XSSHcDVVXVHkm3A1qratkA/qqoyv204H6oF\nc+maoeQBw8llKHmAuXTRUPKA4eQS5r/3r1ik5PuAf1dV/yLJR+jAOGeufeWPZVu/L23EGUqMtuIY\no3txhhKjrThDidFWHGN0L85kY51JZ/AA/CTwG0neCfwF8KPAOcCeJNuBQ8CtAFV1IMkeRt9EcRLY\nUacqSzuAh4HzGH0r19NN+0PA7iQHgWPAWwY9kiRJLZgbs/wSjnMkSVJHTTyDpyucwdMnQ8llKHnA\ncHIZSh5gLl00lDxgOLm0N4Nn2pzBY4zpxDFG9+IMJUZbcYYSo604xuhenMnGOpNeZFmSJEmSJEkd\nYYFHkiRJkiSp5yzwSJIkSZIk9ZwFHkmSJEmSpJ6zwCNJkiRJktRzFngkSZIkSZJ6zgKPJEmSJElS\nz1ngkSRJkiRJ6jkLPJIkSZIkST1ngUeSJEmSJKnnLPBIkiRJkiT1nAUeSZIkSZKknrPAI0mSJEmS\n1HMWeCRJkiRJknrOAo8kSZIkSVLPWeCRJEmSJEnqOQs8kiRJkiRJPWeBR5IkSZIkqecs8EiSJEmS\nJPWcBR5JkiRJkqSes8AjSZIkSZLUcxZ4JEmSJEmSes4CjyRJkiRJUs9Z4JEkSZIkSeq5JRV4kpyT\n5AtJPt6sr0uyL8nLSZ5JsnZs27uSHEzyUpIbx9qvTfJCc999Y+3nJnm8aX8uyeVL6askSZIkSdJQ\nLXUGzweBA0A16zuBfVV1FfBss06STcBtwCZgC/BAkjT7PAhsr6qNwMYkW5r27cCxpv1e4J4l9lWS\nJEmSJGmQJi7wJLkUuAn4VWCuWHMzsKtZ3gVsbZZvAR6rqhNVdQh4Bdic5CLg/Kra32z3yNg+44/1\nBHDDpH2VJEmSJEkasqXM4LkX+Bnga2Nt66vqaLN8FFjfLF8MHB7b7jBwyQLts007zc9XAarqJHA8\nybol9FeSJEmSJGmQ1kyyU5IfAF6vqi8kuX6hbaqqktRC9y23JHePrc60EVOSJLVlBt/eJUmS3t5E\nBR7gu4Gbk9wEfB3wjUl2A0eTXFhVR5rTr15vtp8FNoztfymjmTuzzfL89rl9LgNeS7IGuKCq3lio\nM1V19/j6qcv7SJKk/ru+uc358HS6IUmS1GETnaJVVT9XVRuq6kpgG/Dpqno/8BRwe7PZ7cCTzfJT\nwLYk70xyJbAR2F9VR4AvJ9ncXHT5/cDvjO0z91jvY3TRZkmSJEmSJM0z6Qye+eZOxfolYE+S7cAh\n4FaAqjqQZA+jb9w6Ceyoqrl9dgAPA+cBe6vq6ab9IWB3koPAMUaFJEmSJEmSJM2TU3WWfkpSVZX5\nbadqTn0XzKVrhpIHDCeXoeQB5tJFQ8kDhpNLmP/eP1QLjXPm2lf+WLb1+9JGnKHEaCuOMboXZygx\n2oozlBhtxTFG9+JMNtZZyrdoSZIkSZIkqQMs8EiSJEmSJPWcBR5JkiRJkqSes8AjSZIkSZLUcxZ4\nJEmSJEmSes4CjyRJkiRJUs9Z4JEkSZIkSeo5CzySJEnzJNmQ5PeS/HmSP0tyZ9O+Lsm+JC8neSbJ\n2rF97kpyMMlLSW4ca782yQvNffeNtZ+b5PGm/bkkl7ebpSRJGhILPJIkSW91AvhQVX0b8F7gx5O8\nB9gJ7Kuqq4Bnm3WSbAJuAzYBW4AHkqR5rAeB7VW1EdiYZEvTvh041rTfC9zTTmqSJGmILPBIkiTN\nU1VHquqPm+WvAC8ClwA3A7uazXYBW5vlW4DHqupEVR0CXgE2J7kIOL+q9jfbPTK2z/hjPQHcsHIZ\nSZKkobPAI0mS9DaSXAF8B/A8sL6qjjZ3HQXWN8sXA4fHdjvMqCA0v322aaf5+SpAVZ0EjidZt/wZ\nSJKk1cACjyRJ0mkk+QZGs2s+WFV/PX5fVRVQU+mYJEnSPGum3QFJkqQuSvIORsWd3VX1ZNN8NMmF\nVXWkOf3q9aZ9FtgwtvuljGbuzDbL89vn9rkMeC3JGuCCqnpjgX7cPbY6U1UzS0pMkiR1zExzW5qM\nPnzqryRVVZnfNpwP1IK5dM1Q8oDh5DKUPMBcumgoecBwcgnz3/uXPcLoAsm7GF0E+UNj7R9p2u5J\nshNYW1U7m4ssPwpcx+jUq08B766qSvI8cCewH/gEcH9VPZ1kB3B1Vd2RZBuwtaq2zevHW8Y5c+0r\nfyzb+n1pI85QYrQVxxjdizOUGG3FGUqMtuIYo3txJhvrWODpvKEMxmE4uQwlDxhOLkPJA8yli4aS\nBwwnl1YKPN8DfAb4U049aXcxKtLsYTTz5hBwa1V9qdnn54APACcZndL1yab9WuBh4Dxgb1XNfeX6\nucBuRtf3OQZsay7QPN4PCzzGmEIcY3QvzlBitBVnKDHaimOM7sWxwPOmtmEMYGE4g3EYTi5DyQOG\nk8tQ8gBz6aKh5AHDyWXlCzxdYYHHGNOJY4zuxRlKjLbiDCVGW3GM0b04k411vMiyJEmSJElSz1ng\nkSRJkiRJ6jkLPJIkSZIkST1ngUeSJEmSJKnnLPBIkiRJkiT1nAUeSZIkSZKknrPAI0mSJEmS1HMT\nFXiSbEjye0n+PMmfJbmzaV+XZF+Sl5M8k2Tt2D53JTmY5KUkN461X5vkhea++8baz03yeNP+XJLL\nl5KoJEmSJEnSUE06g+cE8KGq+jbgvcCPJ3kPsBPYV1VXAc826yTZBNwGbAK2AA8kSfNYDwLbq2oj\nsDHJlqZ9O3Csab8XuGfCvkqSJEmSJA3aRAWeqjpSVX/cLH8FeBG4BLgZ2NVstgvY2izfAjxWVSeq\n6hDwCrA5yUXA+VW1v9nukbF9xh/rCeCGSfoqSZIkSZI0dEu+Bk+SK4DvAJ4H1lfV0eauo8D6Zvli\n4PDYbocZFYTmt8827TQ/XwWoqpPA8STrltpfSZIkSZKkoVlSgSfJNzCaXfPBqvrr8fuqqoBayuNL\nkiRJkiTpzNZMumOSdzAq7uyuqieb5qNJLqyqI83pV6837bPAhrHdL2U0c2e2WZ7fPrfPZcBrSdYA\nF1TVG6fpy91jqzOT5iRJkrpoBt/eJUmS3l5GE23OcqfRBZJ3MboI8ofG2j/StN2TZCewtqp2NhdZ\nfhS4jtGpV58C3l1VleR54E5gP/AJ4P6qejrJDuDqqrojyTZga1VtW6AvVVWZ3zacyUPBXLpmKHnA\ncHIZSh5gLl00lDxgOLmE+e/9Q7XQOGeufeWPZVu/L23EGUqMtuIYo3txhhKjrThDidFWHGN0L85k\nY51JCzzfA3wG+FNOZXYXoyLNHkYzbw4Bt1bVl5p9fg74AHCS0Sldn2zarwUeBs4D9lbV3Feunwvs\nZnR9n2PAtuYCzfP7YoGnN4aSy1DygOHkMpQ8wFy6aCh5wHByscBjgWe1xmgrjjG6F2coMdqKM5QY\nbcUxRvfitFjg6RILPH0ylFyGkgcMJ5eh5AHm0kVDyQOGk4sFHgs8qzVGW3GM0b04Q4nRVpyhxGgr\njjG6F2eysc6Sv0VLkiRJkiRJ02WBR5IkSZIkqecs8EiSJEmSJPWcBR5JkiRJkqSes8AjSZIkSZLU\ncxZ4JEmSJEmSes4CjyRJkiRJUs9Z4JEkSZIkSeo5CzySJEmSJEk9Z4FHkiRJkiSp5yzwSJIkSZIk\n9ZwFHkmSJEmSpJ6zwCNJkiRJktRzFngkSZIkSZJ6zgKPJEmSJElSz1ngkSRJkiRJ6jkLPJIkSZIk\nST1ngUeSJEmSJKnnLPBIkiRJkiT1nAUeSZIkSZKknrPAI0mSJEmS1HMWeCRJkiRJknrOAo8kSZIk\nSVLPWeCRJEmSJEnquc4XeJJsSfJSkoNJfnba/ZEkSVoujnMkSdJy6XSBJ8k5wH8FtgCbgB9K8p7p\n9kqLMzPtDmhBM9PugE5rZtod0IJmpt0BDZjjnL6bmXYHtKCZaXdAC5qZdgd0WjPT7oCWUacLPMB1\nwCtVdaiqTgC/Cdwy5T5pUWam3QEtaGbaHdBpzUy7A1rQzLQ7oGFznNNrM9PugBY0M+0OaEEz0+6A\nTmtm2h3QMup6gecS4NWx9cNNmyRJUt85zpEkScum6wWemnYHJEmSVojjHEmStGzWTLsDZzALbBhb\n38Do0603SbLAACkr1qn29TWXDy/Q1tdc5utzHvOPS59zGTeEPOaOzRBymTOUXIaSBwwrl95bwjgH\n2jmWbf2+9DWXab2n9vX5aivG+HHxmHQnRtvHZSgx2opjjG7GOTup6u6HR0nWAP8TuAF4DdgP/FBV\nvTjVjkmSJC2R4xxJkrScOj2Dp6pOJvkJ4JPAOcBDDnokSdIQOM6RJEnLqdMzeCRJkiRJknRmXb/I\nMgBJtiR5KcnBJD97mm3ub+7/kyTf0XYfV6szHZsk/yDJHyb52yQ/PY0+rkaLOC7/qnmt/GmS/5Hk\n26fRz9VmEcfllua4fCHJHyX5J9Po52q0mPeZZrt/mORkkn/ZZv9Wq0W8Zq5Pcrx5zXwhyX+YRj+X\ng2OdbnKc012OdbrJsU43Oc7prmUf61RVp2+Mpiy/AlwBvAP4Y+A987a5CdjbLG8Gnpt2v1fDbZHH\n5u8C3wn8J+Cnp93n1XBb5HH5LuCCZnmLr5nOHJevH1u+Gnhl2v1eDbfFHJux7T4N/HfgB6fd76Hf\nFvmauR54atp9bSlXxzrdPC6Oc7p7bBzrdPO4ONbp4HEZ285xTseOzdmOdfowg+c6Ri/8Q1V1AvhN\n4JZ529wM7AKoqueBtUnWt9vNVemMx6aqvlhVnwNOTKODq9RijssfVtXxZvV54NKW+7gaLea4fHVs\n9RuA/91i/1azxbzPAPwk8DHgi212bhVb7HHp5tdYnB3HOt3kOKe7HOt0k2OdbnKc013LPtbpQ4Hn\nEuDVsfXDTduZtvGP+MpbzLFR+872uGwH9q5ojwSLPC5JtiZ5Efhd4M6W+rbanfHYJLmE0Rvug02T\nF7BbeYt5zRTw3c10/71JNrXWu+XlWKebHOd0l2OdbnKs002Oc7pr2cc6nf4WrcZif7nmV7X8pVx5\nPsfdtOjjkuQfAx8A/tHKdUeNRR2XqnoSeDLJ9wK7gb+/or0SLO7Y/BdgZ1VVkjCMWSNdt5jj8nlg\nQ1X9TZJ/DjwJXLWy3VoRjnW6yee3uxzrdJNjnW5ynNNdyz7W6cMMnllgw9j6BkaVrbfb5tKmTStr\nMcdG7VvUcWkuNvhR4Oaq+quW+raandXrpap+H1iT5JtXumNa1LG5FvjNJP8L+EHggSQ3t9S/1eqM\nx6Wq/rqq/qZZ/l3gHUnWtdfFZeNYp5sc53SXY51ucqzTTY5zumvZxzp9KPB8DtiY5Iok7wRuA56a\nt81TwI8AJHkv8KWqOtpuN1elxRybOVaB23PG45LkMuC3gR+uqlem0MfVaDHH5VubT01Icg1AVR1r\nvaerzxmPTVX9vaq6sqquZHR++h1Vdbq/d1oei3nNrB97zVwHpKreaL+rS+ZYp5sc53SXY51ucqzT\nTY5zumvZxzqdP0Wrqk4m+Qngk4yuMv1QVb2Y5N829/+3qtqb5KYkrwBfBX50il1eNRZzbJJcCHwW\n+Ebga0k+CGyqqq9MreMDt5jjAvw88E3Ag83fixNVdd20+rwaLPK4/CDwI0lOAF8Btk2tw6vIIo+N\nWrbI4/I+4I4kJ4G/oaevGcc63eQ4p7sc63STY51ucpzTXSsx1knz1VuSJEmSJEnqqT6coiVJkiRJ\nkqS3YYFHkiRJkiSp5yzwSJIkSZIk9ZwFHkmSJEmSpJ6zwCNJkiRJktRzFngkSZIkSZJ6zgKPJEmS\nJElSz1ngkSRJkiRJ6rn/BwDM5Baqvv4fAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axs = plt.subplots(2, 2, sharex=True, figsize=(16, 9))\n", "\n", "xs, ys = zip(*approx.items())\n", "axs[0, 0].plot(xs, ys, '.')\n", "\n", "xs = list(approx.keys())\n", "xs.sort()\n", "ys = [approx[x] for x in xs]\n", "axs[0, 1].bar(xs, ys, 0.005)\n", "\n", "def get_bins(my_dict, nbins):\n", " acumulator = [0] * nbins\n", " xmin, xmax = xs[0], xs[-1]\n", " interval = (xmax - xmin) / nbins\n", " bin_xs = [xmin + i * interval + interval / 2 for i in range(nbins)]\n", " #the starting point, midpoint woudl take more interval / 2\n", " curr_bin = 0\n", " for x in xs:\n", " y_cnt = approx[x]\n", " while curr_bin + 1 != nbins and abs(x - bin_xs[curr_bin]) > abs(x - bin_xs[curr_bin + 1]):\n", " curr_bin += 1\n", " acumulator[curr_bin] += y_cnt\n", " return bin_xs, acumulator, interval\n", "\n", "bin_xs, acumulator, interval = get_bins(approx, 10)\n", "axs[1, 0].bar(np.array(bin_xs) - interval / 2, acumulator, 0.5 / 10)\n", "\n", "bin_xs, acumulator, interval = get_bins(approx, 20)\n", "axs[1, 1].bar(np.array(bin_xs) - interval / 2, acumulator, 0.5 / 20)\n", "axs[1, 1].set_xlim(0, 0.5)\n", "\n", "fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def compute_median_from_dictionary(my_dict):\n", " xs = list(my_dict.keys())\n", " xs.sort()\n", " x_cnt = [my_dict[x] for x in xs]\n", " start = 0\n", " end = len(xs) - 1\n", " while start != end:\n", " if start == end - 1 and x_cnt[start] == x_cnt[end]:\n", " return (xs[start] + xs[end]) / 2\n", " if x_cnt[start] > x_cnt[end]:\n", " x_cnt[start] -= x_cnt[end]\n", " end -= 1\n", " elif x_cnt[start] < x_cnt[end]:\n", " x_cnt[end] -= x_cnt[start]\n", " start += 1\n", " else:\n", " start += 1\n", " end -= 1\n", " return xs[start]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.2216\n" ] } ], "source": [ "print(compute_median_from_dictionary(approx))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6 0.2\n", "50 0.22\n", "326 0.222\n" ] } ], "source": [ "#do not put this, but discuss\n", "for my_round in [1, 2, 3]:\n", " my_approx = defaultdict(int)\n", " for val, cnt in approx.items():\n", " my_approx[round(val, my_round)] += cnt\n", " print(len(my_approx), compute_median_from_dictionary(my_approx))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Subsampling" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Mean #1Mean #2Mean #3Max #1Max #2Max #3
10.233000.39770.488600.23300.39770.4886
100.221600.19320.227250.38070.36930.3977
1000.264200.17330.215900.49430.48860.4830
10000.224450.21590.221600.50000.50000.5000
100000.221600.21590.221600.50000.50000.5000
\n", "
" ], "text/plain": [ " Mean #1 Mean #2 Mean #3 Max #1 Max #2 Max #3\n", "1 0.23300 0.3977 0.48860 0.2330 0.3977 0.4886\n", "10 0.22160 0.1932 0.22725 0.3807 0.3693 0.3977\n", "100 0.26420 0.1733 0.21590 0.4943 0.4886 0.4830\n", "1000 0.22445 0.2159 0.22160 0.5000 0.5000 0.5000\n", "10000 0.22160 0.2159 0.22160 0.5000 0.5000 0.5000" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import random\n", "import pandas as pd\n", "arr = np.ndarray(shape=(5, 6), dtype=float)\n", "samp_sizes = [1, 10, 100, 1000, 10000]\n", "for rep in range(3):\n", " for si, samp_size in enumerate([1, 10, 100, 1000, 10000]):\n", " my_vals = random.sample(mafs, samp_size)\n", " arr[si, rep] = np.median(my_vals)\n", " arr[si, rep + 3] = max(my_vals)\n", "df = pd.DataFrame(arr, index=samp_sizes,\n", " columns=['Mean #1', 'Mean #2', 'Mean #3', 'Max #1', 'Max #2', 'Max #3'])\n", "df" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "#bill gates effect" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }