Clutter_Problem.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets try to implement the clutter problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import simps\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets initiliaze the hyper parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 355,
   "metadata": {},
   "outputs": [],
   "source": [
    "w=0.5;\n",
    "clutter_mean=0;\n",
    "clutter_var=130;\n",
    "prior_mean=0;\n",
    "prior_var=30;\n",
    "var=1;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets generate now the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 356,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the value of the true mean is:  2\n"
     ]
    }
   ],
   "source": [
    "n_samples=20\n",
    "x=np.zeros(n_samples)\n",
    "#mean=np.random.normal(loc=prior_mean, scale=prior_var)\n",
    "mean=2\n",
    "print(\"the value of the true mean is: \",mean)\n",
    "for i in range(n_samples):\n",
    "    coin=np.random.random()\n",
    "    if coin<w:\n",
    "        x[i]=np.random.normal(loc=clutter_mean, scale=clutter_var)\n",
    "    else:\n",
    "        x[i]=np.random.normal(loc=mean, scale=var)\n",
    "x;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 357,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculating_cavities(m_n, v_n, m, v):\n",
    "    ## m, v are the current estimates of the posterior, \n",
    "    ## m_n, v_n are the current ones we wanna update\n",
    "    ##this creates an indecision form, i have to right it such that it gives no problem \n",
    "    #cavity_v = v*v_n/(v_n-v)\n",
    "    aux = 1/v - 1/v_n\n",
    "    cavity_v = 1/aux\n",
    "    cavity_m = m + cavity_v*(m-m_n)/v_n\n",
    "    return cavity_m, cavity_v\n",
    "\n",
    "##################################################################################################\n",
    "\n",
    "def normal_pdf(x, mean, var):\n",
    "    if var>0:\n",
    "        norm=(2*np.pi*var)**(1/2)\n",
    "    else:\n",
    "        #norm=-1*(-2*np.pi*var)**(1/2)\n",
    "        #print(\"i have a negative variance!\\n\")\n",
    "        ##it's gonna be negative....not sure is gonna work out\n",
    "        norm=(2*np.pi)**(1/2)*(-1*np.sqrt(-var))\n",
    "    exp=np.exp((-0.5/var) * ((x-mean)**2))\n",
    "    return exp/norm\n",
    "\n",
    "##################################################################################################\n",
    "\n",
    "def normalizer(w, x, cavity_m, cavity_v, clutter_var):\n",
    "    z = (1-w)*normal_pdf(x, cavity_m, cavity_v+1)+w*normal_pdf(x, 0, clutter_var)\n",
    "    return z\n",
    "\n",
    "##################################################################################################\n",
    "\n",
    "def update(w, x, cavity_m, cavity_v, m_n, v_n, clutter_var):\n",
    "    z = normalizer(w, x, cavity_m, cavity_v, clutter_var)\n",
    "    ro = 1-w*normal_pdf(x, 0, clutter_var)/z\n",
    "    \n",
    "    update_m = cavity_m+ro*(cavity_v/(cavity_v+1))*(x-cavity_m)\n",
    "    update_v = cavity_v - ro*(cavity_v)**(2)/(cavity_v+1) + ro*(1-ro)*(cavity_v)**(2)*(x-cavity_m)**(2)/((cavity_v+1)**(2))              \n",
    "    return update_m, update_v\n",
    "\n",
    "##################################################################################################\n",
    "\n",
    "def new_marginal(update_m, update_v, cavity_m, cavity_v, z):\n",
    "    aux=1/update_v-1/cavity_v\n",
    "    if aux!=0:\n",
    "        new_v=1/aux\n",
    "    else:\n",
    "        ##this remains the same\n",
    "        new_v=np.inf\n",
    "    if new_v==np.inf and (1/cavity_v == 0 or update_m==cavity_m):\n",
    "        new_m=cavity_m\n",
    "    else:\n",
    "        new_m=cavity_m+(new_v+cavity_v)/(cavity_v)*(update_m-cavity_m)\n",
    "    \n",
    "    denominator=((2*np.pi*new_v)**(1/2)*normal_pdf(new_m, cavity_m, new_v+cavity_v));\n",
    "    if denominator>0:\n",
    "        new_s=z/denominator\n",
    "    else:\n",
    "        new_s=np.inf\n",
    "    return new_v, new_m, new_s\n",
    "\n",
    "##################################################################################################\n",
    "\n",
    "def convergence(m_old, v_old, s_old, m_new, v_new, s_new, thresh):\n",
    "    #diff_m=np.abs(m_old-m_new)\n",
    "    #diff_v=np.abs(v_old, v_new)\n",
    "    #diff_s=np.abs(s_old, s_new)\n",
    "    #val=np.max(np.concatenate([diff_m, diff_v, diff_s]))\n",
    "    val_old=np.concatenate([m_old, v_old, s_old])\n",
    "    val_new=np.concatenate([m_new, v_new, s_new])\n",
    "    val=0\n",
    "    for i in range(len(val_old)):\n",
    "        if val_old[i]!=np.inf or val_new[i]!=np.inf:\n",
    "            ##this way i have no form of indecision\n",
    "            aux=np.abs(val_old[i]-val_new[i])\n",
    "            if aux>val:\n",
    "                val=aux\n",
    "                \n",
    "    #print(\"val: \", val)\n",
    "    if val<thresh:\n",
    "        return (True, val)\n",
    "    else:\n",
    "        return (False, val)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 358,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.,  0.,  0.,  0.,  0., 12.,  3.,  0.,  2.,  2.]),\n",
       " array([-297.59417975, -246.23402598, -194.87387222, -143.51371845,\n",
       "         -92.15356469,  -40.79341092,   10.56674284,   61.92689661,\n",
       "         113.28705037,  164.64720414,  216.00735791]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 358,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAANB0lEQVR4nO3df6xk9VnH8fdHtlShNCxyWyk/XJoQImlMIDeK1lRTWrqlBGqiCcTqWkj2r1pqTOySTeQPYwLWNGo0kk3BYkQwARqoRGHFEmJS0IUiXbpQfhTL2pW9DaZUm5Ru+vjHPehwuXvv3Zm5P57L+5VM5pzvnJnzPLsznz17fsykqpAk9fMj612AJGk8BrgkNWWAS1JTBrgkNWWAS1JTW9ZyZaeeempt27ZtLVcpSe09+uij366qmYXjaxrg27ZtY9++fWu5SklqL8m/LzbuLhRJasoAl6SmDHBJasoAl6SmDHBJasoAl6Smlg3wJDcnOZxk/8jYZ5I8leSJJF9IcvLqlilJWmglW+CfB7YvGNsLvKeqfhr4OnDtlOuSJC1j2QCvqoeAlxeM3V9VR4bZh4EzVqE2SdISpnEl5lXA3x7twSQ7gZ0AZ5111hRWJ03ftl33rtu6X7j+I+u2bvU20UHMJLuBI8CtR1umqvZU1WxVzc7MvOFSfknSmMbeAk+yA7gUuKj8XTZJWnNjBXiS7cCngV+squ9NtyRJ0kqs5DTC24AvA+cmOZjkauDPgJOAvUkeT3LjKtcpSVpg2S3wqrpykeGbVqEWSdIx8EpMSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWpq2QBPcnOSw0n2j4ydkmRvkmeG+62rW6YkaaGVbIF/Hti+YGwX8EBVnQM8MMxLktbQsgFeVQ8BLy8Yvhy4ZZi+BfjolOuSJC1j3H3g76yqQwDD/TuOtmCSnUn2Jdk3Nzc35uokSQut+kHMqtpTVbNVNTszM7Paq5OkN41xA/ylJKcBDPeHp1eSJGklxg3we4Adw/QO4O7plCNJWqmVnEZ4G/Bl4NwkB5NcDVwPfDDJM8AHh3lJ0hrastwCVXXlUR66aMq1SJKOgVdiSlJTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNWWAS1JTBrgkNTVRgCf57SRPJtmf5LYkPzqtwiRJSxs7wJOcDnwSmK2q9wDHAVdMqzBJ0tIm3YWyBfixJFuAE4BvTV6SJGklxg7wqvoP4I+AbwKHgO9U1f0Ll0uyM8m+JPvm5ubGr1SS9DqT7ELZClwOnA28CzgxyccWLldVe6pqtqpmZ2Zmxq9UkvQ6k+xC+QDwjaqaq6ofAHcBPz+dsiRJy5kkwL8JXJjkhCQBLgIOTKcsSdJyJtkH/ghwB/AY8NXhtfZMqS5J0jK2TPLkqroOuG5KtUiSjoFXYkpSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUxMFeJKTk9yR5KkkB5L83LQKkyQtbcuEz/8T4B+q6leSHA+cMIWaJEkrMHaAJ3k78D7gNwGq6lXg1emUJUlaziS7UN4NzAF/meQrST6X5MSFCyXZmWRfkn1zc3MTrE6SNGqSAN8CXAD8RVWdD/wPsGvhQlW1p6pmq2p2ZmZmgtVJkkZNEuAHgYNV9cgwfwfzgS5JWgNjB3hV/SfwYpJzh6GLgK9NpSpJ0rImPQvlt4BbhzNQngc+PnlJkqSVmCjAq+pxYHZKtUiSjoFXYkpSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDVlgEtSUwa4JDU1cYAnOS7JV5L83TQKkiStzDS2wK8BDkzhdSRJx2CiAE9yBvAR4HPTKUeStFKTboH/MfC7wA+PtkCSnUn2Jdk3Nzc34eokSa8ZO8CTXAocrqpHl1quqvZU1WxVzc7MzIy7OknSApNsgb8XuCzJC8DtwPuT/PVUqpIkLWvsAK+qa6vqjKraBlwB/FNVfWxqlUmSluR54JLU1JZpvEhVPQg8OI3XkiStjFvgktSUAS5JTRngktSUAS5JTRngktSUAS5JTRngktSUAS5JTRngktSUAS5JTRngktTUVL4LRdL4tu26d13W+8L1H1mX9Wp63AKXpKYMcElqygCXpKYMcElqygCXpKYMcElqygCXpKYMcElqygCXpKYMcElqygCXpKYMcElqygCXpKbGDvAkZyb5UpIDSZ5Mcs00C5MkLW2Sr5M9AvxOVT2W5CTg0SR7q+prU6pNkrSEsbfAq+pQVT02TH8XOACcPq3CJElLm8oPOiTZBpwPPLLIYzuBnQBnnXXWNFYnqbn1+hGL9bQaP6Ax8UHMJG8D7gQ+VVWvLHy8qvZU1WxVzc7MzEy6OknSYKIAT/IW5sP71qq6azolSZJWYpKzUALcBByoqs9OryRJ0kpMsgX+XuDXgfcneXy4XTKluiRJyxj7IGZV/TOQKdYiSToGXokpSU0Z4JLUlAEuSU0Z4JLUlAEuSU0Z4JLUlAEuSU0Z4JLUlAEuSU0Z4JLUlAEuSU0Z4JLU1FR+kWctrOcveKzGL2lI0qTcApekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpiYK8CTbkzyd5Nkku6ZVlCRpeWMHeJLjgD8HPgycB1yZ5LxpFSZJWtokW+A/AzxbVc9X1avA7cDl0ylLkrScSX7Q4XTgxZH5g8DPLlwoyU5g5zD730meHmNdpwLfHuN5U5Eb1nR169rrGnsz9QobrN9Vfl9vqF5X2Yp6nfDP+ycXG5wkwLPIWL1hoGoPsGeC9ZBkX1XNTvIaXdjr5vVm6tde18Yku1AOAmeOzJ8BfGuyciRJKzVJgP8rcE6Ss5McD1wB3DOdsiRJyxl7F0pVHUnyCeA+4Djg5qp6cmqVvd5Eu2CasdfN683Ur72ugVS9Ybe1JKkBr8SUpKYMcElqakMFeJLfT/JEkseT3J/kXcN4kvzpcMn+E0kuGHnOjiTPDLcd61f9sUnymSRPDf18IcnJI49dO/T6dJIPjYy3/eqCJL+a5MkkP0wyu+CxTdfvqM3Sx2uS3JzkcJL9I2OnJNk7fA73Jtk6jB/1s9tBkjOTfCnJgeH9e80wvjH6raoNcwPePjL9SeDGYfoS4O+ZP/f8QuCRYfwU4PnhfuswvXW9+1hhrxcDW4bpG4AbhunzgH8D3gqcDTzH/EHi44bpdwPHD8uct959HEO/PwWcCzwIzI6Mb8p+R/rbFH0s6Ol9wAXA/pGxPwR2DdO7Rt7Pi352u9yA04ALhumTgK8P79kN0e+G2gKvqldGZk/k/y8Muhz4q5r3MHByktOADwF7q+rlqvovYC+wfU2LHlNV3V9VR4bZh5k/jx7me729qr5fVd8AnmX+awtaf3VBVR2oqsWuwt2U/Y7YLH38n6p6CHh5wfDlwC3D9C3AR0fGF/vstlBVh6rqsWH6u8AB5q9C3xD9bqgAB0jyB0leBH4N+L1heLHL9k9fYrybq5j/Vxs2f68LbfZ+N0sfy3lnVR2C+dAD3jGMb5r+k2wDzgceYYP0O8ml9GNJ8o/ATyzy0O6quruqdgO7k1wLfAK4jqNftr+iy/nXy3K9DsvsBo4At772tEWWLxb/x3bD9Aor63expy0y1qLfFdrQ79E1sCn6T/I24E7gU1X1SrJYW/OLLjK2av2ueYBX1QdWuOjfAPcyH+BHu2z/IPBLC8YfnLjIKVmu1+Gg66XARTXsQGPpryjY0F9dcAx/t6Pa9rtCb5avnHgpyWlVdWjYZXB4GG/ff5K3MB/et1bVXcPwhuh3Q+1CSXLOyOxlwFPD9D3AbwxHeC8EvjP8t+U+4OIkW4ejwBcPYxteku3Ap4HLqup7Iw/dA1yR5K1JzgbOAf6FzfvVBZu9383Sx3LuAV47C2wHcPfI+GKf3RYyv6l9E3Cgqj478tDG6He9j/IuOOJ7J7AfeAL4InD6MB7mfzziOeCrvP4shquYP/D1LPDx9e7hGHp9lvl9ZY8PtxtHHts99Po08OGR8UuYPwr+HPO7Jda9j2Po95eZ3zr5PvAScN9m7ndB75uij5F+bgMOAT8Y/k6vBn4ceAB4Zrg/ZVj2qJ/dDjfgF5jfBfLEyGf1ko3Sr5fSS1JTG2oXiiRp5QxwSWrKAJekpgxwSWrKAJekpgxwSWrKAJekpv4XHawn8ap4v7cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I should make v_n start at infinity, start with a big value and then check how it goes from here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 359,
   "metadata": {},
   "outputs": [],
   "source": [
    "#s_n=(np.zeros(n_samples)+1e9*2*np.pi)**(1/2)\n",
    "s_n= np.ones(n_samples)\n",
    "m_n=np.zeros(n_samples)\n",
    "v_n=np.zeros(n_samples)+np.inf\n",
    "#v_n=np.zeros(n_samples)+1e6\n",
    "\n",
    "#s_old=np.copy(s_n)\n",
    "#v_old=np.copy(v_n)\n",
    "\n",
    "##i create the prior\n",
    "m_0=prior_mean\n",
    "s_0=(2*np.pi*prior_var)**(-1/2)\n",
    "v_0=prior_var\n",
    "\n",
    "##i start with the prior as the posterior\n",
    "m_posterior=m_0\n",
    "v_posterior=v_0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:54: RuntimeWarning: invalid value encountered in double_scalars\n"
     ]
    }
   ],
   "source": [
    "thresh=1e-9\n",
    "conv=False\n",
    "iteration=0\n",
    "values_convergence=np.zeros(100)\n",
    "posterior_means=[]\n",
    "posterior_variance=[]\n",
    "\n",
    "posterior_means.append(m_posterior)\n",
    "posterior_variance.append(v_posterior)\n",
    "\n",
    "\n",
    "while not conv and iteration<100:\n",
    "    #print(\"\\n\");print(\"i am iter :\", iteration);print(\"\\n\")\n",
    "    \n",
    "    iteration+=1\n",
    "    m_old=np.copy(m_n)\n",
    "    v_old=np.copy(v_n)\n",
    "    s_old=np.copy(s_n)\n",
    "    ##with s i still have to understand if i have to divide for the normalization constant\n",
    "    for i in range(len(x)):\n",
    "        #if i%10==0:\n",
    "            #print(\"i am at example: \", i);\n",
    "        \n",
    "        #m_old, v_old = m_posterior, v_posterior\n",
    "        #if iteration==8:\n",
    "        #print(\"i is equal to: \", i)\n",
    "        if i%10==0 and iteration>10000:\n",
    "        ##i plot now, just for the first iteration to see how it goes\n",
    "            plt.figure(figsize=(15,10))\n",
    "            plt.title('i am ad datapoint: '+str(i))\n",
    "            xs=np.linspace(-10,10,10000)\n",
    "            ys=normal_pdf(xs, m_posterior, v_posterior)\n",
    "            plt.plot(xs,ys, color=\"blue\", label=\"current posterior estimate\")\n",
    "            plt.axvline(x=mean, color=\"red\", label=\"true generating value\")\n",
    "            aux=x[0:i]\n",
    "            aux2=np.zeros(len(aux))\n",
    "            plt.scatter(aux, aux2, marker=\"*\", color=\"green\", s=80, label=\"realization used\")\n",
    "            plt.legend()\n",
    "        \n",
    "        \n",
    "        if v_n[i]!=v_posterior:\n",
    "            #print(\"\\nv_n is: \", v_n[i])\n",
    "            #print(\"v_posterior: \", v_posterior)\n",
    "            #print(\"m_posterior: \", m_posterior)\n",
    "            cavity_m, cavity_v=calculating_cavities(m_n[i], v_n[i], m_posterior, v_posterior)\n",
    "            #print(\"cavity_m: \", cavity_m);print(\"cavity_v:\", cavity_v)\n",
    "        else:\n",
    "            #print(\"skipping i: \", i)\n",
    "            continue\n",
    "            \n",
    "        z=normalizer(w, x[i], cavity_m, cavity_v, clutter_var)\n",
    "        update_m, update_v = update(w, x[i], cavity_m, cavity_v, m_n[i], v_n[i], clutter_var)\n",
    "        new_v, new_m, new_s = new_marginal(update_m, update_v, cavity_m, cavity_v, z)\n",
    "        \n",
    "        v_n[i]=new_v\n",
    "        m_n[i]=new_m\n",
    "        s_n[i]=new_s\n",
    "        #print(\"new_w: \", new_v);print(\"new_m: \", new_m);print(\"new_s: \", new_s);\n",
    "        \n",
    "        m_posterior, v_posterior = update_m , update_v\n",
    "        posterior_means.append(m_posterior)\n",
    "        posterior_variance.append(v_posterior)\n",
    "        \n",
    "        #if i%10==0:\n",
    "            #print(\"mean posterior: \", m_posterior);print(\"variance posterior: \", v_posterior);print(\"\\n\")\n",
    "        #conv=convergence(m_posterior, v_posterior, m_old, v_old, thresh)\n",
    "    #if iteration==3:\n",
    "    #print(\"i am at iteration: \", iteration)\n",
    "    #print(\"m_old: \", m_old);print(\"m_n: \", m_n);print(\"\\n\");print(\"v_old: \", v_old);print(\"v_n: \", v_n);print(\"\\n\");\n",
    "    #print(\"s_old: \", s_old);print(\"s_n: \", s_n);print(\"\\n\")\n",
    "    ##now i calculate the convergence\n",
    "    aux=convergence(m_n, v_n, s_n, m_old, v_old, s_old, thresh)\n",
    "    conv=aux[0]\n",
    "    values_convergence[iteration-1]=aux[1]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I wanna plot now the evolution of the posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 361,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([           inf,            inf, 2.46977023e+37, 2.03137231e+37,\n",
       "       1.08838512e+35, 5.17723004e+32, 1.03782904e+31, 5.36144912e+28,\n",
       "       2.93798253e+27, 2.46743649e+23, 1.21600937e+23, 2.84217094e-14,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,\n",
       "       0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00])"
      ]
     },
     "execution_count": 361,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "values_convergence\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1198fb0f0>"
      ]
     },
     "execution_count": 362,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAb40lEQVR4nO3df3DU9b3v8ed7NwkhSCAExB9AE3s4IPLTBvRCpWitwuhRa+29UmY6FnuY1uttnc44pTO31TvTFm8vc+utpVruKVfr+PvHUI6HM1qvRaBwK3gFBw6owIk1154Cu+FHNiS72f3cPzYJIWzCbnazm3w/r8cMo/vd7373Hb7ktZ988v28v+acQ0REgi9U6gJERKQ4FPgiIp5Q4IuIeEKBLyLiCQW+iIgnFPgiIp4oK3UB/Rk/fryrq6srdRkiIsPKu+++e9w5N6H39iEd+HV1dezevbvUZYiIDCtm9nGm7ZrSERHxhAJfRMQTCnwREU8o8EVEPKHAFxHxhAJfRMQTQ/qyTJFiau9I0tLWwbhRFZhZQY8bNqMsfOHxVTLl+CTaysfRVmpHVTCltorqyvKC1SJ+G9qB/8EHsGRJqauQQZRIOk6eiZNIOqory6gaUUbhojazlHOcauugORanNZ4kkUyRSKZIptL3hqgoCzO2qpwxI9N/ykK5VdSRcpxojXO6rYPTbR20xjsAMDNCZvT3WZJMOZxzjABagH8Z4NfonQJ9QA/2v71czLi0mtGVhY3ooR34UjRHT7dz7HQ7I8pCVJaHqSzv+m8458DrjwNi7R2caE1wojVOS3vHOc+XhUPdQTtmZDkjyrKfdXRAR9LR3pHsDu+e4h0pmlvjnDiTIJVyhMwYNSL9IVMRNspCIcIh41RbgkhLO0dPtQGcM9o3g8ryMFXlYaoqwlSUhbuf60ilaG5NcOpMAuccoZAxekQZl9eMxDCSzpHKUFdPoZAxsvPvvyPpaOtIkkxl/VdQBEPwhkkFKmmofWUVOfzbz9bQDvxp02DLllJX4YXVT+5i5+EIY6vK+cvJtnOeq6kqp278KOpqR/GZ2irqx4/iM7WjmFwzkvIs/lG2JZL8nyNRthw8ytsfHiMSi2MGcyaNZcm0CVw/7WIurh7BjkMRth86zraPjnO8pR2Av7n4IhZ+tpYxIzNPa5xu66CpuZVPomdoam4lFk/2W8vFo0fwxSsv5sYrJ7Lws+MZWRE+b59LgUQyxXt/PsHOwxFaE2c/lNoTKY4cj/HBv53ir6faz3vtlHFVLJt5CTdddQlzJ48lXMAPS5Gs9fETjw3lWxw2NDQ4tVYojtvX/ZHqyjKevvca2hJJPo600hiJ8XEkxr8eb+XjSIyPI618evIMA/0nU1NVzuK/TQf8dVPHU3vRiIz7Oec4+G+n2f7RcbYdOs7uxihticxBPrI8zORxVUyqqWLyuJFMrqliUs1IxowsP28e/qIRZUy/ZDShAoXwidZ49wcTwIiyMJNqRhZ0/l9kIMzsXedcQ+/tQ3uEL0UTjbVTX1sFpKcspl0ymmmXjD5vv7ZEkk+irTRGWmlqbs04ddJbOGTMmTyWOZOyG/GaGVdeWs2Vl1bz94uvyP2LKZKxVRWMraoodRkiWVPgCwDRljjjRmUecfdUWR5m6sTRTJ14/oeBiAxtug5faEskicWT1F6k0apIkCnwhebWOAA1mp4QCTQFvhBpSQf+uFEKfJEgU+AL0Vg68DWlIxJsCnzpDnyN8EWCTYEvRLpG+Ap8kUBT4AvRWDvhkKlJl0jAKfCFaCxBTVV5wVagisjQpMAXorF2zd+LeECBL0RjcQW+iAcU+EIkFqc2i7YKIjK8KfBFI3wRTyjwPdeRTHGiNaHAF/GAAt9zJ84kAK2yFfGBAt9zXats1ThNJPgU+J7rapymVbYiwafA91x3Hx1N6YgEXtEC38yuMLPfmNnLxXpPubBoLH1PVv3SViT4sgp8M9tgZkfNbF+v7UvN7AMzO2Rmq/s7hnPuiHPu3nyKlcKLaA5fxBvZ3tP2SeCXwG+7NphZGFgHfAloAnaZ2SYgDKzp9fqVzrmjeVcrBReNxRkzspzysGb3RIIuq8B3zm01s7pemxcAh5xzRwDM7HngdufcGuDWQhYpg0eLrkT8kc+w7nLgkx6Pmzq3ZWRmtWb2BDDPzH7Qz36rzGy3me0+duxYHuVJNhT4Iv7Idkonk0y9dF1fOzvnIsC3LnRQ59x6YD1AQ0NDn8eTwojG4kweV1XqMkSkCPIZ4TcBk3s8ngR8ml85Umzpxmka4Yv4IJ/A3wVMNbN6M6sA7gY2FaYsKQbnHM2a0hHxRraXZT4H7ASmmVmTmd3rnOsA7gdeBw4ALzrn9g9eqVJop8500JFyCnwRT2R7lc7yPrZvBjYXtCIpmmhr5ypbBb6IF3Txtce0ylbELwp8j51tnKa7XYn4QIHvMTVOE/GLAt9jXX10dFmmiB8U+B6LxuJUVYSpLA+XuhQRKQIFvsd0Db6IXxT4Hoso8EW8osD3mBqnifhFge8xBb6IXxT4HovE2nWFjohHFPieao130JZIMU6LrkS8ocD31NlVthrhi/hCge+p5s7GaTUKfBFvKPA91bXKVr+0FfGHAt9TUU3piHhHge8pNU4T8Y8C31ORWJzysDF6RD73sReR4USB76lorJ1xoyows1KXIiJFosD3VDSW0DX4Ip5R4HsqPcIvL3UZIlJECnxPpfvoaIQv4hMFvqcisbguyRTxjALfQ/GOFKfbOrToSsQzCnwPdbVVUOCL+EWB7yE1ThPxkwLfQ2qcJuInBb6HuhqnaYQv4hcFvoeiLe2A5vBFfKPA91A0FscMxlYp8EV8UrTAN7MrzewJM3vZzL5drPeV80VicWqqKgiH1EdHxCdZBb6ZbTCzo2a2r9f2pWb2gZkdMrPV/R3DOXfAOfct4N8DDQMvWfKVXmWr0b2Ib7Id4T8JLO25wczCwDpgGTADWG5mM8xslpm91uvPxZ2vuQ3YDvzvgn0FkjMFvoifsmqG7pzbamZ1vTYvAA45544AmNnzwO3OuTXArX0cZxOwycz+CXg20z5mtgpYBTBlypRsypMcRWNxPjvholKXISJFls8c/uXAJz0eN3Vuy8jMlpjZL8zs18DmvvZzzq13zjU45xomTJiQR3nSl2gsrjtdiXgon9sdZfqNn+trZ+fcFmBLHu8nBZBKOZpb1ThNxEf5jPCbgMk9Hk8CPs2vHBlsJ84kSDldgy/io3wCfxcw1czqzawCuBvYVJiyZLBEY1p0JeKrbC/LfA7YCUwzsyYzu9c51wHcD7wOHABedM7tH7xSpRDONk7TzU9EfJPtVTrL+9i+mX5+AStDz9nGabq9oYhv1FrBM2cbp2mEL+IbBb5noi0a4Yv4SoHvmUgszugRZYwoC5e6FBEpMgW+Z7ToSsRfCnzPqI+OiL8U+J6JxrTKVsRXCnzPRDt74YuIfxT4HnHOaQ5fxGMKfI+0tHcQT6Y0pSPiKQW+R6Kdi67GadGViJcU+B45u8pWI3wRHynwPdLcPcJX4Iv4SIHvkYgCX8RrCnyPRBX4Il5T4HskGoszoixEVYX66Ij4SIHvkUhLepWtWabbEYtI0CnwPRKNtWvRlYjHFPgeSTdO0zX4Ir5S4Hsk2qrGaSI+U+B7JNqixmkiPlPge6ItkSQWT1KrOXwRbynwPaFr8EVEge8JBb6IKPA9ocZpIqLA94Qap4mIAt8TapwmIgp8T0Rj7YRDRnVlealLEZESUeB7ouvm5aGQ+uiI+KpogW9mS8xsm5k9YWZLivW+ktbVOE1E/JVV4JvZBjM7amb7em1famYfmNkhM1t9gcM4oAWoBJoGVq4MVLqPjgJfxGdlWe73JPBL4LddG8wsDKwDvkQ6wHeZ2SYgDKzp9fqVwDbn3NtmNhH478CK/EqXXERjca68rLrUZYhICWUV+M65rWZW12vzAuCQc+4IgJk9D9zunFsD3NrP4ZoBtWwsMjVOE5FsR/iZXA580uNxE3BNXzub2Z3AzcBY0j8t9LXfKmAVwJQpU/IoT7p0JFOcaE2ocZqI5/IJ/EyXe7i+dnbOvQq8eqGDOufWA+sBGhoa+jyeZK+5NQGgxmkinsvnKp0mYHKPx5OAT/MrRwaD+uiICOQX+LuAqWZWb2YVwN3ApsKUJYUUibUDCnwR32V7WeZzwE5gmpk1mdm9zrkO4H7gdeAA8KJzbv/glSoDFe1unKbflYv4LNurdJb3sX0zsLmgFUnBqXGaiIBaK3ihq3Ha2Cr10RHxmQLfA9FYnDEjyykP63SL+EwJ4IFITIuuRESB74Voi/roiIgC3wtqnCYioMD3QrQ1rlW2IqLADzrnHM0a4YsICvzAO3Wmg46UU+M0EVHgB11XWwVN6YiIAj/gzjZOU1sFEd8p8AMu0t1HRyN8Ed8p8ANOrZFFpIsCP+AU+CLSRYEfcNFYnKqKMJXl4VKXIiIlpsAPOK2yFZEuCvyAU+M0EemiwA+4aKxdI3wRART4gZfulKlr8EVEgR94apwmIl0U+AHWGu+gLZHSlI6IAAr8QIu0dF6Dr8ZpIoICP9C06EpEelLgB1h34GsOX0RQ4AeaGqeJSE8K/ACLdvbC15SOiIACP9CisQQV4RAXjSgrdSkiMgQo8AMsGmunZlQ5ZlbqUkRkCFDgB1i6cZpW2YpImgI/wNQ4TUR6KtrkrpldB6zofM8ZzrmFxXpvX0VjcSbXVJW6DBEZIrIa4ZvZBjM7amb7em1famYfmNkhM1vd3zGcc9ucc98CXgOeGnjJkq104zSN8EUkLdsR/pPAL4Hfdm0wszCwDvgS0ATsMrNNQBhY0+v1K51zRzv//2vAN/OoWbIQ70hxur1DUzoi0i2rwHfObTWzul6bFwCHnHNHAMzseeB259wa4NZMxzGzKcBJ59ypvt7LzFYBqwCmTJmSTXmSQXNretFVjQJfRDrl80vby4FPejxu6tzWn3uB/9XfDs659c65Budcw4QJE/Ioz29djdM0wheRLvn80jbTxd2uvxc45x7K4/0kB2qcJiK95TPCbwIm93g8Cfg0v3KkUCKdbRV08xMR6ZJP4O8CpppZvZlVAHcDmwpTluTr7AhfC69EJC3byzKfA3YC08ysyczudc51APcDrwMHgBedc/sHr1TJRXMsTshg7MjyUpciIkNEtlfpLO9j+2Zgc0ErkoKIxOLUVFUQCqmPjoikqbVCQEVjcV2SKSLnUOAHVCSmVbYici4FfkBF1ThNRHpR4AdUVCN8EelFgR9AyZSjuVUjfBE5lwI/gE6eSeCcVtmKyLkU+AHUdfNyXaUjIj0p8APobOM0rbIVkbMU+AGkxmkikokCP4AinYGvxmki0pMCP4C6Rvg1VQp8ETlLgR9A0Vic0ZVlVJTp9IrIWUqEANIqWxHJRIEfQGqcJiKZKPADKKIRvohkoMAPoGisXZdkish5FPgB45zrbJymRVcici4FfsCcbu8gkXSa0hGR8yjwA6ZZq2xFpA8K/ICJKPBFpA8K/ICJtijwRSSzslIXIIWlxmmSq0QiQVNTE21tbaUuRXJUWVnJpEmTKC8vz2p/BX7AqHGa5KqpqYnRo0dTV1eHmZW6HMmSc45IJEJTUxP19fVZvUZTOgETjbVTWR6iqkKf5ZKdtrY2amtrFfbDjJlRW1ub009mCvyAicYSuvGJ5ExhPzzlet4U+AGjVbYi0hcFfsCocZpI8TU2NvLss88O6LULFy4scDV9U+AHjBqniRTfQAI/mUwCsGPHjpxfM1BFC3wzm2FmL5rZ42Z2V7He1zfpPjoKfBk+GhsbmT59Ot/85jeZOXMmK1as4M0332TRokVMnTqVd955B4BYLMbKlSuZP38+8+bN43e/+13366+77jquvvpqrr766u4A3bJlC0uWLOGuu+5i+vTprFixAufcee+/ZMkSHnjgARYuXMjMmTO73y8ajXLHHXcwe/Zsrr32Wt5//30A3n77bebOncvcuXOZN28ep0+fZvXq1Wzbto25c+fy85//nGQyyYMPPsj8+fOZPXs2v/71r7truv766/na177GrFmzALjooouA9FU3Dz74IDNnzmTWrFm88MILfb5moLK6lMPMNgC3AkedczN7bF8K/A8gDPyDc+6Rfg6zDHjMObfNzDYBLw+8bMmkLZGkNZ5U4MuA/Zd/3M+/fHqqoMeccVk1D/3dVf3uc+jQIV566SXWr1/P/PnzefbZZ9m+fTubNm3ipz/9KRs3buQnP/kJN9xwAxs2bODEiRMsWLCAG2+8kYsvvpjf//73VFZW8tFHH7F8+XJ2794NwHvvvcf+/fu57LLLWLRoEX/84x/5/Oc/f977x2IxduzYwdatW1m5ciX79u3joYceYt68eWzcuJG33nqLr3/96+zZs4e1a9eybt06Fi1aREtLC5WVlTzyyCOsXbuW1157DYD169czZswYdu3aRXt7O4sWLeKmm24C4J133mHfvn3nXUr56quvsmfPHvbu3cvx48eZP38+ixcv7vc1ucr22r0ngV8Cv+3aYGZhYB3wJaAJ2NUZ5GFgTa/XrwSeBh4ys9uA2ryqloy6Fl1pSkeGm/r6+u7R61VXXcUXv/hFzIxZs2bR2NgIwBtvvMGmTZtYu3YtkL6c9M9//jOXXXYZ999/P3v27CEcDvPhhx92H3fBggVMmjQJgLlz59LY2Jgx8JcvXw7A4sWLOXXqFCdOnGD79u288sorANxwww1EIhFOnjzJokWL+N73vseKFSu48847u4/f0xtvvMH777/Pyy+nx7UnT57ko48+oqKiggULFmQM7u3bt7N8+XLC4TATJ07kC1/4Art27aK6urrP1+Qqq8B3zm01s7pemxcAh5xzRwDM7HngdufcGtI/DWTyHzs/KF4dWLnSH62ylXxdaCQ+WEaMOHspcSgU6n4cCoXo6OgA0lMer7zyCtOmTTvntQ8//DATJ05k7969pFIpKisrMx43HA53H6u33pc3mlnG6R8zY/Xq1dxyyy1s3ryZa6+9ljfffPO8/ZxzPPbYY9x8883nbN+yZQujRo3KWEOm9+vS12tylc8c/uXAJz0eN3Vuy8jM6sxsPemfEv5bP/utMrPdZrb72LFjeZTnHzVOkyC7+eabeeyxx7qD8b333gPSo+dLL72UUCjE008/PaBfbHbNl2/fvp0xY8YwZswYFi9ezDPPPAOkg3r8+PFUV1dz+PBhZs2axfe//30aGho4ePAgo0eP5vTp0+fU+vjjj5NIJAD48MMPicVi/dawePFiXnjhBZLJJMeOHWPr1q0sWLAg56+lP/ksx8x0xX+fH1HOuUZg1YUO6pxbD6wHaGho6PsjT84TjbUDCnwJph/+8Ic88MADzJ49G+ccdXV1vPbaa9x333185Stf4aWXXuL6668f0Gi4pqaGhQsXcurUKTZs2ACkf3L4xje+wezZs6mqquKpp54C4NFHH+UPf/gD4XCYGTNmsGzZMkKhEGVlZcyZM4d77rmH7373uzQ2NnL11VfjnGPChAls3Lix3xq+/OUvs3PnTubMmYOZ8bOf/YxLLrmEgwcP5v6X1Qfr78eIc3ZMT+m81vVLWzP7d8DDzrmbOx//AKBzSqcgGhoaXNcvX+TC/mHbEX78TwfY+6ObGFOVXTMlkQMHDnDllVeWuoySWbJkCWvXrqWhoaHUpQxIpvNnZu865877gvKZ0tkFTDWzejOrAO4GNuVxPMlTNBanLGRUj1QfHRE5X7aXZT4HLAHGm1kT8JBz7jdmdj/wOukrczY45/YPWqVyQV2rbNUXRSR7W7ZsKXUJRZPtVTrL+9i+Gdhc0IpkwKJaZSsi/VBrhQDRKlsR6Y8CP0DUOE1E+qPADxA1ThOR/ijwAyKRTHHyTEJTOjLsnDhxgl/96lelLsMLCvyAaG5VHx0ZnvoL/HzbAcu5FPgB0RxLL+Eep9sbyjCzevVqDh8+zNy5c3nwwQfPawfc2NjIzJndTXpZu3YtDz/8MACHDx9m6dKlfO5zn+O6667rXpV6zz338J3vfIeFCxdyxRVXdDcx66nUbZlLQSt0AiKitgpSCA88AHv2FPaYc+fCo4/2+fQjjzzCvn372NP5vlu2bDmnHXBXt8xMVq1axRNPPMHUqVP505/+xH333cdbb70FwF/+8he2b9/OwYMHue2227jrrvNvw1HqtszFpsAPCHXKlCDJph1wS0sLO3bs4Ktf/Wr3tvb29u7/v+OOOwiFQsyYMYO//vWvGY9R6rbMxabADwgFvhREPyPxYurZAK2srIxUKtX9uK2tDYBUKsXYsWO7fzLorWdr5L6mVErdlrnYNIcfEJGWdODXqGmaDDO9Wwv3NnHiRI4ePUokEqG9vb37rlLV1dXU19fz0ksvAelg3rt3b8HrG8y2zMWmEX4RJVOOn24+QHPnaLyQ3v9/JxlbVU5ZWJ/hMrzU1tayaNEiZs6cybJly7jlllvOeb68vJwf/ehHXHPNNdTX1zN9+vTu55555hm+/e1v8+Mf/5hEIsHdd9/NnDlzClrfYLZlLras2yOXQtDaI//r8RjXr93CuFEVjBoRLvjxP/8341lz5+yCH1eCzff2yMNdLu2RNcIvoq4P14f+bga3z+3z5mAiIoNCP/8X0dD9WUpEfKDAF5EhszBIcpPreVPgF5G+p2QoqqysJBKJKPSHGecckUjknMtBL0Rz+CWgO1LJUDJp0iSampo4duxYqUuRHFVWVnYv8MqGAr+oNIKSoae8vPyCq1olGDSlUwIa34tIKSjwi0hTpCJSSgr8EtAUvoiUwpBeaWtmx4CPOx+OAU5m2K2v7eOB44NUWr76qrnUx8319bnsf6F9B/r8cDv/g3XuC3HswTr/2ezX3z4Dec738/8Z59yE8/Zwzg2LP8D6HLfvLnXNuX4tpT5urq/PZf8L7TvQ54fb+R+scz+Uz382+/W3z0Ce0/nP/Gc4Ten8Y47bh7LBqjnf4+b6+lz2v9C+A31+uJ3/wax3qJ7/bPbrb5+BPjcUlfT8D+kpnXyY2W6XoXmQ+EHn3286/5kNpxF+rtaXugApKZ1/v+n8ZxDYEb6IiJwryCN8ERHpQYEvIuIJBb6IiCe8CXwzG2VmT5nZ/zSzFaWuR4rLzK4ws9+Y2culrkWKz8zu6Pze/52Z3VTqekplWAe+mW0ws6Nmtq/X9qVm9oGZHTKz1Z2b7wReds79PXBb0YuVgsvl/Dvnjjjn7i1NpTIYcjz/Gzu/9+8B/kMJyh0ShnXgA08CS3tuMLMwsA5YBswAlpvZDGAS8Ennbski1iiD50myP/8SPE+S+/n/z53Pe2lYB75zbisQ7bV5AXCoc0QXB54HbgeaSIc+DPOvW9JyPP8SMLmcf0v7r8A/O+f+b7FrHSqCGHyXc3YkD+mgvxx4FfiKmT3O8FuOLdnLeP7NrNbMngDmmdkPSlOaFEFf3///CbgRuMvMvlWKwoaCIN7xKlPzYeeciwHfKHYxUnR9nf8I4O03ukf6Ov+/AH5R7GKGmiCO8JuAyT0eTwI+LVEtUnw6/37T+e9HEAN/FzDVzOrNrAK4G9hU4pqkeHT+/abz349hHfhm9hywE5hmZk1mdq9zrgO4H3gdOAC86JzbX8o6ZXDo/PtN5z93ap4mIuKJYT3CFxGR7CnwRUQ8ocAXEfGEAl9ExBMKfBERTyjwRUQ8ocAXEfGEAl9ExBMKfBERT/x/KhAF9oIAYsgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(posterior_means, label=\"mean posterior\")\n",
    "plt.axhline(2, color=\"red\", label=\"truen mean\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 363,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x120e73940>"
      ]
     },
     "execution_count": 363,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de1hVZd438O+9OSpnOYgIchA5iiAHKzVFKg+F+XQwtbHRJjWtZnqaamyueabsmWZ6m5yx8almXhtT62l0zLcxm7Q0xdTSFE1LUZQAFTNFUJSTsNm/9w+UUPcGNrD3Wuz9/VwXl+51/LEW+8viXve+lxIREBGR4zNoXQAREdkHA5+IyEkw8ImInAQDn4jISTDwiYicBAOfiMhJuGpdQFuCgoIkKipK6zKIiHqMvXv3nhORYHPzdB34UVFRyM/P17oMIqIeQyl13NI8NukQETkJBj4RkZNg4BMROQldt+GT42lsbERZWRnq6+u1LoWoR/P09ER4eDjc3Nw6vA4Dn+yqrKwMPj4+iIqKglJK63KIeiQRQUVFBcrKyhAdHd3h9dikQ3ZVX1+PwMBAhj1RFyilEBgYaPVfyrq+wr9Y34hNBWe0LqPHuCmmD3w9O/7nnVYY9kRd15n3kS4DXyk1EcBE99BYzH6H/fA7KqW/Hz58fAQMBgaqPXh7e6O6urrNZV577TXMmTMHvXv3BgDceeed+Mc//gF/f3+b1FReXo7c3Fw0NDRg8eLFuPXWW22yH73uXytr165FXFwckpKSrFpv3bp1KCgowHPPPWejyq6l9PwAlOTUofLP9Vu1LqNH2FVcgZc+Poy/TE3DpLT+Wpdj0eHDh5GYmKh1Gd2iI4F/9cODQUFBdqlp1apV2LBhA1asWNGt221qaoKLi4tm+7eko3XZ2syZM5Gbm4v777+/w+sYjUa4ulp3zX39OubeT0qpvSKSaXYDIqLbr4yMDKGOMTaZZNyiz2XkK5vlcmOT1uVYVFBQoHUJ8u6770pWVpakpqbKnDlzxGg0yptvvinPPvtsyzLLli2TJ554QkRE/vSnP0lycrIkJyfLokWLWpbx8vISEZG8vDy56667WqY//vjjsmzZMvnLX/4ibm5uMnjwYMnOzhYRkcjISCkvL7e43ZKSEklISJBZs2ZJUlKS3HHHHVJbW3vD91BaWio5OTmSkpIiOTk5cvz4cfn6668lIiJCgoKCJDU19Zr11q9fL5MnT255nZeXJ7m5uSIiMnfuXMnIyJCkpCR5/vnnW5aJjIyUF198UUaMGCErV67s0v5FRF588UXJzMyU5ORkmT17tphMJhERGT16tDz55JNyyy23SHJysnz11VciIvLCCy/I9OnTZcyYMRIbGytLlixpqT07O1umTZsmiYmJFo/l7t27JSUlRerq6qS6ulqSkpLk22+/vaamkpISiY+Pl5/+9KeSkpIi9913n9TU1IiIyGeffSZpaWkyePBgefjhh6W+vl5ERObPny+JiYmSkpIiTz/9tHzxxRcSEBAgUVFRkpqaKkVFRVJUVCTjxo2T9PR0GTlypBw+fFhERGbMmCFPPfWUZGdnyy9/+UtZtmyZPP744xaPqbl1WjP3fgKQLxYyVfNQb+uLgW+dvCNnJHL+v+XtHcVal2JR6x/QBesOygN/+7JbvxasO9ju/nNzc6WhoUFERObNmycrVqyQs2fPysCBA1uWGz9+vGzfvl3y8/Nl8ODBUl1dLZcuXZKkpCTZt2+fiLQf+CLXBnzr15a2W1JSIi4uLvL111+LiMjkyZPl3XffveH7yM3NleXLl4uIyNKlS2XSpEkiItcESGuNjY0SEREh1dXVItIc8le3W1FRISIiRqNRRo8eLQcOHGip9ZVXXjF7HK3df+v9iIhMnz5d1q1bJyLNgT9r1iwREfn8888lOTlZRJoDf8iQIVJbWyvl5eUSHh4up06dkry8POndu7cUFzf/nLd1jn7zm9/I008/LY899pj84Q9/uKGmkpISASA7duwQEZGHH35YXn31Vamrq5Pw8HApLCwUEZGHHnpIFi1aJBUVFRIXF9fyy+r8+fMi0hzK77//fst2c3Jy5OjRoyIismvXLhkzZkzLcnfddZcYjcYbjpelY3r9Oq1ZG/jspeNARscF45aYQPzPliJcqm/Uuhxd2rx5M/bu3YusrCykpaVh8+bNKC4uRnBwMGJiYrBr1y5UVFSgsLAQI0aMwI4dO3DPPffAy8sL3t7euPfee7F9+/Yu19HWdqOjo5GWlgYAyMjIQGlp6Q3r79y5Ew8++CAA4KGHHsKOHTva3J+rqyvGjx+Pjz76CEajER9//DEmTZoEAFi9ejXS09MxdOhQHDp0CAUFBS3rTZkyxez2rN0/AOTl5eGmm25CSkoKtmzZgkOHDrXMmzZtGgBg1KhRuHjxIi5cuAAAmDRpEnr16oWgoCCMGTMGu3fvBgAMGzaspTtiW8fy+eefx6ZNm5Cfn49f/epXZuuKiIjAiBEjAADTp0/Hjh07UFhYiOjoaMTFxQEAZsyYgW3btsHX1xeenp6YNWsWPvjgg5Z7M61VV1fjyy+/xOTJk5GWloZHH30Up0+fbpk/efJks81QbR1TS+tYS5c3balzlFL49Z0JuPv1L7BkWzGeHhuvdUltemFist33KSKYMWMGXn755RvmTZkyBatXr0ZCQgLuueceKKWa/wxuh6urK0wmU8vrjnSVa2u7Hh4eLf93cXFBXV1du9vrSI+NKVOm4I033kCfPn2QlZUFHx8flJSUYOHChdizZw8CAgIwc+bMa+r38vJqd7sd2X99fT0ee+wx5OfnIyIiAgsWLLhmP9evf/W1pemt62rrWFZWVqK6uhqNjY2or683+/2Y24elbbq6umL37t3YvHkzVq1ahddffx1btmy5ZhmTyQR/f3/s37/f7DY6c0w7uk57eIXvYIaE+yN3SD/8fXsJzl7kp1mvd9ttt2HNmjU4e/YsgOZAOH68eXDBe++9F2vXrsXKlStbrmxHjRqFtWvXora2FjU1NfjXv/51Q8+TyMhIFBQU4PLly6iqqsLmzZtb5vn4+ODSpUs31NGR7bZl+PDhWLVqFQDgvffew8iRI9tdJzs7G/v27cNbb73V8v1dvHgRXl5e8PPzw5kzZ7Bhwwab7P9quAcFBaG6uhpr1qy5Zv4///lPAM1X635+fvDz8wMAfPjhh6ivr0dFRQW2bt2KrKysG7bd1rGcM2cOfve73+EnP/kJ5s+fb7a2EydOYOfOnQCAlStXYuTIkUhISEBpaSmKiooAAO+++y5Gjx6N6upqVFVV4c4778Rrr73WEuqtz7Ovry+io6Px/vvvA2j+hXTgwIE2jw/QuXNqLV7hO6Bnxsbjk4M/4LXNx/CHe1K0LkdXkpKS8NJLL2Hs2LEwmUxwc3PDG2+8gcjISAQEBCApKQkFBQUYNmwYACA9PR0zZ85seT1r1iwMHTr0mm1GRETggQcewJAhQzBo0KBr5s+ZMwcTJkxAv379kJeX1zLd0nbNNd+Ys3jxYvzsZz/Dq6++iuDgYCxbtqzddVxcXJCbm4vly5e39KJJTU3F0KFDkZycjJiYmJamje7ev7+/P2bPno2UlBRERUXdENwBAQEYPnw4Ll68iLfffrtl+rBhw3DXXXfhxIkT+O1vf4uwsDAcPXr0mnUtHct33nkHrq6uePDBB9HU1IThw4djy5YtyMnJuWb9xMRErFixAo8++igGDRqEefPmwdPTE8uWLcPkyZNhNBqRlZWFuXPnorKyEpMmTUJ9fT1EBIsWLQIATJ06FbNnz8bixYuxZs0avPfee5g3bx5eeuklNDY2YurUqUhNTe3WY9oZuu6WmZmZKRwPv3Ne+PAg/verE9j41CgMDPbWupwWjtQtk7pHdnY2Fi5ciMzMa3sSLliwAN7e3njmmWdstu/S0lLk5ubi4MGDNtuHLVnbLZNNOg7q57cNgqerAa9+Uqh1KUSkE2zScVBB3h6YM2ogFn12FHuPn0dGZIDWJRGZtXXrVrPTFyxYYPN9R0VF9dir+87gFb4Dm3VrNIK8PfDKhiMd6m1CRI6Nge/AvDxc8eTtg7C7tBKbD5/VupwW/OVD1HWdeR8x8B3c1KwIxAR54ZVPjsDYZGp/BRvz9PRERUUFQ5+oC+TKePienp5Wrcc2fAfn5mLAs+PiMe+9ffhg3yk8kBWhaT3h4eEoKytDeXm5pnUQ9XRXn3hlDV0G/tXhkWNjY7UuxSGMHxyKtAh//HnTUUxMDUMvd+1GF3Rzc7PqCT1E1H102aQjIh+JyJyrn7ajrlFK4bkJCfjhYj2Wf1mqdTlEpBFdBj51v5tjApGTEII3txbhfE2D1uUQkQYY+E5k/vgE1Fw24s2tRVqXQkQaYOA7kfhQH9yXHo4VXx5H2flarcshIjtj4DuZp+6Ig1LAnzcebX9hInIoDHwnE+bfCzNHROFf+0+h4PuLWpdDRHbEwHdCj42Oha+nG1755IjWpRCRHTHwnZBfbzc8PmYgPj9aji+LzmldDhHZCQPfSf30liiE+Xni/3xyBCYThzkgcgYMfCfl6eaCX46NxzdlVfj429Ptr0BEPR4D34ndM7Q/EkJ9sHBjIRqM2g+sRkS2xcB3Yi4GhfkTEnC8ohYrd5/QuhwisjEGvpPLjgvGzTF9sHjzMVRfNmpdDhHZEAPfySml8OsJiaioacDS7SVal0NENsTAJ6RG+GNEbCA++uZ7rUshIhti4BMAICehL4rOVuNkJcfYIXJUDHwCAOQkhAAA8gr18+xbIupeDHwCAEQHeSEqsDfyjjDwiRwVA59aZMeH4MvvKlDX0KR1KURkAwx8apGTEILLRhN2FVdoXQoR2YDdAl8pFaOUWqqUWmOvfZJ1hkX3QS83F7bjEzmoLgW+UuptpdRZpdTB66aPV0oVKqWKlFLPAYCIFIvII13ZH9mWp5sLRsQGYcuRsxDhgGpEjqarV/jLAYxvPUEp5QLgDQATACQBmKaUSurifshOxiQEo+x8Hb4rr9a6FCLqZl0KfBHZBqDyusnDABRduaJvALAKwKSu7IfsJzu+uXvmFvbWIXI4tmjD7w/gZKvXZQD6K6UClVJ/AzBUKfVrSysrpeYopfKVUvnl5eU2KI/a0t+/FxJCfZB3hMeeyNHYIvCVmWkiIhUiMldEBorIy5ZWFpElIpIpIpnBwcE2KI/akx0fgj2llbhU36h1KUTUjWwR+GUAIlq9DgfAQVp6kJyEEBhNgh3H+PhDIkdii8DfA2CQUipaKeUOYCqAdTbYD9lI+gB/+Hi6snsmkYPparfMlQB2AohXSpUppR4RESOAJwB8CuAwgNUicsjK7U5USi2pqqrqSnnUSa4uBoyKC0ZeYTmfd0vkQLraS2eaiPQTETcRCReRpVemrxeRuCvt9b/vxHY/EpE5fn5+XSmPuiAnPgTlly6j4PRFrUshom7CoRXIrNHxwVAKHEyNyIEw8MmsIG8PDAn3xxa24xM5DF0GPtvw9WFMfDD2n7yAypoGrUshom6gy8BnG74+jIkPgQjw+VFe5RM5Al0GPulDSn8/BHm781O3RA6CgU8WGQwKo+NC8PnRcjSxeyZRj8fApzblJISgqq4RX584r3UpRNRFugx83rTVj5GDguBiUPzULZED0GXg86atfvj1ckNGZADb8YkcgC4Dn/QlJyEEBacv4oeqeq1LIaIuYOBTu8ZceSjKVjbrEPVoDHxqV1xfb/T378WnYBH1cLoMfN601RelFLLjg/FF0TlcNjZpXQ4RdZIuA583bfVnTHwIahqakF/K7plEPZUuA5/0Z3hsINxdDWzWIerBGPjUIb3dXXFzTCD74xP1YAx86rCc+GAUl9fgeEWN1qUQUScw8KnDsq90z+RDUYh6JgY+dVhUkBdigryQV8hP3RL1RLoMfHbL1K8xCSHYWVyB2gaj1qUQkZV0GfjslqlfY+JD0GA0Yed3FVqXQkRW0mXgk35lRQfAy92F3TOJeiAGPlnFw9UFI2KDsLWwHCJ8KApRT8LAJ6uNSQjBqQt1OHa2WutSiMgKDHyy2tXRM9msQ9SzuGpdAPU8oX6eSOzni42HfsDtiX3bXd7H0xV9fT3tUBkRtYWBT51ye2II/mdLEW7/8+ftLmtQQN4z2YgM9LJDZURkiS4DXyk1EcDE2NhYrUshCx4dPRAJob5oaufG7YXaBjz/4SHsKq5g4BNpTOm5p0VmZqbk5+drXQZ1gYgg/XebcEdSX/zx/lStyyFyeEqpvSKSaW4eb9qSTSmlkBEZgL3HOY4+kdYY+GRz6ZEB+K68BudrGrQuhcipMfDJ5jIGBAAAvj7Jq3wiLTHwyeaGhPvD1aDYrEOkMQY+2Vwvdxckh/ky8Ik0xsAnu0iPDMCBk1VobDJpXQqR02Lgk11kRAagrrEJR05f0roUIqfFwCe7yIhsvnG793ilxpUQOS9dBj6feOV4+vn1QpifJ/aeuKB1KUROS5eBzydeOab0yADs441bIs3oMvDJMWVEBuDUhTqcrqrTuhQip8TAJ7u52o6/7zibdYi0wMAnu0ns5wtPNwP74xNphIFPduPmYkBquD/2nmDgE2mBgU92lREZgEOnqlDf2KR1KUROh4FPdpURGQCjSfBNGbvcEtkbA5/sauiAqx/AYrMOkb0x8Mmu+ni5IybYi4FPpAEGPtldxoAA7DtxHnp+vCaRI2Lgk91lRAagsqYBJedqtC6FyKkw8MnufhxIjc06RPbEwCe7GxjsDV9PV+xjf3wiu2Lgk90ZDArpkQG8wieyM10GPodHdnwZAwJw9Ew1quoatS6FyGnoMvA5PLLju9qO/zWbdYjsRpeBT44vNcIfBgWOj09kRwx80oSXhysS+/lyIDUiO2Lgk2YyIgOw/8QFGJtMWpdC5BQY+KSZjMgA1DQ0ofDMJa1LIXIKDHzSTPqAq0/AYrMOkT0w8Ekz4QG9EOLjwf74RHbCwCfNKKWQERnAG7dEdsLAJ01lRAbgZGUdzl6s17oUIofHwCdNpV/5ABbH1SGyPQY+aWpwmB/cXQ1sxyeyAwY+acrd1YC0CH98UVShdSlEDo+BT5q7PTEEBacvoux8rdalEDk0Bj5p7o6kUADApoIzGldC5NgY+KS56CAvxPX1xqeHftC6FCKHxsAnXRibFIrdJZU4X9OgdSlEDouBT7owLjkUJgE2HzmrdSlEDouBT7owuL8v+vl5slmHyIYY+KQLSimMTeqL7cfKUdfQpHU5RA6JgU+6MS45FPWNJmw7Vq51KUQOyW6Br5TyUkqtUEq9pZT6ib32Sz1HVnQf+PVyw8ZD7J5JZAtdCnyl1NtKqbNKqYPXTR+vlCpUShUppZ67MvleAGtEZDaAu7uyX3JMbi4G3JYQgs1HzvApWEQ20NUr/OUAxreeoJRyAfAGgAkAkgBMU0olAQgHcPLKYmykJbPGJvfFhdpG7C6t1LoUIofTpcAXkW0Arn9nDgNQJCLFItIAYBWASQDK0Bz6Xd4vOa5RccHwcDWwWYfIBmwRvP3x45U80Bz0/QF8AOA+pdRfAXxkaWWl1BylVL5SKr+8nDfvnE1vd1fcOigYmwrOQES0LofIodgi8JWZaSIiNSLysIjME5H3LK0sIktEJFNEMoODg21QHund2OS+OHWhDoe+v6h1KUQOxRaBXwYgotXrcADf22A/5KBuSwiBQQEb+SEsom5li8DfA2CQUipaKeUOYCqAddZsQCk1USm1pKqqygblkd4FensgK6oPNnL0TKJu1dVumSsB7AQQr5QqU0o9IiJGAE8A+BTAYQCrReSQNdsVkY9EZI6fn19XyqMebGxyKI78cAnHK2q0LoXIYXS1l840EeknIm4iEi4iS69MXy8icSIyUER+3z2lkjMZm9QXANhbh6gbsXsk6VJEn95I6ueLjQVsxyfqLgx80q2xyX2Rf/w8zlVf1roUIoegy8DnTVsCmh+KIgJ8xpu3RN1Cl4HPm7YEAIn9fBAe0Iu9dYi6iavWBRBZopTCuORQvLOzFNOW7DK7zH0Z4bg/I9zsPCK6li6v8ImumjZsAIZF90GTSW74Kj5XjYWfFsJk4hAMRB2hyyt8pdREABNjY2O1LoU0Fhvijfdm3Wx23roD3+MXK7/G7tJK3BwTaOfKiHoeXV7hsw2fOuL2xBD0cnPBh/s5cgdRR+gy8Ik6ore7K8Ym98X6b0+jwcgHphC1h4FPPdrdqWGoqmvEdj4Hl6hdDHzq0W4dFAz/3m5s1iHqAAY+9WjurgbcmdIPmwrOoLbBqHU5RLqmy8DnJ23JGnenhqGusQmb+AEtojbpMvDZS4esMSyqD0J9PfHRATbrELVFl4FPZA2DQWFiaj9sLSzH+ZoGrcsh0i0GPjmESWn9YTQJNhzkcMpEljDwySEkh/kiJsgL6w6c0roUIt1i4JNDUErh7rQwfFVSiR+q6rUuh0iXdBn47KVDnXF3ahhEgH9/w5u3ROboMvDZS4c6IybYGyn9/fghLCILdBn4RJ01KS0M356qQtHZSzA2ma75EuEwyuTcdDk8MlFn5Q4Jw+/XH8btf952w7zRccFYNjMLBoPSoDIi7THwyaGE+nni9WnpKC6vvmb6qQt1WLXnJNYfPI3cIWEaVUekLQY+OZy7hvS7YVqTSbD/5AUs/LQQ45JD4ebC1kxyPvypJ6fgYlCYPz4BpRW1WLXnpNblEGmCgU9OIzs+GDdF98FfPjuGmsscWZOcjy4Dn/3wyRaUUpg/IQHnqi9j6Y4SrcshsjtdBj774ZOtpA8IwPjkUPzfz79DRfVlrcshsitdBj6RLT07Ph71RhNezyvSuhQiu2Lgk9MZGOyNBzLD8b+7juNkZa3W5RDZDQOfnNKTt8XBxaCw6LOjWpdCZDcMfHJKoX6eeHBYJNbt/x6nLtRpXQ6RXTDwyWnNujUaAPD37cUaV0JkHwx8clph/r1wd1oYVu0+yUcjklNg4JNTmzt6IOoam/DOzuNal0Jkcwx8cmpxfX1wW0IIVuwsRV1Dk9blENkUA5+c3tzsgaisacDqfI6xQ45Nl4HPoRXInrKi+iAjMgBvbS+GscmkdTlENqPLwOfQCmRvc0cPRNn5Onz87WmtSyGyGV0GPpG93ZYQgtgQb/x163eo5kia5KAY+EQADAaF/7x9EI78cAmj/piHt7YVo76RN3HJsfCJV0RX5A4JQ3hAb/xpYyF+v/4w/r6jGLNvjUGYf6+WZUJ8PJAZ1UfDKok6T4mI1jVYlJmZKfn5+VqXQU5o53cVWLixEHuPn79h3r9/PhKD+/P+EumTUmqviGSancfAJzJPRFBaUYsGY3PPncYmE6Yv/QpDwv3xzs+GaVwdkXltBT6bdIgsUEohOsjrmmnzRg/EyxuOYFdxBW6OCdSoMqLO4U1bIivMGB6Fvr4e+OMnR6Dnv46JzGHgE1nB080Fv7htEPaduIDNh8+aXab6shE7jp3jLwTSHQY+kZUeyIxAVGBvLNxYCJPpx1A3mQTv55/EmIVbMX3pV/hw//caVkl0I7bhE1nJzcWAp+6Iw5Or9uPp9w8gxNcDQHPPnm/KqjB0gD/69HbHHz85gnHJoejl7qJxxUTNGPhEnTBxSPM4+utbDcUQ7OOBRVNSMSm1P3aXVmLqkl1YuqMYT+QM0rBSoh8x8Ik6wWBQWDnnZovzb44JxLjkvnhz63d4IDMCIb6edqyOyDxdtuFztExyBL+ekIjGJhP+tJEPSid90GXgc7RMcgRRQV6YcUsUVu89icIfLmldDpE+A5/IUTyREwsvd1cs3nxM61KIGPhEtuTf2x0Pj4jCx9+e5lU+aY6BT2Rjj4yMhrcHr/JJewx8IhvjVT7pBbtlEtnBIyOjseyLUrz0cQHuzwi3uFxGZADCA3rbsTJyJgx8Ijvw7+2OWbdG47XPjmH7sXMWlwvy9sDWZ7Ph7cG3JnU//lQR2ckvcgbhP9L6w2RhULWSczV4ZEU+/rb1OzwzLt7O1ZEzYOAT2YnBoBB13fj6rcUEe2NSWhje2l6MaTcNQP9Wj1Yk6g68aUukI78anwAAePWTIxpXQo6IgU+kI/39e2HWrdFYu/97bD9WrnU55GAY+EQ688SYQRgU4o2n/nkA5Zcua10OORAGPpHO9HJ3wesPpuNSfSN+uXo/KmsacKG2gU/Qoi5j4BPpUHyoD16YmIztx84h/XebkPbfmzD7nXzUNzZpXRr1YOylQ6RT04ZFINDbHacv1OGHi5fxt8+/w+x38vHWTzPh6da5p2g1NplQ29AED1dDp7dBPRcDn0inlFIYlxza8jom2Avz/983SPvvjejn1ws+nq5QAEJ8PTGkvx8CvNzhYlAYPjAQkYHN3T+r6hrx3lfHce5SA85eqsfWwnJUXzbC1aAwd/RAzLo1Gl4ernBz6dgf+yICkzT/KwBE0PK5giaToEkEJlPzMma/p64cECfj39sNSnXvEVN6bhfMzMyU/Px8rcsg0o1tR8ux7Wg5TlfVo6bBCBHg5PlaFJfXXLPckHA/hPp6Yk9pJc7XNsLHwxW9PVwwJj4EsSHeOHiqCmtbPWTd3cUAg6E5wJuDXK75v6UAJ9s5+tIEuLta3+qulNorIpnm5vEKn6gHGRUXjFFxwTdMr75sRG2DETWXm7Dh4GlsO1qO4xW1SI3wxzNj4zG4/40PE5o6bAAOfX8RtZeNqL7yy0Oh+S8Lpa7+HzAo1XxlrhQMClD4cb7B0HwFqhTgohRcDAqGq8tdd3Vq6eJSwCt/c1wM3X9UeIVPRORA2rrCZy8dIiInwcAnInISdgt8pVSMUmqpUmqNvfZJREQ/6lDgK6XeVkqdVUodvG76eKVUoVKqSCn1XFvbEJFiEXmkK8USEVHndbSXznIArwN45+oEpZQLgDcA3AGgDMAepdQ6AC4AXr5u/Z+JyNkuV0tERJ3WocAXkW1KqajrJg8DUCQixQCglFoFYJKIvAwgt7MFKaXmAJgDAAMGDOjsZoiI6DpdacPvD+Bkq9dlV6aZpZQKVEr9DcBQpdSvLS0nIktEJFNEMoODb+xvTEREndOVD16Z+1SAxU79IlIBYG4X9kdERF3QlcAvAxDR6nU4gO8tLNspe/fuPaeUOn7lpR+AKjOLWc7CesQAAAJvSURBVJoeBMDy06K1Y6lePWzb2vU7unxHlmtrmc7M4/m3/fo8/9azx/mPtLhE85gZ7X8BiAJwsNVrVwDFAKIBuAM4ACC5o9uz9gvAEiun59uqFlt8H3rYtrXrd3T5jizX1jKdmcfzz/Ov9bnW4/nvaLfMlQB2AohXSpUppR4RESOAJwB8CuAwgNUicqgj2+ukj6ycrle2rLer27Z2/Y4u35Hl2lqms/P0iOff+mV4/rtp27oeS6crlFL5YmE8CXJ8PP/OjeffPEceWmGJ1gWQpnj+nRvPvxkOe4VPRETXcuQrfCIiaoWBT0TkJBj4REROwikCXynlpZRaoZR6Syn1E63rIfvi0NzOTSn1H1fe+x8qpcZqXY+WemzgWzlk870A1ojIbAB3271Y6nbWnH/h0NwOx8rzv/bKe38mgCkalKsbPTbw0Txk8/jWE1oN2TwBQBKAaUqpJDQP+3B1oLcmO9ZItrMcHT//5HiWw/rz/19X5jutHhv4IrINQOV1k1uGbBaRBgCrAExC87g/4VeW6bHfM/3IyvNPDsaa86+avQJgg4jss3eteuJo4WdpyOYPANynlPoret5HsanjzJ7/jg7NTT2epff/zwHcDuB+pZRTj9jbldEy9cjskM0iUgPgYXsXQ3Zn6fxzaG7nYOn8Lwaw2N7F6JGjXeHbfMhm0jWef+fG898ORwv8PQAGKaWilVLuAKYCWKdxTWQ/PP/Ojee/HT028HUyZDNphOffufH8dw4HTyMichI99gqfiIisw8AnInISDHwiIifBwCcichIMfCIiJ8HAJyJyEgx8IiInwcAnInISDHwiIifx/wEOEnmuh6jLDgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(posterior_variance, label=\"evolution of var of approx posterior\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 364,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unnormalized(x, mu, prior_var, clutter_var, var):\n",
    "    n=len(x)\n",
    "    #exponent= -1*(mu*mu)/(2*prior_var) + mu/clutter_var*(sum(x)) - n*mu*mu/(2*clutter_var);\n",
    "    first=np.exp(-1*(mu*mu)/(2*prior_var))\n",
    "    res=1;\n",
    "    for i in range(n):\n",
    "        factor=(2*np.pi*var)**(-1/2)*np.exp(-0.5*(x[i]-mu)**2) + (2*np.pi*clutter_var)**(-1/2)*np.exp(-x[i]*x[i]/(2*clutter_var))\n",
    "        res=res*factor\n",
    "        \n",
    "    return first*res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in true_divide\n",
      "  import sys\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x124c0a358>"
      ]
     },
     "execution_count": 365,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIEAAAI/CAYAAADgJsn+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZRdZZkv/u/OTAaqQgYGoQ14VYZMZCINZEIZVEQFEdGL0DSC0Pxa7ZZ16WsztvfKalmXvnhtFaTFAewADtBqX6c2QAAlCSYOAREkSoQmY1VIQsbavz8O1E1ChkqoqjPU57NWrdQ5e+93P+ecWquKL8/7vkVZlgEAAACgsfWqdgEAAAAAdD0hEAAAAEAPIAQCAAAA6AGEQAAAAAA9gBAIAAAAoAcQAgEAAAD0AH2qdePhw4eXo0aNqtbtAQCoZ7/9beXfN7+5unUAQI1ZsGDBirIsR+zsWNVCoFGjRmX+/PnVuj0AAPVs5szKv3PmVLMKAKg5RVH8YVfHTAcDAAAA6AGEQAAAAAA9gBAIAAAAoAeo2ppAO7N58+YsXbo0GzZsqHYpsM8GDBiQQw89NH379q12KQAAANCupkKgpUuXZsiQIRk1alSKoqh2ObDXyrLMypUrs3Tp0hx++OHVLgcAAADa1dR0sA0bNmTYsGECIOpWURQZNmyYbjYAAABqTk2FQEkEQNQ9P8MAAADUopoLgdh7zz33XN773vcmSRYuXJjvf//77ceuvfba3HjjjdUqbTv33Xdfbrjhht2es2TJktx5553dUs93vvOdLF68uP3x1VdfnR//+Mfdcm8AAADobkKgbrBly5YuHf+QQw7JPffck+TVIVAtOeOMM3LllVfu9pzODoG2bt26y2M7hkDXX3993vrWt3bavQEAAKCWCIF28O53vzsTJ07MMccck1tuuaX9+cGDB+dv//ZvM2HChLzlLW/J8uXLkyQzZ87Mxz72sRx//PEZPXp0Hn300SSVDpyLL744p5xySj70oQ9lw4YN+Yu/+IuMGTMmxx57bH76058mSf7X//pfufDCC5Mkv/rVrzJ69OisX79+u5re/va355e//GWS5Nhjj83111+fJLnqqqvypS99KUuWLMno0aOzadOmXH311Zk9e3bGjx+f2bNnJ0kWL16cmTNn5ogjjsjNN9+809c9ePDgfPKTn8y4ceMyderUvPDCC0mSCy64oD1geuW8JJkzZ05mzJiR973vfXnTm96UK6+8MnfccUemTJmSMWPG5Omnn37VPW6//fZcfvnl7eP+9V//dY4//vgcccQR7fe48sor8+CDD2b8+PG56aabsnXr1lxxxRWZPHlyxo4dmy9+8YtJkra2tlx22WU55phjcvrpp+ftb397+xijRo3K9ddfnxNPPDF33313br311kyePDnjxo3LWWedlfXr1+fhhx/OfffdlyuuuCLjx4/P008/vd1r/clPfpJjjz02Y8aMyYUXXpiNGze2j33NNddkwoQJGTNmTJ544old/zABAABADREC7eBf/uVfsmDBgsyfPz8333xzVq5cmSRZt25dJkyYkMceeywzZszIdddd137NunXr8vDDD+ef//mf2wOdJFmwYEHuvffe3Hnnnfnc5z6XpBL0fOMb38j555+fDRs25GMf+1ieeuqpfPvb385f/MVf5Itf/GIGDhy4XU3Tp0/Pgw8+mDVr1qRPnz556KGHkiRz587NtGnT2s/r169frr/++pxzzjlZuHBhzjnnnCTJE088kR/84Ad59NFHc91112Xz5s2vet3r1q3L1KlTs2jRokyfPj233nrrHt+rRYsW5X//7/+dX/3qV/na176WJ598Mo8++mguuuiifPazn93j9c8//3zmzp2b7373u+0dQjfccEOmTZuWhQsX5uMf/3huu+22NDU1Zd68eZk3b15uvfXWPPPMM/nWt76VJUuW5Fe/+lW+9KUv5ZFHHtlu7AEDBmTu3Ll5//vfnzPPPDPz5s3LokWLctRRR+W2227L8ccfnzPOOCOf+cxnsnDhwrzhDW9ov3bDhg254IILMnv27PzqV7/Kli1b8vnPf779+PDhw/PYY4/l0ksvrZmpdgAAALAnQqAd3Hzzze3dMM8++2x+97vfJUl69erVHqr81//6XzN37tz2a84999wklbBmzZo1aWlpSVKZ/rTffvslqQQ25513XpLkyCOPzOtf//o8+eST6dWrV26//facd955mTFjRk444YRX1TRt2rQ88MADmTt3bt7xjndk7dq1Wb9+fZYsWZI3v/nNe3xN73jHO9K/f/8MHz48I0eObO/y2Va/fv1y+umnJ0kmTpyYJUuW7HHcyZMn5+CDD07//v3zhje8IaecckqSZMyYMR26/t3vfnd69eqVo48+eqc1JckPf/jDfPWrX8348eNz3HHHZeXKlfnd736XuXPn5uyzz06vXr1y0EEHZdasWdtd98pnlSS//vWvM23atIwZMyZ33HFHfvOb3+y2rt/+9rc5/PDD86Y3vSlJcv755+eBBx5oP37mmWcm6fj7BAAAALWgT7UL2J2L/+3iPPfic5023iFDDskt77xll8fnzJmTH//4x3nkkUcycODAzJw5c5dbfW+7A9SOu0G98njQoEHtz5Vlucv7/u53v8vgwYPz3HM7f62TJ0/O/Pnzc8QRR+Tkk0/OihUrcuutt2bixIm7HHNb/fv3b/++d+/eO12jqG/fvu11b3tOnz590tbW1v4aNm3atNNxe/Xq1f64V69eHVoHadvrd/X+lGWZz372szn11FO3e/573/vebsfe9r2/4IIL8p3vfCfjxo3L7bffnjlz5uz22t19VtvWvav3EgAAAGpRTYdAuwtsukJra2uGDh2agQMH5oknnsjPfvaz9mNtbW2555578v73vz933nlnTjzxxPZjs2fPzqxZszJ37tw0NTWlqanpVWNPnz49d9xxR0466aQ8+eST+eMf/5g3v/nNaW1tzUc/+tE88MADufzyy3PPPfe07/T1in79+uWwww7LXXfdlauuuirLly/PJz7xiXziE5941X2GDBmSF198sdPek1GjRmXBggV53/vel3vvvXenU8k60471n3rqqfn85z+fk046KX379s2TTz6Z173udTnxxBPzla98Jeeff36WL1+eOXPm5AMf+MBOx3zxxRdz8MEHZ/Pmzbnjjjvyute9bqf3esWRRx6ZJUuW5Kmnnsp/+S//JV/72tcyY8aMrnnBAAAA0E1MB9vGaaedli1btmTs2LG56qqrMnXq1PZjgwYNym9+85tMnDgx//Ef/5Grr766/djQoUNz/PHH5yMf+Uhuu+22nY592WWXZevWrRkzZkzOOeec3H777enfv38+/vGP57LLLsub3vSm3HbbbbnyyiuzbNmyV10/bdq0HHjggRk4cGCmTZuWpUuXbrce0CtmzZqVxYsXb7cw9Gvx4Q9/OPfff3+mTJmSn//859t12HSFsWPHpk+fPhk3blxuuummXHTRRTn66KMzYcKEjB49Opdcckm2bNmSs846K4ceemj7c8cdd9xOw7ck+Yd/+Iccd9xxOfnkk3PkkUe2P//+978/n/nMZ3Lsscdut5D1gAED8uUvfzlnn312xowZk169euUjH/lIl75uAAAA6GrFnqa+dJVJkyaV8+fP3+65xx9/PEcddVRV6tmTwYMHZ+3ata96fubMmbnxxhszadKkKlTVs61duzaDBw/OypUrM2XKlDz00EM56KCDql1Wktr+WQaAhjBzZuXfPUzzBoCepiiKBWVZ7jSkqOnpYLA7p59+elpaWrJp06ZcddVVNRMAAQAAQC0SAnXQzrqAkuxxkWG6jvceAAAAOs6aQAAAAAA9gBAIAAAAoAcQAgEAAAD0AEIgAAAAgB5ACLSNlpaW/PM//3O1y+gWS5YsyZ133tn+eP78+fnrv/7rqtRy7bXX5sYbb6zKvQGAxvT888lXvpKsX1/tSgCgdgiBtrG7EGjr1q3dXM1rt2XLll0e2zEEmjRpUm6++ebuKAsAoEtt2JB88IPJ2rXJpZdWuxoAqB1CoG1ceeWVefrppzN+/PhcccUVmTNnTmbNmpUPfOADGTNmTJYsWZLRo0e3n3/jjTfm2muvTZI8/fTTOe200zJx4sRMmzYtTzzxxKvGX758eU4++eRMmDAhl1xySV7/+tdnxYoVSZKvf/3rmTJlSsaPH59LLrmkPXQaPHhwPvnJT2bcuHGZOnVqXnjhhfaxzjrrrEyePDmTJ0/OQw89lKTSVXPxxRfnlFNOyYc+9KEsWbIk06ZNy4QJEzJhwoQ8/PDD7a/1wQcfzPjx43PTTTdlzpw5Of3009vHuPDCCzNz5swcccQR24VD//AP/5AjjzwyJ598cs4999xXdfC0trZm1KhRaWtrS5KsX78+hx12WDZv3pxbb701kydPzrhx43LWWWdl/U7+19zMmTMzf/78JMmKFSsyatSoJJUQ7oorrsjkyZMzduzYfPGLX9yLTxYA6Em++tXkQx9K/uqvkoEDk4ULq10RANQGIdA2brjhhrzhDW/IwoUL85nPfCZJ8uijj+Z//I//kcWLF+/22osvvjif/exns2DBgtx444257LLLXnXOddddl5NOOimPPfZY3vOe9+SPf/xjkuTxxx/P7Nmz89BDD2XhwoXp3bt37rjjjiTJunXrMnXq1CxatCjTp0/PrbfemiT56Ec/mo9//OOZN29evvnNb+aiiy5qv8+CBQty77335s4778zIkSPzox/9KI899lhmz57dPuXrhhtuyLRp07Jw4cJ8/OMff1WtTzzxRH7wgx/k0UcfzXXXXZfNmzdn/vz5+eY3v5lf/OIX+da3vtUe1myrqakp48aNy/33358k+bd/+7eceuqp6du3b84888zMmzcvixYtylFHHZXbbrttj5/JK2677bY0NTVl3rx5mTdvXm699dY888wzHb4eAOg5vvnN5JxzKt9/5CPJl79c3XoAoFb0qXYBu/Sxj3X+/7YZPz75p3/aq0umTJmSww8/fLfnrF27Ng8//HDOPvvs9uc2btz4qvPmzp2bb3/720mS0047LUOHDk2S/OQnP8mCBQsyefLkJMlLL72UkSNHJkn69evX3qEzceLE/OhHP0qS/PjHP94umFqzZk1efPHFJMkZZ5yR/fbbL0myefPmXH755e3h0pNPPtmh1/2Od7wj/fv3T//+/TNy5Mi88MILmTt3bt71rne1j/3Od75zp9eec845mT17dmbNmpV//dd/bQ/Efv3rX+fv//7v09LSkrVr1+bUU0/tUC1J8sMf/jC//OUvc8899ySpdBz97ne/2+NnAwD0LEuXJocckrz850rGjUsWL07KMimK6tYGANVWuyFQjRg0aFD793369Gmf5pQkGzZsSJK0tbWlubk5C/cQWpVlucvnzz///Hz6059+1bG+ffumePkvlt69e7ev89PW1pZHHnmkPZDZVc033XRTDjzwwCxatChtbW0ZMGDAbmt8Rf/+/du/f+W+u6p/R2eccUb+7u/+LqtWrcqCBQty0kknJUkuuOCCfOc738m4ceNy++23Z86cOa+6dtv3+JX3N6m8R5/97Gf3KjgCAHqeH/0oOfnk7Z87+ujkiSeSo46qTk0AUCtqdzrYP/1TMmdO537toQtoyJAh7d00O3PggQdm2bJlWblyZTZu3Jjvfve7SZL9998/hx9+eO6+++4klcBi0aJFr7r+xBNPzF133ZWk0tmyevXqJMlb3vKW3HPPPVm2bFmSZNWqVfnDH/6w21pPOeWU/J//83/aH+8qgGptbc3BBx+cXr165Wtf+1r7WkN7eq07c+KJJ+bf/u3fsmHDhqxduzbf+973dnre4MGDM2XKlHz0ox/N6aefnt69eydJXnzxxRx88MHZvHlz+3S3HY0aNSoLFixIkvaunyQ59dRT8/nPfz6bN29Okjz55JNZt27dXtUPADS+H/84eetbt3/urW+thEMA0NPVbghUBcOGDcsJJ5yQ0aNH54orrnjV8b59++bqq6/Occcdl9NPPz1HHnlk+7E77rgjt912W8aNG5djjjkm995776uuv+aaa/LDH/4wEyZMyL//+7/n4IMPzpAhQ3L00UfnU5/6VE455ZSMHTs2J598cp5//vnd1nrzzTdn/vz5GTt2bI4++uh84Qtf2Ol5l112Wb7yla9k6tSpefLJJ9u7hMaOHZs+ffpk3Lhxuemmmzr0/kyePDlnnHFGxo0blzPPPDOTJk1KU1PTTs8955xz8vWvfz3nvDIhP5VFpY877ricfPLJ27132/rEJz6Rz3/+8zn++OPbF81OkosuuihHH310JkyYkNGjR+eSSy7Z7e5nAEDPtHJl8vKs+nYnnJD87GfVqQcAaknR0Sk+nW3SpEnljgsLP/744zmqgft0N27cmN69e6dPnz555JFHcumll+5xClmtWbt2bQYPHpz169dn+vTpueWWWzJhwoRql1VzGv1nGQCqbubMyr/bTC9fvjz5m79Jvva1V59+2mnJ//2/3VIZAFRVURQLyrKctLNj1gTqRn/84x/zvve9L21tbenXr1/7Tl/15OKLL87ixYuzYcOGnH/++QIgAKBmLFiQTJy482NDhyarV1f+BYCeSgjUjd74xjfmF7/4RbXLeE3uvPPOapcAALBTCxYk06fv/NikSZXjO64XBAA9iTWBAABoCIsWJePH7/zY2LHJr3/dvfUAQK0RAgEA0BDWrUuGDNn5saOOShYv7t56AKDWCIEAAKh7mzcnvXvv+vjrXpcsXdp99QBALRICAQBQ9555JjniiF0fL4rKv1XaGBcAaoIQqMbMmTMnDz/8cJfe43/+z/+5T9dddNFFWdxNfdQLFy7M97///fbH9913X2644YZOGfuf/umfsn79+k4ZCwCoDU88kRx55O7POeig5IUXuqceAKhFQqBOsnXr1u0el2WZtra2vR6nVkOgrVu35ktf+lKOPvrovbpmX+0YAp1xxhm58sor93m8bQmBAKDxdCQEOvLI5PHHu6ceAKhFQqAdfPWrX83YsWMzbty4nHfeeUmSCy64IPfcc0/7OYMHD05SCWxmzZqVD3zgAxkzZkyWLFmSo446KpdddlkmTJiQZ599Nj/84Q/z53/+55kwYULOPvvsrF27NkkyatSoXHPNNZkwYULGjBmTJ554IkuWLMkXvvCF3HTTTRk/fnwefPDB7Wq79tprc9555+Wkk07KG9/4xtx6661JKoHTFVdckdGjR2fMmDGZPXt2kuT555/P9OnTM378+IwePToPPvhgrrzyyrz00ksZP358PvjBDyZJvv71r2fKlCkZP358LrnkkvbwZvDgwbn66qtz3HHH5ZFHHsnMmTMzf/78JMk3vvGNjBkzJqNHj85/+2//bbv3ZttrtvX000/ntNNOy8SJEzNt2rQ88cQTSZK77747o0ePzrhx4zJ9+vRs2rQpV199dWbPnp3x48dn9uzZuf3223P55Ze3fx6XXnppZs2alSOOOCL3339/Lrzwwhx11FG54IIL2u936aWXZtKkSTnmmGNyzTXXJEluvvnmPPfcc5k1a1ZmzZqVJLv8jACA+tGREOgNb0h+//vuqQcAalJZllX5mjhxYrmjxYsXv+q57vTrX/+6fNOb3lQuX768LMuyXLlyZVmWZXn++eeXd999d/t5gwYNKsuyLH/605+WAwcOLH//+9+XZVmWzzzzTFkURfnII4+UZVmWy5cvL6dNm1auXbu2LMuyvOGGG8rrrruuLMuyfP3rX1/efPPNZVmW5ec+97nyL//yL8uyLMtrrrmm/MxnPrPT+q655ppy7Nix5fr168vly5eXhx56aPmnP/2pvOeee8q3vvWt5ZYtW8r//M//LA877LDyueeeK2+88cbyU5/6VFmWZblly5ZyzZo129VflpX3/PTTTy83bdpUlmVZXnrppeVXvvKVsizLMkk5e/bs9nNnzJhRzps3r/zTn/5UHnbYYeWyZcvKzZs3l7NmzSq//e1v7/SabZ100knlk08+WZZlWf7sZz8rZ82aVZZlWY4ePbpcunRpWZZluXr16rIsy/LLX/5y+Vd/9Vft1277+Pzzzy/POeecsq2trfzOd75TDhkypPzlL39Zbt26tZwwYUL5i1/8YrvPb8uWLeWMGTPKRYsWtb/3r3zGu/uMXotq/ywDQMObMaPy9bK3va0s29p2f8kvflGW//2/d2lVAFB1SeaXu8hi+lQ5g6op//Ef/5H3vve9GT58eJLkgAMO2OM1U6ZMyeGHH97++PWvf32mTp2aJPnZz36WxYsX54QTTkiSbNq0KX/+53/efu6ZZ56ZJJk4cWK+9a1vdajGd73rXdlvv/2y3377ZdasWXn00Uczd+7cnHvuuendu3cOPPDAzJgxI/PmzcvkyZNz4YUXZvPmzXn3u9+d8ePHv2q8n/zkJ1mwYEEmT56cJHnppZcycuTIJEnv3r1z1llnveqaefPmZebMmRkxYkSS5IMf/GAeeOCBvPvd797lNWvXrs3DDz+cs88+u/25jRs3JklOOOGEXHDBBXnf+97X/p7syTvf+c4URZExY8bkwAMPzJgxY5IkxxxzTJYsWZLx48fnrrvuyi233JItW7bk+eefz+LFizN27NjtxtnTZwQA1I9XFn/elcMP1wkEQM9W0yHQxRcnzz3XeeMdckhyyy27Pl6WZYqd/PXQp0+f9vV9yrLMpk2b2o8NGjRou3O3fVyWZU4++eR84xvf2On9+vfvn6QStmzZsqVDr2HH+oqiSLmLbS6mT5+eBx54IN/73vdy3nnn5YorrsiHPvSh7c4pyzLnn39+Pv3pT7/q+gEDBqT3TvZa3dX9dndNW1tbmpubs3Dhwlcd+8IXvpCf//zn+d73vpfx48fv9JwdvfLe9erVq/37Vx5v2bIlzzzzTG688cbMmzcvQ4cOzQUXXJANGzbs9LXs7jMCAGrfpk1Jv357Pq+pKWlt7fp6AKBW1fSaQLfcknz3u533tbsAKEne8pa35K677srKlSuTJKtWrUpSWb9nwYIFSZJ77703mzdv7lD9U6dOzUMPPZSnnnoqSbJ+/fo8+eSTu71myJAhefHFF3d5/N57782GDRuycuXKzJkzJ5MnT8706dMze/bsbN26NcuXL88DDzyQKVOm5A9/+ENGjhyZD3/4w/nLv/zLPPbYY0mSvn37tr+Gt7zlLbnnnnuybNmy9tf8hz/8Ybc1Hnfccbn//vuzYsWKbN26Nd/4xjcyY8aM3V6z//775/DDD8/dd9+dpBK+LFq0KEllraDjjjsu119/fYYPH55nn312j+/DnqxZsyaDBg1KU1NTXnjhhfz7v/97+7Ftx96XzwgAqC3PPpv82Z9VuwoAqH01HQJ1t2OOOSaf/OQnM2PGjIwbNy5/8zd/kyT58Ic/nPvvvz9TpkzJz3/+81d1/+zKiBEjcvvtt+fcc8/N2LFjM3Xq1PbFkHflne98Z7797W/vdGHopDL97B3veEemTp2aq666Koccckje8573tC9mfdJJJ+Uf//Efc9BBB2XOnDkZP358jj322Hzzm9/MRz/60STJxRdfnLFjx+aDH/xgjj766HzqU5/KKaeckrFjx+bkk0/O888/v9saDz744Hz605/OrFmzMm7cuEyYMCHvete79vh+3HHHHbntttsybty4HHPMMbn33nuTJFdccUX7ItPTp0/PuHHjMmvWrCxevLh9Yei9NW7cuBx77LE55phjcuGFF7ZP93rl9b/tbW/LrFmz9ukzAgBqy5IlyahRHTtXNxAAPVmxu6k9XWnSpEnlKztNveLxxx/PUUcdVZV66sG1116bwYMH5xOf+ES1S2EP/CwDQBebObPy75w5+dKXkmHDkve8Z8+XffKTyXvfmxx7bJdWBwBVUxTFgrIsJ+3smE4gAADq2t50Ar3+9ckf/9iV1QBA7arphaHZ3rXXXlvtEgAAas7ehECHHZY8/XRXVgMAtUsnEAAAdW316qS5uWPnHnposnRp19YDALWq5kKgaq1RBJ3FzzAAdL+i6Nh5hx1W2U0MAHqimgqBBgwYkJUrV/qPaOpWWZZZuXJlBgwYUO1SAKBH2LIl6bMXCxw0NSUtLV1XDwDUsppaE+jQQw/N0qVLs3z58mqXAvtswIABOfTQQ6tdBgD0CC+8kBx4YMfP72jHEAA0opoKgfr27ZvDDz+82mUAAFAn/vSn5HWv27triiJpa0t61VRPPAB0Pb/6AACoW889t/ch0MiRicZzAHoiIRAAAHXrT39KDjlk766xQxgAPZUQCACAurUvnUCHHSYEAqBnEgIBAFC39rUTyDbxAPREQiAAAOrWihXJ8OF7d80hh1Q6iACgpxECAQBQ1/Z22/eDD07+8z+7phYAqGVCIAAAepRhwyodRADQ0wiBAACoS+u2DsjAgXt/Xe/eSVtb59cDALVOCAQAQF16buOwvV4UGgB6MiEQAAB16blNw/c5BOrVK9m6tXPrAYBaJwQCAKAuvbBpaA46aN+uHTHCukAA9DxCIAAA6tILm4dm5Mh9u/bgg5Pnn+/cegCg1gmBAACoS8s2Dc2BB+7btQcdZJt4AHoeIRAAAHXphc1CIADYG0IgAADq0rJNQzNixL5dazoYAD2REAgAgLq0peydvn337VqdQAD0REIgAADqUpFyn68VAgHQEwmBAACoS2WKfb520KBk3bpOLAYA6oAQCACAurNu64AM6v1StcsAgLoiBAIAoO68sGloDuy3utplAEBdEQIBAFB3Xtg0NAf2fW0hUJ8+yZYtnVQQANQBIRAAAHVn2eahGfkaO4GGDUtWruykggCgDgiBAACoO53RCTRiRLJ8eScVBAB1QAgEAEDdeWHza18TSAgEQE8jBAIAoO68sOmATgmBli3rpIIAoA7sMQQqiuJfiqJYVhTFr3dxvCiK4uaiKJ4qiuKXRVFM6PwyAQDg/1m2uTkjTQcDgL3SkU6g25Octpvjb0vyxpe/Lk7y+ddeFgAA7Nr6rQMysPfG1zSGEAiAnmaPIVBZlg8kWbWbU96V5Ktlxc+SNBdFcXBnFQgAADsqUr7mMYRAAPQ0nbEm0OuSPLvN41TjO8gAACAASURBVKUvPwcAADVLCARAT9MZIVCxk+d2+r9miqK4uCiK+UVRzF/uNy4AAPugfO1NQEmSgQOTl17qnLEAoB50Rgi0NMlh2zw+NMlzOzuxLMtbyrKcVJblpBEjRnTCrQEA6GlaW5OmPuuqXQYA1J3OCIHuS/Khl3cJm5qktSzL5zthXAAAeJWVK5PhfVurXQYA1J0+ezqhKIpvJJmZZHhRFEuTXJOkb5KUZfmFJN9P8vYkTyVZn+QvuqpYAABYsSIZ1kkhUK9eydatSe/enTIcANS0PYZAZVmeu4fjZZK/6rSKAABgN1as6LxOoAMOSFatqiwSDQCNrjOmgwEAQLfpzOlgdggDoCcRAgEAUFdWrEiG9VnTKWMJgQDoSYRAAADUlc6cDiYEAqAnEQIBAFBXumJNIADoCYRAAADUlZUrk2F9O2c62LBhlfEAoCcQAgEAUFc2bkz699rcKWMJgQDoSYRAAAD0WMOGmQ4GQM8hBAIAoMcaOlQIBEDPIQQCAKBulGXnjte3b7JlS+eOCQC1SggEAEDdaG1NmpurXQUA1CchEAAAdWPlyso6PgDA3hMCAQBQN1asSIYP7/xxO3uaGQDUIiEQAAB1oytCoP33T158sXPHBIBaJAQCAKBudMV0sAMOqIwLAI1OCAQAQN3oik6gYcOEQAD0DEIgAADqhhAIAPadEAgAgLqxcmVl+lZnGjYsWbWqc8cEgFokBAIAoG6sWtX5IZA1gQDoKYRAAADUjZdeSvbbr3PHNB0MgJ5CCAQAQF0pis4dz3QwAHoKIRAAAD2a6WAA9BRCIAAAerSmpqS1tdpVAEDXEwIBAFAXtm5NenXBX6+9eiVl2fnjAkCtEQIBAFAXWlqSoUOrXQUA1C8hEAAAdWH1aiEQALwWQiAAAOrCqlVdFwL17l2ZbgYAjUwIBABAXVi9urKTV1cYOrQyPgA0MiEQAAB1oSs7gWwTD0BPIAQCAKAudHUnUEtL14wNALVCCAQAQF3oyoWhTQcDoCcQAgEAUBe6cjpYc7MQCIDGJwQCAKAuWBgaAF4bIRAAAHWhKzuBhEAA9ARCIAAA6sKWLUnfvl0ztoWhAegJhEAAAPR4OoEA6AmEQAAA9HgWhgagJxACAQDQ4w0YkGzcWO0qAKBrCYEAAKh5Gzcm/fpVuwoAqG9CIAAAat7q1V23MxgA9BRCIAAAat7q1ckBB1S7CgCob0IgAABq3qpVXd8J1Lt3snVr194DAKpJCAQAQM3rjk6gpqaktbVr7wEA1SQEAgCg5nXHmkC2iQeg0QmBAACoed0xHWzoUCEQAI1NCAQAQM3rjulgQiAAGp0QCACAmtcd08GEQAA0OiEQAAA1b9Wq7ukEamnp2nsAQDUJgQAAqHmtrZXdu7qSTiAAGp0QCACAmtfWlvTq4r9c7Q4GQKMTAgEAQHQCAdD4hEAAABAhEACNTwgEAABJBg5M1q+vdhUA0HWEQAAA1LQNG5IBA7r+PkXR9fcAgGoSAgEAUNO6Y2ewVwiCAGhkQiAAAGpaS0tl567uUJbdcx8AqAYhEAAANa27O4Ha2rrnXgDQ3YRAAADUtO7sBBoyJHnxxe65FwB0NyEQAAA1raWl+zqBbBMPQCMTAgEAUNNaW7uvE2jo0EroBACNSAgEAEBN687pYDqBAGhkQiAAAGpady4MLQQCoJEJgQAAqGnd2QnU3CwEAqBxCYEAAKhp3dkJ1NxcuR8ANCIhEAAANW3NmmT//bvnXk1NFoYGoHEJgQAAqGlbtya9e3fPvXQCAdDIhEAAAPCy5madQAA0LiEQAAC8bP/9dQIB0LiEQAAA8LK+fZMtW6pdBQB0DSEQAAA1a+vWpJe/WAGgU/iVCgBAzXrxxe7bGQwAGp0QCACAmtXSUlmsGQB47YRAAADUrNbWpKmpe+9ZFElbW/feEwC6gxAIAICaVY1OoCFDkrVru/eeANAdhEAAANSsanQCNTdXwicAaDRCIAAAalY1OoGamirhEwA0GiEQAAA1q6VFJxAAdBYhEAAANau1tTqdQEIgABqREAgAgJpVjelgzc2mgwHQmDoUAhVFcVpRFL8tiuKpoiiu3MnxPyuK4qdFUfyiKIpfFkXx9s4vFQCAnsbC0ADQefYYAhVF0TvJ55K8LcnRSc4tiuLoHU77+yR3lWV5bJL3J/nnzi4UAICex8LQANB5OtIJNCXJU2VZ/r4sy01J/jXJu3Y4p0yy/8vfNyV5rvNKBACgp9qwIRkwoHvvqRMIgEbVpwPnvC7Js9s8XprkuB3OuTbJD4ui+P+SDEry1k6pDgAAupkQCIBG1ZFOoGInz5U7PD43ye1lWR6a5O1JvlYUxavGLori4qIo5hdFMX/58uV7Xy0AAHQx08EAaFQdCYGWJjlsm8eH5tXTvf4yyV1JUpblI0kGJBm+40BlWd5SluWksiwnjRgxYt8qBgCALrTffsn69dWuAgA6X0dCoHlJ3lgUxeFFUfRLZeHn+3Y4549J3pIkRVEclUoIpNUHAIB9Vu7Ye95Nip31wQNAA9hjCFSW5ZYklyf5QZLHU9kF7DdFUVxfFMUZL5/2t0k+XBTFoiTfSHJBWVbr1zYAAI1gw4ZKVw4A0Dk6sjB0yrL8fpLv7/Dc1dt8vzjJCZ1bGgAAPVk1tocHgEbWkelgAADQ7VpbK4s0AwCdQwgEAEBNqmYn0IABleloANBIhEAAANSkanYC2SYegEYkBAIAoCZVsxOoublyfwBoJEIgAABqUjU7gYRAADQiIRAAADWpmp1ApoMB0IiEQAAA1KSWFp1AANCZhEAAANSk1lZrAgFAZxICAQBQk0wHA4DOJQQCAKAmvfhiMmRIde6tEwiARiQEAgCgJrW1Jb2q9NeqTiAAGpEQCAAAdqATCIBGJAQCAIAd7L9/smZNtasAgM4lBAIAgB306lWZjgYAjUQIBABAzdmyJendu9pVAEBjEQIBAFBz1qypTMkCADqPEAgAgJrT2lpZnBkA6DxCIAAAak5LS/VDoN69k61bq1sDAHQmIRAAADWntTVpaqpuDXYIA6DRCIEAAKg5tdAJ1NxcqQMAGoUQCACAmlMLnUBCIAAajRAIAICaUwudQE1NlTAKABqFEAgAgJqjEwgAOp8QCACAmqMTCAA6nxAIAICa09KiEwgAOpsQCACAmmM6GAB0PiEQAAA1Z+PGZMCA6tZgOhgAjUYIBAAAO6ETCIBGIwQCAICd0AkEQKMRAgEAwE7071+ZlgYAjUIIBABATSnLalcAAI1JCAQAQE156aVk4MBqVwEAjUcIBABATWlpqSzKDAB0LiEQAAA1paWlsihzrTA9DYBGIQQCAKCmtLbWTifQwIHJ+vXVrgIAOocQCACAmlJL08Gam20TD0DjEAIBAFBTWltrZzpYc3MllAKARiAEAgCgptRSJ1BTk04gABqHEAgAgJqiEwgAuoYQCACAmqITCAC6hhAIAICaohMIALqGEAgAgJqiEwgAuoYQCACAmvLii8ngwdWuokInEACNRAgEAEBNKcukV438laoTCIBGUiO/XgEAoPboBAKgkQiBAABgFwYPrkxPA4BGIAQCAIBd6NWrMj0NABqBEAgAgJqxZUvSp0+1qwCAxiQEAgCgZrS2VhZjBgA6nxAIAICaIQQCgK4jBAIAoGa0tFR25KolffpUpqkBQL0TAgEAUDNaWmqvE6ipqdKhBAD1TggEAEDNaG2tvU4gIRAAjUIIBABAzajF6WDNzZW6AKDeCYEAAKgZtbgwtE4gABqFEAgAgJqhEwgAuo4QCACAmqETCAC6jhAIAICaoRMIALqOEAgAgJpRq51AQiAAGoEQCACAmrFpU9KvX7Wr2F5zs+lgADQGIRAAAOyG6WAANAohEAAA7IaFoQFoFEIgAADYjX79KtPUAKDeCYEAAKgJZZkURbWrAIDGJQQCAKAmrFuXDBpU7SoAoHEJgQAAqAm1uD08ADQSIRAAADWhpaWyE1etKstqVwAAr40QCACAmlDLIdCgQcn69dWuAgBeGyEQAAA1oZangzU3V0IqAKhnQiAAAGpCLXcCNTUJgQCof0IgAABqQi2HQM3NlU4lAKhnQiAAAGqC6WAA0LWEQAAA1IRa7gRqatIJBED9EwIBAFATajkE0gkEQCMQAgEAUBNqeTqYTiAAGoEQCACAmrB2bTJoULWr2DmdQAA0AiEQAAA1oyiqXcHO6QQCoBEIgQAAYA90AgHQCDoUAhVFcVpRFL8tiuKpoiiu3MU57yuKYnFRFL8piuLOzi0TAACqZ/DgynQ1AKhnffZ0QlEUvZN8LsnJSZYmmVcUxX1lWS7e5pw3Jvm7JCeUZbm6KIqRXVUwAAB0t6JIyrLaVQDAa9ORTqApSZ4qy/L3ZVluSvKvSd61wzkfTvK5sixXJ0lZlss6t0wAABrZxo1Jv37VrgIAGltHQqDXJXl2m8dLX35uW29K8qaiKB4qiuJnRVGc1lkFAgDQ+FpbK+vuAABdZ4/TwZLsbI+GHZth+yR5Y5KZSQ5N8mBRFKPLstxu+byiKC5OcnGS/Nmf/dleFwsAQGNqba3swAUAdJ2OdAItTXLYNo8PTfLcTs65tyzLzWVZPpPkt6mEQtspy/KWsiwnlWU5acSIEftaMwAADaalpfY7gfr0STZvrnYVALDvOhICzUvyxqIoDi+Kol+S9ye5b4dzvpNkVpIURTE8lelhv+/MQgEAaFz1EAI1N1c6lgCgXu0xBCrLckuSy5P8IMnjSe4qy/I3RVFcXxTFGS+f9oMkK4uiWJzkp0muKMtyZVcVDQBAY6mH6WBNTUIgAOpbR9YESlmW30/y/R2eu3qb78skf/PyFwAA7JWWlmT48GpXsXvNzZU6AaBedWQ6GAAAdKl6mA6mEwiAeicEAgCg6uphOphOIADqnRAIAICq0wkEAF1PCAQAQNW1tOgEAoCuJgQCAKDq1qxJ9t+/2lXsXlOTEAiA+iYEAgCg6rZuTfp0aN/a6mluNh0MgPomBAIAgA7QCQRAvRMCAQBAB1gYGoB6JwQCAIAO6Ncv2by52lUAwL4TAgEAUFVtbUkvf5UCQJfz6xYAgKqqh53BAKARCIEAAKiqlpbKzlsAQNcSAgEAUFWtrZVFl+tBUSRlWe0qAGDfCIEAAKiqeuoEGjQoWbeu2lUAwL4RAgEAUFX1FAI1N1fqBYB6JAQCAKCq6mk6WFOTEAiA+iUEAgCgquqtE6i1tdpVAMC+EQIBAFBVLS06gQCgOwiBAACoqtZWnUAA0B2EQAAAVFW9TQfTCQRAvRICAQBQVfW2MLROIADqlRAIAICq2rAhGTCg2lV0jE4gAOqZEAgAADpIJxAA9UwIBAAAHaQTCIB6JgQCAIAOGjQoWbeu2lUAwL4RAgEAUDVlWe0K9k5R1F/NAPAKIRAAAFXz0kvJwIHVrgIAegYhEAAAVdPSUllnBwDoekIgAACqRggEAN1HCAQAQNW0tla2Xa8nffsmmzZVuwoA2HtCIAAAqqYeO4GamirhFQDUGyEQAABV09JSf51Azc1CIADqkxAIAICqaWlJhg6tdhV7p7m5UjcA1BshEAAAVbN6df2FQKaDAVCvhEAAAFTN6tX1tyaQTiAA6pUQCACAqtEJBADdRwgEAEDV1GMIpBMIgHolBAIAoGo2bkwGDKh2FXtHJxAA9UoIBAAAe6G5udLBBAD1RggEAAB7wXQwAOqVEAgAAPaCEAiAeiUEAgCgKrZuTXrV4V+jvXsnbW3VrgIA9l4d/toFAKARrFlTWWQZAOgeQiAAAKqiHreHB4B6JgQCAKAqVq+urK8DAHQPIRAAAFVRz51A/fsnGzZUuwoA2DtCIAAAqqKlpX5DoKFDKyEWANQTIRAAAFVRz51ABxwgBAKg/giBAACoinpeE0gnEAD1SAgEAEBV1HMnkBAIgHokBAIAoCrqPQRataraVQDA3hECAQBQFRaGBoDuJQQCAKAq1q5NBg2qdhX7RggEQD0SAgEAUDVFUe0K9o0QCIB6JAQCAIC9JAQCoB4JgQAAYC81N1fWNAKAeiIEAgCg25Vl/U4FS5LevZO2tmpXAQB7RwgEAEC3W7s2GTy42lUAQM8iBAIAoNutXl2ZUgUAdB8hEAAA3W716sriyvWsnqezAdAzCYEAAOh2jRAC9e+fvPRStasAgI4TAgEA0O1aWuo/BLJNPAD1RggEAEC3a4ROICEQAPVGCAQAQLdrhIWhhUAA1BshEAAA3U4nEAB0PyEQAADdTggEAN1PCAQAQLezMDQAdD8hEAAA3W7NmmT//atdxWsjBAKg3giBAADodm1tSa86/0tUCARAvanzX70AAFAdBxwgBAKgvgiBAABgHzQ1Ja2t1a4CADpOCAQAQLcqy2pX0Dl6965MawOAeiEEAgCgW61blwweXO0qOkejBFoA9AxCIAAAutXKlcmwYdWuAgB6HiEQAADdatWqyqLKAED3EgIBANCtVq5snBBowIDkpZeqXQUAdIwQCACAbrVqVeNMBxs61DbxANQPIRAAAN2qkTqBhEAA1JMOhUBFUZxWFMVvi6J4qiiKK3dz3nuLoiiLopjUeSUCANBIGqkT6IADKq8HAOrBHkOgoih6J/lckrclOTrJuUVRHL2T84Yk+eskP+/sIgEAaByNtDC0EAiAetKRTqApSZ4qy/L3ZVluSvKvSd61k/P+Ick/JtnQifUBANBgGmk62LBhldcDAPWgIyHQ65I8u83jpS8/164oimOTHFaW5Xc7sTYAABrQ6tWVtXQagRAIgHrSkRCo2MlzZfvBouiV5KYkf7vHgYri4qIo5hdFMX/58uUdrxIAgIaxdWvSp0+1q+gcw4aZDgZA/ehICLQ0yWHbPD40yXPbPB6SZHSSOUVRLEkyNcl9O1scuizLW8qynFSW5aQRI0bse9UAAFADdAIBUE86EgLNS/LGoigOL4qiX5L3J7nvlYNlWbaWZTm8LMtRZVmOSvKzJGeUZTm/SyoGAIAaccABQiAA6sceQ6CyLLckuTzJD5I8nuSusix/UxTF9UVRnNHVBQIA0DjKcs/n1JMBA5KNG6tdBQB0TIdmY5dl+f0k39/huat3ce7M114WAACN6MUXk/33r3YVANAzdWQ6GAAAdIpG2h4eAOqNEAgAgG6zalVlMeVG02jT3ABoTEIgAAC6zapVjdcJNGhQsn59tasAgD0TAgEA0G0acTqYbeIBqBdCIAAAuk0jTgcTAgFQL4RAAAB0G51AAFA9QiAAALpNI3YCHXCAEAiA+iAEAgCg2zTiwtA6gQCoF0IgAAC6zerVSXNztavoXMOGVcItAKh1QiAAALpNW1vSu3e1q+hcOoEAqBdCIAAAeA2EQADUCyEQAAC8Bs3NSUtLtasAgD0TAgEA0C3a2pKiqHYVna9Xr6Qsq10FAOyZEAgAgG6xZk3S1FTtKrqGEAiAeiAEAgCgW6xc2XjbwwNAPRECAQDQLVasSIYPr3YVXaNXr2Tr1mpXAQC7JwQCAKBbrFiRjBhR7Sq6xtChFocGoPYJgQAA6BbLlzduJ5Bt4gGoB0IgAAC6RSN3AgmBAKgHQiAAALpFI68JJAQCoB4IgQAA6BamgwFAdQmBAADoFitXVsKSRiQEAqAeCIEAAOgWW7YkfftWu4quMXx4ZbobANQyIRAAALxGI0ZUprsBQC0TAgEAwGukEwiAeiAEAgCgy23enPTpU+0quk7fvpXpbgBQy4RAAAB0uVWrGndRaACoF0IgAAC63IoVjbs9PADUCyEQAABdbvnyyuLJjaxPn8q0NwCoVUIgAAC6XE/oBLI4NAC1TggEAECX6wkhkG3iAah1QiAAALpcT5gOJgQCoNYJgQAA6HI6gQCg+oRAAAB0OZ1AAFB9QiAAALrcmjXJkCHVrqJrjRiRLFtW7SoAYNeEQAAAdLmiqHw1Mp1AANQ6IRAAAF2uLKtdQdcTAgFQ64RAAAB0qZ4QACXJgAHJxo3VrgIAdk0IBABAl1q/Phk8uNpVAABCIAAAutTy5Y2/PTwA1AMhEAAAXWrZsmTkyGpX0T169Uq2bq12FQCwc0IgAAC61Asv9JwQaNiwZOXKalcBADsnBAIAoEstW5YceGC1q+gedggDoJYJgQAA6FIvvCAEAoBaIAQCAKBLCYEAoDYIgQAA6FI9aU0gIRAAtUwIBABAl1q7Nhk8uNpVdI8RIyprIAFALRICAQDQ5Yqi2hV0j4MOEgIBULuEQAAA0ElGjkz+8z+rXQUA7JwQCACALrN5c9KnT7Wr6D4DBiQbN1a7CgDYOSEQAABdZvnyyjo5PUlPmfoGQP0RAgEA0GV60vbwryjLyhcA1BohEAAAXWbZsv+/vbsPlqyuzwT+fGVmgOV1AsjwNkAQX9CoGAJLUkYirK/Et6iLKXeNa8qXVcvEaKlrSi21KsbE3dJFE11jRVMaX6OCooYoaKWCKCjyKjDK2zgIo8wAo+Aw8Ns/TiOX4V7mMtP3nr59Pp+qLrr7nHvuczlz+nQ//TvnDK8E2mOP7opoADBplEAAACyYG27oTpY8JKtWdX83AEwaJRAAAAtmiIeDrVrlCmEATCYlEAAAC2aIJdD++yuBAJhMSiAAABbMjTc6HAwAJoUSCACABbNhQ7JyZd8pFpfDwQCYVEogAAAWTGvJgwb2jtPhYABMqoHtkgEAWEyt9Z1g8e23X7J+fd8pAOC+lEAAACyIu+4a3iigJFm+PNmype8UAHBfA9wtAwCwGNav70bFAACTQQkEAMCCWLcuOfDAvlP0Z4iHwgEw2ZRAAAAsiOuvTw44oO8U/Vi5Mtm4se8UAHBvSiAAABbEkEcC7b9/csMNfacAgHtTAgEAsCCuv364JdCqVS4TD8DkUQIBALAg1q0b7uFgSiAAJpESCACABfHTn3ZlyBCtWtWNhAKASaIEAgBgQWzenKxY0XeKfhx4YDcSCgAmiRIIAADG7KCDkp/8pO8UAHBvSiAAAMburruSBw34neaeeya33NJ3CgC4twHvmgEAWCjr1yf77dd3iv5U9Z0AAO5LCQQAwNgN+fLwM7XWdwIAuIcSCACAsRvy5eHvtnJlsmFD3ykA4B5KIAAAxs5IICeHBmDyKIEAABg7I4GUQABMHiUQAABjZyRQVwKtW9d3CgC4hxIIAICxu/76ZNWqvlP068ADjQQCYLLMqwSqqqdU1eVVtaaq3jjL9NdW1aVVdWFVfb2qDh1/VAAAloo77khWrOg7Rb8cDgbApNlmCVRVOyV5f5KnJjkqyQuq6qitZvt+kmNaa49O8tkk7x53UAAAWEpWrUp++tO+UwDAPeYzEujYJGtaaz9urW1O8skkz5w5Q2vtrNbaL0cPv53k4PHGBABgqfjVr4wCSpLly5MtW/pOAQD3mE8JdFCS62Y8Xjt6bi4vSfKVHQkFAMDStXZtcrCvBAFg4iybxzw1y3Nt1hmrXpjkmCRPmGP6S5O8NElWr149z4gAACwl112XHHJI3ykmw/LlyebNRkYBMBnmMxJobZKZu/GDk9znYpdVdVKSNyd5RmvtV7MtqLX2odbaMa21Y/bbb7/tyQsAwIRTAt3jgAO6K6UBwCSYTwn03SRHVtXhVbUiySlJTps5Q1UdneSD6QqgG8cfEwCApUIJdA9XCANgkmyzBGqtbUnyqiRfS3JZkk+31i6pqrdX1TNGs/1Nkt2TfKaqLqiq0+ZYHAAAU+6665wT6G6HHJJce23fKQCgM59zAqW1dkaSM7Z67i0z7p805lwAACxR69YlBx7Yd4rJcOihybnn9p0CADrzORwMAADm7c47uxMi05VA11zTdwoA6CiBAABggRx0ULJ2bd8pAKCjBAIAYGx++ctk1137TjE5li9PtmzpOwUAdJRAAACMjSuD3VdV3wkAoKMEAgBgbJRA97XXXsnGjX2nAAAlEAAAY3Tttcnq1X2nmCyrVzs5NACTQQkEAMDYXHVVcvjhfaeYLIce2pVjANA3JRAAAGNz1VXJYYf1nWKyuEw8AJNCCQQAwNhs2JCsXNl3isnicDAAJoUSCACAsXI1rHtzOBgAk0IJBADAWNx+e7Lzzn2nmDy7755s2tR3CgBQAgEAMCbXXNONegEAJpMSCACAsXBlsLntskty2219pwBg6JRAAACMhRJobr/5m8mPf9x3CgCGTgkEAMBYXH21EmguRxyR/OhHfacAYOiUQAAAjMVVVyWHHdZ3isn0kIcka9b0nQKAoVMCAQAwFr/4RXclLO7LSCAAJoESCACAsWit7wST65BDkmuv7TsFAEOnBAIAYIdt3JjsvXffKSbXsmXJnXf2nQKAoVMCAQCww668MjnyyL5TTLZly5ItW/pOAcCQKYEAANhhV1yRPPShfaeYbKtXOyQMgH4pgQAA2GFGAm3bEUe4QhgA/VICAQCww664Qgm0LQ95iCuEAdAvJRAAADvs5puTlSv7TjHZjAQCoG9KIAAAdohLw8/PEUcYCQRAv5RAAADskPXrkwc/uO8Uk2/nnZPNm/tOAcCQKYEAANghzgc0fzvvnNx+e98pABgqJRAAADvkyitdHn6+HvrQ7v8XAPRBCQQAwA754Q+Thz2s7xRLw8Mfnlx2Wd8pABgqJRAAADvk8suNBJqvRzxCCQRAf5RAAADskM2bu3PdsG1KIAD6pAQCAGC7bdqU7LZb3ymWjr32Sm65pe8UAAyVEggAgO126aXJUUf1nWJpedCDkjvv7DsFAEOkBAIAYLtdcknyyEf2nWJpOfTQ5Jpr+k4BwBApgQAA2G5KoAfukY9MLr647xQADJESCACA7XbFFcmRR/adYml57GOTH/yg7xQADJESCACA7XbHHcmKFX2nWFp+67eSCy/sOwUAQ6QEAgBgu2zYkOy9d98plp499uiuqgYAi00JBADAdvn+95Ojj+47xdK0++7Jrbf2nQKAoVECAQCwXb7//eRxj+s7xdL06EcnF13UdwoAhkYJBADAJSYnvgAAD1tJREFUdvne94wE2l6PeUxywQV9pwBgaJRAAABsl5//PNlnn75TLE2uEAZAH5RAAAA8YJs2Jbvt1neKpeuQQ5Jrr+07BQBDowQCAOABu/DCbjQL26cqWbEiuf32vpMAMCRKIAAAHrDzz3dS6B119NHdybUBYLEogQAAeMDOOSc57ri+Uyxtxx2XnHtu3ykAGBIlEAAAD9jPf57su2/fKZa2Y49NvvOdvlMAMCRKIAAAHpB165IDD+w7xdK3zz7Jhg19pwBgSJRAAAA8IOeckxx/fN8ppsO++yY33th3CgCGQgkEAMADcs45ye/+bt8ppsPxxyf/8R99pwBgKJRAAAA8IJdemhx1VN8ppsMTnpCcfXbfKQAYCiUQAADztmlTsssuyYO8ixyLo47qSjUAWAx23wAAzNu3vpU8/vF9p5geVd0Jon/2s76TADAESiAAAObt619PTjyx7xTT5QlPSL75zb5TADAESiAAAObt4ouTRz2q7xTT5YlPTM46q+8UAAyBEggAgHlZv747dMn5gMbryCOTyy/vOwUAQ2AXDgDAvJx5ZnLSSX2nmD5VyRFHJGvW9J0EgGmnBAIAYF6+9KXk6U/vO8V0Ovnk7v8vACwkJRAAANu0eXOycWOy//59J5lOJ56YfOMbfacAYNopgQAA2Kazz05OOKHvFNNr112T5cuTm2/uOwkA00wJBADANn3xi8kzntF3iun29Kcnp5/edwoAppkSCACA+3XHHckVVyQPf3jfSabbc56TfO5zfacAYJopgQAAuF9nnpk86Ul9p5h+e++drFiR3Hhj30kAmFZKIAAA7tcnPpH88R/3nWIYnv/85DOf6TsFANNKCQQAwJw2bEhuvTU56KC+kwzD05+efOELSWt9JwFgGimBAACY00c+kvzJn/SdYjh22SV53OOSc87pOwkA00gJBADArO68M/nyl5M//MO+kwzLK16RfOADfacAYBopgQAAmNXppydPfnKybFnfSYblsMOSzZuT667rOwkA00YJBADAfbSWnHpqNyqFxff61yfvelffKQCYNkogAADu44tfTP7gD5I99+w7yTD9zu8kN9xgNBAA46UEAgDgXjZvTt73vuTVr+47ybC9+c3JW9/adwoApokSCACAe3nve5MXv9gooL4dfXSy667Jt77VdxIApoUSCACAX7vyyuSb30xe+MK+k5Ak73xn8ra3Jbfd1ncSAKaBEggAgCTdYWCvfGXy/vcnVX2nIUlWrkxe97rkL/6i7yQATAMlEAAAaS15zWuSl70sOfTQvtMw09Oeluy+e/Kxj/WdBIClTgkEAEDe+c7k8MOTP/qjvpMwm7/6q+SMM7obAGwvJRAAwIC1lrzpTd2hYK9/fd9pmMtOOyX/+I/Jhz+cfPrTfacBYKlSAgEADNRNNyUveEFy8MHJO97hPECTbpddkk99Kvn615O//Mvkjjv6TgTAUjOvEqiqnlJVl1fVmqp64yzTd66qT42mn1tVh407KAAA43HXXclHP5o897nJa1/bnQyapWH58uSDH0we9ajkqU9Nzj6770QALCXbLIGqaqck70/y1CRHJXlBVR211WwvSbKhtfaQJP8nyV+POygAADtm48bk1FOTk05Kbr45+epXk2OP7TsV2+OUU7pRQaefnpx8cvKZz3SH9AHA/Vk2j3mOTbKmtfbjJKmqTyZ5ZpJLZ8zzzCRvG93/bJJTq6paa22MWQEAeABuvDG54ILk/POTf//37nCi5z+/K39WrOg7HTtqn32S97wn2bAh+fjHk2c9K9ltt+SEE5LHPjZ59KOTPfboOyUAk2Q+JdBBSa6b8XhtkuPmmqe1tqWqbk6yT5KfjSPkJLvyyuTaa8ezrHFXZuNcnmz9L2vcyxtKtqH8neNenmz9L2vcy5vUZY17eUPKtmVL8otfJJs23XO79dZk/fp7ftc++3RlwPHHJ3/+510JxPRZuTJ51au628aNXeF31lnJe9/b/RtJutJv332TPffsiqE99+wKo+XLk2XLutvM+zvtdO9zRM12f1vT7+8+wKTZZ5/kMY/pO8XCm08JNNvL9dZvY+YzT6rqpUlemiSrV6+ex6+efNdfn1x00fiWN+6d4ziXJ1v/yxr38oaSbSh/57iXJ1v/yxr38iZ5WZOabdzLG+eydtqpO6Hzbrslu+/e3fbYo/ug/yCX/hisvffuDg87+eR7P/+rX3UnAr/llq4svOWWriDasqW73XFHcttt9zzesuWen51ZYN59f7bn5nsfYBIdfrgS6G5rkxwy4/HBSdbNMc/aqlqWZK8kN229oNbah5J8KEmOOeaYqdgV/P7vdzcAAJhUO++cHHBAdwNguObzPdF3kxxZVYdX1YokpyQ5bat5TkvyotH95yb5hvMBAQAAAEyObY4EGp3j51VJvpZkpyQfaa1dUlVvT3Jea+20JP+Q5J+qak26EUCnLGRoAAAAAB6Y+RwOltbaGUnO2Oq5t8y4f3uS5403GgAAAADj4rSBAAAAAAOgBAIAAAAYACUQAAAAwAAogQAAAAAGQAkEAAAAMABKIAAAAIABUAIBAAAADIASCAAAAGAAlEAAAAAAA6AEAgAAABgAJRAAAADAACiBAAAAAAZACQQAAAAwAEogAAAAgAFQAgEAAAAMgBIIAAAAYACUQAAAAAADoAQCAAAAGAAlEAAAAMAAKIEAAAAABqBaa/384qr1Sa7p5ZeP375JftZ3CHph3Q+T9T5c1v1wWffDZd0Pl3U/XNb9ME3Tej+0tbbfbBN6K4GmSVWd11o7pu8cLD7rfpis9+Gy7ofLuh8u6364rPvhsu6HaSjr3eFgAAAAAAOgBAIAAAAYACXQeHyo7wD0xrofJut9uKz74bLuh8u6Hy7rfris+2EaxHp3TiAAAACAATASCAAAAGAAlEDzVFXPq6pLququqjpmq2lvqqo1VXV5VT15jp8/vKrOraorq+pTVbVicZIzTqN1d8HodnVVXTDHfFdX1UWj+c5b7JyMV1W9rap+MmPdP22O+Z4yeh1YU1VvXOycjF9V/U1V/bCqLqyqz1fV3nPMZ5ufEtvajqtq59G+YM1ov37Y4qdk3KrqkKo6q6ouG73fe80s85xQVTfP2Be8pY+sjN+2XsOr877Rdn9hVT2uj5yMT1U9bMa2fEFV3VJVf7bVPLb5KVFVH6mqG6vq4hnP/UZVnTn6fH5mVa2c42dfNJrnyqp60eKlXjgOB5unqnpEkruSfDDJ61pr542ePyrJPyc5NsmBSf4tyUNba3du9fOfTvIvrbVPVtXfJ/lBa+3vFvNvYLyq6j1Jbm6tvX2WaVcnOaa19rNFD8bYVdXbkmxqrf3t/cyzU5IrkvyXJGuTfDfJC1prly5KSBZEVT0pyTdaa1uq6q+TpLX2hlnmuzq2+SVvPttxVf3PJI9urb28qk5J8uzW2n/tJTBjU1UHJDmgtfa9qtojyflJnrXVuj8h3XvAk3uKyQLZ1mv46MufVyd5WpLjkry3tXbc4iVkIY1e+3+S5LjW2jUznj8htvmpUFW/n2RTko+11h41eu7dSW5qrb1r9KXPyq3f41XVbyQ5L8kxSVq6fcNvt9Y2LOofMGZGAs1Ta+2y1trls0x6ZpJPttZ+1Vq7KsmadIXQr1VVJXliks+OnvpokmctZF4W1midPj9dAQhJt92vaa39uLW2Ockn070+sIS11v61tbZl9PDbSQ7uMw8Lbj7b8TPT7ceTbr9+4mifwBLWWru+tfa90f1bk1yW5KB+UzFBnpnuw2NrrX07yd6j4pDpcGKSH80sgJgurbVvJblpq6dn7s/n+nz+5CRnttZuGhU/ZyZ5yoIFXSRKoB13UJLrZjxem/u+adgnycYZHyRmm4el5fFJbmitXTnH9JbkX6vq/Kp66SLmYuG8ajQE/CNzDBedz2sBS9v/SPKVOabZ5qfDfLbjX88z2q/fnG4/z5QYHeJ3dJJzZ5l8fFX9oKq+UlWPXNRgLKRtvYbbx0+3UzL3F7u2+em1f2vt+qT7IiDJg2eZZyq3/WV9B5gkVfVvSVbNMunNrbUvzvVjszy39TF285mHCTHPfwcvyP2PAvq91tq6qnpwkjOr6oejBpoJdX/rPcnfJXlHuu32HUnek64QuNciZvlZ2/kSMJ9tvqrenGRLko/PsRjb/HSwTx+4qto9yeeS/Flr7ZatJn8vyaGttU2jw4O+kOTIxc7IgtjWa7jtfkpVd57WZyR50yyTbfNM5bavBJqhtXbSdvzY2iSHzHh8cJJ1W83zs3TDRpeNvjWcbR4mxLb+HVTVsiTPSfLb97OMdaP/3lhVn093iIEPhBNsvtt/Vf2/JF+aZdJ8XguYQPPY5l+U5OQkJ7Y5TqRnm58a89mO755n7Wh/sFfuO8ScJaiqlqcrgD7eWvuXrafPLIVaa2dU1Qeqal/nAlv65vEabh8/vZ6a5HuttRu2nmCbn3o3VNUBrbXrR4d33jjLPGuTnDDj8cFJzl6EbAvK4WA77rQkp4yuFnJ4unb4OzNnGH1oOCvJc0dPvSjJXCOLmHwnJflha23tbBOrarfRSSVTVbsleVKSi2ebl6Vhq+P+n53Z1+d3kxxZ3ZUAV6QbWnzaYuRj4VTVU5K8IckzWmu/nGMe2/z0mM92fFq6/XjS7de/MVc5yNIxOq/TPyS5rLX2v+eYZ9Xd53+qqmPTvY/++eKlZCHM8zX8tCT/vTr/Od2FQa5f5KgsjDlH99vmp97M/flcn8+/luRJVbVydDqIJ42eW9KMBJqnqnp2kv+bZL8kX66qC1prT26tXTK68tel6Q4VeOXdVwarqjOS/Ono24U3JPlkVb0zyffTvdFgabrPccNVdWCSD7fWnpZk/ySfH+0zliX5RGvtq4ueknF6d1U9Nt3wz6uTvCy593ofXT3qVel2DDsl+Uhr7ZK+AjM2pybZOd3hAUny7dFVoWzzU2iu7biq3p7kvNbaaen23/9UVWvSjQA6pb/EjNHvJflvSS6qqgtGz/2vJKuTpLX29+lKv1dU1ZYktyU5RQE4FWZ9Da+qlye/XvdnpLsy2Jokv0zy4p6yMkZV9Z/SXQ3yZTOem7nebfNToqr+Od2Inn2ram2StyZ5V5JPV9VLklyb5HmjeY9J8vLW2p+21m6qqnek+5IoSd7eWlvyo39dIh4AAABgABwOBgAAADAASiAAAACAAVACAQAAAAyAEggAAABgAJRAAAAAAAOgBAIAAAAYACUQAAAAwAAogQAAAAAG4P8DGPGJXdJ0GiIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mu=np.linspace(-10,10,10000)\n",
    "y_mu=np.zeros(10000)\n",
    "for i in range(10000):\n",
    "    y_mu[i] = unnormalized(x, mu[i], prior_var, clutter_var, var)\n",
    "#y_mu=unnormalized(x, mu, prior_var, clutter_var)\n",
    "normalizer = simps(y_mu, mu)\n",
    "y_mu_normalized=y_mu/normalizer\n",
    "ys=normal_pdf(mu, m_posterior, v_posterior)\n",
    "\n",
    "plt.figure(figsize=(20,10))\n",
    "plt.plot(mu, y_mu_normalized, color=\"green\", linewidth=0.7, label=\"approx with num integration\")\n",
    "plt.axvline(x=2, color=\"red\", label=\"true generating value\")\n",
    "plt.plot(mu,ys, color=\"blue\", linewidth=0.7, label=\"current posterior estimate\")\n",
    "plt.legend(loc=\"upper left\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}