src/rl_tb_lidar/src/utils/sensormodel/lidar_sensor_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LIDAR Sensor Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import necessary libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Data\n",
    "### Samples from real turtlebot LIDAR\n",
    "Load data from real turtlebot into a pandas DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Program_Files\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:10: YAMLLoadWarning: calling yaml.load_all() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "  # Remove the CWD from sys.path while we load stuff.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(565, 361)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lidar0</th>\n",
       "      <th>lidar1</th>\n",
       "      <th>lidar2</th>\n",
       "      <th>lidar3</th>\n",
       "      <th>lidar4</th>\n",
       "      <th>lidar5</th>\n",
       "      <th>lidar6</th>\n",
       "      <th>lidar7</th>\n",
       "      <th>lidar8</th>\n",
       "      <th>lidar9</th>\n",
       "      <th>...</th>\n",
       "      <th>lidar351</th>\n",
       "      <th>lidar352</th>\n",
       "      <th>lidar353</th>\n",
       "      <th>lidar354</th>\n",
       "      <th>lidar355</th>\n",
       "      <th>lidar356</th>\n",
       "      <th>lidar357</th>\n",
       "      <th>lidar358</th>\n",
       "      <th>lidar359</th>\n",
       "      <th>pos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.7490</td>\n",
       "      <td>2.66950</td>\n",
       "      <td>3.74000</td>\n",
       "      <td>3.31125</td>\n",
       "      <td>3.02250</td>\n",
       "      <td>2.73900</td>\n",
       "      <td>2.49850</td>\n",
       "      <td>2.50225</td>\n",
       "      <td>2.50850</td>\n",
       "      <td>2.50575</td>\n",
       "      <td>...</td>\n",
       "      <td>3.74150</td>\n",
       "      <td>3.73525</td>\n",
       "      <td>3.73525</td>\n",
       "      <td>3.74550</td>\n",
       "      <td>3.74000</td>\n",
       "      <td>3.73875</td>\n",
       "      <td>3.73650</td>\n",
       "      <td>3.74550</td>\n",
       "      <td>3.65625</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.7470</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.73800</td>\n",
       "      <td>3.34500</td>\n",
       "      <td>2.99675</td>\n",
       "      <td>2.76375</td>\n",
       "      <td>2.49900</td>\n",
       "      <td>2.50175</td>\n",
       "      <td>2.50475</td>\n",
       "      <td>2.50375</td>\n",
       "      <td>...</td>\n",
       "      <td>3.74075</td>\n",
       "      <td>3.74075</td>\n",
       "      <td>3.73450</td>\n",
       "      <td>3.74275</td>\n",
       "      <td>3.74000</td>\n",
       "      <td>3.73725</td>\n",
       "      <td>3.73650</td>\n",
       "      <td>3.77975</td>\n",
       "      <td>3.69900</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.7505</td>\n",
       "      <td>3.70375</td>\n",
       "      <td>3.69350</td>\n",
       "      <td>3.30850</td>\n",
       "      <td>2.97150</td>\n",
       "      <td>2.71600</td>\n",
       "      <td>2.50600</td>\n",
       "      <td>2.50025</td>\n",
       "      <td>2.50725</td>\n",
       "      <td>2.50850</td>\n",
       "      <td>...</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.73725</td>\n",
       "      <td>3.74350</td>\n",
       "      <td>3.73800</td>\n",
       "      <td>3.73375</td>\n",
       "      <td>3.73875</td>\n",
       "      <td>3.78400</td>\n",
       "      <td>3.69775</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.7470</td>\n",
       "      <td>2.67125</td>\n",
       "      <td>3.78825</td>\n",
       "      <td>3.57000</td>\n",
       "      <td>3.11550</td>\n",
       "      <td>2.60650</td>\n",
       "      <td>2.50050</td>\n",
       "      <td>2.50200</td>\n",
       "      <td>2.50425</td>\n",
       "      <td>2.50625</td>\n",
       "      <td>...</td>\n",
       "      <td>3.73025</td>\n",
       "      <td>3.73300</td>\n",
       "      <td>3.74150</td>\n",
       "      <td>3.77775</td>\n",
       "      <td>3.74475</td>\n",
       "      <td>3.73650</td>\n",
       "      <td>3.74150</td>\n",
       "      <td>3.78050</td>\n",
       "      <td>3.65550</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.7870</td>\n",
       "      <td>16.00000</td>\n",
       "      <td>3.83575</td>\n",
       "      <td>3.27550</td>\n",
       "      <td>2.96825</td>\n",
       "      <td>2.71450</td>\n",
       "      <td>2.50175</td>\n",
       "      <td>2.50000</td>\n",
       "      <td>2.50725</td>\n",
       "      <td>2.51300</td>\n",
       "      <td>...</td>\n",
       "      <td>3.77700</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.73525</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.74275</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.74475</td>\n",
       "      <td>3.75325</td>\n",
       "      <td>3.74550</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 361 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   lidar0    lidar1   lidar2   lidar3   lidar4   lidar5   lidar6   lidar7  \\\n",
       "0  3.7490   2.66950  3.74000  3.31125  3.02250  2.73900  2.49850  2.50225   \n",
       "1  3.7470   3.73925  3.73800  3.34500  2.99675  2.76375  2.49900  2.50175   \n",
       "2  3.7505   3.70375  3.69350  3.30850  2.97150  2.71600  2.50600  2.50025   \n",
       "3  3.7470   2.67125  3.78825  3.57000  3.11550  2.60650  2.50050  2.50200   \n",
       "4  3.7870  16.00000  3.83575  3.27550  2.96825  2.71450  2.50175  2.50000   \n",
       "\n",
       "    lidar8   lidar9  ...  lidar351  lidar352  lidar353  lidar354  lidar355  \\\n",
       "0  2.50850  2.50575  ...   3.74150   3.73525   3.73525   3.74550   3.74000   \n",
       "1  2.50475  2.50375  ...   3.74075   3.74075   3.73450   3.74275   3.74000   \n",
       "2  2.50725  2.50850  ...   3.73925   3.73925   3.73725   3.74350   3.73800   \n",
       "3  2.50425  2.50625  ...   3.73025   3.73300   3.74150   3.77775   3.74475   \n",
       "4  2.50725  2.51300  ...   3.77700   3.73925   3.73525   3.73925   3.74275   \n",
       "\n",
       "   lidar356  lidar357  lidar358  lidar359  pos  \n",
       "0   3.73875   3.73650   3.74550   3.65625    1  \n",
       "1   3.73725   3.73650   3.77975   3.69900    1  \n",
       "2   3.73375   3.73875   3.78400   3.69775    1  \n",
       "3   3.73650   3.74150   3.78050   3.65550    1  \n",
       "4   3.73925   3.74475   3.75325   3.74550    1  \n",
       "\n",
       "[5 rows x 361 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MAX_RANGE = 16.0\n",
    "TOTAL_POSITIONS = 6\n",
    "data_files = [\"mea_\"+str(i)+\".txt\" for i in np.arange(1, TOTAL_POSITIONS+1)]\n",
    "cols = ['lidar'+str(i) for i in range(360)]\n",
    "data = pd.DataFrame(columns = cols+[\"pos\"])\n",
    "for i, file in enumerate(data_files, 1):\n",
    "    stream = open(file, \"r\")\n",
    "    docs = yaml.load_all(stream)\n",
    "    file_data = pd.DataFrame(columns = cols)\n",
    "    for doc in docs:\n",
    "        if doc is not None:\n",
    "            for k, v in doc.items():\n",
    "                if k == 'ranges':\n",
    "                    data_line = pd.DataFrame(columns = cols, data = [v], dtype=np.float)\n",
    "                    file_data = file_data.append(data_line, ignore_index = True)\n",
    "    file_data[\"pos\"] = i\n",
    "    data = data.append(file_data, ignore_index = True)\n",
    "data = data.replace(np.inf, MAX_RANGE)\n",
    "data[data>MAX_RANGE] = MAX_RANGE\n",
    "print(data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data from Stage simulator\n",
    "Load simulation data from Stage Turtlebot into a pandas DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Program_Files\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: YAMLLoadWarning: calling yaml.load_all() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(60, 361)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lidar0</th>\n",
       "      <th>lidar1</th>\n",
       "      <th>lidar2</th>\n",
       "      <th>lidar3</th>\n",
       "      <th>lidar4</th>\n",
       "      <th>lidar5</th>\n",
       "      <th>lidar6</th>\n",
       "      <th>lidar7</th>\n",
       "      <th>lidar8</th>\n",
       "      <th>lidar9</th>\n",
       "      <th>...</th>\n",
       "      <th>lidar351</th>\n",
       "      <th>lidar352</th>\n",
       "      <th>lidar353</th>\n",
       "      <th>lidar354</th>\n",
       "      <th>lidar355</th>\n",
       "      <th>lidar356</th>\n",
       "      <th>lidar357</th>\n",
       "      <th>lidar358</th>\n",
       "      <th>lidar359</th>\n",
       "      <th>pos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.67619</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.57681</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.67619</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.57681</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.67619</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.57681</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.67619</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.57681</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.67619</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.57681</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 361 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   lidar0    lidar1    lidar2    lidar3    lidar4    lidar5    lidar6  \\\n",
       "0    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "1    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "2    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "3    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "4    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "\n",
       "     lidar7   lidar8    lidar9  ...  lidar351  lidar352  lidar353  lidar354  \\\n",
       "0  2.851372  2.67619  2.683219  ...  3.585084   3.57681  3.569664  3.563636   \n",
       "1  2.851372  2.67619  2.683219  ...  3.585084   3.57681  3.569664  3.563636   \n",
       "2  2.851372  2.67619  2.683219  ...  3.585084   3.57681  3.569664  3.563636   \n",
       "3  2.851372  2.67619  2.683219  ...  3.585084   3.57681  3.569664  3.563636   \n",
       "4  2.851372  2.67619  2.683219  ...  3.585084   3.57681  3.569664  3.563636   \n",
       "\n",
       "   lidar355  lidar356  lidar357  lidar358  lidar359  pos  \n",
       "0  3.558717  3.554899  3.552176  3.550544      3.55    1  \n",
       "1  3.558717  3.554899  3.552176  3.550544      3.55    1  \n",
       "2  3.558717  3.554899  3.552176  3.550544      3.55    1  \n",
       "3  3.558717  3.554899  3.552176  3.550544      3.55    1  \n",
       "4  3.558717  3.554899  3.552176  3.550544      3.55    1  \n",
       "\n",
       "[5 rows x 361 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_files = [\"result_pos\"+str(i)+\"_f.txt\" for i in np.arange(1, TOTAL_POSITIONS+1)]\n",
    "cols = ['lidar'+str(i) for i in range(360)]\n",
    "gt_ = pd.DataFrame(columns = cols+[\"pos\"])\n",
    "for i, file in enumerate(data_files, 1):\n",
    "    stream = open(file, \"r\")\n",
    "    docs = yaml.load_all(stream)\n",
    "    file_data = pd.DataFrame(columns = cols)\n",
    "    for doc in docs:\n",
    "        if doc is not None:\n",
    "            for k, v in doc.items():\n",
    "                if k == 'ranges':\n",
    "                    data_line = pd.DataFrame(columns = cols, data = [v], dtype=np.float)\n",
    "                    file_data = file_data.append(data_line, ignore_index = True)\n",
    "    file_data[\"pos\"] = i\n",
    "    gt_ = gt_.append(file_data, ignore_index = True)\n",
    "gt_ = gt_.replace(np.inf, MAX_RANGE)\n",
    "gt_[gt_>MAX_RANGE] = MAX_RANGE\n",
    "print(gt_.shape)\n",
    "gt_.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ground truth datasets\n",
    "Average the simulation data for each position to use as ground truth dataset 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "E:\\Program_Files\\Anaconda3\\lib\\site-packages\\pandas\\core\\frame.py:6692: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version\n",
      "of pandas will change to not sort by default.\n",
      "\n",
      "To accept the future behavior, pass 'sort=False'.\n",
      "\n",
      "To retain the current behavior and silence the warning, pass 'sort=True'.\n",
      "\n",
      "  sort=sort)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lidar0</th>\n",
       "      <th>lidar1</th>\n",
       "      <th>lidar2</th>\n",
       "      <th>lidar3</th>\n",
       "      <th>lidar4</th>\n",
       "      <th>lidar5</th>\n",
       "      <th>lidar6</th>\n",
       "      <th>lidar7</th>\n",
       "      <th>lidar8</th>\n",
       "      <th>lidar9</th>\n",
       "      <th>...</th>\n",
       "      <th>lidar351</th>\n",
       "      <th>lidar352</th>\n",
       "      <th>lidar353</th>\n",
       "      <th>lidar354</th>\n",
       "      <th>lidar355</th>\n",
       "      <th>lidar356</th>\n",
       "      <th>lidar357</th>\n",
       "      <th>lidar358</th>\n",
       "      <th>lidar359</th>\n",
       "      <th>pos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.55</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.545198</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.187559</td>\n",
       "      <td>2.851372</td>\n",
       "      <td>2.676190</td>\n",
       "      <td>2.683219</td>\n",
       "      <td>...</td>\n",
       "      <td>3.585084</td>\n",
       "      <td>3.576810</td>\n",
       "      <td>3.569664</td>\n",
       "      <td>3.563636</td>\n",
       "      <td>3.558717</td>\n",
       "      <td>3.554899</td>\n",
       "      <td>3.552176</td>\n",
       "      <td>3.550544</td>\n",
       "      <td>3.55</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.05</td>\n",
       "      <td>4.050621</td>\n",
       "      <td>4.052483</td>\n",
       "      <td>4.055589</td>\n",
       "      <td>4.059945</td>\n",
       "      <td>3.623867</td>\n",
       "      <td>3.629996</td>\n",
       "      <td>3.637263</td>\n",
       "      <td>3.645677</td>\n",
       "      <td>3.655253</td>\n",
       "      <td>...</td>\n",
       "      <td>4.090025</td>\n",
       "      <td>4.080585</td>\n",
       "      <td>4.072433</td>\n",
       "      <td>4.065557</td>\n",
       "      <td>4.059945</td>\n",
       "      <td>4.055589</td>\n",
       "      <td>4.052483</td>\n",
       "      <td>4.050621</td>\n",
       "      <td>4.05</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.04</td>\n",
       "      <td>1.040159</td>\n",
       "      <td>1.040637</td>\n",
       "      <td>1.041435</td>\n",
       "      <td>1.042554</td>\n",
       "      <td>1.043995</td>\n",
       "      <td>1.045761</td>\n",
       "      <td>1.047854</td>\n",
       "      <td>1.050278</td>\n",
       "      <td>1.053037</td>\n",
       "      <td>...</td>\n",
       "      <td>1.050278</td>\n",
       "      <td>1.047854</td>\n",
       "      <td>1.045761</td>\n",
       "      <td>1.043995</td>\n",
       "      <td>1.042554</td>\n",
       "      <td>1.041435</td>\n",
       "      <td>1.040637</td>\n",
       "      <td>1.040159</td>\n",
       "      <td>1.04</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.95</td>\n",
       "      <td>1.950299</td>\n",
       "      <td>1.951195</td>\n",
       "      <td>1.952691</td>\n",
       "      <td>1.954788</td>\n",
       "      <td>1.957490</td>\n",
       "      <td>1.960801</td>\n",
       "      <td>1.964726</td>\n",
       "      <td>1.969272</td>\n",
       "      <td>1.974444</td>\n",
       "      <td>...</td>\n",
       "      <td>1.969272</td>\n",
       "      <td>1.964726</td>\n",
       "      <td>1.960801</td>\n",
       "      <td>1.957490</td>\n",
       "      <td>1.954788</td>\n",
       "      <td>1.952691</td>\n",
       "      <td>1.951195</td>\n",
       "      <td>1.950299</td>\n",
       "      <td>1.95</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.95</td>\n",
       "      <td>4.950758</td>\n",
       "      <td>4.953034</td>\n",
       "      <td>4.956831</td>\n",
       "      <td>4.962155</td>\n",
       "      <td>4.969014</td>\n",
       "      <td>4.977419</td>\n",
       "      <td>4.987382</td>\n",
       "      <td>4.998920</td>\n",
       "      <td>5.012050</td>\n",
       "      <td>...</td>\n",
       "      <td>4.998920</td>\n",
       "      <td>4.987382</td>\n",
       "      <td>4.977419</td>\n",
       "      <td>4.969014</td>\n",
       "      <td>4.962155</td>\n",
       "      <td>4.956831</td>\n",
       "      <td>4.953034</td>\n",
       "      <td>4.950758</td>\n",
       "      <td>4.95</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.10</td>\n",
       "      <td>0.100015</td>\n",
       "      <td>0.100061</td>\n",
       "      <td>0.100138</td>\n",
       "      <td>0.100246</td>\n",
       "      <td>0.100384</td>\n",
       "      <td>0.100554</td>\n",
       "      <td>0.100755</td>\n",
       "      <td>0.100988</td>\n",
       "      <td>0.101254</td>\n",
       "      <td>...</td>\n",
       "      <td>0.100988</td>\n",
       "      <td>0.100755</td>\n",
       "      <td>0.100554</td>\n",
       "      <td>0.100384</td>\n",
       "      <td>0.100246</td>\n",
       "      <td>0.100138</td>\n",
       "      <td>0.100061</td>\n",
       "      <td>0.100015</td>\n",
       "      <td>0.10</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>6 rows × 361 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "   lidar0    lidar1    lidar2    lidar3    lidar4    lidar5    lidar6  \\\n",
       "0    3.55  3.550544  3.552176  3.554899  3.545198  3.563636  3.187559   \n",
       "1    4.05  4.050621  4.052483  4.055589  4.059945  3.623867  3.629996   \n",
       "2    1.04  1.040159  1.040637  1.041435  1.042554  1.043995  1.045761   \n",
       "3    1.95  1.950299  1.951195  1.952691  1.954788  1.957490  1.960801   \n",
       "4    4.95  4.950758  4.953034  4.956831  4.962155  4.969014  4.977419   \n",
       "5    0.10  0.100015  0.100061  0.100138  0.100246  0.100384  0.100554   \n",
       "\n",
       "     lidar7    lidar8    lidar9  ...  lidar351  lidar352  lidar353  lidar354  \\\n",
       "0  2.851372  2.676190  2.683219  ...  3.585084  3.576810  3.569664  3.563636   \n",
       "1  3.637263  3.645677  3.655253  ...  4.090025  4.080585  4.072433  4.065557   \n",
       "2  1.047854  1.050278  1.053037  ...  1.050278  1.047854  1.045761  1.043995   \n",
       "3  1.964726  1.969272  1.974444  ...  1.969272  1.964726  1.960801  1.957490   \n",
       "4  4.987382  4.998920  5.012050  ...  4.998920  4.987382  4.977419  4.969014   \n",
       "5  0.100755  0.100988  0.101254  ...  0.100988  0.100755  0.100554  0.100384   \n",
       "\n",
       "   lidar355  lidar356  lidar357  lidar358  lidar359  pos  \n",
       "0  3.558717  3.554899  3.552176  3.550544      3.55  1.0  \n",
       "1  4.059945  4.055589  4.052483  4.050621      4.05  2.0  \n",
       "2  1.042554  1.041435  1.040637  1.040159      1.04  3.0  \n",
       "3  1.954788  1.952691  1.951195  1.950299      1.95  4.0  \n",
       "4  4.962155  4.956831  4.953034  4.950758      4.95  5.0  \n",
       "5  0.100246  0.100138  0.100061  0.100015      0.10  6.0  \n",
       "\n",
       "[6 rows x 361 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gt1 = pd.DataFrame(columns = cols)\n",
    "for pos in gt_[\"pos\"].unique():\n",
    "    data_ = gt_[gt_['pos']==pos].replace(np.inf, MAX_RANGE)\n",
    "    data_[data_ > MAX_RANGE] = MAX_RANGE\n",
    "    gt1 = gt1.append(pd.DataFrame(data_.mean(axis=0)).T, ignore_index = True)\n",
    "gt1 = gt1[cols+[\"pos\"]]\n",
    "gt1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take mode for each angle of real turtlebot data for each position to use as ground truth dataset 2. Fill in any remaining infinity (out of range) values with the corresponding values from ground truth dataset 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>lidar0</th>\n",
       "      <th>lidar1</th>\n",
       "      <th>lidar2</th>\n",
       "      <th>lidar3</th>\n",
       "      <th>lidar4</th>\n",
       "      <th>lidar5</th>\n",
       "      <th>lidar6</th>\n",
       "      <th>lidar7</th>\n",
       "      <th>lidar8</th>\n",
       "      <th>lidar9</th>\n",
       "      <th>...</th>\n",
       "      <th>lidar351</th>\n",
       "      <th>lidar352</th>\n",
       "      <th>lidar353</th>\n",
       "      <th>lidar354</th>\n",
       "      <th>lidar355</th>\n",
       "      <th>lidar356</th>\n",
       "      <th>lidar357</th>\n",
       "      <th>lidar358</th>\n",
       "      <th>lidar359</th>\n",
       "      <th>pos</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3.749000</td>\n",
       "      <td>2.669250</td>\n",
       "      <td>3.769417</td>\n",
       "      <td>3.332000</td>\n",
       "      <td>2.969125</td>\n",
       "      <td>2.625750</td>\n",
       "      <td>2.499000</td>\n",
       "      <td>2.502000</td>\n",
       "      <td>2.504500</td>\n",
       "      <td>2.505750</td>\n",
       "      <td>...</td>\n",
       "      <td>3.739250</td>\n",
       "      <td>3.738250</td>\n",
       "      <td>3.73725</td>\n",
       "      <td>3.741500</td>\n",
       "      <td>3.73925</td>\n",
       "      <td>3.740750</td>\n",
       "      <td>3.74625</td>\n",
       "      <td>3.749000</td>\n",
       "      <td>3.660875</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.080000</td>\n",
       "      <td>4.081625</td>\n",
       "      <td>4.084500</td>\n",
       "      <td>4.080500</td>\n",
       "      <td>3.604125</td>\n",
       "      <td>3.585750</td>\n",
       "      <td>3.626750</td>\n",
       "      <td>3.626875</td>\n",
       "      <td>3.627500</td>\n",
       "      <td>3.660000</td>\n",
       "      <td>...</td>\n",
       "      <td>4.087000</td>\n",
       "      <td>4.086250</td>\n",
       "      <td>4.08575</td>\n",
       "      <td>4.084500</td>\n",
       "      <td>4.07950</td>\n",
       "      <td>4.078000</td>\n",
       "      <td>4.08050</td>\n",
       "      <td>4.079500</td>\n",
       "      <td>4.078375</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.966125</td>\n",
       "      <td>0.967500</td>\n",
       "      <td>0.967750</td>\n",
       "      <td>0.968375</td>\n",
       "      <td>0.971188</td>\n",
       "      <td>0.972250</td>\n",
       "      <td>0.976000</td>\n",
       "      <td>0.977750</td>\n",
       "      <td>0.982250</td>\n",
       "      <td>0.984875</td>\n",
       "      <td>...</td>\n",
       "      <td>0.969875</td>\n",
       "      <td>0.967750</td>\n",
       "      <td>0.96650</td>\n",
       "      <td>0.966875</td>\n",
       "      <td>0.96650</td>\n",
       "      <td>0.966000</td>\n",
       "      <td>0.96525</td>\n",
       "      <td>0.964500</td>\n",
       "      <td>0.965625</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.730000</td>\n",
       "      <td>1.728250</td>\n",
       "      <td>1.729750</td>\n",
       "      <td>1.730500</td>\n",
       "      <td>1.768250</td>\n",
       "      <td>1.957490</td>\n",
       "      <td>1.960801</td>\n",
       "      <td>1.370500</td>\n",
       "      <td>1.969272</td>\n",
       "      <td>1.370500</td>\n",
       "      <td>...</td>\n",
       "      <td>1.701500</td>\n",
       "      <td>1.677125</td>\n",
       "      <td>1.70150</td>\n",
       "      <td>0.991375</td>\n",
       "      <td>0.99325</td>\n",
       "      <td>0.998333</td>\n",
       "      <td>2.04375</td>\n",
       "      <td>1.778250</td>\n",
       "      <td>1.748750</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4.183250</td>\n",
       "      <td>5.446000</td>\n",
       "      <td>4.953034</td>\n",
       "      <td>4.956831</td>\n",
       "      <td>4.962155</td>\n",
       "      <td>4.969014</td>\n",
       "      <td>4.977419</td>\n",
       "      <td>5.447250</td>\n",
       "      <td>5.447250</td>\n",
       "      <td>5.488625</td>\n",
       "      <td>...</td>\n",
       "      <td>5.444500</td>\n",
       "      <td>5.443000</td>\n",
       "      <td>5.44450</td>\n",
       "      <td>5.180000</td>\n",
       "      <td>5.18125</td>\n",
       "      <td>5.178750</td>\n",
       "      <td>5.17475</td>\n",
       "      <td>4.950758</td>\n",
       "      <td>4.850500</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.180500</td>\n",
       "      <td>0.181000</td>\n",
       "      <td>0.181250</td>\n",
       "      <td>0.181500</td>\n",
       "      <td>0.182250</td>\n",
       "      <td>0.182500</td>\n",
       "      <td>0.182750</td>\n",
       "      <td>0.183500</td>\n",
       "      <td>0.184250</td>\n",
       "      <td>0.185000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.181250</td>\n",
       "      <td>0.181000</td>\n",
       "      <td>0.18050</td>\n",
       "      <td>0.180000</td>\n",
       "      <td>0.18025</td>\n",
       "      <td>0.180250</td>\n",
       "      <td>0.18050</td>\n",
       "      <td>0.180250</td>\n",
       "      <td>0.180500</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>6 rows × 361 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     lidar0    lidar1    lidar2    lidar3    lidar4    lidar5    lidar6  \\\n",
       "0  3.749000  2.669250  3.769417  3.332000  2.969125  2.625750  2.499000   \n",
       "1  4.080000  4.081625  4.084500  4.080500  3.604125  3.585750  3.626750   \n",
       "2  0.966125  0.967500  0.967750  0.968375  0.971188  0.972250  0.976000   \n",
       "3  1.730000  1.728250  1.729750  1.730500  1.768250  1.957490  1.960801   \n",
       "4  4.183250  5.446000  4.953034  4.956831  4.962155  4.969014  4.977419   \n",
       "5  0.180500  0.181000  0.181250  0.181500  0.182250  0.182500  0.182750   \n",
       "\n",
       "     lidar7    lidar8    lidar9  ...  lidar351  lidar352  lidar353  lidar354  \\\n",
       "0  2.502000  2.504500  2.505750  ...  3.739250  3.738250   3.73725  3.741500   \n",
       "1  3.626875  3.627500  3.660000  ...  4.087000  4.086250   4.08575  4.084500   \n",
       "2  0.977750  0.982250  0.984875  ...  0.969875  0.967750   0.96650  0.966875   \n",
       "3  1.370500  1.969272  1.370500  ...  1.701500  1.677125   1.70150  0.991375   \n",
       "4  5.447250  5.447250  5.488625  ...  5.444500  5.443000   5.44450  5.180000   \n",
       "5  0.183500  0.184250  0.185000  ...  0.181250  0.181000   0.18050  0.180000   \n",
       "\n",
       "   lidar355  lidar356  lidar357  lidar358  lidar359  pos  \n",
       "0   3.73925  3.740750   3.74625  3.749000  3.660875  1.0  \n",
       "1   4.07950  4.078000   4.08050  4.079500  4.078375  2.0  \n",
       "2   0.96650  0.966000   0.96525  0.964500  0.965625  3.0  \n",
       "3   0.99325  0.998333   2.04375  1.778250  1.748750  4.0  \n",
       "4   5.18125  5.178750   5.17475  4.950758  4.850500  5.0  \n",
       "5   0.18025  0.180250   0.18050  0.180250  0.180500  6.0  \n",
       "\n",
       "[6 rows x 361 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gt2 = pd.DataFrame(columns = cols)\n",
    "for pos in gt1[\"pos\"].unique():\n",
    "    data_ = data[data['pos']==pos].replace(np.inf, np.NaN).replace(16.0, np.NaN).mode().mean()\n",
    "    data_[data_.isnull()] = gt1[gt1['pos']==pos].mean()\n",
    "#     data_[data_ == np.NaN] = pd.DataFrame(gt1[gt1['pos'==pos]].mean()).T[data_ == np.NaN]\n",
    "    gt2 = gt2.append(pd.DataFrame(data_).T, ignore_index = True)\n",
    "gt2 = gt2[cols+[\"pos\"]]\n",
    "gt2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit Sensor Model\n",
    "Start with a reasonable initialization, like a uniform distribution with unit parameters for the truncated Gaussian and exponential components of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.asarray([0.25, 0.25, 0.25, 0.25, 1.0, 1.0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use expectation maximization to do maximum likelihood estimation of the true parameters. We first do this for ground truth dataset 1 and then for ground truth dataset 2. We shall test both resulting models on the real turtlebot to see which one generalizes better.\n",
    "### Ground truth 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration:  0\n",
      "likelihood:  7.909989209134689\n",
      "params:  [0.52024238 0.05069461 0.34746964 0.08159337 0.34435953 0.64557225]\n",
      "iteration:  1\n",
      "likelihood:  8.005912180006494\n",
      "params:  [0.55412964 0.03156688 0.34746964 0.06683384 0.25478053 0.68047807]\n",
      "iteration:  2\n",
      "likelihood:  8.026942039582943\n",
      "params:  [0.54316022 0.03305634 0.34746964 0.07631381 0.22652867 0.63887563]\n",
      "iteration:  3\n",
      "likelihood:  8.035660499382475\n",
      "params:  [0.53645494 0.03097894 0.34746964 0.08509648 0.21404841 0.64097671]\n",
      "iteration:  4\n",
      "likelihood:  8.039677976940599\n",
      "params:  [0.53195115 0.03097217 0.34746964 0.08960704 0.20738666 0.63519072]\n",
      "iteration:  5\n",
      "likelihood:  8.042067007072502\n",
      "params:  [0.52933452 0.0306394  0.34746964 0.09255644 0.2035232  0.63427821]\n",
      "iteration:  6\n",
      "likelihood:  8.043465558903357\n",
      "params:  [0.52766141 0.03056049 0.34746964 0.09430846 0.20117095 0.63295979]\n",
      "iteration:  7\n",
      "likelihood:  8.044346235349673\n",
      "params:  [0.52661523 0.0304757  0.34746964 0.09543943 0.19970151 0.63240028]\n",
      "iteration:  8\n",
      "likelihood:  8.044899256848176\n",
      "params:  [0.52594139 0.03043516 0.34746964 0.09615381 0.19876841 0.63196588]\n",
      "iteration:  9\n",
      "likelihood:  8.045254305628687\n",
      "params:  [0.5255086  0.03040557 0.34746964 0.09661619 0.19817018 0.63171539]\n",
      "iteration:  10\n",
      "likelihood:  8.045482607543805\n",
      "params:  [0.52522809 0.03038791 0.34746964 0.09691436 0.19778423 0.63154583]\n",
      "iteration:  11\n",
      "likelihood:  8.045630466443669\n",
      "params:  [0.52504619 0.03037613 0.34746964 0.09710804 0.19753427 0.631439  ]\n",
      "iteration:  12\n",
      "likelihood:  8.045726370570064\n",
      "params:  [0.52492789 0.03036863 0.34746964 0.09723384 0.19737197 0.63136886]\n",
      "iteration:  13\n",
      "likelihood:  8.04578873175838\n",
      "params:  [0.52485091 0.03036373 0.34746964 0.09731573 0.19726641 0.63132357]\n",
      "iteration:  14\n",
      "likelihood:  8.045829316751355\n",
      "params:  [0.52480075 0.03036055 0.34746964 0.09736905 0.19719769 0.63129401]\n",
      "iteration:  15\n",
      "likelihood:  8.045855754342405\n",
      "params:  [0.52476807 0.03035848 0.34746964 0.09740381 0.19715292 0.63127479]\n",
      "iteration:  16\n",
      "likelihood:  8.04587298341155\n",
      "params:  [0.52474676 0.03035713 0.34746964 0.09742646 0.19712374 0.63126226]\n",
      "iteration:  17\n",
      "likelihood:  8.045884215500513\n",
      "params:  [0.52473287 0.03035625 0.34746964 0.09744123 0.19710472 0.63125409]\n",
      "iteration:  18\n",
      "likelihood:  8.045891539416887\n",
      "params:  [0.52472381 0.03035568 0.34746964 0.09745087 0.19709231 0.63124876]\n",
      "iteration:  19\n",
      "likelihood:  8.0458963157051\n",
      "params:  [0.5247179  0.03035531 0.34746964 0.09745715 0.19708422 0.63124529]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt0W+WZ7/HvY/keJ7GTgCEOJdyvLQFnKLSlQ0iB0mmhlFLCmg6sDpyUznQodDgHaM8Uhk57FlPans4A5TCFXqAllEB6KJ1COMUUhksgF3MJIZBAyMUmITi2cWLZlvWcP7TtKIpsy5JsOdq/z1pakvb7bu1HO8pP26/2xdwdEREJj5JCFyAiIuNLwS8iEjIKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCprTQBaQzY8YMnz17dlbz7ty5k0mTJuW3oDxSfblRfblRfbmZyPWtWLFiu7vvl1Fnd59wt8bGRs9WU1NT1vOOB9WXG9WXG9WXm4lcH7DcM8xYDfWIiISMgl9EJGQU/CIiIaPgFxEJGQW/iEjIKPhFREJGwS8iEjIT8gAukWLj7sQd+uNO3J1Y3OlPvbkTD9rjDvFgn+v+OME0x4PpidfZ/bqp7YP3DByrA3GHV7bGiL76LjAw/979nODe2T0dILkN9uy7+40OPvagvmDW3dOS1knSbACse7uPN0rWk9RE0sPB6c7e8w613vd8ntK+V/80r5HUa8OGXppjbwy5zBEvZDvCpW6rK0q54i8PG+lVcqbgl6Lg7kT74nT39dMT66enL05PLE60r5+eWJyeWD/RvviIbb2xOLF+py8ep6/fifUn7vv648TiSY8Hp+35vDvaQ8mflxKPB+Huu4N9wli1otAVDG/t64WuYHjr3sxpdrOh22bUVCj4pfi5O53RGG07e2nb2UtntI+dPbHg1s/Onhhdvbufb9gS5c43n2dnb/9gv66eGLt6+3MK1/JICRWlJZSXllAaMcoiJZRFSigtMUojJZRHEvelJUZNRWnS9ET/0pISyiLGtndbOWjWTCIlJURK2PPeLM00iER2t5WYESkxSswwSzwvCdoseFwyML1k9+OBvgaD0yFxb8E8hrFy5XLmzp271zwWPB98DMFzGwwqs/RtwaJItLB7WtJrDUxPdN3dIbXv008/zSdPO22P/oPzpLxOsqH6puufOruldEiXywNdnnzySU4//fS08+1LFPySVz2xfnbs7BsM8rZdvbR19dC2q48dA9OS2nbs7CU2QmCXl5ZQU1HKpIoI9Dn1VXFqq8poqK1kUnkpkypKg/ZSqssjVJSWUFFWQkVphMrgvqI05XlZIugryyKUR0ooKcnPf+Inn2zj9NOPz8trjYX310U4bubUQpcxpKpSY1LFxI2lxBfkvhv4AybuGpYJaVdvjM07utm8Yxeb2pLu2xP3Hd19aeczg9qqMuomlTN9UjkHT6/mpINrqasuZ9qkxK1uUjlTq8oGQ3xSeYRJFaWURXbvg5DY4vrYeL1dkaKk4Jc99PY7b73XxaakcN+0Y1ci7Nt28f7O3j36V5SWMKuuioOmVTPnoFoOmFI5GO511eVMr0nc11aXE8nTVrWI5Caj4Dezq4HLSfxo/QrwFXePJrVXAL8CGoH3gYvcfUNS+4eA14Ab3f2WvFUvOYn29bO6pZPmTe2s2riDlza3s6mtGx7/82CfsojRUFvFrLpqzjqunll11YNBP6uuiv1qKoriT1+RMBkx+M2sAbgSONbdu83st8AC4BdJ3S4Ddrj74Wa2ALgZuCip/cfAH/NWtYyau7OxbRerNrYPBv1rrZ309SfG12dOreTED9Uxd1qM0046lll11Rw0rYr9J1dqS12kyGQ61FMKVJlZH1ANtKS0nwfcGDxeDNxqZububmafB94CduahXslQR3cfL29uHwz65k3ttAXDNFVlET4yayqXfeJQ5hxUy4kfqqV+SiUQjKGfNKuQpYvIGBsx+N19i5ndAmwEuoGl7r40pVsDsCnoHzOzDmC6mXUD1wJnAtfktXLZwwfRPh5bvZXn33qf5k3trNvWNdh2+P41zD96f+Z8qJYTD6rjyPoaSiM6aFskrCz1yLa9OpjVAQ+SGLppBx4AFrv7vUl9VgNnu/vm4Pl64GTgeuAFd/+tmd0IdA01xm9mC4GFAPX19Y2LFi3K6g11dXVRU1OT1bzjIZ/19ced1e/388yWGCu39dMXh8llcGhthMNqSzh0aoRDppYwqSzzoZowrb+xoPpyo/qyN2/evBXuPjejziNdogu4ELgr6fklwO0pfR4DTg0elwLbSRwH8TSwIbi1A23A10dapi69OLR4PO6vbG73m36/2hu/+7gffO0j/pEbH/NvL3nZl29o83g8XtD6xprqy43qy81Ero9RXHoxkzH+jcApZlZNYqhnPrA8pc/DwKXAc8AXgSeCQk4b6JC0xX9rRt9Isod3O6L8rnkLD63czBtbuyiLGGccvT/nnziLeUfvR0VppNAlisg+IpMx/mVmthhYCcSAVcCdZnYTiW+Yh4G7gHvMbB2JrfoFY1hzaOzsifHoq+/y0KrNPLv+fdzhpA/V8t3PH89nP3wgdZPKC12iiOyDMtqrx91vAG5ImfydpPYoiSGh4V7jxtEWF0b9ceeZddtZsmoLj776Lt19/Rw0rYp/OOMIzj+xgUNmTCp0iSKyj9ORuxNErD/ObU3r+fWyd9j2QQ9TKkv5/IkNXHBSA40H1+kgKRHJGwX/BLC9q4ev/2Ylz7/VxhlH78+FjbOYd/T+VJZp3F5E8k/BX2AvbWrnintX0Lazlx9eeAIXNOrgKREZWwr+Arr/xY380+9Ws9/kCh782sc4vmHini5XRIqHgr8AemL9/PPvX+M3yzZy2hEz+LcFJ2oPHREZNwr+cbYjGmfBnc+zamM7Xzv9MK456yidBE1ExpWCfxy98HYbNzwbJUYvP/3rkzjnwwcWuiQRCSEF/zhwd37x7Aa+94c1zKiEe776cY6on1zoskQkpBT8Y6y7t59vL3mFh1Zt4VPH1POFmR8o9EWkoHRu3jG0qW0XF/z0WZY0b+EfzzySO/+mkepRnClTRGQsaIt/jDz1xntcuWgV8bhz96V/wbyj9y90SSIigII/79yd259czy1L13JU/WTu+HIjs3V+HRGZQBT8edTVE+Oa377Eo6vf5XMnzOTmCz5MdblWsYhMLEqlPIn1x7n4zud5rbWT//lXx3DZJw7RidVEZEJS8OfJklVbeGVLBz9ZMIfz5jQUuhwRkSFpr5486In187//35t8ZNZUzj1hZqHLEREZloI/Dxa9sIkt7d3897OP0vCOiEx4Cv4c7eqN8e9PrOOUQ6fxicNnFLocEZERKfhz9PNnNrC9q0db+yKyz1Dw56BjVx//58/rmX/0/jQePK3Q5YiIZCSj4Dezq81stZm9amb3mVllSnuFmd1vZuvMbJmZzQ6mn2lmK8zsleD+jPy/hcK58+n1dEZj/ONZRxW6FBGRjI0Y/GbWAFwJzHX344EIsCCl22XADnc/HPgxcHMwfTvwOXf/MHApcE++Ci+09z7o4e7/2sC5J8zk2JlTCl2OiEjGMh3qKQWqzKwUqAZaUtrPA34ZPF4MzDczc/dV7j7QdzVQaWYVuRY9EdzWtI7e/jhXn3lkoUsRERmVEYPf3bcAtwAbgVagw92XpnRrADYF/WNABzA9pc8FwCp378m16ELbvGMXv1m2kS/NncUhOg+PiOxjzN2H72BWBzwIXAS0Aw8Ai9393qQ+q4Gz3X1z8Hw9cLK7vx88Pw54GDjL3dcPsZyFwEKA+vr6xkWLFmX1hrq6uqipqclq3kzd9UoPz7XGuPm0KqZXje738fGoLxeqLzeqLzeqL3vz5s1b4e5zM+rs7sPegAuBu5KeXwLcntLnMeDU4HEpibH9gS+VWcAbwMdHWtbArbGx0bPV1NSU9byZWLftAz/kukf8pt+vzmr+sa4vV6ovN6ovN6ove8ByzzBjM9lc3QicYmbVlthRfT6wJqXPwyR+vAX4IvCEu7uZ1QJ/AK5392cy+iaa4H70+BtUlUX4u9MPK3QpIiJZyWSMfxmJH2xXAq8E89xpZjeZ2blBt7uA6Wa2DvgmcF0w/evA4cA/mVlzcNtnr0jy6pYO/vByK5d94hCm1xTFb9QiEkIZnZ3T3W8AbkiZ/J2k9iiJIaHU+f4F+JdcCpxIfrh0LVOryrj8k4cWuhQRkazpyN0Mvbihjaa17/G10w9jSmVZocsREcmagj8D7s6/Pvo6+02u4NJTZxe6HBGRnCj4M/DnN97jxQ07uPKMw6kqjxS6HBGRnCj4RxCPOz94bC2z6qq46C8+VOhyRERypuAfwaOr32V1SydXf+pIyku1ukRk36ckG0asP84Pl67l8P1r+PyJuo6uiBQHBf8wlqzawvr3dnLNWUcSKdFFVkSkOCj4hzBwAfUPN0zl7OMOKHQ5IiJ5o+Afgi6gLiLFSsGfxsAF1D96yDROO0IXUBeR4qLgT+MXzyYuoP4/Pq2tfREpPgr+FB3dfdzxpC6gLiLFS8Gf4s6ndAF1ESluCv4kAxdQ/5wuoC4iRUzBn2TwAuqfOqLQpYiIjBkFf8DdeWD5Js49YSaH7jcxr6kpIpIPCv5A+64+dvb2c5yGeESkyCn4Ay0d3QDMrK0qcCUiImNLwR9obY8CCn4RKX4K/sDgFv/UygJXIiIytjIKfjO72sxWm9mrZnafmVWmtFeY2f1mts7MlpnZ7KS264Ppa83s7PyWnz8t7VHKIsaMmopClyIiMqZGDH4zawCuBOa6+/FABFiQ0u0yYIe7Hw78GLg5mPfYoO9xwKeB281sQl67sKW9mwOmVlKi0y+LSJHLdKinFKgys1KgGmhJaT8P+GXweDEw3xInuTkPWOTuPe7+NrAOODn3svOvtaObA6dqfF9Eit+Iwe/uW4BbgI1AK9Dh7ktTujUAm4L+MaADmJ48PbA5mDbhtLRHadAPuyISAqUjdTCzOhJb7ocA7cADZvZld783uVuaWX2Y6emWsxBYCFBfX8+TTz45UmlpdXV1jXreuDutHd30dWzLermZyqa+8aT6cqP6cqP6xom7D3sDLgTuSnp+CXB7Sp/HgFODx6XAdhKhfz1wfbp+w90aGxs9W01NTaOep7W92w++9hH/1XMbsl5uprKpbzypvtyovtyovuwBy32EbB24ZTLGvxE4xcyqg3H7+cCalD4PA5cGj78IPBEU8jCwINjr5xDgCOCFLL6fxtSW9sSunA212pVTRIrfiEM97r7MzBYDK4EYsAq408xuIvEN8zBwF3CPma0D2gj2+nH31Wb2W+C1YN6/d/f+sXkr2WsN9uHXj7siEgYjBj+Au98A3JAy+TtJ7VESQ0Lp5v0e8L1sCxwPOmpXRMJER+6SGOqZVB5hSmVG34MiIvs0BT/BPvy1Vbq+roiEgoKfxD78GuYRkbBQ8JPY4tfJ2UQkLEIf/NG+frZ39WqLX0RCI/TB/25HYo+eA7XFLyIhEfrg15W3RCRsFPzah19EQib0wd/aPnDUroZ6RCQcQh/8LR1Rpk8qp7JsQl4fRkQk7xT87d0cqJOziUiIhD74deUtEQmb0Ae/rrwlImET6uDvjPbR1RPTD7siEiqhDn6djllEwijUwd/SPnDwlrb4RSQ8wh38uvKWiIRQuIO/vZtIibH/5IpClyIiMm5CHfyt7VHqJ1dQGgn1ahCRkAl14rV0dOuHXREJnRGD38yOMrPmpFunmV2V0qfOzJaY2ctm9oKZHZ/UdrWZrTazV83sPjObML+ktrRHOVDBLyIhM2Lwu/tad5/j7nOARmAXsCSl27eAZnf/CHAJ8BMAM2sArgTmuvvxQARYkMf6sxaPO+92RHXlLREJndEO9cwH1rv7OynTjwX+BODurwOzzaw+aCsFqsysFKgGWnKoN2+27+yhtz+uoR4RCZ3RBv8C4L40018CvgBgZicDBwOz3H0LcAuwEWgFOtx9afbl5s/AwVs6aldEwsbcPbOOZuUkttaPc/etKW1TSAzvnAi8AhwNXE4i8B8ELgLagQeAxe5+b5rXXwgsBKivr29ctGhRVm+oq6uLmpqaEfstfzfGrc09/PPHKjl4yvidkjnT+gpF9eVG9eVG9WVv3rx5K9x9bkad3T2jG3AesDSDfgZsAKYAFwJ3JbVdAtw+0ms0NjZ6tpqamjLq97On3/KDr33E3+/qyXpZ2ci0vkJRfblRfblRfdkDlnuGeT6aoZ6LST/Mg5nVBn8RQGJL/yl37ySxxX+KmVWbmZH4jWDNKJY5Zlrbu6koLaGuuqzQpYiIjKuMgt/MqoEzgYeSpl1hZlcET48BVpvZ68A5wDcA3H0ZsBhYSWIIqAS4M2/V56Clo5uG2ioS30ciIuFRmkknd98FTE+ZdkfS4+eAI4aY9wbghhxqHBOJffj1w66IhE9oj9xt7ehmpk7OJiIhFMrg743F2fZBj47aFZFQCmXwb+2M4o6O2hWRUApl8Ld26MpbIhJeoQx+XXlLRMIsnMGvK2+JSIiFM/jbu5laVcakioz2ZhURKSqhDP7W9qhOziYioRXK4G/piNKgH3ZFJKTCGfzt3TpqV0RCK3TBv7MnRkd3n3blFJHQCl3wtwZ79Oh0DSISVqEL/hZdeUtEQi50wT+4xa+hHhEJqdAF/5b2KGZwgLb4RSSkQhf8re3d7FdTQVkkdG9dRAQIYfC3dHRrmEdEQi10wd/aHtXJ2UQk1EIV/O6e2OLXrpwiEmKhCv4du/qI9sV15S0RCbVQBf/gefi1R4+IhNiIwW9mR5lZc9Kt08yuSulTZ2ZLzOxlM3vBzI5Paqs1s8Vm9rqZrTGzU8fijWRi9wVYtMUvIuE14gnp3X0tMAfAzCLAFmBJSrdvAc3ufr6ZHQ3cBswP2n4CPOruXzSzcqA6X8WP1sAlF3WCNhEJs9EO9cwH1rv7OynTjwX+BODurwOzzazezKYAnwTuCtp63b09x5qz1tLRTXmkhBmTKgpVgohIwZm7Z97Z7G5gpbvfmjL9+0Clu3/TzE4GngU+CvQDdwKvAScAK4BvuPvONK+9EFgIUF9f37ho0aKs3lBXVxc1NTVp237aHOWtjjg/+MuC/dExbH0TgerLjerLjerL3rx581a4+9yMOrt7RjegHNgO1KdpmwL8HGgG7gFeJBH0c4EY8NGg30+A7460rMbGRs9WU1PTkG0X3P6Mf+mOZ7N+7XwYrr6JQPXlRvXlRvVlD1juGeb5aC46ew6Jrf2tab48OoGvAJiZAW8Ht2pgs7svC7ouBq4bxTLzqqW9m1MOnV6oxYuITAijGeO/GLgvXUOw50558PRy4Cl373T3d4FNZnZU0DafxLDPuOuPO1s/6NEPuyISehlt8ZtZNXAm8NWkaVcAuPsdwDHAr8ysn0SwX5Y0+z8Avw6+GN4i+MtgvG37IEp/3LUrp4iEXkbB7+67gOkp0+5IevwccMQQ8zaTGOsvqN0Hbyn4RSTcQnPk7uCVtzTUIyIhF6Lg11G7IiIQouBv7YhSU1HKlMqyQpciIlJQoQn+lvZunYdfRIQwBX9HNwfqh10RkfAEv668JSKSEIrgj/b18/7OXu3KKSJCSIJ/9+mYFfwiIuEI/sFdOTXUIyISiuDfoqN2RUQGhSL4B4Z6DtC1dkVEwhL83cyoKaeyLFLoUkRECi4Uwb+lPap9+EVEAqEI/lYdtSsiMqjog9/daWnXUbsiIgOKPvg7ozF29vZri19EJFD0wd/aodMxi4gkK/rgHzgPv4Z6REQSQhD8iX34G7TFLyIChCL4uyktMfabXFHoUkREJoQRg9/MjjKz5qRbp5ldldKnzsyWmNnLZvaCmR2f0h4xs1Vm9ki+38BIWjui1E+pJFJi471oEZEJqXSkDu6+FpgDiQAHtgBLUrp9C2h29/PN7GjgNmB+Uvs3gDXAlHwUPRq68paIyJ5GO9QzH1jv7u+kTD8W+BOAu78OzDazegAzmwX8FfCzHGvNiq68JSKyJ3P3zDub3Q2sdPdbU6Z/H6h092+a2cnAs8BH3X2FmS0G/hcwGbjG3T87xGsvBBYC1NfXNy5atCirN9TV1UVNTQ0AcXf+29JdnD27jC8dVZ7V6+Vbcn0TkerLjerLjerL3rx581a4+9yMOrt7RjegHNgO1KdpmwL8HGgG7gFeBE4APgvcHvQ5HXgkk2U1NjZ6tpqamgYfb+3o9oOvfcR/+ezbWb9eviXXNxGpvtyovtyovuwByz3DPB9xjD/JOSS29rem+fLoBL4CYGYGvB3cFgDnmtlngEpgipnd6+5fHsVys9YycOUtDfWIiAwazRj/xcB96RrMrNbMBsZSLgeecvdOd7/e3We5+2wSXwJPjFfog668JSKSTkZb/GZWDZwJfDVp2hUA7n4HcAzwKzPrB14DLst/qaOnK2+JiOwto+B3913A9JRpdyQ9fg44YoTXeBJ4ctQV5qC1I0pVWYTa6rLxXKyIyIRW1EfutrR3c2BtJYmfHUREBIo9+DuiGuYREUlR1MGvK2+JiOytaIO/Nxbnva4e7copIpKiaIN/a2cUd+3KKSKSqmiDf3BXTp2HX0RkD0Ub/AOXXNRQj4jInoo2+AeuvKWhHhGRPRVx8HdTW11GdfloTkckIlL8ijb4WzuiGuYREUmjaIO/pb2bBg3ziIjspaiDX1v8IiJ7K8rg7+qJ0RmNaVdOEZE0ijL4dR5+EZGhFWXw68pbIiJDK87g1xa/iMiQijL4W9u7MYP6KQp+EZFURRn8LR1R6idXUhYpyrcnIpKTokzGgStviYjI3ooy+Ft15S0RkSGNGPxmdpSZNSfdOs3sqpQ+dWa2xMxeNrMXzOz4YPpBZtZkZmvMbLWZfWOs3sgAd6dFV94SERnSiGcwc/e1wBwAM4sAW4AlKd2+BTS7+/lmdjRwGzAfiAH/6O4rzWwysMLMHnf31/L5JpJ90Ac9sbh25RQRGcJoh3rmA+vd/Z2U6ccCfwJw99eB2WZW7+6t7r4ymP4BsAZoyLHmYbV1xwFdgEVEZCjm7pl3NrsbWOnut6ZM/z5Q6e7fNLOTgWeBj7r7iqQ+s4GngOPdvTPNay8EFgLU19c3Llq0aPTvBnjmnS7+Y41xw6mVHDI1ktVrjKWuri5qamoKXcaQVF9uVF9uVF/25s2bt8Ld52bU2d0zugHlwHagPk3bFODnQDNwD/AicEJSew2wAvhCJstqbGz0bH3rF0v94Gsf8W2d0axfYyw1NTUVuoRhqb7cqL7cqL7sAcs9wzwfzVVKziGxtb81zZdHJ/AVADMz4O3ghpmVAQ8Cv3b3h0axvKy0RZ3y0hKmTyof60WJiOyTRjPGfzFwX7oGM6s1s4GkvRx4yt07gy+Bu4A17v6j3ErNzPvdcQ6cWklJiY3H4kRE9jkZBb+ZVQNnAg8lTbvCzK4Inh4DrDaz10n8ZTCw2+bHgb8BzkjaHfQzeas+jbaoax9+EZFhZDTU4+67gOkp0+5IevwccESa+f4LGNdN77ao8xHtwy8iMqSiOnI31h9nh7b4RUSGVVTBv+2DHhztwy8iMpyiCv6B8/DrBG0iIkMrruAPrrzVoC1+EZEhFVfwD2zxT9UWv4jIUIoq+Fvbu6kqhcmVZYUuRURkwiqq4G/piDK9UgduiYgMp7iCv72baZVF9ZZERPKuqFKytSPKtCpt8YuIDGc0J2mb0OJx5/Qj92NG//ZClyIiMqEVzRZ/SYnxo4vm8LGZRfNdJiIyJoom+EVEJDMKfhGRkFHwi4iEjIJfRCRkFPwiIiGj4BcRCRkFv4hIyCj4RURCxty90DXsxczeA97JcvYZwEQ+fFf15Ub15Ub15WYi13ewu++XSccJGfy5MLPl7j630HUMRfXlRvXlRvXlZqLXlykN9YiIhIyCX0QkZIox+O8sdAEjUH25UX25UX25mej1ZaToxvhFRGR4xbjFLyIiw9hng9/MPm1ma81snZldl6a9wszuD9qXmdnscaztIDNrMrM1ZrbazL6Rps/pZtZhZs3B7TvjVV+w/A1m9kqw7OVp2s3M/i1Yfy+b2UnjWNtRSeul2cw6zeyqlD7juv7M7G4z22ZmryZNm2Zmj5vZm8F93RDzXhr0edPMLh3H+n5gZq8H/35LzKx2iHmH/SyMYX03mtmWpH/Dzwwx77D/18ewvvuTattgZs1DzDvm6y/v3H2fuwERYD1wKFAOvAQcm9Ln74A7gscLgPvHsb4DgZOCx5OBN9LUdzrwSAHX4QZgxjDtnwH+CBhwCrCsgP/W75LYR7lg6w/4JHAS8GrStH8FrgseXwfcnGa+acBbwX1d8LhunOo7CygNHt+crr5MPgtjWN+NwDUZ/PsP+399rOpLaf8h8J1Crb983/bVLf6TgXXu/pa79wKLgPNS+pwH/DJ4vBiYb2bjckFed29195XB4w+ANUDDeCw7j84DfuUJzwO1ZnZgAeqYD6x392wP6MsLd38KaEuZnPwZ+yXw+TSzng087u5t7r4DeBz49HjU5+5L3T0WPH0emJXv5WZqiPWXiUz+r+dsuPqC3PgScF++l1so+2rwNwCbkp5vZu9gHewTfPg7gOnjUl2SYIjpRGBZmuZTzewlM/ujmR03roWBA0vNbIWZLUzTnsk6Hg8LGPo/XCHXH0C9u7dC4sse2D9Nn4myHv+WxF9w6Yz0WRhLXw+Gou4eYqhsIqy/04Ct7v7mEO2FXH9Z2VeDP92We+ruSZn0GVNmVgM8CFzl7p0pzStJDF+cAPw78LvxrA34uLufBJwD/L2ZfTKlfSKsv3LgXOCBNM2FXn+Zmgjr8dtADPj1EF1G+iyMlZ8ChwFzgFYSwympCr7+gIsZfmu/UOsva/tq8G8GDkp6PgtoGaqPmZUCU8nuT82smFkZidD/tbs/lNru7p3u3hU8/k+gzMxmjFd97t4S3G8DlpD4kzpZJut4rJ0DrHT3rakNhV5/ga0Dw1/B/bY0fQq6HoMfkz8L/LUHA9KpMvgsjAl33+ru/e4eB/5jiOUWev2VAl8A7h+qT6HWXy721eB/ETjCzA4JtgoXAA+n9HkYGNiD4ovAE0N98PMtGBO8C1jj7j8GSNQPAAABZ0lEQVQaos8BA785mNnJJP4t3h+n+iaZ2eSBxyR+BHw1pdvDwCXB3j2nAB0DwxrjaMgtrUKuvyTJn7FLgf+bps9jwFlmVhcMZZwVTBtzZvZp4FrgXHffNUSfTD4LY1Vf8m9G5w+x3Ez+r4+lTwGvu/vmdI2FXH85KfSvy9neSOx18gaJX/y/HUy7icSHHKCSxBDBOuAF4NBxrO0TJP4cfRloDm6fAa4Argj6fB1YTWIvheeBj41jfYcGy30pqGFg/SXXZ8Btwfp9BZg7zv++1SSCfGrStIKtPxJfQK1AH4mt0MtI/Gb0J+DN4H5a0Hcu8LOkef82+ByuA74yjvWtIzE+PvAZHNjLbSbwn8N9FsapvnuCz9bLJML8wNT6gud7/V8fj/qC6b8Y+Mwl9R339Zfvm47cFREJmX11qEdERLKk4BcRCRkFv4hIyCj4RURCRsEvIhIyCn4RkZBR8IuIhIyCX0QkZP4/ja18iPUu+fIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from maximumlikelihood import ExpectationMaximization as EM\n",
    "theta1, likelihood_log1, fig1, ax1 = EM(theta.copy(), data.copy(), gt1.copy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ground truth 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iteration:  0\n",
      "likelihood:  8.6426753215856\n",
      "params:  [0.45039944 0.13601098 0.34746964 0.06611994 0.11043395 1.36137458]\n",
      "iteration:  1\n",
      "likelihood:  10.765529369611052\n",
      "params:  [0.56336279 0.07168543 0.34746964 0.01748214 0.04050285 1.75811821]\n",
      "iteration:  2\n",
      "likelihood:  12.527122795827342\n",
      "params:  [0.58099373 0.04593928 0.34746964 0.02559735 0.02558128 1.0515376 ]\n",
      "iteration:  3\n",
      "likelihood:  14.47923998843127\n",
      "params:  [0.57629668 0.03171166 0.34746964 0.04452202 0.01740276 1.12667585]\n",
      "iteration:  4\n",
      "likelihood:  16.50546274737689\n",
      "params:  [0.56716026 0.02875284 0.34746964 0.05661726 0.01271706 0.7477079 ]\n",
      "iteration:  5\n",
      "likelihood:  18.31106197763992\n",
      "params:  [0.55942098 0.02020952 0.34746964 0.07289987 0.01009555 0.94112917]\n",
      "iteration:  6\n",
      "likelihood:  20.120244846864004\n",
      "params:  [0.54886275 0.02405653 0.34746964 0.07961108 0.0082119  0.60580942]\n",
      "iteration:  7\n",
      "likelihood:  21.645393899749298\n",
      "params:  [0.54233279 0.01583485 0.34746964 0.09436272 0.00705399 0.89791961]\n",
      "iteration:  8\n",
      "likelihood:  22.908790439621747\n",
      "params:  [0.53423608 0.0245455  0.34746964 0.09374878 0.00624256 0.53949476]\n",
      "iteration:  9\n",
      "likelihood:  23.812982912676272\n",
      "params:  [0.53040221 0.01445691 0.34746964 0.10767124 0.00576762 0.88350728]\n",
      "iteration:  10\n",
      "likelihood:  24.65861462717362\n",
      "params:  [0.52425243 0.02613996 0.34746964 0.10213797 0.00534173 0.51748248]\n",
      "iteration:  11\n",
      "likelihood:  25.240954895958794\n",
      "params:  [0.52220501 0.01448143 0.34746964 0.11584392 0.00509785 0.8571942 ]\n",
      "iteration:  12\n",
      "likelihood:  25.839812795539615\n",
      "params:  [0.51759636 0.02690605 0.34746964 0.10802795 0.00483748 0.5167883 ]\n",
      "iteration:  13\n",
      "likelihood:  26.208597149931055\n",
      "params:  [0.51663095 0.01508779 0.34746964 0.12081162 0.00470542 0.82158066]\n",
      "iteration:  14\n",
      "likelihood:  26.62519309021752\n",
      "params:  [0.51321461 0.02683189 0.34746964 0.11248387 0.00453919 0.52615736]\n",
      "iteration:  15\n",
      "likelihood:  26.85384638935868\n",
      "params:  [0.51287132 0.01593357 0.34746964 0.12372547 0.00446644 0.78531062]\n",
      "iteration:  16\n",
      "likelihood:  27.146037219982283\n",
      "params:  [0.51034366 0.0262935  0.34746964 0.1158932  0.00435566 0.53958039]\n",
      "iteration:  17\n",
      "likelihood:  27.28812478845914\n",
      "params:  [0.51032402 0.01683564 0.34746964 0.1253707  0.00431477 0.75300503]\n",
      "iteration:  18\n",
      "likelihood:  27.494841109833313\n",
      "params:  [0.50845189 0.02558632 0.34746964 0.11849215 0.00423874 0.55391206]\n",
      "iteration:  19\n",
      "likelihood:  27.5828873402268\n",
      "params:  [0.50858697 0.01769387 0.34746964 0.12624952 0.00421572 0.72616816]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VPXZ//H3TcKWhD0Q9k2QKihLAqJWa6oiRa1L3ahatxat2tan2l+1rY9eWqutddfHaiuiFY1atSrighgX6haC7FsQEFlkC0tCCCHJ/fsjg41xQoaZZGYy83ld11xzlu8385mTyT0n3zlzjrk7IiKSPFrEOoCIiESXCr+ISJJR4RcRSTIq/CIiSUaFX0Qkyajwi4gkGRV+EZEko8IvIpJkVPhFRJJMaqwDBJOZmen9+/cPq++uXbtIT09v3ECNSPkio3yRUb7IxHO+wsLCLe7eNaTG7h53t+zsbA9Xfn5+2H2jQfkio3yRUb7IxHM+YLaHWGM11CMikmRU+EVEkowKv4hIklHhFxFJMir8IiJJpsHDOc2sD/Ak0B2oBh519/vM7FlgSKBZR2C7u48I0n81UAJUAZXuntNI2UVEJAyhHMdfCVzr7nPMrB1QaGYz3P3cfQ3M7C5gx35+Rq67b4kwq4iINIIGC7+7bwA2BKZLzGwJ0AtYDGBmBpwDfL8Jc4qINGuVVdWU7qmkpLySneV7KS2vmS7ZUzO9s7ySlBbGFd87qMmzHNA3d82sPzAS+KTW4mOAje5eVE83B94yMwcecfdHw8gpIhJXqqudjSXlrNlaxhfFZXxZXMbGneWUlFdSuqemkJfUKvC791Y1+DO7tmsdlcJvHuLF1s0sA3gPuM3dX6y1/GFghbvfVU+/nu6+3sy6ATOAX7j7+0HaTQImAWRlZWXn5eUd8JMBKC0tJSMjI6y+0aB8kVG+yCjfgamocjbvdjaVVbO5zFm3Yw/bK1Nr5nc7ldX/bWtAh9ZGWiq0TbWaW8ua6W8sC0yntfzv9L7lrVIs7Ky5ubmFoX6GGlLhN7OWwDTgTXe/u9byVGAdkO3ua0P4OTcDpe7+1/21y8nJ8dmzZzeYK5h3332X4447Lqy+0aB8kVG+yCjff1VWVbN1VwWbS/awuXQPm0v2sGF7OWuKy1hTvIs1xWVs3LnnG33apMCAbu3p27kt/bqk06dzGv06p9G3cxq9OrWlZUrsDpQ0s5ALfyhH9RjwGLCkdtEPOAFYWl/RN7N0oEXgs4F0YBxwSyjBREQOVHW1s62s4utCvrlkD1u+Mf3fQr+trIJg+73d27ehb5c0jhnclb6d0+jXJe3rAj+/4ENyc4+J/hNrZKGM8R8NXAgsMLO5gWW/c/fpwHnAM7Ubm1lP4B/uPgHIAl6qee8gFXja3d9orPAikrw2l+xh8YadLFq/g0Xrd7Jk/U6+KC6jqvrb1bx1agu6tmtN13at6dcljez+neiaUTOfGbjvFljfpmVKvY8ZqGXNXihH9cyiZvgq2LqLgyxbD0wITK8EhkcWUUSSmbuzpriMRet3snj9fwv9ppL/DsP07tSWoT3bM35Y90ABb0NmRquvi31G69SEKdqNIS7Pxy8iyamy2r9R3BdvqNmTL9lTCUBKC2NQ1wy+OyiTQ3u2Z2jPDhzaoz0d0lrGOHnzosIvIlHl7mwr28uqLaWs3LyLVVv+e1uxsYzKtz4AoG3LFL7Tox2njezJ0J4dGNqzPQdntdvvUIyERoVfRJpEWUUlq7eUBYp6KSu37Pq60O/YvffrdqktjL5d0hiYmc6AtuWMP2IoQ3t2YEBmOiktNDzTFFT4RSRi67bv5oPlm1mwbsfXe+8bdpR/o02PDm0YkJnOKYf3YEBmOgO7pjMgM4PetQ6DfPfddzluRK9YPIWkosIvIgespHwvH68sZlbRZj4o2sLKLbsAaN8mlYFdMzjyoC41e/CZGQzITKd/ZhpprVRu4oV+EyLSoMqqauav28Gsoi18ULSZz9Zsp7LaadsyhbEDO3P+2H4cOziTQd0ydPRMM6DCLyJBrdlaxvtFm5lVtIUPP9/CzvJKzOCwXh2YdOxAjhnclVH9OtI6VR+2Njcq/CICwI6yvXy0cgsfFNXc1hSXAdCrY1t+MKwHxxycydEHZdIpvVWMk0qkVPhFktSmneV8urqYT1fV3JZtLMEd0lulcORBXbjsuwP47uBMBmama/gmwajwiySBfd9+/WDtXl57fh4Fq4tZvbVmjz6tVQrZ/Tpx8mE9OGJgF0b27RjTk41J01PhF0lA1dVO0aZSPl21lU9WFVOwuvjrM012TNvI6P6duWBsP0b378zQnu1JVaFPKir8IgmiaGMJ+cs28emqbRSsLv76S1Ld27fhiAFdGDOgM7b5cyaenEsLfTEqqanwizRju/ZUMm3+evIKvuSzNdsBGJCZzvih3RkzoDNjBnSmd6e2X4/Rv/vuKhV9UeEXaW7cnc++3M5zBV/y6rz17KqoYnC3DP5w8iH8cHhPurVvE+uIEudU+EWaieJdFbz02TqeLVjD8o2ltG2ZwqnDe3Du6L6M6ttRR95IyFT4ReJYdbXz4edbyStYw1uLNlJRVc3wPh25/czDOOXwHrRro9MRy4FT4ReJQxt27Ob52Wt5bvaXrN22m45pLTl/bF/OHd2H73RvH+t40syp8IvEib1V1cxcsolnC9bw3vLNVDscPagL/2/8dxh3aJbOQy+NJpSLrfcBngS6A9XAo+5+n5ndDPwM2Bxouu86vHX7jwfuA1KouRbvHY2UXSQhLPuqhOdnf8m/565jS2kFWe1bc+Vxgzgnpw99u6TFOp4koFD2+CuBa919jpm1AwrNbEZg3T3u/tf6OppZCvAQcCKwFigws1fcfXGkwUWasx1le3ll3jqeL1zL/LU7SG1hHH9IN87O7sNxQ7rqC1XSpEK52PoGYENgusTMlgChXilhDLAicNF1zCwPOA1Q4ZekU1Xt/GfFFp4vXMubi76iorKa73Rvx42nHMrpI3rSJaN1rCNKkjB3D72xWX/gfWAY8GvgYmAnMJua/wq21Wl/FjDe3X8amL8QOMLdrw7ysycBkwCysrKy8/LyDvzZAKWlpWRkZITVNxqULzLNMd/GXdXMWlfJf9ZXUlzupLeEsT1SOaZXKv3at4jqYZjNcfvFk3jOl5ubW+juOSE1dveQbkAGUAicGZjPombcvgVwGzA5SJ+zqRnX3zd/IfBAQ4+VnZ3t4crPzw+7bzQoX2SaS77S8r3+bMEaP/vhD73fb6f5gOun+U8e+8RfnbfOd1dUxjxfvFK+8AGzPcR6HtJRPWbWEngBmOruLwbeMDbWWv93YFqQrmuBPrXmewPrQ3pHEmlm3J1lxVW89vw8XluwgbKKKgZkpvObk4bwo1G96d5B36iV+BDKUT0GPAYscfe7ay3v4TXj/wBnAAuDdC8ABpvZAGAdcB7w44hTi8QRd+e95Zu54/WlLP2qnPRWGzj18J6cndOb7H6d9I1aiTuh7PEfTc0QzQIzmxtY9jtgopmNABxYDVwOYGY9qRnemeDulWZ2NfAmNcNCk919USM/B5GYWbR+B7dPX8qsFVvo2zmNy4a14tpzcnVhcYlroRzVMwsItsvyrWP2A+3XAxNqzU+vr61Ic7V++27++tYyXvpsHR3atuTGUw7lgrF9+WjWByr6Evf0ChU5ADvL9/Lwu58zedYqHJh0zECuzB1Eh7Y6Z440Hyr8IiHYW1XN05+s4b6ZRRTvquD0ET257qQh9O6kb9ZK86PCL7If7s6bi77iz28sY9WWXRw5sAu/m3AIh/XuEOtoImFT4RepR+EX2/jT9CUUfrGNwd0ymHxxDrlDuukoHWn2VPhF6li9ZRd/eXMp0xd8Rdd2rbn9zMM4O7u3zp8jCUOFXySgeFcF988sYuonX5DaogXXnDCYnx0zkPTW+jORxKJXtCS9yqpqpn6yhrveWkbpnkrOHd2H/znhYF27VhKWCr8ktY9XbuXmVxax9KsSjh7UhZtOHcrBWe1iHUukSanwS1LasGM3t722hGnzN9CrY1v+dsEoThraXR/cSlJQ4ZekUr63isdmreLBd1ZQ7c6vjh/MFd87iLatdFlDSR4q/JIU3J2ZSzZx62uL+WJrGeOHduf3Jx9Cn876ApYkHxV+SXgrN5dyy7TFvLtsMwd1Teefl43hmMFdYx1LJGZU+CVhle6p5MF3VvDYrJW0Tk3hDycfwkVH9aeljseXJKfCLwnH3Xl57npuf30JG3fu4azs3vy/8UPo1k6HZ4qACr8kmEXrd3DzK4soWL2Nw3t34OELshnVt1OsY4nEFRV+SQh7q6q54/WlPP6fVXRMa8UdZx7GOTl9aNFCh2eK1KXCL83etl0VXDl1Dh+t3MoFY/vym3HfoUOazo8vUp9QrrnbB3gS6A5UA4+6+31mdidwKlABfA5c4u7bg/RfDZQAVUClu+c0XnxJdss3lnDZEwVs3LGHu84ezo+ye8c6kkjcC+XwhkrgWnc/BBgLXGVmhwIzgGHufjiwHLhhPz8j191HqOhLY3p78UbOeOg/lO+tJu/ysSr6IiEK5Zq7G4ANgekSM1sC9HL3t2o1+xg4q2kiinyTu/Pwe59z55vLGNazA4/+JJseHdrGOpZIs2HuHnpjs/7A+9Ts6e+stfxV4Fl3fypIn1XANsCBR9z90Xp+9iRgEkBWVlZ2Xl5e6M+iltLSUjIyMsLqGw3KF5niHaU8tzqVjzdUcUT3FC49rDWtU+LnA9x4337KF5l4zpebm1sY8qiKu4d0AzKAQuDMOst/D7xE4E0kSL+egftuwDzg2IYeKzs728OVn58fdt9oUL7wbdi+24/703Tv99tp/sDM5V5dXR3rSN8Sz9vPXfkiFc/5gNkeYj0P6ageM2sJvABMdfcXay2/CDgFOD7wwMHeWNYH7jeZ2UvAGGr+axAJ2dwvtzPpydnsKKvm0QuzGTe0e6wjiTRbDX64azXnqX0MWOLud9daPh74LfBDdy+rp2+6mbXbNw2MAxY2RnBJHv/+bB3nPPIRrVJb8IexbVX0RSIUylE9RwMXAt83s7mB2wTgQaAdMCOw7G8AZtbTzKYH+mYBs8xsHvAp8Jq7v9H4T0MSUVW1c8frS7nm2bmM6NORl686mj7tdJ4dkUiFclTPLCDYp2fTgyzbN7QzITC9EhgeSUBJTiXle7kmby4zl27ix0f05eZTh9IqVUVfpDHom7sSd77YuoufPjGblVt2cctpQ7lwbD9dGUukEanwS1z5cMUWrnx6Du7w5KVjOHpQZqwjiSQcFX6JG09/soYbX17IgMx0Hrsoh35d0mMdSSQhqfBLXHhs1ipunbaY3CFduX/iSNq10UnWRJqKCr/E3D8+WMkfX1vCD4Z15/6JI3WFLJEmpr8wial9RX/CYSr6ItGiPX6JmUff/5w/TV/KyYf14N7zRqjoi0SJCr/ExN/e+5w7Xl/KyYf34L5zR5Cqoi8SNSr8EnUPv/s5f35jKacO78k95wxX0ReJMhV+iaqH8ldw55vL+OHwntytoi8SE/qrk6jZV/RPG6GiLxJL2uOXqHhgZhF3zVjOGSN78dezh5PSQqdgEIkVFX5pcve9XcQ9by/nzJG9uFNFXyTmVPilSd379nLufbuIH43qzV/OOlxFXyQOqPBLk3B37nm7iPtnFnFWdm/+/CMVfZF4ocIvjc7duWfGcu5/ZwXn5PTmjjMPp4WKvkjcUOGXRuXu3D1jOQ+8s4Jzc/pw+5mHqeiLxJlQrrnbx8zyzWyJmS0ys18Flnc2sxlmVhS471RP/4sCbYoCF2eXBOXu/PWtZTzwzgrOG62iLxKvQjmQuhK41t0PAcYCV5nZocD1wEx3HwzMDMx/g5l1Bm4CjgDGADfV9wYhzd9dby3nofzPmTimL386Q0VfJF41WPjdfYO7zwlMlwBLgF7AacATgWZPAKcH6X4SMMPdi919GzADGN8YwSW+PPPpGh7Mr9nTv+30YSr6InHM3D30xmb9gfeBYcAad+9Ya902d+9Up/11QBt3/2Ng/kZgt7v/NcjPngRMAsjKysrOy8s74CcDUFpaSkZGRlh9oyER8y3eWsVds8s5tEsK14xq3aRH7yTi9osm5YtMPOfLzc0tdPeckBq7e0g3IAMoBM4MzG+vs35bkD6/Af5Qa/5GaoaN9vtY2dnZHq78/Pyw+0ZDouUr2ljiw256w8fd/Z7v3F3RNKFqSbTtF23KF5l4zgfM9hDreUgnSzGzlsALwFR3fzGweKOZ9Qis7wFsCtJ1LdCn1nxvYH1I70gS94p3VXDplAJap7bgsYtzdLlEkWYilKN6DHgMWOLud9da9Qqw7yidi4CXg3R/ExhnZp0CH+qOCyyTZm5PZRWX/3M2X+0s59Gf5NC7U1qsI4lIiELZ4z8auBD4vpnNDdwmAHcAJ5pZEXBiYB4zyzGzfwC4ezFwK1AQuN0SWCbNmLtz/QsLKFi9jbvOHs6ovjpQS6Q5afALXO4+C6jv07rjg7SfDfy01vxkYHK4ASX+PPjOCl76bB3Xnngwpw7vGes4InKAdEJ0OSCvzlvPXTNqzrR59fcHxTqOiIRBhV9CNmfNNq59fh6j+3fi9h8dRs3HPyLS3KjwS0i+LC5j0pOz6d6+DY9cmEPr1JRYRxKRMOkkbdKgkvK9/PSJ2eyprCZv0mg6p7eKdSQRiYAKv+xXZVU1Vz/9GSs2l/LEJWMY1C0+v7UoIqHTUI/s163TFvPe8s388fRhfHdwZqzjiEgjUOGXek35zyqe+OgLfnbMACaO6RvrOCLSSFT4Jaj8pZu4ZdpiTjgki+t/cEis44hII1Lhl29Z+tVOfvHMZ3yne3vuO2+ErpUrkmBU+OUbtu+p5rIps0lvncJjF+eQ3lqf/4skGv1Vy9fK91Zx/5w9FJcZz11+JD06tI11JBFpAtrjFwCqq51rn5/Hqh3V3HveCA7r3SHWkUSkiajwCwB3z1jOa/M3cPaQlpw0tHus44hIE9JQj/CvwrVfXy/3pM5bYx1HRJqY9viT3Ccrt3LDi/M56qAu3Hr6MJ14TSQJqPAnsVVbdnH5U4X06ZzGw+dn0zJFLweRZKC/9CS1vayCy6YUYMDjF4+mQ5qulyuSLBoc4zezycApwCZ3HxZY9iwwJNCkI7Dd3UcE6bsaKAGqgEp3z2mk3BKBispqrniqkLXbdjP1Z0fQr0t6rCOJSBSF8uHuFOBB4Ml9C9z93H3TZnYXsGM//XPdfUu4AaVxuTu/f2kBH68s5p5zhzO6f+dYRxKRKAvlmrvvm1n/YOus5pPAc4DvN24saSp/e28lzxeu5ZfHD+aMkb1jHUdEYsDcveFGNYV/2r6hnlrLjwXurm8Ix8xWAdsABx5x90f38xiTgEkAWVlZ2Xl5eSE+hW8qLS0lIyN+zxkfy3wFX1Xy0Nw9HNE9hSuGtw56BI+2X2SULzLKF77c3NzCkIfT3b3BG9AfWBhk+cPAtfvp1zNw3w2YBxwbyuNlZ2d7uPLz88PuGw2xyjd3zTYf8ofpfsZDs3x3RWW97bT9IqN8kVG+8AGzPYT66u7hH9VjZqnAmcCz+3lTWR+43wS8BIwJ9/EkfOu27+anT84mM6M1j/4khzYtdb1ckWQWyeGcJwBL3X1tsJVmlm5m7fZNA+OAhRE8noShdE8ll00poLyiiscvHk1mRutYRxKRGGuw8JvZM8BHwBAzW2tmlwVWnQc8U6dtTzObHpjNAmaZ2TzgU+A1d3+j8aJLQyqrqvnF03Mo2lTKQ+ePYnBWu1hHEpE4EMpRPRPrWX5xkGXrgQmB6ZXA8AjzSQT++NoS8pdt5rYzhnHswV1jHUdE4oS+uZugnvxoNVM+XM1l3x3A+Uf0i3UcEYkjKvwJ6N1lm7j5lUWccEg3fjdB18sVkW9S4U8wy74q4eqn910vd6Sulysi36LCn0A2l+zh0ikFul6uiOyXKkOCqKis5vJ/zqZ4V4Wulysi+6XCnyDueH0pc9Zs5//OH6Xr5YrIfmmoJwG8sXADk/+zikuO7s+Ew3rEOo6IxDkV/mbui627+M3z8xnepyM3/EBH8IhIw1T4m7HyvVVcOXUOLVoYD/14JK1S9esUkYZpjL8Zu3XaYhat38ljF+XQu1NarOOISDOhXcRm6uW565j6yRqu+N5BHH9IVqzjiEgzosLfDK3YVMoNLy5gTP/OXDfu4FjHEZFmRoW/mdldUcWVUwtp2zKF+yeOJDVFv0IROTAa429mbnx5IUWbSnny0jF079Am1nFEpBnS7mIz8tzsL/lX4Vp+8f3BHDNYp1kWkfCo8DcTSzbs5MZ/L+Sog7rwq+MHxzqOiDRjKvzNQOmeSq6aOof2bVvqjJsiErFQLr042cw2mdnCWstuNrN1ZjY3cJtQT9/xZrbMzFaY2fWNGTxZuDs3vLiA1Vt38cDEkXRtp2vmikhkQtnjnwKMD7L8HncfEbhNr7vSzFKAh4AfAIcCE83s0EjCJqOnPlnDq/PWc+24IYwd2CXWcUQkATRY+N39faA4jJ89Bljh7ivdvQLIA04L4+ckrQVrd3Drq4s5bkhXfv69g2IdR0QSRCRj/Feb2fzAUFCnIOt7AV/Wml8bWCYh2LF7L1c+XUhmRivuOWcELTSuLyKNxNy94UZm/YFp7j4sMJ8FbAEcuBXo4e6X1ulzNnCSu/80MH8hMMbdf1HPY0wCJgFkZWVl5+XlhfWESktLycjICKtvNISSz915cO4e5m6q4oYxbRjUKSVK6RJj+8WS8kVG+cKXm5tb6O45ITV29wZvQH9g4YGsA44E3qw1fwNwQyiPl52d7eHKz88Pu280hJLvHx+s9H6/neZ/f//zpg9URyJsv1hSvsgoX/iA2R5CfXX38IZ6zKz21T7OABYGaVYADDazAWbWCjgPeCWcx0smc9Zs4/bpSzjx0Cwu++6AWMcRkQTU4CkbzOwZ4Dgg08zWAjcBx5nZCGqGelYDlwfa9gT+4e4T3L3SzK4G3gRSgMnuvqhJnkWC2LargqunzqFHxzb89azhmGlcX0QaX4OF390nBln8WD1t1wMTas1PB751qKd8W3W18+vn5rKltIJ//fxIOqS1jHUkEUlQ+uZunLj/nSLyl23mxlMO4fDeHWMdR0QSmAp/HHhn6UbufbuIM0f14oKx/WIdR0QSnAp/jH2xdRfX5M3l0B7t+dMZh2lcX0SanAp/DJVVVHL5PwsxMx65MJs2LaN3vL6IJC9diCVGPHDytWUbS5hyyRj6dNbF0kUkOrTHHyNTPlzNy3PX8+sTDuZ7B+uiKiISPSr8MVCwupjbXlvCCYdkcVXuoFjHEZEko8IfZdvKq7ly6hz6dE7j7nOH6+RrIhJ1GuOPoorKah6au4fScuOpy46gfRt9SUtEok+FP4pue20xK7ZX88DEkQzp3i7WcUQkSWmoJ0pe+mwtT3z0BSf1T+XU4T1jHUdEkpgKfxQsWr+DG15cwBEDOnPOwa1iHUdEkpwKfxPbXlbBFU8V0rFtKx788ShS9GGuiMSYxvibUHW1c82zc/lqRznPXn4kXdu1jnUkERHt8Tele2cW8e6yzdx06lBG9Q12WWIRkehT4W8iby/eyP0zizgruzfnH9E31nFERL6mwt8EVm/Zxf88N5dhvdrzx9OH6YybIhJXGiz8ZjbZzDaZ2cJay+40s6VmNt/MXjKzoFcOMbPVZrbAzOaa2ezGDB6v9p1xM6WF8fD5OuOmiMSfUPb4pwDj6yybAQxz98OB5cAN++mf6+4j3D0nvIjNh7tz/QsLWL6phPvPG6kzbopIXGqw8Lv7+0BxnWVvuXtlYPZjoHcTZGt2Jv9nNa/MW89144ZwrM64KSJxqjHG+C8FXq9nnQNvmVmhmU1qhMeKW/lLN3Hba4sZd2gWP//eQbGOIyJSL3P3hhuZ9QemufuwOst/D+QAZ3qQH2RmPd19vZl1o2Z46BeB/yCCPcYkYBJAVlZWdl5e3gE+lRqlpaVkZGSE1Tdcq3dUcfun5XRPb8ENY9rQJrX+D3Njke9AKF9klC8yyhe+3NzcwpCH1N29wRvQH1hYZ9lFwEdAWog/42bgulDaZmdne7jy8/PD7huOtdvKPOePM/yo22f6xh27G2wf7XwHSvkio3yRUb7wAbM9hPrq7uEN9ZjZeOC3wA/dvayeNulm1m7fNDAOWBisbXO1Y/deLnn8U8r3VvH4JaPp1r5NrCOJiDQolMM5n6Fmz36Ima01s8uAB4F2wIzAoZp/C7TtaWbTA12zgFlmNg/4FHjN3d9okmcRAxWV1fz8qUJWbdnFIxdkc3CWTrMsIs1Dg+fqcfeJQRY/Vk/b9cCEwPRKYHhE6eKUu3P9i/P58POt3HX2cI4alBnrSCIiIdM3d8Nwz9tFvDhnHf9zwsH8KFtHsopI86LCf4Cem/0l988s4uzs3vzyeF0oXUSaHxX+AzCraAu/e3EBxwzO5E9nHqZz8IhIs6TCH6KlX+3k508VMqhbBv93/ihapmjTiUjzpOoVgq92lHPJ4wWktU7h8UtG065Ny1hHEhEJmwp/A0r3VHLJlAJ27t7L5ItH06ND21hHEhGJiC69uB97q6q5auoclm8sYfLFoxnas0OsI4mIREx7/PVwd/735YW8t3wzt50+jO/pbJsikiBU+Ovxf+9+zjOffslVuQdx3hhdOlFEEocKfxAvz13HnW8u47QRPblu3JBYxxERaVQq/HV8snIrv3l+PmMHduYvZx2uY/VFJOGo8NeyYlMpk/5ZSN8uaTxyQQ6tU3W9XBFJPCr8ATvL93LplAJaprTg8YtH0yFNx+qLSGLS4ZzUHMHz+5cWsm77bp67/EhdJF1EEpr2+IF/Fa7l1Xnr+fWJB5Pdr1Os44iINKmkL/wrN5dy0yuLGDuwM1foIukikgSSuvBXVFbzy7zPaJXagnvPHUlKCx3BIyKJL6TCb2aTzWyTmS2stayzmc0ws6LAfdAxEjO7KNCmyMwuaqzgjeHON5eycN1O/vKjw+neQdfLFZHkEOoe/xRgfJ1l1wMz3X0wMDMw/w1m1hm4CTgCGAPcVN8bRLS9t3wzf/9gFReM7cu4od1jHUdEJGpCKvzu/j5QXGfxacD/Q5zqAAAIAElEQVQTgekngNODdD0JmOHuxe6+DZjBt99Aom5L6R6ufW4eB2dl8IeTD411HBGRqIpkjD/L3TcABO67BWnTC/iy1vzawLKYqa52rnt+HjvL93L/xJG0aakvaYlIcjF3D62hWX9gmrsPC8xvd/eOtdZvc/dOdfr8Bmjt7n8MzN8IlLn7XUF+/iRgEkBWVlZ2Xl5eWE+otLSUjIyMete/uXovzyyt4IJDWnFCv+h/SauhfLGmfJFRvsgoX/hyc3ML3T0npMbuHtIN6A8srDW/DOgRmO4BLAvSZyLwSK35R4CJDT1Wdna2hys/P7/edQvWbvfBv5vul00p8Orq6rAfIxL7yxcPlC8yyhcZ5QsfMNtDrOeRDPW8Auw7Suci4OUgbd4ExplZp8CHuuMCy6KurKKSX+Z9Rqf0ljr5mogktVAP53wG+AgYYmZrzewy4A7gRDMrAk4MzGNmOWb2DwB3LwZuBQoCt1sCy6LullcXs2rLLu45ZwSd01vFIoKISFwI6Vw97j6xnlXHB2k7G/hprfnJwOSw0jWS6Qs2kFfwJT8/7iCOGpQZyygiIjGX8N/cXbd9N9e/MJ/hfTry6xMPjnUcEZGYS+jCX1lVzTV5n1HtcP95I2iZktBPV0QkJAl9WuYH81dQsHob9547gn5d0mMdR0QkLiTsLnDB6mLun1nEmSN7cfrImH5nTEQkriRk4d9Rtpdr8ubSp3Mat5w+LNZxRETiSsIN9bg7v3tpARt3lvPCz48io3XCPUURkYgk3B7/++sqeW3BBq4dN4ThfTo23EFEJMkkVOFfsamUqUsqOHpQFy4/dmCs44iIxKWEKfx7Kqv45TOf0aoF3H3OCFroaloiIkElzAB4VbXzne7tOLHHHrLa62paIiL1SZg9/rRWqdx97ghGdkuY9zIRkSaRMIVfRERCo8IvIpJkVPhFRJKMCr+ISJJR4RcRSTIq/CIiSUaFX0Qkyajwi4gkGXP3WGf4FjPbDHwRZvdMYEsjxmlsyhcZ5YuM8kUmnvP1c/euoTSMy8IfCTOb7e45sc5RH+WLjPJFRvkiE+/5QqWhHhGRJKPCLyKSZBKx8D8a6wANUL7IKF9klC8y8Z4vJAk3xi8iIvuXiHv8IiKyH8228JvZeDNbZmYrzOz6IOtbm9mzgfWfmFn/KGbrY2b5ZrbEzBaZ2a+CtDnOzHaY2dzA7X+jlS/w+KvNbEHgsWcHWW9mdn9g+803s1FRzDak1naZa2Y7zeyaOm2iuv3MbLKZbTKzhbWWdTazGWZWFLjvVE/fiwJtiszsoijmu9PMlgZ+fy+ZWdCLUDf0WmjCfDeb2bpav8MJ9fTd7996E+Z7tla21WY2t56+Tb79Gp27N7sbkAJ8DgwEWgHzgEPrtLkS+Ftg+jzg2Sjm6wGMCky3A5YHyXccMC2G23A1kLmf9ROA1wEDxgKfxPB3/RU1xyjHbPsBxwKjgIW1lv0FuD4wfT3w5yD9OgMrA/edAtOdopRvHJAamP5zsHyhvBaaMN/NwHUh/P73+7feVPnqrL8L+N9Ybb/GvjXXPf4xwAp3X+nuFUAecFqdNqcBTwSm/wUcb2ZRuRCvu29w9zmB6RJgCdArGo/diE4DnvQaHwMdzaxHDHIcD3zu7uF+oa9RuPv7QHGdxbVfY08ApwfpehIww92L3X0bMAMYH4187v6Wu1cGZj8Gejf244aqnu0XilD+1iO2v3yBunEO8ExjP26sNNfC3wv4stb8Wr5dWL9uE3jx7wC6RCVdLYEhppHAJ0FWH2lm88zsdTMbGtVg4MBbZlZoZpOCrA9lG0fDedT/BxfL7QeQ5e4boObNHugWpE28bMdLqfkPLpiGXgtN6erAUNTkeobK4mH7HQNsdPeietbHcvuFpbkW/mB77nUPTwqlTZMyswzgBeAad99ZZ/UcaoYvhgMPAP+OZjbgaHcfBfwAuMrMjq2zPh62Xyvgh8DzQVbHevuFKh624++BSmBqPU0aei00lYeBg4ARwAZqhlPqivn2Ayay/739WG2/sDXXwr8W6FNrvjewvr42ZpYKdCC8fzXDYmYtqSn6U939xbrr3X2nu5cGpqcDLc0sM1r53H194H4T8BI1/1LXFso2bmo/AOa4+8a6K2K9/QI27hv+CtxvCtImptsx8GHyKcD5HhiQriuE10KTcPeN7l7l7tXA3+t53Fhvv1TgTODZ+trEavtForkW/gJgsJkNCOwVnge8UqfNK8C+IyjOAt6p74Xf2AJjgo8BS9z97nradN/3mYOZjaHmd7E1SvnSzazdvmlqPgRcWKfZK8BPAkf3jAV27BvWiKJ697Riuf1qqf0auwh4OUibN4FxZtYpMJQxLrCsyZnZeOC3wA/dvayeNqG8FpoqX+3PjM6o53FD+VtvSicAS919bbCVsdx+EYn1p8vh3qg56mQ5NZ/4/z6w7BZqXuQAbagZIlgBfAoMjGK271Lz7+h8YG7gNgG4Argi0OZqYBE1Ryl8DBwVxXwDA487L5Bh3/arnc+AhwLbdwGQE+Xfbxo1hbxDrWUx237UvAFtAPZSsxd6GTWfGc0EigL3nQNtc4B/1Op7aeB1uAK4JIr5VlAzPr7vNbjvKLeewPT9vRailO+fgdfWfGqKeY+6+QLz3/pbj0a+wPIp+15ztdpGffs19k3f3BURSTLNdahHRETCpMIvIpJkVPhFRJKMCr+ISJJR4RcRSTIq/CIiSUaFX0Qkyajwi4gkmf8P4tKGWPkG1gYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "theta2, likelihood_log2, fig2, ax2 = EM(theta.copy(), data.copy(), gt2.copy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model Comparison\n",
    "We plot both resulting models here for a test ground truth to do a comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VOW9x/HPTwIIQYSIICAVa61KtXVBtLX2UlFE2wtuVVwqbnVfsLVq9Xq1WFvcq9cFES1qFa0KStUqiKJSFUUEEQFBRJuAgCBrZAn53T+eExjCJJmQSc7JzPf9ep3XzJzznDk/hkl+eZ7zLObuiIiIJM02cQcgIiKSjhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgkUkHcAaSzzTbbeIsWLeIOQ0Qkr5SWlrq7J6biksgE1aJFC1avXh13GCIiecXMvo07hlSJyZQiIiKplKBERCSRlKBERCSRlKBERCSRlKBERCSRlKBERCSRlKBERCSREjkOSkREEsisC/AosBNQDgzF/S7MngL2iEq1AZbhvm+a8+cBK4ENQBnu3au7nBKUSANzh7KysG3YsHXPN2yA8vItH9Ptq6qMe3he8Zj6vKbHiufZ2io+l619XvE69THdvprKV1WmNuWqKrM1ZQG2KS+j5frlFK5fRst1y8JjpdcVW7MN3/Lj4qfTv1F2lAG/w30yZtsBH2A2FveTNpYwux1YXs17/Bz3rzO5mBKU5IXycigthW+/rX5bs6b642vXwvr1sG7dpi31dVXPU1+Xl8f9aWSfWd231PfZmucVr1Mf0+2rqXxVZWpTrqoyqfu2L1tC13Wfssu62eyybjad18+jdfk3tN6wjO3Kl7HdhmVsV76cwvJV6d84Uo6xapvWrGzShpVN2oS/Ppo0qfacrea+AFgQPV+J2QygM/BJ9A8z4ETgsGxcLpEJqqioiPHjx8cdhiREefmmWkPlraI2kMm+2ioogO22g9atYZttwma26THdVtWxdPsh/fN0r6sqX5fHTJ/XdEyqVrBqFS2Ki8NWUkLLiufFJTRduXJjOd9mG9buuCPrW7emrFUryloVUdbqOywtLGRRq1bRvmgrLNzs9YaWLcMXLPLVW2/VKWQzm5Tyeqi7D01b0qwrsB8wMWXvocBC3GdX8f4OjMHMgQeo6r0rLuFV1StjVFhY6JqLL3etXQuLF8OiRZlta9dW/37bbBOSSMVWkVQqb61aQYsWsO224bGqLfX4tttu9rMvsrny8vBlLimBOXNg9uxN26efwtcpLVlm8J3vwO67b7ntuis0axbfv2NjiFbq7oUZFGwFvAHchPvIlP33A3Nwv72K8zrhPh+z9sBY4BLc36zyMkpQkk3r18OXX8Jnn23aPv8cFi7clHCWV9E63bw5dOgA7dtvubVrB9tvnz4JtWypv+olyzZsCIln/nxYsKDqx6++CmVTde6cPgnttlv4iyfBMkpQZk2BF4BXcL8jZX8BUAIcgHtxBhe7AViF+21VFUlkE58k26pVMHfu5kmoYvvii81/XrfdFrp2hY4dYf/9q05A7duHGo4SjTSYJUvggw/C9sUXmyefhQu3TDwAO+4YvsydOsE++4THite77Ra2wporII1WuMf0EDBjs+QUHA7MrDI5mRUC20T3rgqB3sCgai+nGpSk4x5qQu+8A7NmbZ6EFi7cvGxR0aafzcpbx45qIpMEWLECJk+G99+HSZPCNnfupuM77rh5skn3uNNOiWiGq0811qDMfgq8BUwjdDMHuAb3lzAbDryL+5CU8p2AYbgfjdl3gVHRkQLgCdxvqjYeJSiB8MfitGkwYQL8+9/hsTj6O8gMdt656iTUpk28sYtsprQUpkwJSagiIc2atakfd9eu0L07HHhgeDzggNB+LJnfg2ogGTXxmXE5cA6hB8Y04Ex31qQcb04YvHUAsAQ4yZ15Kce/Q+iGeIM7VbY3SsMpLYWJEzclo7ffhopORZ07w09/GrZDDoG99kp807nkq7Vr4eOPNyWi99+H6dM3Nc916hSS0KmnbkpGO+4Yb8ySsRoTlBmdgUuBbu58a8Y/gP7A8JRiZwPfuPM9M/oDNwMnpRy/E/hX1qKWWlu0aFMymjAhtHaUlYXa0Q9+EH5+K5LSd76je0GSMCtXwsyZMGPGpu2TT0IzXUUy2mGHUCvq23dT7ahTp3jjljrJtJNEAdDCjPVAS2B+peP9gBui588A95hh7rgZxwBzAbXZNaDSUnjuOXj11ZCQZkejEpo3hx494Pe/D8noxz+Gtm3jjVUECE1wixdvnoQqtuKU++5Nm4ZecT/8IZx0Ung88EDYZRf9ZZVjakxQ7pSYcRvwJfAtMMadMZWKdQb+E5UvM2M5sIMZ3wJXAUcAV1R3HTM7FzgXoFmO34isT1OmwLBh8Pe/h+7cO+wQmul+85uQkPbfPyQpkVht2BDamN99d/NEtHTppjKFhaF9+ec/D48V23e/G5KU5LxMmvjaEmpIuwLLgKfNOM2dv6cWS3OqA38E7nRnVU1/2ESjlYdC6CSRUfQChA5KTz4JDz4YmuGbN4fjj4dzzoH/+i/1opOEWLoUXnkFXnwRXn45dPOGMMhtr73ghBPCY7du4XHnnVUjynOZNPEdDnzuzmIAM0YCP4HNElQx0AUoNqMA2B5YChwEnGDGLYQZbsvNWOPOPVn8N+Ql9/DH57BhITmVlsLee8Ndd8Fpp4Wu3yKxcg9dQ198MWzvvBNmXmjXDo4+Gn7xi1A7at8+7kgloTJJUF8CB5vRktDE1wuYVKnMaGAA8A5wAvCaO06YlwkAM24AVik51c2SJfDYYyExTZ8eWkFOPjk04fXooT84JWarV8O4cSEhvfTSpntH++8P114bEtOBB9bfZKaSUzK5BzXRjGeAyYSp1j8EhpoxCJjkzmjCyOLHzJhDqDn1r8eY8055OYwfH5LSs8+GGbF79IChQ6F//zDtj0hsPvssJKMXXwxf1LVrw7QgvXvDDTfAUUepN51sFQ3UTbAFC2D4cHjoofA7oE2b0Hz3m9+EjksisVi9OnQNHTMmJKVZs8L+738/NNv94hdw6KE5P+tCLkraQF0lqARatw4uuyx0etiwIXR0OOec0PGhRYu4o5O8s25duOH52mthe/fdMCtws2bQs2dISEcfDd/7XtyRSh0lLUFpstiEWboUjjsO3ngDLroILr00/GEq0mA2bAgjuSsS0ltvhdUat9km3Ev67W/hsMPC+IVcnhhVYqcElSBz5oQ/RufNC+OYTj017ogkL7iHHjcVCWn8+E1rouy9d2hTPuww+NnPNKpbGpQSVEK89RYcc0zohffqq6EJX6TezJ0bettVJKVFi8L+3XaDE08MCennPw/ro4jERAkqAR57DM4+Oyyq+eKLasqXerJuHYwcCffcEyZmhLCMRO/eISEddliYLkgkIZSgYuQO118PN94Y/lh99lm1oEg9mD8fHnggjEv46qswVdAtt8B//zfssYcGz0liKUHFZM0aOOssGDEiPN5/v3rlSha5h3bje+6BUaNCx4ejjoKLL4Yjj9T8V9IoKEHFYPFi6NcvzPwyeDBceaX+iJUsWbUKHn8c7r03TDPUtm0Ys3DBBeH+kkgjogTVwGbMCD31FiyAp58O82OK1Nmnn8J994WR3cuXw777hqlHTj4ZWraMOzqRraIE1YBefTUkpG23DeOcevSIOyJp1DZsCFMM3XtvmCW8adPwBbv44rDQl6rl0sgpQTWQBx8MrSx77QUvvKDOUlIHS5aE+a/uvz8MmuvcGQYNCuOVdtop7uhEskZ3SutZeXm4x3TuuXDEEaF3r5KT1FrF+ioDBoR1kq66Crp2De3En38O112n5CT1z6wLZq9jNgOz6ZhdFu2/AbMSzKZE29FVnN8Hs1mYzcHs6hovp7n46s/q1fDrX4dOVBdeGNZqKlCdVWpj5Up44olQW5o6NUxdf9pp4Qu1995xRyc5psa5+Mw6Ah1xn4zZdsAHwDHAicAq3G+r5twmwKeEFdaLgfeBk3H/pKpT9OuynsyfD337hinN/vrXMKeebglIxj76CIYMCXNerVwZOj088EDo9KD1VSQu7guABdHzlZjNADpneHYPYA7ucwEwe5KwWrsSVEOaMSMMzv/mG3j++TAeUqRGa9aEJrshQ+Dtt0NvmpNOCjcvtRqlJI1ZV2A/YCJwCHAxZqcTFrT9He7fVDqjM/CflNfFhFXXq5TIBFVUVMT48ePjDmOrffopXHFFmIW8RYsw96ZIVVoUF9Ppn/9kp5dfpumKFZR26cL8Cy/kqyOPpKx16zCT+BtvxB2m5IcCM0tdMX2ouw/dopRZK+BZYCDuKzC7H7gR8OjxduCsymeluV6195gSmaCWLl1Kz5494w5jq0yYAOedB3feGQbui6RVVgajR4d7S6++Gm5OHnMMnH8+LQ87jO+ZoSkZJQZl7t692hJmTQnJ6XHcRwLgvjDl+IPAC2nOLAa6pLzeGZhf7aXUSSK7+vQJ953mzdP4SEmjpCTMiTdsWLhR2aVL6OJ59tlh4laRGGXQScKAR4CluA9M2d8xuj8FZpcDB+Hev9K5BYROEr2AEkIniVNwn17V5RJZg2qs3nsvjJccPFjJSSopLQ1fjFtuCbOK9+kTak9HH62undKYHAL8GpiG2ZRo3zXAyZjtS2iymwecB4BZJ2AY7kfjXobZxcArQBPg4eqSE6gGlVX9+oX5Ob/4Qh2tJOIexhlcfjl8+WXohfenP4UZxUUSJmlLvmugbpZMnRpuKQwcqOQkkZkzw8zhxx8P228fess88YSSk0iGlKCy5KaboHXrMN5J8tzKlWH6kH32Ce2+d98dbkz+13/FHZlIo6LG7yyYMQOeeQb+8Ado0ybuaCQ27qGG9Pvfh+nqzzoL/vIXaN8+7shEGqWMalBmXG7GdDM+NmOEGdtWOt7cjKfMmGPGRDO6RvuPMOMDM6ZFj4fVw78hdn/+cxjvdPnlcUcisZk6NdSQTjstTN767rthQlclJ5GtVmOCMqMzcCnQ3Z29Cb0v+lcqdjbwjTvfA+4Ebo72fw38tzv7AAOAx7IVeFLMmRP+aL7gAmjXLu5opMF98w1ccgnsvz988kmYtn7iRDio2gHyIpKBTO9BFQAtzCgAWrLl4Kp+hL7xAM8Avcwwdz5031h2OrCtGc3rGnSSDB4cluG54oq4I5EGVV4eakjf/35YKPCCC8IUIueco+XURbKkxp8kd0qA24AvCZMELndnTKViG+dYcqcMWA7sUKnM8cCH7qyta9BJ8cUX8MgjWoYn77z3Hhx8cEhGe+wBH3wA99wDRUVxRyaSUzJp4mtLqCHtCnQCCs04rXKxNKduHGBlxg8IzX7nVX0dO9fMJpnZpLKyskxij90tt4T5O6+8Mu5IpEEsXhyS0kEHwX/+A489Fga+7btv3JGJ5KRM2iIOBz53Z7E764GRwE8qldk4x1LUDLg9sDR6vTMwCjjdnc+quoi7D3X37u7evaARjKyfPz+08JxxRpitRnKYOwwfDnvuGarMV1wBs2aFDhGaYVyk3mSSoL4EDjajpRlGmEdpRqUyowmdIABOAF5zx81oA7wI/MGdf2cr6CS47bYw3+fVNa4JKY3a7NnQqxeceSZ06xZ66916axj0JiL1KpN7UBMJHR8mA9Oic4aaMciMvlGxh4AdzJgD/Bao+LV9MfA94DozpkRbo+93u2hRWLLn1FM1KUDOWr8+jB/YZ58wyPaBB8KSF926xR2ZSN7QXHxb4eqrw/2nGTPCPXLJMRMnhp4v06bBCSeEmSA007jkAc3F18gtXQr33gsnnqjklHNWrgxzVf34x+E/+vnnwwq3Sk4isVCCqqW77oJVq+Daa+OORLLqn/8MzXf33AMXXxwG3fbtW/N5IlJvlKBqYfny0NpzzDHh1oTkgAUL4Fe/CsmoTRt4++3wn6xOECKxU4KqhXvvhWXL4H/+J+5IpM7Ky8PKtnvtFWpPf/5z6Axx8MFxRyYiEXWSyNDq1bDLLtCjB7z0UtzRSJ3MmBGWWZ8wAQ47LHTJ3H33uKMSiZ06STRSQ4bAkiVw3XVxRyJbbe1a+OMfw8wP06fD3/4Gr76q5CSSUKpBZeDbb8N4p27dYNy4uKORrTJpEpx+eqg9nXIK3HmnlsIQqSRpNajkzymUAA89BF99BSNGxB2JbJXRo6F/f9hxx9A+e9RRcUckIhlQDaoG69bBbruF+09vvaWp1xqdIUPgoovggAPghRdUaxKpRtJqULoHVYNHHoHi4nDvScmpEXEP3S0vuCDUmF5/XclJpK7MumD2OmYzMJuO2WXR/lsxm4nZR5iNwqxNFefPw2waZlMwm1Tj5VSDqlpZWViPrl27MPuNElQjsX59mKqoYrGu++6DRjBDvkjcaqxBmXUEOuI+GbPtgA+AY4CdgddwL8MsrKjuflWa8+cB3XH/OpN4VIOqxhNPwOefhz/ElZwaiZUr4Re/CMlp0KAwyauSk0h2uC/AfXL0fCVhZYvOuI/BvWIhv3cJCavOEvmTW1RUxPjx4+MOg5Urw+Dc1q0hAeFIDZotWcI+V19Nq7lzmXXllXx16KFhBnIRyVSBbd70NtTdh6YtadYV2A+YWOnIWcBTVby/A2Mwc+ABqnrvikuoiS+9p54KHb/+8Y8wE44k3IwZ4V7T11/DM89Anz5xRyTS6GTcScKsFfAGcBPuI1P2Xwt0B44jXXIx64T7fMzaA2OBS3B/s8rLKEFtqbwcfvQj2LABPv4YtlFDaLJNmBDm0mvWDF58MfTYE5FayyhBmTUFXgBewf2OlP0DgPOBXriXZnCxG4BVuN9WVRH96k1j9OiQmK65Rskp8Z59Fg4/PIxxeucdJSeR+mRmhAVqZ1RKTn2Aq4C+VSYns8KoY0V4Dr2Bj6u9nGpQm3OHAw8Mk8LOnKn764l2990wcGCY4HX06NDdUkS2Wga9+H4KvEVYXb082nsNcDfQHFgS7XsX9/Mx6wQMw/1ozL4LjIqOFwBP4H5TtfEoQW3u9dfD/KHDhsHZZ8cSgtSkvDwsa3zrrWHtkyeegBYt4o5KpNFL2kBd1Q8qmTIlPB57bLxxSBXWroUzzwzzTl14YahFNWkSd1QiUg+UoCopKQl/jLdtG3cksoVly+C440I1d/BguPJKDVATyWFKUJUUF0Pnzvq9lzjFxaEb+axZ8NhjcNppcUckIvVMCaqSkpKQoCRBZs6EI46A5cvhX/+CXr3ijkhEGoA6UVdSUgI7Z2WSDsmKefNCN/L168N08kpOInlDNagU7qpBJcpXX4XkVFoapizaZ5+4IxKRBpRRDcqMy82YbsbHZowwY9tKx5ub8ZQZc8yYaEbXlGN/iPbPMuPILMefVV9/HdZ/UoJKgG++gd69Q5J66SUlJ5E8VGOCMqMzcCnQ3Z29gSZA/0rFzga+ced7wJ3AzdG53aKyPwD6APeZkdg+wSUl4VEJKmarVoUZyWfNgueeCwNxRSTvZHoPqgBoYUYB0BKYX+l4P+CR6PkzQC8zLNr/pDtr3fkcmAP0qHvY9aMiQekeVIzWrg1dySdOhCefDE18IpKXakxQ7pQAtwFfAguA5e6MqVSsM/CfqHwZsBzYIXV/pDjal0jFxeFRNaiYlJXBKafA2LHw8MMaLS2S5zJp4mtLqAntCnQCCs2oPAgl3aghr2Z/muvYuWY2ycwmlZWVpStS70pKwuSwO+0Uy+XzW3l5WP125Ej4619hwIC4IxKRmGXSxHc48Lk7i91ZD4wEflKpTDHQBSBqBtweWJq6P7IzWzYPAuDuQ929u7t3L4hphtaSkpCcNEFsA3OH3/0Ohg+H66+Hyy6LOyIRSYBMEtSXwMFmtIzuK/UiLPObajRQ8SfvCcBr7ni0v3/Uy29XYHfgveyEnn0Vs0hIA7vxxlBruuyykKBERMhgHJQ7E814BpgMlAEfAkPNGARMcmc0YX2Qx8yYQ6g59Y/OnW7GP4BPonMvcmdD/fxT6q6kBHbfPe4o8szdd4ekdMYZcMcdmmNKRDbSchsp2raFU0+Fe+5p8Evnp0ceCYnp2GPhH/9Q26pIzJK23IamOoqsXh0my1YX8wby3HNhwa3DDw9LZyg5iUglSlARDdJtQOPGwUknhaWLR42C5s3jjkhEEkgJKqIE1UAmToR+/eD734cXX4RWreKOSEQSSgkqolkkGsC0aWFNp512gjFjoKgo7ohEJMGUoCKaRaKeffZZmPy1RYswU0THjnFHJCIJpwQVKSmB7beHwsT0X8khJSVhwcH160Ny2nXXuCMSka1h1gWz1zGbgdl0zC6L9hdhNhaz2dFj2yrOHxCVmY1ZjdPFKEFFtFBhPalYNmPx4rAabrducUckIluvDPgd7nsBBwMXYdYNuBoYh/vuwLjo9ebMioDrgYMIk4ZfX2UiiyhBRTSLRD0oLw8Dy2bPhtGjQ689EWm83BfgPjl6vpIwq1BnNl/R4hHgmDRnHwmMxX0p7t8AYwnLMFVJCSqilXTrwR//GGpNd98NP/953NGISDaZdQX2AyYCHXBfABA9tk9zRq1Xt0jk6MiioiLGjx/foNf87W/DffsGvmzO2uHtt9ln0CAW9OnDrD320Acr0jgUmNmklNdD3X3oFqXMWgHPAgNxX5HhFGUZr26xMZhM3rWhLV26lJ49ezbY9UpKwh/4Q4ZAA142d82eDTffDPvvT8eRI+nYokXcEYlIZsrcvXu1JcyaEpLT47iPjPYuxKwj7gsw6wgsSnNmMdAz5fXOwPjqLqUmPtTFPKtWrQor4jZtGtZ2UnISyR1mRpgcfAbud6QcSV3RYgDwfJqzXwF6Y9Y26hzRO9pXpUTWoBqaZpHIEnc45xz45BN4+WXYZZe4IxKR7DoE+DUwDbMp0b5rgMHAPzA7m7BE068AMOsOnI/7ObgvxexG4P3ovEG4L63uYkpQaBaJrLnzTnjqKRg8OIx7EpHc4j6B9PeSIKwVWLn8JOCclNcPAw9nejk18RGa+Jo1g3bt4o6kERs/Hq68MjTvXXll3NGISA5QgiLUoDp10lp5W624GE48Maz2+Le/6YMUkaxQgkKzSNTJ2rVw/PGwZk1YOqN167gjEpEcoQSFZpGok8sug/feC6vj7rln3NGISA7J+wTlrlkkttpDD8EDD8DVV4dl20VEsijvE9SyZfDtt0pQtTZpElx0UViy/U9/ijsaEclBeZ+g1MV8KyxeHHrr7bQTjBgBTZrEHZGI5KC8HwelWSRqqawMTj45JKl//1t980Wk3uR9gtIsErV07bUwbhwMHw777x93NCKSw9TEFyWoTp3ijaNReOYZuOUWuOACGFDjYpgiInVSY4IyYw8zpqRsK8wYWKlMWzNGmfGRGe+ZsXfKscvNmG7Gx2aMMGPb+viHbK3iYmjfPswkIdX45BM480w4+GD461/jjkZE8kCNCcqdWe7s686+wAFAKTCqUrFrgCnu/BA4HbgLwIzOwKVAd3f2BpoA/bMYf52pi3kGli8P3cgLC0MtStlcRBpAbZv4egGfufNFpf3dCOvQ485MoKsZHaJjBUALMwqAlsD8OsSbdZpFogbl5XDGGTB3Ljz9tLK5iDSY2iao/sCINPunAscBmNED2AXY2Z0S4DbC9OsLgOXujEn3xmZ2rplNMrNJZWVltQxr62kWiRrcfDM89xzcdhscemjc0YhIHsk4QZnRDOgLPJ3m8GCgrRlTgEuAD4EyM9oC/YBdgU5AoRmnpXt/dx/q7t3dvXtBQcN0LlyzBpYsUYKq0pgxodfeKafApZfGHY2I5JnaZIKjgMnuLKx8wJ0VwJkAZhjwebQdCXzuzuLo2EjgJ8Df6xh3VsyPGhuVoNL4/PMw3mmffWDoUM1QLiINrjZNfCeTvnkPM9pENSwIi1O9GSWtL4GDzWgZJa5ewIy6BJxNmkWiCqWlYaaI8vKwbHthYdwRiUgeyqgGZUZL4AjgvJR95wO4MwTYC3jUjA3AJ8DZ0bGJZjwDTAbKCE1/Q7P5D6gLzSKRhjucfz5MnQovvgi77RZ3RCKSpzJKUO6UAjtU2jck5fk7wO5VnHs9cH0dYqw3mkUijXvvhcceg0GD4Kij4o5GRPJYXs8kUVICrVppjb2NJkyAyy+Hvn1D5wgRkRjldYKq6GKu+/+EHiO/+hXsuis8+ihsk9dfDRFJgLyeLFazSETWrQvJaeVKePVV2H77uCMSkSQyexj4JbAI972jfU8Be0Ql2gDLcN83zbnzgJXABqAM9+41XS7vE1TPnnFHkQCXXw5vvw3/+Af84AdxRyMiyTUcuAd4dOMe95M2Pje7HVhezfk/x/3rTC+WtwmqvDy0auV9DWr4cLjvPvj970MtSkSkKu5vYtY17TEzA04EDsvW5RKZoIqKihg/fny9XqOsDAYPhu98B+r5UonV6tNP2f/ii1m+3358dOSReL5+ECJSocDMJqW8HurumQ4NOhRYiPvsKo47MAYzBx4gg/c1d8/w2g2nsLDQV69eXa/X+OAD6N4dRo2CY46p10sl09dfhw+gvDx8GDvuGHdEIhIzMyt19+pH5oca1Asb70Ft2n8/MAf326s4rxPu8zFrD4wFLsH9zeoulcgaVEOoGKSbl7NIbNgQpjH66qvQtVzJSUTqwqyAMGH4AVWWcZ8fPS7CbBTQA6g2QeVtX+K8HqT7P/8Teuvdd1+oRYmI1M3hwEzci9MeNSvEbLuNz6E38HFNb5rXCapJk7Cabl559tlw8+288+Css+KORkQaE7MRwDvAHpgVY3Z2dGTLpZjMOmH2UvSqAzABs6nAe8CLuL9c4+Xy9R7UGWfAa6/Bl1/W62WS5ZNP4KCDQlfyN96A5s3jjkhEEiSje1ANKG9rUHm3UGHFsu0tW4Zl25WcRCTh8raTRElJHo1JLS+HAQPgs89CtTEve4aISGOTtzWokpI8+j09eDA8/zzcfjv87GdxRyMikpG8TFArVoRp5/Kiie+VV0KvPS3bLiKNTF4mqLzpYq5l20WkEVOCylVr18Lxx4cVcrVsu4g0QnnZSSIvZpG47jr48MNw70nLtotII5TXNahOneKNo968/jrcdlsYjNu3b9zRiIjvkKzAAAAViUlEQVRslbxNUEVF0KJF3JHUg2++gdNPh+99L/TaExFppPKyiS+nu5hfdBEsWADvvKP7TiLSqOVlgsrZWSSeeAJGjIAbb4QDD4w7GhGROsnbJr6cS1BffAEXXgg/+QlcfXXc0YiI1FneJah162DRohxr4tuwIUxlVF4Of/87FORlxVhEckyNCcqMPcyYkrKtMGNgpTJtzRhlxkdmvGfG3inH2pjxjBkzzZhhxo/r4x+SqQULwtCgnKpB3XZbmJ38//4Pdt017mhERLKixj+13ZkF7AtgRhOgBBhVqdg1wBR3jjVjT+BeoFd07C7gZXdOMKMZ0DJbwW+NnBukO3lyGPN0wgmh956ISI6obRNfL+Azd76otL8bMA7AnZlAVzM6mNEa+BnwUHRsnTvL6hhzneRUgiothVNPDUu2DxmiqYxEJKfUNkFtuWpiMJWwHj1m9AB2AXYGvgssBv5mxodmDDMjbd9nMzvXzCaZ2aSysrJahpW5nJpF4qqrYOZMGD4cdtgh7mhERLIq4wQVNc/1BZ5Oc3gw0NaMKcAlwIdAGaEJcX/gfnf2A1YDabuYuftQd+/u7t0L6vEmf0kJbLsttG1bb5doGP/6F9xzDwwcCEccEXc0IiJZV5tMcBQw2Z2FlQ+4swI4E8AMAz6PtpZAsTsTo6LPUEWCaigVXcwbdWvY4sVw1llhxcW//CXuaERE6kVtmvhOJn3zXkVPvWbRy3OAN91Z4c5XwH/M2CM61gv4ZKujzYJGP4uEO5x7LixdCo8/HqqDIiI5KKMEZUZL4AhgZMq+8804P3q5FzDdjJmEmtZlKadfAjxuxkeE3oB/zkbgW6vRzyLx8MPw3HPw5z/Dj34UdzQikk/MHsZsEWYfp+y7AbMSzKZE29FVnNsHs1mYzcEso5Y0c/fsBJ5FhYWFvnr16qy/r3uYIPbSS+GWW7L+9vVvzhzYd1846CAYOxa2ybtx1iJSj8ys1N2rnsTT7GfAKuBR3PeO9t0ArML9tmrOawJ8SqjoFAPvAyfjXm2LWl79hluyJKzj1yib+MrK4LTToGlTeOQRJScRaXjubwJLt+LMHsAc3Ofivg54EuhX00l59Vuuoot5o2zi+9OfYOJEeOCBRpphRSSHXYzZR1ETYLo+0p2B/6S8Lo72VSuRk7YVFRUxfvz4rL/v8uVhVqBtt4V6ePt603r6dPa78UYWHnEEM9u3b1zBi0hjUmBmk1JeD3X3oTWccz9wI+DR4+3AWZXKpOs3XeP9pby6BzV0aFhk9ssvoUuXrL99/Vi5EvbbLzTxTZ0K228fd0QikqNqvAcVCnUFXth4DyqTY2Y/Bm7A/cjo9R8AcK92nEwia1D1pbg43LrZaae4I6mFyy+HuXPDZLBKTiKSNGYdcV8QvToW+DhNqfeB3THblTCfa3/glJreOq8SVEkJdOgQ+hk0CqNGwUMPwR/+AIceGnc0IpLvzEYAPYF2mBUD1wM9MduX0GQ3DzgvKtsJGIb70biXYXYx8ArQBHgY9+k1Xi6fmvj69Ak9+d5/P+tvnX0LFsA++8Auu4Tl25s1q/kcEZE6yKiJrwHlVS++RjOLRHk5nHFGmK388ceVnEQkL+VVgmo0s0jcfjuMGQN33AF77hl3NCIiscibBFVaCsuWNYIE9d57cM01cPzxocuhiEieypsEVbFQYaKb+FasgJNPhk6d4MEHG/mU6yIidZM3vfgSP4uEO5x/PnzxRehS3ugXrBIRqZu8SVCJX+p9+HAYMQJuvBEOOSTuaEREYpd3TXyJTFAzZ8LFF0PPnmHMk4iI5E+CKi4OEzG0ahV3JJWsWQP9+4d1QP7+d2jSJO6IREQSIa+a+BJZe7rqqjDH3j//mdAARUTikTc1qEQmqH/+E+6+Gy67DH75y7ijERFJlLxKUInqYl5SAmeeGWYqv/nmuKMREUmcvEhQZWVharvE1KA2bIBTTw33n558Epo3jzsiEZHEyYt7UAsXhuntEpOg/vznMNZp+HD4/vfjjkZEJJHyogaVqC7mEybADTfAKafA6afHHY2ISGLlRYKqmEUi9ntQS5eGxLTrrnD//ZrKSESkGnnRxJeIGpQ7nHNOuBn29tvQunWMwYiIJF/eJKimTaFduxiDGDIkrJB7661w4IExBiIi0jjU2MRnxh5mTEnZVpgxsFKZtmaMMuMjM94zY+9Kx5uY8aEZL2T7H5CJinWgtomrQXPaNLj8cjjySPjtb2MKQkSkcamxBuXOLGBfCIkGKAFGVSp2DTDFnWPN2BO4F+iVcvwyYAYQS7tWrIN0S0vhpJOgTRt45JEYs6SISONS29+WvYDP3Pmi0v5uwDgAd2YCXc3oAGDGzsAvgGF1jHWrxZqgBg6EGTPgscegQ4eYghARaXxqm6D6AyPS7J8KHAdgRg9gF6Ciz9xfgSuB8ure2MzONbNJZjaprKyslmFVzT008cXSg+/pp8PCg1ddBUccEUMAIiKNV8YJyoxmQF/g6TSHBwNtzZgCXAJ8CJSZ8UtgkTsf1PT+7j7U3bu7e/eCguz13Vi2DL79NoYa1Lx58JvfwEEHhTWeREQaO7OHMVuE2ccp+27FbCZmH2E2CrM2VZw7D7NpmE3BbFIml6tNDeooYLI7CysfcGeFO2e6sy9wOrAj8DlwCNDXjHnAk8BhZvy9Ftess1i6mK9fH5Zudw+LEDZt2oAXFxGpN8OBPpX2jQX2xv2HwKdAdYva/Rz3fXHvnsnFapOgTiZ98x5mtIlqWADnAG9GSesP7uzsTldC8+Br7pxWi2vWWSwJ6tpr4d13YejQMChXRCQXuL8JLK20bwzuFfdl3mXT7Z06y6gtzYyWwBHAeSn7zg+xMQTYC3jUjA3AJ8DZdQmqqKiI8ePH1+UtNlq5Em67LXSmy9JbVqv9uHF0u/VWSvr2ZXaHDg1zURGR7CiwzZvfhrr70FqcfxbwVBXHHBiDmQMPkMH7mrvX4toNo7Cw0FevXp2V9xo0CK6/HtauhWbNai5fJ5MnwyGHQPfuMG5cA1xQRCR7zKzU3QtrKNQVeAH3vSvtvxboDhxHusRi1gn3+Zi1JzQLXhLVyKqU84NySkpgxx0bIFcsWgTHHBMu9swzSk4ikj/MBgC/BE5Nm5wA3OdHj4sIY2l71PS2OZ+gGqSL+fr18KtfweLFYTojjXcSkXxh1ge4CuiLe2kVZQox227jc+gNfJy2bIqcT1ANMkh34EB480146CE44IB6vpiISEzMRgDvAHtgVozZ2cA9wHbA2KgL+ZCobCfMXorO7ABMwGwq8B7wIu4v13S5nJ8stqQEDj64Hi8wbBjcdx9ccUVYSkNEJFe5n5xm70NVlJ0PHB09nwv8qLaXy+ka1Jo18PXX9djE9/bbcOGF0Ls3DB5cTxcREclPOZ2g5odbcvXTxFdcDMcdB7vsAk8+CU2a1MNFRETyV0438dXbIN01a0JyWr06dCdv2zbLFxARESWo2nKH886D99+H556DH/wgi28uIiIVcrqJr7g4PGb1HtRdd8Gjj8INN0C/fll8YxERSZXTCaqkBAoLoXW2lkl89dXQW++YY+C667L0piIikk7OJ6jOncEsC282d25YGXfPPUMNSivjiojUq5z+LZu1WSRWrQrNee7w/POw3XZZeFMREalOTieorMwi4Q5nnAGffAJPPQW77ZaN0EREpAY524uvvDyMg6pzgrrpJnj22bBmh5ZtFxFpMDlbg1q0CMrK6tjEN3p06Axx2mnw299mLTYREalZziaoOo+BmjEjJKYDDggr42alp4WIiGRKCSqdZctCp4gWLcLyGS1aZDU2ERGpWc7eg9rqBLVhA5x8MsybB6+9Bl26ZDs0ERHJQM4mqOLiMH9rrdcOvPlmePllGDIEfvrTeolNRERqltNNfB071nKS8WnTwhRGv/pVmG9PRERik9MJqlbNe+vXw4ABYWby++6rt7hERCQzOd3EV6uJxv/yF/jwQxg5Etq1q7e4REQkM6pBAUyZAjfeGJZsP/bYeo1LREQyk5MJasUKWLkywwS1bl2YyqhdO7j77voOTUREMlRjgjJjDzOmpGwrzBhYqUxbM0aZ8ZEZ75mxd7S/ixmvmzHDjOlmXFZf/5BUFV3MM5pF4qabYOpUeOAB2GGHeo1LRKRRM3sYs0WYfZyyrwizsZjNjh7TLzFuNiAqMxuzAZlcrsYE5c4sd/Z1Z1/gAKAUGFWp2DXAFHd+CJwO3BXtLwN+585ewMHARWZ0yySwush4DNTkySFBnX469O1b32GJiDR2w4E+lfZdDYzDfXdgXPR6c2ZFwPXAQUAP4PoqE1mK2jbx9QI+c+eLSvu7RYHhzkygqxkd3FngzuRo/0pgBpDNBdjTyihBrV0beu116AB//Wt9hyQi0vi5vwksrbS3H/BI9PwR4Jg0Zx4JjMV9Ke7fAGPZMtFtobYJqj8wIs3+qcBxAGb0AHYBNmtgM6MrsB8wMd0bm9m5ZjbJzCaVlZXVMqzNZZSgBg2Cjz+GBx8MXctFRGRrdMB9AUD02D5Nmc7Af1JeF5NBZSXjbuZmNAP6An9Ic3gwcJcZU4BpwIeE5r2Kc1sBzwID3VmR7v3dfSgwFKBLly4+fvz4TEPbws47h0rRxLSpELabOZP9Bw/mq6OOYlbLllCHa4mI5JACM5uU8npo9Lu5rtLNtu01BlOLCxwFTHZn4RZXCUnnTAAzDPg82jCjKSE5Pe7OyEwutHTpUnr27FmL0DbXrx98/jl89FGag2vWwIUXQqdOdBwxgo7bb7/V1xERyTFl7t69lucsxKwj7gsw6wgsSlOmGOiZ8npnYHxNb1ybJr6TSd+8hxltohoWwDnAm+6siJLVQ8AMd+6oxbXqpNoxUNdfH5bSGDYMlJxEROpqNFDRK28A8HyaMq8AvTFrG3WO6B3tq1ZGCcqMlsARsKkGZMb5ZpwfvdwLmG7GTEJNq6I7+SHAr4HDUrqpH53JNeuiuLiKLubvvhtWxv3Nb+DII+s7DBGR3GI2AngH2AOzYszOJtziOQKz2YQ8MTgq2x2zYQC4LwVuBN6PtkHRvuov515jM2CDKyws9NWrV2/VuevXQ/Pm8L//G+Z93ejbb2G//cLjtGnQunVWYhURyRVmVuruhXHHUSHn5uJbsADc0zTxXXcdzJoFY8cqOYmINAI5N9VRcXF43KyJ79//hjvugAsugMMPjyUuERGpnZxLUFuMgSotDXPt7bIL3HJLXGGJiEgt5VwT3xYJ6pprYM4ceP11aNUqtrhERKR2crIG1bw5FBUBb7wBd90Fl1wCdRhXJSIiDS/nElRFF3NbvQrOOgt22y0sRigiIo1KTjbxde4MXH11mE7ijTegMDG9JkVEJEM5V4MqKYEjCl6He++FgQPh0EPjDklERLZCTg3UdYcdt13Jpy1+SFGHZvDhh9CyZT1EKCKSezRQtx4tWQI3rruStuu/gJcmKDmJiDRiOdXEt+yFCVzAED795e/gJz+JOxwREamDnEpQbY48iImn3k3h7YPiDkVEROoop+5BiYjI1kvaPaicqkGJiEjuUIISEZFEUoISEZFEUoISEZFEUoISEZFEUoISEZFEUoISEZFEUoISEZFEUoISEZFESuRMEmZWDnxbh7coAMqyFE59UpzZpTizS3FmV2OIs4W7J6biksgEVVdmNsndu8cdR00UZ3YpzuxSnNnVWOJMksRkShERkVRKUCIikki5mqCGxh1AhhRndinO7FKc2dVY4kyMnLwHJSIijV+u1qBERKSRU4ISEZFEarQJysz6mNksM5tjZlenOd7czJ6Kjk80s64xxNjFzF43sxlmNt3MLktTpqeZLTezKdH2vw0dZ0os88xsWhTHpDTHzczujj7Tj8xs/xhi3CPls5piZivMbGClMrF8pmb2sJktMrOPU/YVmdlYM5sdPbat4twBUZnZZjYghjhvNbOZ0f/rKDNrU8W51X5HGiDOG8ysJOX/9ugqzq3290MDxPlUSozzzGxKFec22OfZKLl7o9uAJsBnwHeBZsBUoFulMhcCQ6Ln/YGnYoizI7B/9Hw74NM0cfYEXoj7M41imQe0q+b40cC/AAMOBiYm4HvwFbBLEj5T4GfA/sDHKftuAa6Onl8N3JzmvCJgbvTYNnretoHj7A0URM9vThdnJt+RBojzBuCKDL4X1f5+qO84Kx2/HfjfuD/Pxrg11hpUD2COu89193XAk0C/SmX6AY9Ez58BepmZNWCMuPsCd58cPV8JzAA6N2QMWdYPeNSDd4E2ZtYxxnh6AZ+5+xcxxrCRu78JLK20O/V7+AhwTJpTjwTGuvtSd/8GGAv0acg43X2Mu1fMcvAusHN9XT9TVXyemcjk90PWVBdn9DvnRGBEfV0/lzXWBNUZ+E/K62K2/MW/sUz0g7cc2KFBoksjamLcD5iY5vCPzWyqmf3LzH7QoIFtzoExZvaBmZ2b5ngmn3tD6k/VP/hJ+Uw7uPsCCH+wAO3TlEna53oWoaacTk3fkYZwcdQU+XAVTaZJ+jwPBRa6++wqjifh80ysxpqg0tWEKveXz6RMgzCzVsCzwEB3X1Hp8GRCE9WPgP8Dnmvo+FIc4u77A0cBF5nZzyodT9Jn2gzoCzyd5nCSPtNMJOlzvZYwX9zjVRSp6TtS3+4HdgP2BRYQms8qS8znCZxM9bWnuD/PRGusCaoY6JLyemdgflVlzKwA2J6tay6oEzNrSkhOj7v7yMrH3X2Fu6+Knr8ENDWzdg0cZkUs86PHRcAoQlNJqkw+94ZyFDDZ3RdWPpCkzxRYWNEMGj0uSlMmEZ9r1Dnjl8CpHt0gqSyD70i9cveF7r7B3cuBB6u4flI+zwLgOOCpqsrE/XkmXWNNUO8Du5vZrtFf0v2B0ZXKjAYqekOdALxW1Q9dfYnanx8CZrj7HVWU2ani3piZ9SD8nyxpuCg3xlFoZttVPCfcNP+4UrHRwOlRb76DgeUVzVcxqPIv06R8ppHU7+EA4Pk0ZV4BeptZ26jJqne0r8GYWR/gKqCvu5dWUSaT70i9qnTP89gqrp/J74eGcDgw092L0x1MwueZeHH30tjajdCj7FNCb51ro32DCD9gANsSmn/mAO8B340hxp8SmhY+AqZE29HA+cD5UZmLgemEnkbvAj+J6fP8bhTD1Cieis80NVYD7o0+82lA95hibUlIONun7Iv9MyUkzAXAesJf8WcT7nuOA2ZHj0VR2e7AsJRzz4q+q3OAM2OIcw7hvk3F97SiB2wn4KXqviMNHOdj0XfvI0LS6Vg5zuj1Fr8fGjLOaP/wiu9kStnYPs/GuGmqIxERSaTG2sQnIiI5TglKREQSSQlKREQSSQlKREQSSQlKREQSSQlKREQSSQlKREQS6f8BgYC+HPFAjDQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "ax1.plot(likelihood_log1, c='b')\n",
    "ax1.tick_params(axis='y', labelcolor='b')\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(likelihood_log2, c='r')\n",
    "ax2.tick_params(axis='y', labelcolor='r')\n",
    "fig.tight_layout()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model for ground truth 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.5247179  0.03035531 0.34746964 0.09745715 0.19708422 0.63124529]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFPxJREFUeJzt3X+cVXWdx/H3Z4YBhPEHikwqJLoJ6VqOjlpqhZPomkm2a5u55QN3bemRbsWGrT/6uaspbVlUu1aWBKvmyKKl8lBLccZ20zSwQSEkXDVBSESFGBBhmM/+cc7QiINz751zzvfcc1/Px+P7OPfeuT/eXsbPfO/3+z33a+4uAED1qwsdAACQDAo6ABQEBR0ACoKCDgAFQUEHgIKgoANAQVDQAaAgKOgAUBAUdAAoiCFZvtjo0aN9/PjxFT128+bNGjlyZLKBEpBYrhUrouPEiYN/LtXA+5UwcpWHXOUbTLbFixevd/f9B7yju2fWWlpavFLt7e0VPzZNieWaNClqCSn8+5UwcpWHXOUbTDZJi7yEGsuQCwAURKZDLngDX/hC6AQAqhwFPS8mTw6dAECVY8glLzo7owYAFaKHnhfTp0fHjo6gMQBUL3roAFAQFHQAKAiGXAAgRWs2vKK2R57Vgdt7Un8teugAkKK1G7fqO/c/qRe2pF/Q6aHnxVVXhU4AoMpR0PPixBNDJwCQCs/slRhyyYsHH4wagEIyS/816KHnxeWXR0fWoQOF4tl10OmhA0A20u+iU9ABoCAo6ACQogxHXCjoAJCFDOZEmRTNjVmzQicAkIIsJ0Up6HnR3Bw6AYAUZbFskSGXvLjvvqgBQIUGLOhmNs7M2s1suZktM7PPxLfva2b3mtnK+Dgq/bgFduWVUQNQKJ7hmEspPfRuSTPc/XBJ75R0kZkdIelSSQvd/TBJC+PrAIBABizo7r7W3R+NL2+StFzSQZLOkjQ3vttcSR9MKyQAVKvcLls0s/GSjpb0sKQmd18rRUVf0pikwwFAUWSxbNFKHd8xs0ZJD0j6qrvfZmYb3H2fPj9/2d1fN45uZtMkTZOkpqamlra2toqCdnV1qbGxsaLHpimpXM3xnqKdCS1fLPr7lTRylYdcpXvipR2a+chWffpI1zFjK8vW2tq62N2PHfCO7j5gk9Qg6eeSPtvnthWSDogvHyBpxUDP09LS4pVqb2+v+LFpSizXE09ELSGFf78SRq7ykKt0Dz653g++ZIFfO/++ip9D0iIvoVaXssrFJF0vabm7f7PPj+6QNDW+PFXS7aX+tUE/Jk6MGgBUqJQTi06SdJ6kx82sM77tckkzJc0zswskPSvpb9OJWCPuvDM6TpkSNgeARHmG06IDFnR3/1/tfjz/lGTj1LBrromOFHSgkDhTFABQMgo6AKSJHYsAoFiyWIdOQQeAFGV5pihfn5sXN9wQOgGAKkdBz4tx40InAFDlGHLJi1tuiRqAQmHHolr0ve9Fx3POCZsDQCpYhw4AVS7LM0Up6ABQEBR0AMgA69ABoMoxKVqL5s8PnQBAlaOg58Xo0aETAEhBbvcURYrmzIkaAFSIgp4XFHSg0JgUBYAq5xnOilLQASALnCkKACgVBR0AUsT3odeiu+4KnQBAirKYFKWg58WIEaETAEgDe4rWoGuvjRqAQmLZYi2ZNy9qAFAhCjoApIjvQweAomEdOgBUtyy/PpeCDgAZYNliLenoCJ0AQJWjhw4AKWLIpRZ94xtRA1BIrEOvJQsWRA1AobBjEQCgbBR0ACgICjoApCjLHYtYtpgXe+wROgGAFFkGs6IU9Ly4++7QCQCkIFeTomY228zWmdnSPrd9xcyeM7POuJ2RbkwAwEBKGUOfI+n0fm7/lrs3x43tdgbriiuiBgAVGrCgu/svJb2UQZbatnBh1AAUSpZnilopM7BmNl7SAnc/Mr7+FUnnS/qTpEWSZrj7y7t57DRJ0ySpqamppa2traKgXV1damxsrOixaUoqV/P06ZKkzlmzBv1cUvHfr6SRqzzkKt3i57v13d++qkubXW99U2XZWltbF7v7sQPe0d0HbJLGS1ra53qTpHpFPfyvSppdyvO0tLR4pdrb2yt+bJoSyzVpUtQSUvj3K2HkKg+5Snf342v84EsW+Nzb76v4OSQt8hJqbEXr0N39eXff4e49kn4o6fhKngcAkJyKli2a2QHuvja++teSlr7R/VGC/fYLnQBAiiyDhegDFnQzu1nSyZJGm9lqSV+WdLKZNStaYvmMpE+kmLE23Hpr6AQAUpDlpOiABd3dz+3n5utTyAIAGAS+yyUvLrssagAKJcszRTn1Py8eeih0AgBVjh46AGSAHYsAoMqxpygAFEwWPXTG0PNi7NjQCQCkwDOcFqWg58WNN4ZOAKDKMeQCAFnIYMyFgp4X06dHDUCh5OpMUWSkszN0AgApYtkiAKBkFHQASFGuNokGAFQHxtDzYsKE0AkApMAznBWloOfFddeFTgAgRUyKAgBKRkHPi2nTogYAFWLIJS9+//vQCQCkKIMtRemhA0Ca+PpcAEDZKOgAUBCMoedFc3PoBABSwPeh16JZs0InAJAi1qEDQJVjUrQWfexjUQOACjHkkherV4dOACBFrEMHgCrHkAsAoGwUdABIUZYbXDCGnhcnnBA6AYAqR0HPi6uvDp0AQIpYhw4AVS7LHYso6Hlx9tlRA4AKMeSSFy++GDoBgBRkOSlKDx0ACoKCDgAZyMWZomY228zWmdnSPrfta2b3mtnK+Dgq3ZgAUKVydqboHEmn73LbpZIWuvthkhbG1zEYp5wSNQCFlMWyxQEnRd39l2Y2fpebz5J0cnx5rqQOSZckmKv2fPGLoRMASEGWG1xUOobe5O5rJSk+jkkuEgCgElbKove4h77A3Y+Mr29w9336/Pxld+93HN3MpkmaJklNTU0tbW1tFQXt6upSY2NjRY9NU1K53nZJ9AHn8a99bdDPJRX//UoaucpDrtI9sHq7frx0m644zjVuv8qytba2Lnb3Ywe8o7sP2CSNl7S0z/UVkg6ILx8gaUUpz9PS0uKVam9vr/ixaUos16RJUUtI4d+vhJGrPOQq3c0P/8EPvmSBz79rYcXPIWmRl1BjKx1yuUPS1PjyVEm3V/g8AFAT8rJs8WZJD0maaGarzewCSTMlnWpmKyWdGl8HAARUyiqXc3fzI9bYAcAA+D70WnTmmaETAEhRLtahIyMXXxw6AYAUsKcoAKBsFPS8OPnkqAFAhSjoAJCiajj1HwBQBvYUBYAqx6QoAKBsLFvMiw9/OHQCAGnKYMyFgp4XF14YOgGAFLBJdC3asiVqAArJMuii00PPizPOiI4dHUFjAEhYhrOi9NABoCAo6ACQAdahA0CVY1IUAIqGZYs15PzzQycAkIIszxSloOcFBR3AIDHkkhfr10cNQCGxY1Et+dCHoiPr0IFCcdahA0CxsGwRAKocyxYBAGWjoANAQTApmhef/GToBABSwDr0WnTOOaETAEiRZTArypBLXqxaFTUAhZLlpCg99Lw477zoyDp0ABWihw4ABUFBB4AUcaYoABQMZ4oCAErGpGhezJgROgGAKkdBz4spU0InAJAi1qHXkhUrogagUDhTtBZ94hPRkXXoACpEDx0ACmJQPXQze0bSJkk7JHW7+7FJhAKAovAMT/5PYsil1d3ZDBMA3gDr0AGgymU5KWqDOS3VzJ6W9LKiLxT7gbtf1899pkmaJklNTU0tbW1tFb1WV1eXGhsbK86alqRyjVq8WJL0ckvLoJ9LKv77lTRylYdcpbvr6W2at2K7rjnBtd/elWVrbW1dXNKQtrtX3CQdGB/HSFoi6T1vdP+WlhavVHt7e8WPTRO5ykOu8pCrPHnM9f2OJ/3gSxb4PffeX/FzSFrkJdTkQQ25uPua+LhO0k8lHT+Y56tpnZ1RA1AoVfF96GY2UlKdu2+KL58m6d8SS1Zrpk+PjqxDB4opg1nRwaxyaZL0U4vOZx0i6Sfufk8iqQCgIKriTFF3f0rSUQlmAQAMAssWASADrEMHgCpXbWeKIglXXRU6AYAqR0HPixNPDJ0AQAqynBRlyCUvHnwwagBQIXroeXH55dGRdehAIbFjEQCgZBR0AMgAyxYBoMp5hrOiFHQAKAgmRfNi1qzQCQBUOQp6XjQ3h04AIAWsQ69F990XNQCFlMWkKD30vLjyyug4eXLYHAASleUGF/TQAaAgKOgAkAHOFAWAKsekKACgbEyK5sUPfhA6AYAUsMFFLZo4MXQCAFWOIZe8uPPOqAEoJNah15JrromOU6aEzQEgUUyKAgDKRkEHgBT1dtAtg4XoFHQAKAgKOgAUBJOieXHDDaETAEhDhrOiFPS8GDcudAIAVY4hl7y45ZaoASicLL6YS6KHngvbd/Ro49e/rW07etQ9eYrevN+I0JEAJITvQ68hr2zboY/+6GE9uW6Tnn1xs06b9YAWLn8+dCwACcqog05BD+3qu5frkadf0qH7N6r5zaM0oWlPXfSTR7VszcbQ0QAkgDNFa0THinX6r4f+oAvedYjG7DlMw4fUafb5x2mv4Q2aMW+JtnX3hI4IoIpQ0AN5afM2fW7+Y5rYtKc+91d//qbF0Y3DdPXfvE1P/HGTvnv/yoAJASQli7NEJSZFg3B3XX7b49q4Zbvm/v3xGt5QL82fv/PnpxzepLOPGatrO/5Pkw9v0lHj9gmYFsBgZPl96PTQA7jlN6t0z7I/asZpE3TEgXtFN44eHbXYl6YcoTF7DtOM/16irdt3BEoKIAlVMSlqZqeb2Qoze9LMLk0qVJEtWbVBX7pjmd71ltH6+LsP/fMP5syJWmzvPRr0tbPfrifXdenSWx+TZzmzAiAxVTEpamb1kv5T0vskHSHpXDM7IqlgRdS5aoPO//Ej2r9xmL5z7tGqr+vzd3uXgi5J75mwvy4+bYJ+1rlGX75jmbp3MEkKYPcGM4Z+vKQn3f0pSTKzNklnSfpdEsH6Wvn8Ji1dv0PDn3pRQ4fUaWh9nYYOqVNDfZ0a6u01t9XXmeqst2U3GbE7r2zbod+t3ag7OtfopoefVdNew/WTf3yH9h05tKTHX9T6Fm18Zbt++D9Pa9EzL+vj7z5EJ71ltMbsOSz4fxuA0lTDmaIHSVrV5/pqSe8YXJz+zX3oGd24aKu06NdlP9ZMqo8LvJl2Fvy+l3uLf10/73p/Exp9P0Jt27ZNQ391X7/THlu379Cmrd2SpIZ604daxuqy9x2uvUc0lJHf9Pn3H6Gj3zxKM+9+Qp+dt0SSNLyhTo3DGrTH0Do11PX5oBX/J2zZskUjF3e85nn6/DiYzVu2aOSjDwRO8XrkKg+5Sre+69XMXmswBb2/2vC6umZm0yRNk6SmpiZ1dHSU/UJHDe3R/m9zNQzfQ909ru4eqdsVHXuvx5d7XOpRVHT7Xvadl6P7uKQe9+h6n/uX+pe0927bt/eooWHHa27rVV8njRrWoKaRdTp833o1Dn1Jv33kV/0+X/OGDZKkzt28PyMk/etx0lMbh+vpjT164ZUevdq9Q6/u6FZP/K73ffO79+jRkPqtr7k9D8PwjcN6VG+vhI7xOuQqD7lKt/ee0rgDh6irq6ui+lcWj4tauU3SCZJ+3uf6ZZIue6PHtLS0eKXa29srfmyaEss1aVLUElL49yth5CoPuco3mGySFnkJdXkwPfTfSDrMzA6R9Jykj0j6u8H9ealhd90VOgGAKldxQXf3bjP7J0k/l1Qvaba7L0ssWa0ZwTcsAhicQZ0p6u53SaJrmYRrr42OF14YNgeAqsWZonkxb17UAKBCFHQAKAgKOgAUBAUdAAqCgg4ABWGe4emDZvaCpD9U+PDRktYnGCcp5CoPucpDrvLkNZc0uGwHu/v+A90p04I+GGa2yN2PDZ1jV+QqD7nKQ67y5DWXlE02hlwAoCAo6ABQENVU0K8LHWA3yFUecpWHXOXJay4pg2xVM4YOAHhj1dRDBwC8gaoo6HncjNrMxplZu5ktN7NlZvaZ0Jn6MrN6M/utmS0InaWXme1jZvPN7In4fTshdCZJMrN/jv8Nl5rZzWY2PFCO2Wa2zsyW9rltXzO718xWxsdROcn19fjf8TEz+6mZ7ZOHXH1+drGZuZmNzksuM/tUXMeWmdm/p/HauS/oOd6MulvSDHc/XNI7JV2Uk1y9PiNpeegQu/i2pHvc/a2SjlIO8pnZQZI+LelYdz9S0VdBfyRQnDmSTt/ltkslLXT3wyQtjK9nbY5en+teSUe6+9sl/V7RBjdZm6PX55KZjZN0qqRnsw4Um6NdcplZq6I9l9/u7n8p6RtpvHDuC7r6bEbt7tsk9W5GHZS7r3X3R+PLmxQVp4PCpoqY2VhJ75f0o9BZepnZXpLeI+l6SXL3be6+IWyqnYZI2sPMhija7W9NiBDu/ktJL+1y81mS5saX50r6YKah1H8ud/+Fu3fHV38taWwecsW+Jelf1M+WmFnYTa5PSprp7q/G91mXxmtXQ0HvbzPqXBTOXmY2XtLRkh4Om2SnWYp+oXtCB+njUEkvSPpxPBT0IzMbGTqUuz+nqLf0rKS1kja6+y/CpnqNJndfK0WdCEljAufpzz9Iujt0CEkysw9Ies7dl4TOsosJkt5tZg+b2QNmdlwaL1INBb2kzahDMbNGSbdKmu7uf8pBnjMlrXP3xaGz7GKIpGMkfc/dj5a0WWGGD14jHpM+S9Ihkg6UNNLMPhY2VfUws88rGn68KQdZRkj6vKQvhc7SjyGSRikanv2cpHlmpW5JX7pqKOirJY3rc32sAn0k3pWZNSgq5je5+22h88ROkvQBM3tG0fDUe83sxrCRJEX/jqvdvfdTzHxFBT60yZKedvcX3H27pNsknRg4U1/Pm9kBkhQfU/moXgkzmyrpTEkf9Xysf/4LRX+Yl8S//2MlPWpmbwqaKrJa0m3xns+PKPr0nPiEbTUU9J2bUZvZUEUTVncEzqT4r+v1kpa7+zdD5+nl7pe5+1h3H6/ovbrf3YP3ON39j5JWmdnE+KZTJP0uYKRez0p6p5mNiP9NT1EOJmv7uEPS1PjyVEm3B8yyk5mdLukSSR9w9y2h80iSuz/u7mPcfXz8+79a0jHx715oP5P0XkkyswmShiqFLxHLfUGPJ156N6NeLmleTjajPknSeYp6wJ1xOyN0qJz7lKSbzOwxSc2SrgqcR/EnhvmSHpX0uKL/J4KcbWhmN0t6SNJEM1ttZhdIminpVDNbqWjlxsyc5PoPSXtKujf+3f9+TnIFt5tcsyUdGi9lbJM0NY1PNZwpCgAFkfseOgCgNBR0ACgICjoAFAQFHQAKgoIOAAVBQQeAgqCgA0BBUNABoCD+H/0mFU8a93guAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sensormodel import Likelihood\n",
    "\n",
    "TEST_GROUND_TRUTH = 3.0\n",
    "\n",
    "x = np.linspace(0.0, 16.0, 1000)\n",
    "x_ = np.ones_like(x)*TEST_GROUND_TRUTH\n",
    "\n",
    "likelihood1 = np.vectorize(Likelihood(theta = theta1).P)\n",
    "y1 = likelihood1(x, x_)\n",
    "\n",
    "print(theta1)\n",
    "plt.plot(x, y1)\n",
    "plt.axvline(TEST_GROUND_TRUTH, c='r', linestyle='--')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model for ground truth 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.50858697 0.01769387 0.34746964 0.12624952 0.00421572 0.72616816]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAFdZJREFUeJzt3X+0XWV95/H3NwSEEET8dY0ki4hjFZrWyKWoZQmmwUqVJchSl3REWHWMy9aIDlUQuyxDi6XWajt0ZEaEko4MWQxCdTKUFuLFtErBEKEEYgUxQjBDoIZC5Fd+fOePsy9e4g33nnP32Xuffd6vtc46++yzzz6fde6937vP8zz72ZGZSJIG36y6A0iSymFBl6SWsKBLUktY0CWpJSzoktQSFnRJaokpC3pELIiIsYjYEBF3RsQZxfpzI+KBiLituL2t/3ElSXsSU41Dj4h5wLzMXBcRBwC3AicB7wG2Zebn+x9TkjSV2VNtkJmbgc3F8mMRsQE4uN/BJEndmfII/VkbRywE1gCLgP8MnA48CqwFzszMrZO8ZhmwDGC//fYbXbBgQU9Bd+3axaxZzWvyLyvXnPvvB+DxHj+f3bX98yqbubpjru7NJNsPfvCDhzPzJVNumJnTugFz6TS3nFw8HgH2otMOfz5w6VT7GB0dzV6NjY31/Np+Ki3Xscd2biVp/edVMnN1x1zdm0k2YG1Oo05P699FROwNfA24PDOvLv4RPJiZOzNzF3AxcFS3/3UkSeWZsg09IgK4BNiQmV+YsH5edtrXAd4JrO9PxCHxB39QdwJJA27Kgg4cDZwK3BERtxXrzgFOiYjFQAIbgQ/1JeGwOO64uhNIGnDTGeXyT0BM8tS15ccZYrcV/ysXL643h6SBNZ0jdFXhYx/r3N94Y60xJA2uZo7vkSR1zYIuSS1hk0tDPLF9J0/v2MWBdQeRNLAs6A1x+/2PAPCGmnNIGlwW9Ib43DGnAXB1zTkkDS4LekOsm39Y3REkDTgLekMcsWlDsfT2WnNIGlwW9Ib45JoVxdLv15pD0uBy2KIktYQFXZJawoIuSS1hQZeklrBTtCHOW7oMcApLSb2zoDfEXSOH1h1B0oCzoDfE0RvHrx3iOHRJvbGgN8Ty76wslj5daw5Jg8tOUUlqCQu6JLWEBV2SWsI2dEnqo0cef5rvbtzKE09l39/Lgt4Q57z1IwB8s+Ycksr1w4d+xgf/Zi1njj6v7+9lQW+Ie180v+4IkgacBb0hlt5zc7HkOHSpXfrf1DLOgt4QH7zlmmLpvFpzSOqPiP6/h6NcJKmPsroDdAu6JFWj/4foFnRJ6qMKD9At6JJUhQqa0O0UbYqPn3AmADfVnENSuapsQ7egN8Tm57+k7giS+qiKUS4W9IY4YcOaYslx6JJ6Y0FviPd9b/zic39aaw5J5coK21ym7BSNiAURMRYRGyLizog4o1j/woi4PiLuLu4P6n9cSdKeTGeUyw7gzMw8DHgD8HsRcThwNrA6M18FrC4eS5ImaNSwxczcnJnriuXHgA3AwcCJwIpisxXASf0KKUmDrophi9FN+05ELATWAIuA+zLzBROe25qZv9DsEhHLgGUAIyMjoytXrtx9k2nZtm0bc+fO7em1/VRWrn3fvxyAJ//mwhnvC9r/eZXNXN0x1/Rt+Led/Ol3n+Sji5Ij5veWbcmSJbdm5pFTbpiZ07oBc4FbgZOLx4/s9vzWqfYxOjqavRobG+v5tf1UVq7Fyy/PxcsvL2Vfme3/vMpmru6Ya/q+c8/DechZq/JLV93Q8z6AtTmNOj2tUS4RsTfwNeDyzLy6WP1gRMzLzM0RMQ/Y0u1/Hf3c1jkH1h1B0oCbsqBHRACXABsy8wsTnvoGcBpwQXH/9b4kHBLvuuOGYslx6FKbZMPmQz8aOBW4IyJuK9adQ6eQXxkRHwDuA97dn4jD4ecFXVIbNeJM0cz8J/bcQbu03DiS1DLOhy5J7VLFsEULuiT1UaNOLJIkDQYn52qI0999LgDfrzeGpAFmQW+IJ/fet+4IkvrAC1wMofet+7/FkuPQpTZqxLBFVeOE7/9j3REk9UGVJxbZKSpJFXDYoiQNuCrb0C3oktQSFnRJagk7RRvivb99AQAb640hqWSeKSpJLVNFp6hH6A3xwZvHrxviOHSpTbLCXlELekMs/eEtdUeQ1E8VHKLb5CJJfWQbuiS1jCcWSZKmzTb0hnhy9vPqjiCpH5xtcfic/p7/AjgOXWorm1wkacBVOduiR+gNsfzbVxRLjkOXWsn50IfH0T++ve4IkvrA2RYlqWVsQ5ekAecRuiSpa7ahN8TW/Z5fdwRJfeRsi0Pkw+88B3AcutQ2zuUiSeqaR+gN8clvXVYsOQ5dahPnQx9CRzzwfaDzw4+oorVNUpWq+LO2yUWS+qhRbegRcWlEbImI9RPWnRsRD0TEbcXtbf2NKUmaynSO0C8Djp9k/Rczc3Fxu7bcWMOrypMQJLXLlG3ombkmIhb2P8pw23zAi+uOIKkPqjxIi+n0wBYFfVVmLioenwucDjwKrAXOzMyte3jtMmAZwMjIyOjKlSt7Crpt2zbmzp3b02v7qaxcp1/3MwAufescZpXQe9L2z6ts5uqOuabv1gd3cOH3nuLsxclrXtZbtiVLltyamUdOuWFmTnkDFgLrJzweAfai02RzPnDpdPYzOjqavRobG+v5tf1UVq5DzlqVh5y1Knfs3FXK/tr+eZXNXN0x1/T93R0/yUPOWpUrvn5Dz/sA1uY0amxPwxYz88Hx5Yi4GFjVy370c5+54csA5Pm/RTUnCUuqUhXDkXsq6BExLzM3Fw/fCax/ru01tcO33Ft3BEl9UGUb+pQFPSKuAN4MvDgiNgF/CLw5IhbTGWK5EfhQHzNKkqZhOqNcTplk9SV9yCKqPQlBUrt4pqgk9VGVB2nO5dIQ977wYACO8BBdaiXnQx8i5xy/HIB31ZxDUrm8BJ0ktYxH6EPks9ddCED+8WTT5kgaVFlhK7oFvSEO/ekDdUeQNOBscmkYZ1uUWsoLXEjSYLNTVJJaxk7RIXLXSw8FYHHNOSSVyxOLhtB5xy0D4LdrziFpcNnk0jB2ikrtkhX+UXuE3hBf/D+f7yyc5zh0Sb2xoDfEvMceBqo9CUFSdaroFLXJRZJawoLeMLahS+1UwRXoLOiS1E+NugSdqrHu4NcA8Cs155A0uCzoDfG5Y08H4P31xpBUsioHOtjkIkkt4RF6Q1x0zWcByHN/s+Ykkso03obuXC5D5KAnHq07gqQBZ5NLwzhqUWonhy1K0oBzPnRJUtdsQ2+Ibx/yWgAOs81FahXnQx9CFx59CgD/qeYckvrDybmGkUfoUqs4H/oQuuzKP+wsfMZx6JJ6Y0FviH13PFV3BEkDziaXhvECF1K7VPkXbUGXpAp4YtEQ8gIXUss06cSiiLg0IrZExPoJ614YEddHxN3F/UH9jdl+q195FKtfeVTdMST1SVOGLV4G7H4p+rOB1Zn5KmB18VgzcPHrT+bi159sC7rUMo2aDz0z1wA/3W31icCKYnkFcFLJuSRJXYrpDHqPiIXAqsxcVDx+JDNfMOH5rZk5abNLRCwDlgGMjIyMrly5sqeg27ZtY+7cuT29tp/KyrXv+5cD8NBX/isH7DPzL2dt/7zKZq7umGv6vrVpO3+9/mn+6NeSBS/qLduSJUtuzcwjp9wwM6e8AQuB9RMeP7Lb81uns5/R0dHs1djYWM+v7aeyct20YFHetGBRPvzYk6Xsr+2fV9nM1R1zTd8VN/84DzlrVV517eqe9wGszWnU2F5HuTwYEfMAivstPe5HkoZCk4ctfgM4rVg+Dfh6OXFkp6jULo06sSgirgBuAl4dEZsi4gPABcBbIuJu4C3FY0nSHjTimqKZecoenlpacpahtuo1bwLgDA/RpVap8mRBJ+dqiK8e8XYAzqg5h6TBZUFviH23P1l3BEkDzoLeEJf973MByHPfXm8QSaVq1JmikqSZa8pcLqqSnaJSq1TZKWpBl6QqNPjEIvWJB+hSu1T5N22naENc9SvHAfCJmnNIGlwW9IawoEvtFhW0uVjQG+Kgx/8d8BJ0UutU+EdtQW+Ii/72TzoLf3RyvUEk9YXDFodQlSchSOq/Rs22KEkqgcMWJWmweWLRELNTVFKv7BRtiK++7m0AnFNzDkn90YgLXKgaqw47BoBP1ZxDUrnSYYvDZ96jD9UdQVIfeYQ+RL646s8ByPPfU3MSSWVy2KIkqWsWdEnqI4ctDjGHLUrqlQVdkvpo/BgtKugVtVO0IS4+6p0AnFtvDEkDzILeEKv/w+sBC7qk3lnQG+LQf9sE2IYutY0nFg2hz/79X3UWLji13iCS+sL50CVJ02ZBbxgvcCGpVxZ0Seqj8Sb0KoYtWtAbxk5RSb2yU7QhLvz19wJwfs05JA2uGRX0iNgIPAbsBHZk5pFlhBpG3164GKh2ZjZJ/Vdlv1gZR+hLMvPhEvYz1A5/8N5i6c11xpDUJ86HPkQ+s/rLxdLv1JpDUrmq7BeLmZzFFBE/ArbSaSn4H5n55Um2WQYsAxgZGRlduXJlT++1bds25s6d23PWfikr177vXw7Axov+kpftP/O+6rZ/XmUzV3fMNX3X/uhprvzX7fz5G5MXHdhbtiVLltw6rSbtzOz5Bry8uH8pcDtwzHNtPzo6mr0aGxvr+bX9VFaumxYsypsWLMofbnmslP21/fMqm7m6Y67pu+jGe/KQs1blddd/s+d9AGtzGjV5RoeCmfmT4n4LcA1w1Ez2JztFJfWu54IeEftHxAHjy8BvAuvLCiZJrdLw+dBHgGuic/rTbOB/ZeZ1paQaQp875jQA/qzmHJLKVWWnaM8FPTPvBV5bYpahtm7+YYBnikpt5bDFIXLEpg3F0rG15pBUrkE7sUgl+OSaFcXSh2rNIWlwOTmXJPVRlc2oFvSGsQ1dUq8s6JJUAedDlyRNm52iDXHe0mUA/EXNOST1h8MWh8hdI4fWHUFSH2SFHWMW9IY4euNtAGQeU3MSSYPKgt4Qy78zPq3wR2vNIalcDlscYlWeVSapXSzoklSBKjpFLeiS1EdVfue2oDeMZ4pK7VTFiUV2ijbEOW/9CAD/reYckso1EPOhq1z3vmg+4BG6pN5Z0Bti6T03F0tvqjWHpHI5H/oQ+uAt1wCQnFlzEkmDyk5RSeqj8WZUhy1KkqbNgt4wdopK7RQVjFu0oEtSH1V5jGanaEN8/IROZ+jFNeeQNLgs6A2x+fkvAWDHLttcpFZxPvThc8KGNQBs3/nGmpNIGlQW9IZ43/euBWD7jk/UnERSmZJq5nEBO0Ub5+mdu+qOIGlAWdAbZvtO29CltqnoAN2C3jTbPUKXWsVL0A0xC7rUPlWcVAR2ijbGh0/6FADn7LCgS23ibItDaOucAwHb0CX1zoLeEO+64wYAtr/jl2tOIqlMmQPSKRoRx0fEv0bEPRFxdlmhhtG716/mXXfcYBu6pJ71XNAjYi86l8D8LeBw4JSIOLysYBM9tWMnT+1Mnty+k6d37GLHzl1ky6YlHO8zcRy61D5VnVg0kyaXo4B7MvNegIhYCZwI3FVGsIn+eNUG/uc/Pw7XX/cLz80KmBXRuc36+XIU6/eaFcyKTi/zrICY8OVn4occz6yb/JN/1rYTlp984kn2u2Vs0n3tvr+YZIPxxfOLOVwu+ccfcc26BybN0I2fPf44+6/71oz3UzZzdcdc3Wliroe3PVXZe82koB8M3D/h8Sbg9btvFBHLgGUAIyMj3HjjjV2/0ciOnZy4MNl7n33IhF3ZOZ02E3aN3z+zLp9ZvyuTpPN4fN24yQ7wc4rnf2Ebkh177WL27Kcmff5Zr8vnfv6AfYKdCa+YuxN4Yg9bTd/c5+1ir5j5fspmru6YqztNzHXgAbDg5bPZtm1bT/WvK50C2P0NeDfwlQmPTwUufK7XjI6OZq/GxsZ6fm0/lZbr2GM7t5K0/vMqmbm6Y67uzSQbsDanUZdncoS+CVgw4fF84Ccz2N9wu/bauhNIGnAzGeXyXeBVEfGKiNgHeC/wjXJiDaE5czo3SepRz0fombkjIj4C/D2wF3BpZt5ZWrJh86Uvde5/93frzSFpYM3oxKLMvBawraAMV17ZubegS+qRk3NJUktY0CWpJSzoktQSFnRJaonICudEiYiHgB/3+PIXAw+XGKcs5uqOubpjru40NRfMLNshmfmSqTaqtKDPRESszcwj686xO3N1x1zdMVd3mpoLqslmk4sktYQFXZJaYpAK+pfrDrAH5uqOubpjru40NRdUkG1g2tAlSc9tkI7QJUnPwYIuSS0xEAW9iRejjogFETEWERsi4s6IOKPuTBNFxF4R8b2IWFV3lnER8YKIuCoivl98bm+sOxNARHy8+Bmuj4grImLfmnJcGhFbImL9hHUvjIjrI+Lu4v6ghuT6s+Ln+C8RcU1EvKAJuSY89/sRkRHx4qbkiojlRR27MyI+14/3bnxBr/Ji1F3aAZyZmYcBbwB+ryG5xp0BbKg7xG7+ErguM18DvJYG5IuIg4GPAkdm5iI6U0G/t6Y4lwHH77bubGB1Zr4KWF08rtpl/GKu64FFmfmrwA+AT1UdislzERELgLcA91UdqHAZu+WKiCV0rrn8q5n5y8Dn+/HGjS/oTLgYdWY+DYxfjLpWmbk5M9cVy4/RKU4H15uqIyLmA28HvlJ3lnER8XzgGOASgMx8OjMfqTfVM2YD+0XEbGAONV15KzPXAD/dbfWJwIpieQVwUqWhmDxXZv5DZu4oHv4znSuW1Z6r8EXgk+z5Er59tYdcHwYuyMynim229OO9B6GgT3Yx6kYUznERsRB4HXBzvUme8Rd0fqF3TbVhhQ4FHgL+umgK+kpE7F93qMx8gM7R0n3AZuDfM/Mf6k31LCOZuRk6BxHAS2vOM5nfAf6u7hAAEfEO4IHMvL3uLLv5JeBNEXFzRHwrIn6tH28yCAU9JlnXmLGWETEX+Brwscx8tAF5TgC2ZOatdWfZzWzgCOCizHwd8DPqaT54lqJN+kTgFcDLgf0j4n31phocEfFpOs2Plzcgyxzg08Bn6s4yidnAQXSaZz8BXBkRk9W2GRmEgt7Yi1FHxN50ivnlmXl13XkKRwPviIiNdJqnfiMivlpvJKDzc9yUmePfYq6iU+Drdhzwo8x8KDO3A1cDv15zpokejIh5AMV9X76q9yIiTgNOAP5jNuOEllfS+cd8e/H7Px9YFxEvqzVVxybg6uy4hc6359I7bAehoDfyYtTFf9dLgA2Z+YW684zLzE9l5vzMXEjns/pmZtZ+xJmZ/w+4PyJeXaxaCtxVY6Rx9wFviIg5xc90KQ3orJ3gG8BpxfJpwNdrzPKMiDgeOAt4R2Y+XncegMy8IzNfmpkLi9//TcARxe9e3f4W+A2AiPglYB/6MCtk4wt60fEyfjHqDcCVDbkY9dHAqXSOgG8rbm+rO1TDLQcuj4h/ARYDn605D8U3hquAdcAddP4majl9PCKuAG4CXh0RmyLiA8AFwFsi4m46IzcuaEiuvwIOAK4vfvf/e0Ny1W4PuS4FDi2GMq4ETuvHtxpP/Zeklmj8EbokaXos6JLUEhZ0SWoJC7oktYQFXZJawoIuSS1hQZeklvj/qv0cfm4o4PAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "likelihood2 = np.vectorize(Likelihood(theta = theta2).P)\n",
    "y2 = likelihood2(x, x_)\n",
    "\n",
    "print(theta2)\n",
    "plt.plot(x, y2)\n",
    "plt.axvline(TEST_GROUND_TRUTH, c='r', linestyle='--')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test Model for Correctness\n",
    "We now sample from the fitted model and plot the histogram so that we can compare with the likelihood to see if our implementation is correct.\n",
    "### Samples for ground truth 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAEKdJREFUeJzt3X+QVeV9x/HPN2hGYRvFELdEiEs6SmOZhGSp9cc0cQO2aPBHmxjiFAcndtYxiZWObUFi2kx1DNNqQmfaWJ1IYdS6OkSroqZButR01LRAUbCIZoyDEBqCP4irpqntt3+cs3bdLN57nnvvPs95eL9mnjn39/24u3722eecezB3FwCg/t4VOwAAoD0odADIBIUOAJmg0AEgExQ6AGSCQgeATFDoAJAJCh0AMkGhA0AmDhvPN5syZYr39PQEPfe1117TpEmT2huoDdqWa+fOYjtzZuuvpUPg69Vm5KqGXNW1km3z5s373f19DR/o7uM2ent7PdTg4GDwczupbbk+8YlitEn2X682I1c15KqulWySNnkTHcuSCwBkYlyXXPAOrr46dgIANUehp2LevNgJANQcSy6p2Lq1GAAQiBl6KpYsKbYbN0aNAaC+mKEDQCYodADIBIUOAJmg0AEgE+wUTcV118VOAKBDepY9oNXzO39KAgo9FaedFjsBgJpjySUVjz5aDAAIxAw9FcuXF1uOQwcQiBk6AGSCQgeATFDoAJAJCh0AMsFO0VSsXBk7AYCao9BTMXt27AQAao4ll1Q8/HAxACAQM/RUXHttseVfLgIQiBk6AGSiYaGb2XQzGzSzHWb2lJldUd5+jJmtN7Nny+3kzscFABxMMzP0NyVd6e4fknSKpC+a2UmSlkna4O4nSNpQXgcARNKw0N19r7tvKS+/KmmHpOMknSdpTfmwNZLO71RIAEBj5u7NP9isR9IjkmZJ2uXuR4+472V3/4VlFzPrl9QvSd3d3b0DAwNBQYeGhtTV1RX03E5qV64jd+2SJL3xgQ+0/FpS/l+vdiNXNeSqZtueA5px1ITgbH19fZvdfU7DB7p7U0NSl6TNkn63vP7KqPtfbvQavb29HmpwcDD4uZ1ErmrIVQ25qkk11/FL17WUTdImb6KnmzrKxcwOl/RtSbe7+93lzT82s6nl/VMl7av2Owdvc//9xQCAQM0c5WKSbpG0w92/PuKu+yQtLi8vlnRv++MdQm64oRgAEKiZDxadLukiSdvMbGt523JJKyTdZWaXSNol6YLORAQANKNhobv7v0iyg9w9t71xAACh+KQoAGSCQgeATHByrlTcemvsBABqjkJPxfTpsRMAqDmWXFJx553FAIBAzNBTceONxXbhwrg5ANQWM3QAyASFDgCZoNABIBMUOgBkgp2iqVi7NnYCADVHoadiypTYCQDUHEsuqVi9uhgAEIhCTwWFDqBFFDoAZIJCB4BMUOgAkAkKHQAywWGLqXjwwdgJANQchZ6KiRNjJwBQcyy5pOKb3ywGAASi0FNx113FAIBAFDoAZIJCB4BMUOgAkAkKHQAywWGLqdi4MXYCADXHDB0AMkGhp+L664sBAIEo9FSsW1cMAAhEoQNAJih0AMgEhQ4AmeCwxVQceWTsBABqjkJPxUMPxU4AoOZYcgGATFDoqbjmmmIAQCAKPRUbNhQDAAI1LHQzW2Vm+8xs+4jbvmpme8xsaznO7mxMAEAjzczQV0uaP8bt33D32eXgXzgGgMgaFrq7PyLppXHIAgBoQStr6F8ysyfLJZnJbUt0qHrve4sBAIHM3Rs/yKxH0jp3n1Ve75a0X5JLukbSVHf//EGe2y+pX5K6u7t7BwYGgoIODQ2pq6sr6LmdRK5qyFUNuapJNde2PQc046gJwdn6+vo2u/uchg9094ZDUo+k7VXvGz16e3s91ODgYPBzO4lc1ZCrGnJVk2qu45euaymbpE3eRMcGLbmY2dQRV39H0vaDPRZNuuqqYgBAoIYf/TezOySdIWmKme2W9GeSzjCz2SqWXJ6XdGkHMx4aHnssdgIANdew0N39wjFuvqUDWQAALeCTogCQCQodADLB6XNTMW1a7AQAao5CT8Vtt8VOAKDmWHIBgExQ6KlYsqQYABCIJZdUbN0aOwGAmmOGDgCZoNABIBMUOgBkgjX0VJx4YuwEAGqOQk/FzTfHTgCg5lhyAYBMUOip6O8vBgAEYsklFc88EzsBgJpjhg4AmaDQASATFDoAZII19FTMnh07AYCao9BTsXJl7AQAao4lFwDIBIWeikWLigEAgVhyScXu3bETAKg5ZugAkAkKHQAyQaEDQCZYQ0/FqafGTgCg5ij0VHzta7ETAKg5llwAIBMUeio+/eliAEAgllxS8eKLsRMAqDlm6ACQCQodADJBoQNAJlhDT8XcubETAKg5Cj0VX/lK7AQAao4lFwDIBIWeirPOKgYABGLJJRVvvBE7AYCaazhDN7NVZrbPzLaPuO0YM1tvZs+W28mdjQkAaKSZJZfVkuaPum2ZpA3ufoKkDeV1AEBEDQvd3R+R9NKom8+TtKa8vEbS+W3OBQCoKHQNvdvd90qSu+81s2PbmOnQtGBB7AQAas7cvfGDzHokrXP3WeX1V9z96BH3v+zuY66jm1m/pH5J6u7u7h0YGAgKOjQ0pK6urqDndhK5qiFXNeSqJtVc2/Yc0IyjJgRn6+vr2+zucxo+0N0bDkk9kraPuL5T0tTy8lRJO5t5nd7eXg81ODgY/NxOIlc15KqGXNWkmuv4petayiZpkzfRsaHHod8naXF5ebGkewNfB8POOKMYABComcMW75D0mKSZZrbbzC6RtELSmWb2rKQzy+sAgIga7hR19wsPchdnkwKAhPDRfwDIBIUOAJngXC6p+OxnYycAUHMUeiq+8IXYCQDUHEsuqXj99WIAQCBm6Kk4++xiu3Fj1BgA6osZOgBkgkIHgExQ6ACQCQodADLBTtFUXHxx7AQAao5CTwWFDqBFLLmkYv/+YgBAIGboqfjMZ4otx6EDCMQMHQAyQaEDQCYodADIBIUOAJlgp2gqLrssdgIANUehp2LhwtgJANQcSy6peOGFYgBAIGboqbjoomLLcegAAjFDB4BMUOgAkAkKHQAyQaEDQCbYKZqKK6+MnQBAzVHoqTjnnNgJANQcSy6p2LmzGAAQiBl6Ki69tNhyHDqAQMzQASATFDoAZIJCB4BMUOgAkAl2iqbi6qtjJwBQcxR6KubNi50AQM2x5JKKrVuLAQCBmKGnYsmSYstx6AACMUMHgEy0NEM3s+clvSrpfyS96e5z2hEKAFBdO5Zc+tx9fxteBwDQApZcACAT5u7hTzb7oaSXJbmkm9z95jEe0y+pX5K6u7t7BwYGgt5raGhIXV1dwVk7pV253rN9uyTpp7NmtfxaUv5fr3YjVzXkqmbbngOacdSE4Gx9fX2bm1rSdvfgIen95fZYSU9I+vg7Pb63t9dDDQ4OBj+3k8hVDbmqIVc1qeY6fum6lrJJ2uRNdHJLSy7u/qNyu0/SPZJObuX1DmmPPloMAAgUvFPUzCZJepe7v1pe/i1Jf962ZIea5cuLLcehAwjUylEu3ZLuMbPh1/l7d/9OW1IBACoLLnR3f07SR9qYBQDQAg5bBIBMUOgAkAlOzpWKlStjJwBQcxR6KmbPjp0AQM2x5JKKhx8uBgAEYoaeimuvLbb8y0UAAjFDB4BMUOgAkAkKHQAyQaEDQCbYKZqKm26KnQBAzVHoqZg5M3YCADXHkksq7r+/GAAQiBl6Km64odiec07cHABqixk6AGSCQgeATFDoAJAJCh0AMsFO0VTcemvsBABqjkJPxfTpsRMAqDmWXFJx553FAIBAzNBTceONxXbhwrg5ANQWM3QAyASFDgCZoNABIBMUOgBkgp2iqVi7NnYCADVHoadiypTYCQDUHEsuqVi9uhgAEIhCTwWFDqBFFDoAZIJCB4BMUOiJePy5F/X4cy+qZ9kDsaMAqCkKPQGjS5xSBxCCwxYjGy7viy/4atwgAGqPGXpEI2fiPzv8CP3s8CPGvA8AmkGhJ2LRlge0aAtLLwDCUeiRjC7rBU9/Twue/l6kNABy0FKhm9l8M9tpZj8ws2XtCoX/xywdQLOCC93MJkj6G0lnSTpJ0oVmdlK7guWsaklT6gCa0coM/WRJP3D359z955IGJJ3Xnlh56ln2QHA5U+oAGmnlsMXjJL0w4vpuSb/RWpyD27bngC5+h1J7fsWnOvXWLWtXGQ+/Tsr/rQDiMXcPe6LZBZJ+291/v7x+kaST3f3yUY/rl9RfXp0paWdg1imS9gc+t5PIVQ25qiFXNanmklrLdry7v6/Rg1qZoe+WNH3E9WmSfjT6Qe5+s6SbW3gfSZKZbXL3Oa2+TruRqxpyVUOualLNJY1PtlbW0P9N0glmNsPM3i3pc5Lua08sAEBVwTN0d3/TzL4k6R8lTZC0yt2falsyAEAlLZ3Lxd0flPRgm7I00vKyTYeQqxpyVUOualLNJY1DtuCdogCAtPDRfwDIRC0KPcVTDJjZdDMbNLMdZvaUmV0RO9NIZjbBzP7dzNbFzjLMzI42s7Vm9nT5dTs1diZJMrM/LL+H283sDjM7ovGzOpJjlZntM7PtI247xszWm9mz5XZyIrn+svw+Pmlm95jZ0SnkGnHfH5mZm9mUVHKZ2eVljz1lZn/RifdOvtATPsXAm5KudPcPSTpF0hcTyTXsCkk7YocY5a8kfcfdf1XSR5RAPjM7TtIfSJrj7rNU7OD/XKQ4qyXNH3XbMkkb3P0ESRvK6+NttX4x13pJs9z9w5KekXTVeIfS2LlkZtMlnSlp13gHKq3WqFxm1qfik/Qfdvdfk3R9J944+UJXoqcYcPe97r6lvPyqinI6Lm6qgplNk/QpSd+KnWWYmb1H0scl3SJJ7v5zd38lbqq3HCbpSDM7TNJEjfF5ivHg7o9IemnUzedJWlNeXiPp/HENpbFzuft33f3N8urjKj6HEj1X6RuS/kRSlB2EB8l1maQV7v5f5WP2deK961DoY51iIIniHGZmPZI+Kun7cZO8ZaWKH+j/jR1khA9K+omkvyuXgr5lZpNih3L3PSpmS7sk7ZV0wN2/GzfV23S7+16pmERIOjZynrF8XtJDsUNIkpmdK2mPuz8RO8soJ0r6TTP7vpn9s5n9eifepA6FbmPclsyhOWbWJenbkpa4+08TyLNA0j533xw7yyiHSfqYpBvd/aOSXlOc5YO3Kdekz5M0Q9L7JU0ys0VxU9WHmX1ZxfLj7QlkmSjpy5L+NHaWMRwmabKK5dk/lnSXmY3VbS2pQ6E3dYqBGMzscBVlfru73x07T+l0Seea2fMqlqc+aWa3xY0kqfg+7nb34b9i1qoo+NjmSfqhu//E3f9b0t2SToucaaQfm9lUSSq3HflTPYSZLZa0QNLveRrHP/+Kil/MT5Q//9MkbTGzX46aqrBb0t1e+FcVfz23fYdtHQo9yVMMlL9db5G0w92/HjvPMHe/yt2nuXuPiq/VP7l79Bmnu/+npBfMbGZ501xJ/xEx0rBdkk4xs4nl93SuEthZO8J9khaXlxdLujdilreY2XxJSyWd6+6vx84jSe6+zd2Pdfee8ud/t6SPlT97sf2DpE9KkpmdKOnd6sBJxJIv9HLHy/ApBnZIuiuRUwycLukiFTPgreU4O3aoxF0u6XYze1LSbEnXRc6j8i+GtZK2SNqm4v+JKJ82NLM7JD0maaaZ7TazSyStkHSmmT2r4siNFYnk+mtJvyRpffmz/7eJ5IruILlWSfpgeSjjgKTFnfirhk+KAkAmkp+hAwCaQ6EDQCYodADIBIUOAJmg0AEgExQ6AGSCQgeATFDoAJCJ/wN0Fe0Q1U4HbQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sensormodel import SampleLIDAR\n",
    "\n",
    "N_SAMPLES = 10**5\n",
    "ground_truth = np.ones(N_SAMPLES)*TEST_GROUND_TRUTH\n",
    "\n",
    "sampler1 = SampleLIDAR(theta1)\n",
    "\n",
    "samples1 = sampler1.sample(ground_truth)\n",
    "plt.hist(samples1, bins=1000, density=True);\n",
    "plt.axvline(TEST_GROUND_TRUTH, c='r', linestyle='--')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Samples for ground truth 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAErtJREFUeJzt3X+QXWV9x/H3t4EqsBbQyBohQ7CFVEs1utRBmDYsP5xIQbRakakYR2bW8VdJC5YfMuqMVLEFoTOtVEYQRMrCRCyaAgp0kVGoLcEIwRSxlGIQjShQFx0t7bd/nBNnjYl777nn5px79v2aeebsPffcu5+5ST777HPPPYnMRJI0+n6t6QCSpHpY6JLUERa6JHWEhS5JHWGhS1JHWOiS1BEWuiR1hIUuSR1hoUtSR+yyM7/Z4sWLc9myZZUe+9RTT7HHHnvUG6gGteW6//5iu3z54M/FAni9amau/pirf4NkW79+/WOZ+dx5D8zMnTYmJiayqpmZmcqPHabacq1cWYyadP71qpm5+mOu/g2SDbgre+hYl1wkqSN26pKLfoVzzmk6gaQRZ6G3xdFHN51A0ohzyaUtNmwohiRVNO8MPSKeCdwOPKM8fm1mvj8iDgCmgWcDdwMnZ+bPhhm209asKba33dZoDEmjq5cZ+k+BIzPzJcAKYFVEHAp8BLgwMw8EHgdOGV5MSdJ85i308qyZ2fLmruVI4Ehgbbn/CuA1Q0koSepJT2voEbEoIjYAW4Cbgf8AnsjMp8tDNgP7DieiJKkXkX38n6IRsRfwWeB9wCcz87fK/UuBGzLzd7fzmClgCmB8fHxienq6UtDZ2VnGxsYqPXaY6sq1olxD33DRRQM/F3T/9aqbufpjrv4Nkm1ycnJ9Zh4y74G9fPpo7gDeD7wHeAzYpdz3CuAL8z3WT4r+Cl/5SjFq0vnXq2bm6o+5+teKT4pGxHPLmTkRsRtwNLAJmAFeXx62Gri+3586muOww4ohSRX18sGiJcAVEbGIYs392sxcFxHfAKYj4lzga8ClQ8zZfXfcUWwtdUkVzVvomXkP8NLt7H8QePkwQi1IZ59dbD0PXVJFflJUkjrCQpekjrDQJakjLHRJ6ggvn9sWNX2gSNLCZaG3xYoVTSeQNOJccmmLW24phiRV5Ay9Lc49t9j6PxdJqsgZuiR1hIUuSR1hoUtSR1joktQRvinaFh//eNMJJI04C70tli9vOoGkEeeSS1t8/vPFkKSKnKG3xQUXFNvjj282h6SR5QxdkjrCQpekjrDQJakjLHRJ6gjfFG2LK69sOoGkEWeht8XSpU0nkDTiXHJpi2uuKYYkVeQMvS0uvrjYnnhiszkkjSxn6JLUERa6JHXEvIUeEUsjYiYiNkXEfRFxarn/AxHxSERsKMexw48rSdqRXtbQnwZOy8y7I+JZwPqIuLm878LMPH948SRJvZq30DPzUeDR8usfRcQmYN9hB1tw1q5tOoGkEdfXGnpELANeCny13PWuiLgnIi6LiL1rzrawLF5cDEmqKDKztwMjxoAvAX+ZmddFxDjwGJDAB4ElmfnW7TxuCpgCGB8fn5ienq4UdHZ2lrGxsUqPHaa6cj3vppsA+O6qVQM/F3T/9aqbufpjrv4Nkm1ycnJ9Zh4y74GZOe8AdgW+APz5Du5fBmyc73kmJiayqpmZmcqPHabacq1cWYyadP71qpm5+mOu/g2SDbgre+jqXs5yCeBSYFNmfnTO/iVzDnstsLH3nzeSpLr1cpbL4cDJwL0RsaHcdzZwUkSsoFhyeQh421ASSpJ60stZLl8GYjt33VB/HElSVX5SVJI6wotztcUN/sIjaTAWelvsvnvTCSSNOJdc2uJjHyuGJFVkobfFtdcWQ5IqstAlqSMsdEnqCAtdkjrCQpekjvC0xba47bamE0gacc7QJakjLPS2OP/8YkhSRRZ6W6xbVwxJqshCl6SOsNAlqSMsdEnqCE9bbIvddms6gaQR5wy9LW68sRiSOmfZmf+0U76PhS5JHWGht8UHP1gMSarIQm+LW28thiRVZKFLUkdY6JLUERa6JHWE56G3xXOe03QCSSPOQm+Lz3ym6QSSRty8Sy4RsTQiZiJiU0TcFxGnlvufHRE3R8QD5Xbv4ceVJO1IL2voTwOnZeYLgUOBd0bEi4AzgVsz80Dg1vK2qjrrrGJIUkXzLrlk5qPAo+XXP4qITcC+wAnAEeVhVwC3AWcMJeVCcOedTSeQNOIiM3s/OGIZcDtwMPBwZu41577HM/OXll0iYgqYAhgfH5+Ynp6uFHR2dpaxsbFKjx2munKtWLMGgA0XXTTwc0H3X6+6mas/5urPvY88yQF7LqqcbXJycn1mHjLvgZnZ0wDGgPXAH5W3n9jm/sfne46JiYmsamZmpvJjh6m2XCtXFqMmnX+9amau/pirP/ufsW6gbMBd2UNP93QeekTsCnwGuCozryt3fy8ilpT3LwG29PczR5JUp17OcgngUmBTZn50zl2fA1aXX68Grq8/3gKy337FkKSKejkP/XDgZODeiNhQ7jsbOA+4NiJOAR4G/ng4EReIT3+66QSSRlwvZ7l8GYgd3H1UvXEkSVV5LZe2WLOmGJJUkR/9b4sNG+Y/RpJ+BWfoktQRFrokdYSFLkkd4Rp6Wxx0UNMJJI04C70tLrmk6QSSRpxLLpLUERZ6W0xNFUOSKnLJpS2++c2mE0gacc7QJakjLHRJ6ggLXZI6wjX0tlixoukEkkachd4WNf1fopIWLpdcJKkjLPS2eNObiiFJFbnk0habNzedQNKIc4YuSR1hoUtSR1joktQRrqG3xSte0XQCSSPOQm+LD3+46QSSRpxLLpLUERZ6W7zudcWQpIrmLfSIuCwitkTExjn7PhARj0TEhnIcO9yYC8APflAMSaqolxn65cCq7ey/MDNXlOOGemNJkvo1b6Fn5u3AD3dCFknSAAZZQ39XRNxTLsnsXVsiSVIlkZnzHxSxDFiXmQeXt8eBx4AEPggsycy37uCxU8AUwPj4+MT09HSloLOzs4yNjVV67DDVlWv/T30KgP9685sHfi7o/utVN3P1x1z9ufeRJzlgz0WVs01OTq7PzEPmPTAz5x3AMmBjv/dtOyYmJrKqmZmZyo8dJnP1x1z9MVd/2ppr/zPWDZQNuCt76NhKSy4RsWTOzdcCG3d0rCRp55j3k6IRcTVwBLA4IjYD7weOiIgVFEsuDwFvG2LGheFVryq2N97YbA5JI2veQs/Mk7az+9IhZFnYfvKTphNIGnF+UlSSOsJCl6SOsNAlqSO8fG5bHHdc0wkkjTgLvS1OP73pBJJGnEsuktQRFnpbHHFEMSSpIgtdkjrCQpekjrDQJakjLHRJ6ghPW2yLN7yh6QSSRpyF3hbveEfTCSSNOJdc2uLHPy6GJFXkDL0tjj222N52W6MxJI0uZ+iS1BEWuiR1hIUuSR1hoUtSR/imaFu85S1NJ5A04iz0trDQJQ3IJZe2eOyxYkhSRc7Q2+L1ry+2nocuqSJn6JLUERa6JHWEhS5JHTFvoUfEZRGxJSI2ztn37Ii4OSIeKLd7DzemJGk+vczQLwdWbbPvTODWzDwQuLW8rUG8/e3FkKSK5j3LJTNvj4hl2+w+ATii/PoK4DbgjBpzLTwnnth0AkkjLjJz/oOKQl+XmQeXt5/IzL3m3P94Zm532SUipoApgPHx8Ynp6elKQWdnZxkbG6v02GGqK9cztmwB4Kf77DPwc0H3X6+6mas/5urPvY88yQF7LqqcbXJycn1mHjLvgZk57wCWARvn3H5im/sf7+V5JiYmsqqZmZnKjx2m2nKtXFmMmnT+9aqZufpjrv7sf8a6gbIBd2UPHVv1LJfvRcQSgHK7peLzSJJqUrXQPwesLr9eDVxfTxxJUlW9nLZ4NXAnsDwiNkfEKcB5wDER8QBwTHlbktSgXs5yOWkHdx1VcxZJ0gC8OFdbnHZa0wkkjTgLvS2OP77pBJJGnNdyaYv77y+GJFXkDL0t3va2Yuv10CVV5AxdkjrCQpekjrDQJakjLHRJ6gjfFG2Lc85pOoGkEWeht8XRRzedQNKIc8mlLTZsKIYkVeQMvS3WrCm2nocuqSJn6JLUERa6JHWEhS5JHWGhS1JH+KZoW3zoQ00nkDTiLPS2OOywphNIGnEuubTFHXcUQ5IqcobeFmefXWw9D11SRc7QJakjLHRJ6ggLXZI6wkKXpI7wTdG2uOiiphNIGnEDFXpEPAT8CPhf4OnMPKSOUAvSihVNJ5A04uqYoU9m5mM1PM/Cdsstxdb/6EJSRS65tMW55xZbC11SRYO+KZrAFyNifURM1RFIklRNZGb1B0c8PzO/ExH7ADcD787M27c5ZgqYAhgfH5+Ynp6u9L1mZ2cZGxurnHVY6sq1ovwfizbU9OZo11+vupmrP+bqz72PPMkBey6qnG1ycnJ9T+9RZmYtA/gAcPqvOmZiYiKrmpmZqfzYYaot18qVxahJ51+vmpmrP+bqz/5nrBsoG3BX9tDDlZdcImKPiHjW1q+BVwIbqz6fJGkwg7wpOg58NiK2Ps8/ZOZNtaRaiD7+8aYTSBpxlQs9Mx8EXlJjloVt+fKmE0gacX70vy0+//liSFJFnofeFhdcUGyPP77ZHJJGljN0SeoIC12SOsJCl6SOsNAlqSN8U7Qtrryy6QSSRpyF3hZLlzadQNKIc8mlLa65phiSVJEz9La4+OJie+KJzeaQNLKcoUtSR1joktQRFrokdYSFLkkd4ZuibbF2bdMJJI04C70tFi9uOoGkEeeSS1tcfnkxJKkiC70tLHRJA7LQJakjLHRJ6ggLXZI6wkKXpI7wtMW2uOGGphNIGnEWelvsvnvTCSSNOJdc2uJjH+OcV76j6RSSRpgz9La49lqOe/AHTaeQNMIGmqFHxKqIuD8ivhURZ9YVSpLUv8qFHhGLgL8DXgW8CDgpIl5UVzBJUn8GmaG/HPhWZj6YmT8DpoET6oklSerXIIW+L/DtObc3l/skSQ0Y5E3R2M6+/KWDIqaAqfLmbETcX/H7LQYeq/jYYao3V2zvZa1kYbxe9TFXf8zVp8mPDJRt/14OGqTQNwNL59zeD/jOtgdl5iXAJQN8HwAi4q7MPGTQ56mbufpjrv6Yqz9tzQU7J9sgSy7/BhwYEQdExK8DbwQ+V08sSVK/Ks/QM/PpiHgX8AVgEXBZZt5XWzJJUl8G+mBRZt4A7KyLkAy8bDMk5uqPufpjrv60NRfshGyR+UvvY0qSRpDXcpGkjhiJQm/jJQYiYmlEzETEpoi4LyJObTrTXBGxKCK+FhHrms6yVUTsFRFrI+Lfy9ftFU1nAoiIPyv/DDdGxNUR8cyGclwWEVsiYuOcfc+OiJsj4oFyu3dLcv11+ed4T0R8NiL2akOuOfedHhEZEYvbkisi3l322H0R8VfD+N6tL/QWX2LgaeC0zHwhcCjwzpbk2upUYFPTIbbxN8BNmfnbwEtoQb6I2Bf4U+CQzDyY4g3+NzYU53Jg1Tb7zgRuzcwDgVvL2zvb5fxyrpuBgzPzxcA3gbN2dii2n4uIWAocAzy8swOVLmebXBExSfFJ+hdn5u8A5w/jG7e+0GnpJQYy89HMvLv8+kcU5dSKT8pGxH7AHwKfaDrLVhHxG8AfAJcCZObPMvOJZlP93C7AbhGxC7A72/k8xc6QmbcDP9xm9wnAFeXXVwCv2amh2H6uzPxiZj5d3vwXis+hNJ6rdCHwF2zng447ww5yvR04LzN/Wh6zZRjfexQKvfWXGIiIZcBLga82m+TnLqL4C/1/TQeZ4wXA94FPlktBn4iIPZoOlZmPUMyWHgYeBZ7MzC82m+oXjGfmo1BMIoB9Gs6zPW8Fbmw6BEBEvBp4JDO/3nSWbRwE/H5EfDUivhQRvzeMbzIKhd7TJQaaEhFjwGeANZn53y3IcxywJTPXN51lG7sALwMuzsyXAk/RzPLBLyjXpE8ADgCeD+wREW9qNtXoiIj3Uiw/XtWCLLsD7wXe13SW7dgF2JtiefY9wLUR9V3nY6tRKPSeLjHQhIjYlaLMr8rM65rOUzoceHVEPESxPHVkRHy62UhA8ee4OTO3/hazlqLgm3Y08J+Z+f3M/B/gOuCwhjPN9b2IWAJQbofyq3oVEbEaOA74k2zH+c+/SfGD+evl3//9gLsj4nmNpipsBq7Lwr9S/PZc+xu2o1DorbzEQPnT9VJgU2Z+tOk8W2XmWZm5X2Yuo3it/jkzG59xZuZ3gW9HxPJy11HANxqMtNXDwKERsXv5Z3oULXizdo7PAavLr1cD1zeY5eciYhVwBvDqzPxx03kAMvPezNwnM5eVf/83Ay8r/+417R+BIwEi4iDg1xnCRcRaX+jlGy9bLzGwCbi2JZcYOBw4mWIGvKEcxzYdquXeDVwVEfcAK4APNZyH8jeGtcDdwL0U/yYa+bRhRFwN3Aksj4jNEXEKcB5wTEQ8QHHmxnktyfW3wLOAm8u/+3/fklyN20Guy4AXlKcyTgOrh/FbjZ8UlaSOaP0MXZLUGwtdkjrCQpekjrDQJakjLHRJ6ggLXZI6wkKXpI6w0CWpI/4fDfqWVVBB9loAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sampler2 = SampleLIDAR(theta2)\n",
    "\n",
    "samples2 = sampler2.sample(ground_truth)\n",
    "plt.hist(samples2, bins=1000, density=True);\n",
    "plt.axvline(TEST_GROUND_TRUTH, c='r', linestyle='--')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}