Method of Moments - distribution parameter estimation from moments of various powers

HW2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [Uniform](#Uniform-Distribution)\n",
    "## [Exponential](#Exponential-Distribution)\n",
    "## [Graphs](#graphs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniform Distribution\n",
    "* $\\large{ X \\sim U_{\\left[0,~\\theta\\right]} }$\n",
    "* $\\large{ E(X^k) = \\frac{\\theta^k}{k+1} }$\n",
    "* $\\large{ \\hat{\\theta} = \\sqrt[k]{\\overline{X^k} \\cdot (k+1)} }$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hats (means):\n",
      "\t1-th: 9.99778713205253\n",
      "\t2-th: 9.997965651217353\n",
      "\t3-th: 9.998193060547539\n",
      "\t4-th: 9.998396828610742\n",
      "\t5-th: 9.99856103363193\n",
      "\t6-th: 9.998690750021565\n",
      "\t7-th: 9.998794718628666\n",
      "\t8-th: 9.99888038000482\n",
      "\t9-th: 9.998953085178524\n",
      "SKOs:\n",
      "\t1-th: 0.0036147499599368926\n",
      "\t2-th: 0.0021838578795396494\n",
      "\t3-th: 0.0015600307160600502\n",
      "\t4-th: 0.001212244185115949\n",
      "\t5-th: 0.0009897439189595588\n",
      "\t6-th: 0.0008345643457651421\n",
      "\t7-th: 0.0007199455323093284\n",
      "\t8-th: 0.0006318172447374726\n",
      "\t9-th: 0.0005620230699197349\n",
      "Logs of SKOs:\n",
      "\t1-th: -5.622732593139269\n",
      "\t2-th: -6.126662296687705\n",
      "\t3-th: -6.463049768132449\n",
      "\t4-th: -6.715281938760002\n",
      "\t5-th: -6.9180643160135205\n",
      "\t6-th: -7.088600710866559\n",
      "\t7-th: -7.236334998386123\n",
      "\t8-th: -7.36691037535423\n",
      "\t9-th: -7.483967659233855\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = 10000\n",
    "theta = 10  # uniform dist theta parameter\n",
    "moments = 9  # number of moments to analyze\n",
    "\n",
    "hats = [[] for _ in range(moments)]  # 2d array. hats[k-1] = hats of k-th moment\n",
    "\n",
    "for q in range(1000):\n",
    "    data = np.random.uniform(0, theta, size=(n, 1))\n",
    "\n",
    "    mean = np.mean(data)\n",
    "\n",
    "    for k in range(1, moments + 1):\n",
    "        theta_hat = np.power((k + 1) * np.mean(data**k), 1 / k)\n",
    "\n",
    "        hats[k - 1].append(theta_hat)\n",
    "\n",
    "print('Hats (means):')\n",
    "print('\\n'.join(map(lambda k_hat: '\\t{}-th: {}'.format(k_hat[0] + 1, np.mean(k_hat[1])), enumerate(hats))))\n",
    "\n",
    "devs = np.array([np.mean([(theta - theta_hat)**2 for theta_hat in hat]) for hat in hats])\n",
    "print('SKOs:')\n",
    "print('\\n'.join(map(lambda k_sko: '\\t{}-th: {}'.format(k_sko[0] + 1, k_sko[1]), enumerate(devs))))\n",
    "\n",
    "print('Logs of SKOs:')\n",
    "print('\\n'.join(map(lambda k_sko: '\\t{}-th: {}'.format(k_sko[0] + 1, np.log(k_sko[1])), enumerate(devs))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Среднеквадратичные отклонения оценок **уменьшаются** с ростом порядка использованного для построения оценки момента\n",
    "\n",
    "Значит, для **равномерного распределения** _\"лучшей\"_ оценкой будет построенная по моменту как можно большего порядка"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Distribution\n",
    "* $\\large{ X \\sim Exp(\\lambda = \\theta) }$\n",
    "* $\\large{ E(X^k) = \\frac{k!}{\\theta^k} }$\n",
    "* $\\large{ \\hat{\\theta} = \\sqrt[k]{\\frac{k!}{\\overline{X^k}}} }$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hats (means):\n",
      "\t1-th: 9.998919732876985\n",
      "\t2-th: 9.997258219963548\n",
      "\t3-th: 9.997539613346781\n",
      "\t4-th: 10.001976836759138\n",
      "\t5-th: 10.016431671249041\n",
      "\t6-th: 10.050871767416885\n",
      "\t7-th: 10.115556158207653\n",
      "\t8-th: 10.21636266672021\n",
      "\t9-th: 10.353972823886743\n",
      "SKOs:\n",
      "\t1-th: 0.010187041653512217\n",
      "\t2-th: 0.012778194760757698\n",
      "\t3-th: 0.02113168699656047\n",
      "\t4-th: 0.04122588014889843\n",
      "\t5-th: 0.08505426358031817\n",
      "\t6-th: 0.16430323920577938\n",
      "\t7-th: 0.2832065901864233\n",
      "\t8-th: 0.44554380299922913\n",
      "\t9-th: 0.6634041490057427\n",
      "Logs of SKOs:\n",
      "\t1-th: -4.58663879249357\n",
      "\t2-th: -4.360015095036964\n",
      "\t3-th: -3.8569816116125697\n",
      "\t4-th: -3.1886890608820653\n",
      "\t5-th: -2.464465831116698\n",
      "\t6-th: -1.8060415389428768\n",
      "\t7-th: -1.261578647001039\n",
      "\t8-th: -0.8084597136740922\n",
      "\t9-th: -0.4103708983347736\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "\n",
    "n = 10000\n",
    "theta = 10  # exp dist theta parameter\n",
    "moments = 9  # number of moments to analyze\n",
    "\n",
    "hats = [[] for _ in range(moments)]  # 2d array. hats[k-1] = hats of k-th moment\n",
    "\n",
    "for q in range(1000):\n",
    "    data = np.random.exponential(1. / theta, size=(n, 1))\n",
    "\n",
    "    mean = np.mean(data)\n",
    "\n",
    "    for k in range(1, moments + 1):\n",
    "        theta_hat = np.power(math.factorial(k) / np.mean(data**k), 1 / k)\n",
    "\n",
    "        hats[k - 1].append(theta_hat)\n",
    "\n",
    "print('Hats (means):')\n",
    "print('\\n'.join(map(lambda k_hat: '\\t{}-th: {}'.format(k_hat[0] + 1, np.mean(k_hat[1])), enumerate(hats))))\n",
    "\n",
    "devs = np.array([np.mean([(theta - theta_hat)**2 for theta_hat in hat]) for hat in hats])\n",
    "print('SKOs:')\n",
    "print('\\n'.join(map(lambda k_sko: '\\t{}-th: {}'.format(k_sko[0] + 1, k_sko[1]), enumerate(devs))))\n",
    "\n",
    "print('Logs of SKOs:')\n",
    "print('\\n'.join(map(lambda k_sko: '\\t{}-th: {}'.format(k_sko[0] + 1, np.log(k_sko[1])), enumerate(devs))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Среднеквадратичные отклонения оценок **увеличиваются** с ростом порядка использованного для построения оценки момента\n",
    "\n",
    "Значит, для **экспоненциального распределения** _\"лучшей\"_ оценкой будет построенная по моменту 1-го порядка (**выборочное среднее**)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "<a id='graphs'></a>\n",
    "Наглядное сравнение зависимости СКО от порядка момента для равномерного и экспоненциального распределений:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAGHCAYAAAAXwu53AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcTfUbwPHPM2MsY5c1ZKssLdayZa+IokXKvmUtpUTI\nLlv9RqGIJDS2KKXixy+hEi2jtCElSVlSWbLFzPP749yZ7ow7Y+a6c8+dmef9ep3XzP2e7/me5+DV\nffqe7yKqijHGGGOMCZ4wtwMwxhhjjMlqLAEzxhhjjAkyS8CMMcYYY4LMEjBjjDHGmCCzBMwYY4wx\nJsgsATPGGGOMCTJLwIwxxhhjgswSMGOMMcaYILMEzBhjjDEmyCwBM8YYY4wJMkvAjDEBIyLXicgK\nEdkrIqdFZL+IrBORh7zq7BWRVT6u7Swi50VktYhk9yqPFJGRIrJdRE6KyFER+UBEOgfruYwxJtDE\n9oI0xgSCiNQD3gd+BhYAB4HSQB2ggqpe7an3E/C1qrb2uraj55p1wJ2q+o+nvKinzYrAEuADICdw\nD9AIWAp0VPsPmTEmg8nmdgDGmEzjSeAoUEtVT3ifEJHCyV0kIvcD84H38Eq+PBbiJF93quq7XuXP\ni8jTwOPAF8AzAXkCY4wJEnsFaYwJlPLAt0mTLwBVPeLrAhFpB7yK08vVxjv5EpHawK3AK0mSr3jD\ngN3AEyKSIwDxG2NM0FgCZowJlJ+BmiJyTWoqi8jdQDSwEWitqmeTVLkDUJwE7QKqGgssBgoC9f2M\n2RhjXGEJmDEmUP4DRAJfishmEZksIreIiK+hDjVwxm99CNzhI/kCqOL5uT2Fe24HBKh8CXEbY0zQ\nWQJmjAkIVX0PqAu8BVwPDAbWAr+KyB1JqhcEwoH9qnommSbzen5e8ErTS/y5fH4FbYwxLrEEzBgT\nMKoao6ptcRKsG4GJQB5guYhU8qq6HpgFdBaRZ5NpLj65ypvMee9zKSVpxhgTciwBM8YEnKqe9yRj\nI4D+QHbg3iR1HsJ5DfmIiIzy0cwOz8/rU7hVVZxxYt9detTGGBM8loAZY9Lb556fJXyc6wKsAcZ4\nL9bq8Q7O+K4uvhoVkTCgA/AXsDkwoRpjTHBYAmaMCQgRaZzMqVaenzuTnlDV80BbnARqmmdB1vhz\nW3DWBusuIq2SXovzevNKYEoyg/iNMSZk2Ur4xpiAEJGvcWZBrsRJtrLjLA/RDmeJihqqejyZlfDz\nAZuAa4B7VPVtT3lxnCQsfiX8D4EcwN38uxJ+J1sJ3xiT0VgCZowJCBG5FWecVz2gFE4Ctg9YDUyI\nX4xVRPbgJGBtklxfFCfBKgXcpqofeMojgUGetssD54GvgDmqGh2ERzPGmICzBMwYY4wxJshCZgyY\niDwoIj+JyGkR2SoiN6Si/ncickpEdohI52DFaowxxpjE/PgebywiMSJyRkS+F5GuPurc6/mOPy0i\n20XktiTnh4nIpyJyXEQOichKEbk6SZ1XRCQuybE6ME/tv5BIwETkPiAKGA1Ux1ndem1yG/iKSD9g\nAjAKZ7XsMcALyQzUNcYYY0w68uN7vCzOTOf1OMvJTAPmisgtXnXq4Ww39hJQDWeR5zdFpIpXUw2A\nGUBt4GYgAlgnIrmS3HINUAwo7jna+/+0gRESryBFZCvwiao+4vkswC/AdFV92kf9zcBHqvqEV9l/\ngBtVtWGQwjbGGGMMfn2PT8EZ63m9V9kSIL+qtvR8XgpEJpmwswX4QlX7JxNHYeAw0FBVP/KUveJp\n9+7APG1guN4DJiIRQE2cLBgAz4ym+G1NfMkBJN2+5Axwo4iEp0ecxhhjjLmQn9/jdTznva1NUr9u\nKuokVQBnceY/k5Q39ryi3CkiM0WkUAptBIXrCRhQGGdPuENJyg/hdBP6shZ4QERqAIhILaAnTtej\nz+5OY4wxxqQLf77HiydTP5+I5LhIHZ9tenrdnsN5Q+a9O8YanAWdmwJDcJawWe2p75psbt78EozH\neZe7xbMa9kFgPs4fbJyvC0TkSmAgzvT1v4MTpjHGGJMp5MHZFuw5Vf3B7WCSMRNnXHh970JVfc3r\n47eeNQt/BBoDG4IWXRKhkIAdAWJxEipvxXASqwuo6hmcHrA+nnoHgD7ACVX9PZn7DAQeDEjExhhj\nTNaVdNuwNH+Pe8p91T/utbNFcnUuaFNEngdaAg1U9UBKwavqTyJyBGcnjaybgKnqORGJAZoBqyCh\nG7EZMP0i18YCv3muuR94O4XqXwH069eP+vXrp1At44uKimLQoEFuh5Hu7DkzF3vOzCWrPCdkjWfd\nvHkzs2bNAs93qTc/v8e3ALclKbvVU+5dJ2kbtySpE598tQEaqeq+iz2LiJQCLsPpvHGPqrp+4GxV\ncgrnHW0lYDbwB1DEc34SsMCr/lVAR5zs9Uac7Uh+B65I4R4dAI2OjtbM7o477nA7hKCw58xc7Dkz\nl6zynKpZ41mjo6MVZ3B7Bw3M93hZ4AQwBWersf7AP8DNXnXqAmeBxzx1xuBMuKviVWcm8BfOchTF\nvI6cnvO5gadxlqkog5PQfQ7sACJ8PUuwDtd7wMB5P+uZOjoO5w/uS6C5/vs6sThQ2uuScJytSa4G\nzuF0IdbTVGS+xhhjjAmstH6Pq+pez9qdzwIPA/uBnqr6nledLSLSAWfdzwnAbqCNJh5g3xcnMdyY\nJKTuwEKcV6PX4ySGBXDemq0FRqnquQA8ut9CIgEDUNWZOJmsr3Pdk3zeCdQIRlzGGGOMubi0fI97\nyj7AWb4ipTZfB15P4XyKqzmoM2a8RUp13BIKy1AYY4wxxmQploBlQu3bu77DQlDYc2Yu9pyZS1Z5\nTshaz2oCJyS2IgoGz3vkRdHR0XTs2NHtcIwxxqd9+/Zx5MgRt8MwWUzhwoW54oorkj2/aNEiOnXq\nBNBRVRcHLbBMLGTGgBljTFa3b98+KleuzKlTp9wOxWQxkZGR7NixI8UkzASWJWDGGBMijhw5wqlT\np4iOjqZy5cpuh2OyiB07dtCpUyeOHDliCVgQWQJmjDEhpnLlytSoYRO9jcnMbBC+McYYY0yQWQJm\njDHGGBNkloAZY4wxxgSZJWDGGGOMMUFmCZgxxhgTwho3bkyTJk38ujYsLIxx48YFOCITCJaAGWOM\nCYoFCxYQFhbm8wgPD+fTTz91O0TX7Nixg7Fjx7Jv374LzokIYWH2dZ3Z2DIUxhhjgkZEGD9+PGXL\nlr3g3JVXXhn8gELEd999x9ixY2nSpMkFa3H973//cykqk54sATPGGBNULVq0sHXOklBVRMTnuWzZ\n7Ks6M7I+TWOMMSFjzJgxhIeHs2HDhkTlvXv3JkeOHHz99dcAbNq0ibCwMF577TWGDx9OiRIlyJMn\nD23atGH//v0XtLt8+XJq1apFZGQkRYoUoXPnzvz222+J6nTr1o28efPy22+/ceedd5I3b16KFi3K\n4MGDSbpvsqry3HPPce2115IrVy6KFy9O3759OXr0aKJ6ZcuWpXXr1mzevJnatWuTK1cuKlSowKuv\nvppQZ8GCBbRr1w5wxnvFv5L94IMPEsqaNm2aUP/cuXOMGjWKWrVqUaBAAfLkyUPDhg3ZuHFjGv+0\njZssATPGGBNUx44d448//kh0/PnnnwCMGDGCatWq0bNnT06ePAnA2rVrmTt3LmPGjOG6665L1NaE\nCRNYs2YNQ4cO5ZFHHuF///sft9xyC2fPnk2oM3/+fO677z4iIiKYPHkyvXv35o033qBBgwYcP348\noZ6IEBcXR/PmzSlSpAhRUVE0btyYqVOnMmfOnET37d27N0888QQNGjRg+vTp9OjRg0WLFtGiRQti\nY2MTtbl7927uvfdebr31VqZOnUqhQoXo3r07O3bsAKBhw4Y8/PDDCc8fHR3Nq6++mrAdVdKesePH\njzNv3jyaNGnC008/zdixYzly5AgtWrTgq6++uqS/GxNEqpolDqADoNHR0WqMMaEoJiZGAY2JiXE7\nlHQxf/58FRGfR65cuRLqffPNN5ojRw7t3bu3Hj16VEuWLKm1a9fW2NjYhDobN25UEdHSpUvryZMn\nE8qXL1+uIqIzZsxQVdVz585psWLFtGrVqnr27NmEeu+++66KiI4ZMyahrFu3bhoWFqYTJkxIFHeN\nGjX0hhtuSPj84Ycfqojo0qVLE9Vbt26dioguWbIkoaxs2bIaFhammzdvTij7/fffNWfOnDp48OCE\nshUrVmhYWJhu2rTpgj+3xo0ba5MmTRI+x8XF6blz5xLVOXbsmBYvXlwfeOCBROUiomPHjr2gTW+p\n+XcXHR2tgAIdNAS+0zPDYS+WjTEmgzp1CnbuTN97VKoEkZGBa09EmDlzJldddVWi8vDw8ITfr7nm\nGsaOHcuwYcPYvn07f/75J+vXr/c5E7Br165EegXYtm1bSpQowerVq3nooYf47LPPOHz4MOPGjSN7\n9uwJ9Vq2bEmlSpV49913GT16dKI2+/Tpk+hzgwYNiI6OTvi8YsUKChQoQLNmzfjjjz8SyqtXr06e\nPHnYsGED999/f0J5lSpVqFevXsLnwoULU7FiRfbs2XPRPy9fRCRhXJiqcvToUWJjY6lVqxbbtm3z\nq00TfJaAGWNMBrVzJ9Ssmb73iImBQI+Xv+GGGy46CH/w4MEsXbqUzz77jIkTJ1KxYkWf9XzNnLzy\nyivZu3cvAPv27UNEuPrqqy+oV6lSJTZv3pyoLGfOnFx22WWJygoWLMhff/2V8Hn37t0cPXqUokWL\nXtCmiHD48OFEZUlnNfpqM60WLFjA1KlT2blzJ+fOnUsoL1++vN9tJic2FhYuDHizWZ4lYMYYk0FV\nquQkSOl9Dzf8+OOP7N69GyBh4H0wePfEJScuLo5ixYqxePHi+CEuiRQpUiRVbfq6NjWio6Pp3r07\nd999N0OGDKFo0aKEh4czceJEv3vVUjJzJqxbF/BmszxLwIwxJoOKjAx871QoUFW6detG/vz5efTR\nR5kwYQJt27blzjvvvKBufJLm7YcffqBq1aoAlClTBlVl165dNG7cOFG9Xbt2UaZMmTTHV6FCBdav\nX0+9evXIkSNHmq/3JbklKHx5/fXXqVChAitWrEhUPmrUqIDE4u3nn2HYMLj5ZnjvvYA3n6XZLEhj\njDEhJSoqiq1bt/LSSy8xbtw46tWrR79+/RJmSnpbuHAhf//9d8Ln5cuXc+DAAVq2bAlArVq1KFq0\nKC+++GKiV3Vr1qxhx44d3H777WmOr127dpw/f97nFj+xsbEcO3YszW3mzp07YTzXxfjqUfvkk0/Y\nsmVLmu+bElXo2xcKFoT77gto0wbrATPGGBNEqsrq1asTlmDwVq9ePc6cOcOoUaPo3r17QhI1f/58\nqlWrRr9+/Vi2bFmiawoVKsRNN91E9+7dOXjwINOmTePqq6/mgQceAJxFTKdMmUKPHj1o2LAh7du3\n5+DBg0yfPp3y5cszcODAND9Dw4YN6dOnD5MnT+bLL7/k1ltvJSIigu+//54VK1Ywffp07r777jS1\nWa1aNcLDw5kyZQpHjx4lR44cNGvWjMKFC19Q9/bbb+eNN97gzjvvpFWrVuzZs4fZs2dzzTXXJEpG\nL9XixfDf/8Lbb4MfOaW5iJBJwETkQeBxoDiwHRigqp+lUL8jMBi4CjgGrAEGq+qF/4tkjDEmJIjI\nBbMO482dO5cXX3yRokWL8uyzzyaUX3nllUyaNImBAweyYsUK2rZtm9DW8OHD+eqrr5g8eTInTpzg\nlltu4YUXXiBnzpwJ13ft2pXcuXMzefJkhg4dSu7cubnnnnuYPHky+fLluyC+5OL2NmvWLGrVqsXs\n2bN58sknyZYtG2XLlqVLly7Ur18/0XWpabNYsWLMnj2bSZMm8cADDxAbG8uGDRto2LDhBXW7devG\noUOHmD17NuvWraNKlSosWrSI1157LWHx1tTcPyW//w6PPAL33w+33w6LFqW5CXMxbq+D4RmEeB9w\nBugCVAJmA38ChZOpXx84DzwIlAHqAV8DK1K4h60DZowJaZl9HbBAil8H7PXXX3c7lAzP17+7Dh1U\nCxVSPXTI+WzrgAX+CJUxYI8Cs1V1oaruBPoCp4AeydSvA/ykqi+o6s+q+jFO0nZjcMI1xhhjMqfV\nq53Xj88+Cz5W2jAB4noCJiIRQE1gfXyZqirwHlA3mcu2AKVF5DZPG8WAe4F30zdaY4wxJvM6ccIZ\neH/rrdC5s9vRZG6uJ2BAYSAcOJSk/BDOeLALeHq8OgHLROQf4ADwF/BQOsZpjDEmhPgztsmkbPhw\n+OMPmD0b7I83fYXMIPy0EJEqwDRgDLAOKAH8B+c15APuRWaMMSYYGjVqlGjTa3Pptm+HF16AqVOh\nbFm3o8n8QiEBOwLEAsWSlBcDDiZzzVBgs6pO9Xz+RkT6Ax+KyJOqmrQ3LUFUVNQF05jbt29P+/bt\n/QreGGOMyQzGjYMbb4QBA9yOJGtwPQFT1XMiEgM0A1YBiNOv3AyYnsxlkcA/ScricGZopNhp2rxH\ncyY9NOmSYjbGGGMym19+cdb8SsVuTCYAQmEMGMBUoJeIdBGRSsCLOEnWfAARmSQiC7zqvw3cIyJ9\nRaSciNTHeSX5iaom12sGwLKflsUvS2GMMcYYj+7d4dpr3Y4i6wiJBExVX8NZhHUc8AVwPdBcVX/3\nVCkOlPaqvwB4DGcdsK+BZcAO4J6L3eunv3/ijR1vBDR+Y4wxJqPr2dPtCLIW119BxlPVmcDMZM51\n91H2AvBCWu9zXcHrGLFhBG0qtSFbWMg8vjHGGOOq7NndjiBrCYkesGBqV7YdO4/sZOH2hW6HYowx\nxpgsKsslYOXylqPdNe0YvXE0Z86fcTscY4wxxmRBWS4BAxjfZDwHThxg5mc+33gaY4zJQLp160a5\ncuUSlZ08eZIHHniAEiVKEBYWxmOPPeZSdMb4liUTsKsvu5qe1Xsy8cOJHD973O1wjDEmSxgzZgxh\nYWH8+eefPs9fe+21NG3aNM3tighhYYm/ziZMmMDChQt58MEHiY6OprPtq2NCTJYdhT6q0SgWfrWQ\nqI+jGNtkrNvhGGNMpiciKW4f5O/WQnPnziUuLi5R2YYNG6hTpw4jRozwq01j0luW7AEDKJmvJANu\nHEDUligOnzzsdjjGGGP8FB4eTkRERKKyw4cPU6BAgYDdIzY2lnPnzgWsPWOybAIGMPSmoWQLy8aE\nDya4HYoxxhgvmzZtIiwsjOXLlzNhwgRKly5Nrly5uPnmm/nxxx8T1fUeAxZ/3d69e3nnnXcICwsj\nPDycffv2AfD777/Ts2dPihcvTq5cuahWrRoLFyaeFf/zzz8TFhbG1KlTmTZtGldeeSU5c+Zkx44d\nieIaO3YspUqVIl++fNx7772cOHGCf/75h4EDB1KsWDHy5s1Ljx49LHEzPmXZV5AAhXIVYkj9IYzd\nNJZH6z5K2QJl3Q7JGGOMl8mTJxMeHs7gwYM5duwYU6ZMoVOnTmzZsiWhjverzcqVKxMdHc3AgQMp\nXbo0gwYNAqBIkSKcOXOGRo0asWfPHgYMGEDZsmVZvnw53bp149ixYwxIsgnivHnzOHv2LH369CFH\njhwUKlSIv/76C4BJkyYRGRnJsGHD+OGHH5gxYwYRERGEhYVx9OhRxo4dy9atW1mwYAHly5e3V6Hm\nAlk6AQN4pPYjTP9kOmM2jmH+nfPdDscYY4yXs2fPsn37dsI9GxQWKFCAgQMH8t1331GlSpUL6hct\nWpQOHTrw5JNPUrJkSTp06JBwbtq0aezatYtFixZx//33A9C3b18aNmzIiBEj6NGjB7lz506o/+uv\nv/Ljjz9SqFChhLL43rfY2Fg2bdqUENfhw4dZunQpt912G++8805C27t372bevHmWgJkLZPkELHf2\n3IxsOJIBawbweL3HubaobYRljMkYTp07xc4jO9P1HpUKVyIyIjJd75GSHj16JCQ5AA0aNEBV2bNn\nj88ELCVr1qyhePHiCckXOOPHHn74YTp06MCmTZto2bJlwrm2bdsmSr68de3aNVFctWvXZunSpfTo\n0SNRvdq1azNjxgzi4uIumKlpsrYsn4AB9KrZi6gtUYx4fwRv3v+m2+EYY0yq7Dyyk5pzaqbrPWJ6\nx1CjRI10vYe3pDMhS5cunehzwYIFARJeBabFzz//zFVXXXVBeeXKlVFVfv7550TlZcuWTbatpHHl\nz58/2fK4uDiOHTuWELsxYAkYANnDszO+yXg6rezE1v1bqVOqjtshGWPMRVUqXImY3jHpfo9AyZkz\nJwCnT5/2ef7UqVMJdeJ59zJ5U9WAxZWcXLlyJXsuubjcjNdkLJaAebS/rj1TNk9h6HtD2dB1g9/r\n0RhjTLBERkQGtXfqUpUpUwaAXbt2UbJkyUTnTp8+zS+//ELz5s3T9f5ff/31BeU7duxIFJ8xwWAv\npD3CJIyJzSay6edNrPtxndvhGGNMptOsWTMiIiKYNWvWBT1Cs2fPJjY2NtEYrEBr2bIlBw8eZNmy\nZQllsbGxzJgxg7x589KoUaN0u7cxSVkPmJdWV7Wifun6DFs/jFsq3EKYWH5qjDGBUqRIEUaNGsXI\nkSNp2LAhrVu3JjIyks2bN7N06VJatGjB7bffnm737927N7Nnz6Zbt258/vnnCctQbNmyhWnTpiWa\nAekPe81o0sISMC8iwuSbJ9PglQas+G4F7a5p53ZIxhiTqQwfPpxy5crx/PPPM378eM6fP0+5cuUY\nP348Q4YMSVQ3uaEgvsqTlvna9ihnzpxs2rSJoUOHsnDhQo4fP07FihWZP3/+BXtFprRtUlrLjfFF\nskrGLiIdgEXR0dF07NgxxbqtFrdi9x+7+bb/t0SER6RY1xhjAmXbtm3UrFmTmJgYatTIOGO7TMaW\nmn93ixYtolOnTgAdVXVxUAPMpOwdmw8Tm05k95+7eeXLV9wOxRhjjMkQRORBEflJRE6LyFYRueEi\n9RuLSIyInBGR70Wkq48694rIDk+b20XktiTnh4nIpyJyXEQOichKEbnaRzvjROQ3ETklIv8TkSsv\n/YkvjSVgPlQtXpUO13Vg7KaxnD7ne7q0McYYYxwich8QBYwGqgPbgbUiUjiZ+mWBd4D1QFVgGjBX\nRG7xqlMPWAy8BFQD3gLeFBHvFXgbADOA2sDNQASwTkRyebXzBPAQ0Bu4ETjpiS37pT73pbAELBnj\nGo/j8MnDPP/p826HYowxxoS6R4HZqrpQVXcCfYFTQI9k6vcD9qjqEFXdpaovACs87cR7GFijqlM9\ndUYB23CSKQBUtaWqvqqqO1T1a6AbcAXgvULxI8B4VX1HVb8BugCXA3cG4Ln9ZglYMioUqkCvGr2Y\n9NEkjp456nY4xhhjTEgSkQichGd9fJk6A8zfA+omc1kdz3lva5PUr5uKOkkVABT40xNbOaB4ktiO\nA59cpJ10ZwlYCkY2HMmZ82d4ZvMzbodijDHGhKrCQDhwKEn5IZzkx5fiydTPJyI5LlLHZ5viTEN9\nDvhIVb/zakPTGFtQWAKWghJ5SzCwzkCe++Q5Dv590O1wjDHGGJO8mUAV4P6LVQwFtg7YRQypP4QX\nP3+Rpz54iudb2ngwY4wxmduSJUtYsmRJorL9+/endMkRIBYolqS8GJBc78XBZOofV9WzF6lzQZsi\n8jzQEmigqgeS3Ec813n3ghUDvkgmtqAImR6wtExfFZFXRCRORGI9P+OPCzf5ukQFchbgifpPMDtm\nNnv+2hPo5o0xxpiQ0r59e1atWpXoGDRoULL1VfUcEAM0iy/zvA5sBnyczGVbvOt73OopT6nOLUnq\nxCdfbYAmqrovSWw/4SRh3rHlw5k1mVxsQRESPWBe01d7A5/izIJYKyJXq+oRH5c8DDzh9Tkb8BXw\nWnrEN6D2AKZ9Mo1RG0YRfXd0etzCGGMSxG8ObUwwBOjf21RgvojE8O/3eCQwH0BEJgGXq2r8Wl8v\nAg+KyBRgHk6C1BanFyveNGCjiDwGvAu0xxns3yu+gojM9JS3Bk6KSHyP2TFVPeP5/TlghIj8AOwF\nxgP7cZa1cE1IJGB4TV8FEJG+QCuc6atPJ62sqieAE/GfReROnJkP89MjuMiISEY3Gk2/d/sxpP4Q\nri92fXrcxhiTxRUuXJjIyMj4FceNCZrIyEgKF/a5ZFeqqOprnjW/xuG83vsSaK6qv3uqFAdKe9Xf\nKyKtgGdxOlX2Az1V9T2vOls8u9hM8By7gTZeA+zBWe5CgY1JQuoOLPS087SIRAKzcXKFD4HbVPUf\nvx84AFzfisgzffUUcI+qrvIqnw/kV9W7UtHGKiC7qrZIoU6qtyLy5VzsOarMrEKlwpV4u/3bab7e\nGGNSY9++fRw54qvj35h/RUfDs8/C/Plw3XWX3l7hwoW54oorkj1vWxEFXij0gKU0fbXixS4WkRLA\nbaTzrIeI8AieavIU979+Px/t+4ibrrgpPW9njMmirrjiihS/CI3ZswdefBEefhi6XrB5j8koQiEB\nu1TdgL9I5bvcqKgoli1blqisffv2tG/f/qLX3nvNvUzePJmh7w3lw+4f4owxNMYYY4JDFfr0gaJF\nYcIEt6MxlyIUEjB/pq966w4sVNXzqbnZoEGD/HoFCRAmYUxqNonbFt3G6t2raXV1K7/aMcYYY/yx\nYAG89x6sWQN58rgdjbkUri9D4ef01fh6jYEKwMvpGGIizSs0p1GZRgx/fzhxGhes2xpjjMniDh2C\nxx6DTp2gRbIjnk1G4XoC5jEV6CUiXUSkEs701ETTV0VkgY/regKfqGrQ5myLCJOaTeKrQ1+x9Jul\nwbqtMcaYLG7AAAgPdwbfm4wvJBIwVX0NeBxn+uoXwPWkMH0VEhZSuwuYG8RQAahbui6tK7Zm5IaR\n/BPr6ixWY4wxWcBbb8Hy5TB9OlzCahEmhIREAgagqjNVtayq5lLVuqr6ude57qraNEn946qaR1Xn\nBT9amNB0Aj/99RNztwU9/zPGGJOFHDsG/ftDy5Zwf4bY5dCkRsgkYBnNtUWvpXPVzoz/YDwn/znp\ndjjGGGMyqSeegOPHYdYssMn3mYclYJdgbOOx/HHqD6Z/Mt3tUIwxxmRCH3wAs2fD5Mlgy8NlLpaA\nXYKyBcoxMiNkAAAgAElEQVTSt1Zfpmyewp+n/3Q7HGOMMZnImTPQqxfUqwf9+rkdjQk0S8Au0ZMN\nnuR83HmmfDTF7VCMMcZkIuPGwd69MHcuhNm3daZjf6WXqFieYjxW9zGmfzqdX4//6nY4xhhjMoHt\n2+Hpp2HECKhc2e1oTHqwBCwABtUdRO6I3Iz/YLzboRhjjMngzp+Hnj2dxOuJJ9yOxqQXS8ACIH/O\n/Ay7aRhzt81l9x+73Q7HGGNMBvbcc7Btm/PqMXt2t6Mx6cUSsADpf0N/SuQtwcgNI90OxRhjTAb1\n448wahQ88gjUru12NCY9WQIWILkicjGm0RiWfbuMLw584XY4xhhjMhhV6N0bihaF8TaiJdOzBCyA\nulbrSsXLKjL8/eFuh2KMMSaDeeUVeP99mDMH8uRxOxqT3iwBC6BsYdl4qulT/PeH/7Jx70a3wzHG\nGJNBHDgAgwZBly5w661uR2OCwRKwALun8j3ULFGTYeuHoapuh2OMMSYDePhhiIiAqVPdjsQEiyVg\nASYiTL55Mlv3b+Xt7992OxxjjDEh7s03YcUKmD4dLrvM7WhMsFgClg5uLn8zzco1Y/j64cTGxbod\njjHGmBB19Cj07w+33w733ed2NCaYLAFLJxObTeTb379l0deL3A7FGGNMiHriCfj7b5g1C0TcjsYE\nkyVg6eTGkjdyd+W7GbVhFGfPn3U7HGOMMSFm0yZnxuPkyVCqlNvRmGCzBCwdPdXkKX45/gtzYua4\nHYoxxpgQcvo0PPAA3HQT9O3rdjTGDZaApaPKRSrTrWo3xn8wnhNnT7gdjjHGmBAxbhzs2wcvvQRh\n9k2cJdlfezob3Xg0x84e47mtz7kdijHGmBDwxRfwzDMwciRUquR2NMYtloClsyvyX8GDNzzIMx8/\nw5FTR9wOxxhjjIvOn3dePVapAkOGuB2NcZMlYEEw7KZhAEz+aLLLkRhjjHHTs8/Cl1/C3LmQPbvb\n0Rg3WQIWBEVyF+Hxeo/z/KfP88uxX9wOxxhjjAt++AFGjYJHHoEbb3Q7GuM2S8CC5NE6j5IvRz7G\nbhrrdijGGGOCTBV694YSJWD8eLejMaEgZBIwEXlQRH4SkdMislVEbrhI/ewiMkFE9orIGRHZIyLd\nghRumuXNkZcnGzzJK1++ws4jO90OxxhjTBDNmwcbNsDs2ZA7t9vRmFAQEgmYiNwHRAGjgerAdmCt\niBRO4bLlQBOgO3A10B7Ylc6hXpK+tfpSKl8pRm4Y6XYoxhhjguTAARg0CLp2hVtucTsaEypCIgED\nHgVmq+pCVd0J9AVOAT18VRaRFkADoKWqblDVfar6iapuCV7IaZcjWw7GNR7Hiu9W8Nmvn7kdjjHG\nmCB46CHIkQOmTnU7EhNKXE/ARCQCqAmsjy9TVQXeA+omc9kdwOfAEyKyX0R2icgzIpIz3QO+RJ2u\n70SVIlUY/v5wt0MxxhiTzt54wzlmzIBChdyOxoSSbKmpJCJvpLZBVb07jTEUBsKBQ0nKDwEVk7mm\nPE4P2BngTk8bs4BCQM803j+owsPCmdB0Anctu4v1e9bTrHwzt0MyxhiTDo4edXq/7rgD7r3X7WhM\nqElVAgYcS9co0i4MiAM6qOrfACLyGLBcRPqrakjvft2mYhtql6zNsPXD+KTcJ4iI2yEZY4wJsMGD\n4e+/YeZMsP/Mm6RSlYCpavd0jOEIEAsUS1JeDDiYzDUHgF/jky+PHYAApYAfk7tZVFQUy5YtS1TW\nvn172rdvn8aw/SciTL55Mk0WNGHlzpXcXTmtnYbGGGNC2YYNzmKrs2ZBqVJuR2NCUWp7wNKNqp4T\nkRigGbAKQJwuoWbA9GQu2wy0FZFIVT3lKauI0yu2P6X7DRo0iI4dOwYk9kvRuGxjbq1wK0++/ySt\nK7YmW5jrfxXGGGMC4PRp6NULGjRw1v4yxhe/BuGLSFsRec2zXtc278PPOKYCvUSki4hUAl4EIoH5\nnvtNEpEFXvUXA38Ar4hIZRFpCDwNvBzqrx+9TWw6kZ1HdrJw+0K3QzHGGBMgY8bA/v3w0ksQ5vpU\nNxOq0vxPQ0QeBl7BGSRfHfgUJxkqD6zxJwhVfQ14HBgHfAFcDzRX1d89VYoDpb3qnwRuAQoAnwGv\nAm8Bj/hzf7fUvLwm7a5px5iNYzhz/ozb4RhjjLlE27ZBVJSz5VDF5KaRGYN/PWD9gd6qOgD4B3ha\nVW/BeV2Y399AVHWmqpZV1VyqWldVP/c6111Vmyap/72qNlfVPKpaRlWHZKTer3jjm4zntxO/Meuz\nWW6HYowx5hKcOwc9e8I11zgD8I1JiT8J2BXAx57fTwN5Pb+/irMavUmDqy+7mh7VezDhwwkcP3vc\n7XCMMcb4aepU+OorePlliIhwOxoT6vxJwA7irLcFsA+o4/m9HM4sRJNGoxqN4uS5k0R9HOV2KMYY\nY/ywe7cz9uvRR6FWLbejMRmBPwnY+0Brz++vAM+KyP+AZcDKQAWWlZTKV4oBNw5g6tapHD552O1w\njDHGpEFcnDPrsUQJGDvW7WhMRuFPAtYbmACgqi/g7Ne4AxgF9AtcaFnL0JuGEi7hTPxwotuhGGOM\nSYOXX4ZNm2DOHMid2+1oTEaR5gRMVeNU9bzX56Wq+rCqzlDVfwIbXtZRKFchBtcbzKzPZ7H36F63\nwzHGGJMKv/3mDLjv3h1uvtntaExGkqoETESuF5Ewr9+TPdI33MztkTqPUDBnQcZsHON2KMYYY1Lh\noYcgZ074z3/cjsRkNKntAfsSZ8Pr+N+/8PxMenwR6ACzkjzZ8zCy4Uhe/epVvj38rdvhGGOMScGc\nObByJcyYAYUKXby+Md5Sm4CVA373+r2852fSo3ygA8xqetXsRZn8ZRixYYTboRhjjEnGBx/Agw9C\n//5w771uR2PSm4hcJiIviMh3InJERP70PvxpM7Wbcf/s9bEM8LH3ODBPcNmAeoB3XZNG2cOzM67J\nODqv7MzW/VupU6rOxS8yxhgTNHv3wj33OHs9Pvec29GYIHkVuBJ4GWcnIL3UBv3ZAXoDUAJIul5C\nfs+58EsNKqtrf217nt78NMPWD+P9Lu/j7E1ujDHGbSdOQOvWkC8fLF9uC65mIQ2Am1R1e6Aa9GcZ\nCsF35ncZcPLSwjEA4WHhTGw2kY17N/K/Pf9zOxxjjDE463117uz0gK1aBZdd5nZEoUVEHhSRn0Tk\ntIhsFZEbLlK/sYjEiMgZEfleRLr6qHOviOzwtLldRG5Lcr6BiKwSkV9FJE5EWvto4xXPOe9jdRof\nbyeQK43XpCjVCZiIvCEib+AkX/PjP3uOt4C1/LtFkblEra5qRf3S9Rm2fhhxGud2OMYYk+WNHu0k\nXosXO/s9mn+JyH1AFDAaqA5sB9aKSOFk6pcF3gHWA1WBacBcEbnFq049YDHwElANeAt4U0SqeDWV\nG2cSYH9Sfi24BigGFPccad06sT8wQUQaecaD5fM+0tgWkLYesGOeQ4ATXp+P4WxPNAfo5E8Q5kIi\nwqRmk9h2YBsrvlvhdjjGGJOlLVsGTz0FkybB7be7HU1IehSYraoLVXUn0Bc4hbNYuy/9gD2qOkRV\nd3kWdl/haSfew8AaVZ3qqTMK2AY8FF9BVf+rqqNU9S1S3g7xrKr+rqqHPcexND7fUSAfzm5Ah4G/\nPMdRz880S/UYMFXtDiAie4H/qKq9bkxnDco0oOVVLRnx/gjuqnQXEeE22MAYY4ItJga6dYOOHWHI\nELejCT0iEgHUBBK2clFVFZH3gLrJXFYHeC9J2VrgWa/PdXF61ZLWaeNHmI1F5BBOsvQ+MEJV0zJ7\ncRFwDuiAW4PwVdV2ugqiiU0nUm12NeZ/OZ9eNXu5HY4xxmQpBw5AmzZw3XXw0ktgc6J8KowzAe9Q\nkvJDQMVkrimeTP18IpJDVc+mUKd4GuNbA7wO/ARUACYBq0WkrqqmNpG6FqiuqrvSeO9kpXkQvogU\nE5FXReQ3ETkvIrHeR6ACM46qxavS4boOjNk0htPnTrsdjjHGZBlnzsBdd4EqvPkm5AroEGwTLKr6\nmqq+o6rfquoq4HbgRqBxGpr5HCgdyLj8WYZiPnAFMB44QAC64UzKxjUeR6UXKvH8p88zuP5gt8Mx\nxphMTxX69IHt251FVy+/3O2IgmfJkiUsWbIkUdn+/ftTuuQIEIszyN1bMZwx4r4cTKb+cU/vV0p1\nkmszVVT1JxE5grOu14ZUXjYDmCYizwBf47yO9G7zq7TG4U8CdhPQQFW/9ONa44cKhSrQq0YvJn00\niV41e1EgZwG3QzLGmEwtKgoWLoRFi+CGFBdTyHzat29P+/aJJwkuWrSITp18z7NT1XMiEgM0A1YB\niLOAZTNgejK32QLclqTsVk+5d52kbdySpE6aiUgpnKWzDqThsmWen/O8ypR/l+ZK8xqo/qwD9gsp\nzzQw6WBkw5GcOX+G/3xsO74aY0x6Wr3aGWw/dCh06OB2NBnGVKCXiHQRkUrAi0AkzlszRGSSiCzw\nqv8iUF5EpohIRRHpD7T1tBNvGtBCRB7z1BmDM9j/+fgKIpJbRKqKSDVPUXnP59Je558WkdoiUkZE\nmgFvAt/jDOhPreS2X/R7G0Z/ErCBwGTPGh4mSErkLcHAOgN5duuzHPz7knpfjTHGJGPHDmjf3llq\nYsIEt6PJOFT1NeBxYBzwBXA90FxV4/eRLo7XGCpV3Qu0Am7GWcfrUaCnqr7nVWcLzqzD3p46dwNt\nVPU7r1vX8twvBqcnKgpnqYr4CYOxnljeAnbhrCn2GdBQVRO9RkyOZ5bnaCBMVX/2daSmnQvaTf0E\ngIRA/sLJarPhrPGR9D1oSO4JLyIdgEXR0dF07NjR7XD8cvTMUcpPK0+H6zrwfMvnL36BMcaYVPvz\nT6hdG3LkgC1bIG9etyMKHV6vIDuq6mK34wk2ETkGVFPVnwLVpj9jwAYG6uYmbQrkLMAT9Z9gxIYR\nPFb3McoX9KvX0xhjTBLnz8N99zlJ2GefWfJlLvAmcCeJ1ym7JP6sA7bg4rVMehlQewDTPpnG6I2j\nefWuV90OxxhjMoXHHoONG2HdOihv/29rLrQbGCUi9XFedyZajF5Vk5tskCx/xoAhIhVE5CkRWSIi\nRT1lt4mI7Y6VziIjIhndaDSLvlrEV4fSPOvVGGNMEi+9BDNmwPTp0KSJ29GYENUTZ9uhmjhj0h71\nOvx6M+jPQqyNcNbAqI0zIC6P51RV/h30lmZp2UXdsxlm0p3NY+OTwcyuR/UeVChUgSfff9LtUIwx\nJkP78EN48EHo1885jPFFVculcARtFuRknD2UbgH+8Sp/H2dvpzRL6y7qHgpcxb87m5dQ1cP+3D+j\niQiPYHyT8bzz/Tts3rfZ7XCMMSZD2rsX7r4b6teHadPcjsZkFOJxqe34k4BdB6z0UX4YZz8of6R1\nF/V43jubZ4nkK167a9pRrXg1hq4fSlpnshpjTFb399/OHo9588Ly5RAR4XZEJtR51jj7GjgNnBaR\nr0Sks7/t+ZOAHQVK+CivDvya1sa8dlFfH1/m2RwzpV3UwVkM9kvPnpTrRKReWu+dkYVJGJOaTeKj\nfR+x5oc1bodjjDEZRlwcdOkCe/bA229DYX+7DkyWISKPAbOA1UA7z/Ff4EURedSfNv1JwJYCU0Sk\nOM5rwDDPrID/AAv9aC+lXdST2/H8ANAHuAdnHNovwEavlXCzhOYVmtOoTCOGrR9GnMa5HY4xxmQI\nY8Y4m2svXgzX2NQxkzoDgH6q+oSqrvIcQ4D+wMP+NOjPOmDDgRdwkp5w4DvPz8XAU/4EkVaq+j3O\nNgLxtopIBZxXmV1TujYqKoply5YlKvO171VGICJMajaJevPqsfSbpXS4zvbMMMaYlLz2GowfD5Mm\nwR13uB2NyUBKAB/7KP8Y328FL8qfdcD+wdnvaTxwLc4syC9Udbc/AeDfLuq+fArUv1ilQYMGZdiV\n8H2pW7ourSu2ZuSGkbSt0pbs4dndDskYY0LStm3QrZuzv+MTT7gdjclgfsB57TgxSfl9OGuEpZk/\nPWAAqOo+YJ+/13u1488u6r5UI207m2caE5pO4PpZ1/Pytpfpd4PNozbGmKQOHnQG3V97LcydC5c+\nh81kMaOBZSLSEIhffqA+Tq7Szp8G05yAiUg40M1z06IkGUemqk39iGMqMN+TiH2K8yox0S7qwOWq\n2tXz+RHgJ+BbICfQC2gC3OLHvTO8a4teS+eqnRn3wTi6VO1C7uy53Q7JGGNCxtmzznITsbGwciXk\nyuV2RCajUdXXRaQ2Tn5yp6d4B3Cjqn7hT5v+DMKf5jnCgW9w1uzyPtIsrbuoA9lx1g37CtiIszRG\nM1Xd6M/9M4Oxjcfyx6k/mP5JmndDMMaYTEsV+vRxXj+++SaULOl2RCajEJGpIpLb83tDYLuqdlLV\nmp6jk7/JF/j3CvJ+oJ2qrvb3pr6o6kxgZjLnuif5/AzwTCDvn9GVLVCWvrX68vTHT9O3Vl8K5iro\ndkjGGOO6qVNhwQKIjoYbb3Q7GpPBDACm4Oz7uAFnsH3A1hz1pwfsH5zBaCbEPNngSc7HnafV4lbs\nP77f7XCMMcZVa9bAkCHOgPtMNPfKBM9e4GHPFowC1BWRhr4Ofxr3JwGLAh4JxDL8JrCK5SnGe53f\n45fjv1Bjdg02/LTB7ZCMMcYVO3fC/fdDy5YwYYLb0ZgMajDOJtwbcNY9XYkz7Cnp4deXrT8J2E1A\nR+BHEXlbRN7wPvwJwgRO7VK12dZ7G1WLV+XmV29mykdTbKsiY0yW8tdf0Lo1lCoFixZBeLjbEZmM\nSFXfVNXiQD6cHrCKQEEfRyF/2vdnDNhRfO8FaUJEkdxF+G/H/zJ642iGrh/K1l+3Mr/NfPLnzO92\naMYYk67On4d27eCPP+DTTyFfPrcjMhmdqv4tIk2An1T1fKDa9Wch1u4Xr2XcFh4WzlNNn6J2ydp0\nXtmZWi/V4vV2r3N9sevdDs0YY9LN44/Dhg2wbh1UqOB2NCYjE5F8qnrc8/ELIDK50Vde9VLNn1eQ\n8YEVEZGbPEcRf9sx6euOincQ0zuG3BG5qTO3DtFfRbsdkjHGpIuXX4Zp02D6dGjqz4qUxiT2l4gU\n9fx+FPjLxxFfnmb+LMSaG5gBdOHfBC5WRBYCA1T1lD+BmPRToVAFPu75Mf3f7U/nlZ3Z8ssWpjaf\nSo5sOdwOzRhjAuKjj6BfP+jbF/r3dzsak0k0Bf70/N4k0I37MwZsKtAIuIN/l+O/CWfboCjA9sIJ\nQZERkbzS5hXqla7HgDUDiDkQw/J7l1M6f+mLX2yMMSHs55+dle7r1XN6v4wJBFXd5Ov3QPEnAbsH\naJtk1fnVInIaeA1LwEKWiNC7Zm+qF69O2+VtqTGnBkvuWcLN5W92OzRjjPHL3387Mx7z5IEVKyAi\nwu2ITGYlIgWAG/G9DePCtLbnTwIWCRzyUX7Yc86EuBtK3kBM7xg6vtGR5tHNGd9kPENvGkqY+D0k\n0Bhjgi4uDrp2hT17YMsWKFzY7YhMZiUidwCLgDzAcZx1weIpkOYEzJ9v3C3AWBHJ6RVYLpydwrf4\n0Z5xQeHIwqzusJoRDUbw5PtPcufSOzl65qjbYRljTKqNHetsrr1oEVx7rdvRmEwuCpgH5FHVAqpa\n0Ovwax0wfxKwR4D6wH4RWS8i64FfgHqecyaDCA8LZ2yTsbzb4V0+3PchNefU5MuDX7odljHGXNTy\n5TBuHDz1lPMK0ph0VhKYHsiJhmlOwFT1G+AqYBjwpecYClylqt8GKjATPC2vasm23tvInyM/dV+u\ny4IvF7gdkjHGJOuLL5xXj+3bw7Bhbkdjsoi1QK1ANujPGDA8GeBLgQzEuKtcwXJs7rGZh1Y/RLe3\nurFl/xamtZhmS1UYY0LKoUPQpg1cc42z7pftSmyC5F3gGRGpAnwNnPM+qaqr0tqgXwmYiFQEBgCV\nPUU7gOdVdac/7ZnQkCsiFy+3eZm6pevy0OqHiDkQw4p7V1CmQBm3QzPGGM6ehbvugnPn4M03IVcu\ntyMyWUh8p9MoH+cUSPOOo2l+BSki9wDfADWB7Z6jBvC155zJ4B6o8QCbe2zmyKkj1JhTg3U/rnM7\nJGNMFqfqLLK6bZuTfJUs6XZEJitR1bAUDr+2e/dnEP7TwCRVrauqj3mOesBEzzmTCdS8vCYxvWOo\nXbI2LaJbMH7TeOI0zu2wjDFZ1HPPwfz5MHcu1K7tdjTGXDp/XkGWwPd6F9HA4EsLx4SSQrkK8U6H\nd3jqg6cYvXE0W3/dyqt3vUqhXH7NuDXGGL+sXetssj1kCHTq5HY0JisSEV+vHhOo6ri0tulPArYR\naAD8kKT8JuBDP9ozISxMwhjVaBQ3lryRjm90pOacmrze7nVqlKjhdmjGmCxg1y647z647TaYONHt\naEwWdleSzxFAOeA88CMQlARsFTBFRGoCWz1ldYB7gdEikrAiiz+zAkxoanFlC2J6x9D2tbbUe7ke\nM1vNpEf1Hm6HZYzJxP76C+64Ay6/HBYvhnC/RtoYc+lUtXrSMhHJB8wHVvrTpj8J2EzPz/6ew9c5\n8HNWgAldZQuU5aMeH/HImkfouaonW37ZwoyWM8iZLefFLzbGmDQ4f97p+TpyBD79FPLlczsiYxJT\n1eMiMhp4G3g1rdf7sxBrSjMBLnlWgAltObPlZPYds5nXeh7RX0dTf159fvrrJ7fDMsZkMoMHw/vv\nOxtsX3ml29EYk6z8niPN/FoHLDkiEhnIZfpN6OpevTvVilfjntfuoeacmiy6exG3XXWb22EZYzKB\nefOcWY/PPw9Nm7odjTEgIg8nLcKZlNgZWONPm/6sA7ZeRC5YgUVEauNsS+QXEXlQRH4SkdMislVE\nbkjldfVF5JyIbPP33sY/1UtUJ6Z3DPVK16PV4laM2TjGlqowxlySjz5y1vvq3Rv6Jx3kYox7Hk1y\nPAw0BhYAffxp0J91wM4AX4nIfQAiEiYiY3BmQK72JwhPW1HAaKA6zuKua0Wk8EWuy4/z8O/5c19z\n6QrmKsiq9qsY12Qc4zaNo9XiVvxx6g+3wzLGZED79sHdd0PdujBjhm0zZEKHqpZLclRQ1TqqOlxV\nT/jTpj9jwFrhLMU/T0QWAx8BvYDbVXWgP0HgZJOzVXWhZzujvsAp4GLT7F4EFvHvbEzjgjAJY0TD\nEfy303/57NfPqDmnJp//9rnbYRljMpCTJ6F1a8idG15/HbJndzsiY9KXPz1gqOoLwHTgfpzdwe9V\nVb/2qxGRCJxtjdZ7ta84vVp1U7iuO84aHGP9ua8JvFsr3Mq2Ptsomrso9efV56WYl3D+Ko0xJnlx\ncdC1K/zwA6xaBYVTfPdhTObgzxiwgiLyOtAP573na8A6EfH3bX1hnOUqDiUpPwQUTyaGq3C2Puqo\naoOOQskV+a/gw+4f0qNaD3q/05ueq3py+txpt8MyxoSwceOcXq/oaLjuOrejMSY4/OkB+wYoBlRX\n1ZdUtRPQExgvIu8GNDofRCQM57XjaFX9Mb44ve9rUi9HthzMun0WC+5cwJJvllBvXj32/LXH7bCM\nMSFoxQoYOxaeegruvNPtaIwJHn8SsBeBhqqasPiTqi4DqgL+vLU/AsTiJHXeigEHfdTPi/Pa83nP\n7MdzwEigmoj8IyKNU7pZVFQUrVu3TnQsWbLEj7DNxXSp2oWtPbdy4uwJas6pyTvfv+N2SMaYEPLl\nl86rx/vvh+HD3Y7GmOSJSMBXHE/zOmCqOj6Z8v3ALX60d05EYoBmONscISLi+TzdxyXHgWuTlD0I\nNAHuAfamdL9BgwbRsWPHtIZp/FS1eFU+7/05Xd/syh1L7mBEgxGMaTyG8DBbp9eYrOzQIWfQfeXK\n8PLLNuPRhLzDIvIGzhu49YEY/uTXIHwRaSAi0SKyJX5NMBHpLCI3+RnHVKCXiHQRkUo4vWyROHss\nISKTRGQBOAP0VfU77wM4DJxR1R2qagOOQkyBnAVYed9KJjadyMSPJnLbots4cuqI22EZY1xy9qyz\n3MS5c/DmmxAZ6XZExlxUVyA38Bbwq4g8JyK1LqVBfwbh3wOsBU7jrNmVw3MqP+BXJ7KqvgY8jrOb\n+BfA9UBzVf3dU6U4UNqftk1oCJMwhjUYxrpO6/jy4JfUmF2DT3/91O2wjDFBpgr9+kFMDKxcCaVK\nuR2RMRenqitV9V6c4VHDgSrAVhH5XkRG+dOmPz1gI4C+qtoLOOdVvhmo4U8QAKo6U1XLqmouVa2r\nqp97neuuqsluSKGqY1XV73ub4GlWvhnb+myjZL6S3DTvJl78/EVbqsKYLGTaNHjlFXjpJahTx+1o\njEkbVT2hqq+o6q04nUUncRaRTzN/ErCKwAc+yo8BBfwJwmQtpfKVYlO3TfSu2Zt+7/aj21vdOHXO\nthA1JrNbuxYGDXI22u7c2e1ojEk7EckpIu1E5E1gG1AIeMaftvxJwA4CvvamvwmwtQZMqmQPz87z\nLZ8n+q5oln+7nLov1+WHP39wOyxjTDrZtQvuuw9atIBJk9yOxpi0EZHmnrHoh4BZnp+3qmoZVR3q\nT5v+JGAvAdM8m28rcLmIdAT+4wnKmFTreH1HPnngE06fO02tObVYtWuV2yEZYwLs6FFnxuPll8Pi\nxRBuk6BNxrMSyAV0AYqrah9V9fU2MNX8ScAmA4txtg7Kg/M6ci7OXo4zLiUYkzVdV+w6Puv1GU3L\nNaXN0jYMXz+c83Hn3Q7LGPP/9u40uooq+/v4dydAgISAQBLmGUVGURFpJ2RSoVtanDpgazs2Imoj\nKq1/51mfBkHRdmwVRRRwQhBREBQRREGjqKACyiSTIIQZkvO8ODchhCQkl9xUcvP7rFUr91ZOVe0K\nrNydc07tUwz27fN1vjZs8MsMVa8edEQSKWZ2jZktN7OdZjbPzDodon1XM1tgZrtCk9kvyaPN+Wb2\nQ+icaWZ2Vq7vn2Jmk8xstZllmtnZ+VzrHjNbY2Y7zOxDM8trJK8gKc65C5xz7zjn9h66+aGFsxi3\nc7zcp9kAACAASURBVM7djx/3bAucCCQ5524vjoCkfKpeuTpvXPAGD/d4mIfnPMwZr5zB+u3rgw5L\nRA7TzTfD9OkwYQK0KOpHnpQZZnYhMBw/Ib0jkAZMM7M8V/Y0sybAZHxnTgdgFPCcmfXM0eZP+A6f\nZ4Fj8CUg3jaz1jlOFQ98DQzCj8rlda1hwGDgKuAE/MT5aWZW6OLxzrl0M2tuZveZ2TgzSw6d+ywz\na1PY8+QUVh2wUDB7QnW45jvntoV7HpEsZsbNJ93MjItnsGj9Io59+ljmrZoXdFgiEqb//Q8efdRv\n3bsHHY1E2BD8SNgY59xiYCCwA7gsn/ZXA8ucczc755Y4554AJobOk+U6YKpzbkSozR34ie+Dsxo4\n5953zt3hnHuH/JclvB641zk32Tm3CD+MWA8o9OJXZnYa8C3QGeiHHwEEnzzeXdjz5BR2AiYSKV2b\ndGXhVQtpXKMxp75wKqPnj1apCpEyZs4cGDgQrrwSBg8+dHspu8ysInAcvjcL8KNlwHSgSz6HnRj6\nfk7TcrXvUog2h4qtKb6WaM7YtgKfF+U8+OlXtznnegJ7cuz/CH8vRaYETEql+on1mXnJTK4+/mqu\nnXotf3/r72zfsz3osESkEJYt85Xuu3SB0aO1zFA5UBuIxT8ZmNM6fPKTlzr5tE80s7hDtMnvnPld\nxxXDedrhJ+Lnth5//0WmBExKrUqxlRh11ijGnTuOtxe/zYnPn8iPv/8YdFgiUoBFi+Dkk/1k+4kT\noVKhZ9mIlGp/AHXz2N8RWB3OCYu8GLeZVXbO7QrnYiLh+Fvbv9EuuR3njj+X4585npf++hLnHH1O\n0GGJSC5z50KfPtC4Mbz/PiQlBR2RhGPcuHGMGzfugH2rVq0q6JCNQAZ+mZ6cUvC1Q/OyNp/2W51z\nuw/RJr9z5ncdCx2XsxcsBb/0YWG9BjxsZufje9RizOwkfAmuMUU4T7ZwesDWm9mLZtbTzNSDJiWi\nTXIb5l85n17Ne9FvfD+GfThMpSpESpEPPoAePaBtW5g1C1Jyf2xKmZGamsqkSZMO2IYOHZpv+1BZ\nhgVA9qMWZmah95/lc9jcnO1DeoX2F9SmZ642BXLOLccnYTljS8RPps8vtrzcCiwGVuIn4H+PL8P1\nGXBfEc6TLZwEqthXBBcpjMS4RCacP4H/9PwPw+cOp/UTrRk1bxRbdm0JOjSRcm38ePjzn+H0033P\nl2p9lUsjgCvN7GIzawU8BVQFXgQwswdDleSzPAU0M7OHzewoMxsEnBc6T5ZRwJlmdkOozV34yf6j\nsxqYWbyZdTCzY0K7moXeN8xxnpHAbWb2FzNrh++xWoXPYwolVPnhSqA58GfgIqCVc+7vzrmMwp4n\np3DqgBX7iuAihWVmDP3TUD6/4nOOr3c8N354I/VH1OfqyVezaP2ioMMTKXeeecYXWr3gAnjrLaha\nNeiIJAjOufHAjcA9+KG99sAZzrkNoSZ1gIY52v8C9AF64Ot4DQEud85Nz9FmLtAfX7/ra3z5h77O\nue9zXPr40PUW4IcGh+NLVdyd4zyPAI8DT+OffqwCnOWcy/k0Y2Hvc4Vz7j3n3Hjn3E9FPT4nK47H\n+0NF0cYC7Z1zpXKRCTPrD4x95ZVXGDBgQNDhSDH5Lf03nlnwDE8teIq129bStUlXBncaTN9WfakQ\nU+QpjiJSSM7BQw/BrbfCtdfCyJEQo0kpUWvs2LFcdNFFAAOcc68GHU9JMLMRh27lOeduKOr5w/6E\nMrPKwNn47PRM/OS2sFYEFwlX3Wp1ubPrndxyyi289cNbjP5iNOdNOI8GiQ0YeNxArjzuSpLjk4MO\nUySqOAc33QTDh8Ndd8Edd6jUhESljrneH4vPm5aE3h+Jf/hgQTgnL/LfK5FYEVzkcFWKrcSFbS9k\n9qWz+eqfX3Fm8zO5f/b9NHy0IRe/dTHzV88POkSRqLBvH1x+uU++HnsM7rxTyZdEJ+fc6Vkb8C7w\nMdDAOXesc+5Y/JDqTGBKOOcPp8O42FcEFylOx9Q5hmfPfpZVN6zigW4P8OmKT+n8XGdOePYExqSN\nYdc+VVERCceuXXD++TBmDLzyih96FCknhgK3OOc2Z+0Ivb4t9L0iCycBK/YVwUUioWaVmgz901B+\nuvYnJqdOplbVWlzy9iU0fLQht864lRVbVgQdokiZkZ7ua3y9/z688w5oKq2UM4lAXpXtkoBq4Zww\nnKcgi31FcJFIio2Jpc+RfZg6YCpLBi9hQLsBPPHFEzQd1ZR+r/fjo+Ufaa1JkQJs3AjdusGXX/p6\nX336BB2RSIl7C3jBzPqZWYPQdi7wPPBmOCcMZw5Ysa8ILlJSjqx1JCPPHMnqG1bzRO8n+PH3H+k+\npjttnmzDk188Sfru9KBDFClVVq6EU06BFSvg44/9a5FyaCAwFXgV+DW0vQq8DwwK54ThDEEW+4rg\nIiUtoVICA48fyLdXf8vMS2bSOqk11029jvoj6nPd1OtYsnHJoU8iEuWWLIGTToKdO+HTT+GYYw59\njEg0cs7tcM4NAmrhn47sCNR0zg1yzm0P55zhJGDFviK4SFDMjK5NujLxgoksv34513W+jtcWvUar\nJ1rR6+VeTFoyiYzMsIoci5RpCxf63q5q1WDOHGjZMuiIRILnnNvunPsmtIWVeGUJJwEr9hXBRUqD\nhtUbcl+3+1g5ZCUvn/MyW3dvpe9rfWn+WHMemfMIv+/4PegQRUrExx9D167QtCl88gnUrx90RCLR\nJ5wELGtF8DoU04rgIqVJXIU4Lmp/EfOumMf8K+ZzWpPTuH3m7TR4tAGXv3M5C39bGHSIIhEzaRKc\ncQZ07gwzZkCtWkFHJBKdwknAin1FcAAzu8bMlpvZTjObZ2adCmh7kpl9amYbzWyHmf1gZv8K99oi\n+elUvxMv/fUlVg1ZxR2n3sGHyz7kuGeO46T/ncS4b8exJ6PIS4mJlFpjxkC/fn5h7cmTISHh0MeI\nSHjCKUNR7CuCm9mF+AU078QPZaYB08wsvzll2/ELa54CtALuBe4zsyvCub7IoSTFJ3HLKbew7Ppl\nvHnBm1SuUJn+b/an0aONuHPmnaxJXxN0iCKHZeRIuOQSuPRSeP11iIsLOiKR6Bb20qnFuSI4fhX0\np51zY5xzi/GPe+4ALsvn2l875153zv0QiuNVYBo+IROJmAoxFTjn6HOYcfEMvhv0HecefS7D5w6n\n8cjGXDjxQmb/Ols1xaRMcQ5uvx2GDIFhw+CZZyA2NuioRKJfoRbjjuSK4GZWETgOeCDHOZyZTQe6\nFPIcHUNt/68o1xY5HK2TWvNEnyd4oPsDjEkbw+gvRnPqi6fSPqU9gzsNpn+7/sRXig86TJF8ZWb6\n5YSefBIefhhuvjnoiETKj0IlYER2RfDaQCx+Ue+c1gFHFXSgma3ELwMQC9zlnHshjOuLHJbqlatz\nbedrueaEa5ixbAajvxjNPyf/k5un38xlx1zGoE6DaF6zedBhihxgzx4/5Dh+PDz7LFyhCRwiJapQ\nCVhoJXAAzOwGIB24JGtRSjM7AngBmB2JIAtwMv5BgBPxT2b+7Jx7vaADhg8fzuuvH9gkNTWV1NTU\nyEUp5UKMxdCzeU96Nu/J8s3LeerLp3juq+d4dN6jnNXyLAZ3GswZLc4gxsIe+RcpFjt2wHnn+acc\nx4+Hc88NOiKR8qewPWA5DQV65V4R3MxuAz7AT6Yvio343rOUXPtTgLUFHeic+zX08rtQWYy7gAIT\nsKFDhzJAq8hKhDU9oikP93yYu7rexbhF43h8/uP0frU3LWq2YNDxg7i046XUqFwj6DClHNq82T/l\nmJYGU6ZAjx5BRyRSPoXzp3ixrgjunNuLH7rsnrXPzCz0/rMinCoW0HM7UqpUqViFyzpexsKrFjLn\nsjl0qteJm6ffTP0R9Rk4eSDfrvs26BClHPntN19gdfFi+OgjJV8iQQonASv2FcGBEcCVZnaxmbUC\nngKqAi8CmNmDZvZSVmMzG2RmfzazFqHtcnzP3MthXl8kosyMPzX8E6+e+yor/rWCYScNY9KSSbR/\nqj1dX+zKxO8nsjdjb9BhShRbtgxOPhl+/x1mz4YTTgg6IpHyLZwErNhXBHfOjQduBO4BvgLaA2c4\n5zaEmtQBGuaK+8FQ2y+Aq4GbnHN3hnN9kZJUt1pd7jjtDn7916+8ft7rZLpMzp9wPk1HNeW+T+5j\n3bbcz6OIHJ5vv/XJV2ysX1S7deugIxIRC7dmkZnF44uxAiw93EUpI83M+gNjX3nlFc0Bk1InbW0a\no+ePZuy3Y9mXuY8L2lzA4BMG07l+Z/yIvEh45s6F3r2hSRN4/31IyT3bVqQQxo4dy0UXXQQwIFR7\nUw7T4RRiLbYVwUXKuw51OvDs2c+y6oZVPNj9QeaumkuX57vQ6dlOjJ4/miUbl6jAqxTZtGl+nle7\ndjBrlpIvkdJEz8OLlCI1q9Rk6J+G8uPgH5mcOpmk+CSGTBtCqyda0WhkIy5951Je+eYVfkv/LehQ\npZQbPx7+8hfo1s0nYtWrBx2RiOQUThkKEYmw2JhY+hzZhz5H9mHbnm188usnzFg2g+nLp/Pi1y8C\n0CapDT2a9aB70+6c1uQ0EuMSgw1aSo2nn4arr4YBA+B//4OKFYOOSERyUwImUsolVEqgd8ve9G7Z\nG4B129bx0fKPmL5sOm8tfotRn48i1mLp3KAzPZr2oEezHnRu0JlKsZUCjlxKmnPw0ENw661+iaGR\nIyFG4xwipZISMJEyJiUhhdR2qaS2S8U5x9LNS5m+bDrTl03n8fmPc88n9xBfMZ5TG59Kj2Y+IWub\n3FYV+KOcc3DTTTB8ONx9t19gW89viJReSsBEyjAzo0XNFrSo2YKBxw8kIzODr9Z+lT1ceeuMWxn6\nwVCS45Pp3rQ73Zt2p0ezHjSu0Tjo0KUY7dsHV10FL7wAjz8OgwcHHZGIHIoSMJEoEhsTy/H1juf4\nescz7ORh7Ny7k89WfsaM5TOYvmw6ry16DYejRc0W2cOVpzc9nZpVagYduoRp1y5ITYXJk+GVV/y8\nLxEp/ZSAiUSxKhWr0L1Zd7o3684D3R9g085NzPplFtOXTWfG8hk8teApDOPYusdmD1ee1PAkqlSs\nEnToUgjp6dC3r6/19fbb0KdP0BGJSGEpARMpR2pWqUm/o/vR7+h+AKzYsiJ7uPKFr1/g4TkPExcb\nx0mNTsruITu27rHExsQGHLnktmGDL7D644/wwQdwyilBRyQiRaEETKQca1S9EZd2vJRLO16Kc45F\n6xdlD1feP/t+bv3oVmpUrkG3pt2y54+1rNlS1fkDtnIl9OoFmzbBxx/DMccEHZGIFJUSMBEB/IT+\ndintaJfSjn+d+C/2ZOxh/ur52cOV179/Pfsy99EwsWH2cGW3pt2ok1An6NDLlSVLoGdPX17i00+h\nZcugIxKRcCgBE5E8VYqtxMmNTubkRidzV9e7SN+dzie/fuJLXoSGLAHaJrfNHq48tfGpVIurFnDk\n0WvhQjjzTEhK8sOO9esHHZGIhEsJmIgUSrW4atnV+QHWblubXRD2jR/eYOTnI6kQU4HO9Ttn95B1\nrt+ZirEqw14cPv7YLy3UujVMmQK1agUdkYgcDiVgIhKWOgl16N+uP/3b9cc5x0+bfsoerhz1+Sju\n/vhuEiolcFrj07Lnj7VNbqv5Y2GYNAkuuMBPtH/rLUhICDoiETlcSsBE5LCZGUfWOpIjax3JoE6D\nyMjMYOFvC7MTsltm3MINH9xASnwK3Zt1p0fTHnRv1p1G1RsFHXqpN2YMXHYZnHOOr/MVFxd0RCJS\nHJSAiUixi42JpVP9TnSq34lbTrmFnXt3MmflnOwlk8Z9Ow6Ho2XNltnDlac2PpXaVWsHHXqpMnIk\nDBkCV1wBTz0FsaoGIhI1lICJSMRVqVglO9EC2LRzEzOXz2T6sul8uOxD/vvlfwGoV60eHVI6+K2O\n/3pkrSPLXR0y5+COO+C++2DYMHjwQa3rKBJtlICJSImrWaUm57Y+l3NbnwvAL3/8wrxV80hbm0ba\nujRe/uZlHprzEACVK1SmbXLbAxKz9intqVG5RpC3EDGZmXDttfDkk/DII36BbRGJPkrARCRwTWo0\noUmNJvyt7d+y923csZFv1n2TnZQt+G0BL3/zMnsy9gDQuHpjn4wlt8/uLWteszkxFhPUbRy2PXvg\nkktg/Hh47jm4/PKgIxKRSFECJiKlUu2qtenWtBvdmnbL3rc3Yy+LNy4mbV0aaWvT+Gb9Nzy78FnW\nbV8HQHzFeNqltDugt6xdcrsyUZtsxw447zyYMQMmTIB+/YKOSEQiSQmYiJQZFWMrZlfrv6j9Rdn7\n121bl52Upa1LY87KOTz/1fPsy9wHQPMjmmf3krVPaU+HlA40qdGk1JTE2LwZ/vxnSEvzNb569Ag6\nIhGJNCVgIlLmpSSk0CuhF72a98ret3vfbr7f8D1p69L8UOa6NB77/DF+3/k7AIlxidnJWFZvWdvk\ntlStWLVEY//tN1/dftUq+OgjOOGEEr28iARECZiIRKW4CnF0rNuRjnU7Zu9zzrEmfc0BvWUzls/g\nv1/+l0yXSYzF0LJmy4N6yxokNohIb9myZX5dx927YfZsX+VeRMqHUpOAmdk1wI1AHSANuNY590U+\nbc8BrgaOAeKA74C7nHMflFC4IlIGmRn1E+tTP7E+vVv2zt6/Y+8Ovlv/XXZPWdq6NKb9PI0tu7cA\n/qnN3L1lrZNaU7lC5bBj+fZbOOMMX9V+zhxo3Piwb09EypBSkYCZ2YXAcOAqYD4wBJhmZkc65zbm\nccipwAfALcAfwGXAu2Z2gnMurYTCFpEoUbVi1ezCsVmcc6zYsuKA3rIpP03hsc8fw+GItVha1W51\nUG9ZnYQ6h+wtmzsXeveGJk3g/fchJSXCNygipU6pSMDwCdfTzrkxAGY2EOiDT6weyd3YOTck167/\nM7O+wF/wvWciIofFzGhcozGNazTm7KPOzt6/bc82Fq1flJ2Upa1LY9KSSWzbsw2ApKpJ2UlZVm9Z\nq9qtqBRbCYBp0/wTjscdB+++C9WrB3J7IhKwwBMwM6sIHAc8kLXPOefMbDrQpZDnMKAasCkiQYqI\nhCRUSuDEBidyYoMTs/dlukyWb15+QG/Zmz+8yfC5wwGoGFOR1kmtqbS5AwumtqVjvxb8577mxFZp\nBmhlbZHyKPAEDKgNxALrcu1fBxxVyHPcBMQD44sxLhGRQomxGJrXbE7zms3pd/T+Al5bdm3h2/Xf\n8u4XafxvShobY9Oo2O0NFth2Or/o2yTHJ9P8CH9s8yP81uyIZjSv2ZyU+JRSUypDRIpXaUjADouZ\n9QduB87OZ76YiEggbE91Jgw/mdGjT6ZtW3j3aejc2bFhxwaWblrK0s1Ls78u27yM6cums3bb2uzj\n4yvGZydj2YlZKFlrXL0xFWMrBnh3InI4SkMCthHIAHJPQ00B1h7cfD8z+xvwDHCec25mYS42fPhw\nXn/99QP2paamkpqaWuiARUQK4hy88QZcfz388Ydf0/H666FCBQAjOT6Z5PhkujQ8eJbF9j3bWbZ5\nWXZylvX6nSXv8Msfv2QXl421WBpVb3RAUpbza1mo/i9SngWegDnn9prZAqA7MAmy53R1Bx7L7zgz\nSwWeAy50zr1f2OsNHTqUAQMGHF7QIiL5WL4cBg+G996Ds8+Gxx+HRo0Kf3x8pfjsav+57cvcx8ot\nK7N7zLJ6z+avmc+4ReNI35Oe3bZ21doHDW02r+l70eom1NXQphS7opSTCrXviq+A0AZYAdzvnHsp\nV5vzgXuAJsCPwL+dc1OLcl0zewG4JNfl33fO9SZAgSdgISOAF0OJWFYZiqrAiwBm9iBQzzl3Seh9\n/9D3rgO+MLOs3rOdzrmtJRu6iAjs3QvDh8M990Dt2vD229C3b/Feo0JMBZoe0ZSmRzQ96HvOOX7f\n+ftBQ5tLNy9l1i+zWJO+JrttlQpV8h3abFKjSfYTmyKFVdRyUmbWBJgMPAn0B3oAz5nZGufch6E2\nfwJeBYYBU4ABwNtm1tE5930RrzsV+AeQ9ZfH7mK7+TCVigTMOTfezGrjs9wU4GvgDOfchlCTOkDD\nHIdciZ+4/0Roy/ISvnSFiEiJ+fRTGDgQFi/2Q4133+0LrJYkM6N21drUrlqbzg06H/T9HXt3sHzz\n8oOGNqf8NIXlm5ezN3Mv4B8oaJjY0PeW1Wh20NBm9cqqmyF5KlI5KXwx9WXOuZtD75eY2cmh83wY\n2ncdMNU5NyL0/g4z6wkMBgYV8bq7c+QUpUKpSMAAnHNP4jPhvL53aa73p5dIUCIiBfj9dxg2DJ5/\nHjp3hi+/hGOOCTqqvFWtWJU2yW1ok9zmoO9lZGawauuqg4Y2F65dyITvJ2SvCAB+VYD8hjbrVatH\njMWU5G1JKRBmOakTgem59k0DHs3xvgu+dyt3m75hXLerma0DNgMfAbc55wItXVVqEjARkbLCOXj5\nZRg61A89PvkkXHUVxMYGHVl4YmNis4vOdmva7YDvOefYtHPTAQ8GZA1tfrriU1ZtXZXdtnKFyjSt\n0TS796x+Yn0aJDagfjW//FP9avWpUrFKSd+eRF445aTq5NM+0czinHO7C2hTp4jXnQq8ASwHmgMP\nAu+ZWRfnnCvgviJKCZiISBEsXgxXXw2zZkFqKowYAXXqHPKwMsvMqFW1FrWq1jpgqaYsu/btOmBo\nM6sX7cNlH7I6fTVbdx84LfeIykf4pCyUkGUlZzkTtVpVaukhASk2zrmcNUK/M7NvgaVAV6BQFRQi\nQQmYiEgh7NwJDz4IDz3kn2qcNg169Qo6quBVrlCZo5OO5uiko/P8/rY921i9dTWr01ezauuq7Ner\n01fzzbpveO+n91i3fR2ZLjP7mLjYOOpVq1dgola3Wl09LBAh48aNY9y4cQfsW7VqVT6tgfDKSa3N\np/3WUO9XQW2yzhlWGSvn3HIz2wi0QAmYiEjp9eGHvtdrxQr497/hllugikbSCiWhUgJH1T6Ko2rn\nv7DJvsx9rN22dn9ytjWUrIUStQVrFrA6fTU79u444Ljk+OT9iVm1HMla6GuDxAYkxiWqN62I8qqN\nOXbsWC666KI824dZTmoucFaufb1C+3O2yX2OnlltDqOMVQOgFvBbfm1KghIwEZF8rF0LN9wA48ZB\n164weTK0ahV0VNGnQkwFGiQ2oEFig3zbOOf4Y9cf2Qla7kRt3up5rP5hNRt2HPigW3zF+AMSs7wS\ntToJdYiNKaMT+EqPIpWTAp4CrjGzh4H/4ZOm84CctblGAbPM7AZ8GYpU/KT7K4tw3XjgTvwcsLX4\nXq+H8TXFphXb3YdBCZiISC6ZmfD0076nq2JFeOkl+PvfQR0pwTEzjqhyBEdUOYK2yW3zbbd7327W\npK85KFFbnb6a5ZuX8+mKT1mTvoY9GXuyj4mxGOom1C1wXlr9avWJrxRfErdaJhW1nJRz7hcz64N/\n6vE6YBVwuXNueo42c0N1P+8PbT8BfbNqgBXyuhlAe+BioAawBp943eGc21vMP4YiUQImIpLD11/7\nml6ffw5XXOHnfNWqFXRUUlhxFeLyLVabxTnHxh0bD56XFvo685eZrE5fzR+7/jjguOpx1fOcl5YS\nn5K9vFRyfHK5HfYsSjmp0L5P8D1aBZ3zDXzvVbjX3QWcWdDxQVECJiICbNsGd94Jo0b5YcbZs+Hk\nk4OOSiLBzEiKTyIpPolj6uRfuG37nu3ZvWm5HyD4fsP3fLD0A9ZuW0uGyzjguEqxlQ5IyJLjk0mu\nmnzwvvhkkuKTqFyhcqRvWUohJWAiUu69/TZce60vrHr//TBkCFTSA3blXnyleFrWaknLWi3zbZPp\nMtm0cxPrt6/Pd1u6aSlzV85l/fb1BxS1zZIYl1ioZC05PpmaVWpqvlqUUAImIuXWihU+8Zo0CXr3\nhtGjoWn+I1ciB4mxmOwloFontT5k+937drNhx4YCE7YFvy3Ifr0748AlC7OuV9iELaFSQrkcDi0L\nlICJSLmzd68farzzTqhRAyZMgHPP1SR7iby4CnGHfOIzi3OO9D3pBSZr67at49t137J++3o27tiI\n48DC7pUrVC7ScKhqq5UcJWAiUq7Mmwf//CcsWgSDB8O990JiYtBRiRzMzEiMSyQxLpEWNVscsn1G\nZga/7/y9wIRtycYlzN4+m/Xb15O+J/2gc9SoXCPPZO3XNb9G4hbLNSVgIlIubN4Mt97qy0sceyzM\nnw/HFfj8lUjZEhsTm504FcbOvTsPORz6+ebPWb99PWuX5ltYXsKkBExEoppzvpDqkCF+OaFRo2DQ\noLK7cLZIcalSsQqNqjeiUfVGh2w7duxYLno670r4Ep6YoAMQEYmUn37y6zUOGACnnQY//OAn3Sv5\nEpGgKQETkaizezfccw+0awc//wxTpsD48VC/ftCRiYh4GoIUkagyc6avZL9sGdx0E9x2G1StGnRU\nIiIHUg+YiESF9evh4ouhWzdITvZLCj3wgJIvESmd1AMmImVaZiY8/zwMG+breD3/PPzjHxCjPy9F\npBTTrygRKbMWLYJTT4WrroKzz4bFi+Gyy5R8iUjpp19TIlLmbN/ue7w6dvTrN86cCS++CElJQUcm\nIlI4GoIUkTJlyhS45hpYu9YvJXTTTRAXF3RUIiJFowRMRMqEVavg+uvhzTd9ba/p06HFoVdnEREp\nlUrNEKSZXWNmy81sp5nNM7NOBbStY2ZjzWyJmWWY2YiSjFVESs6+fb56/dFHw2efwWuvwfvvK/kS\nkbKtVCRgZnYhMBy4E+gIpAHTzKx2PofEAeuBe4GvSyRIESlxX3wBnTv7ZYQuvthXsr/wQv+0o4hI\nWVYqEjBgCPC0c26Mc24xMBDYAVyWV2Pn3K/OuSHOuVeArSUYp4iUgC1bYPBgn3xlZsK8efDEE1Cj\nRtCRiYgUj8ATMDOrCBwHzMja55xzwHSgS1BxiUjJc84vGXT00f6pxuHDfS/YCScEHZmISPEKapu4\nXAAAD5dJREFUPAEDagOxwLpc+9cBdUo+HBEJwrJl0Lu3H2I88UQ/3DhkCFTQo0IiEoX0q01EApOR\nAXPmwIQJ8NxzfgmhSZPgL38JOjIRkcgqDQnYRiADSMm1PwVYW9wXGz58OK+//voB+1JTU0lNTS3u\nS4lIHjIyYPZsn3S9+aav51W/PgwdCrfcAvHxQUcoIhJ5gSdgzrm9ZrYA6A5MAjAzC71/rLivN3To\nUAYMGFDcpxWRAuzbBx9/DBMn+qRr/Xpo2BBSU+H88/1key0fJCLlSeAJWMgI4MVQIjYf/1RkVeBF\nADN7EKjnnLsk6wAz6wAYkAAkhd7vcc79UMKxi0ge9u6FWbN8T9dbb8HGjdC4sS8ncd55fmK9ykmI\nSHlVKhIw59z4UM2ve/BDj18DZzjnNoSa1AEa5jrsK8CFXh8L9Ad+BZpFPmIRycuePfDRR76n6623\nYNMmaNbML5B9/vlw3HFKukREoJQkYADOuSeBJ/P53qV57NOAhUgpsGePXxZowgR45x3YvNlXqf/n\nP31PV8eOSrpERHIrNQmYiJQdu3fDBx/4nq533vGFU486yi+Sfd550L69ki4RkYIoARORQtm5E6ZN\n80nXpEmQnu4Lpl5/vR9ebNNGSZeISGEpARORfO3Y4Re+njABJk+GbdugbVu48Ubf09W6ddARioiU\nTUrAROQA27fDe+/5nq4pU/z7Dh1g2DCfdLVqFXSEIiJlnxIwEWHbNp9sTZjgk6+dO/3k+f/7Pzj3\nXDjyyKAjFBGJLkrARMqprVv9sOLEiTB1KuzaBccfD3fd5ZOu5s2DjlBEJHopARMpR7ZsgXff9T1d\n06b5pxk7d4Z77/XDi02aBB2hiEj5oARMJMr98YcvFTFxoi8dsWcPdOkCDz7oe7oaNQo6QhGR8kcJ\nmEgU2rTJJ10TJvgiqfv2wUknwSOP+KSrQYOgIxQRKd+UgIlEiY0b4e23fU/XjBmQkQGnnAIjRkC/\nflCvXtARiohIFiVgImXY+vV+zcWJE2HmTHAOTjsNRo3ySVedOkFHKCIieVECJlLGrF27P+maNctX\nnz/9dHjiCfjrXyElJegIRUTkUJSAiZQBv/0Gb7zhk65PPoGYGOjeHZ56yiddSUlBRygiIkWhBEyk\nlNm8GX76CX7+2X+dPh3mzIHYWOjZE557Dvr2hVq1go5URETCpQRMpIQ5559SzEqyshKtrNebNu1v\nm5Lii6O+8AKcfTYccURwcYuISPFRAiYSAc75pxJzJlY5E60//tjftk4daNkS2rTxPVstW0KLFr4S\nfWJicPcgIiKRowRMJEzO+acQc/dgZb3eunV/23r1fGLVoYOvw5UzyUpICO4eREQkGErARArgnH/q\nMHeSlbWlp+9v26CBT6qOOw4uvNC/btkSmjWD+Pjg7kFEREofJWBS7jkHa9bkPR/r559h+3bfzgwa\nNvSJ1QknQP/+/nVWT1aVKsHeh4iIlB1KwKRcyMyE1avzT7J27vTtzKBxY59UdekCF1+8P8lq1gwq\nVw72PkREJDooAZOokZkJq1blPR9r6VLYtcu3i4nxSVbLln6pnksv3T8nq2lTiIsL9j5ERCT6KQGT\nUm/3bj+hPT3df83aVq48MNFatsy3BV8zq0kTn1idfjpceeX+OVlNmkClSkHekYiIlHdKwCQi9u3b\nnzDlTJzCeb1nT97XqFDB91i1aOELlGYNFbZs6Xu4KlYs2XsWEREpLCVgki0z0084P5xkKet11pyq\n/CQkQLVqvs5VYuL+182a5b0/r9dJST4JExERKWtKzceXmV0D3AjUAdKAa51zXxTQviswHGgDrADu\nd869VAKhHjbn/JaR4ZOezMz8Xxf0vbxeZ2TAu++Oo1On1CInTtu2+bjyU7ly3olQvXrQqlXhE6eE\nBD9EeLjGjRtHamrq4Z+olNN9RhfdZ/QpT/dakEh8jpvZ+cA9QBPgR+DfzrmpRb2umd0DXAHUAOYA\nVzvnfg73XotDqUjAzOxC/D/CVcB8YAgwzcyOdM5tzKN9E2Ay8CTQH+gBPGdma5xzHxZ0rf/8B8aO\nLd7kp6hJUkFJTvEYB/hfBhUq5J0I1azp50IdKlnK+bq0DemVl196us/oovuMPuXpXvMTic9xM/sT\n8CowDJgCDADeNrOOzrnvC3tdMxsGDAYuBn4B7gu1Odo5l88kl8grFQkY/gf2tHNuDICZDQT6AJcB\nj+TR/mpgmXPu5tD7JWZ2cug8BSZgMTG+JycmxvfCxMQU7+tInbco1xgyBMaN84lTXJwvrSAiIhJB\nkfgcvw6Y6pwbEXp/h5n1xCdTg4pw3euBe51zk0NtLgbWAX8Fxh/ujYcr8ATMzCoCxwEPZO1zzjkz\nmw50yeewE4HpufZNAx491PVuuAEGDAgz2DIiPh6Sk4OOQkREyoMIfo53wfdu5W7Tt7DXNbOm+KHJ\nGTnabDWzz0NtAkvAYoK6cA61gVh8NprTOvwPLS918mmfaGaq4iQiIlJyIvU5nl+brHMW5rp1AFfE\n2EpE4D1gJSgBYM6cOUHHEXGrVq1i7NixQYcRcbrP6KL7jC7l5T6hfNxrjs/OhCDjiCrOuUA3oCKw\nFzg71/4XgbfyOeZjYESuff8ANhdwndH4LFibNm3atGnTFt42uqQ+x4FfgetytbkL+Kqw1wWaAplA\n+1xtZgGPBpn/BN4D5pzba2YLgO7AJAAzs9D7x/I5bC5wVq59vUL78zMy9PUbYFvYAYuIiJQ/CUB7\n9n+WZovg5/jcPM7RM6vNIa77eKjNcjNbG9r3TahNItAZeKIQ9x0xgSdgISOAF0M/yKzHSKvis1jM\n7EGgnnPuklD7p4BrzOxh4H/4H+x5QO/8LuB8vY/BkboBERGRciwSn+OjgFlmdgO+DEUqftL9lYW4\n7gs52owEbjOzn/FlKO4FVgHvFMeNh6tUJGDOufFmVhtfbC0F+Bo4wzm3IdSkDtAwR/tfzKwP/mmJ\n6/A/yMudc7mfqBAREZEIi8TnuHNurpn1B+4PbT8BfV2oBlghr4tz7hEzqwo8jS/EOhs4ywVYAwzA\nXOSrgoqIiIhIDqWhDIWIiIhIuaIETERERKSERX0CZmanmNkkM1ttZplmdnbQMRU3M7vFzOab2VYz\nW2dmb5nZkUHHFQlmNtDM0sxsS2j7zMzODDquSDKzf4f+7444dOuyxczuDN1bzu37Qx9Z9phZPTN7\n2cw2mtmO0P/jY4OOqziZ2fI8/j0zzezxoGMrTmYWY2b3mtmy0L/lz2Z2W9BxRYKZJZjZSDP7JXSv\nn5rZ8UHHFQ2iPgED4vGT8gbha5hEo1Pwj9x2xi9oWhH4wMyqBBpVZKzEL8x6LP5pmI+Ad8zs6ECj\nihAz64RfZDYt6FgiaBF+8myd0HZysOEUPzOrAcwBdgNnAEcDQ4HNQcYVAcez/9+xDr5kgCPA5V4i\n5N/AP/GfK62Am4GbzSwan7R/Hv+E4gCgLX6dxulmVjfQqKJAuZqEb2aZwF+dc5OCjiWSQk+ErAdO\ndc59GnQ8kWZmvwM3OudeOGTjMsTMEoAF+EVrb8cXH7wh2KiKl5ndiX+qKap6gnIzs4eALs6504KO\npSSZ2Uigt3MuqnrkzexdYK1z7soc+yYCO5xzFwcXWfEys8pAOvAX59z7OfZ/CbznnLsjsOCiQHno\nASuPauD/6twUdCCRFBoG+Bu+5ktBRXjLqieAd51zHwUdSIS1DE0RWGpmr5hZw0MfUub8BfjSzMaH\npgksNLMrgg4qkkILJQ/A96BEm8+A7mbWEsDMOgAnAe8FGlXxq4Bfa3F3rv07icKe6pJWKuqASfEJ\nVQEeCXyas1ZKNDGztviEK+uvs3Occ4uDjap4hRLLY/BDOtFsHn75kSVAXfwyI5+YWVvn3PYA4ypu\nzfA9mcPx9YxOAB4zs93OuZcDjSxyzgGqAy8FHUgEPAQkAovNLAPfmfF/zrnXgg2reDnntpnZXOB2\nM1uMX8C6P9AFX5NLDoMSsOjzJNAa/9dYtFoMdMD/cj8PGGNmp0ZLEmZmDfBJdA/n3N6g44kk59y0\nHG8Xmdl8/PpvF3BgJeuyLgaY75y7PfQ+LfSHxEAgWhOwy4Cpzrm1QQcSARfiE5G/Ad/j/1gaZWZr\nojChvghfqX41sA9YCLyKn4Mrh0EJWBQxs9H4ZRxOcc79FnQ8keKc2wcsC739ysxOAK7H9zBEg+OA\nJGBhqEcT/DDAqaFJvnEuSidvOue2mNmPQIugYylmvwE/5Nr3A9AvgFgizswa4R8I+mvQsUTII8CD\nzrkJofffmVkT4BaiLKF2zi0HTg891JXonFtnZq+x/3ewhElzwKJEKPnqC5zunFsRdDwlLAaICzqI\nYjQdaIf/q7pDaPsSeAXoEK3JF2Q/eNACn7BEkznAUbn2HYXv7YtGl+GHq6JtTlSWqkBGrn2ZRPFn\nqnNuZyj5OgL/JO/bQcdU1kV9D5iZxeN/oWf1JDQLTZjc5JxbGVxkxcfMnsQvUno2sN3MUkLf2uKc\n2xVcZMXPzB4ApgIrgGr4Sb6nAb2CjKs4heY+HTB/z8y2A78753L3opRpZvb/gHfxiUh94G5gLzAu\nyLgi4FFgjpndgi/J0Bm4ggMXFY4KoV7bfwAvOucyAw4nUt7FL+68CvgOXxZnCPBcoFFFgJn1wn9+\nLgFa4nv/vie0yLaEL+oTMPwk5pn4pwIdfhIs+ImhlwUVVDEbiL+3Wbn2XwqMKfFoIisZ/29XF9gC\nfAP0KgdPCkZrr1cD/HySWsAG4FPgROfc74FGVcycc1+a2Tn4ydu3A8uB66Nt0nZID/yiy9E0hy+3\nwcC9+CeVk4E1wH9D+6JNdeBB/B9Im4CJwG3Oudw9gFJE5aoOmIiIiEhpELXj1SIiIiKllRIwERER\nkRKmBExERESkhCkBExERESlhSsBERERESpgSMBEREZESpgRMREREpIQpARMREREpYUrARCQwZjbT\nzEYEHYeISElTAiYiIiJSwpSAiYiIiJQwJWAiUmqYWR8z+8PMUoOORUQkkioEHYCICICZ9QeeBFKd\nc1ODjkdEJJLUAyYigTOzQcBo4M9KvkSkPFAPmIgE7XwgCTjJObcg6GBEREqCesBEJGgLgQ3A5UEH\nIiJSUpSAiUjQlgKnA33N7PGggxERKQkaghSRwDnnfjaz04GZZrbPOTck6JhERCJJCZiIBMllv3Du\nRzPrzv4k7KYA4xIRiShzzh26lYiIiIgUG80BExERESlhSsBERERESpgSMBEREZESpgRMREREpIQp\nARMREREpYUrAREREREqYEjARERGREqYETERERKSEKQETERERKWFKwERERERKmBIwERERkRKmBExE\nRESkhP1/KyZ5P/zzCpQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x8f52978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 10000\n",
    "theta = 10\n",
    "moments = 9  # number of moments to analyze\n",
    "\n",
    "# hats[k-1] = hats of k-th moment\n",
    "hats_exponential = [[] for _ in range(moments)]\n",
    "hats_uniform = [[] for _ in range(moments)]\n",
    "\n",
    "for _ in range(100):\n",
    "    data = np.random.exponential(1. / theta, size=(n, 1))\n",
    "    mean = np.mean(data)\n",
    "\n",
    "    for k in range(1, moments + 1):\n",
    "        theta_hat = np.power(math.factorial(k) / np.mean(data**k), 1 / k)\n",
    "        hats_exponential[k - 1].append(theta_hat)\n",
    "\n",
    "for _ in range(100):\n",
    "    data = np.random.uniform(0, theta, size=(n, 1))\n",
    "    mean = np.mean(data)\n",
    "\n",
    "    for k in range(1, moments + 1):\n",
    "        theta_hat = np.power((k + 1) * np.mean(data**k), 1 / k)\n",
    "        hats_uniform[k - 1].append(theta_hat)\n",
    "\n",
    "devs_exponential = np.array([np.mean([(theta - theta_hat)**2 for theta_hat in hat]) for hat in hats_exponential])\n",
    "devs_uniform = np.array([np.mean([(theta - theta_hat)**2 for theta_hat in hat]) for hat in hats_uniform])\n",
    "\n",
    "plt.title('SKO')\n",
    "plt.xlabel('k')\n",
    "\n",
    "ax1 = plt.axes()\n",
    "plt.ylabel('dev exponential')\n",
    "g1, = ax1.plot(range(1, moments + 1), devs_exponential, 'b', label='Exponential')\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "plt.ylabel('dev uniform')\n",
    "g2, = ax2.plot(range(1, moments + 1), devs_uniform, 'g', label='Uniform')\n",
    "\n",
    "plt.legend(handles=[g1, g2])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}