mixture_distributions_in_sympy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mixture Distributions in Sympy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying to figure out how to specify mixture distributions (e.g. mixture of Gaussians) in sympy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import numpy as np,pandas as pd\n",
    "\n",
    "from sympy.stats import E,Normal,sample_iter\n",
    "from sympy.stats.crv_types import _value_check,rv,SingleContinuousDistribution\n",
    "from sympy import exp,sqrt,pi\n",
    "\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define 2 Gaussian distributions, and a 'mixture' as their sum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gaussian 1\n",
    "G1 = Normal('G1',5,1)\n",
    "\n",
    "# Gaussian 2\n",
    "G2 = Normal('G2',20,4)\n",
    "\n",
    "# Sympy G1G2 mixture\n",
    "G1G2 = G1 + G2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw samples from these three distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "G1_samps = np.array(list(sample_iter(G1,numsamples=1000)))\n",
    "G2_samps = np.array(list(sample_iter(G2,numsamples=1000)))\n",
    "G1G2_samps = np.array(list(sample_iter(G1G2,numsamples=1000)))\n",
    "\n",
    "df_samps = pd.DataFrame([G1_samps,G2_samps,\n",
    "                         G1G2_samps],index=['G1', 'G2', 'G1G2']).T.astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concatenate the G1 and G2 sample vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mixsamps = pd.DataFrame(np.concatenate([G1_samps,G2_samps])\n",
    "                           ,columns=['G1G2_mix']).astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot these four distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9oAAAEnCAYAAABWhDupAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X28nHV94P1PzgnZhCRHQhCwrImr0i+KQAPCzb1b69a6\nbVdu67rl3nVd77Vii1ZbERBEbatbS4uwRvrAAqU+3y1qBbtV3LK1WutKKYSEEE7IVyXyZAhCCDk5\nSeDkPOwf10yYDOckM3PmeT7v1yuvMzPXb67f95rJb2a+1/V7WDAzM4MkSZIkSWqOoU4HIEmSJElS\nPzHRliRJkiSpiUy0JUmSJElqIhNtSZIkSZKayERbkiRJkqQmMtGWJEmSJKmJTLQlSZIkSWoiE21J\nkiRJkprIRFuSJEmSpCZa2OkA1H0i4jXAu4F/CawE9gKbgM8Df5aZ01XlLwSuAG7KzDe3OVxJs6in\nHUfEecBvAi8FdgD/C/hQZj7e7rilQdCK79mIeAlwCfBvgJ8AZoD7gZuAP8zMnVXlXwV8BPgpYBJY\nB3wwMzc25yil7jXobTAiVgM/BN6ZmX/a6voGlVe0dZCI+F3gVuAh4A3Ai4HXAN8E/gj4q4qyKyLi\nr4GLKD6gJHWBOtvxRcCfAp8FTgPeAfwi8OX2Ri0NhlZ8z0bEOcA9wAuBdwI/CbwS+GPgV4E7IuK4\nivJnA98AHgR+Gvi3wFLgGxFxbLOOVepGtkGgOPbjKb771SJe0dYBpQ+J3wLekZk3VGx6BFgfEXcD\nfxYR/yozvwu8GTgSWAPc2faAJT1HA+34fcBnMvPqUrmtpR8h10fEKZm5qa0HIPWxVnzPRsTxwJ8D\nX8nMt1Rt3hwRtwJ3AP8JKLfz9wIPZOZ5Ffv5VeB7wH+kSA6kvmMbLGTmDPDjVtYhE20d7CJgc9UH\nzwGZ+ZWI+B8V3Wm+Bvz3zJyJiLYFKemQ6m3HLwemqor9qPR3WYtilAZVK75nf40iEbh4jn0+GBEv\nqOoK+zaKq2eVbPcaBH3XBiNiGng/cCzwVmAxRXf1dwK/DbwdWERxpf4dmTlZ2XUc+BTFiYDdmfnq\niv1+hOJk/E9l5g/qiUkFE20BEBHDFONUDnkGrfJDIjMfbHVckmrXYDt+apYibwD2APc2NUBpgLXw\ne/bVwN2Z+Vgt+yzd3wfsqyr2Booxpf9YQ51Sz+nzNng+8GngbOB1FF3gTwS+Q9E1/RcpjvvvKcah\nV8YyGRH/BbgzIt6WmZ+OiBMpkveLTLIb5xhtlR0D/DOKMRuSetO823FEvJ7i7Pzlmbm7WYFJatn3\n7Anz3WdEvIjiR/itmfn3TYhJ6kb93AYfy8zfz8ytmfknwBhwTGZ+IDPvz8xrSo+tme3JmXkvxcRs\nV0bE84FrgX/IzGsbiEUlXtFWWflM23Dlg6XGtpXiDNuC0sP/kJnntDE2SbWZVzuOiP+X4kz35zPz\nYy2OVRo0rfqena7eZ2m/64DKvq4zmTkyS7mXU0wM9QjFeFSpX/VzG7yr6v5TQPXs5U8BzzvEPq4C\nXg/8A3Ac8IoGY1GJV7RVtoOiq+iLZ3n8NIqlB06j+CBY0t7QJNWo4XYcEb8J3Ahcm5m/0vJIpcHT\nqu/ZB2fZJ8C/K+3vNOB3ee54UCLipym6lt4P/Os5hpJI/aKf2+Ceqvszczy2gDmUurdfT3Fy4H9m\n5rYGY1GJibaAA43rm8AvRcRQ5eOlbihbM3MrYFdSqUs12o4j4p0UM6G+PzMvbGvQ0oBo4ffs3wIn\nR8RLq+p7pGKfz5ldOCJeCfzPUkz/JjN31Vmv1FNsg4cWESPA5cBXgf8QEa8+zFN0GCbaqnQlxfp/\nvzPbxohYBLykrRFJqldd7TgiXgP8CXBhZn68LRFKg6sV37OfBp4AromIuYYEvryqnudTzKb8N8B/\nyMz9ddYp9Srb4Nz+CBgHzqWYifwzEfGcq/CqnWO0dUBmfjciLgY+HhEvAW6gmPr/ecD/TbHm3yrg\nVwAiYgXFcgELKMamLI6I40q725WZT7f3CCTV244pkuzvAl+saL9l45lZ3fVMUoNa8T2bmU+V5lf4\nH8B3IuJjwN3AEcApFJMb/hug8kTaR0v7vQw4tmrZoonM3NnsY5e6gW1wdqWJUN8CvCoz90fEJcD/\nA/wh8KvtjKWfeEVbB8nMP6RYGmAI+HPg+xRLAbwd+BLwksy8qVT8ZmAbxbp/J1AsS7Ct9O8/tDVw\nSQfU2o4jYhXFWKyf5tm2W/lv1jVBJTWuFd+zmfkd4GSKZYGuAEaBdRQ/5rcCazLz0oowfp4isfge\nz233NyH1sT5sgzOlf/U+NgPMlE4mXE8xR8s/lo5nDPgN4G0R8bo641HJgpmZ6vdAkiRJkiQ1qu6u\n4xGxGriG4kzQbuCLmXnZLOU+DPw2MFF6aAHFmZPVmfl4wxFLOqxa22mp7FKKM5lvBk7KzO9VbDsa\n+ATFmdeFFEs+XJCZj7T2CCSBbVnqF3W25XdSdGH+CeAHwEcy868rtl8OvAk4Cvgn4N2Z+cPWHoGk\nejUyRvsm4E6KBn4c8PWI2J6ZV89S9nOZed58ApTUkJraaUS8APgWRVen2bq3fIbic6I8icfnKSbI\n+PnWhC2pim1Z6g+1tuV/D/w+8LpS+bcCX4qIkzLzgdJSjG8C/i1Fd+Y/AL5CsTSV+kREXEsxZvpQ\nZl2XW92jrjHapWnoT6VYAmY8M+8H1gLntyI4SfWrs50+H3gf8BFmX1vxYeB9mbmzNDHHtRTjeSW1\nmG1Z6g91tuUlwAcy8/bMnMrMT1FcAT+7tP18YG1mfq80WeUHgZdHxFmtPxK10W/z7Brcc/3z5EqX\nq/eK9unAA6UB8mXrgYiIZZk5XlX+tIj4LvAK4CHgosz828bDlVSDmttpZt4D3FPq0vYcmfnuqodW\nAY82O2BJs7ItS/2hnrb855VPjIijgOXAIxGxmKJXyoaK8uMR8X3gTOCOFh6D2igzn6BYMkw9rN5Z\nx1cC1dPNP1mxrdIjFONK3kLRReaTwNci4sR6g5RUl3raac0i4kXA71LMoCmp9WzLUn+YT1u+AfjH\nzPzfwAqKHiuz7euY+QYpqbmasY52uYvaQWPCMvOTFMl12dUR8SaKxPvDTahXUu1mbae1ioiTgFuB\nT2fmZ5oVlKS62Zal/nDIthwRC4HPAi8DfraGfbmMkNRl6k20H+e5Z8yOpmjctXRveIBiBsWazMzM\nzCxYMNtQM0kVqhvJfNvpQUrjvm4BrsrMKxsJ0LYs1aSr27LtWKrJbI2krrZc6iL+18Bi4FWleRWg\nuHI9Pce+al7Rx7Ys1WTejaTeRHsdsDoijs7McpeXs4DNmbm3smBEfAi4LTO/VfHwy4Av1FrZggUL\nGBvbx9TUdJ1h1mZ4eIiRkSUtraNd9VhH99XTzjqq1NxOqzznbHhpqMfXKOZX+HyjcfZDW/b/5mDW\n0a56eqEtt7odQ//8v+mXOtpVT7/VMYt62/IXgKeBczJzf/nBzHwmIu4FzgC+AwfGcL+UYpmvmvid\n3F319Esd7aqnw225LnUl2pl5d0TcAVwRERcDJwAXAlcBRMQW4LzMvI1izMk1EfHvgAeB3wBeQtEN\npmZTU9NMTrbuP0S76mhXPdbRffW061jKamin9wFvL7XTsgXMfubuGuBP55Nkl/XLa+3/zcGso531\nlHVjW+6n19o6uq+efqmjWj2/nyPiPwMnA6dUJtkVrgUui4i/oVje62PAXZm5vp6Y+uW19v9/99XR\nrno60Zbr1cgY7XMpJmbYDuwCrs3M60rbTgSWlW5fRnFW/e8ourSMAq/JzG3zilhSLQ7VTn+SUjst\n9Tz5rdLjM8DGiJgBfg/4HPBa4FWlHwYzPDsO7OdLE7NIai3bstQfDvf7eWnp9tuA1cCTEQHPttXP\nZ+Y7MvP6iDge+HuK9v8t4JfbdRCSald3ol1KlM+ZY9twxe0J4OLSP0ltVEc7vRy4/BC7qndlAklN\nZFuW+kMdbfm1NezrvwL/tXnRSWqFZsw6LkmSJGnATUxMcM89owfG0J500sksWrSo02FJHWGiLUmS\nJGneRkc3cenam1m+chW7dzzElRdNs2bNGZ0OS+oIE21JkiRJTbF85SqOOv7ETochdZxjtiRJkiRJ\naiITbUmSJEmSmshEW5IkSZKkJjLRliRJkiSpiUy0JUmSJElqIhNtSZIkSZKayERbkiRJkqQmMtGW\nJEmSJKmJTLQlSZIkSWqihZ0OQJLUPBMTE4yObgLg5JNPYdGiRR2OSJIkafCYaEtSHxkd3cSla28G\n4MqLYM2aMzockSRJ0uAx0ZakPrN85apOhyBJkjTQHKMtSZIkSVITmWgD69ffxYYNd3U6DEmSJElS\nHzDRliRJkiSpiUy0JUmSJElqIhNtSZIkSZKayERbkiRJkqQmMtGWJEmSJKmJXEdbkvrQ9NQkmVsO\n3D/55FNYtGhRByOSJEkaHCbaktSH9jz1KJ+8ZRvLbx9n946HuPIiWLPmjE6HJUmSNBBMtCWpTy1f\nuYqjjj+x02FIkiQNHMdoS5IkSZLURCbakiRJkiQ1kYm2JEmSJElN5BhtqU9FxGrgGuBsYDfwxcy8\nbI6yS4HrgTcDJ2Xm9yq2HVXa9mpgCvg68BuZ+Uxrj0CS7VitMDExwejoJsAVCSSpVbyiLfWvm4CH\ngRcBrwXeGBHvrS4UES8A7gL2AzOz7OeTwBLgZcAZpb8fa03IaqeJiQk2bLjrwL+JiYlOh6Tnsh2r\n6UZHN3Hp2pu5dO3NBxJuSVJzmWhLfSgiXgmcCrw/M8cz835gLXD+LMWfD7wP+AiwoGo/xwJvAD6Q\nmTszczvwUeBtETHcwkNQG9x7b/Fj+6OfXecP7i5kO1YrLV+5iuUrV3U6DEnqWybaUn86HXggM8cq\nHlsPREQsqyyYmfdk5tfm2M9PAZOZOVq1n+XASc0MWJ1RXgLMH9xdyXYsSVKPMtGW+tNKYGfVY09W\nbKtnP7vm2M8xDcQlqXa2Y0mSepSToUmDo9yddLbxm42oaz/Dw607r1fed6/X0Yx65nre8PAQCxce\nvO+hoQVzlpkv35PG6qhB37bjyv33+v+bbq+j8jmHa/fdfizdWIckgYk2ExMTbNlyHy960Ys7HYrU\nTI/z3CtVR1P8qH6izv0cFRELMrP8g3xlxbaajYwsqad4Q/qljvnUM9fzRkaWsGLF0oMeW7Zs8WHL\nzJfvybwMZDtuVz2DXEflc2pt9916LN1YhyRBA4l2PUuNVDznBOA+4L9l5u82EmgrrF9/F4888kN+\n/5obuezX38RZZ53d6ZCkZlkHrI6IozOz3EX0LGBzZu49xPOqr25toLiCdhpwd8V+dgJZT0BjY/uY\nmpqu5yk1Gx4eYmRkSc/X0Yx6xsb2zfn4zp17DqpjfPzpOcvMl+9JY3VUGah2DP3z/6bb66j8nDhc\nu+/2Y+nGOiQJGruifRNwJ/Am4Djg6xGxPTOvPsRz/giYbKCutli8dEWnQ5CaKjPvjog7gCsi4mLg\nBOBC4CqAiLgPeHtm3lbxtAVUzVacmTsi4svA70XEWymWB/pt4IbMrOuXytTUNJOTrfuB3k91zKee\nuX5Azra/6emZw5aZL9+Txg1qO25XPYNcR+XnRK3P79Zj6cY6JAnqnAytzqVGys95HcWspnPNhiqp\nNc6l+GG+Hfgm8JnMvK607SeBZQAR8aGI2EfR62QG2BgReyPig6Wy7wTGgB9SXA27Hfitth2FNNhs\nx5Ik9aB6r2gfcqmRzByvLBwRi4E/Bs4DfmU+gUqqT2ZuA86ZY9twxe3LgcsPsZ8x4M1ND1DSYdmO\nJfWq6alJMrccuH/yyaewaNGiDkYktVe90yPWu9TIh4HvZua36w1MkiRJUm/a89SjfPKWzXz0s+u4\ndO3NjI5u6nRIUls1Y9bxWZcaiYiXU1zJfsV8dt7KpRIql7YZGlrQtGVtqvXbshW9Xke76nEpEUlS\nr5iYmDgoEfLqo5ph+cpVHHX8iZ0OQ+qIehPtepYa+e/ARzKzrqVDqrVy9sbKpW2WLVvc9GVtqvXL\nshX9Uke76nEGUklStxsd3cSla29m+cpV7N7xEFdeBGvWnNHpsCSpZ9WbaNe01EhErAJeBbw8IsrL\neS0DpiPilzLzlbVW2MplGCqXthkff7ppy9pU67dlK3q9jnbV41IikqRe4tVHSWqeuhLtGpYa2ULR\nXfwfgRdWPf0TwMPAlfXU2cplGCqXtpmenumLJSWso/vqcSkRSZIkabA0Mkb7XOAGiqVGdgHXViw1\nciKwLDNngG2VT4qIvcBYZv54HvFKkiSpTpVjsCtngq7n+Rs3bgQcvy1Jtag70a51qZFZtr2t3rpa\nbf/+/WzdurXTYUiSJLVU5Rjsx7beyXEvPrOu5997b/F8wPHbklSDZsw63rPuv/8H/NlffpvFS1d0\nOhRJkqSWKo/B3r3j4YafL1Wqnq2+kd4SUr8a6EQbMMmWJElSS0XEauAa4GxgN/DFzLxsjrJLgeuB\nNwMnZeb3KrY9ALwAmKJYYncG+F+Z+e9aGf9cKntKAA31lpD61cAn2pIkSVKL3QTcCbwJOA74ekRs\nz8yrKwtFxAuAb1FMLDzznL0Uj702M7/T4nhrVjlbfaO9JaR+ZKItST1uvpMcSZJaJyJeCZwKvCYz\nx4HxiFgLXABcXVX8+cD7gE3AW+fY5YJWxSqpeUy0JanHzXeSI0lSS50OPJCZYxWPrQciIpaVkm8A\nMvMe4J5SV/O5XBARnwKOBW4F3pWZj7cicEmNM9GWpD4w30mOJPWnco+XRpf0uueeUUZGlrBly30A\nTE9NHtjX/v37ATjiiCNc8uvQVgI7qx57smLbOLVbD9wB/H/ACuBzwJeAn51njJKabGAT7YmJCR58\n8IFOhyFJktQy5R4ve3c9Vndvl9l6y+x56lE+ecs2lt8+zmNb7+TI5x0HuORXA8rdv2cbhz2nzPzl\nirt7I+LdwOaI+BeZ+cNa9zM8PFRPtU3Zz/DwEAsXNrfeZh1HJ+vplzraVU8765ivgU20R0c38em/\nvJWRYw7VM0eSJKm3FTNCNzasd7beMpWPLV/5wiZF2dceB46peuxoiiT7iXnu+4HS358Aak60R0aW\nzLPa+vczMrKEFSuWNqXeRurv9nr6pY521dOuY5mPgU20ARYtGel0CJIkSepv64DVEXF0Zpa7jJ8F\nbM7MvYd43kFXuyNiFXAZcEFm7i89/PJSua31BDQ2to+pqel6njLnfuopu3PnnnnXCcUVx5GRJU07\njk7W0y91tKuedtYxXwOdaEuSJA2qyvHWrljQOpl5d0TcAVwRERcDJwAXAlcBRMR9wNsz87aKpy3g\nud0Qfgz8EjAZEZcBRwFrgb/OzEfriWlqaprJyfknKfUkOs2qs9X77FQ9/VJHu+pp17HMh4m2JEnS\nAKoeb+2KBS11LnADsB3YBVybmdeVtv0ksAwgIj4E/Fbp8RlgY0TMAL+Xmb8fEb9AkVz/qLT9ZuCi\nth2FpJqZaEuSJA0oVyxoj8zcBpwzx7bhituXA5cfYj+jwC80PUBJTdfaqeckSZIkSRowXtGWJB00\nVhNwTVxJkqR5MNGWJB00VnP3jodcE1eSJGkeTLQlScCzYzUlSZI0PybaktTnqruFn3baacDSzgUk\nSZLU50y0JanPVXcL//glQxx33M90OixJXar65JwkqX4m2pI0AOwWLqlW5ZNze3c95traktQgE21J\nkiQdZPnKVcCCTochST3LRFvqQxGxGrgGOBvYDXwxMy+bo+x7gHcBxwH3ABdm5vrStpXA1cBrgUXA\nBuCSzNzQ8oNQV5uYmGB0dNOB+y4H1hq2ZUmSetNQpwOQ1BI3AQ8DL6L4Yf3GiHhvdaGIeD3wYeAt\nwPHALcDXImJJqci1wPOBl5W2/xNwS0R4mWPAjY5u4tK1N/PRz67j0rU3H5R0q6lsy5Ik9aCBS7Q3\nbLiLDRvu6nQYUstExCuBU4H3Z+Z4Zt4PrAXOn6X4+cCnM3NdZj4DXAXMAK8vbT8d+EpmPpWZ+4HP\nUVwte0Grj0Pdrzzuu+hiqmazLUuS1LsGLtGWBsDpwAOZOVbx2HogImJZVdkzStsAyMwZ4G6gPPvN\n14D/FBHHR8RS4FeADZm5rVXBSzrAtixJUo9yjLbUf1YCO6see7Ji23gNZY8p3b6EogvqNoqrYw8C\nv9jMYCXNybYsSVKPMtEumZycZMOGu5zQR/2qPA5zpsay5XLXlm7/c2AMuAD424h4WWburSeA4eHW\ndaAp77vX62i0nnpjGhpacNDfufa5cOHc+62us7r8oL8njdZRg4625X56rQeljla/Z5X1HOoz43DP\nrfzbCl3WjiUNABPtkh/96BH+5C/+jmsvv4A1a87odDjSfDzOs1exyo6m+JH9RI1lN0XEkcDbgH9Z\n0b308oi4CPh54K/qCWpkZMnhC81Tv9RRSz0TExNs3LgRgEce+WFd+162bPFBf+eqf8WKpTXHN1f5\nQXpPWqDr2nI/vdaDUkc737NDfWbUuo9W60A7ljSgTLQrHDlybKdDkJphHbA6Io7OzHI307OAzbNc\nuVpHMbbz8wARMUQxLvQGijkcFlDxOVGaofiIRoIaG9vH1NR0I089rOHhIUZGlvR8HfXUs379XVx8\n1ZdZvnIVj229k+NefOacZauNjz990N/ZjI3tY+fOPYfcfqjyg/ieNKOOKl3XlvvptR6UOqrbaqsc\n7jPjULrp9WpGHZIEJtpS38nMuyPiDuCKiLgYOAG4kGIWYiJiC3BeZt5G0Z30xoi4kWLd3UuAp4Gv\nZ+YzEfEt4Lci4q3ALuB9wATw7XrjmpqaZnKydT/Q+6mOWuqZmpo+MOv37h0P17Xv6emZg/42Wn8t\n5QfpPWm2bmzL/fRa93MdlevcZ25paf3VcVTWDdQ1JK9f3hNJAhNtqV+dS3ElazvFj+prM/O60rYT\ngWUAmXlrRHwA+BLFGrt3Aq8rLQ8E8Cbg4xSzF/8zYBPwi5lZPemSpNawLatu5XXuG+nx0sy6d+94\niCsvwiF5kgaSibbUh0rjMM+ZY9tw1f3rgevnKPs48F+aHqB6Uieukg0627Ia1WiPl2bWLUmDzERb\nklSTTl4lkyRJ6iUDuQ7B/v37D7oaMzM9zfbtj3YwIknqDeUrVUc+7/hOhyJJktS1BvKK9tat9/O7\na/+URUtGAHhm3y6+tX4Xi5euIHOLY4kkdSW7bkuSJPWGgUy0gQNJdtnipSs6FIkk1cau25IkSb1h\nILuOS1Kvsuu2JElS96v7inZErAauAc4GdgNfzMzL5ij7YeBtwNHAg8DHMvP/bzzc+du/fz8PPvhA\nJ0OQJEmSJPWxRq5o3wQ8DLwIeC3wxoh4b3WhiLgAeEupzPOAjwCfiYjTGg22GbZuvZ8bv357J0OQ\nJEmSJPWxuhLtiHglcCrw/swcz8z7gbXA+bMUvxt4c2b+IDNnMvMmYBfw8vkGPV+Ox5YkSZIktUq9\nXcdPBx7IzLGKx9YDERHLMnO8/GBmfrt8OyIWA78KTAJ/N494JUmSJEnqavUm2iuBnVWPPVmxbbxq\nGxHxp8DbgQeAN2Tmj+uss2kmJiYcny1poE1PTbJly32MjCxhy5b7Oh2OJElSX2rG8l4LSn9nZtuY\nmedHxG8C/wn4ekT8bGZurHXnw8PNmxj9nntG+exN/4uRY1bPWWZoaAELFzZ3MvbyMTTzWKyj++tp\nZx1SrfY89Sg3fHUbX7htt0uESZIktUi9ifbjwDFVjx1NkWQ/MdeTMvMZionQ3kRxdfs9tVY4MrKk\nzhAPva/q9bOrLVu2mBUrljatzur6W806uq+edh2LVKvyEmG7dzzc6VAkSZL6Ur2J9jpgdUQcnZnl\nLuNnAZszc29lwYj4a+BvMvO/Vzw8Deyvp8KxsX1MTU3XGebc+zqc8fGn2blzT1PqKxseHmJkZElT\nj8U6ur+edtYhSVKnTE9NkrkF4MBfSRp0dSXamXl3RNwBXBERFwMnABcCVwFExBbgvMy8DfjfwPsj\n4jZgE/A64OeAj9VT59TUNJOTzUlSakl2pqdnmlbfbPW3at/W0b31tOtYJEnqhD1PPconb9nG8tvH\nHZIiSSWNjNE+F7gB2E6xXNe1mXldaduJwLLS7f8GHAHcQrGO9g+Bt1fORi5JkqTe55AUSTpY3Yl2\nZm4Dzplj23DF7Wng8tI/SZIkSZIGglMWS5IkSZLURCbakiRJkiQ1kYm2JEmSJElN1MhkaJKkNpmY\nmGB0dBPQvmVzKpfqATj55FNYtGhRW+qWJEnqBybaktTFRkc3cenam1m+clXbls2pXKpn946HuPIi\nWLPmjJbXK6m/VJ6084SdpEFjoi1JXa4Ty+aU65SkRpVP2nHLZk/YSRo4JtqSJElqieUrV3U6BEnq\niIGZDG1iYuKw4xtnpqd58MEHmJiYaFNUkiRJkqR+MzCJ9ujoJq649guHLPPMvl185dt5YOIhSZIk\nNc/ExAQbNtzFhg13eWFDUl8bqK7ji5eu4Ok9Ow9Z5siRY9sUjSR1v8rJjObqFTTbLOULFy5uS3yS\nekt5gkfAcduS+tpAJdqSpPpUzkA+16zns81SfuaZrZ8dXVJvcty2pEFgoi1JOqRaZj13lnJJkqRn\nmWhLfSgiVgPXAGcDu4EvZuZlc5R9D/Au4DjgHuDCzFxfsf2XgCuAFwHfA96Xmd9o6QFIAmzLUr+o\nsy0vBa4H3gyclJnfq9h2VGnbq4Ep4OvAb2TmM609Akn1GpjJ0KQBcxPwMMUP6tcCb4yI91YXiojX\nAx8G3gIcD9wCfC0ilpS2/xTwaeAC4CjgauAjETHchmOQZFuW+kWtbfkFwF3AfmBmlv18ElgCvAw4\no/T3Y60JWdJ8eEVb6jMR8UrgVOA1mTkOjEfEWoof2FdXFT8f+HRmris996pSudcDXwLeA3w+M/+2\nVP4zpX+SWsy2LPWHOtvy84H3AZuAt1bt51jgDcBpmbmz9NhHgS9FxMWZOdXaI5FUD69oS/3ndOCB\nzByreGw9EBGxrKrsGaVtAGTmDHA3UJ7J6qeBHRHxzYh4KiK+GxFrWhi7pGfZlqX+UHNbzsx7MvNr\nc+znp4Din45KAAAgAElEQVTJzByt2s9y4KRmBixp/ky0pf6zEqhex+7Jim21lD2mdPufA78CXFS6\nfTfw1Yhw7aYWmpiYYP36Yp3ZuZbU0kCwLUv9oZ62fLj97JpjP8cgqavYdVwaDAtKf2cb7zVb2ZmK\n25/LzLsBIuJS4Ncoro7VNYnS8HDrzuuV993rdZT3v3HjRi6+6sssX7lqziW1utnw8FDfvSetrqeO\nfXe0LffTa93PdbT6fWrE8PAQCxcOHRRb9WP98p7UqJ62XIu69tOs16Ge/ZTf72bW62dS99TRrnq6\nsC3PyURb6j+P89wz20dTfAk/UWPZTaXb26k4e56ZeyLiCYrJluoyMrKk3qfUrV/qgNqW1OpWIyNL\nDrxO/fSetKueCl3Xlvvpte7nOjrwf/WwRkaWsGLF0oNiKz9Web8dcXRAPW35cPs5KiIWlIaHwLNX\nxB+vJ6BmvQ717Kf6/W53/d1eT7/U0a56uvFzrpqJttR/1gGrI+LozCx3KTsL2JyZe2cpewbweYCI\nGKIYS3ZDaftmijFhlLYvpfix8GC9QY2N7WNqarrep9VkeHiIkZElPV9HuZ5eNza2j7GxfX31nrTr\nWKp0XVvup9e6l+uYmJhg8+Z7WbZsMePjT/Pyl7+CRYsWHdg+Nrav6XXO19jYPnbu3HNQbOXH+uE9\nqaxjFvW05UrVV6k3UFwJP41i+Ed5PzuBrCfWZr0O9fxfK7/fzdCPn/+9Xke76ulwW66LibbUZzLz\n7oi4A7giIi4GTgAuBK4CiIgtwHmZeRtwLXBjRNxIse7uJcDTFOtyAlwHfDEi/gL4DvAHwFbgu/XG\nNTU1zeRk6z7c+6mOXjc1NX3gy6+f3pN2v/fd2Jb76bXu5To2btzIpWtvZvnKVeze8RBXXjTFmjVn\nHFRvtym/FpWxVb8+vfyeHEoNbfk+4O2ltly2gGe7l5f3syMivgz8XkS8lWKZr98GbsjMug6qWa9D\nPf/XWvHa+5nUfXW0q55e+D3W+5dOJM3mXIov8u3AN4HPZOZ1pW0nAssAMvNW4AMUy//sAH4OeF1m\nPlPa/lWKyZNuKG0/tbS9uz/ZpP5hW9asysNLlq9c1elQVJtDteWfpNSWI+JDEbEPuI/iivbGiNgb\nER8slX0nMAb8kOKq9u3Ab7XtKCTVzCvaUh/KzG3AOXNsG666fz1w/SH2dR3F1TBJbWZblvpDrW05\nMy8HLj/EfsaANzc9QElNZ6ItSZKklpmemjywVKFLFkoaFCbakiRJapk9Tz3KJ2/ZxvLbx3tyyUJJ\naoRjtCVJktRS5THlRz6v7tUhJaknmWhLkiRJktREdh2XJEnqARMTE4yObnKcsyT1ABNtSZKkHjA6\nuolL197M3l2POc5ZkrqcibYkSVKPKNbNXtDpMCRJh2GiLUlqmvIyPsPDQ4yMLGHVqpcyNORXjdQO\ndi2XpO7hrx9JUtNULuOze8dDfPySczn11DWdDksaCHYtl6TuYaItSWqq8jI+ktrPruWS1B1MtCWp\nC5S7fA4PD/HIIz/sdDiSJEmaBxPtOWzYcBcAa9ac0eFIJA2CcpfP5StX8djWO+32KQl49iQc4Nhr\nSeohdSfaEbEauAY4G9gNfDEzL5uj7DuB9wI/AfwA+Ehm/nXj4UpS/yp3ud694+FOhyKpS3gSTt2m\nkZM/5Ykyy04++RQWLVrUkvikbtHIFe2bgDuBNwHHAV+PiO2ZeXVloYj498DvA68rlX8r8KWIOCkz\nH5hX1JKkrjc9NcmWLfcxNTUN+MNKapQn4dRNGjn5Uz1R5pUX2WtU/a+uRDsiXgmcCrwmM8eB8YhY\nC1wAXF1VfAnwgcy8vXT/UxHxMYor4Q/MK+oWmp6eInMLL37xSzjiiCM6HY4k9aw9Tz3KDV/dxvKV\nu/1hJUl9pJGTP06UqUEzVGf504EHMnOs4rH1QETEssqCmfnnmXl9+X5EHAUsB37UaLDt8PT4Dq7+\n3K1s3Xp/p0ORpJ5X/mFVzIQsSZI0GOpNtFcCO6see7Ji26HcAPxjZn6nzjrb7siRYzsdgiRJkiSp\nRzVj1vHyYo0zs22MiIXAZ4GXAT9b786Hh+s9F9Cc/QwNLWB4eIiFC+dff7nuZh2LdfRGPe2sQ5Ik\nSVL3qDfRfhw4puqxoymS7CeqC0fEYuCvgcXAqzKz+mr4YY2MLKn3KU3Zz7JlixkZWcKKFUubUn8j\nMVhH6/XTsUiSJEnqDvUm2uuA1RFxdGaWu4yfBWzOzL2zlP8C8DRwTmbubyTAsbF9B2asnY877lhf\nV/nx8acZG9vHzp175l338PAQIyNLmnYs1tEb9bSzDkmSeknlck/Dw0O86lVndzgiSWquuhLtzLw7\nIu4AroiIi4ETgAuBqwAiYgtwXmbeFhH/GTgZOKXRJBtgamqaycn5JynT07P2bD9k+WbVXdbs/VlH\nb9TTrmORJKlXVC/3dMPIEl760pd3OixJappGxmifSzGx2XZgF3BtZl5X2nYiUO5r/TZgNfBkREAx\nlnsG+HxmvmM+QddjYmKC0dFNTE5O1vyc6ekpHnzwAV784pe0MDJJkqTB5XJPkvpZ3Yl2Zm4Dzplj\n23DF7dfOI66mGR3dxK9/6A9546uj5uc8Pb6DG7/+A1avfhFnnWVXJkmSJElS7Zox63jXa2S5rsVL\nV7QgEkmSpOapHOtc/ttrpqcm2bx5M2Nj+zjppJNZtGhRp0OSpHkbiERbkiSpH1WOdX5s650c9+Iz\nOx1S3fY89SifuHEbcBdXXjTNmjVndDokSZq3gUi0p6en2L790U6HIUmS1HTlsc67dzzc6VAatnzl\nqk6HIElNNdTpANrh6fEd3PrdTbWX31P3ct+SJEmSJAEDkmgDLFoy0ukQJEmSJEkDYCC6jkuDJiJW\nA9cAZwO7gS9m5mVzlH0P8C7gOOAe4MLMXD9LuV8C/gr415n5D62KXdKzbMuSJPWmgbmiLQ2Ym4CH\ngRcBrwXeGBHvrS4UEa8HPgy8BTgeuAX4WkQsqSp3JPAJYLy1YUuqYluWJKkHmWhLfSYiXgmcCrw/\nM8cz835gLXD+LMXPBz6dmesy8xngKmAGeH1VuY8A3wCeaFngkg5iW5YkqXeZaEv953Tggcwcq3hs\nPRARsayq7BmlbQBk5gxwN3BgfZiIOIXiKtkHgAWtClrSc9iWJUnqUSbaUv9ZCVRPnf9kxbZayh5T\ncf9a4Lcy80kktZNtWZKkHuVkaNJgKF+9mqmx7AxARPwasCAzPzXfAIaHW3der7zvXq6jlbF3i+Hh\nIRYubN5xtuN9b1c9dey7o225n17rXqljED4bKjX7c6Jyv5V/W2HQ3itJh2aiLfWfxzn4KhbA0RQ/\nuKvHZc5VdlNErAR+F/iFZgQ1MrLk8IUGuI52xN5pIyNLWLFiaUv22w4deI+6ri3302vdK3UMwmdD\n2fTUJI888sMDx3zaaaexaNGiptYxSK+npM4y0Zb6zzpgdUQcXdFF9Cxgc2bunaXsGcDnASJiiGJc\n6J8B51D8UP9GRJSvoq0A/kdEfC4zL6gnqLGxfUxNTTd0QIczPDzEyMiSnq5jbGxf0/fZbcbG9rFz\n556m7a8d73u76inXUaXr2nI/vda9UscgfDaU7XnqUT5x4zaWr3yc3Tse4uOX7OP0089oyr472I4l\nDSgTbanPZObdEXEHcEVEXAycAFxIMQsxEbEFOC8zb6MYs3ljRNxIse7uJcDTFEsDLaCYnbjS7cB7\ngb+rN66pqWkmJ1v3A73X62hl8tItWvnatfp9b2c9Zd3Ylvvpte6VOgbhs6HS8pWrOOr4E4HWvEft\nbseSBpeJttSfzgVuALYDu4BrM/O60rYTgWUAmXlrRHwA+BLwfOBO4HWl5YEAtlXuNCImgScyc1fr\nD0EStmVJknqSibbUhzJzG0V30dm2DVfdvx64vsb9vnj+0UmqlW1ZkqTe5PSIkiRJkiQ1kYn2IUxO\nTrJhw11MTEx0OhRJkiRJUo/o+0Q7c0vDz/3Rjx7h1z/0h4yObmpiRJIkSZKkftb3ifZ8HTlybKdD\nkCRJkiT1ECdDkyS11cTExIGeQieffAqLFi3qcESSJEnNZaItSWqr0dFNXLr2ZgCuvAjWrDmjwxFJ\nkiQ1l4m2JKntlq9c1ekQJEmSWsZEW5IkSV1lemryoAltHWaiXlY5ZAr8/zwoTLQlqY38spWkw9vz\n1KN88pZtLL99nN07Hur5YSYRsRq4Bjgb2A18MTMvm6Pse4B3AccB9wAXZub60ra/B/4lMAksKD1l\nS2auaekBaF7KQ6aWr1zVF/+fVRsTbUlqI79sJak2y1eu4qjjT+x0GM1yE3An8CaKBPrrEbE9M6+u\nLBQRrwc+DPwCsAm4APhaRLwkM/cBM8DbM/PzbY1e89Zn/59VA5f3kqQ2K3/ZOk5ZkvpfRLwSOBV4\nf2aOZ+b9wFrg/FmKnw98OjPXZeYzwFUUyfXrK8osmOV5krqMV7QlqQVq6SJeOQaxciyiJKmvnA48\nkJljFY+tByIilmXmeMXjZwA3lu9k5kxE3A2cCXyp9PCbIuL9wAuB24F3ZubWlh6BpLr1baI9MTHB\n3XevZ+vW+zsdiqQBVEsX8coxiI9tvZPjXnxmh6KVJLXQSmBn1WNPVmwbr6HsMaXbm0vl30zRM/VP\ngL+JiJdn5mQzg5Y0P32baI+ObuL8Sy5nYt8YI8es7nQ4kgZQLeOxymV273i4TVFJkrpAufv3TI1l\nZwAy892VGyLifIpE/FXAt2qtfHi48dGj83lu5T4WLpx/DM2IpVX1TExMcO+9Rc+2738/n7Pf8vFX\n1lH5HIBXvKI5E6b2wuvVjXXMV98m2gCLl67odAiSJOwmL9WrPPzE9tIXHufZK9JlR1Mkz0/UWHYT\ns8jM8Yh4EviJegIaGVlST/GmPbdyHytWLG3Kftqh1nomJibYuHEjAJs3b+YTN97F8pWrntNrbbbj\nHxlZQua9XHzVlw/0hrvho0s488zm9Xbrtter2+uYr75OtCVJ3cFu8lJ9ysNP9u56zPbS+9YBqyPi\n6Mwsdxk/C9icmXtnKXsG8HmAiBiiGON9Q0QsB64APpqZ20vbjwGeD9Q1RntsbB9TU9MNHczY2L6G\nnle9j5079zT8/OHhIUZGlszrOFpRz/r1dx1IlMvfdbP1Wqs8/so6xsb2HdQbbr6vU6PH0c31tLOO\n+TLRliS1hd3kpfoUKxM4wXSvy8y7I+IO4IqIuBg4AbiQYkZxImILcF5m3gZcC9wYETdSrKF9CfA0\n8PXMfCYizgb+uNRlHIq1ue/OzH+sJ6apqWkmJxtLUpqR3Myn/lbsp1n1TE1N1/RdN9v+pqamD3pt\np6cm2bx584HHZptUtV7d9np1ex3z5fJekiRJUmudS5Fgbwe+CXwmM68rbTsRWAaQmbcCH6CYYXwH\n8HPA60pLfQG8geLsy/eABykump3TpmNQGxU9wTbz0c+u49K1Nx+0kol6g1e0JUmSpBbKzG3MkRBn\n5nDV/euB6+co+whF0q4BcLhJVauXEoXmXPlWc5hoS5IkSVKPqVxKFJhzOVF1Rt2JdkSsphgPcjaw\nG/hiZl42R9mlFGfk3gyclJnfm0eskiRJkqSSWpYSVWc0Mkb7JuBh4EXAa4E3RsR7qwtFxAuAu4D9\n1LZGoCRJkiRJPa+uK9oR8UrgVOA1mTkOjEfEWuAC4Oqq4s8H3kex7t9bmxCrJEmSpDarHAvs2u5S\nbertOn468EBmjlU8th6IiFhWSr4ByMx7gHtKXc0lSZIk9aDKscDl9aElHVq9ifZKYGfVY09WbBun\nz+wd+zGZW5xUQFLDpqcmD1wB8EqAJKkX1bI+dC0qvxNh8GbJrj7+0047DVjauYDUMs2YdXxB6W9L\nxmEPDze21Hejz6s0NLTgwN+FCxvfXzmWZsRkHb1TTzvrUHcr1sLcxvLbx70SIEkaaJXfib0+S3b1\n8lq1nEyvPv6PXzLEccf9TCvDVIfUm2g/DhxT9djRFEn2E02JqMrIyJK2Pq/SkiXF2bVlyxazYsX8\nzzQ1Iybr6L162nUs6m7NuhIgSVKv65eZsquX16r1ZHq/HL8Ord5Eex2wOiKOzsxyl/GzgM2ZufcQ\nz2v4avfY2D6mpqYbet587ds3wcz0NKOjW3jssZ0Nd2sZHh5iZGRJw8diHb1ZTzvrkCRJUvtVJs2e\nTFeluhLtzLw7Iu4AroiIi4ETgAuBqwAi4j7g7Zl5W8XTFvBs9/K6TU1NMzlZf5IyNTXN03uqh5PX\nZ3p6hmf27eKmb23hta/dOO9uLY0ei3X0dj3tOpZKda53/x7gXcBxwD3AhZm5vrRtMXAF8MsUA4ju\nBC7KzNGWH4Qk27Ik9bHpqUm2bLnvwIUZ53HpL40M8DyXIsHeDnwT+ExmXlfa9pPAMoCI+FBE7APu\no7iivTEi9kbEB+cfdnsdOXJsp0OQ6lXrevevBz4MvAU4HrgF+FpElC+TXwn8K4of+ScADwFfaXXw\nkg6wLUtSn9rz1KPc8NVRLrr623zkU3fwx3/xrU6HpCaqezK0zNwGnDPHtuGK25cDlzcemqRG1Lne\n/fnApzNzXem5V5XKvR74EvAU8L7M/FFp+9XAeRFxfGZub8sBSQPKtixJ/c+u5/2rL6csnpiYsOuF\nBtkh17uvKntGaRsAmTkD3A2cWbr/O5n57Yryq4CneXZZP0mtY1uWpC4zMTHBhg13sWHDXeYbOqRm\nLO/VdUZHN3HFtV/odBhSp9Sz3v1cZatXFyAiVgB/CFyVmRPNCVXSIdiWJanLVM403oklOyuXFDPR\n7259mWgDLF66Yt6ToUl9pJ717hdUl4uIFwD/E7gL+K+NBNAva5bXWodrnNdmeHiIhQsbe63a8b63\nq5469t3RttxPr3U31TExMcG99xY/nF/xilNYtGiRnyFV5vNZUX5+5d9W8D0bHJ1csrPTib5q17eJ\ntjTA6lnvfq6ym8p3IuIlwDeArwIXlLqk1q1f1iyvtQ6XXavNyMgSVqxYOu99tEMH3tOua8v99Fp3\nUx133rmZi6/6MgA3fHQJZ555pp8hVZrxWVHej1SviYkJNm7cCHTHVeR6E/2JiQnuvHPzgSVnTz75\nlIaXLVbtTLSl/lPPevfrKMZ2fh4gIoYoxoX+Wen+SuBW4M9KExw2rF/WLK+1jrGxfS2Jo9+Mje1j\n5849DT23He97u+op11Gl69pyP73W3VTH2Ng+lq9cdeD2zp17/AypMp/PCuhoO1YfuPfe3riKPD01\nedCJgHJCfe+9m7j4qi+zfOUqdu94iCsvYt7LFuvwTLSlPlPDevdbgPNK691fC9wYETdSrLt7CcUE\nSbeUdncFcPt8k2zonzXLa62jlclIP2nGe9auterbVU9ZN7blfnqtu6mOys+L8nP8DDlYs96vdrdj\n9Y9Odhev1Z6nHuWTt2xj+e3jz0moK2c3V3uYaEv96VzgBor17ncB11asd38ipfXuM/PWiPgAxfI/\nzwfuBF6Xmc+Uyr4NmIyIX6borloe8/lrmfnn7ToYaYDZlqVZVE4IZTdY6Vkm1N3DRHsOM9PTbN/+\naKfDkBpS63r3pfvXA9fPUdbPCKmDbMvS7MoTQgF2g5XUlfzincMz+3bxrfW7Oh2GJEmSZlEe1y7p\nuSrHa3//+9nhaAaTifYhuESYJEmSpF5TOV67mydw62cm2pKkjqscbwmOuZQkab56YQK3fmaiLUnq\nuPJ4S5cekSRJ/cBEW5LUFZwpVZKk56ocb125Tra6m4m2JEmSJFWYmJjgnntGGRlZwpYt93U0Fsdb\n9yYT7RpMT0+RucUxg5IkqWW8ajW78uuyf/9+AI444ghfH7Vc5ZCmbkhuHW/de0y0a/D0+A6u/tyt\nRJzkmEFJktQSXrWaXfl12bvrWxz5vOO6JvFR/zO51XyYaNfoyJFjOx2CJEnqc/6wn12xZvYClq98\noa+PpJ5goi1J8+CyVJLUHSo/j/0sltRpfZVolz9gy2N4JKnVXJZKkrpD+fMY8LNYUsf1VaI9OrqJ\nX//QH/LGV0enQ5E0QFyWSpK6Q9HFXJI6r68SbXAstaTOccbg5rNrviSpXSq/c/we13z1XaItSZ3i\njMHNZ9d8SVK7dNuSXq1WfTIbPKHdTCbaktREzhjcfHbNlyS1QnWimblloL7HK08sAJ7QbjITbUmS\nJEkDpzrRHISr2NU8md06JtqS1IDyWXDHcDWu3jHtleUBTjvtNGBpq8KTJPWJQ833UZloDsJVbOeT\naZ++S7Snp6fYvv3Rlux3dPRe9u/fz0/91OmOXZAGXPks+N5djw3c2e9mqXdMe2X53Tse4uOXDHHc\ncT/TpmilxlUvP3rEEUccdNsfu42bLWmoPil38smnsHDh4o7E18u+8c2/53//010APPqjB2HhyR2O\nqHHO9/Es55Npn75JtDdsuIvMLTw9voNbv/sgI8esbur+n3rsB/y3G37A4qUruPbyCwa2cUp6VtHV\nbEGnw+hp9Y6Fs4ubelHlibkjn3fcgYmWKm/7Y7cxsyUN1SflrrwIzjzT17det925nocWrAHgR08+\nwtIeX9jH749nDdI49E7qm0S70qIlIy3Z7+KlK1w+TJJazG5t6kflE3PLV77wwA/cyttq3GxJg0mV\n5uJ3jNqlLxNtSVLvslubJKlV/I5Ru5hoS5K6jt3a1EsqJ1rav38/w8NDrFw5wqpVL2VoyJ9anVS+\nejk8PMTIyBLfEwF+x6g9+uqT5sEHH2hrfRs2FBNEOF5bkqTBVTnRUnnsNcDHLzmXU09d0+HoBttz\nJ1H0Pekm1RPXwcEzgku9rK8SbUmSpE6ovEK2fOULOx2OKjheu3tVnggBBn5G8G5yqCXRVBsTbUmS\nJEkd0aoTIZWJopOe1c8l0eavLxLt7Y9tZ9369S1ZP1uSJOmPr/sznnhyjP37n+Z9v/kOjj565SHL\nT09NsmXLfUxNTfsjv8tUj6mHYj3z6it2XtHrbdVDOpz0rH6znQSxXdSuLxLtb3zz23zhL/+KR58c\nb/r62dX2jv2YzC2e0ZGkLuQPALXKPd9/lKnn/yt27Rjlxz/+8WET7T1PPcoNX93G8pW7/ZHfZeYa\nU199xc4rer3PSc+az3ZRu7oT7YhYDVwDnA3sBr6YmZfNUfY9wLuA44B7gAszc33j4c5twfBwy9bP\nlnpNs9ppRCwC/gg4B1gEfBt4Z2Y+2fKDkBrQbz8AbMu9zR/53avWMfXN6tZsW1YvqHWNcec9qM1Q\nA8+5CXgYeBHwWuCNEfHe6kIR8Xrgw8BbgOOBW4CvRcSShqOVVKtmtdM/ANYA/xcQFJ8Zn2518NJ8\nlH8ALF+5qtOhNINtWWqScnf+DnXlty23wcTEBBs23MWGDXc5ZKMBxeR0m/noZ9fxx3/xrU6H0/Pq\nuqIdEa8ETgVek5njwHhErAUuAK6uKn4+8OnMXFd67lWlcq8HvjTfwDtlZnqaBx98gImJCaAY27Nh\nw112T1TXaFY7jYgvA+cBb8nMbaXtHwI2R8Txmbm9PUd0eN/+zneYmZlkz55nOPUVp/CCF7yg7n3M\nNWavfHvx4n/GkiULGRvbx9DQsF/gHVb+wTwysoSxsX2cdNLJDX0Gd3NX80Fsy1Irlbvz7931WFu7\n8tuWa1e93NfhPpMnJibYuHHjgfuZW/jkLZsdlz0Ps/XEqfVKtw5Wb9fx04EHMnOs4rH1QETEstKH\nR9kZwI3lO5k5ExF3A2fSw4n2M/t2cePXb2f16r8k4iS2br2fqz93K9defkFPd09UX2lWO70beB6w\noWJ7RsS+0vNuaeEx1OUTN3yZ4WNfCcCaTV/isosuqHsfs43ZO9xtv8A7p/yD+Qu37S51EZ9u6DO4\ny7uaD1xbllqt6OmyoN3V2pZrVL3u+eE+k++999nPcODAd7NDNpqr8n2p/P1TfWJkeHiIV73q7E6F\n2XXqTbRXAjurHnuyYtt4DWWPqbPOrrN46YqD7h85cmyHIpFm1ax2uhKYmWX7TrqsHS8+coRFRxVX\nsRfwRMP7qR6zd7jb6qxmjRHr4rFmA9eWpT5lW65D+TO5Momr7mk2PDzEypUjbNly30Gf4X43t85s\nV7pnWwf9hpElvPSlL5+zpyDM3VOhm3uZNaIZs46XTwvO1Fi2lnIHDA8ffhj50PACFgAT+4oThdV/\nZ3tsvtsefvhBhoYW8PDDD7J37McMDw+xcOHssZaPoZZjaZR1dF897ayjBs1sp3W3Y2jt67Cg4uLE\nY9u3cc89G+YuPIfvfz/ZveMhAPbu2k75EA93e++ux2ou287b3RpXK2LcveMhvv/95QwPDx30PlY+\nPpfq8sPDZz3ns9y2/KxB/mxeUPFB873v3cfExNMH7s/2+THb/+/D/Z+v9Tnup/n7Kd7Hgz8vavl8\nqFTn/6mub8sLFiw4aK+H+4481LZablfff/zBDVx9/wRHjmzgyUeTxUtXcOTIsc+5/fxVp7akfp9f\n2+3yrP1lmzdvZnz8aTZv3szaz/zNc96zvWM/5qJf+UVOOullVNuy5b4Dz9k79mP+9IoLOf30g3s0\ndNl38iHVm2g/znPPmB1N8WpXX0aaq+wmardgZOTwc6dd8O63c8G7317HbpvryhrL1XIs82Ud3VdP\nu46lQrPa6eMUX97HUEzgUraitK0eNbXlRv3NF9ZW3HtDQ/v4uZ/7Gd797ubEo86p932sp7xtubXt\nuFI3fjZ/5bMfK9167meMnx/9qYnva0+25T+66ncq7jX23arB9bM/+yre/e531PWcos3V9pwOfCfX\nrd50fR2wOiKOrnjsLGBzZu6dpeyBUxARMUQxRuWfGglUUs2a0U5vB7ZSdEer3P4KiuVE1rUmdEkV\nbMtSf7AtSwNowcxMfT1AI+I24F7gYuAEiokXrsrM6yJiC3BeZt4WEb9AMZnDv6VYA/ASipkSIzOf\naeIxSKrSrHYaEX9AaRkSYB/FEiJ7M/NNbT8oaQDZlqX+YFuWBk8jHdDPpfiA2A58E/hMZl5X2nYi\nsJ59TLQAAAkoSURBVAwgM28FPkAxw/gO4OeA15lkS23RrHb6OxRn0TcC9wO7gF9r0zFIsi1L/cK2\nLA2Yuq9oS5IkSZKkubV2+lBJkiRJkgaMibYkSZIkSU1koi1JkiRJUhOZaEuSJEmS1EQm2pIkSZIk\nNZGJtiRJkiRJTbSw0wHMJiJWA9cAZwO7gS9m5mVNrmMaeAaYARaU/t6QmRfMc7+/AHwW+GZmvrlq\n238EPgj8CyCBD2bm3zarjoh4K/ApiuOCZ4/rZzL/T3v3HyNHXYdx/F0bJCGIEY38MFeMP3i0EttU\nqLGAqSX+QfQPjGLU0gQBESOBYNFaWqxYkCK0NBWtBYlt5Ie/QECpoSqICTFiqYqR8DFSLzUqtgLS\nVAtUc/7xnQvL3u7ezu53l53xeSVNezPbee57M89dZu47s7G9ZMYsYD3wTuA54B7gwojYK2lusW4u\n8HdgU0Ss62EcLTOAVwB/Ap5pGsfKHnPmAGuB44H9wP3ABRGxW9Ii4ErgTcAu4MqIuCVTxoXFdu9r\nMZYlEXFb2ZyGvGtJ++MlxcdZxpHTMHpc5LjLnTNq0eU69rjIdJdxj7vIGHiPO+XgLk+XN/I9Bne5\n14yqdbkuPe6QU5kuj+SJNnAb8CvgQ8ARwFZJj0fE+owZE8CxEfHnXBuU9GngLOAPLdbNBTYDp5EO\njA8A35d0bET8NUdG4f6IWFTyU2/lB6R9MEYq5h3ANZIuKNZtAk4lHXjbJO2MiDtyZABXABMRcUi/\ng5D0UtI3mA3F53sY8D1go6RPAncC5wO3AicDd0l6NCJ2ZMj4arFsPCJe1+9YGvLmAktIxzCSjsox\njgEYRo/BXZ5O5btcxx4Xme7y89zjzobR47Y5uMud8qrSY3CXe8ooVKnLle/xNDmV6fLITR2XdDzw\nVmBZROyLiMeAdcC5maNmFH9y2g/MBx5rse5s4O6IuCciniuuhvwOOCNjRhaSXk4q6PKI2F98o9pC\nujL2HuAg4Ipi3a+Br1Ny/0yTkdMhpKucayLiQEQ8AdwOHAcsBiIithT75KfAXcA5GTOykjQD2Ei6\nujcp1ziyGWKPwV1uq0ZdrlWPwV1uwT1uYxg97iInl1p1uSo9Bne5z4ws/DM5a05Wg+ryKP5Gex7p\nCsXehmU7AEk6NCL2Zcy6StIC4GXAd4FPRcS/et1YRFwHIKnV6rcBP2xatgM4IWMGwJikbaQpFk8C\nqyLi5pIZTzP1IBoD/kIax8MRMdGwbkeL1/eSMavIAJghaQvwbmAmcCNwaUT8t2TOP0nTfYB0EAFn\nAt8ijaX5itQO4IOZMm4tFh0m6XbSlbBngHURcW2ZjAbnkX4Y3AJcXiybR4ZxZDbMHoO73C6jFl2u\nYY/BXW7FPW6dMfAed8hxlzurSo/BXe4nAyrS5br0eJqcynR55H6jDbwSeKpp2ZMN63L5BbANeAPw\nDtL9Kl/JuP1m7cb1qowZe0hTXi4mTQlaAXxD0sJ+NlpcBT2fNOWk3TgOz5RxOekemAdIU5zGSFf5\nzgAu7WP7syQ9C/we+CVwGZn3SZuMvcDDpKvGR5GmJa2SdGYP2z8C+DzwiaZVwzi2yhpWj8Fd7lrV\nu1yHHhcZ7vJU7nGXhtHjphx3ufX2q9RjcJf7UdkuV73HHXIq0eVR/I12K5NTUCY6vqqEiDix8UNJ\ny0jz7j8WEQdy5Uxj8sb9LCJiK7C1YdG3Jb0P+Cjws162KelE0jSJZRFxr9IDJ5r1NY6GjM9ExH3F\n4pMbXrJd0heB5aQilBYRu4CDJb0euB74ZpuX9jyWFhk3RcRioPGenh9L+hppn2wuGbEWuDEiQumB\nJp1kPbYyyd5jcJe7VYcu16TH4C5P4R53Zxg9bspxl9ureo/BXe5KVbtchx63yalMl0fxN9p7mHql\n4HDSoP4xwNxx0hSKVw9o++3GtWdAeZPGgaN7+Y+S3gvcTXqK4OTVyHbjeCJjRivjwJG9ZDQq7kla\nAXyY9CTG7PukMUNSq6vE45TcJ5JOARYAq4tFjfc/vVjHVicvVo/BXZ6ibl2uao/BXS5hHPf4BYbR\n4w45rYzzf9rlCvYY3OXcxhnhLtetx805VenyKJ5obweOkdQ4VWI+8EhE/DtHgKS5kq5pWjybNKWi\n6ycUlrSddN9CoxNIUyCykPRxSac3LX4zsLOHbS0gPTzh/U33oGwH5khqPHZ6Gke7DEmLJF3S9PLZ\npAKVzXiXpEebFk8Uf35CutemUemxTJOxUNJ5TetmU36fLCb9kNolaQ/wEOk+m92kh370PY7MBt5j\ncJe73Fblu1yjHoO7PIV73NW2Bt7jTjnu8hRV6zG4yz2rWpfr0OMucirR5ZGbOh4Rv5H0ILBG0lLg\nNcBFwNUZY3YD5xZfxPXAa4EvkN6vblBTe24AHpR0KnAvace+EbgpY8bBwAZJO4HfAqeTHoc/v8xG\nJM0sPt9lxRP2Gm0l3RexUtLVpCdYng18hBKmyXgK+JykceA7pPcTXAp8qUxG4SHSwxLWkO7pOBRY\nBfwcuBm4TNJZxb9PIX293p4x4wCwVtIfSW87sYj0IIclJTMuAlY2fDxGug9qDqnHyzOMI5sh9Rjc\n5Y5q1OW69Bjc5Vbc4w6G0eMuctzlF6pUj8Fd7lNlulyjHk+XU4kuz5iYGL1bRiQdTTpIFgJPAxsj\nYnXH/1Q+4yTSgXUc6Ul1m4EV/dw/Imk/6SrLQcWi/9DwfnWSTgOuIj397xHSdI4HMmdcQnra4JGk\nN6W/OCJ+VDLjJNIbwj/L8/ciTP4t0lMkN5Gu8jxOevP26zNnzCPdL3Is6RvDhojo5UQbSW8BriNd\nhdpH+ka8NCL+VnweXya9X+E48NmIuDNzxjmkB2iMkb5eqyNicy9jacg7BtgZETOLj7OMI6dh9LjI\ncZfbZ9Smy3XscZHpLuMeT5Mx8B53meMut88b+R6Du9xnRiW6XKced5Ez8l0eyRNtMzMzMzMzs6oa\nxXu0zczMzMzMzCrLJ9pmZmZmZmZmGflE28zMzMzMzCwjn2ibmZmZmZmZZeQTbTMzMzMzM7OMfKJt\nZmZmZmZmlpFPtM3MzMzMzMwy8om2mZmZmZmZWUY+0TYzMzMzMzPLyCfaZmZmZmZmZhn5RNvMzMzM\nzMwsI59om5mZmZmZmWX0PzyxA8oLeg6fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5ebbc6e050>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(ncols=4, figsize=(12,3))\n",
    "\n",
    "df_samps['G1'].hist(ax=ax[0],bins=50,normed=True); ax[0].set_title('G1')\n",
    "df_samps['G2'].hist(ax=ax[1],bins=50,normed=True); ax[1].set_title('G2')\n",
    "df_samps['G1G2'].hist(ax=ax[2],bins=50,normed=True); ax[2].set_title('G1G2')\n",
    "df_mixsamps['G1G2_mix'].hist(ax=ax[3],bins=50,normed=True); ax[3].set_title('G1G2_mix')\n",
    "\n",
    "for a in ax.ravel(): a.set_xlim([0,40])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***The sympy mixture distribution (G1G2) doesn't look right...***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom mixture distribution class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This seems like it does the trick. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MixtureOfNormalDistributions(SingleContinuousDistribution):\n",
    "    _argnames = ('mean1', 'std1', 'mean2', 'std2', 'frac1')\n",
    "\n",
    "    @staticmethod\n",
    "    def check(mean1,std1,mean2,std2,frac1):\n",
    "        _value_check(std1 > 0, \"Standard deviation must be positive\")\n",
    "        _value_check(std2 > 0, \"Standard deviation must be positive\")\n",
    "\n",
    "    def pdf(self, x):\n",
    "        pdf1 = exp(-(x - self.mean1)**2 / (2*self.std1**2)) / (sqrt(2*pi)*self.std1)\n",
    "        pdf2 = exp(-(x - self.mean2)**2 / (2*self.std2**2)) / (sqrt(2*pi)*self.std2)\n",
    "        pdf = self.frac1*pdf1+(1-self.frac1)*pdf2\n",
    "        return pdf\n",
    "        #return pdf1 \n",
    "    \n",
    "    def sample(self):\n",
    "        if np.random.rand()<=self.frac1:\n",
    "            return random.normalvariate(self.mean1, self.std1)\n",
    "        else:\n",
    "            return random.normalvariate(self.mean2, self.std2)\n",
    "\n",
    "    #def _cdf(self, x):\n",
    "    #    mean, std = self.mean, self.std\n",
    "    #    return erf(sqrt(2)*(-mean + x)/(2*std))/2 + S.Half\n",
    "\n",
    "    #def _characteristic_function(self, t):\n",
    "    #    mean, std = self.mean, self.std\n",
    "    #    return exp(I*mean*t - std**2*t**2/2)\n",
    "\n",
    "def MixtureOfNormals(name, mean1, std1,mean2,std2,frac1):\n",
    " \n",
    "    return rv(name,  MixtureOfNormalDistributions, (mean1,std1,mean2,std2,frac1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1 = 5.\n",
    "std1 = 1.\n",
    "\n",
    "mu2 = 20.\n",
    "std2 = 1.\n",
    "\n",
    "frac1 = 0.8\n",
    "\n",
    "MoN = MixtureOfNormals('thing', mu1,std1,mu2,std2,frac1)\n",
    "\n",
    "MoN_samp = pd.DataFrame(np.array(list(sample_iter(MoN,numsamples=1000))),\n",
    "                        columns=['MixtureOfNormals']).astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAF0CAYAAACg3QoAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XucXHV9//HXZpNogCyEKASVqCh+qKlCgOCtUov6U36l\nv2q16g+8X2gfUkVRUFDrFS8g8VLR8kMqVMFqW+/gFbzWimBQYSkfEMSg3GIS2CSuLtnd3x9nVsfJ\nXs4kZ87szr6ejwePMOd7Zs5nvjO7897v+X7P9I2PjyNJkjSTBd0uQJIkzQ2GBkmSVIqhQZIklWJo\nkCRJpRgaJElSKYYGSZJUiqFBkiSVYmiQJEmlGBokSVIpC7tdgNQLIuJ84PnA+Zn54in2+Rzwf4C3\nZObbImIMeHdmnlZfpe2LiHsBLweeAwRwL+AW4MvAGZn5q5b9nw58ANgH+AvgYcDHgMsy80mTPP5b\ngBdk5oM7+DR2WUS8gOJ5HJSZ13e7HqkbHGmQqjEObAOeGRG7tTZGxHLg6MY+E1YA7yh7gIh4UUR8\nc1cLbUdE7Al8FzgJ+ChwBEVwOAlYA/w0Io5ouds7gE2N/X7UtP3IiHjGJIcZb/w3F8yVOqWOcKRB\nqs5VwCrgWcD5LW3HAjcAvw8UmXlnm4//GOr/0DobeDBwSMuIwvqI+ArwDeDfI+Jhmfm7Rtsy4CuZ\n+QuAiICi7nOAsyLi4sz8bW3PQFJlDA1SdbYDFwMvYsfQ8ALg8xThAYDm0xMR8UWKv8z/NDNHGu0v\npPjr/gnA24E/b2wfbRzjF8A3gadm5teaHvdbwOLMfGzTcU4F/gx4IrAmMwcj4kDg3RSjB8uBQeCt\nmfmlxv3uT3FK4vWtpyAAMnN7RJwEXAk8OyK+DfycIiC8sDGc/xdNd3nTxOMBb5muIyPiGOANwMHA\nGLCuUduljfY/bzz3ZwH/COyTmSsap4kOA04GzgIOAP6n0V97A2uBA4FrgRdn5jWNx+sH3gwcB+wP\nbAS+B7x2IvxMUuNK4AyK12Uv4FfAx4G3Z+bYdM9Pmqs8PSFV6yLg8RHxkIkNEfGnwGrgk9Pc72UU\nf6G/oXGf5RQfSGdl5veAv6H44Pw+xWmNTzXuN9nIw2TbXkJxmuEg4PqI2Ltx+0EUH7yrG7c/1/hA\nhuLDsA+4ZKqiM3MdcDtwFLC+UduvG/WtaNQ7se9dFMHh5Ih44FSPGRFPoghYVwGHA48CbgW+HBGH\ntOx+GkWfHdr03O8DvIIioBwJ3Jfiw/yNwIspQth+FPMumh/nFOC1FCMrxwAPBP5jqjqBCxvHegpF\nEHkdcCLwmmnuI81pjjRI1fo6cCfFX7ZvbGx7IXBNZl7TGKrfQWbeHhH/AJwfEZ+k+Ev5DooPWTJz\nc0TcA4xk5gb4/bB/X8m67srMMyZuRMTLKD7wHp2ZNzc2nxQRT6AYlfg28IDG9vUzPPbNwP6ZOQ7c\n2RjZGG6pc8I5wPEUf/FPNr8Biud+bWa+vKne5wFPopiQeXzTvl/PzC+23H8f4JWZeWPjvp8FTgAe\nl5lXNW17btN9zgb+LTNvaNz+VUScB3w4IpZn5sZJ6jyUYlLrTxu3fxkRg/zxvBWppxgapApl5mhE\n/Bvw/Ih4E8WH+rHA+0rc91MR8TTgcxR/5T5m4lRFBa5suX0EcGNTYJjwTeB5jf+fGLHon+Gx+yhO\nIcwoM8cj4pXAdyLiSZn5jUl2Oxz4dMv97omIK/nDiMKEH7GjrROBoWFT49+ftGzbs+n27yhes78G\n7g8s5g+/H5dTnK5o9XngLY3TOJcA38nM6ybZT+oZhgapep+gGB5/MsXP2D5Mf2qi2Ycp/sr/fmb+\nuMKa7mq5PQA8JCK2tGxfCCyKiIUUcyagmBdw1TSPfQDwhbKFZOb3GsHqAxHxyEl2GQDunmT7EMWp\ng2atzwt2/Et/vHHc4dZtTS6ieL1OAb4F/IZiJOTdkzz+hOcBf08xD+IEYKQxSnRSZg5Ncz9pzjI0\nSBXLzCsjIoH/CywCvpeZv5zpfo3JeO8DvkYxL+J5mfnxae4y8cHXeopiD+CeGQ63GbgReOok95+Y\n5PgtYBR4GlOEhog4nGLOwFdnOF6rkykmKJ44Sdtd/PEowIQ9mTwk7JKIWEoxh+Fdmfmhpu3T/n7M\nzFGK0xpnR8ReFPNOzqQYmXlR1XVKs4ETIaXOuBD4XxTXZriw5H3eRDFz/1iKax18sDH03az5A37i\nA/S+ExsaH14HlTjWDxrH2pKZN038RxES7oTfLwn9V+DExkqLP9L4UD0TuB74bIlj/l5jNca7KVY+\nrGhpvpxipUfzse5Ncdrih02bq1p+uoiiX3/ddLwFFCMIMEmoiohlEXFcYz8y867M/BeKCZerK6pL\nmnUcaZA64xPA2yj+4v/3mXaOiNUUExCfn5mbIuIM4NkUSzef3NhtM3BwRBxGMUnyZ41tJ0TEVRQf\nfm+nWM0wk49RXKDpPyLiNIrlgo8C/onig++1jf1Oorj2xHcj4h0UEz1/CzyCYunkAcCTG391t+u9\nFKsZXkSxOmLCGcClEXE28EFgCcUSzXs16ptQdhLotBr9fQPFMtFvUPxefAfwHeBPgD+PiIklp31N\n//4z8BcR8UGKORIHUVzxs/SpGmmucaRB6oDG2v7vApc0lhpOaL764TgwHhGLgAuAr2Xmpxr3HwVe\nSvGBNbGKYC3Fh9V3gWdm5m8o/hoeAK6gWOZ4UeP/m+1wxcXM3Aw8nmJU4YsUowWnN45xctN+QxTL\nFt/ZONblFKcVzqRYTnlwZg7OdLwp+miEIpQsat4/M79D8eF7GMUy0+9QXBTrCS2Xb57qGGWXoTZv\nOxYYoXh+n6ZYavlKiuf4AeCZzffJzE0U17zYn2Ly6PUU81E+RTEvQupJfePjXhVVkiTNrO3TExHx\nFIq/ii7LzGNb2vYDPkKxnnoL8C+Z+Yam9tMpLriyF0WiPyEzf77z5UuSpLq0dXoiIk4G3k8xFDeZ\nzwI3USwxOxJ4YuNiMUTExBXajgZWUpyPbWvylCRJ6p52RxqGKS4K80GKSUm/17j07AHAn2Xmdoov\n53l00y7HA2snzkk2Jl9tiogjMrN5RrQkSZqF2hppyMwPZWbrxWAmPA74KfDOiNgQET9rfJnNxHKp\nh9O01jszt1IEizU7VbkkSapVlUsuHwA8luJyqvtTfCnMZyPiZxSzufsoloc120Rx/XtJkjTLVRka\n+oA7M3Nt4/ZXGl8K8yx2XALWfJ/SyzfGx8fH+/oqWZotSdJ8s8sfoFWGhtvZ8RKvN1PMgdhE8YU2\nraMKewMbyh6gr6+PoaFhRkf9qvo69PcvYGBgiX1eI/u8fvZ5/ezz+k30+a6qMjRcC7w2InZrXHQG\n4EHALzLzdxFxDcXFWr4Lv7/c7UMpll6WNjo6xvbtvsnqZJ/Xzz6vn31eP/t87qkyNHyRYs7CmY2l\nmY8G/pri+vtQXL/h9RHxFYpL1r4H+FFmrquwBkmS1CFthYaIGKaYg7CocfvpwHhm7paZv42IpwLn\nUHzxy53A32XmfwFk5jkRsYLia2f3oLj06jOqeiKSJKmz5tplpMc3b97mcFZNFi5cwLJlu2Of18c+\nr599Xj/7vH6NPt/liZB+YZUkSSrF0CBJkkoxNEiSpFIMDZIkqRRDgyRJKsXQIEmSSjE0SJKkUgwN\nkiSpFEODJEkqpcrvnlAPGBkZYXDwamDyb6JbteoRLF68uJslSpK6xNCgPzI4eDWnrP0MS5ev3KFt\ny8b1nHESrF59WBcqkyR1m6FBO1i6fCV7rTiw22VIkmYZ5zRIkqRSDA2SJKkUQ4MkSSrF0CBJkkox\nNEiSpFIMDZIkqRRDgyRJKsXQIEmSSjE0SJKkUgwNkiSpFEODJEkqxdAgSZJKMTRIkqRSDA2SJKkU\nQ4MkSSrF0CBJkkpZ2O4dIuIpwAXAZZl57BT77A5cB3w9M1/ctP104DnAXsDlwAmZ+fOdKVySJNWr\nrZGGiDgZeD9w/Qy7vg3YveW+r6AIDEcDK4GfAZ9t5/iSJKl72j09MQwcAdw41Q4R8UiKcHB+S9Px\nwNrMvD4ztwGnAQ+PiCParEGSJHVBW6EhMz+UmVtm2O0jFIHg7okNEXFv4OHAVU2PtRW4AVjTTg2S\nJKk72p7TMJ2I+DtgNDMviIg3NzUtA/qAzS132QTcp51j9Pc7d7OTZurf/v4FLFzoa9ApE/3v+7w+\n9nn97PP6VdXXlYWGiNgHeCtwVBt36wPG2znOwMCSdnZXm2bq34GBJSxbtvu0+2jX+T6vn31eP/t8\n7qlypOEs4ILMvHaStk3AGDuOKuwNbGjnIENDw4yOju1chZrR0NDwjO2bN2+rqZr5p79/AQMDS3yf\n18g+r599Xr+JPt9VVYaG44DNETGxxHI3YEFEHJOZ+0TENcBhwHcBImIv4KEUSy9LGx0dY/t232Sd\nMtMPsP1fD/u5fvZ5/ezzuafK0PCAltuvAe4PvLpx+yPA6yPiK8CvgPcAP8rMdRXWIEmSOqSt0BAR\nwxRzEBY1bj8dGM/M3TLz1pZ9h4BlmXkbQGaeExErgG8BewDfBJ6xy89AkiTVoq3QkJmlT4hk5lun\n2LbDdkmSNPu53kWSJJViaJAkSaUYGiRJUimGBkmSVIqhQZIklWJokCRJpRgaJElSKYYGSZJUiqFB\nkiSVYmiQJEmlGBokSVIphgZJklSKoUGSJJViaJAkSaUYGiRJUimGBkmSVIqhQZIklWJokCRJpRga\nJElSKYYGSZJUysJuF6C5Y2x0O5nXTbvPqlWPYPHixTVVJEmqk6FBpW276zbOu/hWlv5g66TtWzau\n54yTYPXqw2quTJJUB0OD2rJ0+Ur2WnFgt8uQJHWBcxokSVIphgZJklSKoUGSJJViaJAkSaW0PREy\nIp4CXABclpnHtrT9DfCPwEOAXwJnZeZHm9pfCbwc2Bf4KfDqzFy38+VLkqS6tDXSEBEnA+8Hrp+k\nbQ3wCeCNwJ7AScDZEfHYRvtfAW8GngusAC4GvhQRS3blCUiSpHq0e3piGDgCuHGStr2B0zPzS5k5\nlplfphhNOLLRfjzwscy8MjN/B5wJjAN/tXOlS5KkOrUVGjLzQ5m5ZYq2r2bm6RO3I6If2I/iNAXA\nYcC6pv3HgR8Da9otWpIk1a+TEyHPALYCn27cXg5sbtlnE3CfDtYgSZIq0pErQkbEe4BnA0/IzJFp\ndu2jOEVRWn+/Cz46aVf7t79/AQsX+hrtrIn+931eH/u8fvZ5/arq60pDQ0T0AecDhwOPzcz1Tc0b\n2HFUYW/g6naOMTDgvMlO2tX+HRhYwrJlu1dUzfzl+7x+9nn97PO5p+qRhg8Af0IRGO5uabuSYl7D\nxwEiYgFwKPBR2jA0NMzo6FgFpWoyQ0PDu3z/zZu3VVTN/NPfv4CBgSW+z2tkn9fPPq/fRJ/vqspC\nQ0Q8DjgOOGiSwADwEeCTEfFJilUVJwO/pVh6Wdro6Bjbt/sm65Rd/QH29amG/Vg/+7x+9vnc01Zo\niIhhijkIixq3nw6MZ+ZuwIuAAeAXEdF8t+9k5lMz86sRcSrFxMj7AlcA/7ux/FKSJM1ybYWGzJxy\nbCMzXwq8dIb7nwOc084xJUnS7NCR1ROavUZGRhgcnHruaeZ1NVYjSZpLDA3zzODg1Zyy9jMsXb5y\n0vY7brqCfQ/weluSpB0ZGuahpctXsteKAydt27LxlpqrkSTNFV5ZQ5IklWJokCRJpRgaJElSKYYG\nSZJUiqFBkiSVYmiQJEmlGBokSVIphgZJklSKoUGSJJViaJAkSaUYGiRJUimGBkmSVIqhQZIklWJo\nkCRJpRgaJElSKYYGSZJUiqFBkiSVYmiQJEmlGBokSVIphgZJklTKwm4XoGqNjIwwOHj1lO2Z19VY\njSSplxgaeszg4NWcsvYzLF2+ctL2O266gn0PWFNzVZKkXmBo6EFLl69krxUHTtq2ZeMtNVcjSeoV\nzmmQJEmlGBokSVIpbZ+eiIinABcAl2XmsS1tzwZOAx4MJHBaZn69qf104DnAXsDlwAmZ+fOdL1+S\nJNWlrZGGiDgZeD9w/SRthwDnA6cA9wHeB3w2Iu7XaH8FRWA4GlgJ/Az47C7ULkmSatTu6Ylh4Ajg\nxknaXgJcnJlfzcyRzLwIuBp4bqP9eGBtZl6fmdsoRiQeHhFH7GTtkiSpRm2Fhsz8UGZumaL5MGBd\ny7Z1wJqIuDfwcOCqpsfaCtwAuP5PkqQ5oMqJkMuBzS3bNlGcqlgG9E3TLkmSZrlOX6ehDxjfhfYd\n9Pe74GM63e6f/v4FLFzoa7SzJl6/br+O84l9Xj/7vH5V9XWVoWEDO44a7N3YvgkYm6a9tIGBJTtb\n37zQ7f4ZGFjCsmW7d7WGXtDt13E+ss/rZ5/PPVWGhisp5jU0WwNclJm/i4hrGu3fBYiIvYCHUiy9\nLG1oaJjR0bEKyu1NQ0PDXT/+5s3bulrDXNbfv4CBgSW+z2tkn9fPPq/fRJ/vqipDw7nADyPiaOAy\n4DjgQODCRvtHgNdHxFeAXwHvAX6Uma2TJ6c1OjrG9u2+yabS7R9AX59q2I/1s8/rZ5/PPW2FhogY\nppiDsKhx++nAeGbulpmDEXEcxXUcVgLXAn+ZmXcCZOY5EbEC+BawB/BN4BlVPRFJktRZbYWGzJx2\nbCMzPwd8bpr2twJvbeeYkiRpdnDqqiRJKsXQIEmSSjE0SJKkUgwNkiSpFEODJEkqxdAgSZJKMTRI\nkqRSDA2SJKkUQ4MkSSrF0CBJkkoxNEiSpFIMDZIkqRRDgyRJKsXQIEmSSjE0SJKkUgwNkiSpFEOD\nJEkqxdAgSZJKMTRIkqRSDA2SJKkUQ4MkSSrF0CBJkkoxNEiSpFIMDZIkqRRDgyRJKmVhtwtQ7xgb\n3U7mddPus2rVI1i8eHFNFUmSqmRoUGW23XUb5118K0t/sHXS9i0b13PGSbB69WE1VyZJqoKhQZVa\nunwle604sNtlSJI6oNLQEBEHA2uBQ4Fh4FLgVZm5MSKOAt4FHASsB96VmRdVeXxJktQ5lU2EjIgF\nwCXA94H7AquAfYAPR8QK4PPAhxttrwLOjYhDqzq+JEnqrCpXT9wP2A/4RGZuz8zNwGeA1cBxQGbm\nBZk5kpmXAl8AXlrh8SVJUgdVGRp+BVwFHB8Ru0fEPsAzgS8BhwHrWvZfB6yp8PiSJKmDKpvTkJnj\nEfFM4BsUpx8AvgWcRnFq4paWu2wC7tPucfr7vbTEdGZ7//T3L2DhwtldYzdNvH6z/XXsJfZ5/ezz\n+lXV15WFhohYDHwR+BTwTmAPijkMF05xlz5gvN3jDAws2dkS54XZ3j8DA0tYtmz3bpcx683217EX\n2ef1s8/nnipXTzwReFBmnta4vTUi3gL8GPgyO44q7A1saPcgQ0PDjI6O7UqdPW1oaLjbJUxraGiY\nzZu3dbuMWau/fwEDA0t8n9fIPq+ffV6/iT7fVVWGhn5gQUQsyMyJd8G9KUYTvgG8sGX/NcDl7R5k\ndHSM7dt9k01ltv8A+vqVYz/Vzz6vn30+91QZGr4PbAXeGhHvBHajmM/wbeDjwJsj4sUUpyueCBwN\nPKrC40uSpA6qbBZKZm4CngI8DvglcDXwG+DYzPw1cAzwCuAu4CzguMwcrOr4kiSpsyq9ImRmXgUc\nNUXb9yiu2SBJkuYg17tIkqRSDA2SJKkUQ4MkSSrF0CBJkkoxNEiSpFIMDZIkqRRDgyRJKsXQIEmS\nSjE0SJKkUgwNkiSpFEODJEkqxdAgSZJKMTRIkqRSDA2SJKkUQ4MkSSrF0CBJkkoxNEiSpFIMDZIk\nqRRDgyRJKsXQIEmSSjE0SJKkUgwNkiSpFEODJEkqxdAgSZJKWdjtAtS+kZERBgevnrQt87qaq5Ek\nzReGhjlocPBqTln7GZYuX7lD2x03XcG+B6zpQlWSpF5naJijli5fyV4rDtxh+5aNt3ShGknSfOCc\nBkmSVEpHRhoi4g3ACcBS4L+Bl2XmLyLiKOBdwEHAeuBdmXlRJ2qQJEnVqnykISJOAI4FjgT2A64F\nXh0RK4DPAx8G7gu8Cjg3Ig6tugZJklS9Tow0nASclJk/a9x+FUBEvAbIzLygsf3SiPgC8FLg5R2o\nQ5IkVajS0BAR9wMeDCyPiEFgX+AyilBwGLCu5S7rgGdVWYMkSeqMqkcaHtD495nAUUA/8J/AucBu\nQOvU/k3Afdo5QH+/czfnch/09y9g4cK5W3+nTby2c/k1nmvs8/rZ5/Wrqq+rDg19jX/fk5l3AETE\nm4EvA1+fYv/xdg4wMLBklwrsBXO5DwYGlrBs2e7dLmPWm8uv8Vxln9fPPp97qg4Ntzf+vbtp280U\n4WARO44q7A1saOcAQ0PDjI6O7Wx9PWFoaLjbJey0oaFhNm/e1u0yZq3+/gUMDCzxfV4j+7x+9nn9\nJvp8V1UdGn4JDAGHAD9ubHswMAJcAjy/Zf81wOXtHGB0dIzt2+f3m2wu/5D5+pVjP9XPPq+ffT73\nVBoaMnM0Is4D3hAR3wW2AG8CPg78K/CmiHgxcCHwROBo4FFV1iBJkjqjE7NQTgW+AvwQuAFI4MTM\n3AAcA7wCuAs4CzguMwc7UIMkSapY5ddpyMwRimDwiknavgesrvqYkiSp81zvIkmSSjE0SJKkUgwN\nkiSpFEODJEkqxdAgSZJKMTRIkqRSDA2SJKkUQ4MkSSrF0CBJkkoxNEiSpFIMDZIkqRRDgyRJKsXQ\nIEmSSjE0SJKkUgwNkiSpFEODJEkqZWG3C9D8MTa6nczrpt1n1apHsHjx4poqkiS1w9Cg2my76zbO\nu/hWlv5g66TtWzau54yTYPXqw2quTJJUhqFBtVq6fCV7rTiw22VIknaCcxokSVIphgZJklSKoUGS\nJJViaJAkSaUYGiRJUimGBkmSVIqhQZIklWJokCRJpXTs4k4R8T7gxMxc0Lh9FPAu4CBgPfCuzLyo\nU8eXJEnV6shIQ0QcAjwPGG/c3g/4PPBh4L7Aq4BzI+LQThxfkiRVr/LQEBF9wEeAs5o2HwdkZl6Q\nmSOZeSnwBeClVR9fkiR1RidGGv4eGAaaTz0cCqxr2W8dsKYDx5ckSR1Q6ZyGiNgXeAtwZEvTcuCW\nlm2bgPtUeXxJktQ5VU+EPAs4LzMzIh44w759NOY8tKO/3wUfvdwH/f0LWLiwd5/fTCZe215+jWcb\n+7x+9nn9qurrykJDRDwReCzwssamvqbmDew4qrB3Y3tbBgaW7FR9vaSX+2BgYAnLlu3e7TK6rpdf\n49nKPq+ffT73VDnScBywD7A+IqCYL9EXEXdSjEAc27L/GuDydg8yNDTM6OjYLpY6u42MjHDNNVdP\n2X7ddf9TYzX1GhoaZvPmbd0uo2v6+xcwMLBkXrzPZwv7vH72ef0m+nxXVRkaXg28sen2/sB/Awc3\njnNqRLwYuBB4InA08Kh2DzI6Osb27b39JvvJT37CKWs/w9LlKydtv+OmK9j3gN6cQzofXt8y7If6\n2ef1s8/nnspCQ2beDdw9cTsiFgHjmXlb4/YxwD8BZwM3A8dl5mBVx+81S5evZK8VB07atmVj65xS\nSZI6r2NXhMzMXwD9Tbe/B6zu1PEkSVJnOXVVkiSVYmiQJEmlGBokSVIphgZJklSKoUGSJJViaJAk\nSaUYGiRJUikdu06DJEkjIyMMDv7xZfGbLyN90EGrWLx4cZeqU7sMDZKkjhkcvHrKy+Jv2bieM04a\nY/Xqw7pQmXaGoUGS1FHTXRZfc4tzGiRJUimONEiSdtpkcxaaZV5XYzXqNEODJGmnTTdnAeCOm65g\n3wPW1FyVOsXQIEnaJdPNWdiy8Zaaq1EnOadBkiSVYmiQJEmlGBokSVIphgZJklSKoUGSJJXi6gnN\nGmOj26dd071q1SO8Rr0kdZGhQbPGtrtu47yLb2XpD7bu0FZcox6vUS9JXWRo0KziNeolafZyToMk\nSSrFkQZJUlfMNI/pnnvuAWDRokVT7uNcp3oZGiRJXTHdPCYovrditz33nfJ7LZzrVD9DgySpa2b6\n3oqly/d3ntMs4pwGSZJUiqFBkiSVUunpiYhYCbwfOBIYAb4KnJiZQxFxSKPtEOAO4JzMXFvl8SVJ\nUudUPdLwRWATsD9wOLAKeG9E3LvR9g1gP+A5wKkR8bSKjy9JkjqkstAQEXsCVwCnZuZwZt4KXEAx\n6vCXwCLg9EbbVcBHgeOrOr4kSeqsyk5PZObdwEtbNu8P/Ao4DPhpZo43ta2bZH9JkjRLdWwiZEQc\nDvwDcDqwHNjcsssmYO9OHV+SJFWrI9dpiIjHAV8AXpeZl0XEsyfZrQ8Yn2T7tPr7e3/Bx3x4jjuj\nv38BCxf2dt9MvPa+B+pjn09vZGSEa665esr2G27IGqvZ0Xz4vVCFqt7flYeGiDgG+ARwQmZe2Ni8\nAXhoy657AxvbffyBgSW7VuAcMB+e484YGFjCsmW7d7uMWvgeqJ99PrkrrriW15z5H1NelfGOm65g\n3wPW1FxVYWx0O7/85c+nfO0OPvhgLzFdsaqXXD6WYvLjMzLz0qamK4G/j4gFmTnW2LYGuLzdYwwN\nDTM6OjbzjnPY0NBwt0uYlYaGhtm8eVu3y+io/v4FDAwsmRfv89nCPp/e0NDwjFdt7JZtd93G+z55\nK0uXb9ihbcvG9Zx18jCHHuolpuEP7/NdVVloiIh+4FyKUxKXtjRfAgwBb4yIM4FHAi8Bjm33OKOj\nY2zf3ts/2P7imtx8eO0nzKfnOlvY55Ob7b+Ppgs0vqbVq3Kk4THAQcAHI+KfKOYrTMxbCOAY4Bzg\nVOB24PWZ+ZUKjy9JkjqoyiWX3wP6Z9jt8VUdT5Ik1cspp5IkqRS/GltzwtjodjKvm3afVase4Uxp\nSeogQ4PmhG133cZ5F9/K0h9snbR9y8b1nHESrF7tTGlJ6hRDg+aM6WZJS5I6zzkNkiSpFEODJEkq\nxdAgSZJKcU6DJM1jIyMjDA5O/YVUM61a0vxiaJCkeWxw8GpOWfuZWfmFVJp9DA1dYLKXNJvM1i+k\n0uxjaOgCk70kaS4yNHSJyV6SNNe4ekKSJJXiSIMkzXEzzZMCv5tF1TA0SNIcN9M8Kb+bRVUxNHTI\ndMnf1RHl7ChAAAAHUklEQVSSquZ3s6gOhoYOmS75uzpCkjQXGRo6aKrk7+oISdJc5OoJSZJUiqFB\nkiSVYmiQJEmlOKdBknrc2Oj2KVdtuZpL7TA0SFKP23bXbZx38a0s/cHWHdpczaV2GBokaR5wNZeq\n4JwGSZJUiiMN6gnTnbOd4LX3pfljpt8J99xzDwCLFi2ach9/Z+zI0KCeMN05W4C7N/ycl/3VdUQc\nNGm7vxyk3jLT74Q7brqC3fbc1+/raJOhQT1jumvvb9l4C+ddfO2kv0D85SD1ppl+Jyxdvr/f19Gm\nWkNDRDwQOBt4NLAF+FRmvr7OGjR/+YU+ksqa6fTGfB2drHuk4T+BK4DnAPsCl0TE7Zn5/prrkCRp\nStOd3pjPo5O1hYaIOBx4JHBUZm4FtkbEWuBEwNCgrnESpaTJODq5ozpHGg4Fbs7MoaZt64CIiD0a\nQWLOuPbaQd783nPZY2DZpO133XY9ffsdWXNV2hkzTZiaz39VSFKzOkPDcmBzy7ZNTW2lQkN//+y4\ntMTw8DbG9lxF374PmbR9/Ne3sXXj+knbfnP37cD4lI+9K+1z9bG7eezf3H07u+2575T3Bbjhhqzl\nvbdgQR977HFvtm79LWNjUz8fVacX+vyGG5ItU/y+gfn7c92px96ycT39/UewcOHs+Dwqo6rfX33j\n4/X8kETEqcDTMvNRTdseCiTw4Myc+h0vSZK6rs6YtAG4T8u2vSmi3K9rrEOSJO2EOkPDlcADI2Lv\npm1HANdm5m9qrEOSJO2E2k5PAETE94FrgNcA9wcuBs7MzH+urQhJkrRT6p7F8UyKsHA7cBlwvoFB\nkqS5odaRBkmSNHfNnfUikiSpqwwNkiSpFEODJEkqxdAgSZJKMTRIkqRSDA2SJKmUOr+waqdFxAOB\ns4FHA1uAT2Xm67tbVW+LiDHgdxSX+e5r/HtuZp7Y1cJ6SEQ8BbgAuCwzj21pezZwGvBgiu9nOS0z\nv15/lb1lqj6PiBcA/0Lxnoc/vOePzMwray+0h0TESuD9wJHACPBV4MTMHIqIQxpthwB3AOdk5tqu\nFdsjpupzYBnwc+C3jV0n3udvLNvvcyI0AP8JXAE8B9gXuCQibs/M93e3rJ42DjwsM2/pdiG9KCJO\nBl4MXD9J2yHA+cDTgG9SXBTtsxHxsMy8tc46e8l0fd7w7cw8qsaS5osvUvz+3p/iQ+tzwHsj4pWN\ntnOAo4GDgK9FxE2Z+bluFdsjJu1z4HRgPDN329kHnvWnJyLicOCRwOsyc2tm3gisBY7vbmU9r6/x\nnzpjmOK7V26cpO0lwMWZ+dXMHMnMi4CrgefWWWAPmq7P1QERsSfFh9epmTncCL0XUPwF/JfAIuD0\nRttVwEfxd/sumaHPd9lcGGk4FLg5M4eatq0DIiL2yMytXaprPnhPRDwWWAr8O3BSZm7rck09ITM/\nBBARkzUfBnypZds6YE2Hy+ppM/Q5wP4R8TXgcGAT8ObMvLCm8npSZt4NvLRl8/7Aryje5z/NzObL\nEq+bZH+1YYo+X0nR5wB9EXEB8GSgHzgPeFNmjpZ5/Fk/0gAsBza3bNvU1KbO+G/ga8BDgcdQzCc5\nu6sVzR9Tvedbv1pe1dlAcdritRSnQN8AfCwintDNonpNY+T4HyiGyad6n+/dej/tvKY+fwfFnJ3/\nojjlvz/FaM9zgTeVfby5MNIwmYlhc784o0My83HNNyPidcAXIuJlmXlPt+qaxyYmLKkDMvMS4JKm\nTZ+KiKcDLwK+1ZWiekxEPA74AsWp5ssak31b+T6vUFOfn5KZ32xsfnzTLldGxDuBU4G3lHnMuTDS\nsIEd/8Lam+KN9ev6y5m3bqYYytqny3XMB1O95zd0oZb57Gbgft0uohdExDHAxcArM3NixHKq9/nG\nOmvrVVP0+WRuBlaUfdy5EBquBB4YEc1DVkcA12bmb7pUU0+LiEMi4r0tmx9OMbTl7P3Ou5LifG+z\nNcDlXahlXoiIv4uIv23Z/CfATd2op5c05kVdADyjZY7IlcDBEdH8OeT7vAJT9XlEHBURp7Xs/nCK\n4FDKrD89kZk/jogfAu+OiNcA9wdeDZzZ3cp62p3A8RFxJ8Va3wcBb6NYQ+3QYeedC/wwIo4GLgOO\nAw4EPtHVqnrbvYAPRsRNwE+Av6VYBnhEV6ua4yKin+L9/LrMvLSl+RJgCHhjRJxJsUruJcCxaKfN\n0OebgX+MiJuBT1NcH+M1wBllH79vfHz2fwZExP0oOuEJwN3ARzLz7V0tqsdFxJ9RvJH+lOJCIOcD\nb3A+QzUiYpjiFNuixqbtNK2fjoinAe+hmPV8LcUQ4391o9ZeUaLPT6OYdb6C4gI4r83ML3ej1l7R\n+D3ybYpRyon5ChP/BsXKrHMoVqzcDrwrM/9fd6rtDSX6/FCK+QsPowgRH8zM3goNkiSp++bCnAZJ\nkjQLGBokSVIphgZJklSKoUGSJJViaJAkSaUYGiRJUimGBkmSVIqhQZIklWJokCRJpRgaJElSKYYG\nSZJUyv8H41R0nXe2Y7YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5ef0470c90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "MoN_samp.hist(bins=50);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.92922059961305e-52*sqrt(2)*(8.60712277201052e+51*pi + sqrt(pi)*(112.5 + 8.60712277201052e+51*sqrt(pi)))/pi + 5.54537317456468e-50*sqrt(2)*(-0.561714787011297*sqrt(pi) + 1.0)/sqrt(pi) + 14.1421356237309*sqrt(2)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "E(MoN)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:_jupyter]",
   "language": "python",
   "name": "conda-env-_jupyter-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}