{ "metadata": { "name": "", "signature": "sha256:8b89a9f1e0735686e77b8d9c337173dcf4b33f9119241cd493f3b26ac7ac8b5c" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Fitting power law data**\n", "\n", "Consider a random variable $x$ that takes values governed by a power law probability distribution, $p(x)=C x^{-\\alpha}$.\n", "\n", "First we generate some synthetic data with a power law distribution:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "alpha=2.5\n", "xmin=1\n", "n=1000\n", "xr=xmin*rand(n)**(-1/(alpha-1))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "[This uses the so-called \"transformation method\": if $r$ is a random real number uniformly distributed in the range $0 < r \\le 1$,\n", "then it is straightforward to show that $x = x_{\\rm min}r^{\u22121/(\\alpha\u22121)}$ is a random power-law-distributed real number in the range $x_{\\rm min} \\le x < \\infty$, with probability distribution governed by exponent $\\alpha$, i.e., $p(x)=C x^{-\\alpha}$.]\n", "\n", "Here is a histogram of the above 1000 data points, with number counts in bins of linear size .1:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(6,5))\n", "yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')\n", "xlabel('x')\n", "ylabel('# occurrences in bin of size .1')\n", "ylim(.5,300),xlim(1,200)\n", "xscale('log');" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFHCAYAAACs30uOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHPRJREFUeJzt3XuUZFV96PHvb14wyEUcMWJ4CMo7IRGNQIya9j3LB4gI\nhKveKxIMcCUYTdREHjVRSLxGsvABRCKgYQmIPAQSQaOWkmBAMCACo4wDyQwXEQUZZJxhZvzdP061\n3T3TXX26u6pOVZ3vZ61adc4+r18PRf1qn73P3pGZSJLqa17VAUiSqmUikKSaMxFIUs2ZCCSp5kwE\nklRzJgJJqrkFVQcwGxFhn1dJmoXMjM3LBrZGkJkD8zr99NMH6hqzPddMjiu773T7zXZ7L/6b9Ot/\n315cpxefobL7z3Wf2W7rx9dUBjYRDJKRkZGBusZszzWT48ruO91+c90+KHr1d3TqOr34DJXdf677\nDMtnqJ1olyX6VUTkIMat/tFoNGg0GlWHoQE3aJ+jiCCH6daQNBd1+JWn7huWz5E1AkmqCWsEkqRJ\nmQgkqeZMBJJUcyYCSao5E4Ek1ZyJQJJqzkQgSTVnIpCkmjMRSFLNmQgkqeZMBJJUcyYCSaq5gU0E\njUaDZrNZdRiS1PeazWbb4bIdfVSSasLRRyVJkzIRSFLNmQgkqeZMBJJUcyYCSao5E4Ek1ZyJQJJq\nzkQgSTVnIpCkmjMRSFLNmQgkqeZMBJJUcwObCJYs2fJ1zjlVRyVJg2dB1QHM1ooVE9c/+EFYt66a\nWCRpkA1sIliyZOL64sXVxCFJg25gbw1JkjrDRCBJNWcikKSaMxFIUs2ZCCSp5kwEklRzJgJJqjkT\ngSTVXN89UBYRhwKvA7YDPpOZX604JEkaan2XCDLzS8CXImJ74O8AE4EkdVFPbg1FxAUR8VBE3LlZ\n+dKIWB4R90bE+zc77BTgk72IT5LqrFdtBBcCS8cXRMR8ii/6pcB+wNERsW8UPgJ8OTNvn8lF1q6F\nRx7Z8iVJmlpPbg1l5o0RsdtmxQcCKzLzfoCIuBQ4FHgl8Apgu4jYIzP/ocw1Fi+Gs84qXmPXhU2b\nYM2aOf8JkjS0qmwj2AlYNW59NXBQZp4EfGK6gxuNxq+XR0ZGOOOMEc44Y+I+a9bAzjt3IlRJGjzN\nZpNmszntfpGZ3Y8GaNUIrs3M/VvrhwNLM/O41vpbGUsE050ry8Q9mgisEUgSRASZGZuXV/kcwQPA\nLuPWd6GoFUiSeqjKRHArsGdE7BYRi4CjgGsqjEeSaqlX3UcvAW4C9oqIVRFxTGZuBN4F3ADcDVyW\nmfeUPWej0Sh170uS6q7ZbE5oV91cz9oIOsk2AkmauX5sI5Ak9QETgSTVXN+NNdRpjz8OCxduWf7k\nkxBbVJAkqX6Guo0gEzZu3LJ80SL41a9MBJLqZao2goGtETQaDUZGRhgZGZlyn4jJawOSVCfTPWE8\n1DWCqY+3RiCpfuw1JEma1MDeGuqkBx+EW27Zsvx1r4MF/gtJGnJ+zQE33wzHHw8HHTRWdt118Nhj\nsO221cUlSb0wsImgTGPxTBx8MFx99di6CUDSsLCxeNLjJzYWX301XHTRlongxz82IUgaHh1tLI6I\nfeYekiSpH8z21tBXgF07GUiv3XXXWI3gv/+72lgkqUpTJoKIaDdd5NO6EEvP7LcfHHXUxLKXv7ya\nWCSpalO2EUTE48CfA+uB8TsF8LHMfHr3w5vcXNsIyrCNQNKwmc0QE7cC38/Mf5/kZI0OxiZJqlC7\nRHA4sG6yDZm5W1eimYFOdx+VpGFl99FZ8taQpGHT6e6jy+YekiSpH8x20LlbOxqFJKkys0oEmXlt\npwORJFXDYaglqeZMBJJUcyYCSaq5aRNBROwYEZ+JiOtb6/tFxLHdD629RqPRtl+sJKnQbDZpNBpT\nbp/2OYJWArgQ+GBm/k5ELAT+MzN/u5OBzoTPEUjSzM3lOYIdMvMyYBNAZm4ANnY4PklSRcokgl9E\nxK8HmIuIg4HHuheSJKmXysxH8F7gWuA5EXET8AzgzV2NSpLUM9Mmgsy8LSJeCuxDMQT1DzLzya5H\nJknqiTKNxSuBj2bmuePKrsvM13c7uDYx9aSxeNddYf78sbLrr4eddoK1a+Gggybu/5a3wAc+0NWQ\nJGlOZjMfwagNwEhEHAgcn5nrgZ06HWC/ufnmYoL7Ua95DSxdClttBZs2wYoVcMstxbaLL4YHH6wm\nTkmaqzKJYG1mHhUR7wO+FRFHdjuofvBbvzVx/frr4clxN8TmzYP99y+Wf/M3YeXK3sUmSZ1UevL6\nzPy/EfFdionrl3QvpP70O79TdQSS1B1lEsFpowuZ+a8R8Wrgf3cvpHKcoUySypn1DGURsW9m3hMR\nL2Di5PWjx93WsShnqBeNxTNx9tnFraGzz646Ekma2mwai98DHAd8jC0TAcDLOhSbJKlCUyaCzDyu\n9T7Ss2gkST1XZvTRIyJiu9byqRFxZUQ8v/uhSZJ6ocxYQ6dl5pqIeDHwCuAC4LzuhiVJ6pUyiWBT\n6/31wPmZeR2wsHshSZJ6qUwieCAiPg0cBfxzRGxd8jhJ0gAo84V+JHAD8OrM/DnwNOAvuhqVJKln\nyow++gRwxbj1BwFH1pGkIeEtHkmquSlrBBGxdWau62UwwyAT/uu/JpbtuCNsvXU18UjSdNrVCG4C\niIiLexTLjDQajbZjZ1Rl0ybYfXcYGSlee+8N//mfVUclqc6azSaNRmPK7e3GGroLOBP4EPDnFLOT\njcrMvLJzYc5MP481tHFj8et/48Zi2+//Ppx1VvEuSVWazVhDxwNvAZ4KvGGS7ZUlAklS57Qba+hG\n4MaIuDUz/7GHMUmSeqjMfASfi4iTgZe21pvAeZm5oWtRDanPfQ7uuWfybaefboOypGqU6T56LvB8\n4FPAOcALWmWaoSuugNWrYbvtJr7OOmviNJiS1EtlagQvzMzxEzV+LSK+162ABtUFF8BVVxXdR9t5\n85vh0EMnlv3N33QvLkmaTplEsDEi9sjMFQAR8VxgY3fDGizHHAOHHVZ1FJI0O2USwV8AX4+I+1rr\nuwHHdC2iATR6i0eSBlGZsYa+FhF7AXtTTFn5Q584lqThUaZGQOuL/44uxyJJqoCDzklSzZkIJKnm\nSt0aioidKBqJ51OMOZSZ+a0uxiVJ6pFpE0FEfIRimsq7GZu/GMBEIElDoEyN4DBg78xc3+1g6ux9\n74OttoI99oCTTqo6Gkl1UqaN4EfAom4HUmcf+Qjsu28xl8GVjukqqcfK1Ah+CdweEV8DRmsFmZl/\n2r2w6uWEE4r3b3wD7rqr2lgk1U+ZRHBN6zVe5bPCNBoNRkZGGBkZqToUSeprzWaz7YyOZZ4svqiD\n8XRMu2nX+s0998DChfDoo1VHIqmORn80L1u2bNLt7Savvzwzj4iIOyfZnJuNSKop7LMPnHPO2Pr2\n21cXiyRNpl2N4OTW+2TTVKqkCy+sOgJJam/KXkOZ+f9a7/cD64DfBfYH1rXKJElDYNruoxHxx8At\nwJuANwM3R8Sx3Q5MktQbZXoNvQ84IDN/BhARTwe+DXymm4FJknqjzANlPwV+MW79F60ySdIQaNdr\n6L2txRUUt4Oubq0fCjhnsSQNiXa3hv4HxYNjPwJWMvYQ2ZfogwfKJEmdMWUiyMxGD+OQJFXEiWkk\nqeZMBH3quONg552L19KlVUcjaZiZCPrUI4/AaafB+efDT35SdTSShlmZGcp+AziOYqrK0f0zM9/R\nxbgEPP3p8MxnVh2FpGFX5oGyL1FMS/lV4FetMnsN9ZEf/Qhe9rJieckSuP32auORNFjKJILFmfn+\nrkeiWdu4ERYsgCuugFe/uupoJA2aMm0E10XE67oeieZk4cKiYVmSZqpMIng3cG1ErIuIx1uvNd0O\nTJLUG2VmKNu2F4FIkqrRbqyhfTPznoh4/mTbM/O73QtLktQr7WoE76HoNnoWk/cSellXItKUVq6E\n9euL5Wc8A3bYYct9Nm4s5kgetXAh7LFHb+KTNJjajTV0XOt9pGfRqK03vQkeewzWrIE/+zM45ZSJ\n2+fPhx13LPYDePJJyCwSiCRNpUz3UfWRq66CL35x8m1LlkysDaxcCa98ZW/ikjS4HGJCkmqur2oE\nEbE78EHgqZl5RNXx9LvVq+Eb34AHH6w6EkmDrMzk9S+OiG1by2+LiLMi4tndCCYz78vMP+7GuYfR\n7bfDBz4Ad94JhxxSdTSSBlWZW0PnAk9ExO9S9CT6EfC5sheIiAsi4qGIuHOz8qURsTwi7o0Ih7CY\npQMOgH/6J/joR6uORNKgKpMINmZmAm8EPpWZn6KYxrKsC4EJI+pHxHzgk63y/YCjI2LfGZxTktQh\nZRLB4xHxV8BbKcYdmg8sLHuBzLwReHSz4gOBFZl5f2ZuAC4FDo2IJRFxHvA8awljnngCvvMdWLu2\n6kgkDaMyjcVHAUcD78jMH0fErsBcb0TsBKwat74aOCgzHwGOL3OCRqPx6+WRkRFGRkbmGFJ/2mYb\n2G47OPFEeOpTYfHiqiOSNCiazSbNZnPa/cqMNfRgRFwJjD6f+lPg6jlF14H5DMYngmG2zz5FbUCS\nZmrzH8nLli2bdL8yvYbeCVwO/EOraGfgqjnG9wCwy7j1XShqBZKkHivTRvB/gBcDawAy84fAb8zx\nurcCe0bEbhGxiOL20zVzPKckaRbKJIL1mbl+dCUiFjCDWzsRcQlwE7BXRKyKiGMycyPwLuAG4G7g\nssy8p9156mL5cjj2WLjttvb7XXUVnH12b2KSNNzKNBZ/MyI+CGwTEa8CTgSuLXuBzDx6ivIvA18u\ne57NNRqNoWsk3ntv+PCHi+UXvQieP+kA4PDGN8LuuxfLu+wy+T6SNGq6RuMoHhGYWqu76LHA6Gy4\nNwD/mNMd2EURUeXlB8booHOOPioJICLIzNi8vEyNYGvgM5n56daJ5gOLAXu1S9IQKNNG8HWKL/5R\n2wD/2p1wJEm9ViYRbJWZvxhdyczHKZKBBsCjj8KZZ8LFF1cdiaR+VSYRPBERLxhdiYjfA37ZvZDK\naTQapZ6Yq7Ptt4cTToC77oLzz686GklVaTabbR/CLdNY/EKKsYBGR71/FnBUZt7aoRhnzMbimfnm\nN+G004p3SfU168bizPxOa2TQvSmeH/hBa6A4SdIQKDtD2e8Bu7f2f34rq5Sek0CS1L+mTQQRcTHw\nHOB2YNO4TSYCSRoCZWoELwD286Z8fXzlK3DhhcXyYYcVw1lAUbb11tXFJak7yvQa+j5FA3FfsddQ\n96xYAT/9KVx/fTEf8k9+AldeCRs3Vh2ZpNmYrtdQmRrBM4C7I+IWYHTwuczMSqdLr8t8BFXZc0+4\n++6x5ZtvrjYeSbM3Oi7bVPMRlEkEjdZ7AjFuWZI0BMp0H21GxG7AHpn5rxGxTZnjJEmDoaoZyiRJ\nfaKqGcokSX2i6zOUqf+ccAIsWVK8rrtubPnTny5/jle9qjjm1FMnln/2s0X5H/5hZ2OW1D1lEsHm\nM5RdzgxmKOsWu4/O3hNPFDOhLVwIGzbAc54Dhx8O69dPf+yoNWvgJS+BX242/OD69fDsZxfbJfWH\n6bqPlkkE7wceBu4E/gT4F+CUTgQ3F6NTVWp2tt0W5rX+6y9YMLsHxRYtmrx8gV0JpL4yMjIy++cI\nWreBvp+Z+wAzuHEgSRoUbWsEmbkR+EFEPLtH8UiSeqxMJX4JcFfryeInWmWVP1ksSeqMMongFMae\nKB5lryFJGhJl2gg+nZl79ygeSVKPlWkjWG4bweD71reKXj4nnTSz4844o9x+e+5ZnP9zzlIhDZyB\nbSMY7T5qF9LpveQlRf/+886D5cvLH3fffcX7vHlw8cXt992wAQ4+2KGqpX7UbDbbPndVJhGcOv0u\nvecw1OXNm1f8Wp9p//6pnhNodx1J/WfOw1BnZrPTQUmS+keZOYt/wVgvoUXAQuAXmbldNwOTJPVG\nmRrBtqPLETEPOAQ4uJtBSZJ6Z0Z3dTPzV5l5NbC0S/FIknqszK2hw8etzgNeAPxyit0lSQOmTD+S\nNzDWRrARuB84tFsBSZJ6q0wbwdt7EIcGwI03wjbbwOLFY2Vr18I3vwnr1lUXl6S5KXNr6LPAyZn5\n89b604CPZeY7uh1cOz5Q1lsvfjH87d8WX/w//3kxCxnAqlWwdCkceCBstRU8+WS1cUra0nQPlJVp\nLP7d0SQAkJmPAs+fe2hz48Q0vXX99cUv/89/fsttu+5abHvWs3ofl6TpTTcxTZlEEBGxZNzKEmD+\n3EOTJPWDMo3FHwO+HRFfoBiO+gig5FBkkqR+V6ax+HMRcRvwcoreQ4dl5t1dj0yS1BNlGosPBu7O\nzE+01reLiIMy8+auRydJ6roybQTnAY+PW3+iVSZJGgKlhpjIzBy3vAkbiyVpaJRJBPdFxJ9GxMKI\nWBQRJwMrux2YJKk3yiSC44E/AB4AVlOMPPrObgYlSeqdMr2GHgKO6kEskqQKTFsjiIhdIuKqiHi4\n9boiInbuRXCSpO4rc2voQuAa4Ddbr2tbZRoCZ545cf1Nb4Jzzy1//C23zPyazSYcckhxrUMOgVOn\nmBX74x8vtt9228yv0c6aNfDKV8JnPjOz4y67rDhuw4bOxiNVrcyTxc/IzPFf/BdFxJ91K6CyHHRu\n7i66CNavh6c9Db7whaLsuuvgkktgzz2nP/6Nb4S3vQ2e8pSZXfeBB+ArX4HM4kv1l1PMbnHHHXDt\ntXDiiTM7/3Q2bICvfQ32339mx917b3HcWB86aTBMN+hcmUTws4h4G/B5iiEm/gj4aUeim4N2Ayip\nnNe8Zmx5NBEAvOENsGjR9MfvsQccdNDsrr1gQfGFHDG74zuhymtLvTT6o3nZsmWTbi9za+gdwJHA\nj4EHKcYaOqZjEUqSKlWm19D9FLOUSZKG0Iwmr5ckDR8TgSTVnIlAkmquzANlp4xb3rq74UiSem3K\nRBARH4iIF1H0Ehp1U/dDkiT1UrteQ8spksDuEfFvwD3ADhGxT2Yu70l0kqSua3dr6OfAXwI/AkaA\nj1NMVfn+iPh290OTJPVCuxrBa4BTgedSTGD/PWBtZvowmSQNkSlrBJn5l5n5CuA+4J8oksYOEfHv\nEXFtrwKUJHVXmbGGbsjMW4FbI+L4zPyDiHhGtwOTJPXGtN1HM/N941bf3ip7uFsBSZJ6a0YPlGXm\nHd0KRJJUjTK3hjREfvWr4tVtmbBpE8yfP/HameWuP1Wco+fILIaRnjdv7sNJb9pUvE91rtFrzsbo\nseP/HaR+M7BDTDQajbYTLWhLEXD++cXEM90ei/9739tyToMPfaiYh+Ckk6Y//pRTiolzNnfuucU5\nFi4s3tevn3uso+f67ncn337yycX22Xj962HbbWcfm9QJzWaz7RwuA50InJ1sZk44ofj1u2lTMbNY\ntxxwwNTTOT796eXPc8YZEyfP6abnPW/qbaM1htlYvRrWrZv98VInjIyMDGcikCR1holAkmrORCBJ\nNWcikKSaMxFIUs2ZCCSp5kwEklRzJgJJqjkTgSTVnIlAkmrORCBJNWcikKSaMxFIUs2ZCCSp5kwE\nklRzJgJJqjkTgSTVnIlAkmrORCBJNWcikKSaMxFIUs2ZCCSp5kwEklRzC6oOYLyIeApwDrAeaGbm\n5ysOSZKGXr/VCN4EfCEz3wkcUnUwklQHXU8EEXFBRDwUEXduVr40IpZHxL0R8f5W8U7Aqtbypm7H\nJknqTY3gQmDp+IKImA98slW+H3B0ROwLrAZ26WFsklR7Xf+yzcwbgUc3Kz4QWJGZ92fmBuBS4FDg\nSuDwiDgHuKbbsUmSqvvVPf4WEBQ1gZ0yc21mviMzT8zMSyqKrbY+9jHYsGH6/VatgrvuKnfOI46A\ngw8uXuefP3Hbk08W7//xH8X297538nO8+91j5/joRydue+lL4cwz4cgjYdGiYp8vfanYduKJsP/+\nxfLee0NEsf2008rF/q53wWtfO7Z+3nkTtz/xBDz3ufDFL2557LHHFtefzGGHwRlnwIoVxfHLlxfl\ny5bBW99aLN97b7Ft1aqJxx5xBGy/PZxwwpbnPfVUeNGLxtY/+9nib/7rv27/d45673uLf8uZuv12\nOOCAsfUjj4Qrrih//Fe/Cm94w8yvO513vxsuuqjz5x1GVfUayrmeoNFo/Hp5ZGSEkZGRuZ6y1t7z\nHnjLW4rlBW0+FTvvDM1msbzfftOf94474MMfhl13Ldavvrr4Qn3JS+D00+GWW+Dkk+Hmm2GHHbY8\n/u//Hh57bGJZowE33FBse+ghWLkSvv71IondfDP85CfFfrfcAt//frH8wx+O/W0rVkwfN8C//Avc\nd9/Yl+5rXzsxmW3aVFz7oYe2PPaCC4ov4d/+7S23XX01/OxncMghxfFr1xblV1wBd94JF19clK1c\nCevWTTx2NOmcdx6ce+7Ebd/5Dnz722Pr//ZvxfunPlUu+Z11Fuy4I/zVX02/73gPP1wkg1GXXw7P\nfCYcfni542+7Da67bmbXLOPss4vP2dvf3vlzD4pms0lz9H/YNqpKBA8w1hZAa3n1TE4wPhFo7nbf\nvXhNZ/Hi4lf1TDzvebDXXsXybbcV70uWFOd5/PH2x+677+SxAhx0UFEzefjhyWsyW221ZdmSJbB+\nfbm4Fy6cuL7LLluWtbNo0dTbIsqVzcS8Ker3cz2vBtfmP5KXLVs26X5V3Rq6FdgzInaLiEXAUdgm\nIEmV6EX30UuAm4C9ImJVRByTmRuBdwE3AHcDl2XmPd2ORZK0pa7fGsrMo6co/zLw5dmet9Fo2DYg\nSSVM11bQV0NMzIRtBJJUzuiP5n5rI5Ak9QkTgSTVnIlAkmpuYBNBo9Eo9aCEJNVds9ls265qY7Ek\nDTkbiyVJbZkIJKnmTASSVHMmAkmquYFNBPYakqRy7DUkSTVnryFJUlsmAkmqOROBJNWciUCSas5E\nIEk1N7CJwO6jklSO3UclqebsPipJastEIEk1ZyJQTTWrDkBDYFjaKU0Eqqlm1QFoCJgIVFovPiyd\nvMZsz/XjH5c/7qGHyu07XSw/+MF05yl3nX7Xqy+c9es7c53Vq2d3npn+nWX2n+s+w/Jl346JoAdM\nBFsyEcxMr76MnnyyM9d54IHZncdEUI3IzKpjmLGIGLygJakPZGZsXjaQiUCS1DneGpKkmjMRSFLN\nmQgkqeZMBJJUc0ORCCLiKRHx2Yj4dET8z6rj0eCJiN0j4h8j4vKqY9HgiohDW99Dl0bEq6qOp6yh\n6DUUEW8DHsnMf46ISzPzj6qOSYMpIi7PzCOqjkODLSK2B/4uM/+46ljK6NsaQURcEBEPRcSdm5Uv\njYjlEXFvRLy/VbwTsKq1vKmngapvzfAzJE1qlp+jU4BP9i7KuenbRABcCCwdXxAR8yn+cZcC+wFH\nR8S+wGpgl9Zu/fw3qbdm8hmSplL6cxSFjwBfzszbex/q7PTtl2Zm3gg8ulnxgcCKzLw/MzcAlwKH\nAlcCh0fEOcA1vY1U/Womn6GIWBIR5wHPs5ag8Wb4XfQu4BXAmyPiT3ob6ewN2gxl428BQVETOCgz\n1wLvqCYkDZipPkOPAMdXE5IG0FSfo5OAT1QT0uz1bY1gCoPfsq2q+RlSJwzV52jQEsEDjLUF0Fpe\nXVEsGkx+htQJQ/U5GrREcCuwZ0TsFhGLgKOwTUAz42dInTBUn6O+TQQRcQlwE7BXRKyKiGMycyNF\nY8wNwN3AZZl5T5Vxqn/5GVIn1OFzNBQPlEmSZq9vawSSpN4wEUhSzZkIJKnmTASSVHMmAkmqOROB\nJNWciUCSas5EIEk1ZyKQpJozEUgdEBEvjIg7ImKr1hza34+I/aqOSyrDISakDomIDwFbA4uBVZn5\nkYpDkkoxEUgdEhELKUal/CXw++n/XBoQ3hqSOmcH4CnAthS1AmkgWCOQOiQirgE+DzwHeFZr2kKp\n7w3anMVSX4qI/wWsz8xLI2IecFNEjGRms+LQpGlZI5CkmrONQJJqzkQgSTVnIpCkmjMRSFLNmQgk\nqeZMBJJUcyYCSao5E4Ek1dz/B5HIiMJi2/5AAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "An obvious way to try to fit a power law to this data is to use the `linregress()` function, as described in [lec19](http://nbviewer.ipython.org/url/courses.cit.cornell.edu/info2950_2014fa/resources/lec19.ipynb),\n", "applied to the log-log data (where we're careful only to use the bins with non-zero counts):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from scipy.stats import linregress\n", "slope, intercept, r_value, p_value, std_err = linregress(log(bh[yh>0]),log(yh[yh>0]))\n", "print slope,intercept\n", "\n", "figure(figsize=(6,5))\n", "yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')\n", "xlabel('x')\n", "ylabel('# occurrences in bin of size .1')\n", "ylim(.5,300),xlim(1,200)\n", "xscale('log')\n", "\n", "xl=arange(1,200,.1)\n", "plot(xl,exp(intercept)*xl**(slope),'k-',linewidth=2,label='$\\\\alpha={:.2f}$'.format(-slope))\n", "\n", "legend();" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "-1.27505386327 3.51358457684\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFHCAYAAACs30uOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOX5xvHvE0KQfUeUHUHQH1gUqQmbCYpsAcqmVaSt\nCFS0allaEbfggrUVXCvghqAVwVo2S1JBOwhNkBYrKAUFEdlRQEGDYX1/f0wiAbJMlpmTmbk/1zVX\n5ux3cJwn55z3vK855xARkegV43UAERHxlgqBiEiUUyEQEYlyKgQiIlFOhUBEJMqpEIiIRLlYrwMU\nh5mpzauISDE45+zMeWF7RuCcC5vXAw88EFbHKO6+irJdoOsWtl5xl4fiv0lZ/e8biuOE4jMU6Pol\nXae4y8riKz9hWwjCSWJiYlgdo7j7Ksp2ga5b2HolXR4uQvV7lNZxQvEZCnT9kq4TKZ+hglhBVaKs\nMjMXjrml7EhJSSElJcXrGBLmwu1zZGa4SLo0JFIS0fBXngRfpHyOdEYgIhIl8jsjCMtWQyISGczO\n+k6SUlKUP5ZVCETEUzq7L31FLbC6RyAiEuVUCEREopwKgYhIlFMhEBGJcrpZLCISQh999BGvvfYa\njz/+eJ7LX3/9dXbv3s3q1asZMGAAP//5zwFYvHgxO3bsICsriyZNmjBw4MBSy6RCICISIlOnTmXl\nypVUr149z+WbN29m//79jBs3jn379tGyZUuuuOIKYmNj+fTTTxk/fjwAI0aM4JprrqFKlSqlkkuX\nhkREQmTs2LH0798/3+Xr16/nj3/8IwB16tShRYsWrFmzhn379rFs2TKOHj0KQOXKlYmLiyu1XDoj\nEBEppi1btvDCCy/kuzw+Pv6sL/6Cnpvo3bs3qampP663e/duWrRoQbt27Th58iQdOnRg1KhRXHPN\nNSoEIhL5SvOp4+I+tLZr1y5mzpxJu3bteP/99xk9ejR16tTh+++/p379+jRv3pxHH320SPss6Pcq\nX748bdq0AeDvf/87l19+Oe3atQNgwoQJPProo4wfP54nn3yyWL9PflQIRETykJmZyYABA1iyZAm1\na9emXr16jBkzhmHDhpGcnFzs/QZSlL799lteeeUVXnvtNQA+++wzfD4fS5cuZdmyZdx00020bduW\njh07FjtHbmFbCFJSUkhMTIyY3v9E5HRedz0xd+5c2rdvT+3atQGoW7cu69evx8x+vCxTnEtDhZ3p\nOOf4wx/+wIsvvkiVKlXYunUrixcvZsiQIQBcffXVzJo1i5UrVwZcCHw+Hz6fL9/lYV0IRESC5dix\nY7Ro0eLH6czMTMqVK8eAAQN+nFecS0N5FbjPP/+c5s2bY2Y888wzDBkyhKysLFavXs0PP/xAs2bN\n+OSTT2jbti0AR44cIT4+PuBj5vzRPGnSpDyXqxtqEfFMdrfIXsfI06FDh3jsscfo1KkTx44do3Ll\nyrz00kt0796dn//851SqVKnI+3z22WeZN28e27dv51e/+hVjxoyhWrVqXHbZZbz00ktkZmZy5ZVX\n/vhvYmZs27aNBg0a8NRTT5GZmUnlypWpUaMGv/zlL/M9Tn7/rvl1Q61CICKeKcuFIJwVtRDoOQIR\nkSinQiAiEuVUCEREopwKgYhIlFMhEBGJcioEIiJRToVARCTKhe2TxSISGUqzczkpHhUCEfGMHiYr\nG3RpSEQkyoVtIahV6+zXc895nUpEJPyE7aWhzZtPn77nHsjK8iaLiEg4C9tCUKvW6dMVK3qTQ0Qk\n3IXtpSERESkdKgQiIlFOhUBEJMqpEIiIRDkVAhGRKKdCICIS5VQIRESinAqBiEiUK3MPlJlZf6AP\nUA14yTm31ONIIiIRrcwVAufcQmChmdUAHgdUCEREgigkl4bM7GUz22tmH58xv6eZbTSzTWZ21xmb\n3Qs8G4p8IiLRLFT3CGYCPXPPMLNy+L/oewIXA9eb2UXm9xiQ6pz7qCgHOXwYDhw4+yUiIvkLyaUh\n59wKM2t6xuyfApudc1sBzOwNoD9wNXAVUM3MWjjnZgRyjIoVYepU/+vUceHECTh0qMS/gohIxPLy\nHkEDYHuu6R3AFc6524FnCts4JSXlx/eJiYk88kgijzxy+jqHDkHDhqURVUQk/Ph8Pnw+X6HrWaiG\niss+I1jsnGubPT0I6OmcG5k9fSOnCkFh+3KB5M4pBDojEBHxjw/tnDtrkGgvnyPYCTTKNd0I/1mB\niIiEkJeF4D9ASzNramZxwHXAIg/ziIhEpVA1H50DpAMXmtl2M7vJOXcc+A3wD+B/wFzn3IZA95mS\nkhLQtS8RkWjn8/lOu696ppDdIyhNukcgIlJ0ZfEegYiIlAEqBCIiUa7M9TVU2r77DsqXP3v+0aNg\nZ50giYhEn4i+R+AcHD9+9vy4ODh5UoVARKJLfvcIwvaMICUlhcTERBITE/NdxyzvswERkWhS2BPG\nEX1GkP/2OiMQkeijVkMiIpKnsL00VJp274bVq8+e36cPxOpfSEQinL7mgA8+gFtugSuuODXv7bfh\n4EGoUsW7XCIioRC2hSCQm8VFER8PCxacmlYBEJFIoZvFeW5/+s3iBQvglVfOLgR79qggiEjkKNWb\nxWbWuuSRRESkLCjupaF3gMalGSTU1q8/dUawbZu3WUREvJRvITCzgoaLrBmELCFz8cVw3XWnz+vW\nzZssIiJey/cegZl9B4wHjgC5VzJginOudvDj5a2k9wgCoXsEIhJpitPFxH+AT5xz/8pjZymlmE1E\nRDxUUCEYBGTltcA51zQoaYqgtJuPiohEKjUfLSZdGhKRSFPazUcnlTySiIiUBcXtdO4/pZpCREQ8\nU6xC4JxbXNpBRETEG+qGWkQkyqkQiIhEORUCEZEoV2ghMLP6ZvaSmaVlT19sZjcHP1rBUlJSCmwX\nKyIifj6fj5SUlHyXF/ocQXYBmAnc45y7xMzKA/91zrUpzaBFoecIRESKriTPEdRxzs0FTgA4544B\nx0s5n4iIeCSQQvC9mf3YwZyZxQMHgxdJRERCKZDxCMYBi4HmZpYO1AUGBzWViIiETKGFwDm3xsy6\nAq3xd0H9qXPuaNCTiYhISARys3gL8Cfn3LRc8952ziUHO1wBmUJys7hxYyhX7tS8tDRo0AAOH4Yr\nrjh9/aFDYcKEoEYSESmR4oxHkOMYkGhmPwVucc4dARqUdsCy5oMP/APc5+jRA3r2hAoV4MQJ2LwZ\nVq/2L3vtNdi925ucIiIlFUghOOycu87Mfg+8b2bXBjtUWfB//3f6dFoaHM11QSwmBtq29b8//3zY\nsiV02URESlPAg9c75/5oZh/iH7i+VvAilU2XXOJ1AhGR4AikENyf88Y5t8zMrgF+GbxIgdEIZSIi\ngSn2CGVmdpFzboOZtef0wetztltTaimLKBQ3i4viqaf8l4aeesrrJCIi+SvOzeKxwEhgCmcXAoCk\nUsomIiIeyrcQOOdGZv9MDFkaEREJuUB6Hx1iZtWy399nZn8zs8uCH01EREIhkL6G7nfOHTKzzsBV\nwMvA9ODGEhGRUAmkEJzI/pkMvOCcexsoH7xIIiISSoEUgp1m9jxwHfB3MzsnwO2C6tFHH6UstRwS\nEQlXgXyhXwv8A7jGOfctUBP4XVBTBWDixInceuutHD+uoRFEREqi0ELgnMt0zr3lnNuUPb3bOfdO\n8KMVrEKFCkyfPp0BAwaQmZnpdRwRkbDl+SWe4nr33XepVasWb7/9NklJSezdu9frSCIiYSnfQpB9\nL6DM6tSpE+np6TRr1ox///vfJCQk8Nlnn3kdC+dg69bTX1lZ3mYSESlIQWcE6QBm9lqIshRJSkoK\nu3fvJiMjg/bt2/PFF1/QsWNH0tPTPc114gQ0awaJif5Xq1bw3/96GklEopzP5yMlJSXf5QX1NbQe\nmAw8BIzHPzpZDuec+1vpxSyaM/sa+v7777nuuutYsmQJ55xzDq+//joDBgwIWZ7cfQ0dPw7nnOP/\nCZCQAFOn+n+KiHgpv76GCjojuAXoAlQH+uJ/jiDn1TcYIYurSpUqLFy4kFGjRpGVlcWgQYN4+umn\nvY4lIhIWCupraAWwwsz+45x7MYSZiiU2Npbp06fTpEkT7rnnHu68806+/PJL/vSnPxETE7b3xEVE\ngi6Q8Qhmm9mdQNfsaR8w3Tl3LGipisnMmDhxIo0aNWL48OFMnTqV7du3M3v2bM45x/t737Nnw4YN\neS974AH/JSURkVAL5E/lacBlwJ+B54D22fPKrGHDhpGWlka1atV488036d69OwcOHPA6Fm+9BTt2\nQLVqp7+mTj19GEwRkVAK5Iygg3Mu90CN75rZumAFKi1XXXUVK1eupFevXqxcuZJOnTqRmppK06ZN\ng3K8l1+G+fP9zUcLMngw9O9/+rxHHw1KJBGRgARyRnDczFrkTJjZBUBY9OvQtm1bVq1aRdu2bdm4\ncSPx8fGsWVP6A6vddBOsXw8rV8K//qWB7EUkvARSCH4HvGdmy81sOfAe/uakYaFhw4asWLGCq666\nir1793LllVeyZMmSUj1GtWrQuPHpLxGRcBFIX0PvAhcCdwC3A62cc+8FO1hpql69OkuWLOHGG28k\nMzOTfv368eKLZb4hlIhISATUrtI5l+WcW+ucW+ecC8sOE+Li4pg9ezYTJ07kxIkTjBw5kvvvv19d\nWYtI1IuqBvZmxiOPPMKMGTOIiYnhoYce4qabbuKomuyISBSLqkKQY9SoUSxatIhKlSoxa9Ys+vTp\nw6FDh7yOJSLiiYAKgZk1MLNOZtbVzK40s66Fb1W29enTh+XLl1OvXj2WLVtGly5d2Llzp9exRERC\nrtDnCMzsMfzDVP6PU+MXA7wfrFChcvnll5ORkUGvXr1Yt24d8fHxpKam0qZNG6+jiYiETCAPlA3A\n31LoSLDDeKF58+akp6fTv39//vWvf9G5c2fmz59PUlJSSHP8/vdQoQK0aAG33x7SQ4tIlAvk0tDn\nQFywg3ipdu3aLF26lEGDBnHw4EF69OjBX/7yl5Ad/7HH4KKL/GMZ/M2zzr1FJFoFckbwA/CRmb0L\n5JwVOOfcHcGLFXoVK1Zk3rx5jB8/nieeeIIbb7yRbdu2MWHCBMzO6r67VI0e7f/5z3/6n1AWEQml\nQArBouxXbp43vk9JSSExMZHExMRS22dMTAxTp06lcePGjB07lokTJ7Jt2zaeeeYZYmMD+acSESl7\nfD4fPp8v3+X5jlBWlp05Qlkw/PWvf+XGG2/kyJEj9O3blzlz5lC5cuUi7ychAUaOhEsugbFjYdy4\nszudy/HPf8KDD/p/ioiUtvxGKMv3z1wze9M5N8TMPs5jsTujR9KIM3jwYOrXr0///v1ZvHgxSUlJ\nvP3229SrV69I+2ndGp577tR0jRqlHFREpIQKGrP4fOfcLjNrmtdy59zW4MUqWCjOCHJ8+umn9OzZ\nk61bt9K8eXNSU1O58MILg3IsnRGISDAVecxi59yu7J9bgSzgJ0BbIMvLIhBqrVq1IiMjg/bt27Nl\nyxY6duxIenq617FEREpNoc1HzWwEsBoYCAwGPjCzm4MdrCypX78+Pp+P3r17s3//fq666irmz5/v\ndSwRkVIRyHMEvwcudc790jn3S/zDVt4V3FhlT5UqVVi4cCGjRo0iKyuLQYMG8fTTT3sdS0SkxAIp\nBPuA73NNf589L+rExsYyffp0HnnkEZxz3HnnnYwbN46TJ096HU1EpNgKajU0LvvtZvyXgxZkT/cH\nyvyYxcFiZkycOJFGjRoxfPhwpk6dyvbt25k9ezbnnHOO1/FERIqsoDOCqkAV/F1MLMD/EJkDFgJR\nPyrvsGHDSEtLo1q1arz55pt0796dAwcOeB1LRKTI9EBZCa1bt47evXuzc+dOWrduTWpqKk2bNi3W\nvtR8VESCqcjNRyUwl1xyCatWraJNmzZs3LiRhIQEPvzwQ69jiYgETIWgFDRs2JCVK1fSrVs39uzZ\nQ9euXUlNTS3RPkeOhIYN/a+ePUspqIhIHlQISkn16tVJTU3lxhtvJDMzk759+/Liiy8We38HDsD9\n98MLL8BXX5ViUBGRMwQyQlk9YCTQNNf6zjk3PIi5wlJcXByzZ8+mcePGTJ48mZEjR7Jt2zYmTZpU\nrK6sa9eGc88NQlARkVwCOSNYCFQDlgJ/z/WSPJgZjzzyCDNmzCAmJoaHHnqIm266iaNHjwbtmJ9/\nDo0b+1/t2gXtMCISoQLpZL+icy7qniQuqVGjRtGgQQOuvfZaZs2axc6dO3nrrbeoVq1aqR/r+HGI\njYW33oJrrin13YtIhAvkjOBtM+sT9CQRqE+fPixfvpx69eqxbNkyunTpws6dO4NyrPLl/TeWRUSK\nKpBC8FtgsZllmdl32a9DwQ4WKS6//HIyMjK48MILWbduHfHx8XzyySdexxIR+VGhhcA5V8U5F+Oc\nO8c5VzX7VfrXNyJY8+bNSU9Pp1OnTuzYsYPOnTvzTz01JiJlRL6FwMwuyv55WV6v0EWMDLVr12bp\n0qUMGjSIgwcP0qNHD15//XWvY4mIFHizeCz+ZqNTyXuw+qSgJIpgFStWZO7cuYwfP54nn3ySoUOH\nsm3bNu66666Ampdu2QJHjvjf160Ldeqcvc7x47Bhw6np8uWhRYtS+gVEJCLlWwiccyOzfyaGLE0U\nKFeuHE888QSNGzdm3Lhx3H333Wzbti17bIOCG3ENHAgHD8KhQzBmDNx775n7hvr1/esBHD0KzvkL\niIhIfvRksUfGjBnDvHnzqFChAtOmTWPgwIH88ENmodvNnw+jR+e9rFYt/9lAzmvp0lIOLSIRSYXA\nQ4MHD2bZsmXUqlWLxYsXM3ZsEkePqj8JEQmtQB4oCxkzawbcA1R3zg3xOk8odO7cmfT0dHr27Mmn\nn/6bc85J4LPPUoELC912xw5/l9W7dwc/p4hErkAGr+9sZlWy3w8zs6lm1iQYYZxzXzjnRgRj32VZ\nq1atsp81aE9W1hY6duzI/v3phW730UcwYQJ8/DH06xeCoCISkQK5NDQNyDSzn+BvSfQ5MDvQA5jZ\ny2a218w+PmN+TzPbaGabzCzqu7CoX78+Tzzho1at3uzfv58VK67igw/mF7rdpZfCq6/Cn/4UgpAi\nEpECKQTHs4cD+xnwZ+fcn/EPYxmomcBpPeqbWTng2ez5FwPX5zy3EM0qVqxCmzYLGTVqFCdPZvH4\n44OYM+dpr2OJSIQLpBB8Z2YTgRvx9ztUDigf6AGccyuAb86Y/VNgs3Nuq3PuGPAG0N/MapnZdKBd\ntJ4lxMTEMn36dC6++BGcczz++J1s3TqeDz44yeHDXqcTkUgUyM3i64DrgeHOuT1m1hgo6YWIBsD2\nXNM7gCuccweAWwLZQUpKyo/vExMTSUxMLGGkssPMaN16Ildd1Yhp04bzzTdT6NVrO82azaJixXO8\njiciYcLn8+Hz+Qpdr9BC4JzbbWZ/A3KeT90HLChRuryfVC6S3IUgUl155TD69TuPgQMH8s0382jT\nZjd16y4AankdTUTCwJl/JE+aNCnP9QJpNTQKeBOYkT2rIVD4XcyC7QQa5ZpuhP+sQM5w9dVXs3Ll\nSho0aMCKFSvo1KkT33671etYIhJBArlHcBvQGTgE4Jz7DKhXwuP+B2hpZk3NLA7/5adFJdxnxLrk\nkktYtWoVbdq0YePGjbzySgK7dn3odSwRiRCBFIIjzrkjORNmFksRLu2Y2RwgHbjQzLab2U3OuePA\nb4B/AP8D5jrnNhS0n2ixcSPcfDOsWXP6/IYNG7Jy5Uq6detGZuYeZszoyn33pXoTUkQiSiA3i5eb\n2T1AJTPrDtwKLA70AM656/OZnwoU+5ssJSUl4m4St2oFDz/sf9+xI1x2Rmff1atXJzU1lZ/97GZS\nU19j3bq+JCZOB6LuGTwRKYLCbhqb/xGB/GU3F70ZyBkN9x/Ai66wDYPIzLw8vOecc9x7771MnjwZ\ngPvuu49Jkyad1ZX1li1w9dXqfVRE/MwM59xZfd4HUggqA1nOuRPZ0+WACs45z1q1R3shyPH8888z\nevRoTp48yS9/+Uuef/554uLiflyuQiAiueVXCAK5R/AeUDHXdCVgWWkFk+IbNWoUCxcupFKlSsya\nNYvk5GQOHdJw0iJSNIEUggrOue9zJpxz3+EvBlIGJCcn4/P5qFevHkuXLqVr167s3Lnzx+XffAOT\nJ8Nrr3kYUkTKtEAKQaaZtc+ZMLPLgR+CFykwKSkpAT0xFw06dOiQ3Xvphaxdu5aEhATWr19PjRr+\nQWzWr4cXXvA6pYh4xefzFfgQbiD3CDrg7wsop9f784DrnHP/KaWMRaZ7BHnbv38//fr1Iz09nerV\nqzN//nySkpJYvhzuvx+WL/c6oYh4qdj3CJxz/wYuAkbj7weotZdFQPJXu3Ztli1bxqBBgzh48CA9\nevTg9ddf9zqWiJRxgQ5VeTlwCdAef5fRvwheJCmJihUrMnfuXH77299y7Ngxhg4dyl/+8gd0BiUi\n+Qnk0tBrQHPgI+BEznzn3O3BjVZgJl0aCsATTzzBuHHjcM5x3nmj2bbtaWJjy9TopCISQiV5jmAD\ncHFZ+uZVIQjcX//6V2644UaOHTtC3759mTNnDpUrVy5wm3fegZkz/e8HDID52V0MzpwJ56gXbJGw\nVZLnCD7Bf4O4TFGrocAMHjyYKVOWERtbi8WLF5OUlMRXX31V4DabN8O+fZCW5h8P+auv4G9/g+PH\nQxRaREpVabQa8gHtgNVATudzzjnn2XDpOiMomuXLYfz4T9m3rydbt26lefPmpKWl0bJlyzzXf+45\n+OQTWLgQhg+Hr7/2P4ewZw9UqRLi8CJSavI7IwjkgnFK9k8HWK73EkYqVWpFRkYGycnJrFmzhoSE\nBBYvXkxCQoLX0UTEY4E0H/UBW4Hy2e9XA/8NaioJivr16+Pz+ejduzf79++nW7duzJ9f0jGGRCTc\neTVCmXikSpUqLFy4kJEjR5KVlcWgQYN45plnvI4lIh7yaoQy8VBsbCwzZszg4YcfxjnHHXfcwfjx\n4zl58qTX0UTEA0EfoUzKntGjoXZtY8qUexgzZjYQy5QpU/jpT68nKysroH107w61asF9950+f9Ys\n//wrryz93CISHIEUgjNHKHuTIoxQFixqPlp8mZn+kdDKl4cuXYbRsmUq5ctXZc2aeVxzzTVkZh4o\ndB+HDkGXLvDDGd0PHjkCTZr4l4tI2VBY89FACsFdwNfAx8CvgSXAvaURriRyhqqU4qlSBWKy/+vX\nqnU1gwatpHr1BqxYsYIpUzpx6NDWQveRawyc0+jhZZGyJTExsfiFIPsy0P+cc8875wZnv15QI/7I\nU6fOJYwdu4o2bdqwd+9G5s9P4OjRD72OJSIhUGAhcM4dBz41syYhyiMeqlmzIStXruTCC7tx+PAe\n9u/vyubNqV7HEpEgC+TSUC1gvZm9Z2aLs1+Lgh1MvFG9enVuuy2Vli2H4lwmc+f2ZePGl7yOJSJB\nFMjV3Hs59URxDl0aimCxsXFcddWr7N7dmO+/f5Tly0dQvvw2nEvh7I+CiIS7QO4RPO+c853x0lhX\nEc7MqFZtMr16TccshmPHHmT06OEcO3bM62giUsoCuUewUfcIwt/77/tb+dxexFEk0tJ+TY8eC4FK\n/OUvr9CnTx8O5dE2tGVL//5nzy6NtCISSoFcGsq5R7AayMye52nvo3Cq+aiakBauSxd/+/7p02Hj\nxsC3++IL/8+YmGSqVvVRtWoyS5cupWvXrpj9HWjw47rHjkF8vLqqFimLfD5fgc9dBVII7it8ldAr\nqE2snC4mxv/XelHb9+d+TqBcuQ68+24Ggwf3Yu3atcTFJVC3birwf6cdR0TKnpw/midNmpTn8kK/\nGrJ7HBWhWbPmpKen069fP9LT0/H5OlGjxnwgyetoIlICgfQ++r2ZfZf9OmJmJ81MHQhEqdq1a7Ns\n2TJq1hzEsWMHeeutHrz++utexxKREghkPIIqzrmqzrmqQEVgIPBc0JNJmVWxYkVatJhLixa/5eTJ\nYwwdOpSDB/+AHjgXCU9FuqrrnDvpnFsA9AxSHgkTZuVo1+4JEhOnYmZ8++3dbN58G86d8DqaiBRR\nofcIzGxQrskYoD3wQz6rS5Rp334Mt97akGuvHcauXdP44YedNG06B6jkdTQRCVAg7Uj6cupJ4uP4\nh63sH6xAEn6GDBnCueeex/79/fjmm0UcPZpE06aL0fhFIuEhkFZDvwpBDgkDK1ZApUpQseKpeYcP\nw/LlAJ259NJ01q3rRWbmajZtSmDTpjRatmzpVVwRCVAgrYZmmVmNXNM1zezl4MYqnAamCa3OneEP\nf4Dx42Ho0FPzt2+Hnj2hVSuoWbM1bdpkUKlSe44e3UJCQgIZGRnehRYRoHQGpvmJc+7bnAnn3DfA\nZSWPVjIamCa00tL8f/nn1VK0cWP/svPOg7i4+rRq5aNatd7s37+fbt26MX/+/NAHFpEflWhgmmxm\nZrVyTdQCypU8mkSqcuWq0KzZQkaOHElWVhaDBg3imWee8TqWiOQjkEIwBcgws4fM7GEgA/hTcGNJ\nuDOLZcaMGTz88MM457jjjjsYP348J0+e9DqaiJwhkAfKZuN/iOwrYA8wIHueSIHMjHvuuYfZs2cT\nGxvLlClTuP7668nKyvI6mojkEsjN4nhgu3PuGefcs8AOM7si+NEkUgwbNozU1FSqVq3KvHnzuOaa\nazhw4IDXsUQkWyCXhqYD3+WazsyeJxKwq6++mpUrV9KgQQNWrFhBp06d2Lp1q9exRIQAu5hwuTqR\ncf4+BHSzWIrskksuISMjgzZt2rBx40YSEhL48MMPvY4lEvUCKQRfmNkdZlbezOLM7E5gS7CDSWRq\n1KgRK1asICkpiT179tC1a1fS0tK8jiUS1QIpBLcAnYCdwA4gHhgVzFAS2WrUqEFaWhpDhw4lMzOT\n5ORkXnrpJa9jiUStQFoN7XXOXeecq5f9ut4591UowknkiouL49VXX+Xuu+/mxIkTjBgxggceeEBd\nWYt4IJBWQ43MbL6ZfZ39esvMGoYinEQ2M2Py5MlMnz6dmJgYHnzwQYYPH86xY8e8jiYSVQK5NDQT\nWAScn/1TfhUEAAARjUlEQVRanD1PIsDkyadPDxwI06YFvv3q1UU/ps8H/fr5j9WvH+zY8WsWLlxI\npUqVeOWVV+jTpw+HDh3i6af9y9esKfoxCnLoEFx9NRT1atTcuf7tVKck0gTSDXVd51zuL/5XzGxM\nsAIFKqevIfU3VHyvvAJHjkDNmjBvnn/e22/DnDkQSKehP/sZDBsGlSsX7bg7d8I774Bz/i/VH36A\nhx5KxufzkZyczNKlS+natSutW/+dxYsbcOutRf7VCnTsGLz7LrRtW7TtNm3yb6erVxJufD5fgZ10\nBlII9pvZMOB1wICfA/tKJV0JFNSBkgSmR49T73MKAUDfvhAXV/j2LVrAFcV8tDA21v+FbHZqXocO\nHcjIyKBXr16sXbuWzZsTgFTg/4p3kELkPrZIJMv5o3nSpEl5Lg/k0tBw4Fr83UvsBoYAN5VaQpFc\nmjdvTnp6Oh07diQzczvQiXXrfF7HEologbQa2uqc6+ucq5v96u+c2xaKcBKdateuzbJly2jSZCBw\nkHvu6cGcOXO8jiUSsYo0eL1IqFSsWJGkpHnAnRw/fpQbbriBxx57TM1LRYJAhUDKrJiYcsCTjBo1\nFYAJEyZw2223ceLECW+DiUQYFQIp8wYMGMO8efOoUKEC06ZNY+DAgRw+fNjrWCIRI5AHyu7N9f6c\n4MYRyduQIUNYtmwZNWvWZNGiRSQlJfHVV3rAXaQ05FsIzGyCmXXE30ooR3rwI4nkrXPnzqSnp9O0\naVNWr15NQkICmzZt8jqWSNgr6IxgI/4i0MzMVprZC0AdM2sdmmgiZ2vdujUZGRm0b9+eLVu2kJCQ\nQEZGhtexRMJaQYXgW+Bu4HMgEXgacMBdZqb/88Qz9evXx+fz0bt3b/bv30+3bt2YP3++17FEwlZB\nhaAH8HfgAvwD2P8UOOycu8k5lxCKcCL5qVKlCgsXLmTkyJFkZWUxaNAgnnnmGa9jiYSlfAuBc+5u\n59xVwBfAq/i7o6hjZv8ys8WhCiiSn9jYWGbMmMHDDz+Mc4477riD3/3ud5w8edLraCJhJZDmo/9w\nzv3HOTcD2OGc64S/2wkRz5kZ99xzD7NmzSI2NpbHH3+cG264gaysLK+jiYSNQLqY+H2uyV9lz/s6\nWIFEiuMXv/gFS5YsoWrVqsydO5cePXpw4MABr2OJhIUiPVDmnFsbrCAiJdW9e3dWrFjB+eefz/vv\nv0/nzp358ssvvY4lUubpyeIoc/Kk/xVszsGZPUGcPOmfF8jx88uZs4/jx/0/z+x66Cc/+QmrVq2i\nTZs2bNiwgfj4eD788MMCj3XiRN77OvOYxVGSbUVCJWwLQUpKSoEDLcjZzOCFF/wDzwS7L/51684e\n0+Chh/zjENx+e+Hb33uvf+CcM02b5t9H+fL+n0eOnL1Oo0aNWLFiBUlJSezZs4euXbuSlpaW77Fy\n9pVfvbjzTv/y4khOhipViretSGnx+XwFjuES1oVAo5MVzejRp/76HTYseMe59NL8h3OsXTvw/Tzy\nyOmD5xRFjRo1SEtLY+jQoWRmZpKcnMxLBYxN2a5d/vsqyV/0O3aA7luL1xITEyOzEIgUJi4ujldf\nfZW7776bEydOMGLECB544AF1ZS1yhmKe8IqEBzNj8uTJNGnShFtvvZUHH3yQTZu2Ac8D5b2OJ1Im\n6IxAosKvf/1rFi5cSKVKlZgz5xWgD0eOHPI6lkiZoEIgUSM5ORmfz0edOnWBpbz5Zld27drldSwR\nz6kQSFTp0KEDqakZQEu+/not8fHxOLfe61ginlIhkKjTrNkFQDrnnZfA9u3bgU58953P41Qi3lEh\nkChVh8GD32XgwIHAQbZs6UFa2hyvQ4l4QoVAolZsbEXmzZsH3IlzR7n33ht47LHH1LxUoo4KgUS1\ncuXKYfYk558/FYAJEyZw2223cUL9QkgUUSEQAerVG8Mf/jCPChUqMG3aNAYOHMjhw4e9jiUSEioE\nItmuvnoIy5Yto2bNmixatIikpCS++uorr2OJBJ0KgUgunTt3Jj09naZNm7J69Wo6duzIt99u8jqW\nSFCpEIicoXXr1mRkZHDZZZfx+eefM39+R2CV17FEgkaFQCQP9evXZ/ny5fTq1YusrH1AEv/73wKv\nY4kEhQqBSD6qVKnCokWLuOiiEUAWb7wxEHjW61gipU6FQKQAsbGxXHnl88BD2c8X3M6ECb/jZCiG\neRMJERUCkUKYGXAvAwfOAmJ54onHueGGG8jSiDMSIVQIRAJ06aW/AJZQtWpV5s6dS48ePThw4IDX\nsURKTIVApEi68957Kzj//PN5//336dy5M19++aXXoURKRIVApIguueQnrFq1ijZt2rBhwwbi4+P5\nML+R70XCgAqBSDE0atSIFStWkJSUxJ49e+jatStpaWlexxIpFhUCkWKqUaMGaWlpDB06lMzMTJKT\nk3nppZe8jiVSZCoEIiUQFxfH7Nmzufvuuzlx4gQjRowgJSVFXVlLWClThcDMKpvZLDN73sxu8DqP\nSCBiYmKYPHky06ZNIyYmhkmTJnHzzTdz7Ngxr6OJBKRMFQJgIDDPOTcK6Od1GJGiuOWWW1iwYAGV\nKlVi5syZJCcnc+LEIa9jiRQq6IXAzF42s71m9vEZ83ua2UYz22Rmd2XPbgBsz36vkUEk7PTt2xef\nz0fdunV555132LKlK7DL61giBQrFGcFMoGfuGWZWDn+nLT2Bi4HrzewiYAfQKITZREpdhw4dyMjI\noGXLlmRlrQXiWb9+vdexRPIV9C9b59wK4JszZv8U2Oyc2+qcOwa8AfQH/gYMMrPngEXBziYSLBdc\ncAHp6elUqpQAbKdTp074fD6vY4nkyau/unNfAgL/mUAD59xh59xw59ytzrk5HmWLWlOmQCD3N7dv\nh0D/wB0yBOLj/a8XXjh92dGj/p+rVvmXjxuX9z5++9tT+/jTn05f1rUrTJ4M114LcXH+dRYu9C+7\n9VZo29b/vlUrMPMvv//+wLL/5jfQu/ep6enTT1+emQkXXAB//evZ2958M9x/fx2aNXsXGMjBgwfp\n0aMHc+bMYcAAeOQR2LzZv/3Gjf5tJk2CG2/0v9+0yb9s+/bT9ztkCNSoAaNHn33M++6Djh1PTc+a\n5f+dH3wwsN933Dj/v2VRffQRXHrpqelrr4W33gp8+6VLoW/foh+3ML/9LbzySunvNxLFenTcEret\nS0lJ+fF9YmIiiYmJJd1lVBs7FoYO9b+PLeBT0bAh5Pxhe/HFhe937Vp4+GFo3Ng/vWCB/wu1Sxd4\n4AFYvRruvBM++ADq1Dl7+yeegIMHT5+XkgL/+Id/2d69sGULvPeev4h98AHkjC65ejV88on//Wef\nnfrdNm8uPDfAkiXwxRenvnR79z69mJ044T/23r1nb/vyy/4v4TZtKgLzuPPOcTz11FPccMMNwDb2\n7fs9/foZW7ZAztDIb70FH38Mr73mn7dlC5zZr11O0Zk+HaZNO33Zv/8NGRmnpleu9P/8858DK35T\np0L9+jBxYuHr5vb11/5ikOPNN+Hcc2HQoMC2X7MG3n67aMcMxFNP+T9nv/pV6e87XPh8voDORL0q\nBDs5dS+A7Pc7irKD3IVASq5ZM/+rMBUr+v+qLop27eDCC/3v16zx/6xVy7+f774reNuLLso7K8AV\nV/jPTL7+Ou8zmQoVzp5XqxYcORJY7vLlT59u1OjseQWJi8t5V44nn3ySJk2aMHbsWGACn3++jRMn\nngbK/bi+WeD7zktMPuf3Jd2vhK8z/0ieNGlSnut5dWnoP0BLM2tqZnHAdeiegES4MWPGMHfuXCCO\n3bufY+zYQcBhr2OJhKT56BwgHbjQzLab2U3OuePAb4B/AP8D5jrnNgQ7i4jXrr32WmAZsbE1+ec/\nFwLd+Oabr72OJVEu6JeGnHPX5zM/FUgt7n5TUlJ0b0DCVBfatfsXu3b1YteuD/jVrxJ4771UoKXX\nwSRCFXavIGzb6ucUApFwVKnSRbz22irgMnbs+JyOHTuSmbnK61gSoRITEwu8rxq2hUAk3NWpUx9Y\nTqdOvdi3bx9btiQBC7yOJVFIhUDEU1WYOnURI0aMwLksYCDPPvus16EkyqgQiHgsNjaW559/nnPP\nfQhw3H777Uyd+jvgpNfRJEqEbSFISUnRI/sSMcyMc8+9F5hFbGwss2c/DtzAkSNZhW0qUiifzxeZ\n9wh0s1gi0y9YsmQJlStXBeYyfHgPDhw44HUoCXO6WSwSZrp3787LL68Azuff/36fzp078+WXX3od\nSyKYCoFIGdSq1U+AVbRs+X9s2LCB+Ph4/vvf/3odSyKUCoFImdWIOXNWkpSUxJ49e+jatSv+h/FF\nSpcKgUgZVq1aDVJTUxk6dCjff/890Ad42etYEmHCthCo1ZBEiwoVKjB79mzuvvtu/CO43gyk4FyJ\ne3OXKKFWQyIRICYmhsmTJwPT8P9vO4mbb76ZY4GMJCRRT62GRCLKLfi7oajEzJkzSU5O5tChQ16H\nkjCnQiASdvoCPurWrcs777xD165d2bVrl9ehJIypEIiEpQ5kZGTQsmVL1q5dS3x8PN99F+BA0iJn\nUCEQCVMXXHAB6enpJCQksH37dlat6gT4vI4lYUiFQCSM1alTh3fffZcBAwZw/PhBoAdvvPGG17Ek\nzIRtIVDzURG/ihUr8uabb9KkyR3AUa6//nr++Mc/qnmp/EjNR0WiQLly5bjooieBKQDcddddpKf/\nBv9zBxLt1HxUJEqYGTCWuXPnEhcXx4YNzwGDcO6w19GkjFMhEIkw1157LcuWLaNChZrAQr75phtf\nf/2117GkDFMhkCjl8zpAUHXp0oXk5H8BTTh27AMSEhLYtGmT17EiTqTcp1QhkCjl8zpA0NWseRGw\nitjYy/j888/p2LEjq1at8jpWRFEhkICF4sNSmsco7r727Al8u717A1u3sCyfflrYfgI7TllX/P++\n9alZczm9evVi3759JCUlsWDBgnzXPnKkuMc53Y4dxdtPUX/PQNYv6TqR8mVfEBWCEFAhOJsKQdGU\n5L9vTEwVFi1axIgRI8jKymLgwIE8++yzea579Gjxj5Pbzp3F248KgTcsHNsam1n4hRYRKQOcc3bm\nvLAsBCIiUnp0aUhEJMqpEIiIRDkVAhGRKKdCICIS5SKiEJhZZTObZWbPm9kNXueR8GNmzczsRTN7\n0+ssEr7MrH/299AbZtbd6zyBiohWQ2Y2DDjgnPu7mb3hnPu515kkPJnZm865IV7nkPBmZjWAx51z\nI7zOEogye0ZgZi+b2V4z+/iM+T3NbKOZbTKzu7JnNwC2Z79Xv7sCFPkzJJKnYn6O7gXyfmqvDCqz\nhQCYCfTMPcPMyuH/x+0JXAxcb2YXATuARtmrleXfSUKrKJ8hkfwE/Dkyv8eAVOfcR6GPWjxl9kvT\nObcC+OaM2T8FNjvntjrnjgFvAP2BvwGDzOw5YFFok0pZVZTPkJnVMrPpQDudJUhuRfwu+g1wFTDY\nzH4d2qTFF+t1gCLKfQkI/GcCVzj/yBvDvYkkYSa/z9AB4BZvIkkYyu9zdDvwjDeRiq/MnhHkI/zv\nbIvX9BmS0hBRn6NwKwQ7OXUvgOz3OzzKIuFJnyEpDRH1OQq3QvAfoKWZNTWzOOA6dE9AikafISkN\nEfU5KrOFwMzmAOnAhWa23cxucs4dx38z5h/A/4C5zrkNXuaUskufISkN0fA5iogHykREpPjK7BmB\niIiEhgqBiEiUUyEQEYlyKgQiIlFOhUBEJMqpEIiIRDkVAhGRKKdCICIS5VQIRESinAqBSCkwsw5m\nttbMKmSPof2JmV3sdS6RQKiLCZFSYmYPAecAFYHtzrnHPI4kEhAVApFSYmbl8fdK+QOQ4PQ/l4QJ\nXRoSKT11gMpAFfxnBSJhQWcEIqXEzBYBrwPNgfOyhy0UKfPCbcxikTLJzH4BHHHOvWFmMUC6mSU6\n53weRxMplM4IRESinO4RiIhEORUCEZEop0IgIhLlVAhERKKcCoGISJRTIRARiXIqBCIiUU6FQEQk\nyv0/LSySmB1xThcAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is not very successful, given that the data was generated with $\\alpha=2.5$. The problem appears to be that the number counts are small out in the tail so they drag the best fit to the wrong slope. Suppose we instead fit just the first $x_{\\rm max}$ points, for $x_{\\rm max}$ in (64,32,16,8,4,2):" ] }, { "cell_type": "code", "collapsed": false, "input": [ "figure(figsize=(6,5))\n", "yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')\n", "xlabel('x')\n", "ylabel('# occurrences in bin of size .1')\n", "ylim(.5,300),xlim(1,200)\n", "xscale('log')\n", "\n", "slope, intercept = linregress(log(bh[yh>0]),log(yh[yh>0]))[:2]\n", "xl=arange(1,200,.1)\n", "plot(xl,exp(intercept)*xl**(slope),'k-',linewidth=2,\n", " label='$\\\\alpha={:.2f}$ (all)'.format(-slope))\n", "\n", "for xmax in (64,32,16,8,4,2):\n", " retain=logical_and(yh>0,bh[:-1]<=xmax)\n", " slope,intercept = linregress(log(bh[retain]),log(yh[retain]))[:2]\n", " xl=arange(1,xmax+1,xmax-1)\n", " plot(xl,exp(intercept)*xl**(slope),linewidth=1,\n", " label='$\\\\alpha={:.2f}$ ({})'.format(-slope,xmax));\n", "legend();" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFHCAYAAACs30uOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdcVfUfx/HXl43iQgQUJymKK3CbI0wz0xT3KPXnTisz\n1Mw0FWdprhylVq7UHGlO1DTFwhGKmooDF05ERVFAQMb398eFG+MCFwQu4/t8PM5D7jnfc+73CNwP\nZ33fQkqJoiiKUngZGboDiqIoimGpQqAoilLIqUKgKIpSyKlCoCiKUsipQqAoilLIqUKgKIpSyJkY\nugNZIYRQ97wqiqJkgZRSpJyXb48IpJT5ZpoyZUq+eo+sbisz6+nbNqN2WV2eG9+TvPr9zY33yY2f\nIX3bv2qbrC7Li1Na8m0hyE/c3Nzy1XtkdVuZWU/fthm1e9Xl+UVu7Ud2vU9u/Azp2/5V2xSUn6H0\niPSqRF4lhJD5sd9K3uHp6Ymnp6ehu6Hkc/nt50gIgSxIp4YU5VUUhr/ylJxXUH6O1BGBoihKIZHW\nEUG+vGtIUZScIUSqzwgln8rMH8uqECiKkow62s7/MlvQ1TUCRVGUQk4VAkVRlEJOFQJFUZRCThUC\nRVGUQk4VAkVRlEJOFQJFUZQc9uWXX/Ldd9/p1bZy5cocOnRI+/Wff/6pXda4cWMuXryY7f1Tt48q\nilKonT17lnXr1jF37lydyzds2EBQUBC+vr506dKF3r17A7Br1y7u3r1LVFQUlSpVomvXrjrXf/To\nEb/88gvXr1/Xqz9Jb/0UQiR7PXbsWCZPnsxvv/2m7+7pRRUCRVEKrfnz5+Pj40OJEiV0Lr927Roh\nISGMGTOGx48fU61aNRo3boyJiQlXrlxh7NixAAwZMoS2bdtiZWWVahurV6+mQ4cOmJubv3J/O3bs\nyPDhwwkODsbOzu6Vt5dInRpSFKXQGj16NO7u7mku9/f3Z86cOQDY2NhQtWpV/Pz8ePz4MQcPHuTl\ny5cAFC1aFDMzM53b2LdvH2+++ab29TfffEPVqlUpXrw4tWrVYvv27Xr318LCgvr167N//36919GH\nOiJQFKXAuHHjBj/++GOay5s0aZLqgz+9J6nbt2/P3r17te2CgoKoWrUqLi4uxMfH07BhQ4YNG0bb\ntm3TLATnz5+nevXq2tdVq1bFx8cHe3t7Nm/eTN++fbl27Rr29vZ67aOzszP//vuvXm31lW8Lgaen\nJ25ubgVm9D9FyeuycxyirA5jcf/+fVatWoWLiwt//fUXI0aMwMbGhvDwcOzt7XF0dOTrr7/O1DbT\n2y9TU1Nq164NwJ49e2jQoAEuLi4AjB8/nq+//pqxY8eycOHCNLcRGhpKsWLFtK+7d++u/bpnz558\n/fXX+Pr60qlTJ736W6xYMYKCgvRqm8jb2xtvb+80l+frQqAoSuERERFBly5d8PLyonTp0tja2uLh\n4UG/fv147733srxdfYpSaGgoq1evZt26dQAEBATg7e3NgQMHOHjwIAMHDqROnTq88cYbqdYtVaoU\nYWFh2tdr165lwYIFBAYGAhAeHk5ISIje/X3+/DmlSpXSuz2g/aN56tSpOpfn20KgKEruMvRgdJs2\nbaJ+/fqULl0agDJlyuDv748QQntaJiunhjI60pFS8s033/DTTz9hZWVFYGAgu3btokePHgC0adOG\nNWvW4OPjo7MQ1K1blytXrlC/fn1u3brFsGHDOHToEE2bNkUIgaura6b+by9dukT//v31bq8PVQgU\nRckXYmJiqFq1qvZ1REQExsbGdOnSRTsvK6eGdH0IX79+HUdHR4QQLF68mB49ehAVFYWvry+RkZFU\nqVKFCxcuUKdOHQCio6Np0qSJzu23b9+eI0eO8P777xMREYEQAhsbG+Lj41m7di0XLlzQu69RUVGc\nPn2aX375JVP7mBF115CiKPlCnz59CAkJwcvLix07dhAUFISLiwsrV67kxYsXWdrmkiVLWLlyJd7e\n3kydOpXnz58D0KNHD86ePYuPjw8eHh40bNiQcuXK0bRpU6pWrUrXrl15+PAhs2bN4rvvvuPhw4e0\nbNlS53v0798fLy8voqKiqFmzJmPGjKFp06bY29tz4cIFmjdvrnd/d+3aRatWrfS+sKwvlVCmKIpW\nQoKVobtR4EycOBFbW1tGjRr1Sttp0qQJK1eupGbNmum2S+v7mFZCmSoEiqJoqUJQMGS2EKhTQ4qi\nKIWcKgSKoiiFnCoEiqIohZwqBIqiKIVcvi0E8fGxhu6CoihKgZBvC0Fw8BpDd0FRFKVAyLeFIDDQ\nk7i4SEN3Q1EUJd/Lt4Xgr78aMGTIUqyt0U7ff2/oXimKouQ/+faBsnv3LnD9eitq1AjA2LgkEydC\ntWowerShe6co+Zd6oKxgKDQPlJUrV4syZToQFjYXa2uwtDR0jxRFUXTLTHh9RnIiwD7fFgKAypWn\ncv/+D0RHPzB0VxRFyafOnj2rzR7WZceOHaxfv55p06bxvY7zz76+vumOeJoYXj98+HDtvI0bN+Ls\n7IyVlZU2sSylq1evYmFhQb9+/ZLNTwywz075dhjq48ePJ4zg9z9u3ZoBLDF0lxRFyWcyCq8PDQ2l\nV69ehIaGYm5ujo2NDR06dKBSpUoAxMfHM3nyZJ05BIlShtcfOHCA8ePHs3nzZho1akRQUJDO0zgf\nf/wxjRo1SpWXkBMB9vn2iGDkyJHExcVRseKXPHz4K0WKXDd0lxRFyWcyCq8vWbIkfn5+WFhYIIQg\nNjY22Yf2li1baNOmTbrXVVKG10+ZMoUpU6bQqFEjAMqWLUu5cuWSrbNx40ZKlSpF69atU207JwLs\n8+0RgaWlJT///DPDhg2jfPlR1KgxmYcP1xu6W4qiGFB2h9cD1KpVCwAfHx/c3NyoXLkyoDnlY2xs\nTJkyZYiIiEhz/aTh9XFxcfj5+eHu7k61atWIioqic+fOfPvtt1hYWACaKMopU6Zw+PBhVqxYoXOb\n2R1gn28LweLFi3nnnXfo3r075ct7YGNTjbCws4CLobumKEoOye3w+kTbtm1jy5YtzJs3L9m8YcOG\nsXbt2nTXTRpeHxwcTExMDFu3bsXHxwcTExPc3d2ZMWMGM2bMAGDSpEkMGTKEcuXKpdm3rATYpyff\nFgIXFxe6d+/OpEmTWLp0KVevTqRixYnAHkN3TVEKLDE14w9Nfcgpmb9F1ZDh9V27dqVt27a4urpy\n4MABHjx4QOPGjfW63TZpeL1lwu2NI0eO1J7fHz16tLYQnD17lj///JMzZ86k27esBNinJ98Wgmcn\nnjF9+nScnZ0ZOnQot24Nw8lpPqGhf1GypO7IOEVRXk1WPsCziyHC6/fs2cOsWbM4evQoVlZW2Nra\n8ttvv2Fubs6LFy/Yv38/R48eJTIykp07d9KpU6dU20gaXl+qVCnKly+f5vt5e3sTGBhIxYoVAQgP\nDycuLo5Lly5x6tQpbbvsDrDPt4Ug4MMA6p+qz/Tp0xk5ciQNGvxFUNA0btwYj6vrUb0O9xRFyT8M\nEV5vbGyMm5ubtt2dO3eoW7cubdu21bb19PRECKGzCEDy8HqAgQMHsnjxYtq1a4eJiQkLFiygY8eO\nAHz44YfadlJK5s6dS2BgIMuWLdNuLycC7PPcXUNCCHchxAohxEYhxNtptTOzNePud3cZPHgwL168\n4PLlDTx58j5xcWGEhOzKzS4ripILDBFe365dOxwcHFi8eDGff/45EydOTFYENm/ezM6dO9m5cydb\ntmzR+R5Jw+tBcw2gYcOGODk5UbNmTerXr8/EiRMBzakjW1tbbG1tsbOzw8rKCktLS+1REORMgH2e\nHWJCCFESmCulHKJjmYy4GsHpJqdp4NeAM/fP8Pbb3Zkw4TLDhh3hxo3xNGz4L0IYG6DnipJ/qSEm\nckZ2hdeDfgH2eTK8XgixEugAPJRS1kkyvx2wEDAGfpJSzk6ybC6wTkp5Vsf2pJSSwOmBhJ0Mo/aO\n2tSuPRBbW1sOHZrNmTMtKFduKPb2/8vxfVOUgkQVgoIhr441tApol6JDxmgeB24H1AT6CCGchcZs\nYK+uIpBUxXEVeRHwgsfbH9O8+Tf4+q7in3+uUKbMN1y/PoXHj6N58iSndklRFKVgyJVCIKX8G3ia\nYnYj4JqUMlBKGQNsBNyBT4DWQHchxIfpbdfI3AinZU5c+/QadsVskHICLVp8St26zTh+vDYeHstI\nePZDURRFSYMh7xpyAO4keX0XaCylHAkszmhlT09P7ddValWhf1wZJj37BBeXn5gxYztvvz2Ls2fb\n8scfg4Bi2dx1RVGUvM/b2xtvb+8M2+XaxWIhRGVgV+I1AiFEN6CdlHJowuu+/FcIMtqWTNrvl49f\ncrLWSerurcvJ0JMMHjyYixcvcvXqMJYufY3lyz1zYI8UpeBR1wgKhrx6jUCXe0CFJK8roDkqyDQz\nGzMcv3Yk4MMAWr3ZigYNGjB79mzKlp1Khw6LefnyYbZ0WFEUpSAyZCE4BVQTQlQWQpgBvYCdWd2Y\n/UB7jCyNuPfDPebNm8fixYsJChL89df73Lo1K9s6rSiKUtDkSiEQQvwKHAOchBB3hBADpZSxaC4M\n7wcuApuklJf03aanp2eyc19CCJyWOXFr6i3sTOzw8PBg4sTRbN78FcHBvxAZGZit+6QoipJfeHt7\nJ7uumlKefaAsPSmvESR1Y+INIq9G8tra13B2rkVw8Pf4+x8lKuoWzs5rcrmnipK/qGsEBUN+ukaQ\nIyp9VYkwvzAiDkfwzTcLiYr6FDu7kTx5spfw8AuG7p6iKEqeU+AKgbGlMU7fO3H146u0bfkuRkZV\nWbp0NRUrjufmzYmG7p6iKIWQCq83AOt3rCnWqBgPvr1NXNxCxo2bjZNTZy5ePMvrrx/D1BTU0a+i\nKJBxeH2ilCH1u3bt4ocffmDBggVs27YtzfVShtf37duXsmXLUrx4cRwdHZk5c6a27cuXLxk8eDCV\nK1emePHiuLq6sm/fvmTby4nweqSU+W7SdDt9UfejpI+Nj3x6OkyOG/el7NPnA3nnzkrp59dCQryM\nj89wE4pS6Ojzu1WQzJs3T3bp0kUOGDAg3XZxcXHynXfekVOnTpVSSnn79m357bffapcPHjxYhoWF\n6Vx3zpw5ctiwYdrXFy5ckJGRkVJKKS9fvizt7Ozk3r17pZRSRkRESE9PT3nr1i0ppZS7d++WxYoV\nk4GBgdr1IyMjpbW1tXzw4EGa/U3r+5gwP9Vnar49Ikh511BK5mXNqTy1MjdGBvDVxC/5++8jXL9e\nhdjYxzRuvDfX+qkoSt6VUXh9opQh9Y8ePeLgwYO8fPkSgKJFi2rDcVJKGV5fq1YtbT4xgImJCba2\ntgAUKVKEKVOmaINpOnToQJUqVTh9+rS2fVbC6zO6ayjfBtOkt1OJyn1YjgdrHhC2OYy5c+cyatRn\n7Ns3jaFDv0TKdgiRb+ugoig65ER4va6Q+nr16hEfH0/Dhg0ZNmwYbdu2TbMQJA2vT/TRRx+xZs0a\noqOjWbJkCfXq1dO5bnBwMAEBAdSqVSvZ/MyG17u5ueHm5sbUqVN1Ls+3hUAfwljgtNyJc23P0fl8\nZ5YtW8bWrcEIYcnDhxuxt3/f0F1UFCUTDBFen1ZI/fjx4/n6668ZO3YsCxcuTHP9pOH1ib7//nuW\nLl3KkSNH6N69O/Xq1aNRo0bJ2sTExPDBBx8wYMAAnJycki1T4fWZVMylGHb97Ljx+Q0WLVpE69at\nsbX9kXr1RmNr2x0jIzOCgsDXN/W6HTqASYH/H1KUTMiuCNgs3K1hiPD6EydO6AypDwgIwNvbmwMH\nDnDw4EEGDhxInTp1eOONN1JtI2l4fVJCCNzc3OjRowe//vprskIQHx9Pv379sLCwYMmSJanWVeH1\nWVB5amVO1jxJjQE16NOnD4sWeWFpWY2goB9xcPiYf/6B4cOhceP/1tm9G549Aysrg3VbUfIeA95u\nZ4jw+pMnT6YKqd+xYwfXrl2jR48eALRp04Y1a9bg4+OjsxAkDa/XJSYmJlkUpZSSwYMH8+jRI7y8\nvDA2Tp20qMLrE3h89RXubdpog6XTY2JlQrXF1QgYEcCUv6awaFEtnj6dS0TEWOzs/gdY0aQJbN/+\n3zqqAChK3mKI8PqRI/8bDDkxpN7d3Z1t27Zx4cIF6tTRBC5GR0fTpEkTndtPGl7/6NEj/vzzTzp2\n7IiFhQUHDx5ky5YtHDx4UNt+xIgRXL58mYMHD2Jubp5qe1kJr89wOGpdtxLl9QmQDU+dkpGxsWne\nPqXLOfdz8qbnTQk/yyZNmsjz53vJwMAZ8vffpXR3T962aFEp07gbTFEKLPLw7aPPnj2TEyZMkHv2\n7JHbt2+XBw4ckL1795Y///yzjIiIyNI2Fy9eLFu0aCErV64sPT095bNnz6SUUrq6usrTp09r223a\ntEm6urrKevXqyS1btkgppVy4cKGcOXOmXLhwoVy9enWa7/H48WNZvnx5GRkZKR89eiTffPNNWbJk\nSVmiRAnZsGFDuWPHDm3bwMBAKYSQlpaW0srKSjtt2LBB22bz5s2yW7du6e5XWt9H0rh9NEtjDQkh\nakgpL2d6xWwihJDdL1yglIkJK1JcjU9P1J0oTrmeoneICw6NWjNoUFdq1fqWkJArrFpVOtURwYMH\n6shAKVzUWEM5o0CG1wshbkspK2Z6xWwihJDPY2JofPo0o8uXZ0i5cnqve2fBHbaMDqHshpd8+qk7\nv/76Lk+fWrN+/VxVCJRCTxWCgiGzhSDNawRCiPTiIrPvcnUWFTMx4ffatWlx5gyvW1nRsHhxvdZz\nGOlAmQnB7B9bkdjYDvTvL1ixYhVOTqNInpOjKIpSOKT3RNUA4ALghyZEJnHyA17meM/0UL1IEZY7\nOdHd359HL/XrkpGJEe5HnBgWfwP/f6YRE7MdE5OuDBqk+0ELRVGUgi7NU0NCiMPAV1LKozqWBUop\nK+dw39KUMo9gwo0b/PP8Ofvr1sXESL+nhQM+CUC+lPxR5w+2bdvC9OkXcXX1oWjRGoA6NaQUTurU\nUMGQnXkE3YAzuhYYsggkSjrW0PQqVTAWgok3b+q9vuNMR0L2hPB+3fcJCQnl/Pl2aphqRVEKpEKT\nUPb45Usa+Pkxr2pVupUpo9d2Hm5+yK3ptwhfGM7/BvZn1ao4GjbcTvHijdQRgVIoqSOCgiFXEsqE\nEHnuhLqNmRlba9dmeEAAlxIGhspImR5lMC9vzmunXqN58xbs3OnKjRvj1S+CoiiFSlaH3zyVrb3I\nJvWLFWOOoyNdLlzgeWxshu2FEFRbWo3b395m2shprFt3gmvXbvL06cEM11UURSkoCsypoaSGX7nC\nw5gYttaqleHIggC3vr7FM59neLXw4tChLcyaBa1anSQoyEidGlIKFXVqqGAo9OH1AN9Vq8b96Ghm\n376tV/sKYyoQFRjFB5U+IDAwDB+f5zRr9lsO91JRFCVvKJCFwNzIiN9q1WLRvXscfPIkw/ZGZkY4\nLXfizud3WPDNAhYtekGvXhOIj4/Jhd4qilLQqfB6AylvYcEGZ2f6XrrEraioDNuXbF4S63bWOHk7\nUadOA377DUJCVuVCTxVFMaSMwuszCqlPGWqfUsrw+iVLltCgQQMsLCwYOHCgznU2btyIs7MzVlZW\nVK1aFR8fH+2ynAivz7AQCCHshRA/CyH2JbyuKYQYnK29yIKMMosB3EqVYlzFinS7cIGouLgMt/na\n7Nd4uPkh0wZOY//+x5w5M5m4uBfZ1GNFUfKa+fPnM23aNEJCQnQuv3PnDleuXGHEiBF4eHjg5eVF\neHi4dnl8fDyTJ08mJibtswerV6+mQ4cO2iGlHRwcmDRpEoMGDdLZ/sCBA4wfP541a9YQHh7O33//\njaOjo3Z5x44dOXz4MMHBwXrvZ0bPEehzRLAa+ANIHNntKuChdw9yiKenp15ZBB7ly1PV0pKPr17N\n8CKYaWlTXpvzGjHTYzAWI1mxwpR799IbcklRlPwso/D6x48fpxtSnzLUXpeU4fVdunTB3d09WRhN\nUlOmTGHKlCnaxLKyZctSLsnAmlkJr3dzc3vl8HobKeUmIcR4AClljBAi43sz8wghBD9Vr06T06f5\nMSiIYRmMVGrXz44Hqx/gLvpz4uJKduyYxfDhwzA1Nfg4e4qiZCC7w+tdXV3TDKnXFWqvi67w+rTe\nNy4uDj8/P9zd3alWrRpRUVF07tyZb7/9FgsLC227zIbXZ0SfQhAuhNCWLiFEE+BZtvUgF1iZmLCt\ndm2aJ4xU2jidkUqFEDj94MQHNU7TfuxcZi3+mLZtv6Z69Tm52GNFUXQxRHh9WiH1aYXap6QrvD6t\n9w0ODiYmJoatW7fi4+ODiYkJ7u7uzJgxgxkzZmjbGSK8fgywC3AUQhwDygDds60HWSVlpoK0nYoU\n4afq1enh78+p+vWxTXJ4l1KR6kXYY+rAiIMlqFChNt9/v5Q5cz7D3Fz/3ANFKYhEBtfl9CX1OK2b\nkiHC69MKqTc2NtYZaq9LWuH1utaztLQEYOTIkdjZ2QGa01cpC0Guh9dLKf2EEC2BGoAArkgpDT8M\n9aefwqJFmSoGnWxs8H3+nF4XL3Igg5FKt5hWZODFU0wf8w3vTWxDz57jadYs/cqvKAVdVj7As4sh\nwut37dqlM6Te0tIyVaj9zp076dSpU6ptpBVer+t9S5UqRfny5dPsT6JcD68XQtwAvpVS/pBk3m4p\nZdZLcHY4cQLGjIF58zJVDKZWqUKHc+cYf+MGc5MEYacUI4yZE1ONvqOuUMR6AJ6eK1mxYhJVqlTj\nxQto3Dh5+w8+gPHjs7oziqJkxBDh9VWqVNEZUt+yZUtt28RQe11FAJKH14PmOkBMTAyxsbHExcUR\nHR2NiYkJxsbGAAwcOJDFixfTrl07TExMWLBgAR07dtRuLyvh9RnSFWScdAKuAJuAVYB5wrwzGa2X\nkxMg5ZMnUrq6SvnFF1LGx+sMak5LyMuXssrx43JTcHCabS5ckPLcOSmPdvCXB98/I4sUKSZHjmwj\n69eX0sVFSgsLzfJz56QcN07KTz/NVBcUJU9ChddLKZOH16cXUp801H7z5s063yNpeL2UUk6ZMkUK\nIZJNU6dO1baPiYmRH330kSxZsqS0t7eXo0aNktHR0drlBgmvF0KckVK6CiHGocko6Alsl1K6Zl85\nyhztWEMhIdCqFXTpAlMzNyDqmbAw2p47h7eLC7WKFk2z3cvgl5ysc5KDvU7z64ExrFp1DFPTBhgZ\nQb16mjbffQc3bmj+VZT8TI01lDPyfXh9YiFI+LoNsBSwllLqN+h/Dkg26NzDh5pi0KcPfPVVpraz\n5sEDZt26hW/9+pQwSfss2f0V97m/8j4j4nrRvr0FU6f6J1uuCoFSUKhCUDDkxKBz2meZpZQHgbaA\nwZ+y0j5ZbGsLf/4J69bBnMzd4vk/e3valCrF/y5dIj6dH/6yQ8piZGzE9LbfsXjxFW7e3PGKvVcU\nRck9WU4oE0I4SykvCSHqAykbCSmlX7b1MpN0DkN97x64ucHHH8Nnn+m9rZfx8bidPct7pUszoVKl\nNNuFnw/n39b/srz1ZKLjL7Bx4wPtVX91RKAUFOqIoGDI7BFBencNjQaGAvNIXQgAWmW1kznCwUFz\nZODmBqammoKgBzMjI7bUqkVDPz8aFCtGW2trne2s6lhhP9CeD69PpeOhVnh7L6BVq9HZuAOKoiiG\nkW+DaR6GP6RMUR2XKW7e1BSDr76CoUP13uZfoaH09PfnRL16VE54qCOluIg4TtY+yZ+dfua3v9Zz\n6tQzjI1N1RGBUmCoI4KCIduvEQghegghiid8PUkIsU0IUS9bevsKmq9qzu1nOoJnqlTRHBlMmwar\nV+u9vZYlSzK+YkW6+vsTmcZIpcZFjam2pBpuewYQGWnEsmX6HXUoiqLkZXpdLJZSPhdCNAdaAyuB\nZTnbrYwNrz+c5iubc+nRpdQLq1aFgwdh4kRYv17vbY4qX57qRYowIiAgzb+KSncoTQmXEkxu/SVT\np64kNPRRVndBURQlT9CnECT+efwe8KOUcjdgmnNd0o9HUw+mt5pOqzWtOHX/VOoG1avDgQMwdixs\n3qzXNhNHKvULD2fZ/ftptqv2XTUqbXajUX17vvyyV1Z3QVEUJU/QpxDcE0KsAHoBe4QQFnqul6O+\n/vpr+r/enxUdV9B+fXsO3TyUulHNmrB/v2Zcot9/12u7RY2N+b1WLaYEBnL8me5BVs0dzKk0qRKf\nGE/h11+9efRIRyFSFEXJJ/T5QO8J7AfaSilDgVLA5znaKz1MmDCBjz76iPavtWdLjy30/q03v1/S\n8WFfty54ecHw4bBrl17brlqkCCurV6fnxYsEv9Q9vp7Dxw7YPHBlYLe67N7dS11gUxQl38qwEEgp\nI6SUW6WUVxNeB0kp/8j5rqXP3NycZcuW0aVLFxqUacC+vvv4yOsjVp5ZmbpxvXqwezcMHgz79um1\n/fdsbBhkb09Pf39i4uNTLRfGgurLq+N+4kuePg3kzh01MqmiKLplJrw+J8LpM5Jvbx/18fGhU6dO\nPHnyhIYNG7Jr1y5CjUN5Z907fNLoE8a+oSOM+tgx6NwZNmyANm0yfJ94KXnv/HlqFCnC/DRGKr06\n6ipbnn/Egt3/cPt2sHY8cUXJjwrb7aMbNmwgKCgIX19funTpQu/evVO12bFjB+Hh4Vy/fh0bGxs+\n+uijdOen9OjRI1xdXbl+/Trm5ubcvXuXESNGcOzYMczMzOjevTsLFy7Ujj66ZcsWNm3axG+//Zbl\n/crs7aPpjfBpkdYyQ08kjKx3+fJlWaVKFQnIKlWqyCtXrsjbobdljSU15BcHvpDxukYlPXJEShsb\nKb29Uy/TIeTlS+l4/Lj89cEDnctjnsXI/RX3yObNTOVXX30m4+OlvHkz+ZQw6KCi5Hnk4dFHs9vV\nq1flokWLpJRSPnr0SJYsWVLeuHEjWZunT59Kc3NzGRkZKePj46W1tbUMDAxMc74uc+bMkcOGDdO+\n7tKlixwwYICMjo6WDx48kHXq1NH2Q0opIyMjpbW1tXyQxmeOPtL6PpLG6KPpnRo6llBB1mW5LOUg\nT09PgoJDJNMcAAAgAElEQVSCOH78OPXr1+fmzZu88cYb3PG/w98D/+bQzUN8uPtD4uJTPBPQsiVs\n2gTdu8PRoxm+j7WpKdtq12bktWucDw9PtdykuAn3WzdkSKVuLF78PdevB1KliuaZNjc3zc1LZ85k\nzz4ripJ9/P39mZMwPpmNjQ1Vq1bFzy/5yDklS5bEz88PCwsLhBDExsYipUxzvi4pw+v9/f3p1asX\nZmZm2NnZ0a5dO/z9/xvIMivh9Bl5lbGG/IFZwHRgLJp0skRSSrkt23qZSSnHGgoPD6dXr154eXlh\nYWHBhg0baNO+DV02daGUZSnWdVmHuYl58o388Qf07au5gJwyZUaHdQ8eMPXWLU7Wq0dJ0+R3z363\nUGK//AQHGrXh4dOmeHkdJDZWs6xpU5g/X/OvouR1+f3UUGYSymJiYrhy5Qq1a9dGSkmFChXYvXs3\nLi4uOtf18fHh22+/ZceOHXrNT2Rra8vevXu1CWWffvopoaGhLF++nCdPntCuXTtmzJiRLDlt1KhR\nmJiYMG/evEztf6LsPDXUAs2DYyFoQmmSTWmtlxsTOg57YmJi5LBhwyQghRDyu+++k1ExUbLrpq7y\n7bVvy7DosNTHSbt3S2lrK+WpU2keYiX1SUCA7HjunIxLccpp4UIpxw94If94/2PpUM5SGhn9oV3W\npImUx47ptXlFMThdv1t5yb179+SMGTPk7t275bhx4+TNmzdlWFiYDAoKeqXt7tq1S7q7u6e5fOvW\nrbJ3797y6tWres1PytTUVF65ckX7OiQkRLq6ukoTExMphJADBw5Mtc7EiRPloEGDsrAnGml9H0nj\n1JA+H7pDMmqT21NaOxkfHy9nzpwp0QySJ0ePHi2jY6Ll4B2DZaMfG8nHEY9Tr7R9u5R2dlKePZv2\n/2qC6Lg42czPT067eTPZ/IULNQllgXMC5IzxpaS5eUVtopAqBEp+klEhOMzhbJmyIjw8XDZq1Eg+\nfqz5Pfb19ZWdO3eWW7duTZbglVlPnz6V3bp1k2FhOv5YTCIsLExWrVpV3kzx+5/W/ES2trbyVMIf\nm/Hx8bJBgwZy1qxZ8uXLlzIkJES6u7vLcePGJVtn5MiRcsyYMVnep8wWggwzi4G1QohRQGJIpzew\nTEoZo8e6uUoIwYQJE6hQoQKDBg1i/vz53LlzhzVr1jDFZwotV7fkj75/4FDc4b+V3N0hJgbatdM8\niVy7dprbTxyptEHCSKXvJoRoJ6rw2Wu4fTiSzY7fsWjRIsaOTX7n0tq1cEnHiBgAU6aAhUWWd11R\ncoWbdDPYe+dEeL2Ukm+++YaffvoJKysrbt26RaUkw9Hv2bOHWbNmcfToUaysrLC1teW3337D2dlZ\n5/yUv/OQPLz+8ePH+Pn5cejQIUxNTbG2tmbAgAFMmjSJ2bNna9fJ7nD6jOhTCH5IaLcUzXWCfgnz\nhuRgv15Jv379KFeuHF27dmXLli0EBQWxY8cObIrY0HxVc/7o+wfVSlf7b4Xu3TXFoG1bOHQIatRI\nc9tlzc3ZVLMm3fz9OV6vHo5Jbhc1MjWixv9G4VF1JWNmTeODDz4AymqXb90KxYtrHnhOytMTvvxS\nFQJFSU9OhNcvXryYHj16EBUVha+vL5GRkVSqVEkbXm9sbIybmxugKRp37tyhbt26xMfH65yvS9Lw\nehsbG8qWLcsPP/zAmDFjCAsLY82aNbz++uva9jkSTp8RXYcJMvlpmHP6zMvNCT3PY547d046ODhI\nQNaoUUPevHlT/uj3oyw7t6w8ff906hXWrJHSwUHKgIAMt73ozh35uq+vjIiNlQsXSmllJWWFClKW\nLy/l7De+lx90LSX79eub7NRQp06aM1EpFSsmZUJmtqIYlL6/W4aQ3eH1f//9tzQyMtIGyBsZGcm7\nd+9KKZOH1y9dulQuWrRIjhkzRi5btky7flrzU0oZXn/ixAnZvHlzWbJkSWljYyN79eolHz58qG2v\nTzh9RtL6PvIK1whOA1WTvH4NOJ3Rejk5ZeaH9c6dO7JOnToSkHZ2dvLUqVPyN//fZJk5ZeSRwCOp\nV/jpJ80n+vXr6W43Pj5efuDvL/tdvChDQ+PlrVtSO13/N1p6Lawp7cuUlLVq+ahCoOQbebkQ5GcT\nJkyQCxcu1Ktt48aNpb+//yu9X04UgtbAbeBIwnQLeCuj9XJyyuwPa2hoqGzdurUEZNGiReWePXvk\ngesHZJk5ZeTOyztTr/DDD1JWrixlGg+IJIqIjZV1fX3lkoS/IpK6sXm7nPRZSVm0aF3599+xUkpV\nCJS8TxWCgiGzhUCfsYb+BJyAT4GRQHUppY6hPvOuEiVK4OXlRd++fYmIiKBTp04EHgpk9/u7Gbpr\nKL/8m+Jc3PDh4OEBb70Fd++mud0ixsZsq12baYGBHE0xUmnl7p1o61ybKtZR3LhxISd2S1EUJVvo\nNZy0lDJKSvmvlPKclDIqpzuVE8zMzFi7di0TJkwgLi6OoUOHsnv5bv7s/ycTDk3guxMpBoT69FP4\n6CNNMQgKSnO7r1lasqpGDXr5+/MgOlo7XwhBjeZzWbjgEa8VK5dTu6UoivLKDJ4rkJuEEMycOZPl\ny5djZGTE9OnT+Xbctxzqe4ilJ5cy+fDkxFNPGmPGwMCBmmIQHJzmdtuXLs3QcuXocfFispFKbWo2\n5uHllsSe9szBvVIURXk1haoQJBo2bBg7d+6kSJEirFmzho/e/4i9PfayO2A3I/eOJF4mGXb6yy+h\nd2/NaKWPH6e5zUmVKlHC2Jix168nm7/BZy6y0XqCdqTxAIGiKIqB6VUIhBAOQohmQoiWQog3hRAt\nM14rb+vQoQNHjhzB1taWgwcP0vWdrqxvu57zD8/Td1tfXsYlCaSZPBk6ddIUgydPdG7PSAh+cXZm\nT0gIG5IcPTwOrYrRy55cOz4V08QBiBRFUfKQDAuBEGI2cBSYiCaZbCx5IKEsOzRo0IDjx4/j5OTE\nuXPnaNuyLfNc5xERE0HnjZ15EfNC01AImDED3n5b89BZaKjO7ZVKGKl01LVrnEsyUqlJ2enEu+3j\nzUcZj3aqKIqS2/Q5IuiC5k6h9lLKjolTTncstzg6OnLs2DGaNWvG3bt3afNmG0bajqRM0TK8/cvb\nPI18qmkoBMyZA82aaYajeP5c5/bqWlnxXdWqdLlwgacxMQmr2uFQ/iNqtJmD0Y0wneuNGwejRsHi\nxTmym4qiKGnSpxBcB8xyuiOGVLp0aQ4cOEC3bt149uwZ7du1p014GxqVa8Sbq98kKCzhriEhYOFC\nTfRl+/agI58A4H07OzqWLk3fS5eQQnPxuXKNL4hv6ovl9r3IuOTDw86eDc7OEBcH2ww2uLeiKIWV\nPoUgEjgrhFghhFicMC3K6Y7lNktLSzZv3oyHhwcxMTH079efMqfL0KNmD5qvas6Npzc0DYWAJUs0\n4xG99x68eKFze9++9hphcXHce+sWACYmJTjhNx6jrj9wf9n9ZG1HjNAcDXTrlqO7qCiKopM+hWAn\nmnCaY4BfksmgPD098fb2ztZtGhkZMX/+fBYsWIAQgokTJ3J/4308GnvQclVLzgefT2wIK1ZApUqa\n0UsjI1Nty9TIiM01a/Kw4X2OxoUA8O+5j8H5Gjc27Cb6fnSqdRRFKZgMHV6fUUKZwYaJeJWJXHgM\nfsuWLdLc3FwCsmPHjnLVqVXS9ltbefT20f8axcZK2aePlO3aSRkVpXM7tXqGyuIHfOTv/7yQLVpI\nuWfPz/LY1kbyfM/zqdoeOiSlm1tO7ZGiZCw3frfykvXr18u5c+fKnj17yl9//TXV8ri4OFm8eHFZ\nsmRJ7dSzZ0+91k308OFD6eDgIKNSfEYEBARIc3Nz2bdv32Tz89Sgc8CWhH/P65jyxeijr+rvv/+W\n1tbWEpANGzaUv578VZaZU0buvbr3v0YxMVJ26yZlx45S6gjHGDBAygoj70rL9b7StUmsPHw4Rp44\nXkP6dJovH+9NHpSjCoFiaIWpEOgTXn/jxg25fv16efPmTRkYGCgXLlwoL168qNe6iVKG1yd6++23\nZYsWLWS/fv2Szc9r4fWjEv7tqGPqlOljk3yoefPmHDt2jMqVK3Py5Ekm9prIojcW8b/t/2PjhY2a\nRiYm8OuvmtNFvXtrcg2SWLUKbn1Xju6uVjj/fIU33zTG8bWZGI9aScDHV4h7EWeAPVMURZ/wenNz\nczp37kzlypUpXrw4pqamODs767VuopTh9QAbN26kVKlStG7dOvloBuRMeH1G0gymkVLeT/g3UAhh\nDzQG4oGTUsoHudQ/g6tevTrHjx/nvffew8/Pj0/cP2H++vmM+WMMTyOfMqLhCDA1hU2boGtX6NcP\n1q3TFIgEQgiWOTnR7MwZFt+7x0iHLtwuMRt6H+XWDHscZzkacA8VpeDITEJZ+/bt2bt3L6A5MxIU\nFJQs+AagXLn/xglbvnw5Hh4eeq+b6Pz581SvXl37+vnz50yZMoXDhw+zYsUKnes4Ozvz77//ZrS7\n2SbDhDIhxBBgMnA4YdYSIcQ0KeXPOdqzPMTe3h5vb2969eqFl5cXH3b+kPmr5jP3+FxCIkOY2GIi\nwtxcE0HWqZNmfKLVq8HYWLuNIsbGbKtViyanT1PPyopajt9wuf0Q7ndviN0HdhStVdRwO6go+cT9\n+/dZtWoVLi4u/PXXX4wYMQIbGxvCw8Oxt7fPVEKZqakptROiaffs2UODBg1wcXHR2fbJkyc8fvwY\nc3PzTK8bGhpKsWLFtK8nTZrEkCFDKFeuHEIInesUK1aMoHQGu8xu+kRVjgNcpZQhAEKI0sBxoNAU\nAgArKyt27NjBxx9/zIoVK/j4/Y+ZOn8qm/03E/IihHnvzMPIwgK2b9fcVjp0KPz0k+aUUYIqlpas\nqVGDXhcvcrL+GxSxqorJvGMEDC+ByxEXNEmgipJ3eXtnz8+om5vMuFEKERERdOnSBS8vL0qXLo2t\nrS0eHh7069eP9957L8t9CQ0NZfXq1axbty7NNps2bcLZ2TlL65YqVYqwMM2DpGfPnuXPP//kzJkz\nAKlOCyV6/vw5pUqVysxuvBJ9CsFjIOmTU+EJ8wodExMTli1bRqVKlZg4cSKTPSbz0eiPOGl+koE7\nBvJTx58wLVIEdu2Cd9/VPCCwbJnm2YME7UqXZkS5cvTw92fXazO5HNEJs/iWPFj1ABzLpvPuimJ4\nWfkAzy6GCK9PdPjw4VRh8vqumzS83tvbm8DAQCpWrAhAeHg4cXFxXLp0iVOnTmnXye3wepFWRRJC\njEn48nWgLrA94bU7mruG/pfz3dNNCCHT6ndu+eWXXxg0aBCxsbF06dmFF++9wMzEjE3dN2Fpaglh\nYfDOO1C/PixalKwYxEtJ5wsXqGRhwfCXkzANq84j97a8XNaQaYvMOHw4nTdWlBwkhEjzr1RDW758\nOREREYwePRrQXOzt3r07ly5lfWTfRYsW0axZMxwcHLh9+zaRkZG8+eab2vD6xFM3rq6uzJ07l9at\nW2e4bkoLFizg8uXLLF++nMjISO3RgZSSuXPnEhgYyLJly7QFLioqCgcHB/z9/bG3t8/SfqX1fUyY\nn+qwLr27hooBVmiGmNgOyIRpB3AjS70rQPr168e+ffsoXrw4v2/+nbAfwzDDjHbr2/Es6hkUKwZ7\n98I//2hyDZJ8U4yEYG2NGux78oRzxT7lYexSbAaawbLr6byjohRuffr0ISQkBC8vL3bs2EFQUBAu\nLi6sXLmSF2k84Z8eHx8fPDw8aNiwIeXKlaNp06baC749evTg7Nmz2rbW1tY4ODjotW5K/fv3x8vL\ni6ioKCwtLbG1tcXW1hY7OzusrKywtLTUFgGAXbt20apVqywXgaxI84ggL8sLRwSJzp07R/v27bl3\n7x7Va1Sn4ZSG+D/zZ1/ffdgWtYWnT6F1a82opV9/nezI4EJ4OK3+/ZcdJVZSxrgkd926s7ZsDVad\nzb1zg4qSVF4+IsjPJk6ciK2tLaNGjcqwbZMmTVi5ciU1a9bM8vtl9ohAFYJscPfuXd59910uXLiA\nnb0dneZ3wvuxNwf6HaBSyUoQEgKtWkGXLjB1arJ1Nz18yJzr/7AwbgCxxw4T4hlJ16cNMTIvlJlB\nioGpQlAwZOepIUVP5cuXx8fHh7feeovgB8FsGLqBt6zeosWqFlx8dBFKl4aDB+G33zS5Bkn0srXl\nzTLO+Jh0RjZcxOMiRbg9+zZDh0L58pqpXTsD7ZiiKIWCKgTZpESJEuzdu5e+ffsSERHBT0N/oq1x\nW95a8xa+93zB1hb+/FPzsFnCE4mJZjs6csSsP9H8zok34rm76C5xt18weTL8+CM8fGignVIUpVDQ\nJ6HMVggxUQjxoxBiVcK0Mjc6l9+YmZmxdu1aJkyYQFxcHD+P+hm35268t+E9Dt44CPb2mmKwYgUs\nWKBdz9TIiDW1m7KNXrze/SsqTahEqzMBlLaW2NkZcIcURSkU9Dki2AEUBw4Ae5JMig5CCGbOnMny\n5csxMjJi0/RN1Ltej/e3vs/Wi1vBwQEOHdJEkS1dql3PzswMVyZRofwJnvW5icXLWMz+1u9Q4Pp1\nqFhRM6XxcKOiKEqa9HmgzFJK+UWO96SAGTZsGA4ODvTs2ZP9y/fTyL0Rn4hPeBL5hKH1h2qODNzc\nNOMUDRsGgDNlWXdhLCEmY/Fr8Ae9115EdrEGTNN9r9hYzdBGW7dqbk5SFEXJDH2OCHYLITrkeE8K\noA4dOnDkyBFsbW3x3eFLiW0lmO49ndk+s6FKFU0xmD5dM0Rpgps7PLAXD7k9YCvRjW2IWqzfIxum\nppoLy4qiKJmlTyH4DNglhIgSQoQlTLqT25VUGjRowPHjx3FycuLK8SvErojlR98fGXdgHPK11zR3\nE331leYiMiDjzahf7Ru6283j90GmxP4dQuXwZwbeC0VRCrIMC4GU0kpKaSSltJBSFkuYiudG5woK\nR0dHjh07RrNmzQgKCOLhnIfsPr+bobuGElvtNThwAD7/nDLemwFwsOuJ1TNz7tqs4/6osvS4E0B8\nTLyB90JRlIIqzUIghHBO+Leerin3ulgwlC5dmgMHDtCtWzfCgsO4Oukq/1z5h16/9SLKyRH276fa\n4k9p/uh3hDDihM8cPpFr8Gh2m6eWJtxdcNfQu6AoSgGV3hHB6IR/5wPzdExKJllaWrJp0yY+++wz\nYl/EcuHLCwRcCaDD+g6EVa/Cua+98AgYDrt2cedOG0zjK/G5yXHmjY/j9pzbBPwVyaVLcOkSPE5j\n/NfYWLRtLl2Ca9dydx8VRUnN0OH1GdKVX5nXJwpArur8+fOlEEIikM7jnGX95fXl7/sfyQ/r+UpZ\npoyc9sZeuX37P9Lb20GWWewrF3mckkuK/SsrV4qX1tZSTp/+37YuX5bSyUnKkBApa9T4b3J0lLJK\nFcPto5L/FITfrczQN4BeSin/+ecfOWvWLO1rR0dHaWZmJm1tbeWaNWvSXC9peH10dLQcNGiQrFSp\nkixWrJh0cXGRe/fuTdbeEOH16sliA/Hw8GDz5s2Ym5lzac4lnp1+xsgzzfi3jD1s386nfv2xPv0c\nY+PGfHDnT37oGkMJ43A2ffKYESN0b9PaOvnRwIEDubtPipKfXLt2jZCQEMaMGcPSpUsZMWIEN2/e\n1Nk2Pj6eyZMnE5Mkk3z8+PFcvXqVe/fupZsdsHr1ajp06IC5uTmxsbFUrFiRv/76i+fPnzNjxgx6\n9uzJrVu3tO07duzI4cOHCQ4Ozr6dzYAqBAbUvXt3Dh48iLW1NdeWXyP6+AvOuDTjilNp5jb+jUYL\n+1DqmjutW8xjc53yzB0fR9i3ARhHxxq664qS72UmgH7Lli20adMm2UBuZmZmVKxYEROT9B/HShpe\nX6RIEaZMmaINpunQoQNVqlTh9OnT2vZ5KrzeEIQQVYCJQAkpZQ9D9yc3NG/enGPHjtGuXTsCdwRi\ncr8MLUq1oHaFffzTahNNR/fmqEcTatb8ntCLn/F3gytU+vs699+rzt27cPgw5GK0qaLkadkdXg/w\n6NEjjI2NKVOmDBEREdr5J0+eJDo6mufPn+Pk5ESnTp10vmfK8PqkgoODCQgIoFatWsnm58Xw+ubA\nWSlluBCiH+AKfCelvJXBqpkmpbwJDBFCbMnubedl1atX5/jx47z55nsEnPQjEit83m7N4Sq/U3b6\nOjrN7MO1Fn9R8eZI7k+ww7VPEKGuZTl7tjjjx8Nbb0EaP4OKUqAYIrx+27ZtDBs2jLVr1yab37p1\na7p06QKAi4sLLVu2pGTJkqnWTxlenygmJoYPPviAAQMG4OTklGxZbofX63Nq6AcgQgjxOpo7ia4D\na9Nf5T9CiJVCiGAhxPkU89sJIS4LIa4KIQr9EBb29vYsWOCNtXV7wk+GE7sxgjk33dlVJZKZlr9g\nvyuGXm+NZJhdVTb3tcR+3wWIk7i6wi+/wLffGnoPlMJACJEtU1YkhtcPHz6cDh060L17dzw8PPjj\njz+wtrbO8j6lF0B/4sQJGjdurHN8/6TZx6VKlcLb21vn9pOG1yeKj4+nX79+WFhYsGTJklTr5MXw\n+lgppRRCdAaWSil/EkIMzsR7rAIWk6R4CCGMgSVAG+AecFIIsVNKmfXw0QLA0tKK2rV3UKPGx6xY\nsYIXP8YwWfandO1FzH9zOXXl/7C4vhvb0u8QUvwEYXsDAN2HnIqSE1J+GOYmQ4TXnzx5khcvXrB/\n/36OHj1KZGQkO3bsICwsjJ07d7J5s+Yh0IiIiDSvFSQNr098z8GDB/Po0SO8vLwwNjZOtU5uh9fr\nUwjChBATgL5Ai4QP8fRHQUtCSvm3EKJyitmNgGtSykAAIcRGwF0IEQzMAlyEEF9IKWfr+z4FhZGR\nCcuWLcPHpxIXL04k4vvnRPYfyeeRnlT17o2bQy+q3D/JrfbOtPz+IqKNHZD6cFRRCpqYmJhk5/Aj\nIiIwNjbWnp4BMnVqCGDx4sX06NGDqKgofH19iYyMpFKlStrw+pEjR2rbenp6IoTA3d0dHx8fhg8f\nDsCLFy949OgRb731ls73aN++PUeOHOH9998HYMSIEVy+fJmDBw9ibm6eqn1UVBSnT5/ml19+0Xs/\nXpU+haAX0AcYJKV8IISoCLzqiQgH4E6S13eBxlLKJ8BwfTbg6emp/drNzQ03N7dX7FLeIYSgRo0J\ntG5dgR9+GETsz2EsjvWkwoshvN7Ygp5bWvG729/c7WBNkxvnCY9tilUGdy4oSn7Xp08fZs+ejZeX\nFzExMRQtWlQbXt+7d2+KFCmSqe0lBtAnHuUIIbh9+zagCa//+eefcXV1BWDz5s3s3LkTIQS1atWi\ne/furF+/noULF3Lr1i02btyY5vv3798fFxcXoqKiCA4OZsWKFVhYWCQLp1+xYgV9+vQBsje83tvb\nO81TVknplVmc8Bd9VSnlQSFEEcBESqn3wHMJ6++SUtZJeN0NaCelHJrwui+aQjAyzY0k35405CFq\nTjl8GKZN0/zbrRu8/z6UKHGQrl27EhYXhtVwKya2qkW9xxE0GfWUU7P+5OHUh/h5FmXOCNdU515v\n3IA2bTT/Koo+VGZxzsjr4fX6JJQNA7YAyxNmlQd+z3IPNe4BFZK8roDmqEBJoU2bNvj4+OBQyoHw\npeF8ve8cobY3+KtfT5p6vs2/5UvRdGYYCwNuG7qriqKkYebMmXoVAdBcoH6VIpAV+tw19DHQHHgO\nIKUMAGxf8X1PAdWEEJWFEGZoTj/tfMVtFlh169blxIkT1K5Wm+fLI/npn5c8aPMjF7uP5rPL71Gh\nvgV3pt3C++lTQ3dVUZR8SJ9CEC2ljE58IYQwAfQ+dhRC/AocA5yEEHeEEAOllLHAJ8B+4CKwqbDf\nMZTo8mUYPBhSPuBYvnx5fHx8eKvlWxycGYuICWdcpblsfe1DXv93OB33wbjd/tyNijJMxxVFybcy\nvEYghPgWCAX6o/nw/gi4KKWcmPPdS7NPcsqUKQXuIvH9+5DwoCOgeVCsSpXkbV6+fEnnzoN58nQd\nw8fAuPPluGo+mIjv7uNXdgAzf7TgSP16mBsZqWsESqapawQFQ8rvY+JF46lTp+q8RqBPITAGBgOJ\nabj7gZ8MebW2oF4s1peUkq++moit7dfsuw4nrEsSEDuIO/Nrst2jJg/7W7OsenVVCJRMU4WgYMjs\nxWJ9CkFRIEpKGZfw2hgwl1K+yJ4uZ15hLwSJ1q79AnPzOfRbACZvW3IhcCT3N7Zg9O+lGF7fEbfI\nsqoQKJmiCkHBkO13DQGHAMskr4sAB7PcQyXb9O8/m9KlG9PjNVMit0XiXGEBljWvsWjYGb64do1z\n0SpaWlGUjOlzRHBWSumS0bzcpI4I/hMefoFTp95kwAAjbpk8pqh7EbyWbyC02T98NKEDYf0a8sVw\nMypWhL59Dd1bJa/L6jhASt6T3UcEEUKI+kk21ACIfKUeZgNPT0+9npgr6KysamNn9x7btvXGydiJ\niM0v+PKd2ZQ4Wo8P/vgTm3nnOX9Rks7wK4qipSu9Sk35c0rK29s72WgMKelzRNAQ2AgkjolaFugl\npTz1Kj9wr0IdESQXGRmIn199qlb1oVu3IRy7doypxafSIu4FMyZWo2zFZtyZUYMjRwzdU0VRDCnL\nF4sTVjZDM8ylBK5IKWMyWCVHqUKQ2tWrnwHxlC8/m379+nHkwBF+fPkj0mUJn37+CaX3NuDsj+UN\n3U1FUQzoVQvBG0AVNIPUJabH651JkN1UIUjt5cuH+Po6U7/+KczMKjJ27FjufH+HZjbNiKuzmVkj\nvuJYi6bUeIVx2xVFyd9e5fbRdYAjcBaIS5wv9RwgLieoQqDbzZueREVdx9lZM3ztgnkLKPp5UbY1\n20qlKlb81akPvu3bUyyTozQqilIwvEohuATUzEufvKoQ6BYbG8Y//1Tj9df/wMqqLgDbZ28nfnw8\nQ7oMpKHLhxS1c2TrgIEIHeOgJ/rjD1i1SvN1ly7we8IQg6tWgYVFTu+Foig55VXuGrqA5gJxnqLu\nGmq9k90AACAASURBVErNxKQYlSp9yc2b/43+0fmLzsS3LM2QHSP589xCLpjCzFkzICbtyzzXrsHj\nx7BvH5w/Dw8fwrZtEBubG3uhKEp2y467hrwBF8AXSBx8TkopDRaXro4I0hYfH80//1TH2XkdJUs2\nB8DbK5bQzsf4vvQMvJudxarvD/y0dy9dly4DHYE2338PFy7Ajh0waBA8egTr1sGDB2Blldt7pChK\ndnmVIwJPoDMwE5iXZFLyICMjc6pUmcqNG+O19xKLoibsr1oDz+KeuNyoTPjvMxn67rsc79UV4uIy\n2KKiKAVdhoVASukNBAKmCV/7AmdytFfKK7Gz60tsbCghIXu08y7alMGqqhXrO63n7cdlee6zhYHt\n3+X0Oy0hPt6AvVUUxdAMlVCm5CAhjHF0nMXNm1+SMFYgCEG1pdUIXhLMpgWbGPC8ONeeBTO5WVOO\nN3NRxUBRCjFDJZQpOax06Y4YGxcjOPhX7TzLypZU/LwiN0bdYPmyZUx4aYFX9drscKzA3gbViFen\niRSlUMrxhDLFMIQQODp+Q2DgJJJ8+yg/ujw3fKPpVPwxS+aMp/8lc+b0GsKDImasrlaBqEj9hpF6\n+22wtoZJk5LPX7NGM//NN7NzbxRFyUn6FIIjQoiJQBEhxNtoThPtytluZUzdPpqxkiVbUqRITWCF\ndp6RqRE+DZ0YZXKN4saxuNfrg/264gwZNZVKRi9YW68cEeEhGW77+XNo0QJS1o3oaKhUSbNcUZS8\nIaPbR/UpBF8Aj4DzwIeAF/BVdnTuVXh6ehaomMqc4ug4C5iJmVmYdt5D2xJEu5SmzwtNYk3F262p\nfrciHTym0yDkOXErKvLs+c0Mt21mpnu+jjtSFUUxIDc3t6wXgoTTQBellCuklN0Tph/VTfz5h5XV\n60BrmjZdmGx+aE9HGkU9xviq5k/3VjcaU75cIzp+OJaOT15Q1bc60S/9DNBjRVFyW7qFQEoZC1wR\nQlTKpf4oOWIajRt/x8uXj7Rz4q1M2VDMkSLLAjCS8RghGBzkivG77nRuN5g+D2MYbNuIgGu7Ddhv\nRVFygz6nhqwBfyHEIfH/9u48Lqqqf+D458wMu4i4objhAoqoWT0ZuPCI5q6ZmBnYz720Mi0r01bU\n6tFKfZ4ytcXU0tzK3TQ1wyU1W2xhUTZFVCAVcQFBGO7vj0EF2WbYke/79ZpXzL3n3jnHLvPl3nPO\n9yi1Nfu1pawrJkqPUi0JCXmc06f/k2v7QVsXNEcDnRPOAlDDaM3G9u2JHzuagPajePJ0Fq4nBhEW\nsaQiqi2EKCfmBILXgYHALGRmcZW1f//rJCSsIC3t9O2NSpH6lAd+52Kxu5oGwP2OjgxJbMmpl5+m\nn90knouA9qlPk870PKseCSHuDub0EXyqaVrwHS9Z66qKuXatAa6uT3PqVFCu7Vmu9hxyacw9B6Ju\nbfO53BDX87U48/IEZv5rFtP+gjH3zWXMpEAyCklWJ4SomszpIzgufQRV3/79cN99L3Pq1DZq1QrL\ntS/YtSk1k1Jw/PvCrW1dj7ljdE7ny0YjCfrXEt44Co5n1tBjSA+u5DM21N3dNIroywpbrkgIUVxV\nto9A5hGYr1s30/j+WbOcOH78Fe6//7Vc+406Hce6e+D6TSSZ10y5pvVZOo4P96LhM2cZ/u0wBmjr\neecnRbPMg3Tq34kbN87mOkdGBnh7S6YKISqj0phH8AZ5+wjml0blSkLmEZhPpzP9tW4wQHT0M9Sp\n8yvW1kdylTnf2JmUVrU4FXTq1rYWjrZ83daT0RHhRLkO5NqaHcz7UUcH5xOEJN3P5cuheT5HCFH5\nlGgeAZiyj+b3KsU6inJkNNpx7FgQtWpN585MIfGPtCTxy0Rs4m5PPuvu7My0pk1JfzWE2j0ewmb3\nj3y8W8/AexPZe/JBTp/+sZxbIIQobeZkH72mlLqa/UpXSmUppSSBQBUWFTUKvT6RDh125dpudLSm\n+bvNcV0dgcq6HSReaNwYFW/H1NOROHXrhuPe/XyyzYp+7VP4JqQXX3/9dXk3QQhRisy5I6ihaZqj\npmmOgB3gDywq85qJMqNpBpKT3yEgYDqQ+6F+w7EN0fSKFmHnbm1TSmHzYWt+vXaFz+Ljse3cmZc9\n9vPFVht6exoZ8cEILl+eI8NLhaiiLHqqq2lalqZpm4C+ZVQfUU6uXx9CZqYVVlbrcm1XOkV8QGu8\nfjlFevztrKUqzcCqVu14/eRJfr5yheOO3vy32498vc2ens0h2XMGUVHP3l7/QAhRZZjzaGhojtcw\npdQcwLxcxaISU6xZMwc7u9fR6XLPDUh3deCkZ0OiXojKtd3d1p7PWrdmWGgoGQ43iKjjw+ohO9m6\npwa+LnDOczHHTwzBaEwtz4YIIUrInDuCQZhGDQ0EegNXgcFlWSlRPkJDe5CV1YIHHliaZ1/Y/c24\nevQqSd8n5do+uG5dRrq4EDUijCyVxcnG3bDbtIVvtjjRxV5HctetnIjsTkbGP+XVDCFECZnTRzBa\n07Qx2a8nNU17R9M0+S2/S1y//h969pyFXp+Sa7vRSo/7IncinonAeN30uOfAAdi3Dwacb44uU/G3\nTwypqbBP58ekmuvYssMJb6M1qf1+ISLGm8jIyIpokhDCQuY8GlqhlKqV472zUuqLsq1W0WRCWekw\nGu/n5MlutGr1YZ59dfrWwfFfjsS+HUvXrjBnDrz0Eox8QtFiTVvOtbjAz3b/0LcvJHTozfz2K9m+\ny5FOl2uQ8ehJvH29OXz4cAW0SgiRU2lMKLtH07Tkm280TbsE3FfyqpWMTCgrPbt3z8bdfT56fVKe\nfa3+24r4T+P5dn4K+/bBzZGiVqlWeO/0IqRnJC4Pmvad6dCfuS2X8t1eO/51rj5JA5Pw6+/Hxo0b\ny7lFQoicSjyhDFBKqdo53tQG9CWvmqgsLlzw4OxZf+rVm5tnn01DG9yC3IiYGIGWlXt4qPMFRzz3\nteTc0yFczjSlpjhQezBzmy3mu306/DM8SA9Mx3+0Px999FG5tEUIYTlzAsE84LBSarZS6m3gMPB+\n2VZLlLfw8Ldwdv6ctLQzefa5TnQlKy2LhOUJefY1CWuAfbgzo8LD0bJnKv/oPJT5jebzza6rjNd3\ngTEwOWgyL730ElmSjEiISseczuIvMU0i+wdIAIZkbxN3kbQ0V5KSniQ2dlaefUqv8PjEg5gZMRiT\nbuTZX39dK/7JyODv9rfXOtjpHICaM4fPNpxiRs2hMAbmrZpHQEAAaWlpZdoWIYRlzOks9gbiNE37\nSNO0hcAZpdSDZV81Ud7On3+F5OR96HRX8+xzvNcRlxEuXJkbnWefMupY7+XFcc+zJLvn6GcYORJm\nzuTd5UdZ0GQSaqRi3c/r6N27N0lJefsjhBAVw5xHQ0swzR24KSV7m7jLZGU506lTGFlZjvnud5vl\nRvqRZDxSL+XZ18jGBt99bYl6NJx05xzzDceNg1df5fkFW1nZ4R10j+s4kHiALl26cOrUqTJqiRDC\nEmalmNByJJHRTDkEpLP4LqVUwf9rDTUMOL3eiifOR6DP51l/g8RaNNrXlOiRoWRZ5Ug1MXEiTJ1K\n4Jufs/Pfn2PwN3Dc+jg+Pj78/vvvZdEMIYQFzAkEJ5VSk5VSVkopa6XUFCCmrCsmKie7XvVItLan\ne/zpfPc3ONQY2/P2xD0akTsJ3eTJ8Mwz9JryH44M2oBNfxsSmiXg6+vLzp07y6n2Qoj8mBMIJgJd\ngLPAGcAbeKosKyUqt9V13emSeAanq3lzCikUzda35nrjayw5dy73zhdfhDFjuH/cNMKe+BHHno6k\nPJDCgIEDWLo0b5oLIUT5MBRVQNO0RGB4OdRFVBFJVrb82LAZ/z4WyboHOgAq1359hp7my7x4q9kx\nOtaogY+T0+2dM2bAjRu0GP4kUVsP08HQk0T7RMY/OZ7Tp08TFBSEUrnPJ4QoW+aMGmqilNqolDqf\n/fpWKdW4PConKq+fGjTCLv0Gbc7ln3bK5oI9X7RuzWNhYSSkp+fe+eabMHgw9YeMIGLCEdy7uIM/\nzHpnFmPHjiUjIyPfcwohyoY5j4aWAVsA1+zX1uxt4i7w7ru53/v7w+LFRR+XpXQE39canz+jsTfm\n/8U9sG5dxjZowPCwMDJydC4H71M8/NfbbErpRWLbRwm8tI9OXTuhC9SxfNVyBgwYwJUrV/jwQ3j4\nYfjttxI0MB9XrsBDD4GlT6PWrjUdJ3FK3G3MCQT1NE1bpmlaRvZrOVC/jOtVJEk6V3LLl8Mbb8D7\nOeaJb9sGQUGmgT5F8RlXkxp96jLPs+CxA2+5ueGg1/NKzO0yZ8/Crt2K4bHv8d3lrgz/dAgHR31H\nH98+WI21YveB3fj6+nLo0Fm2boXz54vfxvxkZMAPP0BIiGXHRUaajpOF2ERVUxpJ5y4qpf5PKaVX\nShmUUk8AF0qrgsUlSedKrk8f01/c3brl3j5oEHToUPTxrVqB39fNsf3tIpcPX863jE4pVnp6svnC\nBVYnJt7abjAASjFVLSCy5n1YDRrMtuFreMLvCawnWPNn9J9s2+YDhBa/gUWQrghRXZRG0rmxwGOY\n0kvEA8OAMaVROVH1WdWyouX8lkRMiEAVkEeotpUVG9q1Y3JUFH9fu5Z7p1IsbL0QPD3RDXqYpQM/\n4sW+L2L7jC0p1nFAF/76K7jM2yFEdWZOrqFTmqYN0jStXvZrsKZp+Q8iF9VS/eH1sW5oTduwvAnr\nbrqnRg0WtGyJf2goKSr3Q3ZN6eCTT8DNDfXII7z77zeYM3gO+vG2UP8yr73Wh9WrV5d1M4Sotixa\nvF6I/Cil8FjkQfuQ09RKLzih3BMNGtC3dm2WOB9HU3c8aNfpTL239euDvz9T7ptAl9SlMNKOzIY3\nCAwMZO7cubknqQkhSoUEAlEq7FraEda2CYNORhbamzqvZUtSdBlkDI/Nu1OvhxUroEYNGDaM1tce\nhU0bsB3tAK1g+vTpPPvssxiNxrzHCiGKTQKBKDUhXk2ok3addskFjyWw1ul4LsmLzH7nyPrXxbwF\nDAbTMmh6PROCH8cQ1ZN32u/CaZQTho4GFi9ejL+/P6mpeWc1CyGKx5wJZa/n+Nm2bKsjqrIsvY5N\nLTwYfDoKG2NmgeWcs2ywmduWzJeOozW8nreAlRWsXYsh6wYreQIv+07sf3I/zsOdsfe1Z8uWLfj5\n+fHPP/lPZhNCWKbAQKCUmq6U6oxplNBNh8q+SqIqO1WzFhE1nekXf7LQcvqwWuhXN0MLCsFoyOdR\nj40Ni/y+oRbJtJ83mg51vTgy4Qh1BtWh1sO1OHr0KD4+PkRGRpZRS4SoPgq7IziOKQg0V0odVEp9\nBtRVSrUpn6qJqmp74xbcm/wPV3/Pu8BNTrrNjeCUA+H9IvLtBM402PIIm7BJiofx42nh5MbPE36m\n4UMNqT+yPjEnY/Dx8eHw4cNl1RQhqoXCAkEyMAOIBroDHwIa8IpSSn7zRIFSrazZ1rAFERMi0IwF\ndxwrFGp+a665XOPjs2fzLZOGHb8HbYHoaJg4kYYOLvw0/idadG1Bo2cbcfHSRXr06MHGjRvLqjlC\n3PUKCwR9gO1AS0wL2HcCUjVNG6Npmk95VE5UXb/UboDOQcfZRfl/wd+k0vV0+KYds2Jj+ely/rOT\njbYOsH27KSfE5Mk429bih9E/0P7B9jR9uSlpmWkMHTqUjz76qCyaIsRdr8BAoGnaDE3TegInga8w\npayuq5T6SSm1tbwqKKoopfBY4kHsrFjSz6YXWtQ+2Y7lbdowPDSU+Dszld7k6Ag7dsDRozB1KvYG\nOzYHbKbLA11we90NzVpj8uTJvPzyy2QVMMNZCJE/c4aPfq9p2q+apn0CnNE0rQumtBNCFMqhjQOu\nE12Jej6qyLL969ThSVdXHrsjU2kuTk7w/fewbx9Mn461zoqV/ivp/6/+NHurGfqaej744AMCAwNJ\nSyt4YpsQIjdzUkxMy/F2dPa2Us4HKe5WTV9tytVjV7n4XT5zBu7wRrNmOOn1vBQdXXAhZ2fYvdt0\nd/DWW+iUjoX9FzLKexQNXm2Ag6sDa9eupU+fPiQlJZViS4S4e1k0oUzTtD/LqiLi7qS30+Ox2IPI\nZyMxphY+I1inFF95erL94kW+zpGpNI86dWDPHvj2W5g9G6UUM/1mMq37NBynOFKvbT32799P165d\niY3NZwazECIXmVlczWRlmV5lTdPgZiaI2r1qU9OnJqdmnsJoLPzznbMzlU6JiuKC07U8+7OyTOfN\nrF0f464f0FatgrlzAZj84GTe7/s+jIIW3VoQHh6Ot7c3v//+e6F1NRpNr4IyY9z8zOIoybFClJcq\nGwhkYRrLKQWffQarV5d9Lv6//gJr69vvW85vSfTCBFoZrvHcc4Uf26FGDTr/1oot94dAjdyZShcv\nNmWhsLICQ+MGpG//wdSoBQsAeKLDEywdvJTLAy9zz5B7SEhIwNfXl507dxb4eVZWpnMWFC+mTMle\nP6EYBg40pU4SoiKVxsI0lZIsTGO5p5++/dfv//1f2X3OvffmXc7RpoENp/zceMUQgTIjg6hXggud\nVR3qLggni0LKN2oEe/fCRx/Bxx8DMKj1IDYM38A5n3N0m9CNlJQUBg4cyNJC1qbs2LHgjyjJX/Rn\nzoD0W4uKVhoL0whRKuLvd0UBD2XEm1W+b2xLjDZGVuuLeM7ftKkpGLz3Hnz6KQC+zXz5/onviXKP\not9r/TAajYwfP5633npLUlkLcQcJBKL86BSfOngwIu0kTtqNIovrNR33bGrLDsM5tl8sYtSRm5sp\nGMyeDcuWAXBvw3vZN3of4fXCGfLBEHQ6HbNmzWLy5LGArEAvxE0SCES5ijXU4EerBoxLL2SIaA42\nKTa8esOLscePE309n0ylObVsaRpN9PrrsHIlAO513Dk45iAn7E7wyMePYGdvx+rVy4EBpKdfKVlj\nhLhLSCAQ5W6NrRteWcl0zLpkVvm2mhNvurkxJCSEdFXEA/vWrU3zDKZNg7VrAWhUsxH7R+/nrOEs\nvT7sRZ16dYHdrF/vy7lz50rYGiGqPgkEotylKz2fWLvzXGYEVpp5PbHPuLrSsUYNVruegMI6jwHa\ntjXNQJ4yBTZsAKCOfR32jNxDqnUqHWbfD4aWnD//J97e3mhaaAlbJETVJoFAVIijhrqcUg4EaKfN\nKq+UYomHBwk2qeBfeCI7ANq3N80+fvpp2GpKjVXDugbbArbh6OAIga64NO1EXFwc0IWrV4OL3xgh\nqjgJBKLCLDa04hHOUjvFvGUn7fV6xsd5wYhYaJ9c9AH33gvbtsG4caagANgYbPiszxq41AZtpJEB\nwwYAl4mJ6cPOnatL0Bohqi4JBKLCXFC2rFTN6Hs8/4Vp8lM3ww7mtIE3wqBO4VlNAXjgAdi8GUaN\nMnUkA3qdHrZ+QrPMXsT8OwZqjkXTbvD664HMnTtXhpeKakcCgahQm2iEXWYmiSsLyS10p1/qwBZX\nCArlhjn5Mnx8THmJAgLg1mx0Rdf0/zD23rEw9gfqe04HYPr06Tz77LMYJS+EqEYkEIgKlaV0fNfG\ng5hpMWQkWTC2f1UzuGzFK7HmDUOlWzdYtw6GDcNw5OCtzS91fgn2vUXS4OU89+4cbGxsWLx4Mf7+\n/qSmmvfISoiqTgKBqHDxTjWp92g9Yl6JMf8gTcF/2rD7chJfJSSYd4yfH6xaheNofx7kyK3N6o8x\nNP5zESuZx/vr38fZ2ZktW7bg5+fHP//8Y2FrhKh6JBCISqH52825uOMiyQfN6AS+KcWKtR5eTI2O\n5o+rV807pndvrn20nC08TOOEX29trhU/hHfvW8Ps47N5e/3buLm5cfToUTp37kxycqSFrRGiapFA\nICoFg5OBVgtaETExgqwb5ufJ9rKvwUetWuEfGkrSnZnuCpDRqz9P8hlPbh4Af/xxa/sDdXuwPXA7\ns47N4oUVL3DfffcRHR3Nxo2dIccdhBB3GwkEotKo92g9bJvaEjc/zqLjHndx4ZG6dRkRHo7RzBE/\nWxjMtz0WQb9+eGkht7Y/0OgBgkcHM+/3eQybP4x+/fqRlnYB8CMsbJNF9RKiqpBAICoNpRTuH7sT\n90EcdslF5BW6w9wWLUg1Gpl56pTZx/zlPhTmz+d7etM8LfzW9jZ123BgzAGW/72cjlM70sZzHJDG\nmjX+wEKL6iVEVSCBQFQqds3taPJSE9p8H1nwkmH5sNLpWOflxbKEBLZeuGD+BwYEMIM5LI7uhU1s\nxK3NTZ2acmDMAXbF7EI3yABqZvb8gueYPv1lsspjmTchyokEAlHpNHmxCXZX0qgTct6i41ysrVnf\nti3jTpwg0oKhn1+pkSxpMBOPZx6C6NvDUes51GPvqL0kG07A0DAeGboUMLBgwQcEBgaSJivOiLuE\nBAJR6eisdIT1a43bd1HYZGRadKy3kxMz3dzwDw0lxYJJYZvqjCN+zKvQsyfkWPC+pk1NBlzaAYY0\nQtqvA6tvcXR0ZO3atfTp04ekpCSL6idEZSSBQFRKyY2dSPaog2/USYuPnejqyr8cHRl/4oRF6SIu\nPDoRpk6FHj1Ma0xmM2AL677BUTWEkXPYuHMbrq6u7N+/n65duxIbW8QKakJUchIIRKUV26cFrRPO\no05YtoCMUopF7u5EpKby3xxf6GaZPBmeecYUDHKuVZBl4BG1FOJ8mPLHs2zau4l27doRHh6Ot7c3\nvxe08r0QVYAEAlFpZdpb8WPrFlh9FEFWpmWds3Z6Pd96eTHn9Gn2JVswSQ3gxRdhzBjTY6LE2zmQ\ndEoHuz7g8baBDN85nBVbV+Dn50dCQgK+vr7s3LnTss8RopKQQCAqtdCGLmg1DJxdaMYaBHdws7Pj\nK09PAsLCOJtuRqbSnGbMMCWp69mTGtdzdlorpnWewfSu0xm0YRDvrXiPESNGkJKSwsCBA1m6dKnF\n9RSiokkgEJWbUmRO8iD27VjSzlg+Sqd37dpMatSIR0NDSbd0yOcbb8Ajj/Dc1l44k7tT+Kn7n2JB\nnwX0X9Ofp95+ihkzZmA0Ghk/fjxBQUGSylpUKZUqECilHJRSK5RSnyqlAiu6PqJy0Brb02hSI6Im\nRxXr+OlNm+JiZcULURYerxTMns3xJr3ZTS9srud+xPSY12N8NeQrhq4fStcxXVm8eDE6nY6ZM2cy\nbtw4MsxMeSFERatUgQDwB9ZpmvYU8HBFV0ZUHk2nNyUlJIUaf1swWSybTilWeHryw6VLLI+Pt+xg\npdjoPZcDdCPwyz44krvjuk+rPmwN2MqYzWOo2aUmmzZtwt7enmXLljFw4ECMRss6uoWoCGUeCJRS\nXyilEpVSf9+xva9S6rhSKlIp9Ur25kbAzUQzsjKIuEVvq8djiQcN10diW4xLw8lgYGO7drwcE8Of\n183MVHqTUrzAAuJd72cH/eDatVy7vRt788PIH5i2exqnXU4THBxMvXr12LVrFzExvsC5/M8rRCVR\nHncEy4C+OTcopfSYkrb0BdoCAUopT+AM0KQc6yaqEOcezqS2qsUoThXr+LYODixyd2dMbCjUtPSx\njWLngIWE0RbDIwMgJSXX3nb123FgzAH++/N/2ZG6g0OHDuHu7k5a2p+AN6GhocWqsxDlocy/bDVN\nOwBcumNzJyBK07RTmqZlAGuAwcAGYKhSahGwpazrJqqehCEt6U0CLblWdOF8DKtfn0FO9eD1MLKw\nsENXp2MCn6C5NYfBg+F67sR4zZ2bc3DMQTaEb+DDyA85+NNB7O19gDi6dOlC8K1lMoWoXCrqr+6c\nj4DAdCfQSNO0VE3Txmqa9oymaasrqG7V1rx5YE7/ZlwcmPsH7rBh4O1ten32We59N26Y/nvkiGn/\niy/mf47nn799jjlLrFlKc14gAoWGry+8+y489hhYW5vKbN5sOu6ZZ6B9e9PPrVub+n69vcG4pDno\nNQ63KXzW8qRJ0L//7fdLloCGDuMnS8HFhcyHh+DZPI1vvrldxqWGC8Gjg1m3/xhd503FrcVOwJ/L\nly/Tp08fVq9ezZAh8M47EBUFLVvC8eOmY2fOhCeeMP0cGWnaF3dHRu5hw6BWLXj66bz1feMN6Nz5\n9vsVK0xtnjWr0Gbe8uKLpn9LS/3xB9x77+33jz1mWiLaXLt3w6BBln9uUZ5/HpYvL/3z3o0MFfS5\nJR5bFxQUdOvn7t27071795KeslqbOhVGjDD9bCjkqmjc+Pb6723bFn3eP/+Et9+Gpk1N7zdtMn2h\ndusGb70FR4/ClCnw889Qt27e4xcsgMuXc28LeqshWbsSWPHoOcJaNSImBvbuNQWxn3+Gm6tLHj0K\nIdlLDURE3G5bTKQOVrUlfNVvbL5QE8jng4HvvoOTJ29/6fbvnx3M9HpYsQLt0QDeOzWMM2e/Baxv\nHVfLthaJ874ncdhjOHYJhONrmPJsE/73v/8RGBgInObChWk8/LAiJgZu5sf79lv4+29YudK0LSYG\n7sxrdzPoLFkCixfn3vfLL3D48O33B7OXZv74Y3jzzXybmMv8+dCgAbz6atFlczp/Ptf6PqxfDy4u\nMHSoecf/9hts22bZZ5rjf/8zXWejR5f+uauK4OBgs+5EKyoQnOV2XwDZP1uUCyBnIBAl17y56VUU\nOzvTX9WW6NgRPDxMP//2m+m/tWubzlPUCpOenvnUtYViPq1ZtucPMr3rcv68Tb53MjY2ebfVrg3p\n6UCyNQN+9eJJp7/RGtsD9nnKWlnlft+kSY5tBgPXl36NcfMwei97HJ5Zm/uADHusN25EHzAG/q8/\nQXO20KxZM6ZOnQpMJzr6NEbjh4D+1iFKFf5vURRdAff3JT2vqLru/CN55syZ+ZarqEdDvwLuSik3\npZQ1MBzpExAWOIUD+kENqbsuuujCBWiYXJO3mzeHWSEYrS3LcgqAlRXDWYs+84bpmU5m7nOoLCsa\n//IlJHSk+/LuBD4ZyNq1awFr4uMXMXXqUMD8dNlClJXyGD66GjgEeCil4pRSYzRNywQmAd8DmKzI\nOgAACxhJREFUYcBaTdPCCzuPEHfSj2lGqldtixawudOTDRtCWE3iAizLVHrTDWzYOf4bSE42PYO4\nI/W1Qgc7/oe/pz9dl3WlU+9OwB4MBmd+/HEz0INLlyxbd0GI0lYeo4YCNE1z1TTNRtO0JpqmLcve\nvkPTtNaaprXSNO0/lp43KChIRmFUc8pWz9XODUr07EMpBf9zJ712GqsyLMxUms1oZWvq/IiPh/Hj\nIU8qC8Wb/36TKQ9OoduyblDfmY4df8LVtRnwM6NH+xAZGVnsNghRlODg4EIfp1fZsfpBQUHSQSxK\nhcrQ03y5F1/dOM2Pl+4c6WwmOzvYssW0wtnEiSjy5jWa1GkScx+aCyN7ktngMitXHgHu48yZaDp3\n7kxKypGSNUSIAnTv3v3uDARClCbrS7bMsvUkMDycuOIuQengANu3Q0gIHzI530dWge0DYfMXhHYY\nRKTxb2AfXbr048KFC8TE+AGbStQOIYpDAoEQ2R401GZKcTOV3uToCDt20ImjvGecmn//ReQA2v69\nkRm/PgFtdzB//hbGjx+PpqUB/ixcuLBE7RDCUhIIhMjhlaZNaWRjw5SSPLN3cqIP39Mtax/PJ0wn\nv2kzTpe78kmXXdD3ebac/YJPP/0UF5fZgMZzzz3H/PkvQz6Pl4QoCxU1j6DEbvYRSD+BKE1KKZa3\nacODv//O5c7xENGwWOdJxpkBVrsJvubHTGyAvNN7W9e6B5bvY3mD3jgcvEh9l9dITGyKwTCOL7/8\nAIgjPX05YFuSJglR5MSyKntHIJ3FoqzUNBjY4OXFBf8Y8Ch+GukkVYenmu9hKN/C7NkFFGrF0i4H\n+Trka+LbvQz8H9999x0ODo7AWsaO7UNSUlL+xwphJuksFqIYPB0cqL/KA2aGct36RrHPk2SoT09+\ngFWrYO7cfMvUs3Vl3+h9pNb5CQaPxa+nH198cQBw5Zdf9tO1a1diY2OLXQchiiKBQIgCOB6rBz/W\nZ493OOiKP2ktkQamZEiff25K6JOP2na1aXFwDzieY9j6Ybi1ag0cwd3di/DwcLy9vTl27Fix6yBE\nYSQQCFGYz5ujKQ3GFp6ptEiurqZgsHAhz5L/qCCd0QFWb8Vab82kw/3BxonVqw/i5+dHQkICvr6+\nmCbjC1G6JBAIUZgsHb2OtIWeiYTVKWEqiCZNYO9eXuZ9Bp37JP8yRmu+9v+aZjU8YFQPMq0z2LFj\nByNGjODatWvAAOCLktVDiDtU2UAgKSZEebFLt4YgL7a4R6A1SSn6gMK4udGDvYw8/TbOG/P/Qtfr\n9Lx2z2KI6kvAnm4kpiXy5ZdfMmPGDEwruI4DgoqVG0lUT5JiQojScKImvU42J/OtUK4ai5GpNIcY\nWjK1wx5cPn6DEazMt4xSCva+zeOtJtBtWTcikiJ49913gcWYfm1nMm7cODLMWUlIVHsyakiIUnJ/\noisqxIknI4+X+K/xOPvWnPxkN+8xDeddawssN6bNC8zqPovuy7vz67lfgYmY0lDYs2zZMgYOHMiV\nK8Uf4ioESCAQwiL6j92Ju5HO+3euIVkM6S3b0ofvafzBFNiwocByozqO4pOBn9BvVT9w+xEYBART\nr149du3aha+vL+fOnStxfUT1JYFACAuoDB1rWnux4MwZgq8UM1NpDiG0J+rDHfD00/z7ytYCyw1u\nM5j1w9bDsOHQZhPwAIcPH8bd3Z0///wTb29vrl41cyFpIe4ggUAICzWxsWWVpydPnQqH+sXMVJrD\n9Tb3wrZtBMWNoy87CizX3a07rNwBA56Gjsto2bIlhw4dwsfHh7i4OI4c6QIEl7g+ovqRQCBEMfRw\nduZZl8YwM5QMZSz6gKI88ABT3DazglGwZ0/B5eLvh+XB0H0m8w7No27duvzwww8MGTKEzMzLQB/W\nrFlT8vqIaqXKBgIZPioq2uT6TSDBlvUNo0rlfH85+JjyEgUE4PBLcMEFL7aGLw7w+bHPmbFnBra2\ntqxfv55mzSYDNwgICOC9996T4aXiFhk+KkQZUUrBe62JsbvMZ6XUWXuQbrBuHW7ThtGFgwUXvNKE\nA2MOsOfkHiZsmwAKPD3/C8wD4JVXXuHQoUmY5h2I6k6GjwpRlq4beDKuHa+ePMnR0hrG6edH7Dur\n2IA/tn8UvHxlXfu67B25l+hL0QR8G0CW7gYwlbVr12JtbU14+CJgKJqWWjr1EnctCQRClJDLDXs+\n9fBgWGgo528UP1NpTtc692Y0y2k48WH49dcCyznaOLI9cDuZWZn81noQWF/jscceY8+ePdjYOAOb\nuXSpB+fPlzA9hrirSSAQ1VRwqZ5tSL16jHBx4fGwMDKLu8zlHXbQn3/e/gwGDIA//iiwnK3BlnXD\n1mGX3hRGPsTF1It069aNgQN/ApqRkfEzPj4+RJZk1TWRr7uln1ICgaimgkv9jLObN0evFK+eLGGm\n0hxSHhoMixZBv34QElJgOYPOQLuTn0GsL77LfTl75SzOzp7AEQyG+4iOjqZz584cOVLwoyZhOQkE\nwmzlcbGU5mcU91wJCeYfl5hoXtmi6nLiRFHnMe9zSoNeKb729GTdP/+A7z+ld+KhQwkePx5696YN\n4QUWUyjY/R6j7hlF12VduWyIBBrg7LyPfv36ceHCBfz8/Ni0aVOB50hPDy6VKp85U7zzWHrtmVO+\npGXuli/7wkggKAcSCPK6GwMBQF1ra75t1w5eiCSraQkzleYQrNfD3LnsphfuRBRadlqXabzW7TW2\n1/s3NDiGTleDLVu2MH78eNLS0vD392fhwvzXRLhxI7hU6nv2bPHOI4GgYqiqONZYKVX1Ki2EEJWA\npmnqzm1VMhAIIYQoPfJoSAghqjkJBEIIUc1JIBBCiGpOAoEQQlRzd0UgUEo5KKVWKKU+VUoFVnR9\nRNWjlGqulPpcKbW+ousiqi6l1ODs76E1SqleFV0fc90Vo4aUUv8HJGmatl0ptUbTtMcruk6ialJK\nrdc0bVhF10NUbUqpWsAHmqaNr+i6mKPS3hEopb5QSiUqpf6+Y3tfpdRxpVSkUuqV7M2NgJuLyEre\nXQFYfA0Jka9iXkevA/nP2quEKm0gAJYBfXNuUErpMf3j9gXaAgFKKU/gDNAku1hlbpMoX5ZcQ0IU\nxOzrSJnMBXZomlZwpsBKptJ+aWqadgC4c3XwTkCUpmmnNE3LANYAg4ENwFCl1CJgS/nWVFRWllxD\nSqnaSqklQEe5SxA5WfhdNAnoCTyqlJpQvjUtPkNFV8BCOR8BgelO4EHNtPLG2IqpkqhiCrqGkoCJ\nFVMlUQUVdB09B3xUMVUqvkp7R1CAqt+zLSqaXEOiNNxV11FVCwRnud0XQPbPZyqoLqJqkmtIlIa7\n6jqqaoHgV8BdKeWmlLIGhiN9AsIycg2J0nBXXUeVNhAopVYDhwAPpVScUmqMpmmZmDpjvgfCgLWa\nphW8Uoeo1uQaEqWhOlxHd8WEMiGEEMVXae8IhBBClA8JBEIIUc1JIBBCiGpOAoEQQlRzEgiEEKKa\nk0AghBDVnAQCIYSo5iQQCCFENSeBQAghqjkJBEKUAqXUA0qpP5VSNtlraIcopdpWdL2EMIekmBCi\nlCilZgO2gB0Qp2na3AqukhBmkUAgRClRSllhykp5HfDR5JdLVBHyaEiI0lMXcABqYLorEKJKkDsC\nIUqJUmoL8DXQAmiYvWyhEJVeVVuzWIhKSSk1EkjXNG2NUkoHHFJKddc0LbiCqyZEkeSOQAghqjnp\nIxBCiGpOAoEQQlRzEgiEEKKak0AghBDVnAQCIYSo5iQQCCFENSeBQAghqjkJBEIIUc39P51tILxl\nAoWnAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives a number of different slopes, from too small to too large, with no principled method for deciding where to apply the cut. A more principled method (derived) below, suggests instead to use all data points $x_i$ above some $x_{\\rm min}$ (in this case $x_{\\rm min}=1$), suppose that's\n", "$n$ data points, then calculate:\n", "$$\\alpha = 1 + {n\\over \\sum_{i=1}^n \\ln {x_i\\over x_{\\rm min}}}$$\n", "with the result:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "alpha=1+n/sum(log(xr))\n", "sigma=(alpha-1)/sqrt(n)\n", "print 'alpha=',round(alpha,3),', with stderr sigma=',round(sigma,3)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "alpha= 2.52 , with stderr sigma= 0.048\n" ] } ], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "How does this more principled method work?\n", "\n", "First consider the normalization constant $C$ in the probability $p(x)=C x^{-\\alpha}$,\n", "determined by\n", "\n", "$$1=\\int_{x_{\\rm min}}^\\infty p(x){\\rm d}x=\\int_{x_{\\rm min}}^\\infty C x^{-\\alpha}\n", "={C\\over 1-\\alpha}x^{-\\alpha+1}\\Big|_{x_{\\rm min}}^\\infty={C\\over \\alpha-1}x_{\\rm min}^{-\\alpha+1}$$\n", "\n", "(where $\\alpha>1$ is assumed so that the integral converges). The result is\n", "\n", "$$C=(\\alpha-1)x_{\\rm min}^{\\alpha-1}\\ .$$\n", "\n", "Now consider the probability of the data.\n", "Given a set of $n$ values $x_i$, the probability $P(x|\\alpha)$ that they were generated by the $p(x)$ distribution (likelihood of the dataset) is the product of the probabilities:\n", "\n", "$$P(x|\\alpha)=p(x_1)\\cdots p(x_n)=\\prod_{i=1}^np(x_i)=\\prod_{i=1}^nC x_i^{-\\alpha} \\ .$$\n", "\n", "The $\\alpha$ that provides the best fit for the data is given by maximizing $P(\\alpha|x)$, in turn given by Bayes theorem as\n", "\n", "$$P(\\alpha|x)={P(x|\\alpha)P(\\alpha)\\over P(x)}\\ .$$\n", "\n", "But $P(x)$ is fixed for a given set of data $x$, and for no prior knowledge of $\\alpha$ it is natural to assume $P(\\alpha)$ if flat (i.e., with no $\\alpha$ dependence). Then $P(\\alpha|x)$ is just proportional to $P(x|\\alpha)$ and $\\alpha$ is chosen to maximize the above likelihood $P(x|\\alpha)$ of the dataset.\n", "\n", "Since the logarithm is an increasing function, it is equivalent (and slightly more convenient) to maximize the log of the function :\n", "\n", "$${\\cal L}(\\alpha)=\\ln P(x|\\alpha)\n", "= \\sum_{i=1}^n (\\ln(C) + \\ln x_i^{-\\alpha})= n\\ln C-\\alpha\\sum_{i=1}^n \\ln x_i\n", "= n\\ln \\bigl[(\\alpha-1)x_{\\rm min}^{\\alpha-1}\\bigr] -\\alpha\\sum_{i=1}^n \\ln x_i\\ .\n", "$$\n", "\n", "Setting to zero the derivative with respect to $\\alpha$,\n", "\n", "$$0={d\\over d\\alpha}{\\cal L}(\\alpha)= n{1\\over \\alpha-1} -\\sum_{i=1}^n \\ln {x_i\\over x_{\\rm min}}$$\n", "\n", "gives the value of the exponent $\\alpha$ that maximizes the likelihood of the data (for the assumed $p(x)=C x^{-\\alpha}$):\n", "\n", "$$\\alpha = 1 + {n\\over \\sum_{i=1}^n \\ln {x_i\\over x_{\\rm min}}}$$\n", "\n", "\n", "With a bit more effort (see, e.g., eqns B7-B12 from appendix B [here](http://arxiv.org/pdf/cond-mat/0412004.pdf)), the expected error (stdev) of the estimation of the exponent $\\alpha$ can be shown to be given by the simple formula\n", "\n", "$$\\sigma={(\\alpha -1)\\over\\sqrt{n}}$$\n", "\n", "(which decreases, as expected, as the inverse square root of the number of data points).\n", "\n", "Finally, here is the data again showing the relation to the best fit power law distribution:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "C=(alpha-1)*xmin**(alpha-1)\n", "\n", "figure(figsize=(6,5))\n", "yh,bh,p=hist(xr,bins=arange(xmin,1000,.1),log=True,histtype='step')\n", "xlabel('x')\n", "ylabel('# occurrences in bin of size .1')\n", "ylim(.5,300),xlim(1,200)\n", "xscale('log')\n", "\n", "xl=arange(1,200,.1)\n", "plot(xl,(.1*n)*C*xl**(-alpha),'k-',linewidth=2,label='$\\\\alpha={:.2f}$'.format(alpha))\n", "\n", "legend();" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFHCAYAAACs30uOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FWX+/vH3J5UqXemCGl1gBVYsCCrZVcpXQRfLWtZ1\nFd0V/S0SqiwohDWgESLoiooNVLDAqqDYV4wVRayIgiAgRZGyFGmBJM/vj0noSU6Sc86ccr+ua66c\nmTPlDh7PJzPzzPOYcw4REYlfCX4HEBERf6kQiIjEORUCEZE4p0IgIhLnVAhEROKcCoGISJxL8jtA\nRZiZ2ryKiFSAc84OXRa1ZwTOuaiZRo0aFVXHqOi+yrNdoOuWtV5F3w/Hf5NI/e8bjuOE4zMU6PqV\nXaei70XiVJKoLQTRJD09PaqOUdF9lWe7QNcta73Kvh8twvV7BOs44fgMBbp+ZdeJlc9Qaay0KhGp\nzMxFY26JHJmZmWRmZvodQ6JctH2OzAwXS5eGRCojHv7Kk9CLlc+RzghEROJESWcEUdlqSERig9lh\n30kSJOX5Y1mFQER8pbP74CtvgdU9AhGROKdCICIS51QIRETinAqBiEicUyEQEQmTp59+mpycHC6/\n/HKeffbZI65z/PHHk5qayjHHHMOTTz5Zrm0rSq2GRETCYNmyZWzatIlBgwaxceNG0tLSOOOMM2jZ\nsuVB6w0bNozu3bvTuHFjkpKSyrVtRemMQEQkDBYtWsTdd98NQP369TnhhBP47LPPDlsvJSWF5s2b\n7ysC5dm2onRGICJSQcuXL+eRRx4p8f2OHTty0UUXAXD++efz2muvAd6zEz///DMnnHDCYdt8+umn\n5OXlsW3bNk488UQuvPDCgLetKHUxISK+KeryoMT3gqWi3xc//fQTU6ZMoX379rz33nvcdNNN1K9f\nn+3bt9OwYcMK55kzZw6PPvoos2bNOuy9F198kd69ewPQvn17cnNzqV27dkDbFivp31WdzomIlMOO\nHTvo3bs3ffv25YILLuDSSy9lwIABvPnmm9StW7fC+92yZQtTp05l2rRpR3y/+AwCoE6dOuTm5ga8\nbUVF7aWhzMxM0tPTY6b3PxE5mN9n/c899xwdOnSgXr16ADRo0IBFixZhZqSkpADluzQE3u901113\n8eijj1KjRg1+/PFHjj322H3vT5s2jZdeeokZM2YAXjEqvldQ1ralyc3NPaigHEqXhkTEN6VdGvLb\n5MmT2bFjBwMHDgS8G7aXXnop3333XYX3ed9999G5c2eaNGnCqlWr2LVrF126dOGHH37guOOO48MP\nP2TPnj384Q9/YOfOnbRp04ZFixZRrVq1Erc9kvJeGlIhEBHfRHIh2LZtG9nZ2XTu3Jm9e/dSvXp1\nHnvsMbp27coVV1xBtWrVyrW/Dz74gC5duuz7fc2MVatW0aRJE0455RQee+wxfve73zF9+nQ2bNjA\njz/+yBVXXMEZZ5xR6rZHokIgIlEjkgtBNNPNYhERKRcVAhGROKdCICIS51QIRETinAqBiEicUyEQ\nEYlzKgQiInEuaruY2LBhAw0aNPA7hohUUjA7l5OKidozgpEjR/odQUQqyTmnKURTeUTtk8UJCQl8\n8cUXtG3b1u84IiJRIeaeLC4sLKRDhwzq1HHUrQt168IDD/idSkQk+kRtIahTpy75+e9w772zWLYM\nLr8cdu/2O5WISPSJ2kLwr3+NBmD06MFUr55H1ao+BxIRiVJRWwj69u1LmzZtWL58ORMnTvQ7johI\n1IraQpCUlMSECRMAyMrKYseOdT4nEhGJTlFbCAC6du1Kr1692L59Ox9+ONzvOCIiUSmqCwFATk4O\nycnJLFo0ldWrP/M7johI1In6QpCWlsYtt9wCOGbN6q/RjkREyilqHyg7MPfWrVtp1CiNXbs28Oyz\nz3L55Zf7mE5EJDLF3ANlB6pVqxadO48BYOjQoezcudPnRCIi0SMmCgFAmzZ9aNy4HatWrWL8+PF+\nxxERiRoRVwjM7CIze9jMnjWzroFul5CQyB//eC8A2dnZrFmzJmQZRURiScQVAufcbOfc34G+QLku\n9h9/fBcuueQSdu7cybBhw0ITUEQkxoSlEJjZ42b2i5ktPGR5DzNbbGZLzezWQza7Dbi/vMcaN24c\nqampTJ8+nY8//rgysUVE4kK4zgimAD0OXGBmiXhf9D2A1sCVZtbKPNnAa865L8tzkJ07oVatltx8\n8yAA/t//68/GjYVB+QVERGJV2JqPmlkL4GXn3MlF82cCo5xzPYrmi6/l7AD+CnwKfOmcm3yEfblD\nc48YAQ8+6L12bjtbt56Icz9TpcoT7Np1TUh+JxGRaFJS81E/h6psAqw+YH4NcIZzrh/w77I2zszM\n3Pc6PT2dMWPSGTOmeEkNnnjiTq699lry8oaxffvF1KhRI3jJRUSiQG5uLrm5uWWu5+cZwSVAD+fc\n34rmr2Z/IShrX4edERyqsLCQ007ryOeff8qIESPIysqq7K8gIhLVIvGBsrVAswPmm+GdFQRFQkIC\n2dlec9Lx48ezcuXKYO1aRCSm+FkIFgBpZtbCzFLwmoq+FMwDnH76mSQlXUVeXh5Dhw4N5q5FRGJG\nuJqPPgN8BJxoZqvN7DrnXD7wD+AN4FvgOefcd4HuMzMzM6BrX6mpd1G1alVmzpzJu+++W7FfQEQk\niuXm5h50X/VQMdHpXEm2bYOmTWHQoNFkZmbSvn17FixYQGJiYhhSiohElki8RxA2Q4YMoVmzZnz5\n5ZdMmTLF7zgiIhElLgpBtWrVuPvuuwEYMWIE27Zt8zmRiEjkiPlLQ7VqQVISOOcoKDgb+BCzIRQU\n3I0ddoIkIhK74vLSUM2asGeP1/XErl3GvHkTAXBuIkuXLvU5nYhIZIjaQhBIqyEzSE7eP3XseCrX\nXnstsJchQwaHI6aIiO/iutXQkfz88880bnwisJ233nqL8847L7jhREQiVFxeGjqSRo0aAcMByMjI\nID8/399AIiI+i7szAm/73bRs2ZoVK1YwadIkeve+mfnzD1/vggu8G80iIrGgpDOCOC0E8J//vMCl\nl15CvXr1mDDhe4YMqcsZZ+xfZ84c2LoV1GmpiMSKSOyGulIyMzNJT08nPT29Qtv37t2b9PR0cnNz\nee650XTseC+zZu1/XwVARGJFWd1Rx+0ZQWEhfP31V5xyyimAkZ6+kLffbrVvnRo1YN06FQQRiR1B\nvVlsZr+pfCT/tWvXjhtuuIHCwgK++Wag33FERHxR0VZDbwY1hQ8WLfKmq6/OomrVWqxf/zqvvvqq\n37FERMKuxEtDZlbacJHXOudqhiZS2Sp7aahNm4PnN226h19+GcRJJ53E119/TUpKii4NiUjMKXer\nITP7FRgM5AEHrmRAjnOuXiiCBqKyheBQe/bs4eSTT+b777/nnnvuYcCAASoEIhJzKlII3gFuc859\neIT3VjrnWgQ9ZYCCXQgA5syZQ69evahVqxZLly6lZcsGKgQiElMqcrP4EuCLI73hZxEoFugIZYG6\n4IIL6NatG1u3bmXkyJFB26+IiN/U11A5fPvtt7Rt2xbnHKmpX7B+fVudEYhIzAh289HRlY8UeVq3\nbs3NN99MYWEheXkZRGORFBEpr4o2H10Q1BQRJDMzk7p161JY+A4vvzyr7A1ERKJchQqBc+7lYAeJ\nFHXr1mX0aO+E57bbBpOXl+dzIhGR0Iq7bqgD0bdvX8xas2LFciZOnOh3HBGRkNLN4hJUrfoWu3d3\no0aNGixdupSGDRuG9HgiIqGmgWnKKTGxK//3f73Yvn07I0aM8DuOiEjIlFkIzKyhmT1mZq8Xzbc2\ns+tDH610wX6O4EjGjs0hOTmZKVOm8Nlnn4X0WCIioVLp5wiKCsAUYIRzrq2ZJQNfOOd+G8yg5RGO\nS0PFXUxkZg4mJyeHzp078/7772N22FmViEhUqMylofrOueeAAgDn3F4gbgb6vf3222nQoAEffvgh\nM2bM8DuOiEjQBVIItpvZvg7mzKwjsDV0kSJLrVq1GDNmDABDhw5l165dPicSEQmuQArBIOBl4Dgz\n+wh4CrglpKkiTJ8+fWjXrh2rVq1i/PjxfscREQmqgJqPmlkS8Bu8LqiXOOf2hDpYGXnCdo+guK+h\nd999l/T0dKpVq8aSJUto2rRpSI8vIhJs5e6G+oANlwPjnHMPHrBsjnOuZ/BjBiZchaB5c0hM3L+s\nefNLefXV57n88j+zaNG0g9b/859h2LCQRhIRqZTKFIIlwJfATqCvcy7PzL5wzv0uNFHLFo5CsGiR\nN8B9se7doWbNFSxd2grn8khJmceCBR0BmDYNdu+Ge+8NaSQRkUqpTKuhnc65y4HvgPfM7Nigp4tA\nbdrAySfvn15/HaZPb8lf/+oNcp+W1p82bQo5+WRo3NjnsCIilRDwk8XOubuBEXgD18fdBfK2beHU\nU+G++/5Jo0aNWLRoPtOnT/c7lohIpQVSCPYN1+Wc+y/QDShtYPuwCMeTxUdSs2ZN7rzzTgCGDRvG\n9u3bw55BRKQ8KvxksZm1cs59Z2YdOHjw+uLtfOtzIRz3CEpTWFhIx44d+fTTTxkxYgQNGmSxfLnu\nEYhIZKvI4PWPOOf+Zma5HF4IcM79PugpA+R3IQD46KOP6Ny5M6mpqQwdupitW1uoEIhIRKtwq6FI\nFAmFAOCqq67imWeeoX37yzjnnBkqBCIS0SrcasjMLjOzo4pe325mL5jZKaEIGW2ys7OpWrUqX345\nk7Vr3/M7johIhQR0s9g5t83MzgLOBR4HHgptrOjQrFkzbr31VgDeey+DgoICnxOJiJRfIIWg+Nut\nJ/CIc24OkBy6SNFlyJAh1K7djA0bvmDKlCl+xxERKbdACsFaM3sYuBx4xcyqBLhdXKhWrRq9emUD\nMGLECLZt2+ZzIhGR8gnkC/1PwBtAN+fcFqAOMCSkqaLMKadcQaNGnVi/fj1ZWVl+xxERKRe1GgqC\ne++FefMW8Nxzp5GcnMy3337LCSec4HcsEZGDaPD6EDvmmFO59tpr2bt3L4MHD/Y7johIwEosBEX3\nAqQcxo4dS40aNZg9ezbTpv2XlSth5UqvZ1IRkUhV2hnBRwBmNq2UdXzjV19DpWnUqBHDhg0HoE+f\nDLp0yeekk+CLL3wOJiJxrTJ9DS0CxgJ3AIPxRicr5pxzLwQvZvlE4j2C4r6Gtm/fTc2arYEVTJo0\niaeeupl77oEzz/Q7pYjEu4rcI+gLnA3UAnrhPUdQPPUKRchYUKVKFRISvHGNR44cSX7+Zp8TiYiU\nLqmkN5xz7wPvm9kC59yjYcwUA3qTnp5Obm4uycmjgYl+BxIRKVGJheAAT5pZf+Ccovlc4CHn3N6Q\npYpyZsbEiRM55ZRTWLfuflauvJEzz2zFk0/Cd98deZtRo6CKbs+LiA8CaT76IHAKMAl4AOhQtExK\n0a5dO2644QaggPvu84a3fP55WLMGjjrq4Omee2DPHn/zikj8CmTw+q+dc23LWhZOkXiz+LbboE4d\ncA5+/hny82H9+vU0bpxGQcE2XnnlFSZPPp8+feCiiw7e/qij9hcIEZFQqcwDZflmtu8xWTM7HsgP\nZrhod911sGgRfPABfPih14II4Oijj6Zp01EADBw4kMJCXU0TkcgTyD2CIcBcM1tRNN8CuC5kiaJQ\n8SWeIznmmH/g3GSWLFlCUtIkICOs2UREylLmGYFz7m3gROAWoB9wknNubqiDxYqEhBT69csBYMmS\nTLZu3eBzIhGRgwXU15Bzbrdz7ivn3NfOOXWYUE6dOl1At27dyM/fyjPPjPQ7jojIQdTpXBiYGRMm\nTMAskTfffJivv/7a70giIvuoEIRJ69atadHiJgoLCxkwYACR1OpJROJbQIXAzJqYWWczO8fMupjZ\nOWVvJYc66aTR1KxZl7lz5zJ79my/44iIAAEUAjPLBj4ERuC1IBqMRiirkJSUulxxxWgABg0aRF5e\nns+JREQCOyPojddS6HznXK/iKdTBYlWPHn1p3bo1y5cvZ+LE/X0QDR0K/fvDv//tYzgRiUuBFIIf\ngJRQB4kXiYlJTJgwAYCsrCzWrVtHdja0agUFBfCCb517i0i8CuSBsl3Al2b2NlB8LcM5524JXazY\n1q1bN3r27MmcOXMYMWIEjz32GADvvOM9oSwiEk6BnBG8hDc4zUfAZwdMvorEEcrKIycnh+TkZKZM\nmcLnn3/udxwRiWEVHqEskkVap3OlOfNM+NvfoG1bGDgQBg3a3+nc4MGDycnJ4ayzzuK9994jN9f4\n17+8MwMRkWArd6dzZjaz6OfCI0x6IipAv/kNPPAA9O0LO3dC7dr737v99ttp0KABH3zwATNmzPAv\npIjEtdLGLG7snPvJzFoc6X3n3MrQxSpdNJ0RlOXhhx/mxhtvpHnz5kyevJjs7Ko6IxCRkCj3GYFz\n7qeinyuB3UA74GRgt59FINZcf/31tGvXjlWrVjFjxni/44hIHArkgbIbgPnAxcClwCdmdn2og8WL\nxMTEfc8TPPPMXeTlrfU5kYjEm0BGKPseONM5t6lovh4wzzl3YhjylZQpZi4NFbv00kt5/vnnOeaY\nq1m37im/44hIDKrMCGUbge0HzG8vWiZBNG7cOJKTU/nll2l8/PHHfscRkThSWquhQWY2CFiGdzko\n08wygY+BpWHKFzdatmzJZZd5g9z379+fwsJCnxOJSLwo7YygJlADr4uJWYArmmYDy0MfLf5cddU/\nSUlpyPz585k+fbrfcUQkTuiBsgjyzjtw881PsHjxtTRu3JglS5ZQo0YNv2OJSIyozD0CCaNjjvkL\np512Gj/99BPZ2dl+xxGROKBCEGHMEvY1J83KGkfDhivp0cPnUCIS01QIIlCnTp1o2vRKII+0tKGs\nX+93IhGJZWV2Q21mRwN/A1ocsL5zzvUJYa6499vfZrN+/Sw++GAmaWnvARodVERCI5AzgtnAUcBb\nwCsHTBJC1ao1449/vBWANWsyKCgoKHHdH36A5s29qX37cCUUkVgRSCGo6py71Tk3wzn3n6Lp+ZAn\nEy66aAjHHNOMXbu+YOrUqSWul58PSUkwezasVQ8VIlJOgRSCOWZ2QciTyGFSU6txyy1ey6Hhw4ez\nbdu2EtdNToamTcOVTERiSSCFIAN42cx2m9mvRVPJ30gSVN27X0H16p1Yv349WVlZfscRkRhUZiFw\nztVwziU456o452oWTUeFI5x4D4A0aeI1J504cSLLli3zOZGIxJrS+hpqVfTzlCNN4Yso1aufxrXX\nXsvevXsZPHiw33FEJMaU1nx0IF6z0Xvw+hg61O9DkkiOaOzYscyc+R9mz57NY4/9l4suOo/69Q9f\nLz8fvvtu/3xyMpxwQvhyikj0KbEQOOf+VvQzPWxppESNGjWiZs3h7NgxnBtvzGDNmi8ZNerg/3yJ\nidCwIVx8sTe/Zw84B8vVRaCIlEJPFkeRBg0G0KRJSwoKFrFgwcOHvV+3rnc2UDy99ZYPIUUk6qgQ\nRJGEhCpkZIwDYO7ckWzevNnnRCISCyKqEJhZSzN71Mxm+p0lUp177sU0b96FnTs3MWTIaJ56ynuQ\nTESkogIZvP4sM6tR9PovZnaPmR0bijDOuRXOuRtCse9YYWZ07ToRswSmTp3EoEHfsXAhXHih38lE\nJFoFckbwILDDzNrhtST6AXgy0AOY2eNm9ouZLTxkeQ8zW2xmS83s1nKljnMNG7anQ4cbKCjIBwby\n1FMwbpzfqUQkWgVSCPKLhgP7IzDJOTcJbxjLQE0BDupR38wSgfuLlrcGrix+bkECc+65d1Ct2lFs\n2PA6r776qt9xRCSKBVIIfjWz4cDVeP0OJQLJgR7AOfc+cOhdzdOBZc65lc65vcCzwEVmVtfMHgLa\n6yxhvx074NNPYefO/ctq1DiaK64YCcDAgQPZu3evT+lEJNqVOR4BcDlwJdDHObfOzJoDlb0Q0QRY\nfcD8GuAM59z/gL6B7CAzM3Pf6/T0dNLT0ysZKTJVqwZHHQU33wy1akHVqvvf69mzH889N5klS5Yw\nadIkMjIy/AsqIhEnNzeX3NzcMtcrsxA45342sxeA4udTNwKzKpXuyE8ql8uBhSCW/eY33tnAkSQn\np9Cq1T0sWNCL0aNHc/XVV1P/SI8bi0hcOvSP5NGjRx9xvUBaDf0dmAlMLlrUFHixkvnWAs0OmG+G\nd1Yg5XT00RfQrVs3tmzZwsiRI/2OIyJRKJB7BP8POAvYBuCc+x44upLHXQCkmVkLM0vBu/z0UiX3\nGZfMjAkTJpCYmMjkyZP5+uuv/Y4kIlEmkEKQ55zLK54xsyTKcWnHzJ4BPgJONLPVZnadcy4f+Afw\nBvAt8Jxz7rvS9hMvFi+G66+Hzz4rfb0XX4R77/Vet27dmptuuonCwkIGDBiA18hLRCQwVtaXhpmN\nA7YA1+B9ed8MfOucGxH6eCVmcqNGjYq5m8Q//QSvvbZ//g9/gJYtD19vwQL46ivvdbNm0K0bbNq0\nibS0NDZv3syLL77IH//4R5Yvh/POU6dzIvGu+Kbx6NGjcc7Zoe8HUggSgeuBbkWL3gAedT7+2Wlm\nfh4+Yt1///3069eP448/nkWLFrF2baoKgYjsY2YVLgTVgd3OuYKi+UQg1Tm3s9QNQ0iF4Mjy8/Np\n164d3377LdnZ2Vx66VAVAhHZp6RCEMg9grnAAa3XqQb8N1jBJHiSkpKYMGECAHfccQcbNqzzOZGI\nRINACkGqc2578Yxz7le8YiARqFu3bvTs2ZPt27eTkzOCzZth7FiYNs3vZCISqQIpBDvMrEPxjJmd\nCuwKXaTAZGZmBvTEXDzKyckhOTmZ//xnCr17f86iRfDII36nEhG/5ObmlvoQbiD3CE7D6wvo56JF\njYDLnXMLgpSx3HSPoGyDBw8mJyeHs846izvueI9Ro4x33/U7lYj4qcI3i4s2TgFOwnt+YElRR3G+\nUSEo29atW0lLS2PDhg2MGvUc77zzJxUCkThXmZvFAKcCbYEOeF1GXxPMcBJ8tWrVIisrC4CHHhpC\nQYHvV/NEJEIFcmloGnAc8CVQULzcOdcvtNFKzaQzggAUFBTQoUMHvvrqK4499l+sXHm735FExEeV\neY7gO6B1JH3zqhAELjc3l9///vckJFRj1arvadKkSZnbvPkmTJnive7d2+vOArxlVaqEMKyIhFRl\nLg19g3eDOKKo1VBg0tPTOeecSygs3MmwYcMC2mbZMti4EV5/HRYuhPXr4YUXID8/xGFFJCSC0Woo\nF2gPzAeKO59zzjnfhkvXGUH5PPvsCq66qhXO5TFv3jw6duxY6voPPADffAOzZ0OfPrBhg/ccwrp1\nUKNGmEKLSNCVdEYQyAhlmUU/HWAHvJYo0ahRS5o2Hcjq1XfSv39/5s2bR0JCoO0ERCTWlflt4JzL\nBVYCyUWv5wNfhDSVBF2zZv+kYcOGzJ8/n+nTp/sdR0QiiF8jlEmYJSXV5M477wRg2LBhbN++vYwt\nRCRe+DVCmfjgmmuu4dRTT+Wnn34iOzvb7zgiEiFCPkKZRI6EhATuLRrWLCtrPLVqrWTOHKhb15se\nfjjwfXXt6m1z+yGPJjzxhLe8S5cgBheRkAqkELxrZiOAambWFe8y0cuhjVU2NR+tmE6dOtGixZXA\nbvLyhrJ3Lxx3HFxyCeTllbn5Ptu2wdlnw65DHljOy4Njj/XeF5HIUFbz0UAKwa3ABmAhcCPwKnBb\nMMJVRmZmZkwNUxlOp5ySTUpKVfLyZrJo0XskJVXsQbGUlCMvTwqkLZqIhE16enrFC0HRZaBvnXMP\nO+cuLZoeUSP+6Fa9ejPOP38oAI8/nkHR4HMiEqdKLQTOuXxgiZkdG6Y8Eibnnz+UhISmrFjxBRs2\nTPU7joj4KJBLQ3WBRWY218xeLppeCnUwCa3U1GrUrHk3AKtXD2fPHl3UF4lXgVzNvY39TxQX06Wh\nGFClyhU0bHg/S5Z8xKefZvHb397tdyQR8UEg9wgeds7lHjJpiJMYYGb06TMRgK++msiGDct8TiQi\nfgjkHsFi3SOIfu+957Xy6XfIKBJpaadRv/5fKSzcy6xZgw/bbsyYwPaflubt/8kngxBWRMIqau8R\n6DmCwJ19tte+f/x4KDhCA6FmzcaSnFydhQtn8/bbb+9bvmIF7N4N999f9jH27oWOHaGwMIjBRSQo\nynqOIJB7BBE5rFVpv5QcLCHB+2u9pPb9KSmN6dBhBB9/PJyMjAz+/vcvgKQSnxMo7TgiEnnS09NJ\nT09n9OjRR3y/zEJQ1OOoxLj27QewdOkjfPPNN3z44cPUrXuz35FEJEwC6X10u5n9WjTlmVmhmamt\nYYxJSqrCRReNA2DOnJHk5W32OZGIhEsg4xHUcM7VdM7VBKoCFwMPhDyZhF27dhfTpUsXduzYxKef\nHvkUUkRiT7mu6jrnCp1zs4AeIcojPjIzJk6ciJmxaNEkFi9e7HckEQmDQC4NXXLAdJmZ3QXsKms7\niU7t27enU6e/UViYz8CBA/2OIyJhEMgZQS+gZ9HUDfgVuCiUocRfvXrdQUrKUbz22mu8+uqrfscR\nkRALpNXQtWHIIRGkZs2j6dBhJPPmDWbgwIF07doVSOb996FaNahadf+6O3fCu+96zxuISHQK5NLQ\nE2ZW+4D5Omb2eGhjlU0PlIXWySf3Iy0tjSVLljBp0iTOOgvuugsGD4Y//3n/eqtXQ48ecNJJkJrq\nX14RKVkwBqZp55zbUjzjnNsMnFL5aJWjgWlCKzExhZycHABGjx7NtGkbefddePrpw9dt3tw7K2jU\nKMwhRSQglRqYpoiZWd0DZuoCiZWPJpGuZ8+edOvWjS1btjBy5Ei/44hIiARSCHKAeWZ2h5llAfOA\ncaGNJZHAzJgwYQKJiYlMnjyZhQsX+h1JREIgkAfKnsR7iGw9sA7oXbRM4kDr1q256aabKCwsJCMj\nA41SKhJ7ArlZ3BFY7Zz7t3PufmCNmZ0R+mgSKTIzM6lTpw5z587lv/+d7XccEQmyQC4NPYT37ECx\nHUXLJE7Uq1dvX6+F2dmDcS7P50QiEkwBdTHhDrge4JwrQDeL407fvn1p1aoVq1b9wJYt9/odR0SC\nKJBCsMIx9GeQAAAQBElEQVTMbjGzZDNLMbP+wPJQB5PIkpyczIQJEwD43/+y2LNnnc+JRCRYAikE\nfYHOwFpgDdAR+HsoQ0lk6t69O7//fU+c+5U1a0b4HUdEgiSQVkO/OOcud84dXTRd6ZxbH45wEnmG\nDcsBktm4cQqbN3/udxwRCYJAWg01M7MXzWxD0fS8mTUNRziJPC1bnkjt2v0Ax5df9ldzUpEYEMil\noSnAS0DjounlomUSA8aOPXj+4ovhwQdL36ZOndtJSqrPpk0f8P77M8t9zNxcuPBC71gXXgi3lzAq\n9n33ee9/9lm5D1GqbdvgvPPgscfKt91zz3nb7d0b3Dwifgtk8PoGzrkDv/inmtmAUAUKVHFfQ+pv\nqOKmToW8PKhTB2bM8JbNmQPPPANpaSVvl5hYm6ZNx7By5Y2kpg5hwoReeIPXBWbtWnjzTXDO+1Ld\nVcLoFl99BS+/DDcHefjkvXvh7bfh5JPLt93Spd52OgmSaJObm1tqJ52BnBFsMrO/mFmimSWZ2dXA\nxmAFrCh1Old53bt7f3GfffbBy3v1grZtS9+2QYPrSUtrx4YNq/jss/HlPnZS0Z8gZuXeNGj8PLZI\nOAWj07k+wJ/wupf4GbgMuC4Y4SR6mSWSkTERgLvuuou1a9f6nEhEKiqQVkMrnXO9nHMNiqaLnHOr\nwhFOIluHDulcfPHF7Ny5k2HDhvkdR0QqqFyD14scaty4caSkpDBt2jQ2bPjY7zgiUgEqBFIpxx13\nHIMGDQJg/vwMnCv0OZGIlJcKgVTaP//5Txo2bMjGjZ+wceN0v+OISDkF8kDZbQe8rhLaOBKNatas\nyZ133gnAjz8Oo6Bgu8+JRKQ8SiwEZjbMzDrhtRIq9lHoI0k0uuaaa6hX71T27v2Jdeuy/Y4jIuVQ\n2hnBYrwi0NLMPjCzR4D6Zvab8ESTaJKQkMDpp3vNSdetG8+ePT/6nEhEAlVaIdgC/BP4AUgH7gMc\ncKuZzQt9NIk2Rx/dmXr1rsS53fz001C/44hIgEorBN2BV4Dj8QawPx3Y6Zy7zjl3ZjjCSfQ59ths\nEhKqsmXLDN577z2/44hIAEosBM65fzrnzgVWAE/h9UtU38w+NLOXwxVQoktqajOOOcY7G8jIyKCg\noMDnRCJSlkCaj77hnFvgnJsMrHHOdcbrdkLkiBo2HEpyclO++OILpk6d6nccESlDIF1MHHix99qi\nZRtCFUiiX2JiNRo39loODR8+nG3btvmcSERKU64HypxzX4UqiMSW2rWvpFOnTqxfv54xY8b4HUdE\nSqEni+NMYaE3hZ6Rk+M1J504cSLLli2jsBAKCgI7fkk5i/eRn+/9DMbYAAUFpe+r+JgVUZltRcIl\nagtBZmZmqQMtyOHM4JFHvIFnQt0X/9dfQ+fOp/HXv/6VPXv2MHjwYO64wxuHoF+/sre/7TZv4JxD\nPfigt4/kZO9nXl7lsxbv6/MShmDu33//+Anl1bMn1KhR8WwiwZCbm1vp8QgikgamKb+bbtr/1+9f\n/hK64/zud/uHcxw7dizVq1dn9uzZLF/+NvXqBb6fMWO8wXPCoX37kt+rzF/0a9bA7t0V314kGIIx\nMI1IhTVu3Jjhw4cD8PrrGTiX73MiETmUCoGE3MCBA2nRogXr139DXt4jfscRkUOoEEjIValShfHj\nvXGNd+68Hec2+5xIRA6kQiBhcfHFF3PssV1wbhO7d4/2O46IHECFQMLCzOjRYyJg7NkzicLCxX5H\nEpEiKgQSNo0atSc19QYgnz17BvodR0SKqBBIWFWrlgUcRUHBaxQWvuZ3HBFBhUDCLCHhaKpUGQlA\nfv4AnNvrcyIRUSGQsEtJ6YdZGrAEmOR3HJG4p0IgYWeWQkpKDgDOjWbPno0+JxKJbyoE4ovExJ6Y\ndQW2sHz5SL/jiMQ1FQLxhZmRmDgBSGTNmsksXLjQ70gicUuFQHyTkNAGuAkoJCMjAxeMPqVFpNxU\nCMRXZpkkJdVh7ty5zJ492+84InFJhUB8ZVaP447zupwYPHgwecEYYEBEykWFQHzXtGlfWrVqxQ8/\n/MC9997rdxyRuKNCIL5LSEhmwoQJAGRlZfHLL7/4nEgkvqgQSETo3r07PXv25Ndff2XEiBF+xxGJ\nKyoEEjFycnJITk7m8ccfZ926EgYQFpGgUyGQiHHiiSfSr18/nHO89VZ/NScVCRMVAokot99+O/Xr\n12fNmg/45puZfscRiQsqBBJRateuTVZWFgCvvz6EgoJdPicSiX0qBBJxbrjhBo4+ui1btqxi5coc\nv+OIxLyIKgRmVt3MnjCzh83sKr/ziD8SExM57zzveYIVK+5k48a1PicSiW0RVQiAi4EZzrm/Axf6\nHUb8c+yx6bRpczEFBTuZMmWY33FEYlrIC4GZPW5mv5jZwkOW9zCzxWa21MxuLVrcBFhd9Log1Nkk\nsvXoMQ6zFObOncbHH3/sdxyRmBWOM4IpQI8DF5hZInB/0fLWwJVm1gpYAzQLYzaJYHXrHkeLFt4g\n9xkZGRQWFvqcSCQ2hfzL1jn3PrD5kMWnA8uccyudN2jts8BFwAvAJWb2APBSqLNJ5DvuuOHUqdOQ\nTz75hKefftrvOCIxya+/ug+8BATemUAT59xO51wf59zNzrlnfMoWt3JyYG8AY8mvXg2LFgW2z8su\ng44dvemRRw5+b88e7+fHH3vvDxp0+PZJSTVJTb0TgD59biU7e/tB759zDowdC3/6E6SkePsp7s36\n5pvh5JO91yedBGbe+yMDHBDtH/+A88/fP//QQwe/v2MHHH88/Oc/h297/fXe8Y+kd28YMwaWLfO2\nX7zYWz56NFx9tfd66VLvvdWrD972ssugdm246abD93v77dCp0/75J57wfud//av037PYoEHev2V5\nffkl/O53++f/9Cd4/vnAt3/rLejVq/zHLUtGBkydGvz9xqIkn45b6UdGMzMz971OT08nPT29sruM\nawMHwp//7L1OKuVT0bQp5OZ6r1u3Lnu/X30FWVnQvLk3P2uW94V69tkwahTMnw/9+8Mnn0D9+odv\nP2ECbN58DTfcMInvvlvAuedm8/PPd/DGG957v/wCy5fD3LleEfvkE1i/3tt2/nz45hvv9fff7//d\nli0rOzfAq6/CihX7v3TPP//gYlZQ4B37SH3kPf649yX8298e/t6sWbBpE1x4obf9zp3e8uefh4UL\nYdo0b9ny5bB798HbFhedhx6CBx88+L1PP4V58/bPf/CB93PSpMCK3z33QMOGMHx42eseaMMGrxgU\nmzkTjjkGLrkksO0/+wzmzCnfMQNx773e5+zaa4O/72iRm5tLbvH/sKXwqxCsZf+9AIperynPDg4s\nBFJ5LVt6U1mqVvX+qi6P9u3hxBO915995v2sW9fbz6+/lr5tq1YACTzyyETOOussnn12PJdccgNw\nLGec4Z2ZbNhw5DOZ1NTDl9WtC4EOeZCcfPB8s2aHLytNSkrJ75kFtqw8Eko4v6/sfiV6HfpH8ujR\no4+4nl+XhhYAaWbWwsxSgMvRPQEpRefOnbniiivYvXs3H3881O84IjElHM1HnwE+Ak40s9Vmdp1z\nLh/4B/AG8C3wnHPuu1BnkeiWnZ1N1apV+eGHGcD7fscRiRnhaDV0pXOusXMu1TnXzDk3pWj5a865\nk5xzJzjn7izvfjMzMwO69iWxo3nz5gwdWnw20J+CAj1qIhKI3NzcUi+nR21b/czMTN0gjkNDhw6l\nevWmwBe8+upUv+OIRIX09PTYLAQSn6pVq0bHjtkAPPTQcHbt2uZzIpHop0IgUeeEE64EzmTz5vW8\n8soYv+OIRD0VAok6ZgZ4vZO+/fZEtm37wd9AIlEuaguBbhbHu9M4//y/kp+/h/nzB/sdRiSi6Wax\nxKwbbxxLamp1fvxxFvn5b/sdRyRi6WaxxKwGDRrzf//n9Yewe3cGkO9vIJEopUIgUa1bt4HUqNGC\nwsJvgEfKXF9EDqdCIFEtObkKp58+rmjudg7v8VxEyqJCIFGvRYtLSEzsAmwCAuxzWUT2idpCoFZD\nUszMqFJlImDA/fz882K/I4lEFLUakriQmNgeuAHIZ+bMgX7HEYkoajUkcSQLOIpvvnmN1157ze8w\nIlFDhUBiyNF4N4xh4MCBFBYGMO6miKgQSKy5haOPTmPx4sWsW/eA32FEooIKgcSYFC67LAeA1asz\ngY2+phGJBioEEnPatu1J165dKSjYAgQwartInIvaQqDmo1ISM2PChAlAIjCZhQsX+h1JxFdqPipx\nqU2bNjRs2BcoZMCAAYDzO5KIb9R8VOJWs2ajgTq8/fbbwEt+xxGJWCoEErOSk+sBo4vmBgF5PqYR\niVwqBBLj+tKqVSvgB4pHNfPk+hNHYkqs3KdUIZAYl1x04xggi7y8X4pe5/qUR2KJCoEELBwflmAe\no6L7Wrcu8O1++SWwdcvKsmRJWfvJpXv37sAFwK8sXjwioONGmnB94eTlBec4a9ZUbD/l/T0DWb+y\n68TKl31pVAjCQIXgcOEsBJ57gCRWrXqcbds+D+jYkSRcX0Z79gTnOGvXVmw/KgT+MOeir1mdmUVf\naBGRCOCcs0OXRWUhEBGR4NGlIRGROKdCICIS51QIRETinAqBiEici4lCYGbVzewJM3vYzK7yO49E\nHzNraWaPmtlMv7NI9DKzi4q+h541s65+5wlUTLQaMrO/AP9zzr1iZs86567wO5NEJzOb6Zy7zO8c\nEt3MrDYw3jl3g99ZAhGxZwRm9riZ/WJmCw9Z3sPMFpvZUjO7tWhxE2B10euCsAaViFXOz5DIEVXw\nc3QbcH/4UlZOxBYCYArQ48AFZpaI94/bA2gNXGlmrYA1QLOi1SL5d5LwKs9nSKQkAX+OzJMNvOac\n+zL8USsmYr80nXPvA5sPWXw6sMw5t9I5txd4FrgIeAG4xMweQB3PS5HyfIbMrK6ZPQS011mCHKic\n30X/AM4FLjWzG8ObtOKS/A5QTgdeAgLvTOAM59xOoI8/kSTKlPQZ+h/Q159IEoVK+hz1A/7tT6SK\ni9gzghJE/51t8Zs+QxIMMfU5irZCsJb99wIoer3GpywSnfQZkmCIqc9RtBWCBUCambUwsxTgcnRP\nQMpHnyEJhpj6HEVsITCzZ4CPgBPNbLWZXeecy8e7GfMG8C3wnHPuOz9zSuTSZ0iCIR4+RzHxQJmI\niFRcxJ4RiIhIeKgQiIjEORUCEZE4p0IgIhLnVAhEROKcCoGISJxTIRARiXMqBCIicU6FQEQkzqkQ\niASBmZ1mZl+ZWWrRGNrfmFlrv3OJBEJdTIgEiZndAVQBqgKrnXPZPkcSCYgKgUiQmFkyXq+Uu4Az\nnf7nkiihS0MiwVMfqA7UwDsrEIkKOiMQCRIzewl4GjgOaFQ0bKFIxIu2MYtFIpKZXQPkOeeeNbME\n4CMzS3fO5focTaRMOiMQEYlzukcgIhLnVAhEROKcCoGISJxTIRARiXMqBCIicU6FQEQkzqkQiIjE\nORUCEZE49/8BtN3SqhvM3lMAAAAASUVORK5CYII=\n", "text": [ "" ] } ], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, here is the same thing with `n=` 1 million data points:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "alpha=2.5\n", "xmin=1\n", "n=1000000\n", "xr=xmin*rand(n)**(-1/(alpha-1))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "alpha=1+n/sum(log(xr))\n", "sigma=(alpha-1)/sqrt(n)\n", "print 'alpha=',round(alpha,4),', with stderr sigma=',round(sigma,4)\n", "C=(alpha-1)*xmin**(alpha-1)\n", "\n", "figure(figsize=(6,5))\n", "yh,bh,p=hist(xr,bins=arange(xmin,20000,.1),log=True,histtype='step')\n", "xlabel('x')\n", "ylabel('# occurrences in bin of size .1')\n", "ylim(.5,150000),xlim(1,20000)\n", "xscale('log')\n", "xl=arange(1,200,.1)\n", "plot(xl,(.1*n)*C*xl**(-alpha),'r-',linewidth=2,label='$\\\\alpha={:.4f}$'.format(alpha))\n", "title('n={:,} data points'.format(n))\n", "legend();" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "alpha= 2.5004 , with stderr sigma= 0.0015\n" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFWCAYAAABkVZqwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVPX1//HXYWlKkWoBVJQiGpBiIxF1ExvKV40KKnbs\nDWxRMVGX1V+iYFSsKIJoNBZUTDQWbKxgiahRBBUFEaRXkV73/P64szAsO7szOzN7Z3bez8djHju3\nzL1n787Omfv53Ps55u6IiEjuqhF2ACIiEi4lAhGRHKdEICKS45QIRERynBKBiEiOUyIQEclxSgQi\nITKzfDObHXYcFTGzs8xsbNhxSHooEUiVMbNaZvaSmf1kZsVmdkQF6zcxs1fMbJWZzTSzvqWWH2lm\nU81stZm9b2Z7lFo+2MyWRB53VbCvSm/LzFqb2bjIa78zsyPjOyKJMbPzzWxCOrZdEXf/p7sfG8+6\nYcYplaNEIFVtPHA2sACo6G7Gh4F1wM7AWcAwM9sPwMyaAS8DfwEaA58DL5S80MwuBU4C9o88TojM\n204KtvUc8AXQJLKNlyLbFMkO7q6HHhU+gJnA9cAkYDnwPFAnie3NBg4vZ3k9YD3QNmreU8CdkeeX\nAB9GLdsRWAO0j0x/DFwUtbwf8EmMfVV6W0B7gmRVL2r5B8ClMfa1A/AksAz4BrgBmB21fCAwHVgR\nWf7HyPx9gbXAJmAlsCwyvxfwJfAr8DNQUM4xzQfmADcDi4GfgDOjlu8E/ANYFPl7/wWwyLLzgQlR\n6xYDlwI/AL8AD1UQ5/GR32dFJIbrw35P67H1oTMCiZcDfYBjgb0IvhmfD2Bme5jZL+U8zqjE/toD\nm9x9etS8ScBvIs9/E5kOgnNfQ/ABWrJ8v+jlwNdRy0pLZlu/AWa4++oYcZZWQHD89iY4luex7ZnR\ndKCHuzcECoFnzGwXd/8OuIwgATVw9yaR9VcBZ7v7TgRJ4XIzOynGvgF2AZoCLSL7Hm5m7SPLHgQa\nROI7AjiXIOnF0gs4kOC9cJqZHVtOnCOBSyK/12+A98vZrlQxJQJJxAPuvsDdfwFeA7oAuPvP7t64\nnMfzldhXfYJvj9FWEnxQxVq+otTyX0stqx9jX/WS2FbpZaVfW1of4K/uvtzd5wD3A1ay0N1fcvcF\nkeejgWnAIZHFVnpj7v6Bu38TeT6Z4Eyt3L4X4FZ33+ju44HXCT7E84DTgZvdfbW7zwLuAc4pZzt3\nufsKd58NjCPyfigrTmAD8Bsza+juv7r7lxXEKFVIiUASsSDq+Vpif7CmwiqgYal5O7H1A3tljOUr\nY7x+p8i8RPYVz7bKem0jtk8sJVoQNIuV+Dl6oZmda2ZflpxNAR0JvsGXycwOiXRULzKz5QTNNTHX\nB35x97VR07OA3SKvqRWZjo6tZTnbin4/rCFIqLGcStA8NNPMisyseznrShVTIpDK2tKcEWkaWlnO\no295G4rhB6CmmbWNmteZoJ2ZyM/OUTHUA9qUWt6l1GunxNhXMtv6BtjbzOqXWv4NZZsPRF+RtOW5\nme0JDAeuBJq4e+PIfkq+YZfVuf4s8C+glbs3Ah6l/P/rxma2Y9T0nsA8YAmwEWhdKrY55Wwrlu3i\ndPfP3f2PQPNIvKMrsV1JEyUCqazo5oyfI+3BsR7PbXmRWR0zqxuZjH5ectnhT5FtrgbGALeb2Y5m\n1gM4AXg6svorQEczOyWyjQLgK3f/IbL8H8B1ZtbCzFoC1xF00pbsq8jMCiKT/6rstiLrfAUUmFld\nMzuF4Fv8yzGO22jgZjNrZGatgP5Ry+oRfIguAWqYWb/ItkosBFqZWa2oefUJvuVvMLODgTOp+Gqs\nwsilvIcRtPO/6O7Fkdj+amb1I0npWuCZCrZVwtj6ntgmzsi+zjKzndx9M8GZ1uY4tytVQIlAKsup\n+AOnLN8TNCO0AMYCq6Ou2d8d+DBq3SsIrrJZRPCBdFmkMxJ3X0LQ3PBXgitwDgS2dEq7+2ME/RiT\nCTp3X3P34VHbblWyL3dfnOS2zoi8ZllkG6e6+9IYv38hQfPLT8BbBEnGI/v5lqBd/hOCZpeOpY7H\newRnGgvMbFHUMbrdzFYAtxJ12WsMCwiu8plHkFQvjUp4/YHVwAxgAvBPYFTJYWDbv3fpv3308rLi\nPBv4ycx+JbhK66wK4pQqVHJpmEjoIneuDnD379O8n1bA8+7eI537yTRmlg887e67hx2LZJaaYQcg\nUsLjvHM1BfuZA+RUEhApj5qGRHKLmgBkO2oaEhHJcTojEBHJcVnZR2BmOo0REakEd9/uzu+sPSPw\nevXwn39O+eBLBQUFaX1deesluqz0vESnq/txi3d+dT1uFa2TyHFL9P2XzmOm41b5RyxZmwhYvZo3\n2l/DZZeldrP5+flpfV156yW6rPS8iqbTKROPW7zzq+txq2idRI5bZd5/6aTjlmLpzD5pzGq+ecd6\n7uB/6/G6S3wKCgrCDiEr6bglTsesctJ93IKP/O0/UzPqjMCCsn0TzGxYRdWratxeCEC/L/vD2rXl\nrSoRGfHNIwvpuCVOx6xywjpuGXX5qJkdTlCYYwHBUL0/xljPfcMGlrc9gEY/T4Zbb4Xbb6/SWEVE\nso2Z4WV0Fqc9EZjZEwQDWy1y905R83sCQ4E8YIS7DzYzc3c3s52Be9397BjbdHfn/ds/5A8Fh0Ht\n2jB5MrRvX9bqIlIFzMoqQyBhKeuzPVYiqIqmoVFAz1LB5AEPRebvB/Q1s319a+TLgToVbXhJhx78\nI68fbNjAhM5XMvCmzDm7EclFZbU/6xFKP2pC0p4I3H0CwWiH0Q4Gprv7THffSFBV6SQzO9nMHiUY\nkfHBirbdqxfkTxzC5kZNOGzdu+z9WUUDL4qISGlh3VDWkm2rNM0BDnH3uwjGma/QoEGDtjzPv7Af\n+ffcw2n/vRZ+PQ522imVsYqIZKWioiKKiooqXK9KOovNrDXBGO6dItOnAj3d/eLI9NkEiaB/zI1s\nuz3fJu7iYha2O5RdZvwXrr4ahg5N8W8gIhWJtD+HHYYQ+28RZh9BWeYSFCEpsTsJlsQbNGjQ1kxX\nowYfnjWMzdTAH3yQCQ98yYIF5b5cRCRnFBUVbdOKUlpYZwQ1CSpVHUlQKWki0Ncj1afi2J6Xjvut\nt2Dd5dfyx5lD+aLmIcwZ/TEnnZxRt0mIVGu5eEbw7LPPMn/+fCZOnMjJJ5/MGWecsd06bdq0Yc6c\nOTRq1Ii7776bc889F4B//etffPvtt9SoUYOWLVtyzjnnlDu/xMSJE3nvvfe4+eabY8aV6BlBVfRe\nP0fwYb+eoF+gX2T+cQTJYDpwc4Lb9DL9+qt7ixbu4F9e/ljZ64hIWsT8v6ympk2b5g888IC7uy9e\nvNgbNWrkM2bM2G694cOH+6xZs3zjxo1b5i1fvty7deu2Zbp79+6+ZMmSMucvXrx4y/TmzZv92GOP\n9cLCwnJji/W3IKw7i929r7u3cPc67r67u4+KzH/T3fdx97bufmei292maahEw4Zw330A7Pf0QFi8\nOPlfQESkDN988w1DhgwBoFmzZrRt25Yvvvhiu/Vq167NHnvsQc2aW6/NGT9+PPvtt9+W6c6dO/P+\n+++XOX/cuHFbpl988UWOOuqohM+8KmoaysphqIHYv1SfPnx51Qi6Ln6H9VffyPpHR9GwYZWGJiJZ\nasaMGTz++OMxl3fv3p2TTjoJgOOPP54333wTCFpW5s+fT9u2bbd7zWeffcb69etZsWIF7du358QT\nT9zSVFSiUaNGTJs2jSZNmpQ5H2Dx4sXk5eXRvHlzVq9endDvlZ+fT35+PoWFhWUuz9pEEJMZTx/y\nMPv9pyN1nnuSSydewJPTDws7KhFJ5Z3HleyLmDdvHqNGjaJLly6MHz+eyy+/nGbNmrFq1Sp23XVX\n9t57b+68M74Gilq1atGxY0cAXn/9dQ488EC6dOmy3XpHHnkkJ598MgBdunTh8MMPZ/ny5dStW3fL\nOrVr12bVqlWYWZnzAcaMGcMll1zCP/7xj0r97uXJ2t7UMpuGIu59rR11bhsIwF/mXgEbN1ZhZCKS\niVavXs3JJ5/MZZddRq9evejduzfXXnstb7/9Nk2aNKn0dpcvX86TTz7JM888U+bykjMIgMaNG1NU\nVESDBg22ad5Zu3YtTZo0iTn/008/5ZBDDql0h3zuNQ2VGDiQ9U88Q7s5U+D+++FPf6qSuEQkhpCv\nKHrhhRc44IADaNq0KQDNmzfnm2++wcyoXbs2kFjTEARNQnfddRcjRoygfv36zJo1iz333HPL8mee\neYZXX32V0aNHA0EyqlmzJm3atOHzzz/fst7SpUvp1q0bjRo12m5+165dmThxImvWrGHs2LF89NFH\nrF27lldffZUTTzwxrt+9oqah0MfEqMyDOK9O+PGhN9zBi+vV81+n/OybNsX1MhGphHj/L8Py6KOP\n+j333LNlesqUKd6hQ4ektnn//ff7559/7vPnz/dPP/3Ui4qK3N19+vTpXlxc7BMmTPD33nvP3d1X\nr17trVu39tWrV/uqVau8Y8eOW7az//77+8KFC2POj1ZQUOCDBg0qN65YfwtiXDWUUcNQx8vMvKCg\nYEuWi2XyZJjRrTcnbXqZlzmFTt+/rAFKRdIk0+8jWLFiBYMHD+bQQw9l48aN1KtXj5EjR3L00Udz\nxhlnsOOOOya0vQ8//JAjjjhiy+9sZvz888+0bNmSbt26MXLkSLp27co///lPFi9ezKxZszjjjDM4\n5JBDAHj66aeZNWsWxcXFtGnThrPOOqvc+QCjR4/mrrvuwswYOHAgffr0KTO20n+LkqEmCgsLy7yP\nIGsTQdxxz5kDHTrA6tXMHf46LS8+Pr3BieSoTE8EuSRbhpioOq1aQaRdrPkdqmYmIlJa9U8EAAMG\nMLV2J2rPnsFnp9zJ9OlhByQikjmyNhGUd/nodmrV4uOzhwGw/1uDmf7GD+kLTEQkw2TEoHOpllAf\nQbQLL4QnnmBhp6PYZdLbqb3BRSTHqY8gcyTaR5BbiWDJEla23IcGG5bBc89BGSMFikjlKBFkDnUW\nl6dZM57pOBiABWdeywvDfw05IBGR8OVWIgDO++ACNh3UnV19AZ1evC3scEREQpfVQ0xUdENZWXas\nXwOGD6O42wF0eO8h+PJ86No1LTGK5BpTv1tGqqh2cW71EUQZ+5trOfbboUxvegjL3/iYAw/OuZMj\nEckx6iMopfH9haxp3IK2Sz8l78kRYYcjIhKanD0jAGD0aDj9dDbUb0ztGd9D8+bJb1NEJEPpjKAs\nffrwVfOjqb3qF9494EbmzQs7IBGRqpfbicAMe+RhNteszVGzn6T4gwlhRyQiUuVyOxEAnXu3I+/P\nQTWz5oNUzUxEck/WJoKExhqqyMCBzKrZhjo/TOGNY+9PzTZFRDKExhqKU9HAt8gffByrrR71Zn0H\nu++e0u2LiIRNncUVyL+rJ6uP6009Xw3XXBN2OCIiVUaJIMrSW+5jFfVgzBhevvCNsMMREakSSgRR\ndjuoFesGBtXMfj/mKlUzE5GcoEQQpVYtaHb7ABbt0okmy3/i81PvZM2asKMSEUkvJYLSatVixg1B\nNbNObw5m9ZeqZiYi1ZsSQRm6X38oXHABddhAwz9fCVl4ZZWISLyUCGIZPJhl1oQ6499lyq0vhB2N\niEjaZG0iSOkNZWVp1owPjguqme1+77Xwq6qZiUh20g1lySguZkG7Huw64xNmnTSAVi/fT15e+ncr\nIpIOuqGsMmrU4IdrhrGZGrT690NsnPhl2BGJiKScEkEFDu/fmbxrBpBHMbWvvhyKi8MOSUQkpZQI\n4lFYyDxaUOOzT/n0ElUzE5HqRYkgHg0bUnzPfQB0fm4gLFoUckAiIqmjRBCnVtf24btWR1N3zS98\nfNhNYYcjIpIyumooAfM+mMYuR3Ukb9MGGD8eDjusymMQEaksXTWUAi2OaMfPfYNqZr+ceQXF61XN\nTESynxJBgtZdM5CFDdrQeM4UfKiqmYlI9su4RGBm9czsMzPrFXYsZdm32w7sMvohANbdPIgvX50d\nckQiIsnJuEQA3Ahk9uA+PXuyKlLNrNU9qmYmItkt7YnAzJ4ws4VmNrnU/J5mNtXMppnZTZF5RwPf\nAovTHVey6g+/jzV59Wk+fgwLn1Q1MxHJXlVxRjAK6Bk9w8zygIci8/cD+prZvsARQHfgTOBiM9uu\ndztjtGrFuz2CamZ1/6RqZiKSvdKeCNx9AvBLqdkHA9Pdfaa7bwSeB05y91vc/VrgWWB4KNeIJuDE\nd/ozr2kndlr6E1PP/VvY4YiIVErNkPbbEojuZZ0DHFIy4e5PVbSB6CFV8/Pzyc/PT1108apVi3m3\nDqPFNT1oM2YIKz47m4YH7VP1cYiIlKGoqCiu4fqr5IYyM2sNvObunSLTpwI93f3iyPTZwCHu3j/O\n7WXUycIPh11I+w+fYFKzI+m86B3I4BYtEcldmXZD2Vxg96jp3QnOCuKW9sI0CWj/ymA21G9C5yXv\nwQuZfcGTiOSejChMU8YZQU3ge+BIYB4wEejr7t/Fub2MOiMA+PLKEXR95GKW77ArjeZPhZ12Cjsk\nEZFthHZGYGbPAR8D7c1stpn1c/dNwFXAWILLRV+INwlkqjZ/vYBFbX9Lo7UL4Lbbwg5HRCRulToj\nMLMO7j41DfHEu38vKCgIr5M4hv+NmkTnC7phBjW++By6dg07JBGRLZ3GhYWFZZ4RVDYR/Ozue6Qk\nwkrIxKYhgCVLYObJ13Lgh0P5ofEhtF/yMdTIxJu3RSQXxWoaipkIzOzBcrZ3vrs3SFVwicrURACw\nduEK2HdfdvhlHjz2GFxySdghiYgAlesjOB+YAnwBfB71+ALYkIYYE5JJVw1F22GXhiwcGFQzW3mV\nqpmJSPgqfdWQmY0DbnH3j8pYNtPdW6coxoRl8hkBwLq1zvLf9mTXSW8z5aDz6ThxVNghiYhUqmmo\nCbDO3dekO7hEZXoiAFj632ns1KMjNTdvYMV/xtOwl6qZiUi4Em4acvdlmZgESmRq01CJpt3bMe+c\noJrZ8jMvh42qZiYi4UjLDWVmVujuBUnElZRsOCMAYO1aVu3VifoLf+THy+6mzbA/hR2RiOSwVN9Q\n9nmS8eSGHXZg4a1BNbNWIwfBbFUzE5HMUyVDTKRa1pwRREzr2od2X73EjC6nsPeXL4cdjojkqEwb\ndC5pmd5HEK3V6PtYV6s+e381hp8eVjUzEalaGTHoXKpl2xkBwE/972Wvh65n6U57UfuHb2iw8w5h\nhyQiOabanRFkm73u7c+iXTvR9Nef2FioamYikjkqTARmtquZjTSztyLT+5nZhekPrZqpVYudXxoG\nQINhQ1j+6fchByQiEojnjOBJ4G2gRWR6GnBtugKq1g49lOlHXEAt38CPPa+ELGveEpHqKZ5E0Mzd\nXwA2A0SKzW9Ka1RxyKbO4mhtXxrMxoZNOGC5qpmJSNVIurPYzIqAU4F33b2rmXUHBrv7ESmMMyHZ\n2Fkc7fsbRrDP3y9mVYNdqT9b1cxEpGokPNZQ1AsPAB4EfgN8AzQHerv7pHQEGo9sTwRrVxezsH0P\nWs/7hP8dNoAuRferbIGIpF2lE0HkxTWBDoAB37t7qMNQZ3siAFj49iSaHdsNgE2ffE6d7qpmJiLp\nVenLR81sBnCxu09x98nuvsHM/pOWKHPILsd0Ju+aAeRRzKpzLofi4rBDEpEcFU+DxEYg38xGmVmd\nyLyWaYwpdxQWsqxuC5pO/xRGjAg7GhHJUfEkgjXufjrwHTDezPZMc0xxydarhrbRsCELbgqqmW24\nXtXMRCQ9UnHV0Jfu3jXy/CjgYaCJuzdPYZwJqQ59BCWWLHZmd+xJ10VvM/P359P6fVUzE5H0SOaq\noRPc/bWo6T2B89z99tSHGZ/qlAgAVn05jVoHdKKOr4fx4+EwVTMTkdRLuLPYzPaNPJ1nZt1KHkBT\n4PU0xZmT6ndtx8orgmpmi/uompmIVK3yahY/7u4XR24o224ld/99mmOLqbqdEQBsWrmWX1p1ovmK\nH1l0wxB2HnJD2CGJSDWT1H0EmaY6JgKAFaPfouHpx7HWdmSHWVNh993DDklEqpFk7iPoY2YNI89v\nNbMxkSYiSbGGp/Vk3u96s4OvYWHfq8MOR0RyRDyXj97m7ivMrAdwJPAE8Gh6w8pdOzx6H6usPrt8\n9AqLRqkrRkTSL55EsDny8/+Ax939P0Ct9IWU2xp3asV3ZxQCsPbi/iydszbkiESkuosnEcw1s+HA\n6cDrZlY3ztelVbW4oSyGg57qz+q9O7Hn5p9YfqOqmYlIclJxQ1k9oCfwtbtPM7PdgE7u/nYqA01E\nde0s3sZHH0GPHqynNnPf+Jq9j9sn7IhEJMvpqqEstKDXhez6xhPMancke37/Dth2fz8RkbgpEWSj\nJUtYvfs+1Fu3jE1PP0fNs88IOyIRyWKVubO4bnpDkgo1a8Z35w0GYPE517J0xq8hByQi1VF5nb4f\nA5jZM1UUi5ThwEcuYE2X37IbC5h62m1hhyMi1VB5Q0x8A/wNuAP4E0F1shLu7mPSH17ZcqZpqMSk\nSRR37YY7LHvrM5ofq/v5RCRxlbmz+DLgMGAn4ASC+whKHiekI0iJoXNnph8XVDOb90dVMxOR1Irn\n8tGL3D2jymfl3BkBsPjHFdTusi87rZrHuDMe4/fPXRJ2SCKSZZKpR1AbuBw4PDKrCHjU3UMbKzkX\nEwHA8uGjaXTp6SyjMXPfnUqnI3cOOyQRySLJJIKRQE3gKYJ+gnOATe5+URqC7ABcTVDzYKy7j4yx\nXk4mAtzZeFRPar3/Ngt6ns+ub6qamYjEL5lE8LW771/RvFQysxrA8+5+WozluZkIAKZNY137TtRl\nPd8+Op79LlU1MxGJT6WHoQY2mVnbqA21ATYlsOMnzGyhmU0uNb+nmU01s2lmdlPU/BMIKqA9H+8+\nckq7dkw/Nahm1ujPqmYmIsmL54zgSGAU8FNkVmugn7u/H9cOzA4DVgH/cPdOkXl5wPfAUcBc4DOg\nr7t/F/W6f7v7STG2mbtnBADr1rFw547ssvJHGDIEblA1MxGpWFJDTETuMt6HoGTlD+6+LsGdtwZe\ni0oEvwUK3L1nZHpgZNVPgFOAusB37j40xvZyOxEAnxa+xSGDjmN9zR2p8f1Uau2tamYiUr5kmoZw\n93XuPsndv040CcTQEpgdNT0HaOnuH7j71e5+aawkIIG9r+jJpPa9qbNpDZuuUjUzEam8miHtN+mv\n89Fja+fn55Ofn5/sJrNK8+bQ/P2hrGz1Fg3efAVefx169Qo7LBHJIEVFRXHVbamS0UfLaBrqDgyK\nahq6GSh298Fxbi/nm4ZK3L3bvdyw4HpWNtuLBrOmwI47hh2SiGSopJqGzKylmR1qZoeb2RFmdnjF\nryrX50A7M2sduWHtdODVRDZQnSuUJeLSKf2ZbJ1osOQnFl93Z9jhiEgGSkWFssEEH9TfsrV+Me4e\n13hDZvYccATBTWKLgNvcfZSZHQcMBfKAke4e96eYzgi2tWrsR9TvGVQze/WOr+lzi6qZicj2krmh\n7AeC0pTr0xVcopQItvf9YReyz4dP8L4dyR82q5qZiGwvmaahH4HaqQ8pOWoa2tY+rwxm005N+IO/\nx7/6vqD7zERki1Q0DY0BOgPvASVnBe7uA1IUY8J0RlC2zY+NIO+yi5nPriyZMJVOPXYKOyQRySDJ\nNA2dX8Zsd/enUhRbwpQIYiguhh494JNPeG3vAZzw4/1hRyQiGaTaFa8vKCjIyfsHKjRpEpu6HIDh\n5H3xGXRTNTORXFdyP0FhYWFiicDMXnT3PqUHi4vwdI4+WhGdEZRv5inX0fqV+/h5t4PZY84nUCOu\nq4RFpJpL+IzAzFq4+7zIzWDbcfeZqQwwEUoE5du4dAWr99yXRqvnsfSvj9L0z5eGHZKIZICErxpy\n93mRnzOBdQQdxp2AdWEmgRK6aii2Wk0bMnPAfQDU+MvNsGhRyBGJSJhScdXQRcBtwLjIrHzg9ljV\nw6qCzggq9tabjh3fk2N5m5lHnEfroifDDklEQpbsDWW/dfelkemmwCfu3j4tkcZBiaBiP/0Ej/1p\nGoPGBNXMZj/zAbuflezIICKSzZK5oWwJQWGZEqsi8ySD7bUX3PlSO37sHZR6+PXsK5j4ke4yE5Ht\nxUwEZna9mV0PTAc+NbNBZjYI+C8wrYrii0l9BBUzg988PZCljdrQkW+o9fBQfv017KhEpKpVuo8g\n8qFfstBKP3f3wpRFmSA1DSVmyTNv0eyc41jNjvTrPpXRn6iamUguqnY3lGVj3GFaeEQfdhn/EmM4\nmY7fj6F9aD08IhKWpOoRSPZrNGooK6nPKbzCupdfDzscEckgSgQ5os7eLVn1p6A1r9WQ/qxYsCbk\niEQkU2RtIlBnceJ2u3MA39TcnybLf2LYHncyb17YEYlIVUjFDWU7AxcDrdla7N7d/YLUhJg49RFU\n3oh+H3HRkz3YQC0OqDWZ137Yh9atw45KRKpCMjeUfQKMB74AiiOz3d1fTnmUcVIiSM6yUy6kyStP\n8C5HMv+pd+j1f0aTJmFHJSLplkwi+Mrdu6QtskpQIkjSkiVs2Hsfaq9cxhk8x+v1z2DlyrCDEpF0\nS+aqof+YWa80xCRhadaMFX8eDMB9XEuNVbrLTCSXxXNGsArYEdgAlIxR4O7eMM2xlReTzgiSVVzM\nx3k9+B2fMLLeAC5cpWpmItVdpc8I3L2+u9dw97ru3iDyCC0JlNBVQ0mqUYN5twxjs+Vx/uqH+PCB\n/4UdkYikSTJDTOzr7t+ZWZm1Dt09tE8OnRGkzpLzrqPZP+7jUw7mgPWfULN21l5RLCIVqEyFssfd\n/WIzK2LrOENbuPvvUx5lnJQIUmjlSuY27EBL5rHo9kfZ+VZVMxOprjTWkMT0l31e5K8/nMbyGo1p\nNH8q7Lxz2CGJSBporCGJ6a9TezOWY2hU/AvceGPY4YhIFVMiEDCj/pMPs4468NRTLH91PGs0FJFI\nzlAiEAAOPKMtdxFUM5tz0hXc/TdVMxPJFRUmAjPrYWb1I8/PMbN7zWzP9IcmValOHWgzfCDTCaqZ\nrfrrUCauG5GVAAAaTElEQVROhA0bwo5MRNItnhvKJgP7Rx5PAiOA09z9iLRHFzsmdRanyXm7vMVT\ni4JqZh2YypBnd6dv37CjEpFUSKazeFPkU/ePwMPu/jDQINUBJko3lKVH83N68rL1ph5ruJ+rOfNM\nmD497KhEJBmpGIZ6PPAW0A84DFgMfOXunVIXZmJ0RpBmc+eyeo8O1CteRS/+w35/6sXdd4cdlIgk\nK5kzgtOBdcAF7r4AaAnoY6E6a9mSzbcG1cwepD8P/12XEIlUZ/GMNTQfGAPUicxaAvwrnUFJ+Bre\nMoB5zfZnb37iZu7ko4/CjkhE0iWeq4YuAV4EHovMagW8ks6gJAPUrMm9bR8B4CYGc0GP7/n555Bj\nEpG0iKdp6EqgB7ACwN1/ADQGQQ6Y2vRQRnAhtdnIw1zJsEeczZvDjkpEUi2eRLDe3deXTJhZTcoY\nhE6qn7w8GMhdLKUJR/EeMwc/zwknhB2ViKRaPIngAzP7C7CjmR1N0Ez0WnrDkkzQvz+cc00zbmQI\nAPdyHYunq5qZSHUTz+WjecCFwDGRWWOBEWFev6nLR6uOO2zeWMyMlj1ov+QTHqA/A/yBsMMSkUpI\npnh9PWCdu2+OTOcBddw9tGsKlQiqXmebxBccgOGsn/AZO/Yos16RiGSwZO4jeB/YIWp6R+DdVAVW\nmpmdZGbDzez5SFOUZICv6cwDDCCPYiYfdjkUF4cdkoikSDyJoI67ryqZcPeVBMkgLdz93+5+CXAZ\nwc1skiEKKGQuLTiEiSy56/GwwxGRFIknEaw2swNKJszsQGBtIjsxsyfMbGFkALvo+T3NbKqZTTOz\nm0q97BbgoUT2I+nTqhWsogHXMBSAHW6/mYNbLwo5KhFJhXgSwTXAaDP70Mw+BF4A+ie4n1FAz+gZ\nkb6GhyLz9wP6mtm+FhgMvOnuXyW4H0mTl1+GWrXgJYJqZvXW/8IVs27kggvCjkxEkhVXzWIzqw3s\nQ3D/wPfunnDVEjNrDbxWMlidmf0WKHD3npHpgZFVVwPnAZ8RDG73WBnbUmdxSHr0gAUfTWcKHanL\neg7nAw667nDuuSfsyESkIsnWLD6QoB7BAQTf3M9NQUwtgdlR03OAlu7+oLsf6O6Xl5UEJFwvvQRv\nfL+1mtkjXMED926ke3f4+uuQgxORSqlZ0Qpm9gywN/AVED3AwD+S3HdSX+mjx9bOz88nPz8/yXAk\nHrvuGjz+2WogZ895ho58wzUM5e+f3sDChWFHJyLRioqK4qrbEs99BN8B+yXbFlNG01B3YFBU09DN\nQLG7D45jW2oaCtnrr8OaV96iz8igmtm+fMfIt/fgaF3wK5KxkmkamgLslvqQ+BxoZ2atI30QpwOv\nxvtiVSgLV69e8MdhPXmRoJrZUK7hmGOgWbOwIxOR0lJRoawI6AJMBEoGn3N3PzHeIMzsOeAIoCmw\nCLjN3UeZ2XHAUCAPGOnud8a5PZ0RZICNG6F17blMpQMNCKqZvUEv+vWDJ54IOzoRKS2ZISbyI08d\nKNmAu/sHKY0wAUoEmWHjRqhdG67lXu7lemawFx2Zwlp2RH8ekcxT6aYhdy8CZgK1Is8nAl+mOL6E\nqWkofLVqBYPSPcAAJrG1mpmIZJZUNA1dAlwMNHH3NmbWHhjm7kemMtBE6Iwgs/zvfzD+zo+45qUe\nbKAWnZjMD+zD22+jzmORDJJMZ7EqlEm5unWDk/++bTUzcObPDzsyEYlH1lYoU9NQ5omuZnYGz9O0\nadgRiQikpmnobmA5cC5wFXAF8K27/yV1YSZGTUOZZ8UKuPdemF04kpFcxHx25Z+3TOWP5+1E27Zh\nRycikNxVQzWAi1CFMolDDSvmQ3rwO4JqZlfzgK4gEskQlUoEkWagKe7eIZ3BJcrMvKCgQENLZCAz\n6MTX/I9uGM5BfMYrM7ux555hRyaSu0qGmigsLKz0GcG/gQHuPitdQSZKZwSZyyJvsXu4juu4j085\nmN/xMUf8Po9Ro6BJE2jQINwYRXJVMk1DE4CuBPcPrI7MTujO4lRTIshcJYmgPiuZSgdaMo9LeZTh\nXArAQw/BlVeGGKBIDksmERzB1juKS+jOYimTRb1TevMiL3Iay2hMB6ayOHLV8fjxcNhhIQUoksMq\ndR9BpI9guLsXlXqElgRK6PLRzHTHHfDVV/D441urmTXhF4Zw45Z15s4NMUCRHJSKy0fVRyAJGzEC\nLr4Y2rBtNbMJHM4TT0C/fmFHKJJ7krmzuAnwjZm9b2avRR5xDxctualWreDnj2xbzawmG5k6NcTA\nRGQ7iYw+uo3IAHSh0BlB5pszJ+gLOOssqMM6ptCRtvzIDQyhQeENDBgAO+20bZ+CiKRXpTuLM5ES\nQfYo+aA/hrGMpeeWamaz2YOffoLWrUMNTySnVLppyMxWmdnKyGO9mRWb2Yr0hBk/dRZnh1tuCX6+\nzbGMps+WamYAa9aEGJhIDkm6s3iblYPhJk4Eurv7wKSjqySdEWSPkuI1AENvmMsFd29bzez22+Ga\na3STmUhVSKazeAt3L3b3fwE9UxaZVGslncYAddu05DZuB+AhrmIH1nDbbfDFF7BhQ0gBikhcTUOn\nRj36mNldwNoqiE2qGXd4kP5MYn/2YiZ/5m8A/P73cN11IQcnksPiOSM4Afi/yOMYYCVwUjqDkuql\nfv3gZ716sJmaXMEjANzIENrzPQDz58PNN4cVoUhu01VDknZvvw3HHgvjxsGrr8J998HjXMRFjORd\njuRo3qFkFBP9WUXSJ5mrhp4ys0ZR043N7IlUB5goXTWUPY45Zuvz226D2bO3r2ZW4rzzYNIkePrp\noNiNiCQvFUNMfOXuXSqaV5V0RpB9Lr00aPopuW/ADC5gazWzDkxlBTtt85r//Q+6dq36WEWqq2Su\nGjIzaxI10QTIS2VwUv099ti2N4+99x6Moh8f81t2YwF3cOt2r9m8ueriE8ll8SSCe4BPzOwOM/t/\nwCfA3ekNS6q7P/wBnBpcxqNsIo8reZiu/G+bdX79NaTgRHJMXJ3FZvYb4A+AA++7+7fpDqyCeNQ0\nVA3EqmZWHDnhbNkyGLNIRFIjmcI03YFv3X1FZLohsK+7f5qWSOOgRFA9VFTNDHQVkUgqJdNH8CjB\nvQMlVkfmiSTl0EODn6towDUMBeAuBtKcRVvWmTQpjMhEcktcQ0xEf/12982os1hSIPrq35fozVsc\nS2OWb1PNrEto16aJ5I54EsFPZjbAzGqZWW0zuxqYke7ApPqrWRP22KNkyriKh1hHHc7nKQ5j/Jb1\nzOD558vchIikQDyJ4DLgUGAuMAfoDlySzqDioRvKqoeNG7c+L6uaWYlvvoFWrVTIRqQyUjoMdaZQ\nZ3H10bw5LFmydbp0NbO/c8N2r9GfXqRykhliYncze8XMFkceL5tZq/SEKbmmV69gCIpnnw2m11OX\nK3kYgEEMYnd+DjE6kdwQz+Wj7wL/BJ6JzDoLOMvdj05zbOXFpDOCambVqm2L07zAaZzGi4zhZE5l\nzDbr6k8vUjnJ3Ecwyd07VzSvKikRVE/dusGXXwbPWzCXqWxbzaxEcbH6CkQqI5n7CJaa2Tlmlmdm\nNc3sbGBJha8SSVD0h/s8tq9mVqJGDfjPf6o6OpHqK55EcAFwGrAAmA/0AfqlMyjJTfvtt+10WdXM\nSpxwAnz8cRUGJ1KN6aohyRgLFsAbb8CFF26d91s+5mMOZQO16MRkfmCfLcvOOQeuvBIOOQSmTQvm\ntWtXxUGLZJFK9xFkIiWC6mvNmqCkZbSSambv8QeO4l1KqpmVcIfu3YPn//1v1cQpko2UCCRrdOwY\n3EBWoilL+J59aMoy+vIsz9N3m/Xdt/Yv6G0hElsyncVVxsz2MrMRZvZi2LFI5lhKM25kCAD3ch0N\nUaECkVSK54ayW6Ke101nMO7+k7tflM59SHYqr5rZ3LkhBSVSTcRMBGY20Mx+R3CVUImEr9MwsyfM\nbKGZTS41v6eZTTWzaWZ2U6Lblerv2mu3Pi+vmtlXX4UQnEg1Ut4ZwVSCJLCXmX1oZo8DzcysQ4L7\nGAX0jJ5hZnnAQ5H5+wF9zWzfBLcr1VSXLtC27fbzJ7M/DzCAPIoZxuXUIChq/H//t+16ZrDnnlUQ\nqEg1UV4iWA7cDPwI5AMPEJSqvMnMPol3B+4+Afil1OyDgenuPtPdNwLPAyeZWRMzexToorOE3PX4\n41vvMC6tgELm0oJDmMhFjIi5jZ81RJFI3MpLBMcCrwNtCArYHwyscfd+7v7bJPfbEpgdNT0HaOnu\ny9z9Mndv5+6Dk9yHZKkddoD69cteVl41M4BatdIdnUj1UzPWAne/GYJxhYCngQMImoY+Apa5+wlJ\n7Dfpi/yix9bOz88nPz8/2U1KhnrpJejdO2o6Us2sJ2MZwo3048ktyzZtqvr4RDJVUVFRXHVb4hl0\nboi73xh5/qW7dzWz5u6+ON5gzKw18Jq7d4pMdwcGuXvPyPTNQHG8ZwG6jyA3DBwIgwdve59AiTZM\nZwodqct6DucDJnD4dq+/8UbYay84/vjoSmgiuavS9xGUJIGI8yPz4k4CMXwOtDOz1mZWGzgdeDWR\nDahCWfV37rlBf0FZfqQtd3IzAMO4fJtqZiWGDIHLLw86ju+4A846S6OWSm4KvUKZmT0HHAE0BRYB\nt7n7KDM7DhgK5AEj3f3OBLapM4IcU9YHeHQ1sxsZzN3cuP1KEYWFMHUqPPec7j6W3KUhJiSrxfom\nfwxjGUtPVrMj+/Ids6m4DUhvHclVWTHERCLUNCQAb3Mso+lDPdYwlGvCDkckI4XeNJQOOiPIPeW1\n7ZdXzawseutIrqp2ZwSSW/LyYi8rr5pZLL/+GiSXyZMrXFWk2svaRKCmodyybFnws0WLspeXV82s\nvO2pypnkAjUNSbVhBrvtBvPnl728vGpm0dxh1ixo3RoeeSS4xFQkF6hpSKqFHXeMvewTfscILqQ2\nG3mEK4h1A/v558OxxwbPR46EKVOC5yeeCMOHb7++O6xalVTYIhkta88ICgoKNLREjnnuuaCu8XXX\nxV6nCUv5nn1oxtIyq5mV5aqr4MEHY1c5Ky4O+iiy8F9FBNg61ERhYaHuI5Ds9847cMwx5a9zASMZ\nyUXMZ1c6MJUV7FTu+hddFNzBrEQg1Z2ahiRnjKIfH/G7MquZVYbqIUt1p0Qg1Y5Tg8sZVmY1s7KM\nGBH0G1RkzBiYNCk1MYpkkqxNBLp8VD7/PPayWNXMYnnqqa3P69YNzgLOOWfbdXr3hr9VfGWqSMbR\n5aNSrUT3EXzxBRxwQOx167OSqXSgJfO4lEcZzqVx7aNmzaCuQd++8OyzQZNQjchXpmOOgbFjk/wl\nREKiPgLJORVVM4uluDj4qe8akiuUCKRaK6lm1pjlDClnmOpoSgSSa5QIJKvUiHrHxldkxriKh1hH\nHc7nKQ5jfNz7euEFOPpo+Oc/t85bvjy4uWz16qCCWokFC1QmU7JX1iYCdRbnps6dYdiw4Nv6XnvF\n95p4qpnF8u6723YaT5wI++8PS5cGZTRL7LYbfPpp3JsVqVLqLJZqbcCA4K7gitRhHZPpRDumV1jN\nrCKtW8MHHwQlMEvehmYwYQL06FHpzYqknTqLpVqKtwbxeupyJQ8DUEAhu/NzymPRdxPJVkoEktUS\nKUb/DsdsqWZ2P1enPBYlAslWSgSSU67lPlZSn5P5F734T6W2EesDX4lAspUSgWS1RM4IYNtqZg/S\nP65qZqXNmrW1oM3UqVvnf/118PO772Dt2uD5xInw9NOxtzV3LixcCCtXwrRpCYcikhJZ21msYagF\n4LXX4L77YNy4+F+Txya+4AA68zX/j79wK/8vZfG4B8npnXfgqKMqHrDu0ENh553h+OPhkkt0ViHp\noWGopdpbsgSaN0/sNfFWM0tUSSJ4+WU45ZSKE4EZNG4MV18NgwYpEUh66aohqbYq8+EZbzUzkVyg\nRCA56yYGs4SmHMn7nMHzYYcjEholAsl6lW1OWUZTbiK4PfherqMhv6YwKpHsoUQgWS+ZdvVUVzMT\nyUZKBJLTEq1mJlIdKRFI1kv2SpvJ7M/9XE0exTzKZRVWMxOpbpQIRIBBDGIuLTiYz7iYxyu9nccj\nL120CN5/f+v8M8+Et96Ciy6Cu++GX37Zum5ps2fDiy/CV1/Bmjjvd1u8GKZPr3i9r78OhtEW2Ya7\nZ90D8IKCAh83bpyLrFvnHpwXuPftG/y86qqt8+J99Ga0O/gyGnlzFib8+so+Gjd2HzQoeO7u3rXr\ntstHj674GAwY4N68ecXrgfuLLyZ3vCX7jBs3zgsKCjz4yN/+M1U3lEm1lejwE+C8yXH0ZCxPch79\neDINUW2v9A1lpeMePBhurGDU7KOPDmonVPRvYQaPPRbcxSy5RzeUiVSo8tXM0knfeSTdlAhEoiRT\nzUwkWykRiJQymJuYRlt+w7dcy31hh6MzAkk7JQKRUqqimplIJlEiEClDuquZJUJnBJJuSgQiMaSi\nmplINlAiEIkhFdXMRLKBEoFIOR6kP5PYn72YyZ/5WygxqGlI0i2jEoGZ1TOzp8xsuJmdGXY8Ipup\nyeUMA+BGhtCe79OyH33YS5gyKhEApwCj3f0S4MSwgxGB8KuZKUlIuqU9EZjZE2a20Mwml5rf08ym\nmtk0M7spMrslMDvyXENASsZQNTOpzqrijGAU0DN6hpnlAQ9F5u8H9DWzfYE5wO5VGJtIXMKsZqYz\nAkm3tH/YuvsE4JdSsw8Gprv7THffCDwPnASMAU41s0eAV9Mdm0giVM1Mqq2yhiRN9QNoDUyOmu4N\nPB41fTbwYALbS+kQrVI99eyZ+iGjOzHJN5Lnm6jhXfmiyoaqrurHQQeVPX/nncuef+mlW58PGRL8\nPP5499tuc+/WLZi+777g71Ky3jHHuL//vvupp26dd/DB7sOHu7dosXXeZ58Frxs7Npj+97/dL7zQ\n/eKL3a+/Ppg3Y8b2f/9x49z//Ocqe7tto3dv97lzg+dr1rgfeWTVx/D73wdDtEcjxjDUSX3Ax/so\nIxGcmmwiKCgo2PJQXQIpyy+/uJ93nm8Zg7/k+8NFFyX3Ifl3rnMH/5SDvAabQv/QzqaH+7bTp59e\n8WtKPsxLEvuBB26/zgsvbP/3L6lNEQZwf++94PnCheHEESTNcdt8VsZKBDWr6syjlLls7Qsg8nxO\nIhsYNGhQKuORaqhRo+AB0K5d6rY7iEGczgtbqpk9xmWp27hICvXokc+JJ+ZvmS4sLCxzvbA6ZD8H\n2plZazOrDZxOgn0CgwYNoqioKB2xSTUXnFRW3ioacA1DAbiTm2nOohREJZI+RUVF5X55rorLR58D\nPgbam9lsM+vn7puAq4CxwLfAC+7+XSLbHTRoEPn5+SmPV6qXxKuUxedlTuUtjqUxy7mbG9KzE5EU\nyc/PLzcRpL1pyN37xpj/JvBmuvcvuS3Zb/+xBdXMptCR8/gHT3AB4zkiXTsTSausvVZfTUNSWalK\nDtHVzB7hClUzk4wVetNQuqhpSDJBplUzEylLRU1DWZsIRDKBqplJdZC1iUBNQ1JZqe43yKRqZiJl\nUdOQSBVQNTPJZGoaEqkCqmYm2SxrE4GahhKn4xVIvGmoKK61MqGaWeYoCjuArJSu/1E1DckWSgSV\nVRTXWqWrme3D1DTGlOmKwg4gK6Xrf1RNQwmq7B8i3teVt16iy0rPq2g6nTLxuMWaP39+6fkVTcfv\nE37H41xEbTbyMFdChdXMKruveF5X0Tqxlpc1P555Fe0vlSq3r3jebxWtk8j7rTL/t2FQIiglEz/Q\nYi1TIih/Waz5CxaUnl/RdGIGcteWamZ9ea6CtSu7r3heV9E6sZaXNT+eeRXtL5Uqt6+qTATu2ZMI\nzNN3D37amFn2BS0ikgHcfbsRuLIyEYiISOqoaUhEJMcpEYiI5DglAhGRHKdEICKS46pFIjCzemb2\nlJkNN7Mzw44nW5jZXmY2wsxeDDuWbGFmJ0XeZ8+b2dFhx5MtzKyDmQ0zs9FmdmHY8WSTyOfbZ2bW\nK237qA5XDZnZOcAyd3/dzJ539zPCjimbmNmL7t4n7DiyiZk1Av7u7heFHUs2MbMawPPuflrYsWQL\nMysEVgLfufvr6dhHxp4RmNkTZrbQzCaXmt/TzKaa2TQzuykyuyUwO/J8c5UGmmESPG5CpY/ZLcBD\nVRdl5kn0uJnZCcDrwPNVHWsmSeS4Rc46vwUWpzOmjE0EwCigZ/QMM8sj+OfrCewH9DWzfYE5wO6R\n1TL5d6oKiRw3CcR9zCwwGHjT3b+q+lAzSkLvNXd/zd2PA86r6kAzTCLH7QigO3AmcLGZbXczWCqk\nvXh9Zbn7BDNrXWr2wcB0d58JYGbPAycBDwAPRdrQXq3CMDNOIsfNzBYCfwO6mNlN7j64KmPNFAm+\n144CjgQamllbd3+sCkPNKAm+13YGTgHqAuOqMMyMk8hxc/dbItPnAYs9TW35GZsIYohuAoLgTOAQ\nd18DXBBOSFkh1nFbBlwWTkgZL9Yx6w88GE5IWSHWcfsA+CCckLJCmcetZMLdn0rnzrOtGSX7e7bD\noeOWOB2zytFxq5xQj1u2JYK5bO0LIPJ8TkixZBMdt8TpmFWOjlvlhHrcsi0RfA60M7PWZlYbOJ0c\n7xOIk45b4nTMKkfHrXJCPW4ZmwjM7DngY6C9mc02s37uvgm4ChhLcEnVC+7+XZhxZhodt8TpmFWO\njlvlZOJxqxY3lImISOVl7BmBiIhUDSUCEZEcp0QgIpLjlAhERHKcEoGISI5TIhARyXFKBCIiOU6J\nQEQkxykRiIjkOCUCkRQws4PMbJKZ1YnUmJ1iZvuFHZdIPDTEhEiKmNkdBIVXdgBm52qhH8k+SgQi\nKWJmtQhGkVwL/DZd1aREUk1NQyKp0wyoB9QnOCsQyQo6IxBJETN7FXgW2BvYLVLWUiTjZVvNYpGM\nZGbnAuvd/XkzqwF8bGb57l4UcmgiFdIZgYhIjlMfgYhIjlMiEBHJcUoEIiI5TolARCTHKRGIiOQ4\nJQIRkRynRCAikuOUCEREctz/B8nKk4rQmxaGAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }