baby weight bayes 1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WHO data url\n",
    "https://www.who.int/childgrowth/standards/weight_for_age/en/\n",
    "#### WHO weight for age txt link\n",
    "https://www.who.int/childgrowth/standards/wfa_girls_z_exp.txt\n",
    "<br>https://www.who.int/childgrowth/standards/wfa_girls_p_exp.txt\n",
    "\n",
    "##### Ages 5-19 data\n",
    "https://www.who.int/growthref/who2007_weight_for_age/en/\n",
    "<br>https://www.who.int/growthref/wfa_girls_z_WHO2007_exp.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "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>Age</th>\n",
       "      <th>L</th>\n",
       "      <th>M</th>\n",
       "      <th>S</th>\n",
       "      <th>P01</th>\n",
       "      <th>P1</th>\n",
       "      <th>P3</th>\n",
       "      <th>P5</th>\n",
       "      <th>P10</th>\n",
       "      <th>P15</th>\n",
       "      <th>P25</th>\n",
       "      <th>P50</th>\n",
       "      <th>P75</th>\n",
       "      <th>P85</th>\n",
       "      <th>P90</th>\n",
       "      <th>P95</th>\n",
       "      <th>P97</th>\n",
       "      <th>P99</th>\n",
       "      <th>P999</th>\n",
       "      <th>day</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.3809</td>\n",
       "      <td>3.2322</td>\n",
       "      <td>0.14171</td>\n",
       "      <td>2.002</td>\n",
       "      <td>2.273</td>\n",
       "      <td>2.440</td>\n",
       "      <td>2.532</td>\n",
       "      <td>2.678</td>\n",
       "      <td>2.779</td>\n",
       "      <td>2.932</td>\n",
       "      <td>3.232</td>\n",
       "      <td>3.550</td>\n",
       "      <td>3.729</td>\n",
       "      <td>3.853</td>\n",
       "      <td>4.041</td>\n",
       "      <td>4.166</td>\n",
       "      <td>4.409</td>\n",
       "      <td>4.846</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>0.3259</td>\n",
       "      <td>3.1957</td>\n",
       "      <td>0.14578</td>\n",
       "      <td>1.963</td>\n",
       "      <td>2.231</td>\n",
       "      <td>2.398</td>\n",
       "      <td>2.490</td>\n",
       "      <td>2.635</td>\n",
       "      <td>2.737</td>\n",
       "      <td>2.892</td>\n",
       "      <td>3.196</td>\n",
       "      <td>3.520</td>\n",
       "      <td>3.704</td>\n",
       "      <td>3.831</td>\n",
       "      <td>4.026</td>\n",
       "      <td>4.155</td>\n",
       "      <td>4.408</td>\n",
       "      <td>4.865</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2</td>\n",
       "      <td>0.3101</td>\n",
       "      <td>3.2104</td>\n",
       "      <td>0.14637</td>\n",
       "      <td>1.972</td>\n",
       "      <td>2.240</td>\n",
       "      <td>2.408</td>\n",
       "      <td>2.500</td>\n",
       "      <td>2.646</td>\n",
       "      <td>2.748</td>\n",
       "      <td>2.904</td>\n",
       "      <td>3.210</td>\n",
       "      <td>3.538</td>\n",
       "      <td>3.723</td>\n",
       "      <td>3.853</td>\n",
       "      <td>4.049</td>\n",
       "      <td>4.181</td>\n",
       "      <td>4.438</td>\n",
       "      <td>4.902</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>3</td>\n",
       "      <td>0.2986</td>\n",
       "      <td>3.2315</td>\n",
       "      <td>0.14657</td>\n",
       "      <td>1.986</td>\n",
       "      <td>2.255</td>\n",
       "      <td>2.424</td>\n",
       "      <td>2.516</td>\n",
       "      <td>2.663</td>\n",
       "      <td>2.766</td>\n",
       "      <td>2.923</td>\n",
       "      <td>3.232</td>\n",
       "      <td>3.562</td>\n",
       "      <td>3.749</td>\n",
       "      <td>3.879</td>\n",
       "      <td>4.079</td>\n",
       "      <td>4.212</td>\n",
       "      <td>4.471</td>\n",
       "      <td>4.942</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4</td>\n",
       "      <td>0.2891</td>\n",
       "      <td>3.2558</td>\n",
       "      <td>0.14658</td>\n",
       "      <td>2.004</td>\n",
       "      <td>2.274</td>\n",
       "      <td>2.443</td>\n",
       "      <td>2.536</td>\n",
       "      <td>2.684</td>\n",
       "      <td>2.787</td>\n",
       "      <td>2.945</td>\n",
       "      <td>3.256</td>\n",
       "      <td>3.589</td>\n",
       "      <td>3.778</td>\n",
       "      <td>3.909</td>\n",
       "      <td>4.110</td>\n",
       "      <td>4.245</td>\n",
       "      <td>4.507</td>\n",
       "      <td>4.983</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Age       L       M        S    P01     P1     P3     P5    P10    P15  \\\n",
       "0    0  0.3809  3.2322  0.14171  2.002  2.273  2.440  2.532  2.678  2.779   \n",
       "1    1  0.3259  3.1957  0.14578  1.963  2.231  2.398  2.490  2.635  2.737   \n",
       "2    2  0.3101  3.2104  0.14637  1.972  2.240  2.408  2.500  2.646  2.748   \n",
       "3    3  0.2986  3.2315  0.14657  1.986  2.255  2.424  2.516  2.663  2.766   \n",
       "4    4  0.2891  3.2558  0.14658  2.004  2.274  2.443  2.536  2.684  2.787   \n",
       "\n",
       "     P25    P50    P75    P85    P90    P95    P97    P99   P999  day  \n",
       "0  2.932  3.232  3.550  3.729  3.853  4.041  4.166  4.409  4.846    0  \n",
       "1  2.892  3.196  3.520  3.704  3.831  4.026  4.155  4.408  4.865    1  \n",
       "2  2.904  3.210  3.538  3.723  3.853  4.049  4.181  4.438  4.902    2  \n",
       "3  2.923  3.232  3.562  3.749  3.879  4.079  4.212  4.471  4.942    3  \n",
       "4  2.945  3.256  3.589  3.778  3.909  4.110  4.245  4.507  4.983    4  "
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_girls_percentiles = pd.read_csv('https://www.who.int/childgrowth/standards/wfa_girls_p_exp.txt', sep=\"\\t\") #, header=None)\n",
    "data_girls_percentiles['day'] = data_girls_percentiles['Age']\n",
    "data_girls_percentiles.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_day = np.max(data_girls_percentiles['day'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### LMS z score\n",
    "https://www.cdc.gov/nchs/data/nhsr/nhsr063.pdf\n",
    "<br>https://www.indianpediatrics.net/jan2014/jan-37-43.htm\n",
    "<br> Z score and x-coordinate are same thing\n",
    "#### scipy.stats.norm\n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.039999999999999994\n",
      "0.39862325420460504\n"
     ]
    }
   ],
   "source": [
    "def LMS_z_score(weight_kg, L, M, S):\n",
    "    if L != 0:\n",
    "        return 1/(S * L) * ((weight_kg / M) ** L - 1)\n",
    "    elif L == 0:\n",
    "        return ln(weight_kg / M)/S\n",
    "print(LMS_z_score(weight_kg = 10, L = -1, M = 8, S = 5))\n",
    "\n",
    "def LMS_pdf(weight_kg, L, M, S):\n",
    "    x = LMS_z_score(weight_kg = weight_kg, L = L, M = M, S = S)\n",
    "    return scipy.stats.norm.pdf(x = x, loc=0, scale=1)\n",
    "\n",
    "print(LMS_pdf(weight_kg = 10, L = -1, M = 8, S = 5))\n",
    "\n",
    "def kg_to_lbs(kg):\n",
    "    return kg * 2.2046\n",
    "\n",
    "def lbs_to_kg(lb):\n",
    "    return lb / 2.2046"
   ]
  },
  {
   "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>day</th>\n",
       "      <th>lbs</th>\n",
       "      <th>oz</th>\n",
       "      <th>lbs_2</th>\n",
       "      <th>kg</th>\n",
       "      <th>L</th>\n",
       "      <th>M</th>\n",
       "      <th>S</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "      <td>9.0</td>\n",
       "      <td>7.5625</td>\n",
       "      <td>3.430327</td>\n",
       "      <td>0.3809</td>\n",
       "      <td>3.2322</td>\n",
       "      <td>0.14171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>15.0</td>\n",
       "      <td>6.9375</td>\n",
       "      <td>3.146829</td>\n",
       "      <td>0.2986</td>\n",
       "      <td>3.2315</td>\n",
       "      <td>0.14657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>16</td>\n",
       "      <td>8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8.0000</td>\n",
       "      <td>3.628776</td>\n",
       "      <td>0.2218</td>\n",
       "      <td>3.6438</td>\n",
       "      <td>0.14258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>65</td>\n",
       "      <td>11</td>\n",
       "      <td>6.4</td>\n",
       "      <td>11.4000</td>\n",
       "      <td>5.171006</td>\n",
       "      <td>0.0877</td>\n",
       "      <td>5.2360</td>\n",
       "      <td>0.12934</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   day  lbs    oz    lbs_2        kg       L       M        S\n",
       "0    0    7   9.0   7.5625  3.430327  0.3809  3.2322  0.14171\n",
       "1    3    6  15.0   6.9375  3.146829  0.2986  3.2315  0.14657\n",
       "2   16    8   0.0   8.0000  3.628776  0.2218  3.6438  0.14258\n",
       "3   65   11   6.4  11.4000  5.171006  0.0877  5.2360  0.12934"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "baby_weights = pd.read_csv('baby weights.csv')\n",
    "baby_weights['lbs_2'] = baby_weights['lbs'] + baby_weights['oz'] / 16\n",
    "baby_weights['kg'] = lbs_to_kg(lb = baby_weights['lbs_2'])\n",
    "baby_weights = baby_weights.join(data_girls_percentiles[['day', 'L', 'M', 'S']].set_index('day'), on = 'day')\n",
    "\n",
    "baby_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean z_score of baby's weight: 0.388 sigma.\n",
      "Standard deviation of baby's weight: 0.014\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>day</th>\n",
       "      <th>lbs</th>\n",
       "      <th>oz</th>\n",
       "      <th>lbs_2</th>\n",
       "      <th>kg</th>\n",
       "      <th>L</th>\n",
       "      <th>M</th>\n",
       "      <th>S</th>\n",
       "      <th>z_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "      <td>9.0</td>\n",
       "      <td>7.5625</td>\n",
       "      <td>3.430327</td>\n",
       "      <td>0.3809</td>\n",
       "      <td>3.2322</td>\n",
       "      <td>0.14171</td>\n",
       "      <td>0.364552</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>15.0</td>\n",
       "      <td>6.9375</td>\n",
       "      <td>3.146829</td>\n",
       "      <td>0.2986</td>\n",
       "      <td>3.2315</td>\n",
       "      <td>0.14657</td>\n",
       "      <td>0.392501</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>16</td>\n",
       "      <td>8</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8.0000</td>\n",
       "      <td>3.628776</td>\n",
       "      <td>0.2218</td>\n",
       "      <td>3.6438</td>\n",
       "      <td>0.14258</td>\n",
       "      <td>0.398775</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>65</td>\n",
       "      <td>11</td>\n",
       "      <td>6.4</td>\n",
       "      <td>11.4000</td>\n",
       "      <td>5.171006</td>\n",
       "      <td>0.0877</td>\n",
       "      <td>5.2360</td>\n",
       "      <td>0.12934</td>\n",
       "      <td>0.397088</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   day  lbs    oz    lbs_2        kg       L       M        S   z_score\n",
       "0    0    7   9.0   7.5625  3.430327  0.3809  3.2322  0.14171  0.364552\n",
       "1    3    6  15.0   6.9375  3.146829  0.2986  3.2315  0.14657  0.392501\n",
       "2   16    8   0.0   8.0000  3.628776  0.2218  3.6438  0.14258  0.398775\n",
       "3   65   11   6.4  11.4000  5.171006  0.0877  5.2360  0.12934  0.397088"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "baby_weights['z_score'] = baby_weights[['kg', 'L', 'M', 'S']].apply(axis = 1,\n",
    "                                                                    func = lambda in_row: \n",
    "                                                                    LMS_pdf(weight_kg = in_row['kg'],\n",
    "                                                                            L = in_row['L'],\n",
    "                                                                            M = in_row['M'], \n",
    "                                                                            S = in_row['S']))\n",
    "likelihood_mean = np.mean(baby_weights['z_score'])\n",
    "print(\"Mean z_score of baby's weight: \" + \n",
    "      str(np.round_(likelihood_mean, decimals = 3)) +\n",
    "     \" sigma.\")\n",
    "\n",
    "likelihood_stdev = np.std(baby_weights['z_score'])\n",
    "print(\"Standard deviation of baby's weight: \" + \n",
    "      str(np.round_(likelihood_stdev, decimals = 3)))\n",
    "\n",
    "baby_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max posterior: 0.13423535505552397\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>lbs</th>\n",
       "      <th>day</th>\n",
       "      <th>prior_z_score</th>\n",
       "      <th>likelihood_z_score</th>\n",
       "      <th>prior_pdf</th>\n",
       "      <th>likelihood_pdf</th>\n",
       "      <th>posterior_pdf</th>\n",
       "      <th>prior_cdf</th>\n",
       "      <th>likelihood_cdf</th>\n",
       "      <th>posterior_cdf</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>1856</td>\n",
       "      <td>-51.272834</td>\n",
       "      <td>-3726.946392</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.099099</td>\n",
       "      <td>1856</td>\n",
       "      <td>-48.966315</td>\n",
       "      <td>-3560.548885</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.198198</td>\n",
       "      <td>1856</td>\n",
       "      <td>-46.925274</td>\n",
       "      <td>-3413.303522</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.297297</td>\n",
       "      <td>1856</td>\n",
       "      <td>-45.100721</td>\n",
       "      <td>-3281.676138</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.396396</td>\n",
       "      <td>1856</td>\n",
       "      <td>-43.455616</td>\n",
       "      <td>-3162.994549</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        lbs   day  prior_z_score  likelihood_z_score  prior_pdf  \\\n",
       "0  1.000000  1856     -51.272834        -3726.946392        0.0   \n",
       "1  1.099099  1856     -48.966315        -3560.548885        0.0   \n",
       "2  1.198198  1856     -46.925274        -3413.303522        0.0   \n",
       "3  1.297297  1856     -45.100721        -3281.676138        0.0   \n",
       "4  1.396396  1856     -43.455616        -3162.994549        0.0   \n",
       "\n",
       "   likelihood_pdf  posterior_pdf  prior_cdf  likelihood_cdf  posterior_cdf  \n",
       "0             0.0            0.0        0.0             0.0            0.0  \n",
       "1             0.0            0.0        0.0             0.0            0.0  \n",
       "2             0.0            0.0        0.0             0.0            0.0  \n",
       "3             0.0            0.0        0.0             0.0            0.0  \n",
       "4             0.0            0.0        0.0             0.0            0.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "year_5_pdf = pd.DataFrame({'lbs': np.linspace(1, 100, num = 1000)})\n",
    "year_5_pdf['day'] = max_day\n",
    "year_5_pdf = year_5_pdf.join(data_girls_percentiles[['day', 'L', 'M', 'S']].set_index('day'), on = 'day')\n",
    "year_5_pdf['kg'] = lbs_to_kg(lb = year_5_pdf['lbs'])\n",
    "year_5_pdf['prior_z_score'] = year_5_pdf[['kg', 'L', 'M', 'S']].apply(axis = 1,\n",
    "                                                                    func = lambda in_row: \n",
    "                                                                    LMS_z_score(weight_kg = in_row['kg'],\n",
    "                                                                                L = in_row['L'],\n",
    "                                                                                M = in_row['M'], \n",
    "                                                                                S = in_row['S']))\n",
    "for colname in ['L', 'M', 'S', 'kg']:\n",
    "    del year_5_pdf[colname]\n",
    "\n",
    "year_5_pdf['likelihood_z_score'] = (year_5_pdf['prior_z_score'] - likelihood_mean) / likelihood_stdev\n",
    "year_5_pdf['prior_pdf'] = scipy.stats.norm.pdf(x = year_5_pdf['prior_z_score'], loc=0, scale=1)\n",
    "year_5_pdf['likelihood_pdf'] = scipy.stats.norm.pdf(x = year_5_pdf['likelihood_z_score'], loc=0, scale=1)\n",
    "year_5_pdf['posterior_pdf'] = year_5_pdf['prior_pdf'] * year_5_pdf['likelihood_pdf']\n",
    "year_5_pdf['prior_cdf'] = year_5_pdf['prior_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['likelihood_cdf'] = year_5_pdf['likelihood_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['posterior_cdf'] = year_5_pdf['posterior_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['prior_cdf'] = year_5_pdf['prior_cdf'] / np.max(year_5_pdf['prior_cdf'])\n",
    "year_5_pdf['likelihood_cdf'] = year_5_pdf['likelihood_cdf'] / np.max(year_5_pdf['likelihood_cdf'])\n",
    "year_5_pdf['posterior_cdf'] = year_5_pdf['posterior_cdf'] / np.max(year_5_pdf['posterior_cdf'])\n",
    "\n",
    "### check max posterior\n",
    "print(\"Max posterior: \"+ str(np.max(year_5_pdf['posterior_pdf'])))\n",
    "year_5_pdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.01801801801802"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Location of max posterior\n",
    "max_posterior_lbs = (year_5_pdf\n",
    "                     .loc[year_5_pdf['posterior_pdf'] == np.max(year_5_pdf['posterior_pdf']), 'lbs']\n",
    "                     .iloc[0])\n",
    "max_posterior_lbs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Helpful intro to matplotlib\n",
    "https://bespokeblog.wordpress.com/2011/07/07/basic-data-plotting-with-matplotlib-part-2-lines-points-formatting/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcFNW5//HPMzNssquYqKhgAgIGRESQTTEqonIxLlx3XIJeVIx440+90WiM3uRGE0WuCxIV1EvEqEgw4BIQBhdQBlFQEEII4ijqiIgQ1pl5fn9UzdAM09OzdE33dH/fr1e/pmvpqqd6uuvpc07VOebuiIiIAOSkOgAREUkfSgoiIlJOSUFERMopKYiISDklBRERKaekICIi5ZQUJFJmNsjMVqY6jorMrIOZuZnlpTqWWDV5v8xssJkVRh1TdZjZRWb2WqrjkLpTUsgwZrbWzLaZ2RYz22hmM83skFTF4+5vuPsRUWw7HY7VzA4Mk8v3YubdGmfeK4m2l8z3y8wmm9ndCdZxM/tX+B6WPW5K8Jq9Eqq7T3H3IcmIu5L9zTOzUVFsW/ampJCZ/s3dWwAHAl8C/5vieKKU0mN19/XAauD4mNnHAx9XMm9+PYZWE0e5e4uYxz2pDkhSR0khg7n7duB5oFvZPDM7w8yWmNl3Zvapmf0qZtlMM7sudhtmttTMfhI+72JmfzOzb8xspZn9e8x6p5vZcjPbbGafmdmN4fw9qjjM7BYz+0e43nIzOytm2WVm9qaZ/T785f9PMzstimONcYWZfW5m683s5+Hrvm9mW81sv5htHWNmRWbWqJJtzCdMAGaWCxwNPFBhXr9wPcysSXiM68zsSzObYGbN4rxfvcJj2Gxmz5nZsxV//ZvZz83sq/AYLg/nXQVcBNwU/vp/qTrvY4Xt9jGzgvD9+9LM7os5XoBvw233K/vfxbzWzewaM/t7GPtdZvYDM1sQbu/PZtY4XLetmf01fH83hs/bh8v+GxgEPBju68FwftzPotSRu+uRQQ9gLXBy+Hwf4EngqZjlg4HuBD8IehD8uv5JuOzfgXdi1j0K2AA0BpoDnwKXA3lAL+Br4Mhw3fXAoPB5W6BXzP4KY7Y5Ajgo3P95wL+AA8NllwG7gCuBXOBq4HPAIjjWDoADz4TH1h0oitneLODqmG3dD/xvnDguBT4In/cmOGl2qjBvG9A4nB4HzAD2BVoCLwG/rfh+he/7J8D1QCPgbGAncHfMusXAr8PlpwNbgbbh8sll61bxeXHgh3GWLQAuCZ+3AI6r8N7lxax7GfBmhe3OAFoBRwI7gDnA4UBrYDlwabjufsA54f+wJfAcMD1mW/OAUTHTVX4W9ajjOSTVAeiR5H9ocKLcAnwbnjA+B7pXsf444P7weRPgG6BTOP174OHw+XnAGxVe+yhwR/h8HfAfQKsK65Sf5OLs/33gzPD5ZcDqmGX7hCeX70dwrGUnti4xy+8BHo853rfC57nAF0CfONvtAJQQJMMbgP8O538WM29uOM8IEuEPYl7fD/hnxfeLoKTxGTFJEXiTPZPCNvY8OX/F7pP3ZKqXFL4L38Oyx6nhsvnAncD+lRxvdZLCgJjpxcDNMdN/AMbFiaknsDFmeh57JoUqP4t61O2h6qPM9BN3b0Nwkh8D5JvZ9wHMrK+ZzQ2L6puA0cD+AO6+A/gzcLGZ5QAXAE+H2zwM6Gtm35Y9CKonvh8uP4fgl+onZpZvZv0qC8zMRprZ+zHb+FHZ/kNflD1x963h0xbJPtYYn8Y8/4SgFAPwF6CbmR0OnAJscvd3KwvA3dcChcBAghP5G+GiBTHzyqpc2hEku8Ux78Er4fyKDgI+8/CsV0m8ABvcvThmeitVv1+V6eXubWIer4bzfwp0Bj42s0VmNqyG2/0y5vm2SqZbAJjZPmb2qJl9YmbfEbxXbcJqt8ok+ixKHSgpZDB3L3H3aQS/YgeGs/9EUKw/xN1bAxMIfr2WeZLgC3YSsNXdF4TzPwXyK5w8Wrj71eG+Frn7mcABwHSC5LIHMzsM+CPByXu/8GT+YYX91+exAsRerXQoQWkDD9oo/kzwXlzC7uQYzxsEJ/9+wNsV5g1kd1L4muCEeGTM+9jag8byitYDB5tZbMw1ubqqTl0gu/vf3f0Cgv/p74Dnzax5XbdbiZ8DRwB93b0Vuxvoy4674v6q/CxK3SgpZDALnElQhbEinN0S+Mbdt5tZH+DC2NeESaCUoHgfeyL8K9DZzC4xs0bh41gz62pmjS24Tr21u+8iqI4oqSSkshNKURjf5QQlhZQca+iX4S/VIwnqqJ+NWfYUQbXIcOD/EoQwHxgJfO7u34Xz3gzntSYoNeDupQSJ8X4zOyCM/WAzO7WSbS4geB/HmFleeHx9EsQR60uCOvxaMbOLzaxdGPO34ewSgv9faV22XUFLgkT5rZntC9xRYXnF44j7WUxSPFlNSSEzvWRmWwhOzv9N0KD3UbjsGuDXZrYZuJ1KftETnAy7E3MidPfNwBDgfIJf018Q/HpsEq5yCbA2LP6PBi6uuFF3X06QbBYQfNG7A2/V6Ujrfqz5BJeUzgF+7+7lN2C5+1sEJ7/3wiqiquQT/KJ+M2be+0AzYHFMVRjAzeE+F4bv12yCX8p7cPedBI3LPyU4KV9McELckSCWMo8TVIF9a2bTq1jvA9vzPoVx4fyhwEfh+/sAcL67bw+P5b+Bt8JtH1fNeOIZR/A+fQ0sJKhOi/UAcG54ZdL4anwWpQ5sz+pKkaDeH7jK3QcmXDnDmdnrwJ/c/bFUxwJgZu8AE9x9UqpjkcykkoLswcz2IfiFPTHVsaSamR1LcLnjs4nWjTCGEyy4byLPzC4luLQ24Z3RIrWlpCDlwnrtIoKqnT+lOJyUMrMnCap1xobVFalyBPABsImgQfZcD+6iFomEqo9ERKScSgoiIlIurboNro7999/fO3TokOowREQalMWLF3/t7pXdJLmHBpcUOnToQEFBQarDEBFpUMzsk+qsp+ojEREpp6QgIiLllBRERKRcg2tTqMyuXbsoLCxk+/btqQ5FKtG0aVPat29Po0aVjU8jIukkI5JCYWEhLVu2pEOHDuzZoaSkmruzYcMGCgsL6dixY6rDEZEEIqs+MrMnLBgi8MM4y83MxpvZaguGfOxV231t376d/fbbTwkhDZkZ++23n0pxIg1ElG0Kkwl6WYznNIIhCzsBVwGP1GVnSgjpS/8bkYYjsuojd59vZh2qWOVMgvF0naAL4TZmdqD6dZEoLVm9hFufvZUtO7ekOpS0VOpGaWkOpSU5lJbmBs/Dh7vV6BHLy8bLCUck9d0LKBtLp+K8xOvUVt1+pKSyY6CTO/XliVuvjXQfqWxTOJg9hxYsDOftlRTM7CqC0gSHHnpovQQnmWnUE6N4r8l7SRjrraFKcOBGUH+QEa2NmefVVSVA5iaFyj6dlSZhd59I2JVz7969M7YHv+nTp9O5c2e6detWo9fNmDGD5cuXc8stt0QU2W6TJ0+moKCABx98kKKiIoYNG8bOnTsZP348gwYNinz/dfXF9i9otKMRO/+wM9WhJMXOnfDhh7BkCaxatfuxenWwLJYZtGkDbdsGf2Oft2oF++wDzZrtfjRtuufzxo0hLw9yc4O/ZY+K02XzcnKCfZbtO9GjOutVXKe2VKMZXyqTQiF7jjfbnnB83Gw1ffp0hg0bVqOkUFxczPDhwxk+fHiNXpOXV/d//Zw5c+jSpQtPPvlknbdVXzaVbqKltUx1GLX27bcwbx7MmQMLF8LSpbtP/o0bww9/CJ07wxlnQIcOcNBBux/f+x7oqmBJJJVJYQbB2LNTgb7ApmS0J4wdO5b333+/zsHF6tmzJ+PGjatynbVr1zJ06FD69u3LkiVL6Ny5M0899RQLFizgxhtvpLi4mGOPPZZHHnmEJk2acMsttzBjxgzy8vIYMmQIZ599NjNmzCA/P5+7776bF154AYBrr72WoqIi9tlnH/74xz/SpUsXLrvsMvbdd1+WLFlCr1696N69e/mv908++YQrrriCoqIi2rVrx6RJkzj00EP3es0f/vCHvY5hy5YtXHfddRQUFGBm3HHHHZxzzjlMmjSJ3/72txx44IF07tyZJk2a8P7773PTTTexbds2evbsyYIFC2jWrFlS3/co7MzdyX7sl+owamTNGnjuOZg2DQoKoLQUmjeHPn3g+uuhd2/o1Qs6dgx+oYvURWRJwcyeAQYD+5tZIcFg3I0A3H0CMAs4nWCs2q0Eg6Y3aCtXruTxxx9nwIABXHHFFdx33308+uijzJkzh86dOzNy5EgeeeQRRo4cyYsvvsjHH3+MmfHtt9/Spk0bhg8fzrBhwzj33HMBOOmkk5gwYQKdOnXinXfe4ZprruH1118HYNWqVcyePZvc3FwmT55cHsOYMWMYOXIkl156KU888QQ/+9nPmD59+l6vqcxdd91F69atWbZsGQAbN25k/fr13HHHHSxevJjWrVtz4okncvTRR9OzZ09+/etflyejhsLNyfX0P3Nu2QJTpsBjjwWJAODYY+GXv4STToK+fYOSgUiyRXn10QUJljsRtJgk+kUfpUMOOYQBAwYAcPHFF3PXXXfRsWNHOnfuDMCll17KQw89xJgxY2jatCmjRo3ijDPOYNiwYXtta8uWLbz99tuMGDGifN6OHbvHax8xYkSlJ/cFCxYwbdo0AC655BJuuummhK8pM3v2bKZOnVo+3bZtW6ZPn87gwYNp1y7ocfe8885j1apV1Xo/0pGbk5PGvbt89hncey9MmgTffQc9esA998CIEUF1kEjUdI1BElX3evy8vDzeffdd5syZw9SpU3nwwQfLSwBlSktLadOmTdyqsObNm9c4pkSvcfdKjyHT7jPIsfRLCl98Af/zPzBhApSUwL//O1x7LfTrp0ZRqV/p9+1owNatW8eCBQsAeOaZZzj55JNZu3Ytq1evBuDpp5/mhBNOYMuWLWzatInTTz+dcePGlZ/4W7ZsyebNwXDArVq1omPHjjz33HNAcML+4IMPEsbQv3//8l/7U6ZMYeDAgdWOf8iQIXtUBW3cuJG+ffsyb948NmzYwK5du8rjaajcPK2SQnExjBsXNA4/+CBcdFFw9dCUKdC/vxKC1L/0+XZkgK5du/Lkk0/So0cPvvnmG2644QYmTZrEiBEj6N69Ozk5OYwePZrNmzczbNgwevTowQknnMD9998PwPnnn8+9997L0UcfzT/+8Q+mTJnC448/zlFHHcWRRx7JX/7yl4QxjB8/nkmTJtGjRw+efvppHnjggWrHf9ttt7Fx40Z+9KMfcdRRRzF37lwOPPBAfvWrX9GvXz9OPvlkevWqdW8kaSGdqo+WLg0aiW+4AQYMgBUr4PHHgwZjkVSxoGq/4ejdu7dXHHltxYoVdO3aNUURBdauXcuwYcP48MNKu3rKeunwPwLIuSWHI/1Ilv1uWcpicIfx4+Gmm2DffeGhh+Css1QqkGiZ2WJ3751oPbUpSHYxUlpS2LIFLrkEpk+HYcPgiSegXcJRc0Xqj5JCknTo0KHBlBImTZq0V7XSgAEDeOihh1IUUf1xc3JJzSWphYXwb/8WVBvddx+MHavSgaQfJYUsdPnll3P55Q3+tpDasdRcfbRiBZx8MmzeDC+9BKefXu8hiFSLkoJklxzq/ea1Dz8MbjgzgzffDO49EElXSgqSXYx6rT5asQJOPDHoc+j116FLl3rbtUitKClIdqnH6qPPP4ehQ4PeQvPzoVOnetmtSJ0oKUjWKC0trbfqo+++g9NOg2++CXo1VUKQhiI97uLJIrfffjuzZ8+ul31ddtllPP/88wC88cYbHHnkkfTs2ZNt27bVy/7Tza7iXQDk5kSbFEpLYeRI+OgjeP55OOaYSHcnklQqKdSjkpISfv3rX9f4NVV1YlddU6ZM4cYbb8zeq46AXSVhUrBok8I998Bf/gL33w+nnhrprkSSLuOSwtixkOThFOjZM+ifpirxxlPo1q0bV1xxBa+99hpjxozhlVdeKe8ee86cOZWOtdChQ4c9XnP++efvtb/Vq1czevRoioqKyM3N5bnnnuPwww/nuuuu4/XXX6djx46U3a3+2GOP8ec//5lXX32V2bNnM2XKlOS+QQ3Ezl3BaDRRJoXXX4dbb4XzzgvGOhBpaFR9lEQrV67kqquuYunSpbRq1YqHH34YgKZNm/Lmm2/ucXLfvn07l112Gc8++yzLli2juLiYRx55pHx5Za+JddFFF3HttdfywQcf8Pbbb3PggQfy4osvsnLlSpYtW8Yf//hH3n77bQBGjRrF8OHDuffee7M2IQDsLA6TQkTVRxs3BtVGnTsH4yDoxjRpiDKupJDC4RT2Gk9h/PjxQDAGQUUrV66sdKyFsWPHxn1Nmc2bN/PZZ59x1llnAUECAZg/fz4XXHABubm5HHTQQfz4xz9O3sFlgLKkkJcTzcf+uuvgyy+DLixatIhkFyKRU0khiSqOO1A2Xdk4Bok6Iqxq7IOqXptpYx8kU3FxMRBNSeH554Purn/5y6DnU5GGSkkhiSqOp1DVWAZdunSpdKyF6mjVqhXt27cvH2Zzx44dbN26leOPP56pU6dSUlLC+vXrmTt3bh2PKLOUVx8luU1h06aglHDMMfCLXyR10yL1TkkhiSqOp3D11VfHXbdp06aVjrVQXU8//TTjx4+nR48e9O/fny+++IKzzjqLTp060b17d66++upqJ5lsEdUlqbfdBl99BY8+CnkZVyEr2UYf4STKyclhwoQJe8xbu3btHtOTJ08uf37SSSexZMmSvbZT8TWV6dSp015DeAJ7jJwWb7/ZKoo2hcWL4eGH4ZprdD+CZAaVFCRrlJcUknDfBwSD5YwZE4yHcPfdSdmkSMqppJAkUY2ncO211/LWW2/tMe/666/P6pvQaqvs5rVklRReeAEWLgyG0GzdOimbFEk5JYU0lw0D39SXZCaFnTvhv/4LjjwSLr20zpsTSRtKCpI1ktnQPHEirF4NM2dCkmqjRNKC2hQka5SUlgB1vyR161a46y4YPDjoCVUkk6ikIFmjLCnk5NTtt9DEicElqC+8oK4sJPOopCBZo+xO8Lrc9b19e9AL6uDBUMW9iSINlpJCkrQIO7v5/PPPOffcc4Hg3oAxY8ZUexux60+YMIGnnnoKgMGDB1NQUJDkiJOz3Q4dOvD1118DMH78eLp27cpFF12UjPCSrtRLATBqnxQefxzWrw+6sxDJRKo+SrKDDjqofGCbuqjJ3c3p4uGHH+bll1+mY8eOqQ6lUqWlQVKo7XCcu3bB734H/fsH4y6LZKKMSwpjXxnL+18kd0CFnt/vybih1et+de3atQwbNmyvexZmzpzJ3XffzUsvvYS7M3r0aNatWwfAuHHjyntXLfOrX/2KFi1acOONNwLw3HPPcc011/Dtt9/y+OOPM2jQILZv387VV19NQUEBeXl53HfffZx44olx52/bto3LL7+c5cuX07Vr14QjsL3yyiv84he/oKSkhP333585c+awYcMGLrjgAoqKiujTp095lczo0aNZs2YNw4cP54orruCGG26o1vtVn+paffT88/Dpp8EdzGpLkEyVcUkhHb344ovcd999zJo1i7Zt23LhhRdyww03MHDgQNatW8epp57KihUrqtxGcXEx7777LrNmzeLOO+9k9uzZ5fcwLFu2jI8//pghQ4awatWquPMfeeQR9tlnH5YuXcrSpUvp1atX3P0VFRVx5ZVXMn/+fDp27Mg333wDwJ133snAgQO5/fbbmTlzJhMnTgSC6q5XXnmFuXPnsv/++yfjbUu6suqj2pQU3IOR1Dp1gtNPT3ZkIukj0qRgZkOBB4Bc4DF3/58Ky1sD/wccGsbye3efVJd9VvcXfX2ZO3cuBQUFvPbaa7Rq1QqA2bNns3z58vJ1vvvuOzZv3lzlds4++2wAjjnmmPK+kd58802uu+46IOh19bDDDmPVqlVx58+fP5+f/exnAPTo0YMePXrE3d/ChQs5/vjjy6uC9t13XyAYs2HatGkAnHHGGbRt27ZG70cqlbcp1OJn/sKFsGgRPPgg1PHiJZG0FllSMLNc4CHgFKAQWGRmM9x9ecxq1wLL3f3fzKwdsNLMprj7zqjiqm+HH344a9asYdWqVfQOO9ovLS1lwYIFNGvWrNrbadKkCRD021M2LkC8cRWSMd6Cu8ddt6GO2VCXksK4cUFXFrp7WTJdlL95+gCr3X1NeJKfCpxZYR0HWlpwlmkBfAMURxhTvTvssMOYNm0aI0eO5KOPPgJgyJAhe/Rm+n4tB5U+/vjjy4fXXLVqFevWreOII46o1vwPP/yQpUuXxt12v379yM/P55///CdAefVR7DZefvllNm7cWKvYU6GsobmmSW3duuCehCuv1IhqkvmiTAoHA5/GTBeG82I9CHQFPgeWAde7hz/nYpjZVWZWYGYFRUVFUcUbmSOOOIIpU6YwYsQI/vGPfzB+/HgKCgro0aMH3bp126u77eq65pprKCkpoXv37px33nlMnjyZJk2axJ1/9dVXs2XLFnr06ME999xDnz594m67Xbt2TJw4kbPPPpujjjqqfHjQO+64g/nz59OrVy9ee+01Dj300FrFngq1bWh+9NHdPaKKZDpLNCxkrTdsNgI41d1HhdOXAH3c/bqYdc4FBgD/CfwA+BtwlLt/F2+7vXv39orX1q9YsYKuXbsm/yAkadLhf/TU7Ke49K1Lue3w27jrkruq9ZriYjj00GCshJdeijhAkQiZ2WJ3TzhYbJQlhULgkJjp9gQlgliXA9M8sBr4J9Alwpgki9WmpDBrVnCz2qhRUUUlkl6ivPpoEdDJzDoCnwHnAxdWWGcdcBLwhpl9DzgCWBNhTFKJvn37smPHjj3mPf3003Tv3j1FEUWjNg3Njz0G3/++LkOV7BFZUnD3YjMbA7xKcEnqE+7+kZmNDpdPAO4CJpvZMsCAm93961rur8FeFZNq77zzTqTbj6qKsqZqeknqZ58FXWPffDM0ahRlZCLpI9L7FNx9FjCrwrwJMc8/B4bUdT9NmzZlw4YN7LfffkoMacbd2bBhA02bNk11KDWuPpo8GUpL4ac/jTAokTSTEXc0t2/fnsLCQhrilUnZoGnTprRv3z7VYdSo+sgdnngi6OPoBz+IOjKR9JERSaFRo0Zp2wmbpI+alBQWLIA1a+D226OOSiS96IZ9yRo1KSn86U/QtCmcdVbUUYmkFyUFyRplJYVEI6/t2gXPPgvDh0PYXZVI1lBSkKxRNhxnouqjv/0Nvv4a0nSsIJFIKSlI1igvKST42E+ZAm3bwtCh9RGVSHpRUpCsUZ37FLZsgenTYcQIaNy4viITSR9KCpI1yksKVTQ0v/QSbN2qqiPJXkoKkjXKrz6qoqF52rSgW4uBA+srKpH0oqQgWSPRJalbtwYd4J11lkZXk+ylj75kjUQ3r732WpAYwpFPRbKSkoJkjUQlhRdegH33hRNOqM+oRNKLkoJkjapKCjt3Bo3Mw4erR1TJbkoKkjWqKim8/jps2gTnnFPfUYmkFyUFyRpVlRSmTYMWLeDkk+s7KpH0oqQgWSNeSaGkJLhhbdiwoBM8kWympCBZI15JYeFCKCqCn/wkFVGJpBclBcka8e5onjkTcnPh1FNTEZVIelFSkKwRr/ror3+FQYOgTZtURCWSXpQUJGtUVn20bh0sWwZnnJGqqETSi5KCZI1S9i4pzJwZ/FVSEAkoKUjWqGzktZkz4fDDoUuXVEUlkl6UFCRrlJbuWVLYuhXmzAlKCQkGYxPJGkoKkjUqtinMnQvbtwf3J4hIQElBskbF8RRmzoTmzdUBnkgsJQXJGs7u+xTcg7ETTjoJmjRJcWAiaURJQbJG7M1rq1fDJ5/ohjWRipQUJGvE3rz22mvBvCFDUhiQSBpSUpCsUZYULMf429+gY0f4wQ9SHJRImlFSkKxRVn1UWpzD66/DKafoUlSRiiJNCmY21MxWmtlqM7slzjqDzex9M/vIzPKjjEeyW1lS+Pvfm7F5s6qORCqTF9WGzSwXeAg4BSgEFpnZDHdfHrNOG+BhYKi7rzOzA6KKR6Ss+mjJkubk5MCPf5zigETSUJQlhT7Aandf4+47ganAmRXWuRCY5u7rANz9qwjjkSxXdknqkvf24dhjoW3bFAckkoaiTAoHA5/GTBeG82J1Btqa2TwzW2xmIyvbkJldZWYFZlZQVFQUUbiS6cqqj1aubKqqI5E4okwKlTXheYXpPOAY4AzgVOCXZtZ5rxe5T3T33u7eu127dsmPVLJCWVLw0hxOOSXFwYikqcjaFAhKBofETLcHPq9kna/d/V/Av8xsPnAUsCrCuCRLlbUpNGvmHHdcioMRSVNRlhQWAZ3MrKOZNQbOB2ZUWOcvwCAzyzOzfYC+wIoIY5IsVtamcHTPHTRqlOJgRNJUZCUFdy82szHAq0Au8IS7f2Rmo8PlE9x9hZm9AiwFSoHH3P3DqGKS7LZpU2MAevfemeJIRNJXlNVHuPssYFaFeRMqTN8L3BtlHCIA69e3gv3hmF5KCiLx6I5myRrrv2gFwOGHl6Q4EpH0paQgWeOLL4OkkJebm+JIRNKXkoJkhU8+gX/9Kxg4oWw4ThHZm74dkhXy86HsNhnLUS94IvEoKUhWmDcPGjcJGphVUhCJT98OyQr5+XDAAd8CkJujNgWReJQUJOMVFsKaNbD/ARsBlRREqlLlt8PMJsc8vzTyaEQikB+O0nGAkoJIQom+HUfFPL8+ykBEopKfD61bQ6vWWwAwDbcmEleipFCxV1ORBic/HwYNgqAnFcjJUUlBJJ5E3Vy0N7PxBN1glz0v5+4/iywykSRYvx5WrYIrr4SFYdfZqj4SiS9RUvh/Mc8LogxEJArz5wd/TzgBFsxVUhBJpMqk4O5P1lcgIlGYNw9atoSjjwYvSwqqPhKJK+G3w8wuNbP3zOxf4aMg3rCZIukmPx8GDIC8vN0jr6mkIBJfoktSRwJjgZ8DBxGMsXwTcL0Sg6S7r76CFStg8OBgujwpqKQgEleib8c1wFnuPtfdN7n7t+7+OnBOuEwkbcW2J8Dukdd0SapIfImSQit3X1txZjivVRQBiSRLfj4RNKTRAAAShklEQVQ0bw7HHBNMl43RrOojkfgSfTu21XKZSMrl50P//pSPx1xWUlD1kUh8iS5J7WpmSyuZb8DhEcQjkhQbNsCyZXDeebvnqaFZJLFESeEo4HvApxXmHwZ8HklEIklQsT0BdicF9ZIqEl+in0z3A9+5+yexD2BruEwkLeXnQ7NmcOyxu+fp6iORxBJ9Ozq4+17VR+5eAHSIJCKRJMjPh379oEmT3fPK2xRUfSQSV6JvR9MqljVLZiAiybJxI3zwwZ5VR6CSgkh1JPp2LDKzKyvONLOfAoujCUmkbt58E9z3Tgql6JJUkUQSNTSPBV40s4vYnQR6A42Bs6IMTKS28vODaqO+ffecr4ZmkcQSdYj3JdDfzE4EfhTOnhne1SySlvLzg4TQtELlp+5TEEksUUkBAHefC8yNOBaROtu0Cd57D269de9luk9BJDF9OySjvPUWlJbu7gQvlkoKIonp2yEZJT8/6NbiuOP2XqaSgkhi+nZIRsnPhz59YJ999l6mkoJIYvp2SMbYsgUKCva+FLVMWUlBROKLNCmY2VAzW2lmq83slirWO9bMSszs3Cjjkcz29ttQUlJFUsBBeUGkSpElBTPLBR4CTgO6AReYWbc46/0OeDWqWCQ75OdDbm7QXXZl3JUURBKJsqTQB1jt7mvcfScwFTizkvWuA14AvoowFskC8+YFHeC1aFH5cldGEEkoyqRwMHt2uV0YzitnZgcT3Bk9oaoNmdlVZlZgZgVFRUVJD1Qavq1bYdGi+FVHoJKCSHVEmRQqGwi34ldyHHCzu5dUtSF3n+juvd29d7t27ZIWoGSOBQtg164ESUEZQSShat3RXEuFwCEx0+3Ze2Ce3sDUcCD1/YHTzazY3adHGJdkoPx8yMmBAQPir6Orj0QSizIpLAI6mVlH4DPgfODC2BXcvWPZczObDPxVCUFqIz8fevWCVq3ir6OSgkhikVUfuXsxMIbgqqIVwJ/d/SMzG21mo6Par2Sf7dvhnXeqrjoCtSmIVEeUJQXcfRYwq8K8ShuV3f2yKGORzPXOO7BjR+X9HcVSSUEksUiTgkh9yM8HMxg4sOr11KYgkpi6uZAGb9486NkT2rSpej2VFEQSU1KQBm3HjuBy1ETtCaA2BZHqUFKQBu3dd4OG5hNPTLyuSgoiiSkpSIM2b17QnjBoUOJ11SGeSGJKCtKgzZ0btCe0bZt4XXfHKr3RXkTKKClIg7V9e9CekOhS1DKqPhJJTElBGqyy9oTqJgVA1UciCSgpSINVk/YE0H0KItWhpCANVtn9CdVpTwBVH4lUh5KCNEg1bU8AXX0kUh1KCtIg1aY9QdVHIokpKUiDNHduzdoTICgp6JJUkaopKUiDVNP2hHIqLIhUSUlBGpyy9oTqdG0RS9VHIokpKUiDU93xEyrS1UciiSkpSINT1p6QaPyEitSmIJKYBtmRBmf2bOjdu+btCao+EklMJQVpUL77DhYuhJNPrvlrVX0kkpiSgjQo+flQUgKnnJLqSEQyk5KCNCizZ0OzZtC/f81f6zjmalMQqYqSgjQof/sbHH88NGlS89eqTUEkMSUFaTA++wxWrKhdewKoTUGkOpQUpMGYPTv4W6ekoLwgUiUlBWkwZs+Gdu2gR4/ab0P3KYhUTUlBGgT3ICmcdBLk1PJTq+ojkcSUFKRB+Ogj+OKL2lcdgRqaRapDSUEahJdfDv6eemrdtqPqI5GqKSlIgzBrVtCW0L597behhmaRxJQUJO1t2gRvvgmnn1637ahNQSSxSJOCmQ01s5VmttrMbqlk+UVmtjR8vG1mR0UZjzRMs2dDcTGccUbdt6XqI5GqRZYUzCwXeAg4DegGXGBm3Sqs9k/gBHfvAdwFTIwqHmm4Zs2CNm3guOPqth2VFEQSi7Kk0AdY7e5r3H0nMBU4M3YFd3/b3TeGkwuBOtQYSyZyD5LCqadCXh07etfVRyKJRZkUDgY+jZkuDOfF81Pg5coWmNlVZlZgZgVFRUVJDFHS3fvvB5ei1rU9oYyqj0SqFmVSqOzbV+lPNTM7kSAp3FzZcnef6O693b13u3btkhiipLtZs4K/Q4fWfVuqPhJJLMqR1wqBQ2Km2wOfV1zJzHoAjwGnufuGCOORBmjGDOjTBw44oO7bUlIQSSzKksIioJOZdTSzxsD5wIzYFczsUGAacIm7r4owFmmACgvh3XfhrLOSt02NpyBStchKCu5ebGZjgFeBXOAJd//IzEaHyycAtwP7AQ+bGUCxu/eOKiZpWKZPD/4mKymopCCSWJTVR7j7LGBWhXkTYp6PAkZFGYM0XNOmQbducMQRydme42poFklAdzRLWvr662A85mRWHYlIYkoKkpZeeglKS+Hss5O3TZUURBJTUpC0NG0aHHYYHH108rapNgWRxJQUJO188w28+iqcey5YMn/YKyeIJKSkIGnnuedg1y646KLkblfVRyKJKSlI2pkyBbp2hZ49I9i4SgsiVVJSkLSybh288QZceGGSq45QSUGkOpQUJK0880zw98ILUxuHSLZSUpC04R5UHfXrB4cfHsH2VVIQSUhJQdLGokWwbBmMHJnqSESyl5KCpI1HH4XmzaOrOlJJQSQxJQVJC5s2wdSpcMEF0KpVqqMRyV5KCpIWpkyBrVvhP/4jun2opCCSmJKCpJw7PPxw0KXFMcekOhqR7BZp19ki1fHKK/DRR/DUU8m/NyGW4+SYfgeJVEXfEEm5e++Fgw+G886Ldj+llJKjj7xIlfQNkZRavBjmzoWxY6Fx42j3VWql5JIb7U5EGjglBUmp3/wmuNroyiuj31cppao+EklA3xBJmUWLgnETfv5zaN06+v25uUoKIgmooVlS5he/gHbt4IYb6md/SgoiiSkpSErMnh087r8fWrasn30qKYgkpqQg9W7HDhgzJuj0bvTo+tuvkoJIYkoKUu/uvRdWroSXX4amTetvv0oKIokpKUi9WroU7r4bRoyAoUPrd9+e4+TpIy9SJX1DpN5s2xZ0eNe2LTz0UP3v301JQSQRfUOkXrgH7QfLl8OrrwZXHdW7HFR9JJKA7lOQevGb3wR9G915JwwZkqIgciDP9DtIpCr6hkjk7r8fbrsNLroIfvnLFAaSA3muj7xIVVRSkMiUlgbJ4D//E849FyZNirYX1KoUlxQHSSFXSUGkKvqGSCTWr4dRo2DWLLjiimCozbwUftp27NwBQKOcRqkLQqQBUElBkmrr1qC6qFs3mDMHHnwQHnsstQkBYNvObQDk5eh3kEhVIk0KZjbUzFaa2Wozu6WS5WZm48PlS82sV5TxSDTcYdkyuOUW6NgxqC7q3TuYd+21qasyirV953ZASUEkkci+IWaWCzwEnAIUAovMbIa7L49Z7TSgU/joCzwS/pU0U1oKmzbBxo2wYQOsWQN//3tw4p8/H774AnJzgxvSbr4ZBg1KdcR7+nzD5wA0b9w8xZGIpLcofzb1AVa7+xoAM5sKnAnEJoUzgafc3YGFZtbGzA509/XJDuayu/6XJ796oJIliX7GejW2Hm8b1XltDTaXUNn+arqBRHFa/G1+fxd23g5ycnZgtpVXKeHV6cD0amyyHpXmlEJz+HG3H9fvjkUamCiTwsHApzHThexdCqhsnYOBPZKCmV0FXAVw6KGH1iqY77VtS7PVh8VZWvVJ0fY4g3nsggS80pXiv8yrXMvixblHeOGEJU5INTkv5+SUkJNTXP7Iy9tBo0bbMQOjKdAUqMGgCF7xfa0Zq2mdVCm03dGWkSeNrPU+RbJBlEmhsm9txTNVddbB3ScCEwF69+5dq5/fvxtzMb/j4tq8VEQka0TZ0FwIHBIz3R74vBbriIhIPYkyKSwCOplZRzNrDJwPzKiwzgxgZHgV0nHApijaE0REpHoiqz5y92IzGwO8CuQCT7j7R2Y2Olw+AZgFnA6sBrYCl0cVj4iIJBbpRdvuPovgxB87b0LMcweujTIGERGpPt3RLCIi5ZQURESknJKCiIiUU1IQEZFyFrT1NhxmVgR8UoOX7A98HVE46Swbjzsbjxmy87iz8Zihbsd9mLsnHAi3wSWFmjKzAnfvneo46ls2Hnc2HjNk53Fn4zFD/Ry3qo9ERKSckoKIiJTLhqQwMdUBpEg2Hnc2HjNk53Fn4zFDPRx3xrcpiIhI9WVDSUFERKpJSUFERMpldFIws6FmttLMVpvZLamOJwpmdoiZzTWzFWb2kZldH87f18z+ZmZ/D/+2TXWsyWZmuWa2xMz+Gk5nwzG3MbPnzezj8H/eL0uO+4bw8/2hmT1jZk0z7bjN7Akz+8rMPoyZF/cYzey/wnPbSjM7NVlxZGxSMLNc4CHgNKAbcIGZdUttVJEoBn7u7l2B44Brw+O8BZjj7p2AOeF0prkeWBEznQ3H/ADwirt3AY4iOP6MPm4zOxj4GdDb3X9E0BX/+WTecU8GhlaYV+kxht/x84Ejw9c8HJ7z6ixjkwLQB1jt7mvcfScwFTgzxTElnbuvd/f3wuebCU4SBxMc65Phak8CP0lNhNEws/bAGcBjMbMz/ZhbAccDjwO4+053/5YMP+5QHtDMzPKAfQhGaMyo43b3+cA3FWbHO8YzganuvsPd/0kwJk2fZMSRyUnhYODTmOnCcF7GMrMOwNHAO8D3ykaxC/8ekLrIIjEOuAkojZmX6cd8OFAETAqrzR4zs+Zk+HG7+2fA74F1wHqCERpfI8OPOxTvGCM7v2VyUrBK5mXs9bdm1gJ4ARjr7t+lOp4omdkw4Ct3X5zqWOpZHtALeMTdjwb+RcOvMkkorEc/E+gIHAQ0N7OLUxtVykV2fsvkpFAIHBIz3Z6gyJlxzKwRQUKY4u7TwtlfmtmB4fIDga9SFV8EBgDDzWwtQbXgj83s/8jsY4bgM13o7u+E088TJIlMP+6TgX+6e5G77wKmAf3J/OOG+McY2fktk5PCIqCTmXU0s8YEjTIzUhxT0pmZEdQxr3D3+2IWzQAuDZ9fCvylvmOLirv/l7u3d/cOBP/X1939YjL4mAHc/QvgUzM7Ipx1ErCcDD9ugmqj48xsn/DzfhJB21mmHzfEP8YZwPlm1sTMOgKdgHeTskd3z9gHcDqwCvgHcGuq44noGAcSFBuXAu+Hj9OB/QiuVvh7+HffVMca0fEPBv4aPs/4YwZ6AgXh/3s60DZLjvtO4GPgQ+BpoEmmHTfwDEGbyS6CksBPqzpG4Nbw3LYSOC1ZcaibCxERKZfJ1UciIlJDSgoiIlJOSUFERMopKYiISDklBRERKaekIBnHzO43s7Ex06+a2WMx038ws/9MsI23q7GftWa2fyXzB5tZ/ype9xMzuz18/iszuzF8Ps/MqjUou5m1M7NXqrOuSE0oKUgmepvgjlfMLAfYn6A3yTL9gbeq2oC7xz2pV8Pgsv3HcRPwcB22j7sXAevNbEBdtiNSkZKCZKK32H1SPpLghqfNZtbWzJoAXYElAGb2/8xskZktNbM7yzZgZlvCvzlm9nDYl/9fzWyWmZ0bs6/rzOw9M1tmZl3CTglHAzeY2ftmNig2MDPrDOxw96/jxH6xmb0djhvQJ3zNCeG23g87wmsZrjsduKjW75JIJfJSHYBIsrn752ZWbGaHEiSHBQQ9SPYDNgFL3X2nmQ0h6B6gD0EHYzPM7HgPujAuczbQAehO0EPlCuCJmOVfu3svM7sGuNHdR5nZBGCLu/++kvAGAO9VEX5zd+9vZseH+/kRcCNwrbu/FXZ8uD1ctwC4u7rvi0h1qKQgmaqstFCWFBbETJe1FwwJH0sITtRdCJJErIHAc+5e6kHfQ3MrLC/rgHAxQfJI5ECC7q/jeQbK+9ZvZWZtwmO5z8x+BrRx9+Jw3a8Ieg0VSRolBclUZe0K3QmqjxYSlBRi2xMM+K279wwfP3T3xytsp7IuimPtCP+WUL2S9zagaRXLK/Y74+7+P8AooBmw0My6hMuahtsTSRolBclUbwHDgG/cvcTdvwHaECSGBeE6rwJXhFUymNnBZlZxoJY3gXPCtoXvETQiJ7IZaBln2Qrgh1W89rwwloEEg8lsMrMfuPsyd/8dQZVRWVLoTJDwRJJGSUEy1TKCq44WVpi3qayR14PRu/4ELDCzZQTjE1Q8mb9A0GPlh8CjBKPabUqw75eAsypraAbmA0eHXUBXZmN4OewEgl4yAcaGDc8fEJQMXg7nnwjMTBCLSI2ol1SRBMyshbtvMbP9CPqsHxC2L9R2ew8AL7n77DrGNR8409031mU7IrF09ZFIYn8NG3wbA3fVJSGEfgP0rcsGzKwdcJ8SgiSbSgoiIlJObQoiIlJOSUFERMopKYiISDklBRERKaekICIi5f4/o2xceLJbfygAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(year_5_pdf['lbs'], year_5_pdf['posterior_cdf'], label = 'posterior_cdf',\n",
    "        linestyle = '-', color = 'k')\n",
    "plt.plot(year_5_pdf['lbs'], year_5_pdf['prior_cdf'], label = 'prior_cdf',\n",
    "        linestyle = '-', color = 'b')\n",
    "plt.plot(year_5_pdf['lbs'], year_5_pdf['likelihood_cdf'], label = 'likelihood_cdf',\n",
    "        linestyle = '-', color = 'g')\n",
    "plt.xlabel('Weight (lbs)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Bayesian Baby Weight Estimate')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max posterior: 0.14760105060194184\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4FFXWwOHfSUJYZF90UFBQQVAJCAwMixp3RhkRBz53EdxAUMEFcRZFGZfRGWFQFp1BVAZFcRBBEZElsgUFVxAE2YQgYgg7AbKd749biU1ISCBdqaT7vM/TD9VV1VWnukOfvvfWvVdUFWOMMQYgJugAjDHGlB2WFIwxxuSxpGCMMSaPJQVjjDF5LCkYY4zJY0nBGGNMHksKxlcicr6IrA46jvxEpJGIqIjEBR1LqGN5v0QkUURS/I6pOETkJhGZFXQcpuQsKUQYEdkoIgdEZJ+I7BSRD0WkYVDxqOoCVT3Lj2OXhWsVkfpecjkpZN2fC1k3s6jjhfP9EpHXRORvReyjIrLfew9zH4OLeM0RCVVVJ6rq5eGIu4DzJYnIHX4c2xzJkkJk+oOqVgXqA9uAFwOOx0+BXquqbgXWAheErL4A+L6AdfNLMbRj0VJVq4Y8ngs6IBMcSwoRTFUPAu8CZ+euE5GrROQrEdkjIptFZGjItg9F5N7QY4jItyJyjbfcTEQ+EZEdIrJaRP4vZL8rRWSliOwVkS0i8pC3/rAqDhEZIiLrvP1Wikj3kG23ichCEfmH98t/g4j83o9rDdFHRH4Ska0i8qD3ut+ISLqI1Ak5VhsRSRWRCgUcYz5eAhCRWOA84F/51nXw9kNEKnrXuElEtonIWBGpXMj71dq7hr0iMllE3s7/619EHhSRX7xr6O2tuwu4CRjs/fqfXpz3Md9x24nIMu/92yYiL4RcL8Au79gdcj+7kNeqiNwjIj94sQ8TkTNEJNk73jsiEu/tW0tEPvDe353ecgNv21PA+cBL3rle8tYX+rdoSkhV7RFBD2AjcKm3XAV4HXgjZHsi0AL3gyAB9+v6Gm/b/wGfhezbEkgD4oETgM1AbyAOaA1sB87x9t0KnO8t1wJah5wvJeSYPYGTvfNfB+wH6nvbbgMygTuBWKAf8BMgPlxrI0CBt7xrawGkhhxvBtAv5FjDgRcLiaMX8I233Bb3pdkk37oDQLz3fAQwDagNVAOmA8/kf7+89/1H4H6gAnAtkAH8LWTfLOBJb/uVQDpQy9v+Wu6+R/l7UeDMQrYlA7d4y1WB3+V77+JC9r0NWJjvuNOA6sA5wCFgDnA6UANYCfTy9q0D/NH7DKsBk4GpIcdKAu4IeX7Uv0V7lPA7JOgA7BHmD9R9Ue4DdnlfGD8BLY6y/whguLdcEdgBNPGe/wMY7S1fByzI99qXgce95U3A3UD1fPvkfckVcv6vgW7e8m3A2pBtVbwvl9/4cK25X2zNQrY/B4wLud5F3nIs8DPQrpDjNgKycclwEPCUt35LyLp53jrBJcIzQl7fAdiQ//3ClTS2EJIUgYUcnhQOcPiX8y/8+uX9GsVLCnu89zD3cYW3bT7wBFC3gOstTlLoFPL8C+CRkOf/BEYUElMrYGfI8yQOTwpH/Vu0R8keVn0Uma5R1Zq4L/kBwKci8hsAEWkvIvO8ovpuoC9QF0BVDwHvADeLSAxwAzDBO+ZpQHsR2ZX7wFVP/Mbb/kfcL9UfReRTEelQUGAicquIfB1yjHNzz+/5OXdBVdO9xarhvtYQm0OWf8SVYgDeB84WkdOBy4Ddqvp5QQGo6kYgBeiM+yJf4G1KDlmXW+VSD5fsvgh5D2Z66/M7Gdii3rdeAfECpKlqVsjzdI7+fhWktarWDHl87K2/HWgKfC8iS0Wk6zEed1vI8oECnlcFEJEqIvKyiPwoIntw71VNr9qtIEX9LZoSsKQQwVQ1W1Wn4H7FdvZWv4kr1jdU1RrAWNyv11yv4/6DXQKkq2qyt34z8Gm+L4+qqtrPO9dSVe0GnAhMxSWXw4jIacC/cV/edbwv8xX5zl+a1woQerfSqbjSBuraKN7BvRe38GtyLMwC3Jd/B2BxvnWd+TUpbMd9IZ4T8j7WUNdYnt9W4BQRCY35WO6uKtEQyKr6g6regPtM/w68KyInlPS4BXgQOAtor6rV+bWBPve685/vqH+LpmQsKUQwcbrhqjBWeaurATtU9aCItANuDH2NlwRycMX70C/CD4CmInKLiFTwHr8VkeYiEi/uPvUaqpqJq47ILiCk3C+UVC++3riSQiDX6vmr90v1HFwd9dsh297AVYtcDfy3iBDmA7cCP6nqHm/dQm9dDVypAVXNwSXG4SJyohf7KSJyRQHHTMa9jwNEJM67vnZFxBFqG64O/7iIyM0iUs+LeZe3Ohv3+eWU5Nj5VMMlyl0iUht4PN/2/NdR6N9imOKJapYUItN0EdmH+3J+Cteg95237R7gSRHZCzxGAb/ocV+GLQj5IlTVvcDlwPW4X9M/4349VvR2uQXY6BX/+wI35z+oqq7EJZtk3H/0FsCiEl1pya/1U9wtpXOAf6hqXgcsVV2E+/L70qsiOppPcb+oF4as+xqoDHwRUhUG8Ih3ziXe+zUb90v5MKqagWtcvh33pXwz7gvxUBGx5BqHqwLbJSJTj7LfN3J4P4UR3vouwHfe+/sv4HpVPehdy1PAIu/YvytmPIUZgXuftgNLcNVpof4F9PDuTBpZjL9FUwJyeHWlMa7eH7hLVTsXuXOEE5G5wJuq+p+gYwEQkc+Asao6PuhYTGSykoI5jIhUwf3CfiXoWIImIr/F3e74dlH7+hjDheL6TcSJSC/crbVF9ow25nhZUjB5vHrtVFzVzpsBhxMoEXkdV60z0KuuCMpZwDfAblyDbA91vaiN8YVVHxljjMljJQVjjDF5ytSwwcVRt25dbdSoUdBhGGNMufLFF19sV9WCOkkeptwlhUaNGrFs2bKgwzDGmHJFRH4szn5WfWSMMSaPJQVjjDF5LCkYY4zJU+7aFAqSmZlJSkoKBw8eDDoUU4BKlSrRoEEDKlQoaH4aY0xZEhFJISUlhWrVqtGoUSMOH1DSBE1VSUtLIyUlhcaNGwcdjjGmCL5VH4nIq+KmCFxRyHYRkZEislbclI+tj/dcBw8epE6dOpYQyiARoU6dOlaKM6ac8LOk8BrwEm7EzYL8HjdlYROgPTDG+/e4WEIou+yzKcOSkyEpCRIT3fPc5Q4dCt8WDctl8fo7FDhvVdj5lhRUdb6INDrKLt1w8+kqbgjhmiJS38Z1McZnuV92deqQM/B+OJRBTqyrNIjJzia7QhwfD/g9V7w0g9jMbLJjYxBAsnPIiYLl7AqxfNT/Cn4/6mNiM7PLTEwP9rqF07tez6A/XO7rn0eQbQqncPjUgineuiOSgojcBdwFcOqpp5ZKcMZEpORk9JJL0EMHyUYRhTgFzckBXH1ydkYmcVOnIRkQo6A5br6kWIAoWM7JyKLS+x8SkwGxWnZiqrp2PO98fpLvSSHIW1ILqlMocHQ+VX1FVduqatt69YrspV2mPfbYY8yePbtUznXbbbfx7rvvArBgwQLOOeccWrVqxYEDB0rl/KbsOTh7JjmHDhCT4xKCxMS4X6Tx8cRUjEdjY4mrVJlLHx5NbKXKaGwssfHxxFasCFGyXKFSZa4Y/DJxlSoHGgcV4smSimQSSzaV+eM1i0ke9ozvfyNBlhRSOHy+2QZ48+NGquzsbJ588sljfk1sbGHzlxffxIkTeeihh+jdu3eJj2XKp+ycbAZnzuDZGKgkMcTGV0RGjIC0tCPqr+M6dICEVsHXowe13KEDtGgRyLk3boT5w5IYsyqRE+vB0MQkWt6fSNtOpdOm4OvQ2V6bwgeqesQ8vCJyFW4C9ytxDcwjVbXI+Wfbtm2r+cc+WrVqFc2bu+lZBw6Er78uceiHadUKRow4+j4bN26kS5cutG/fnq+++oqmTZvyxhtvcPbZZ9OnTx9mzZrFgAEDmDlzJl27dqVHjx7MmTOHhx56iKysLH77298yZswYKlasSKNGjQ57zfXXX3/E+dauXUvfvn1JTU0lNjaWyZMnc/rpp3Pvvfcyd+5cGjdujKrSp08fdu3axeDBg6lRowYdO3Zk4sSJ4X2DiiH0MzIBSE5m/utP8kjWTP7U+VH+sLVaqTZemqL98gs88QS8/DJUqQKPPuq+zypXDs/xReQLVW1b1H6+lRRE5C0gEagrIim4ybgrAKjqWGAGLiGsBdJxk6aXa6tXr2bcuHF06tSJPn36MHr0aMB13lq40E3dO3OmmzTr4MGD3HbbbcyZM4emTZty6623MmbMGAYOHHjEawpy0003MWTIELp3787BgwfJycnhvffeY/Xq1Sxfvpxt27blJaQ77riDhQsX5iUjE2WSk9FLLqbjwYMkVYih4u1/gNssGZQVBw64H53PPAPp6XD33fD443DiicHE4+fdRzcUsV2B/uE+b1G/6P3UsGFDOnXqBMDNN9/MyJEjAbjuuuuO2Hf16tU0btyYpk2bAtCrVy9GjRqVlxQKek2uvXv3smXLFrp37w64BAIwf/58brjhBmJjYzn55JO5+OKLw3dxpvxKSkIPHSJOITZbXDWFlRACl5MDEyfCn/8MmzdDt27w7LPQrFmwcdnYR2GU/3783OcnnHDCEfsWVW1X0GuK81rrE2Dyy77gfA7FQlYMSHz8r3XYJjALFkC7dnDrra5EkJQEU6cGnxDAkkJYbdq0ieTkZADeeustOnfuXOi+zZo1Y+PGjaxduxaACRMmcOGFFxbrPNWrV6dBgwZMnToVgEOHDpGens4FF1zApEmTyM7OZuvWrcybN6+EV2QiwdzfHODiW5VV994Ac+ZYKSFA69dDjx5wwQWwbRtMmACffw7F/K9fKiwphFHz5s15/fXXSUhIYMeOHfTr16/QfStVqsT48ePp2bMnLVq0ICYmhr59+xb7XBMmTGDkyJEkJCTQsWNHfv75Z7p3706TJk1o0aIF/fr1K3aSMZFt2uppfNO4Mmc+P84SQkB274ZHHoHmzeGjj+DJJ2H1arj5ZogpY9/Cvt595Iei7j4KysaNG+natSsrVhQ41FPUKwufUTRSVc4YeQbnnHgO02+YHnQ4UScrC8aNg7/+FVJToVcveOopOOWU0o+luHcflbEcZYwJp++3f8+GXRu4qslVQYcSdT75BM47D/r2dW0Fy5bBa68FkxCOhSWFMGnUqJEvpYT+/fvTqlWrwx7jx48P+3lMZFqwaQEAl55+acCRRI/Vq+EPf4DLL4f9++Hdd+HTT6FNm6AjK56ImE8hko0aNSroEEw5lpySTN0qdTmj1hlBhxLxduxwnc9Gj3Ydzv7+d7jvPvDuGC83LCkYE8EWb15Mx4Yd7VZlH2VmwpgxMHSoa1C+807XkBxU57OSsuojYyLU7nkzuXbqGnrsKuOV2OWUKnzwgRsi6f77XfXQ11/D2LHlNyGAJQVjIlNyMlWvvIZhc+HGQa+6ORRM2KxYAVdc4doOVGH6dJg1yyWI8s6SgjGRKCkJycggTiEmM8t1mTUllpYG/ftDy5bubqIRI1yC6NoVIqWGzpJCmFStWhWAn376KW/Quddee40BAwYU+xih+48dO5Y33nAzmSYmJpK/b0Y4hOO4jRo1Yvv27QCMHDmS5s2bc9NNN4UjPFMSiYlkxsWQKTa0RThkZbkG5KZNXfVQv37www+u2qhChaCjCy9raA6zk08+OW9im5I4lt7NZcXo0aP56KOPaNy4cdChmA4d6DvwTDpvyOL2QROsJ3MJJCW5u4iWL4eLLoJ//SsyqokKE3FJYeDMgXz9c3gnVGj1m1aM6FK84VcL69n84Ycf8re//Y3p06ejqvTt25dNmzYBMGLEiLzRVXMNHTqUqlWr8tBDDwEwefJk7rnnHnbt2sW4ceM4//zzOXjwIP369WPZsmXExcXxwgsvcNFFFxW6/sCBA/Tu3ZuVK1fSvHnzImdgmzlzJn/605/Izs6mbt26zJkzh7S0NG644QZSU1Np165d3uB8ffv2Zf369Vx99dX06dOHQYMGFev9Mv7IzslmUvUfqX37PZYQjtOPP8JDD7l+Bqed5v699trIqSYqTMQlhbLovffe44UXXmDGjBnUqlWLG2+8kUGDBtG5c2c2bdrEFVdcwapVq456jKysLD7//HNmzJjBE088wezZs/P6MCxfvpzvv/+eyy+/nDVr1hS6fsyYMVSpUoVvv/2Wb7/9ltatWxd6vtTUVO68807mz59P48aN2bFjBwBPPPEEnTt35rHHHuPDDz/klVdeAVx118yZM5k3bx5169YNx9tmSuDH3T9yMOsgzevZ0CLHKj3dDWH9/PNuXKJhw+DBB8M32U1ZF3FJobi/6EvLvHnzWLZsGbNmzaJ69eoAzJ49m5UrV+bts2fPHvbu3XvU41x77bUAtGnTho0bNwKwcOFC7r33XsCNunraaaexZs2aQtfPnz+f++67D4CEhAQSEhIKPd+SJUu44IIL8qqCateuDbg5G6ZMmQLAVVddRa1atY7p/TCl44e0HwBoUrtJwJGUH6rw9tvw8MOQkgLXXw/PPQcNGxb92kgScUmhrDn99NNZv349a9asoW1bNxZVTk4OycnJVD6Gnx4VK1YEIDY2lqysLKDweRXCMd+Cqha6r3WEKvvW7nBDsjepY0mhOL76yjUaL1jgpt998004//ygowqG3X3ks9NOO40pU6Zw66238t133wFw+eWX89JLL+Xt8/VxTip9wQUX5M23vGbNGjZt2sRZZ51VrPUrVqzg22+/LfTYHTp04NNPP2XDhg0AedVHocf46KOP2Llz53HFbvz1w44fOKHCCdSvWj/oUMq01FQ3/WWbNrBqlZsfedmy6E0IYEmhVJx11llMnDiRnj17sm7dOkaOHMmyZctISEjg7LPPZuzYscd13HvuuYfs7GxatGjBddddx2uvvUbFihULXd+vXz/27dtHQkICzz33HO3atSv02PXq1eOVV17h2muvpWXLlnnTgz7++OPMnz+f1q1bM2vWLE499dTjit3464cdP3Bm7TOtVFeIzEx3F1HTpm5o6/vugzVr4K67IDY26OiCZfMpmFJhn1HpOuuls0g4KYHJPScHHUqZ88knMHAgrFwJl13mOqCdfXbQUfnP5lMwJkrlaA4bdm7gzFpnBh1KmbJxo7ul9PLL4eBBNyfyxx9HR0I4FtbQbGjfvj2HDh06bN2ECRNoEck9dCLYz/t+JjMnk9NqnhZ0KGXCwYPu9tKnn3a3mD71FDzwQPkb0rq0WFIwfPbZZ0GHYMJo027XKfLUGtbe8+GH7q6ideugRw/45z/BmsGOzqqPjIkwuUmhYfUou8E+xPr1bgTTrl3d2ESffAKTJ1tCKA5LCsZEmM27NwPRWVI4cAAef9y1E8yb5zqfffMNXGqzkRabVR8ZE2E27d5E9YrVqVGpRtChlBpVeP99GDTINSjfcINrRzjF5hc6ZlZSMCbCbNqzKapKCT/8AFdeCd27wwknuBLCm29aQjhelhTKkKlTpx42JlJxTZs2jWeffdaHiI4UOudDamoq7du357zzzmPBggWlcn5TtE27N0VFe8L+/fDnP8O558KiRfDCC264Cps6omSiNykkJ8Mzz5SpaQqPJylkZWVx9dVXM2TIkGN6TTjMmTOHZs2a8dVXX3F+NI8LUMZs2h3ZJQVVN4x18+buNtPrrnO9kQcNirwJb4IQnUkhORkuuQT++lf3bxgSw8aNG2nWrBm9evUiISGBHj16kJ6ezpw5czjvvPNo0aIFffr0yesPMGTIEM4++2wSEhJ46KGHWLx4MdOmTePhhx+mVatWrFu3jnXr1tGlSxfatGnD+eefz/fffw/AbbfdxgMPPMBFF13EI488ctiv9x9//JFLLrmEhIQELrnkkrw5G/K/piD79u2jd+/etGjRgoSEBP73v/8BMH78eJo2bcqFF17IokWLADde0+DBg5kxYwatWrUqcm4GUzrSM9PZnr49YpPC99+7zmc9e0KtWm4AuzfegN/8JujIIoiqlqtHmzZtNL+VK1cese6onn5aNTZWFdy/Tz99bK8vwIYNGxTQhQsXqqpq7969ddiwYdqgQQNdvXq1qqrecsstOnz4cE1LS9OmTZtqTk6Oqqru3LlTVVV79eqlkydPzjvmxRdfrGvWrFFV1SVLluhFF12Ut99VV12lWVlZqqo6fvx47d+/v6qqdu3aVV977TVVVR03bpx269atwNcUZPDgwXr//ffnPd+xY4f+9NNP2rBhQ/3ll1/00KFD2rFjx7xzhZ63KMf8GZnjsnr7amUo+sbXbwQdSljt3as6eLBqhQqqNWqovviiamZm0FGVL8AyLcZ3rK8lBRHpIiKrRWStiBxRvyEiNURkuoh8IyLfiUhvP+PJk5gI8fFu5Kswzl/bsGHDvBnUbr75ZubMmUPjxo1p2rQpAL169WL+/PlUr16dSpUqcccddzBlyhSqVKlyxLH27dvH4sWL6dmzJ61ateLuu+9m69atedt79uxJbAEjdyUnJ3PjjTcCcMstt7Bw4cIiX5Nr9uzZ9O/fP+95rVq1+Oyzz0hMTKRevXrEx8fnDYxnyqZ9n37CkAXQfO2uoEMJC1X43/+gWTN3e+nNN7uqogEDIM7unfSFb2+riMQCo4DLgBRgqYhMU9XQSvP+wEpV/YOI1ANWi8hEVc3wKy7ATU84Z46bfDUxMWzTFRZ3RMq4uDg+//xz5syZw6RJk3jppZeYO3fuYfvk5ORQs2bNQofVPuGEE445pqJeo4XMoWAjbZYTyckk3PQACYcgZtFgOLltuZ6Kc906uPde+OgjaNnSdT4rx5dTbvhZUmgHrFXV9d6X/CSgW759FKgm7lunKrADCE8raFE6dIBHHw3rX9mmTZtI9ton3nrrLS699FI2btzI2rVuwpMJEyZw4YUXsm/fPnbv3s2VV17JiBEj8r74q1WrljcDW/Xq1WncuDGTJ7tRLlWVb775psgYOnbsyKRJkwCYOHEinTt3Lnb8+ed52LlzJ+3btycpKYm0tDQyMzPz4jFlUFISMRmZxClIZqb70VMOHTrkpsA891zXZjB8uJvjwBJC6fAzKZwCbA55nuKtC/US0Bz4CVgO3K+qOfkPJCJ3icgyEVmWmprqV7wl1rx5c15//XUSEhLYsWMHgwYNYvz48fTs2ZMWLVoQExND37592bt3L127diUhIYELL7yQ4cOHA3D99dfz/PPPc95557Fu3TomTpzIuHHjaNmyJeeccw7vv/9+kTGMHDmS8ePHk5CQwIQJE/jXv/5V7Pj/8pe/sHPnTs4991xatmzJvHnzqF+/PkOHDqVDhw5ceumlR53X2QQsMZGsCrFkCmGtFi1Nn3wCLVrAY49Bt26uYXngQKsqKlXFaXg4ngfQE/hPyPNbgBfz7dMDGA4IcCawAah+tOOGpaHZBxs2bNBzzjkn6DDKrLLwGUWDJ57pon+/qqbq4sVBh3JMtmxRve46d+/HmWeqfvxx0BFFHspAQ3MKENqDpgGuRBCqNzDFi3mtlxSa+RiTMRFt4SnZTLnmrHJT15KV5WZAa9bMzW/wxBOwfLm77dQEw89C2VKgiYg0BrYA1wM35ttnE3AJsEBETgLOAtb7GJNvGjVqxIoVK4IOo1jGjx9/RLVSp06dGDVqVEARmXDZum8rZ9YuH5PrLFkC/frB119Dly7w0ktwxhlBR2V8SwqqmiUiA4CPgVjgVVX9TkT6etvHAsOA10RkOa4K6RFV3X6c57O7ZIqpd+/e9O5dOnf/ArlVhaYU/LzvZzo3LP7NBUHYscPd4/Hvf8PJJ7veyddeC/bft2zwtflGVWcAM/KtGxuy/BNQ4oJipUqVSEtLo06dOpYYyhhVJS0tjUo2zZXvMrIz2J6+nfrV6gcdSoFycuD112HwYNi5081+9vjjUK1a0JGZUBHRpt+gQQNSUlIoy3cmRbNKlSrRoEGDoMOIeL/s/wWA31Qte2M+LF8O99wDCxdCp04wZoy7y8iUPRGRFCpUqEDjxo2DDsOYQG3d63q8169adkoK+/a5xuPhw6FmTXj1VejVy82VbMqmiEgKxhjXngBlp6Tw/vtuOIqUFLjzTjcocZ06QUdlimL52pgIUVaSwubNbsKba66B2rVh8WJ45RVLCOWFJQVjIsT2dHfjXt0qdQM5f1YWjBjh5keeNctNh2nDU5Q/Vn1kTIRIO5BG5bjKVK5QudTP/cUXcNdd8OWXbmrMUaOgUaNSD8OEgZUUjIkQOw7soE6V0q2j2bvXjU3Urh1s3epGMv3gA0sI5ZmVFIyJEGkH0qhTufSSwtSpbmjrLVtcz+Snn4YaNUrt9MYnVlIwJkKkpaeVSklh82bXiNy9u2tITk521UWWECKDJQVjIoTfJYX8DcnPPecaktu39+2UJgBWfWRMhEhL9y8pWENy9LCSgjERQFV9aWgObUj+6Sd45x1rSI50VlIwJgLsPrSbbM0Oa0nBGpKjk5UUjIkAaelpANSuXLvEx8rfkLx4sTUkRxNLCsZEgLQDLimUpPooO9vNgpa/Ifl3vwtXlKY8sOojYyJAbknheKuPli+HO+6Azz93s6CNHg028HB0spKCMRHgeEsKBw/CX/4CrVvD+vUwcSLMmGEJIZpZScGYCHA8JYX5892Q1mvWwK23wj//CXWDGUvPlCFWUjAmAqQdSEMQalaqWeS+u3bB3XfDhRdCZiZ8/LGbJtMSggFLCsZEhB0HdlCrci1iY2KPut+UKa4h+T//gQcfdG0Jl5d4lnQTSaz6yJgIUNQQF1u2uFnQpk6FVq1g+nRo06YUAzTlhpUUjIkAhQ2Gl5MDY8e60sHMmfD3v7s7jCwhmMJYScGYCJB2II36Vesftu77711D8sKFcPHF8PLLcOaZAQVoyg0rKRgTAUJLChkZMGwYtGwJ330Hr74Ks2dbQjDFYyUFYyJA2oE0aleqTXKyKx189x1cd53roXzSSUFHZ8oTKykYU85lZGewL2MfS+bVoVMn2L3bNSRPmmQJwRw7SwrGlHOTprmOa0vm1aF/f1i5Erp2DTgoU25ZUjCmnEpNhRtvhF59XVIY9mgdXnzKyU0DAAAad0lEQVQRqlULODBTrlmbgjHljKqrGrrvPldVNOzW+exfAFecsy3o0EwEsKRgTDmSO+HN9OluNrQ3703mtDsHwSGIWfQwnNIWOnQIOkxTjvlafSQiXURktYisFZEhheyTKCJfi8h3IvKpn/EYU16pwr//7TqhzZ4N//iHm/zmjM1JxGRkEqcgmVmQlBR0qKac862kICKxwCjgMiAFWCoi01R1Zcg+NYHRQBdV3SQiJ/oVjzHl1fr17jbTuXMhMdElh7w+B4mJZFeIIzsjk7j4eLeDMSXgZ0mhHbBWVderagYwCeiWb58bgSmquglAVX/xMR5jypXsbBgxAlq0gKVL3XAVc+bk64TWoQMv/70nT14ah8yZY1VHpsT8TAqnAJtDnqd460I1BWqJSJKIfCEitxZ0IBG5S0SWiciy1NRUn8I1puxYtQo6d4ZBg9yP/+++c8NdxxTwP/arRpUY3+UkSwgmLPxMClLAOs33PA5oA1wFXAH8VUSaHvEi1VdUta2qtq1Xr174IzWmjMjMhKeeciOZ/vAD/Pe/8MEH0LBh4a9JO1DwYHjGHA8/7z5KAUL/lBsAPxWwz3ZV3Q/sF5H5QEtgjY9xGVMmffkl9OkD33wD//d/8OKLcGIxWtnSDqRRu3Jt/wM0UcHPksJSoImINBaReOB6YFq+fd4HzheROBGpArQHVvkYkzFlzsGD8Oij7hbTbdvgvffg7beLlxDAGwzvGKbhNOZofCspqGqWiAwAPgZigVdV9TsR6ettH6uqq0RkJvAtkAP8R1VX+BWTMWXNwoVw++1unuTevd08ybVqHdsxippgx5hj4WvnNVWdAczIt25svufPA8/7GYcxZc2+fa50MGoUnHoqzJoFl1127MdRVXYc2GFtCiZsbOwjY0rZJ5/Auee6hDBgAKxYcXwJAWDPoT1k5WRZScGEjSUFY0rJzp2uIfnyy6FSJViwAEaOhKpVj/+YaQfcYHhWUjDhYknBmFIwfTqccw688QYMGQJffw2dOpX8uGnpXlKwkoIJE0sKxvhoxw649Va4+mqoWxc++wyeecaVFMLBSgom3CwpGOOTadNc6eCtt+Cxx2DZMmjTJrzn2HFgB2AlBRM+NnS2MWGWlgb33w8TJ0LLljBjBpx3nk/nSreSggkvKykYE0ZTp7rSwdtvw9Ch8Pnn/iUE+LX6qFalY+zcYEwhjpoUROS1kOVevkdjTDm1fbubGrN7d6hf31UVPf44xMf7e9609DRqVqpJbEysvycyUaOokkLLkOX7/QzEmPJqyhRXOnj3XXjySVc6aNmy6NeFg/VmNuFWVJtC/lFNjTGe1FS4915XVdS6teuUlpBQujHYCKkm3IpKCg1EZCRuGOzc5Tyqep9vkRlThr37LtxzD+zaBX/7GwweDBUqlH4caelpnHiCTVhowqeopPBwyPIyPwMxpjz45Rc3NMXkye720rlz3ZAVQUk7kEbzes2DC8BEnKMmBVV9vbQCMaYsU3WJoH9/2LMHnn4aHn4Y4gK+qduGzTbhVuQtqSLSS0S+FJH93mNZYdNmGhOJtm2Dnj3huuugcWM3Gc6jjwafEDKyM9ibsdeSggmro/5Ze1/+A4EHgC9xbQutgedFBFV9w/8QjQmGqmtEHjAA9u6FZ5+FBx8MPhnkyuvNbA3NJoyKKincA3RX1XmqultVd6nqXOCP3jZjItLPP8Mf/wg33ABnnAFffQWPPFJ2EgLYYHjGH0UlheqqujH/Sm9ddT8CMiZIqvDmm67fwYwZ8NxzsGgRnH120JEdyQbDM34o6nfPgePcZky5s20b9O3rhqr43e9g/Hho1izoqApng+EZPxSVFJqLyLcFrBfgdB/iMSYQ77zj+h3s2+dKBw88ALFlfOQIGwzP+KGopNASOAnYnG/9acBPvkRkTCnavt3dZvrOO/Db38Lrr0PzcnLbf271Ue3KtQOOxESSotoUhgN7VPXH0AeQ7m0zptx6/33XdvDee/DUU7B4cflJCOBKCvGx8ZxQ4YSgQzERpKiSQiNVPaL6SFWXiUgjXyIyxmc7d7r5DiZMcAPXzZpVegPYhVPuYHgiEnQoJoIUVVI42qSBlcMZiDGl4aOP3LAUb74Jf/1r6Y5oGm42GJ7xQ1FJYamI3Jl/pYjcDnzhT0jGhN+ePXDnnXDllVCzJixZ4oa59nu+Az/ZEBfGD0VVHw0E3hORm/g1CbQF4oHufgZmTLjMnQu9e0NKiuuANnQoVDpaGbicSDuQRrO6ZfieWVMuFTUg3jago4hcBOSOBfmh16vZmDJt/36XBEaNgqZNYeFC6NAh6KjCx0oKxg/F6rSvqvOAeT7HYkzYLFwIt90G69fDwIHu7qIqVYKOKnxU1WZdM74ocpRUY8qTAwfcoHUXXAA5OZCUBMOHR1ZCANibsZesnCxraDZhV4aG9zKmZD77DHr1gtWroV8/1zO5atWgo/KHDYZn/GIlBVPuHTrk5jfo2BHS091cyaNHR25CABs22/jHSgqmXPvyS1c6WLEC+vSBF16AGjWCjsp/NsSF8YuvJQUR6SIiq0VkrYgMOcp+vxWRbBHp4Wc8JnJkZrpbS9u3h7Q0+OADGDcuOhICWPWR8Y9vSUFEYoFRwO+Bs4EbROSIUem9/f4OfOxXLCayLF/uksETT8D117tSwlVXBR1V6arw+TKGLICTlq8POhQTYfwsKbQD1qrqelXNACYB3QrY717gf8AvPsZiIkBWFjzzDLRpA1u2wJQpbvyi2tFWg5KcTLf+Ixk2F2p17QnJyUFHZCKIn0nhFA4fcjvFW5dHRE7B9Ywee7QDichdIrJMRJalpqaGPVBT9q1eDZ07w5/+BNdc40oH3aO1T31SEjGZ2cQpSEaGu+/WmDDxMykUNHSj5ns+AnhEVbOPdiBVfUVV26pq23r16oUtQFP2qboeyeedBz/8AJMmubkPovrPIDGRrLgYsmJwgzclJgYdkYkgft59lAI0DHnegCMn5mkLTPKG/q0LXCkiWao61ce4TDmxZYu7o2jWLOjSxTUkn3xy0FGVAR06MGRIW85cnkL/hydH1tgdJnB+JoWlQBMRaQxsAa4HbgzdQVUb5y6LyGvAB5YQDMDbb7sOaIcOwZgxcPfdYNMG/Gphg2y+b9qC/pYQTJj5Vn2kqlnAANxdRauAd1T1OxHpKyJ9/TqvKd927oQbb3R3FTVtCl9/DX37WkLIzwbDM37xtfOaqs4AZuRbV2Cjsqre5mcspuz75BM3xPW2bTBsGAwZAnHWvbJAaQfSqFulbtBhmAhkw1yYwKWnw733wuWXQ7VqbgKcv/zFEkJhMrIz2HNojyUF4wtLCiZQS5dC69bw0ktu3uQvv3T9EEzhrDez8ZMlBROIzEzXI7lDBzcZzuzZMGIEVLaZv4u0PX07gJUUjC+sgG5K3erVcMstrpRw002ulFCzZtBRlR+WFIyfrKRgSk1oR7R161wntP/+1xLCscodIdWSgvGDlRRMqbCOaOFjJQXjJyspGN+9/Ta0aOHmTR4zBmbMsIRQErlJwSbYMX6wpGB8Yx3R/LE9fTvVK1YnPjY+6FBMBLKkYHzxySeudDB5suuItnAhNGkSdFSRYXv6dqs6Mr6xpGDCyjqi+W97+nbro2B8Y0nBhI11RCsdVlIwfrKkYErMOqKVLksKxk9WqDclEtoR7eab4cUXrd+B3ywpGD9ZScEcl4I6ok2YYAnBbwezDrI/c78lBeMbKymYY2Yd0YKTOxieJQXjFyspmGMyaZJ1RAuS9WY2frOkYIpl50644Qb3sI5owbGkYPxmScEUKbcj2rvvWke0oFlSMH6zpGAKZR3Ryp7U9FTAJtgx/rGkYApkHdHKpm37thEjMVZSML6xpGAOk5kJQ4daR7Syatv+bdStUpfYmNigQzERyioCTB7riFb2bdu/jZNOOCnoMEwEs5KCsY5o5ci2fds4qaolBeMfSwpRbssW1wFtwAC48EJYvhx69gw6KlMYKykYv1lSiGLWEa18UVVXUrCkYHxkSSEKWUe08mlfxj4OZB2w6iPjK0sKUcY6opVf2/ZvA7CSgvGVJYUoEdoRrXp164hWHm3b5yUFKykYH1lSiAL5O6J98YV1RCuPrKRgSoMlhQhmHdEii5UUTGnwNSmISBcRWS0ia0VkSAHbbxKRb73HYhFp6Wc80WT1aujUyU2TecMN7lbTSy4JOipTEhU+X8aQBXDiN+uCDsVEMN9qlEUkFhgFXAakAEtFZJqqrgzZbQNwoaruFJHfA68A7f2KKRqowujR8PDDrkQweTL06BF0VKbEkpO59YE3iMmEuMVXwJw5rghoTJj5WVJoB6xV1fWqmgFMArqF7qCqi1V1p/d0CdDAx3giXv6OaCtWWEKIGElJxGZlEadARgYkJQUdkYlQfiaFU4DNIc9TvHWFuR34qKANInKXiCwTkWWpqalhDDFyFNQRrX79oKMyYZOYSGaskBUDxMdDYmLQEZkI5WdSKKgrlBa4o8hFuKTwSEHbVfUVVW2rqm3r1asXxhDLvx07rCNaVOjQgR531+L9m9pa1ZHxlZ9JIQVoGPK8AfBT/p1EJAH4D9BNVdN8jCfiWEe06JGZncmMujv5tveVlhCMr/xMCkuBJiLSWETigeuBaaE7iMipwBTgFlVd42MsESW0I1qNGtYRLRr8vO9nFKVBdWt2M/7y7WtEVbNEZADwMRALvKqq34lIX2/7WOAxoA4wWlx9R5aqtvUrpkiwdKmb82D1atcR7ZlnrN9BNEjZkwLAKdWP1ixnTMn5+ttSVWcAM/KtGxuyfAdwh58xRIrMTHjqKfjb31wD8uzZ1u8gmmzZuwXASgrGd1bhUA7YjGgmr6RQzUoKxl82zEUZpurGK8qdEW3yZJsRLVpt2bOFSnGVqF25dtChmAhnJYUyassW6NMHZs1yHdJefdX6HUSzLXu3cEq1UxC719j4zEoKZZB1RDP5pexJsUZmUyosKZQh1hHNFGbzns3WyGxKhSWFMsI6opnCZGZnsnn3Zk6veXrQoZgoYEkhYNYRzRRl857NZGs2jWs1DjoUEwUsKQQodEa0gQNtRjRTsA07NwBwei0rKRj/WVIIQEEzog0fbj2TTcHW71wPQOOaVlIw/rNKilJmHdHMsdqwawNxMXHW0GxKhZUUSol1RDPHa/3O9ZxW4zRiY2KDDsVEASsplILQjmi//z2MG2f9Dkzxrdu5ztoTTKmxkoLP8ndE+/BDSwim+FSV77d/T/O6zYMOxUQJSwo+sY5oJhxS9qSwL2Mfzeo2CzoUEyUsKfhg1izriGbCY9X2VQA0r2clBVM6LCmEUW5HtCuusI5oJjx2zp3BkAWQsG5/0KGYKGFJIUysI5oJu+Rkug94iWFzodYfekJyctARmShgSaGErCOa8U1SEjGZ2cQpSEYGJCUFHZGJAlaxUQLWEc34KeuCzmTEgWQLsfHxkJgYdEgmClhJ4Tjk5FhHNOO/75vU4pJb4dsBPWDOHFccNcZnVlI4Rlu2QO/ebqhr64hm/LTsp2UsaQiV+z8JdkuqKSVWUjgGuR3RFi2yjmjGf0u3LKVqfFWa1mkadCgmilhSKIbQjmhnnWUd0UzpSPoxic6ndiZG7L+pKT3211aE/B3RFiywjmjGfz/v+5mVqSu5uNHFQYdioowlhUKkp8OAAdYRzQRj7oa5AFzc2JKCKV2WFAqwdKm7s2jUKOuIZoIxfc106lapS6vftAo6FBNlLCmECO2IduCAuwvQOqKZ0rY/Yz/TVk+jR/MeNoeCKXVWGeKxjmimrFj0zj+4b246PZudG3QoJgpFfVLIyYHRo2HwYFcimDwZevQIOioTrXTxYi7o8wQXZ0Hs4odhTmvrtGZKVVRXH23ZAl26uJFNExNhxQpLCCZY3//vZeKy1MY7MoGJ2qRgHdFMWbPr4C4ezZpJZpygsbFg4x2ZAPiaFESki4isFpG1IjKkgO0iIiO97d+KSGs/4wHriGbKpl3zPmLSTS1JS09j3eSXkWHDbLwjEwjf2hREJBYYBVwGpABLRWSaqq4M2e33QBPv0R4Y4/0bfsnJ/PDvJB6cnsjOnTD3siTO/2sicU06uHHqk5J+/VUWrcsdovy9KOXr35+xn10zp/Jl5iYuHTmdO7LgjorxxPU5F7pZMjDB8LOhuR2wVlXXA4jIJKAbEJoUugFvqKoCS0SkpojUV9WtYY0kOZmMCy+icdYh3okBESF2tpKZJDz5x5N49H/bqJClZMcAIsRmB70MsdmU6nJmnDD0mpoMnbor5L0o/TiyY0AF4kq4nOUds7jLGXHwpz9U4enp6cRnHfvrs7wyd1xO8ZdjgZNy4AqBWHUPMrNd0rASggmIn0nhFGBzyPMUjiwFFLTPKcBhSUFE7gLuAjj11FOPPZKkJOKyM4lRkBwQVWIAyVauXanEZyuxCjE5gLct2GVKfVmy4eYfKhOfvSvkvSj9OIK7fuGujXWomH2AWNVSObcAoqASg8TFulvhrB3BBMzPpFBQLb0exz6o6ivAKwBt27Y9YnuREhOJqVgRMjKIjY11DQhZWcTGx9Oq/5Ou23JGBjEh26JtOTY+nhb9Ho/a9yI2Pp6z7/5LINcv8fEwYgSkpf1ajWVMQPxMCilAw5DnDYCfjmOfkuvQwTXaFVaP3KJF8PXZZWE52t+LIK/fEoEpI8RV5/twYJE4YA1wCbAFWArcqKrfhexzFTAAuBJXtTRSVdsd7bht27bVZcuW+RKzMcZEKhH5QlXbFrWfbyUFVc0SkQHAx7g2tVdV9TsR6ettHwvMwCWEtUA60NuveIwxxhTN12EuVHUG7os/dN3YkGUF+vsZgzHGmOKLCToAY4wxZYclBWOMMXksKRhjjMljScEYY0we325J9YuIpAI/HufL6wLbwxhOeWDXHB3smqNDSa75NFWtV9RO5S4plISILCvOfbqRxK45Otg1R4fSuGarPjLGGJPHkoIxxpg80ZYUXgk6gADYNUcHu+bo4Ps1R1WbgjHGmKOLtpKCMcaYo7CkYIwxJk9EJwURiRWRr0TkA+95bRH5RER+8P6tFXSM4VbANfcUke9EJEdEIvL2vQKu+XkR+V5EvhWR90SkZtAxhlsB1zzMu96vRWSWiJwcdIzhlv+aQ9Y/JCIqInWDis0vBXzOQ0Vki/c5fy0iV4b7nBGdFID7gVUhz4cAc1S1CTDHex5p8l/zCuBaYH4w4ZSK/Nf8CXCuqibg5vR4NJCo/JX/mp9X1QRVbQV8ADwWTFi+yn/NiEhD4DJgUyAR+e+IawaGq2or7zGjoBeVRMQmBRFpAFwF/CdkdTfgdW/5deCa0o7LTwVds6quUtXVwUXlr0KueZaqZnlPl+Bm9IsYhVzznpBdTqCAaW3Ls0L+PwMMBwYTYdcLR71mX0VsUgBG4P5YckLWnaSqWwG8f08MIjAfFXTNka6oa+4DfFR64ZSKAq9ZRJ4Skc3ATUReSeGIaxaRq4EtqvpNYFH5q7C/7QFeVeGrflSBR2RSEJGuwC+q+kXQsZQWu+YCt/8ZyAImlmpgPjraNavqn1W1Ie56B5R6cD4p6JpFpArwZyIv+QFH/ZzHAGcArYCtwD/DfW5fZ14LUCfgaq8RphJQXUT+C2wTkfqqulVE6gO/BBpleBV4zap6c8Bx+anQaxaRXkBX4BKNrM44xfmc3wQ+BB4PIkAfHHHNwASgMfCNiICrIvxSRNqp6s+BRRo+RX7OIvJvXPtReKlqRD+AROADb/l5YIi3PAR4Luj4/L7mkHVJQNugYyulz7kLsBKoF3RcpXjNTULW3wu8G3R8fl9zvvUbgbpBx1cKn3P9kPWDgEnhPl+klhQK8yzwjojcjrtboWfA8fhORLoDLwL1gA9F5GtVvSLgsPz2ElAR+MT7FblEVfsGG5LvnhWRs3D1zz8CkX690eo5EWmFa1jfCNwd7hPYMBfGGGPyRGRDszHGmONjScEYY0weSwrGGGPyWFIwxhiTx5KCMcaYPJYUTMQRkeEiMjDk+cci8p+Q5/8UkQeKOMbiYpxnY0Ejc4pIooh0PMrrrhGRx7zloSLykLecVNyRbEWknojMLM6+xhwLSwomEi0GOgKISAxQFzgnZHtHYNHRDqCqhX6pF0Ni7vkLMRgYXYLjo6qpwFYR6VSS4xiTnyUFE4kW8euX8jm44cP3ikgtEakINAe+AhCRh0VkqTfA2BO5BxCRfd6/MSIy2puT4gMRmSEiPULOda+IfCkiy0WkmYg0wnUcG+SNd39+aGAi0hQ4pKrbC4n9ZhFZLCIrRKSd95oLQ8bP/0pEqnn7TsUNfmdM2ERbj2YTBVT1JxHJEpFTcckhGTgF6ADsBr5V1QwRuRxoArQDBJgmIheoaujcE9cCjYAWuFF1VwGvhmzfrqqtReQe4CFVvUNExgL7VPUfBYTXCfjyKOGfoKodReQC7zznAg8B/VV1kYhUBQ56+y4D/lbc98WY4rCSgolUuaWF3KSQHPI8t73gcu/xFe6LuhkuSYTqDExW1Rx1A63Ny7d9ivfvF7jkUZT6QOpRtr8F4CWm6t6scYuAF0TkPqCm/jpXxC9AxM2wZoJlScFEqtx2hRa46qMluJJCaHuCAM/or7NYnamq4/IdR4o4zyHv32yKV/I+gBv1sjD5x51RVX0WuAOoDCwRkWbetkre8YwJG0sKJlItwg2dvUNVs1V1B1ATlxiSvX0+Bvp4VTKIyCkikn/ipYXAH722hZNwjchF2QtUK2TbKuDMo7z2Oi+WzsBuVd0tImeo6nJV/Tuuyig3KTTFJTxjwsaSgolUy3F3HS3Jt253biOvqs7CzT2QLCLLgXc58sv8f0AK7sv3ZeAzXLvE0UwHuhfU0IybK/s88YZvLcBO73bYscDt3rqBXsPzN7iSQe5Mchfh5k0wJmxslFRjiiAiVVV1n4jUAT4HOmkJJnIRkX8B01V1dgnjmg90U9WdJTmOMaHs7iNjivaB1+AbDwwrSULwPA20L8kBRKQe8IIlBBNuVlIwxhiTx9oUjDHG5LGkYIwxJo8lBWOMMXksKRhjjMljScEYY0ye/wdwkWUg575QUgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "year_5_pdf = pd.DataFrame({'lbs': np.linspace(40, 45, num = 1000)})\n",
    "year_5_pdf['day'] = max_day\n",
    "year_5_pdf = year_5_pdf.join(data_girls_percentiles[['day', 'L', 'M', 'S']].set_index('day'), on = 'day')\n",
    "year_5_pdf['kg'] = lbs_to_kg(lb = year_5_pdf['lbs'])\n",
    "year_5_pdf['prior_z_score'] = year_5_pdf[['kg', 'L', 'M', 'S']].apply(axis = 1,\n",
    "                                                                    func = lambda in_row: \n",
    "                                                                    LMS_z_score(weight_kg = in_row['kg'],\n",
    "                                                                                L = in_row['L'],\n",
    "                                                                                M = in_row['M'], \n",
    "                                                                                S = in_row['S']))\n",
    "for colname in ['L', 'M', 'S', 'kg']:\n",
    "    del year_5_pdf[colname]\n",
    "\n",
    "year_5_pdf['likelihood_z_score'] = (year_5_pdf['prior_z_score'] - likelihood_mean) / likelihood_stdev\n",
    "year_5_pdf['prior_pdf'] = scipy.stats.norm.pdf(x = year_5_pdf['prior_z_score'], loc=0, scale=1)\n",
    "year_5_pdf['likelihood_pdf'] = scipy.stats.norm.pdf(x = year_5_pdf['likelihood_z_score'], loc=0, scale=1)\n",
    "year_5_pdf['posterior_pdf'] = year_5_pdf['prior_pdf'] * year_5_pdf['likelihood_pdf']\n",
    "year_5_pdf['prior_cdf'] = year_5_pdf['prior_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['likelihood_cdf'] = year_5_pdf['likelihood_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['posterior_cdf'] = year_5_pdf['posterior_pdf'].cumsum(axis = 0)\n",
    "year_5_pdf['prior_cdf'] = year_5_pdf['prior_cdf'] / np.max(year_5_pdf['prior_cdf'])\n",
    "year_5_pdf['likelihood_cdf'] = year_5_pdf['likelihood_cdf'] / np.max(year_5_pdf['likelihood_cdf'])\n",
    "year_5_pdf['posterior_cdf'] = year_5_pdf['posterior_cdf'] / np.max(year_5_pdf['posterior_cdf'])\n",
    "\n",
    "### check max posterior\n",
    "print(\"Max posterior: \"+ str(np.max(year_5_pdf['posterior_pdf'])))\n",
    "\n",
    "plt.plot(year_5_pdf['lbs'], year_5_pdf['prior_cdf'], label = 'prior_cdf',\n",
    "        linestyle = '-', color = 'b')\n",
    "plt.plot(year_5_pdf['lbs'], year_5_pdf['likelihood_cdf'], label = 'likelihood_cdf',\n",
    "        linestyle = '-', color = 'g')\n",
    "plt.plot(year_5_pdf['lbs'].iloc[lambda x: x.index % 10 == 0], \n",
    "         year_5_pdf['posterior_cdf'].iloc[lambda x: x.index % 10 == 0], \n",
    "         label = 'posterior_cdf',\n",
    "        linestyle = 'None', marker = '.', color = 'r')\n",
    "\n",
    "plt.xlabel('Weight (lbs)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Bayesian Baby Weight Estimate')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Start MCMC\n",
    "Approximate same posterior with MCMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "43.01801801801802"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_posterior_lbs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "iteration = [0]\n",
    "MCMC_lbs = [max_posterior_lbs]\n",
    "MCMC_day = max_day\n",
    "MCMC_L = data_girls_percentiles.loc[data_girls_percentiles['day'] == MCMC_day]['L'].iloc[0]\n",
    "MCMC_M = data_girls_percentiles.loc[data_girls_percentiles['day'] == MCMC_day]['M'].iloc[0]\n",
    "MCMC_S = data_girls_percentiles.loc[data_girls_percentiles['day'] == MCMC_day]['S'].iloc[0]\n",
    "MCMC_prior_z_score = [LMS_z_score(weight_kg = lbs_to_kg(lb = MCMC_lbs[0]), \n",
    "                                  L = MCMC_L, M = MCMC_M, S = MCMC_S)]\n",
    "MCMC_likelihood_z_score = [(MCMC_prior_z_score[0] - likelihood_mean) / likelihood_stdev]\n",
    "MCMC_prior_pdf = [scipy.stats.norm.pdf(x = MCMC_prior_z_score[0], loc=0, scale=1)]\n",
    "MCMC_likelihood_pdf = [scipy.stats.norm.pdf(x = MCMC_likelihood_z_score[0], loc=0, scale=1)]  \n",
    "MCMC_posterior_pdf = [MCMC_prior_pdf[0] * MCMC_likelihood_pdf[0]]\n",
    "MCMC_alpha = [MCMC_posterior_pdf]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, 30000):\n",
    "    MCMC_lbs.append(MCMC_lbs[i - 1] + (np.random.uniform(-1.0, 1.0) * 10))\n",
    "    MCMC_prior_z_score.append(LMS_z_score(weight_kg = lbs_to_kg(lb = MCMC_lbs[i]), \n",
    "                                          L = MCMC_L, M = MCMC_M, S = MCMC_S))\n",
    "    MCMC_likelihood_z_score.append((MCMC_prior_z_score[i] - likelihood_mean) / likelihood_stdev)\n",
    "    MCMC_prior_pdf.append(scipy.stats.norm.pdf(x = MCMC_prior_z_score[i], loc=0, scale=1))\n",
    "    MCMC_likelihood_pdf.append(scipy.stats.norm.pdf(x = MCMC_likelihood_z_score[i], loc=0, scale=1))\n",
    "    MCMC_posterior_pdf.append(MCMC_prior_pdf[i] * MCMC_likelihood_pdf[i])\n",
    "    if MCMC_posterior_pdf[i] == 0:\n",
    "        MCMC_alpha.append(0)\n",
    "    elif MCMC_posterior_pdf[i - 1] == 0:\n",
    "        MCMC_alpha.append(2)\n",
    "    else:\n",
    "        MCMC_alpha.append(MCMC_posterior_pdf[i] / MCMC_posterior_pdf[i - 1])\n",
    "    if MCMC_alpha[i] < 1:\n",
    "        if np.random.uniform(0, 1) > MCMC_alpha[i]:\n",
    "            MCMC_lbs[i] = MCMC_lbs[i - 1]\n",
    "            MCMC_prior_z_score[i] = MCMC_prior_z_score[i - 1]\n",
    "            MCMC_likelihood_z_score[i] = MCMC_likelihood_z_score[i - 1]\n",
    "            MCMC_prior_pdf[i] = MCMC_prior_pdf[i - 1]\n",
    "            MCMC_likelihood_pdf[i] = MCMC_likelihood_pdf[i - 1]\n",
    "            MCMC_posterior_pdf[i] = MCMC_posterior_pdf[i - 1]\n",
    "            MCMC_alpha[i] = MCMC_alpha[i - 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XeYFFXa9/HvzZCTStAVUEFRARXUd8S4imlFRFlfdcGAWR8DKyZWfXGNu7quPiuYVnkURVfBrIBgDpgVjCAGBNRRUAQflCAwcL9/nBqmGXtmmpmpqQ6/z3X1VbGr79MDdfepU3WOuTsiIiLVaZB0ACIikhuUMEREJCNKGCIikhElDBERyYgShoiIZEQJQ0REMqKEIXnJzNzMuiYdRyoz29zMlphZUQb7do7K0LA+Yqsmlt+b2WdJxyHJU8IoQGY218xWmlm7Cus/iE5SnVPW9TazSWb2v2a2yMzeMbOTom19ov0fq3CcXtH6l1PWmZmdY2bTzWypmZWY2cNmtkMlMb5sZr9GJ9jFZjalsn3jYmYNo8/vnbLu2KhsFdd9Wt3x3P1rd2/p7qvrILYrzOw/1ewz18yWR2Uoe92SwbHXSbbu/qq7b1vbmCv5rHvM7G9xHFvqnhJG4ZoDHF22EJ2Mm6XuYGa7Ay8CrwBdgbbAmcDBKbstAPYws7Yp604APq/weSOBocA5QBtgG+AJ4JAqYhzi7i2jz30ZuC+zotUNdy8F3gT2SVm9N/BpmnVT6jG09XFolKTKXkOSDkhymLvrVWAvYC5wKfBuyrobgOGAA52jda8Bt1ZxnD5ACXA7cHa0rihadxnwcrRua2A10Hs9YnwZODVluQewMmW5N+Fk/r/APOAWoHHKdickp9nAj8D1hB9ITYBFwA4p+24MLAfap4njr8CElOVPgBPTrDsumm8AXAx8CSwEHgLaRNs6R3E1jJa7EBLNL8DzwK3AfyrsewLwdVSG4dG2vsBKYBWwBPiwir/zAZVs60r4IbA4OvaD0fop0ecujY49sOzvXOG4w4CPov3uAjYBJqeUZaOU/R8G5kefNQXYLlp/elSGldFnTYjWdwAeJfwYmQOck/T/Gb3CSzWMwvUW0NrMukfX1AcCay9xmFlzYHfgkQyOdS9wfDR/EDAD+C5l+/6EE847NQnUzBoDx0Yxl1kNnAe0i+LcHzirwlsPB4qBnYEBwMnuvgIYBxyXst/RwPPuviDNx08B9jSzBtElvBaEJNA7ZV03ymsY5wB/JNRAOgA/ERJBOg8A7xBqUFcAg9PssxewbVS+y8ysu7s/DVxDOMm3dPdelRy/KlcDzwIbAZ2AmwHcfe9oe6/o2A9W8v4jgAMJNcVDCcni/xH+Hg0I30OZyYQfDRsD7wH3R581Kpr/Z/RZh5pZA2AC8CHQMSr3uWZ2UA3KKHVMCaOw3Uc40R9IuMzybcq2jQj/PuZVdxB3fwNoY2bbRse7t8IubTM5Tho3mdn/En59DgGuTPnMae7+lruXuvtc4A7WvUwEcJ27L3L3r4ERlF+CGwMcE52cIJyoK7vc9TbQHNgB+D3wmrsvI/zyLVv3VfQZAP9FqAmURMnpCuDIio3XZrY5sAtwmbuvdPfXgPFpPv9Kd1/u7h8STqLrmxyeiNqfyl6nRetXAVsAHdz91+jz18fN7v69u38LvAq87e7vR2V+HNipbEd3H+3uv6R8H73MbINKjrsLoaZ3VfS9zAb+Bxi0nvFJDJQwCtt9wDGESywVT/I/AWuATdfjWEOAfQknjFQL1+M4qc5x9w2BpkB/4BEz6wlgZtuY2UQzm29mPxN+cber8P5vUua/Ivzix93fJlxK2cfMuhEuz6Q7WePuvxJqAXtHr1ejTa+lrEttv9gCeLzsBA3MJNSGNqlw6A7Aoij5pIu3zPyU+WVAy3RxVuGP7r5hyut/ovV/AQx4x8xmmNnJ63nc71Pml6dZbglgZkVm9g8z+zL6O82N9qn4tyqzBdAhNckRai4Vvz9JgBJGAXP3rwi/lPsBj1XYtozQRnBEhoe7j3BJaFKFkyDAC0AnMyuuYZxr3P1VYBbwh2j1vwm1oq3dvTXhpGIV3rpZyvzmrHuZbAzhstRg4JEoMVRmCiEx/J7yhPFqyrrUhPENcHCFk3TT6Jd4qnmEWlnzSuKtTq26mXb3+e5+mrt3INSKbovpNuRjCJcDDwA2ILTNQPnfqmI5vgHmVPj+Wrl7vxhik/WkhCGnAPu5+9I02/4CnGhmw8rugopumR1XcUd3n0O4JDQ8zbYvgNuAsdGtuI3NrKmZDTKzizMJMrpjqwehfQSgFfAzsCSqJZyZ5m3DzGwjM9uMcIdW6vX4+whtHMfx29pVRVMINafNCA3cEGoYfYAdWTdh3A783cy2iOJub2YDKh4wStZTgSui72N3QltApr4HOqdcVlsvZnaUmXWKFn8inLjLbvf9HtiyJsdNoxWwglDLbE6oCaaq+FnvAD+b2UVm1iyqoWxvZrvUUTxSC0oYBc7dv3T3qZVsewPYL3rNNrNFwChgUiX7v+bu36XbRmgEvYXQAPy/hLuIDic0cFbmlrLnBwgn+EvdfXK07ULCr9dfCNe40zXOPglMAz4AniLczVMWawmhAdYprzVU5g3Cr+O33cNtPO6+kHAXzw9RQiwzknB561kz+4XQUL9rJcc9ltBgvxD4W1SGFdXEUubhaLrQzN6rYr8JFZ7DKLtcuAvwdvTdjgeGRkkfQjvDmOiS0J8yjKcy9xIuB35LSLZvVdh+F9Aj+qwnPDyjcighEc8h3MF1J+H7l4RZ9O9fpOCY2WjgO3e/NOlYAMzsQeBTd7886VhE0lHCkIJk4Wn2D4CdUn5Z13cMuxCeCZlDaJt5Atjd3d9PIh6R6uiSlBQcM7samA5cn1SyiPyO8IDiEuAm4EwlC8lmqmGIiEhGVMMQEZGMJN518vpq166dd+7cOekwRERyyrRp03509/a1OUbOJYzOnTszdWrau0BFRKQSZvZVbY+hS1IiIpIRJQwREcmIEoaIiGQk59ow0lm1ahUlJSX8+mtV/cdJTTRt2pROnTrRqFGjpEMRkYTlRcIoKSmhVatWdO7cGbOKHZZKTbk7CxcupKSkhC5duiQdjogkLLZLUmY22sx+MLPplWw3M7vJzGaZ2UdmtnNNP+vXX3+lbdu2ShZ1zMxo27atam4iAsTbhnEPYezhyhxMGLZxa8LYvv+uzYcpWcRD36uIlIntkpS7T4k6eKvMAODeqLvot8xsQzPb1N1rMpSniEjivvwS5s+HX3+F0lJYvbry1/z50KRJeN+aNZm99twT/vCHqmOIU5JtGB1Zd0jKkmjdbxKGmZ1OqIWw+eab10twIiKZuvpquPxyiLtrvosuKtyEke5aR9qv291HEQbuobi4OCt7SzQzjjvuOO677z4ASktL2XTTTdl1112ZOHEiAJMnT+avf/0rS5cuxd3p378/N9xwA1dccQVXXnklX3zxBV27hlEyb7zxRs4//3zeffddiouLWbJkCRdccAHPP/88TZs2pW3btlx//fXsumtlY/NUr2XLlixZsgSAYcOGMWnSJPr168f1119fy29DJHe9/z488QSsXBleq1bBihXw4Yfw3XehVrBiRahFlE1XrgzvHTIE+vWDFi2gYUMoKqr61aQJNG0KDRpU/8qGq8NJJowS1h3DuBPrjrmcU1q0aMH06dNZvnw5zZo147nnnqNjx45rt0+fPp0hQ4bw1FNP0a1bN0pLSxk1atTa7TvssAPjxo3j0kvDWD6PPPIIPXr0WLv91FNPpUuXLnzxxRc0aNCA2bNnM3PmzDqL/4477mDBggU0Kasji+SJl16Ce+6Bn36CmTOhTZtwSajsklFpafn8ihUhKQA0bhxejRqFk3+DBiERtG4NPXuGE33Zq0mTkCh69060qLFLMmGMB4ZE40PvCiyui/aLc8+FDz6odWzr2HFHGDGi+v0OPvhgnnrqKY488kjGjh3L0UcfzauvhtE///nPfzJ8+HC6desGQMOGDTnrrLPWvvePf/wjTz75JJdeeimzZ89mgw02WPvsw5dffsnbb7/N/fffT4MG4T6FLbfcki23rHzY5XvvvZcbbrgBM6Nnz57cd999zJkzh2OOOYbS0lL69i2/H+Gwww5j6dKl7LrrrlxyySUMHDhwvb8jkWx1883w1FPQtWs4+a9cCR06hCRQ9ioqKp9v3BgOPRQOOSTpyLNPbAnDzMYCfYB2ZlYCXA40AnD32wnjQvcDZgHLgJPiiqW+DBo0iKuuuor+/fvz0UcfcfLJJ69NGNOnT+eCCy6o9L2tW7dms802Y/r06Tz55JMMHDiQu+++G4AZM2aw4447UlRUlFEcM2bM4O9//zuvv/467dq1Y9GiRQAMHTqUM888k+OPP55bb7117f7jx4+nZcuWfFDXmVYkAe7w5pswbx489xw8/nj45f/220lHlvvivEvq6Gq2O3B2XX9uJjWBuPTs2ZO5c+cyduxY+vXrt97vHzRoEOPGjeOZZ57hhRdeWJsw1teLL77IkUceSbt27QBo06YNAK+//jqPPvooAIMHD+aiiy6q0fFFstWqVTBlChxwwLrrUyrzUgt58aR3NjnssMO48MILefnll1m4cOHa9dtttx3Tpk2jV69elb730EMPZdiwYRQXF9O6det13vvhhx+yZs2atZekquLulT4/oecqJFf9+CN89VV5m8OqVWE6fz4MHw7ffhvaIco88ADsvjv87nehnUFqT50P1rGTTz6Zyy67jB122GGd9cOGDeOaa67h888/B2DNmjX861//WmefZs2acd111zF8+PB11m+11VYUFxdz+eWXUzak7hdffMGTTz6ZNob999+fhx56aG3CKrskteeeezJu3DgA7r///lqWVCR+a9bAaafBRhtB+/ZQXAy77QZ77QX77gsHHgiDB8PXX0PLlnDVVeEW16eegoEDoXNnJYu6pBpGHevUqRNDhw79zfqePXsyYsQIjj76aJYtW4aZcUiaVrVBgwalPe6dd97JBRdcQNeuXWnevPna22rT2W677Rg+fDj77LMPRUVF7LTTTtxzzz2MHDmSY445hpEjR3LEEUfUrqAiMfj559DWsHx5uF116NBQg4CQGHbbDbp0Kb9zqezVrl1o1JZ4mcf9pEkdKy4u9ooj7s2cOZPu3bsnFFH+0/cr9WH58nDiX7Zs3fVNmsAPP4TbWaXmzGyauxfX5hi6JCUiiTvvPGjePCSLLl1g6lSYPh3mzAmJRMkiO+iSVA5buHAh+++//2/Wv/DCC7Rt2zaBiETW3+LFMGECbLklHHssXHihEkS2ypuEUdWdQfmqbdu2sT87kWuXLCW3zJ4NW20V5k84ITRaS/bKi0tSTZs2ZeHChTq51bGyAZSa6jYTicnll4fpmWfCjTcmG4tULy9qGJ06daKkpIQFCxYkHUreKRuiVaQu3XADPPwwvPNOWP7Xv3T7ay7Ii4TRqFEjDSEqkuXcw22zZ50VHqqD8DzFTTcpWeSKvEgYIpI9Pv88XGqaPj10yb10KfzyC1S8APDhh6HXV8kdShgiUieWLAkd/R1/fFhu2RI22SQ8bNeyZegFtlUr6NYNBgyADTZINl5Zf0oYIlJrX34ZEkFpaVj+05/gwQeTjUnqnhKGiNTaMceEZNGnDzz6aOj7SfKPEoaI1MrPP4e7nX73O3jxxewYSlTikRfPYYhIcq6+OkzPPlvJIt+phiEiNbJqFYwdG56pgDAmheQ3JQwRWS+ffw7jxsGtt4ZeZAHOPVe1i0KghCEiGXvpJdhvvzDfuDFssQVMmwbq67IwqA1DRDJ2ww0hUQwdCitWwNy5ShaFRAlDRKrlDqNHw6RJ0KsXjBiRdESSBCUMEanWmDFwyilhvn//ZGOR5KgNQ0QqNWkSPP98edfjixbpobxCpoQhIr/x2Wdw991w3XVhuXlzOP98JYtCp4QhIuuYPz/0CwXQvn1IHIcckmxMkh3UhiEi6zj99DC96KLwnIWShZRRwhCRtR54ACZMgCZN4OKLk45Gso0ShogA8MQTcOyxYf6TT2DDDZONR7KP2jBECtzSpfDFF3DUUWF51CjYcstkY5LspIQhUsC+/z4kh2XLwvLFF8NppyUbk2QvXZISKUBLlsB554UxLJYtgx49Qj9RV16ZdGSSzWJNGGbW18w+M7NZZvabJjQz28DMJpjZh2Y2w8xOijMeEQkP351xRnn3HldfDR9/HEbLa9w40dAky8V2ScrMioBbgQOBEuBdMxvv7p+k7HY28Im7H2pm7YHPzOx+d18ZV1wihWz5cth889BuAbBwIbRpk2xMkjvirGH0Bma5++woAYwDBlTYx4FWZmZAS2ARUBpjTCIFbcyYkCzOPBPmzVOykPUTZ8LoCHyTslwSrUt1C9Ad+A74GBjq7msqHsjMTjezqWY2dcGCBXHFK5L3HnssTM8/P7RfiKyPOBNGuvG3vMLyQcAHQAdgR+AWM2v9mze5j3L3Yncvbt++fd1HKlJAdt8dunZNOgrJRXEmjBJgs5TlToSaRKqTgMc8mAXMAbrFGJNIwVqwAJ57DoqKko5EclWcCeNdYGsz62JmjYFBwPgK+3wN7A9gZpsA2wKzY4xJpOD89FOoVWy8cVg+/PBk45HcFdtdUu5eamZDgGeAImC0u88wszOi7bcDVwP3mNnHhEtYF7n7j3HFJFJoFi0KI+SVlECXLnDZZXDiiUlHJbkq1ie93X0SMKnCuttT5r8D/hBnDCKFauVKuPbakCw6d4b334cNNkg6Ksll6hpEJM+UlsJWW8HXX4fljTaCOXOSjUnygxKGSJ65776QLLp2hXPOgZ13TjoiyRdKGCJ5ZPp0OPnkMH/ffbDbbsnGI/lFnQ+K5JEBUV8Kd9yhZCF1TwlDJE+sXAmzZ8MJJ5QPsypSl5QwRPLE6tVh2r17snFI/lLCEMkTN90Ups2bJxuH5C8lDJE88OOPYbQ8gMGDk41F8pcShkiOu+ceKOuT889/hg03TDQcyWNKGCI57I034KRonMozzoCRI5ONR/KbnsMQyVGrV8Pee4f5yZOhb99k45H8pxqGSI5avDgkjTPOULKQ+qGEIZKjzj8/THfaKdk4pHAoYYjkoFmzwvjcAH/6U7KxSOFQwhDJQS+8EKZPP627oqT+KGGI5KCSkjBVf1FSn5QwRHLQvHnQqpUGRJL6pYQhkqNat046Aik0ShgiIpIRJQyRHLNmDdx1F5glHYkUGiUMkRzz6adhWlycbBxSeJQwRHLMvfeG6VlnJRuHFB4lDJEccu21cN11YX777ZONRQqPOh8UyQHucPbZ8O9/h+U5c2DTTZONSQqPahgiWeyVV6BPH2jQoDxZPP88dO6cZFRSqFTDEMlS7rDvvmHasGHolXbECCgqSjoyKVSqYYhkqalTQ7K46CJYtQpuvlnJQpKlhCGShRYvDpeiAPbfP9FQRNZSwhDJMsuXw2WXwbJlcOmlcOCBSUckEqgNQyTL/P73MG0atGgBV12VdDQi5WKtYZhZXzP7zMxmmdnFlezTx8w+MLMZZvZKnPGI5IK5c6FrV3j/fXX/IdklthqGmRUBtwIHAiXAu2Y23t0/SdlnQ+A2oK+7f21mG8cVj0gu+OorWLgQBg2CrbdOOhqRdcVZw+gNzHL32e6+EhgHDKiwzzHAY+7+NYC7/xBjPCJZbeHC8ucryhq8RbJJnAmjI/BNynJJtC7VNsBGZvaymU0zs+PTHcjMTjezqWY2dcGCBTGFK5Ksb78N00GD4Mgjk41FJJ04E0a6q69eYbkh8H+AQ4CDgL+a2Ta/eZP7KHcvdvfi9u3b132kIlnAo/8dRx2VbBwilYnzLqkSYLOU5U7Ad2n2+dHdlwJLzWwK0Av4PMa4RLLSgw+GqR7Ok2wVZw3jXWBrM+tiZo2BQcD4Cvs8CfzezBqaWXNgV2BmjDGJZK3Ro8NU7ReSrWKrYbh7qZkNAZ4BioDR7j7DzM6Itt/u7jPN7GngI2ANcKe7T48rJpFs1bMnfP89DBsGG2yQdDQi6Zl7xWaF7FZcXOxTp05NOgyROjNjRhjbYostwnyLFklHJPnIzKa5e63GaVTXICIJuv768oGQRo9WspDspoQhkpAZM+AvfwnzI0bAfvslG49IddSXlEhCRowI07vvhhNPTDQUkYyohiGSkIkTw/SII5KNQyRTShgiCZgxA+bPh759oVWrpKMRyYwShkgCvvoqTE87Ldk4RNaHEoZIAubPD9OyO6REckGVCcPM7kmZPyH2aEQKTNOmSUcgkrnqahi9UuaHxhmISL5zDyPpDR0Kp5wS1jXUfYqSQ6r755pbj4GLZLHzzy+/lRbgyith002Ti0dkfVWXMDqZ2U2ErsrL5tdy93Nii0wkj6xeXZ4sXn8ddthBd0dJ7qkuYQxLmVcHTiI1sHo17LhjmL/8cthjj2TjEampKhOGu4+pr0BE8tHSpTBmDEyfDu3awTmqk0sOq/a2WjM7wczeM7Ol0WtqZUOpikiwYAH8139By5Zw9tlh3Zgx0KZNsnGJ1EaVNYwoMZwLnA+8R2jL2Bm43sxw93vjD1Ekd/zyC9x1F5x3Xvm6a64JSaN16+TiEqkL1bVhnAUc7u5zU9a9aGZHAOMAJQwRYM0a+Pxz6N69fN0VV8Bll4GlG91eJAdVlzBaV0gWALj7XDPT7yURoLQUttkG5swJywcdFC4/bbJJsnGJ1LXqEsbyGm4TKRjz54dkcdRRcOyxcNhhqlVIfqouYXQ3s4/SrDdgyxjiEck5Y8eG6bHHwoABycYiEqfqEkYvYBPgmwrrtwC+iyUikRyzPKpr9++fbBwicavuttobgZ/d/avUF7As2iZS8F5+GRo1gqKipCMRiVd1CaOzu//mkpS7TwU6xxKRSA6ZMwdeekkN3FIYqksYVXW+3KwuAxHJReeeG6Y335xsHCL1obqE8a6Z/WZMMDM7BZgWT0giuWHFChg/Psyr/UIKQXWN3ucCj5vZsZQniGKgMXB4nIGJZLvPPw/Tq67SuBZSGKrrfPB7YA8z2xcoG0zyKXd/MfbIRLLcAw+E6V57JRuHSH3J6HeRu78EvBRzLCI546OP4B//CP1D9emTdDQi9aPa3mpFZF0PPgi9osGLb7xRT3VL4VDCEFkPkyfDoEFh/rbb4OSTk41HpD6pqU4kQ3PmQL9+Yf7NN2G33ZKNR6S+qYYhkqGyZy0uuUTJQgpTrAnDzPqa2WdmNsvMLq5iv13MbLWZHRlnPCK1MS26sXzYsKr3E8lXsSUMMysCbgUOBnoAR5tZj0r2uw54Jq5YROrC3LlhrIuNNko6EpFkxFnD6A3McvfZ7r6SMEJfus6f/ww8CvwQYywitfLQQ/D117ClOvWXAhZnwujIut2il0Tr1jKzjoQnxm+v6kBmdrqZTTWzqQsWLKjzQEUq4x7G4x44MCyffXay8YgkKc6Eke7udK+wPAK4yN1XV3Ugdx/l7sXuXty+ffs6C1CkKj/9BIMHh9tnAZ56CrbbLtmYRJIU5221JcBmKcud+O2gS8XAOAtPPrUD+plZqbs/EWNcItVaswbatg01DIBvv4UOHZKNSSRpcdYw3gW2NrMuZtYYGASMT93B3bu4e2d37ww8ApylZCHZYMSIkCz22w+WLFGyEIEYaxjuXmpmQwh3PxUBo919hpmdEW2vst1CJElTpoTpo49CixbJxiKSLWJ90tvdJwGTKqxLmyjc/cQ4YxFZH59+Cm3awIYbJh2JSPbQk94iFbz5Jnz2GRxwQNKRiGQXJQyRFO+8A3vsEeaLi5ONRSTbqPNBkRTfRE8O3X03nHhioqGIZB3VMERSlJSE6X77JRuHSDZSwhBJ8euvYdquXbJxiGQjJQyRFDNnhmnjxsnGIZKNlDBECA/nnXQSjBkTnrsoKko6IpHso4QhBc8dnngC7rknjNX95JMap1skHd0lJQWvT5/yJ7tvv12j6YlURglDCtrTT5cni48/hu23TzYekWymS1JSkNxDt+UHHxyWJ09WshCpjhKGFKQPPywfDGniROjbN9l4RHKBEoYUpCVLwvQ//4FDDkk2FpFcoYQhBW2TTZKOQCR3KGFIQbrqqqQjEMk9ShhScMaMgeeeC/O9eycbi0guUcKQgjJvXnkvtN99B61bJxqOSE5RwpCC8uKLYXrBBbDppsnGIpJrlDCkoNx1V5hefnmycYjkIiUMKRi//AIvvQTbbQetWiUdjUjuUcKQgnH11WHav3+ycYjkKiUMKQiLF8PDD4f5a69NNhaRXKWEIQXhqqtg7lzYcUd1XS5SU0oYkvfc4csvw6BIb76ZdDQiuUvdm0veKy6G994LD+k1bZp0NCK5SzUMyVslJXDQQSFZNGsGd9yRdEQiuU01DMlLM2dCjx5hvm1bmDoVOndONCSRnKcahuSdW24pTxbHHw8//qhkIVIXVMOQvDN8eJhOmqSBkUTqkhKG5JXZs+Hnn+Hii8uHXxWRuhHrJSkz62tmn5nZLDO7OM32Y83so+j1hpn1ijMeyW9ffAFbbRXmlSxE6l5sCcPMioBbgYOBHsDRZtajwm5zgH3cvSdwNTAqrngk/40cGabnnAN7751sLCL5KM4aRm9glrvPdveVwDhgQOoO7v6Gu/8ULb4FdIoxHslzkybBttuWJw4RqVtxJoyOwDcpyyXRusqcAkxOt8HMTjezqWY2dcGCBXUYouSDpUvhk0/CgEjbbJN0NCL5K85G73Q99njaHc32JSSMvdJtd/dRRJeriouL0x5DClNpaRgI6ZdfwvIuuyQbj0g+izNhlACbpSx3Ar6ruJOZ9QTuBA5294UxxiN56IILQrLYcku47Ta1XYjEKc6E8S6wtZl1Ab4FBgHHpO5gZpsDjwGD3f3zGGORPPXAA2H66qvQoUOysYjku9gShruXmtkQ4BmgCBjt7jPM7Ixo++3AZUBb4DYLfU6XuntxXDFJfjnhhPAU96BBShYi9cHcc6tJoLi42KdOnZp0GJKwu+6CU08N8z/+GPqLEpHKmdm02v4gV19SknNKS8uTxYQJShYi9UUJQ3LOtGlh+s9/anxukfqkhCE5ZcIE2G23MN+xqqd6RKTOqfNBySllT3GPGweqPyF4AAAMdElEQVQDByYbi0ihUQ1DcsZLL8ELL4RR9JQsROqfEobkhMmTYb/9wvyhhyYbi0ih0iUpyXpz5kC/fmH+449h++2TjUekUKmGIVlr9WqYODF0+wFhuFUlC5HkqIYhWWPpUnjssTBdsQL++tfyTgUPPzw8rCciyVHCkESVlsKCBTB4cGjQrqhXLxgzJkxFJFlKGJKoAw6AV14pX775ZjjiCGjaFJo0gWbNwNJ1lC8i9U4JQxLx9dfQuzd8/31IDHffHZJH+/ZJRyYilVHCkHq3ahVssUWY33DD0NVHWcO2iGQv3SUl9W7mzDAdOBAWLVKyEMkVShgSO3dYsgSuvx523bW8AfvMM9U+IZJLdElKYrFoETz8MIwYAZ9+uu62rl3h5JNhn32SiU1EakYJQ2Jx6aXw73+H+Q4dYOhQ2GCD8PBds2bJxiYiNaOEIbGYMQM23jjcDdWkSdLRiEhdUBuG1Lm33oIpU2DnnZUsRPKJEobUqenTYffdw/yf/pRsLCJSt5QwpNb+9jfo3h022wx22CGsGzUKTjop2bhEpG6pDUNqZOlS6NEjtFGUGTwYWrWCbt3gtNOSi01E4qGEIRmbORPGjoXXXguj3wF06gRnnAGnnAK/+12y8YlIvJQwpFrLloUawwMPhOX27cPlpyFDYNgwPXwnUiiUMKRKK1bAyJHlyeKVV2DvvZONSUSSoYQhaa1aFdokHnywfN2yZXroTqSQKWEUiOXL4ddfw7CnpaXlr9WrQ3J49ln46afQmL10Kdx2W/l7zz0XTjhByUKk0Clh5Im33oLXXw8n/bvvLk8IK1eGjv8y1bQptGgBm28OPXuGIVMbNYovbhHJHUoYeWD+/PKH5SA8Xd28ORxzTJhv3Dgkkm7doGHD376KisJ+ffuG94mIpKOEkcPc4bjjyhukzz8frrsuJAERkbqmU0sWKS0NT00vXhzaFVauDNMvvwwd+a1aFfYpm775ZmiXAHjoITjySN3iKiLxiTVhmFlfYCRQBNzp7v+osN2i7f2AZcCJ7v5enDFlA3cYPRrmzQsn/7Lk8N//Xb5PmzbhUlKjRqHG8N57YRyJsuVGjcJlqAYN4N57QxfiIiJxii1hmFkRcCtwIFACvGtm4939k5TdDga2jl67Av+OpjnDPbPX44+HZxiWL4dPPgkJoEzDhiE5tGgR+mKaMkUNzSKSfeKsYfQGZrn7bAAzGwcMAFITxgDgXnd34C0z29DMNnX3eXUdzEcfwSGHwM8/h0beTE/0Vb1qYsstQwNz797hbqbu3XUZSURyQ5wJoyPwTcpyCb+tPaTbpyOwTsIws9OB0wE233zzGgWzdGnoGG/nncMto2b1/zryyJAwRERyUZwJI93v5oq/yzPZB3cfBYwCKC4urtFv+913D5eCRESkZuIcD6ME2CxluRPwXQ32ERGRLBBnwngX2NrMuphZY2AQML7CPuOB4y3YDVgcR/uFiIjUXmyXpNy91MyGAM8Qbqsd7e4zzOyMaPvtwCTCLbWzCLfVaow2EZEsFetzGO4+iZAUUtfdnjLvwNlxxiAiInVDY3qLiEhGlDBERCQjShgiIpIRJQwREcmIeU37uEiImS0Avqqnj2sH/FhPnxU3lSU7qSzZK5/K0w5o4e7ta3OQnEsY9cnMprp7cdJx1AWVJTupLNkrn8pTV2XRJSkREcmIEoaIiGRECaNqo5IOoA6pLNlJZcle+VSeOimL2jBERCQjqmGIiEhGlDBERCQjBZswzKzIzN43s4nR8vVm9qmZfWRmj5vZhmnes62ZfZDy+tnMzq3/6H8T13qXJdrvPDObYWbTzWysmTWt38jTxlTTsgyNyjEjG/4mZdKU5+qoLB+Y2bNm1qGS9/U1s8/MbJaZXVy/UadXi7KMNrMfzGx6/UZcuZqUxcw2M7OXzGxm9O9saP1H/ls1LEtTM3vHzD6MynJlJp9VsAkDGArMTFl+Dtje3XsCnwOXVHyDu3/m7ju6+47A/yF0yf54fQRbjfUui5l1BM4Bit19e0IX9IPqIdbq1KQs2wOnEcaR7wX0N7Ot6yHWTFQsz/Xu3jP6NzQRuKziG8ysCLgVOBjoARxtZj3qI9hqrHdZIvcAfWOObX3VpCylwAXu3h3YDTg7h/8uK4D93L0XsCPQNxqTqEoFmTDMrBNwCHBn2Tp3f9bdS6PFtwij/1Vlf+BLd6+vp87TqmVZGgLNzKwh0JyERzusRVm6A2+5+7Jo31eAw+OOtzqVlOfnlF1akGZIYkLim+Xus919JTAOGBBnrNWpRVlw9ynAolgDXA81LYu7z3P396L5Xwgn6Y7xRlu1WpTF3X1JtNgoelV7B1RBJgxgBPAXYE0l208GJldzjEHA2LoMqoZqVBZ3/xa4AfgamEcY7fDZuILMUE3/LtOBvc2srZk1JwzKtVma/epb2vKY2d/N7BvgWNL/+usIfJOyXELCJyZqXpZsVOuymFlnYCfg7XhCzFiNyxJdyvoA+AF4zt2rLUvBJQwz6w/84O7TKtk+nFD1vL+KYzQGDgMejiXIDNWmLGa2EeFXaxegA9DCzI6LMdwq1aYs7j4TuI5w+epp4MNo38RUVR53H+7umxHKMiTd29OsS+z+91qWJavURVnMrCXwKHBuhV/z9aq2ZXH31dFlq05A7+jSbtXcvaBewLWEX2xzgfmEdoj/RNtOAN4EmldzjAHAs7lcFuAo4K6U5eOB23KxLGmOdQ1wVrb+bVL22QKYnua9uwPPpCxfAlySi2VJ2d65qu25UhbCpZtngPNzvSwV9rscuLDa/ZIudMJfeB9gYjTfF/gEaJ/B+8YBJyUdf23KAuwKzCC0XRgwBvhz0uWo6d8F2Diabg58CmyUdDkqKc/WKev/DDySZv+GwGxC7a8xoca0XdLlqElZUrZnRcKo5d/FgHuBEUnHXgdlaQ9sGM03A14F+lf3OQV3SaoKtwCtgOei29FuBzCzDma2dlzy6Br5gcBjyYSZkWrL4uF65SPAe8DHhMuT2dgVQkZ/F+BRM/sEmACc7e4/JRBrJv4R3f77EfAHwh0uFf82pYTLCM8QGlYfcvcZSQVchWrLEi2PJdQQtzWzEjM7JZlwq5RJWfYEBgP7Wfmt9f0SircqmZRlU+ClaJ93CW0YE6s7sLoGERGRjKiGISIiGVHCEBGRjChhiIhIRpQwREQkI0oYIiKSESUMyTtmdqOl9FhrZs+Y2Z0py/9tZudXc4w3MvicuWbWLs36Pma2RxXv+6OZXRbNX2FmF0bzL5tZcXWfG+3b3syezmRfkbqihCH56A1gDwAzawC0A7ZL2b4H8HpVB3D3Sk/4GehT9vmV+AtwWy2Oj7svAOaZ2Z61OY7I+lDCkHz0OuUn7O0InRP+YmYbmVkTQu+27wOY2TAzezcaP2DtmABmtiSaNjCz26IxAyaa2SQzOzLls/5sZu+Z2cdm1i3qlO4M4Lzowa7fpwZmZtsAK9z9x0piP87M3ogevOodvWeflAfF3jezVtG+TxA6lxOpFw2TDkCkrrn7d2ZWamabExLHm4TeXncHFgMfuftKM/sDsDWhO3EDxpvZ3h664y7zfwndWuwAbEx48np0yvYf3X1nMzuL0BfPqdHT6Evc/YY04e1JeLq+Mi3cfQ8z2zv6nO2BCwlPr78edXz3a7TvVOBvmX4vIrWlGobkq7JaRlnCeDNluax94g/R633CSbwbIYGk2gt42N3XuPt84KUK28u6iJlGSCzV2RRYUMX2sbB2DInWFkYYfB34l5mdQ+j/p6wn3h8IPQ2L1AslDMlXZe0YOxAuSb1FqGGktl8YcK1Hoyi6e1d3v6vCcdJ1NZ5qRTRdTWY19uVAVUPhVuyrx939H8CphE7i3jKzbtG2ptHxROqFEobkq9eB/sAiD/3+LwI2JCSNN6N9ngFOji7zYGYdzWzjCsd5DTgiasvYhNCgXZ1fCB0mpjMT6FrFewdGsexFGNRqsZlt5e4fu/t1hMtQZQljG0IyFKkXShiSrz4m3B31VoV1i8sanD2MMPgA8KaZfUzovbfiif5RwpgD04E7CCOsLa7msycAh6dr9AamADuZWWU1l5+iW3pvB8p6dT03agT/kFCjKBt1cF/gqWpiEakz6q1WpBpm1tLdl5hZW+AdYM+oPaOmxxsJTHD352sZ1xRgQBZ35S55RndJiVRvYtT43Bi4ujbJInINYQCrGjOz9sC/lCykPqmGISIiGVEbhoiIZEQJQ0REMqKEISIiGVHCEBGRjChhiIhIRv4/jhc/PoHPbpoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "MCMC_results = pd.DataFrame({'MCMC_lbs': MCMC_lbs}).sort_values('MCMC_lbs')\n",
    "MCMC_results['cdf'] = list(range(MCMC_results.shape[0]))\n",
    "MCMC_results['cdf'] = (MCMC_results['cdf'] + 1) / MCMC_results.shape[0]\n",
    "#MCMC_results\n",
    "plt.plot(MCMC_results['MCMC_lbs'], MCMC_results['cdf'], label = 'MCMC_cdf',\n",
    "        linestyle = '-', color = 'b')\n",
    "plt.xlabel('Weight (lbs)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('MCMC Baby Weight Estimate')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## At 10 years old"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "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>Month</th>\n",
       "      <th>L</th>\n",
       "      <th>M</th>\n",
       "      <th>S</th>\n",
       "      <th>SD4neg</th>\n",
       "      <th>SD3neg</th>\n",
       "      <th>SD2neg</th>\n",
       "      <th>SD1neg</th>\n",
       "      <th>SD0</th>\n",
       "      <th>SD1</th>\n",
       "      <th>SD2</th>\n",
       "      <th>SD3</th>\n",
       "      <th>SD4</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>61</td>\n",
       "      <td>-0.4681</td>\n",
       "      <td>18.2579</td>\n",
       "      <td>0.14295</td>\n",
       "      <td>10.742</td>\n",
       "      <td>12.352</td>\n",
       "      <td>13.961</td>\n",
       "      <td>15.899</td>\n",
       "      <td>18.258</td>\n",
       "      <td>21.169</td>\n",
       "      <td>24.817</td>\n",
       "      <td>29.468</td>\n",
       "      <td>34.118</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>62</td>\n",
       "      <td>-0.4711</td>\n",
       "      <td>18.4329</td>\n",
       "      <td>0.14350</td>\n",
       "      <td>10.828</td>\n",
       "      <td>12.456</td>\n",
       "      <td>14.083</td>\n",
       "      <td>16.043</td>\n",
       "      <td>18.433</td>\n",
       "      <td>21.386</td>\n",
       "      <td>25.090</td>\n",
       "      <td>29.823</td>\n",
       "      <td>34.555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>63</td>\n",
       "      <td>-0.4742</td>\n",
       "      <td>18.6073</td>\n",
       "      <td>0.14404</td>\n",
       "      <td>10.915</td>\n",
       "      <td>12.559</td>\n",
       "      <td>14.204</td>\n",
       "      <td>16.187</td>\n",
       "      <td>18.607</td>\n",
       "      <td>21.601</td>\n",
       "      <td>25.363</td>\n",
       "      <td>30.178</td>\n",
       "      <td>34.992</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>64</td>\n",
       "      <td>-0.4773</td>\n",
       "      <td>18.7811</td>\n",
       "      <td>0.14459</td>\n",
       "      <td>11.000</td>\n",
       "      <td>12.662</td>\n",
       "      <td>14.324</td>\n",
       "      <td>16.330</td>\n",
       "      <td>18.781</td>\n",
       "      <td>21.817</td>\n",
       "      <td>25.637</td>\n",
       "      <td>30.535</td>\n",
       "      <td>35.433</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>65</td>\n",
       "      <td>-0.4803</td>\n",
       "      <td>18.9545</td>\n",
       "      <td>0.14514</td>\n",
       "      <td>11.085</td>\n",
       "      <td>12.764</td>\n",
       "      <td>14.444</td>\n",
       "      <td>16.473</td>\n",
       "      <td>18.954</td>\n",
       "      <td>22.032</td>\n",
       "      <td>25.911</td>\n",
       "      <td>30.893</td>\n",
       "      <td>35.876</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Month       L        M        S  SD4neg  SD3neg  SD2neg  SD1neg     SD0  \\\n",
       "0     61 -0.4681  18.2579  0.14295  10.742  12.352  13.961  15.899  18.258   \n",
       "1     62 -0.4711  18.4329  0.14350  10.828  12.456  14.083  16.043  18.433   \n",
       "2     63 -0.4742  18.6073  0.14404  10.915  12.559  14.204  16.187  18.607   \n",
       "3     64 -0.4773  18.7811  0.14459  11.000  12.662  14.324  16.330  18.781   \n",
       "4     65 -0.4803  18.9545  0.14514  11.085  12.764  14.444  16.473  18.954   \n",
       "\n",
       "      SD1     SD2     SD3     SD4  \n",
       "0  21.169  24.817  29.468  34.118  \n",
       "1  21.386  25.090  29.823  34.555  \n",
       "2  21.601  25.363  30.178  34.992  \n",
       "3  21.817  25.637  30.535  35.433  \n",
       "4  22.032  25.911  30.893  35.876  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_girls_5_to_10 = pd.read_csv('https://www.who.int/growthref/wfa_girls_z_WHO2007_exp.txt', sep=\"\\t\") #, header=None)\n",
    "data_girls_5_to_10.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_month = np.max(data_girls_5_to_10['Month'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max posterior: 0.1450449479137154\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>lbs</th>\n",
       "      <th>Month</th>\n",
       "      <th>prior_z_score</th>\n",
       "      <th>likelihood_z_score</th>\n",
       "      <th>prior_pdf</th>\n",
       "      <th>likelihood_pdf</th>\n",
       "      <th>posterior_pdf</th>\n",
       "      <th>prior_cdf</th>\n",
       "      <th>likelihood_cdf</th>\n",
       "      <th>posterior_cdf</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>40.000000</td>\n",
       "      <td>120</td>\n",
       "      <td>-3.874752</td>\n",
       "      <td>-307.541119</td>\n",
       "      <td>0.000219</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>40.100111</td>\n",
       "      <td>120</td>\n",
       "      <td>-3.854516</td>\n",
       "      <td>-306.081230</td>\n",
       "      <td>0.000237</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000004</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>40.200222</td>\n",
       "      <td>120</td>\n",
       "      <td>-3.834360</td>\n",
       "      <td>-304.627146</td>\n",
       "      <td>0.000256</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000006</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>40.300334</td>\n",
       "      <td>120</td>\n",
       "      <td>-3.814284</td>\n",
       "      <td>-303.178828</td>\n",
       "      <td>0.000276</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000008</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>40.400445</td>\n",
       "      <td>120</td>\n",
       "      <td>-3.794288</td>\n",
       "      <td>-301.736241</td>\n",
       "      <td>0.000298</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000010</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         lbs  Month  prior_z_score  likelihood_z_score  prior_pdf  \\\n",
       "0  40.000000    120      -3.874752         -307.541119   0.000219   \n",
       "1  40.100111    120      -3.854516         -306.081230   0.000237   \n",
       "2  40.200222    120      -3.834360         -304.627146   0.000256   \n",
       "3  40.300334    120      -3.814284         -303.178828   0.000276   \n",
       "4  40.400445    120      -3.794288         -301.736241   0.000298   \n",
       "\n",
       "   likelihood_pdf  posterior_pdf  prior_cdf  likelihood_cdf  posterior_cdf  \n",
       "0             0.0            0.0   0.000002             0.0            0.0  \n",
       "1             0.0            0.0   0.000004             0.0            0.0  \n",
       "2             0.0            0.0   0.000006             0.0            0.0  \n",
       "3             0.0            0.0   0.000008             0.0            0.0  \n",
       "4             0.0            0.0   0.000010             0.0            0.0  "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "year_10_pdf = pd.DataFrame({'lbs': np.linspace(40, 130, num = 900)})\n",
    "year_10_pdf['Month'] = max_month\n",
    "year_10_pdf = year_10_pdf.join(data_girls_5_to_10[['Month', 'L', 'M', 'S']].set_index('Month'), on = 'Month')\n",
    "year_10_pdf['kg'] = lbs_to_kg(lb = year_10_pdf['lbs'])\n",
    "year_10_pdf['prior_z_score'] = year_10_pdf[['kg', 'L', 'M', 'S']].apply(axis = 1,\n",
    "                                                                        func = lambda in_row: \n",
    "                                                                        LMS_z_score(weight_kg = in_row['kg'],\n",
    "                                                                                    L = in_row['L'],\n",
    "                                                                                    M = in_row['M'], \n",
    "                                                                                    S = in_row['S']))\n",
    "for colname in ['L', 'M', 'S', 'kg']:\n",
    "    del year_10_pdf[colname]\n",
    "\n",
    "year_10_pdf['likelihood_z_score'] = (year_10_pdf['prior_z_score'] - likelihood_mean) / likelihood_stdev\n",
    "year_10_pdf['prior_pdf'] = scipy.stats.norm.pdf(x = year_10_pdf['prior_z_score'], loc=0, scale=1)\n",
    "year_10_pdf['likelihood_pdf'] = scipy.stats.norm.pdf(x = year_10_pdf['likelihood_z_score'], loc=0, scale=1)\n",
    "year_10_pdf['posterior_pdf'] = year_10_pdf['prior_pdf'] * year_10_pdf['likelihood_pdf']\n",
    "year_10_pdf['prior_cdf'] = year_10_pdf['prior_pdf'].cumsum(axis = 0)\n",
    "year_10_pdf['likelihood_cdf'] = year_10_pdf['likelihood_pdf'].cumsum(axis = 0)\n",
    "year_10_pdf['posterior_cdf'] = year_10_pdf['posterior_pdf'].cumsum(axis = 0)\n",
    "year_10_pdf['prior_cdf'] = year_10_pdf['prior_cdf'] / np.max(year_10_pdf['prior_cdf'])\n",
    "year_10_pdf['likelihood_cdf'] = year_10_pdf['likelihood_cdf'] / np.max(year_10_pdf['likelihood_cdf'])\n",
    "year_10_pdf['posterior_cdf'] = year_10_pdf['posterior_cdf'] / np.max(year_10_pdf['posterior_cdf'])\n",
    "\n",
    "### check max posterior\n",
    "print(\"Max posterior: \"+ str(np.max(year_10_pdf['posterior_pdf'])))\n",
    "year_10_pdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4FFXWwOHfScK+C7gBsoysCiIiyCJuDDqKoCiKomxuIKgwMm6jgsLnhgsgICICLijjAoiKqCCKKKJBZBeGAcSISgSUHZL0+f64ldCEzt6VTnfO+zz9pLqquup0pbtO33ur7hVVxRhjjAGIi3QAxhhjig5LCsYYYzJYUjDGGJPBkoIxxpgMlhSMMcZksKRgjDEmgyUFUyAicq6IrI90HJmJSB0RURFJiHQswfJyvETkfBFJ8jsmY4JZUogyIrJFRA6IyF4R2SUiH4pIrUjFo6pfqmpDP7ZdFN6riJzkJZcTgub9O4t583LaXjiPl4hME5GRuVhPRORaEVkgIttF5DcR+UREuuXitZNEZL2IBESkT4jlQ7zt/SUiU0SkVD7fjikiLClEp8tVtTxwEvA78HyE4/FTRN+rqv4KbAQ6BM3uAPwYYt6iQgwtV0QkHngDuAV4DPgbUAsYDtwqIi+JiGSziRXA7cD3IbZ9MXAfcBFQB6gHPBLG8POsqJUMo5ElhSimqgeBd4Am6fNE5DIRWS4iu0XkZxEZHrTsQxG5I3gbIrJSRK7wphuJyKcistP7dXhN0HqXishaEdkjIr+IyFBv/lFVHCJyn4j8z1tvrYhcGbSsj4gsFpGnvV/+m0XkH3681yD9RGSbiPwqInd7rztRRPaLSNWgbZ0lIskiUiLENhbhJQDvJHsmMCbTvDbeeohIKe89bhWR30VkooiUyeJ4tfDewx4ReVtE/pP517+I3O39wv9VRPp6824FegL3eCWp97M4dP8GDgF/V9UFqrpHVVNU9WvgH0BF4MYsXouqjlfVBcDBEIt7Ay+r6hpV3QWMAPqE2o6IjBeRZzLNe19EBnvTJ4vIu97/YLOI3Bm0XisRWSIif3rHYJyIlAxariIyUET+C/zXKxk95x2zv7zP+OlZvUeTiaraI4oewBagozddFngFeDVo+flAU1zCb4b7dX2Ft+waYGnQumcAO4CSQDngZ6AvkAC0AP4ATvPW/RU415uuArQI2l9S0Da7Ayd7+78W2Aec5C3rA6TgfrXGAwOAbYD48F7rAAq86b23pkBy0PbmAgOCtvUc8HwWcfQGVnjTLXEn//qZ5h0ASnrPRwNzgOOACsD7wOOZj5d33H8C7gJKAN2Aw8DIoHVTgUe95ZcC+4Eq3vJp6etmEXc5b/vlvNdPAX4D5gFTcUmtBpCYi8/dYqBPpnkrgGuDnlfzjnnVEK9v5f2v44LW3Q+c4P3/lgEPe8ekHrAJuNhb9yzgHNznsg6wDhgctG0FPvWOdxngYm97lQEBGuN9Bu2Ri3NMpAOwRx7/Ye5EuRf40zthbAOaZrP+aOA5b7oUsBOo7z1/GpjgTV8LfJnptS8Cw7zprcBtQMVM62Sc5LLY/w9AV2+6D7AxaFlZ7wt9og/vtY637UZBy5/C/bJNf79fedPx3smyVRbbrQOk4ZLhEOD/vPm/BM1b6M0TXCL8W9Dr2wCbMx8v3En5F4KSIu7kG5wUDgAJQcu3A+d409PIPil0BF7xpgcAn+BKBvW8Y3m+t2xzLj53oZLC/4BLgp6X8I55nSy2sQ5XYgEYBMz1plsDWzOtez8wNYvtDAZmBT1X4MKg5xcCG3CJJK4wvpex9LDqo+h0hapWxp3kBwFfiMiJACLSWkQWesXwv4D+uF9lqOoh4C3gBhGJA64DXvO2WRto7RXR/xSRP3HVEyd6y6/C/VL9SUS+EJE2oQITkV4i8kPQNk5P37/nt/QJVd3vTZYP93sN8nPQ9E+4UgzAe0ATEakH/B34S1W/DRWAqm4BkoD2uBP5l96iJUHz0tsTquOS3bKgYzDPm5/ZycAv6p3JQsQLsENVU4Oe7yf74xXseFzSAVdSmq2qu1V1E+4kj4hUwCWx/NiLSzLp0qf3ZLH+K8AN3vQNHP3ZOznTZ+8BXCkCEWkgIh+Ia9DejWsbyfL/rKqfAeOA8cDv4hrLK2JyxZJCFFPVNFWdifsV296b/Qau6qKWqlYCJuJ+vaZ7BXeyvwjYr6pLvPk/A1+oauWgR3lVHeDt6ztV7Yo70czGJZejiEht4CXcybuqdzJfnWn/hflewTWqpjsF9wsZdW0Ub+GOxY0cOUFl5Uvcyb8N8HWmee05khT+wP26Py3oOFZS11ie2a9ADZGjGnrzcnVVTl0c/4FroAdYBVwhIhVEpK4XcxVgAq5aKT/W4Kog050B/K6qO7JY/3Wgq4icgavSme3N/xlXWgn+7FVQ1Uu95S/gGvbrq2pFXMLI/H8+6lio6lhVPQs4DWgA/Ct/b7H4saQQxbwGta64L/c6b3YFYKeqHhSRVsD1wa/xkkAAeIajT4QfAA1E5EYRKeE9zhaRxiJSUkR6ikglVU0BduNOzpmVw305k734+uJKChF5r56HRKSsiJyGay/5T9CyV3FVWl1wJ6zsLAJ6AdtUdbc3b7E3rxKu1ICqBnCJ8TkROd6LvYa4K3UyW4I7joNEJMF7f61yiCPY77iqoKwsAS7wGrlfxpV2NnrTn+Dq8Jfgqt1C8v73pXEn4RIiUtorZYI7fjeJSBMRqQI8iKvSCklVk4DvcJ+7d1X1gLfoW2C3iNwrImVEJF5ETheRs73lFXCfub0i0ghXFZYl73PbWtxFA/twjeShPq8mlEjXX9kjbw9cPfsBXNF9D+6XeM+g5Vfjqkn24E7044DXM23jQdzJu16m+Q2BD3En9R3AZ0BzXOPfPGAX7sv5HdDee835HN3Q/H+4dos/gGeBL4CbvWV9gMWZ9qnAqeF+rxxpU7gVVzr4DbgnxD7+iysh5XTcG3rbGxs0L947HksyrVsaV8WxyVu+Drgzi+PVEtfushd4G5gJPBRq3aBjkt5YXt977Z+4qqFQcT8GvJjFsoRcvO/Pvfcd/Dg/aPk/cclpN67xulQO27vB28YFmeafjLso4Dfvc/ZN0PtMvwR4L6509mjw5yjzZwhXCl7prf8HMB0oH+nvbrQ8xDuIphgRkV7AraraPseVY5yIfAa8oaqTIx0LgIgsBSaq6tQwbS8Bl2zigJG4JFIWd5XYvbiryLJqAwg7EemAK5XVUVeqMkWMVR8VMyJSFncz0qRIxxJpXvVEC46uUirsGM4Td99Egoj0xl1am+Od0bmlrpH6KlzD+nO4X+I/4q7QubKQE0IJ3OW3ky0hFF12918x4tVrzwTm4xppiy0ReQW4ArirME+MITTENXiXx13iebW6u6jDxjsBTyH/DcoFJiKNgUTcvQ19IxWHyZlVHxljjMlg1UfGGGMyRF31UbVq1bROnTqRDsMYY6LKsmXL/lDVUDdRHiXqkkKdOnVITEyMdBjGGBNVROSn3Kxn1UfGGGMyWFIwxhiTwZKCMcaYDFHXphBKSkoKSUlJHDwYahwQE2mlS5emZs2alCgRavwaY0xREhNJISkpiQoVKlCnTh0k25EFTWFTVXbs2EFSUhJ169aNdDjGmBz4Vn0kbhDv7SKyOovlIiJjRWSjN1xei/zu6+DBg1StWtUSQhEkIlStWtVKccZECT/bFKYBl2Sz/B+4Xh7r43qyfKEgO7OEUHTZ/8aY6OFb9ZGqLhKROtms0hU33q4C34hIZRE5Kdz9vpjiY/uu7fR/qT+7DuwiLWDd55uiQQENCIFAHAGNQwNxbjrTQzV9WoLmifca97dj/dZM+fdAX+ONZJtCDY4eejDJm3dMUhCRW3GlCU455ZRCCc5En45PdmRVmVXuSRw5j0tmTL7lsfQb7z0K6OMNaUDsJoVQRzXk11hVJ+F19dyyZcuY/arPnj2bBg0a0KRJkzy9bs6cOaxdu5b77rvPp8iOmDZtGomJiYwbN47k5GQ6d+7M4cOHGTt2LOeee67v+8/Ojyk/csKhE/jlmV+IjwvDN9BEtQMHIDk568eOHfDnn0cef/3lHoFcdOpdujSULw/lyh39yDyvTBm3bqlSR/8NNS/U31KloGRJKFHCPQqjJjaSSSGJo8ejrYk3fm5xNXv2bDp37pynpJCamkqXLl3o0qVLnl6TkFDwf/2CBQto1KgRr7zySoG3VVCBQICUMinU0lqWEGLcoUOwbRv88suRR1LSkelt22D7dti3L/TrExKgWjX3qFIFatWCpk2hcmWoVMn9DX5UquQeFSq4E33ZshAfwx+xSCaFObixaWcArYG/wtGeMHjwYH744YcCBxesefPmjB6d5TC2AGzZsoVLLrmE1q1bs3z5cho0aMCrr77KkiVLGDp0KKmpqZx99tm88MILlCpVivvuu485c+aQkJBAp06d6NatG3PmzOGLL75g5MiRvPvuuwAMHDiQ5ORkypYty0svvUSjRo3o06cPxx13HMuXL6dFixY0bdo049f7Tz/9RL9+/UhOTqZ69epMnTqVU0455ZjXPPPMM8e8h71793LHHXeQmJiIiDBs2DCuuuoqpk6dyuOPP85JJ51EgwYNKFWqFD/88AP33HMPBw4coHnz5ixZsoQyZcqE9bjnxY8//wgloHbZ2hGLwYRHWpo7sW/aFPqxffuxrylbFmrUcI82beCEE6B69dCPSpUK5xd3tPItKYjIm7gxZquJSBIwDCgBoKoTgbnApbiBxPcTAwNvrF+/npdffpl27drRr18/nn32WV588UUWLFhAgwYN6NWrFy+88AK9evVi1qxZ/Pjjj4gIf/75J5UrV6ZLly507tyZq6++GoCLLrqIiRMnUr9+fZYuXcrtt9/OZ599BsCGDRuYP38+8fHxTJs2LSOGQYMG0atXL3r37s2UKVO48847mT179jGvCWXEiBFUqlSJVatcvfyuXbv49ddfGTZsGMuWLaNSpUpccMEFnHnmmTRv3pxHH300IxlF2v9+/R8ANSrXiHAkJrf274f162HdOli79sjfTZvg8OEj68XHwymnQL160KWLm65RA2rWPJII7EQfPn5efXRdDssVH1pMcvpF76datWrRrl07AG644QZGjBhB3bp1adCgAQC9e/dm/PjxDBo0iNKlS3PzzTdz2WWX0blz52O2tXfvXr7++mu6d++eMe/QoUMZ0927dw95cl+yZAkzZ84E4MYbb+See+7J8TXp5s+fz4wZMzKeV6lShdmzZ3P++edTvbrrcffaa69lw4YNuToehenPfX8CUKlMpQhHYjJTddU7y5fD99+7x+rVsGWLWwbuxH/qqdC4sTvx/+1vLgnUq+eqd+xm+MITE3c0FxW5vR4/ISGBb7/9lgULFjBjxgzGjRuXUQJIFwgEqFy5cpZVYeXKlctzTDm9RlVDvodouM8gPSlULFMxwpGY5GRYssQ90pPAH3+4ZSLQqBG0agV9+rgk0KSJSwilSkU0bOOxDvHCaOvWrSxZsgSAN998k44dO7JlyxY2btwIwGuvvcZ5553H3r17+euvv7j00ksZPXp0xom/QoUK7NnjhguuWLEidevW5e233wbcCXvFihU5xtC2bduMX/vTp0+nffv2uY6/U6dOR1UF7dq1i9atW/P555+zY8cOUlJSMuIpanbv3w1ApXJWUihMqvDjj/Dyy9CvHzRsCMcfD127wtNPu/r/Ll1g3Dj4+mvYs8dVEc2YAQ8/DN27w2mnWUIoSqykEEaNGzfmlVde4bbbbqN+/fqMGTOGc845h+7du2c0NPfv35+dO3fStWtXDh48iKry3HPPAdCjRw9uueUWxo4dyzvvvMP06dMZMGAAI0eOJCUlhR49enDGGWdkG8PYsWPp168fo0aNymhozq0HH3yQgQMHcvrppxMfH8+wYcPo1q0bw4cPp02bNpx00km0aNGCtLSid2PY7gMuKVQuVznCkcS+bdtgwQL49FOYPx9+9S4PqVoV2rZ1yaFdOzjrLHdJpokuohpdl/23bNlSM4+8tm7dOho3bhyhiJwtW7bQuXNnVq8O2dVTsef3/2jACwOYuH0in17+KR1bdPRtP8VRaip8+SW8/z588gmsWePmV6sGHTvCRRfBuedCgwbW2FuUicgyVW2Z03pWUjAxYe+hvQBUqVAlwpHEhn374OOPYfZs+PBD2LnTVfF06AC9esHf/w5nnAFxVgEdcywphEmdOnWippQwdepUxowZc9S8du3aMX78+AhFVHDpSaFaxWoRjiR6HTzoEsD06fDRR+55lSpw+eWujaBTJ3fHroltlhSKob59+9K3b9TfFnKU/Yf3A1ZSyKu0NPjiC5cI3n3XdfNw4olw661w5ZXQvr27A9gUH/bvNjHhUNohEChf2n7K5sa2bTBlCkyeDD/95Lpw6NYNevaECy6wRFCc2b/exITUQCoIxFkld5YCAddQ/OKLrtE4Lc01FD/5pLts1K4UMmBJwcSI9KRgjnXgALz2Gjz7rOtWonp1uPtuuOUWd9OYMcEsKZiYYEnhWMnJMGECjB/vplu0cG0HV1/tumM2JhQraxeyhx9+mPnz5xfKvvr06cM777wDwJdffslpp51G8+bNOXDgQKHsvzClBlKRgGUFcHcRDx0KtWvD8OGuS4mFCyExEa6/3hKCyZ6VFApRWloajz76aJ5fk10ndrk1ffp0hg4dGnNXHaVL1dRiP9JacjKMGuVKBgcPukbj++93/QsZk1sxlxQGD4YwD6dA8+aQU+erWY2n0KRJE/r168cnn3zCoEGDmDdvXkb32AsWLAg51kKdOnWOek2PHj2O2d/GjRvp378/ycnJxMfH8/bbb1OvXj3uuOMOPvvsM+rWrUv63eqTJ0/mrbfe4uOPP2b+/PlMnz49vAeoCEgNpEZFx31+2LvXJYNnnnHdUV9/PTz0kOuHyJi8suqjMFq/fj233norK1eupGLFikyYMAGA0qVLs3jx4qNO7gcPHqRPnz785z//YdWqVaSmpvLCCy9kLA/1mmA9e/Zk4MCBrFixgq+//pqTTjqJWbNmsX79elatWsVLL73E119/DcDNN99Mly5dGDVqVEwmBIA0TUO0eCWFtDR3WWmDBvDoo3Dppa4Litdft4Rg8i/mSgoRHE7hmPEUxo4dC7gxCDJbv359yLEWBg8enOVr0u3Zs4dffvmFK6+8EnAJBGDRokVcd911xMfHc/LJJ3PhhReG780VcamaihSjluaFC2HIEFixAs45x9141qZNpKMyscBKCmGUufoi/XmocQxy6ogwu7EPsnttca1CCWigWJQUfv8dbrgBLrzQ3X08Y4brktoSggkXSwphlHk8hezGMmjUqFHIsRZyo2LFitSsWTNjmM1Dhw6xf/9+OnTowIwZM0hLS+PXX39l4cKFBXxH0SNVU4nT2P04BwIwaZIboOatt1ybwbp1cO211jOpCa/Y/RZFQPp4Cs2aNWPnzp0MGDAgy3VLly7N1KlT6d69O02bNiUuLo7+/fvnel+vvfYaY8eOpVmzZrRt25bffvuNK6+8kvr169O0aVMGDBiQ6yQTCwIEYrb66McfXdfUt93meiZdudK1IXi1hsaElY2nECY2nkL2/P4fHT/4ePbJPvY9t8+3fRS2QADGjnWXlZYt6+5I7tXLSgYmf2w8BVOsBAgQF0MF3y1boG9f+Pxz6NwZXnrJ9V5qjN8sKYSJX+MpDBw4kK+++uqoeXfddVfM3oSWX2mkxUxSeOMN6N/fjX/88ssuOVjpwBQWSwpFXDQPfFOYYqGkcOAA3HWXKxW0b+86satTJ9JRmeImur9FxngCBIin4N2BRMr69dC6tUsI99/v7kOwhGAiwUoKJiZEc0nh7bddFVGZMm4YzEsuiXREpjiLzm+RMZkEJPpKCoGAu9/gmmugWTNYvtwSgok8KymYmKCiUVVS2LMHbrwR3nsP+vVz4x6UKhXpqIyxkkLYlC/vxgbetm0bV199NQDTpk1j0KBBud5G8PoTJ07k1VdfBeD8888n870Z4RCO7dapU4c//vgDgLFjx9K4cWN69uwZjvDyRImepLB5M7RtCx98AGPGuHGSLSGYosJKCmF28sknZwxsUxB5ubu5qJgwYQIfffQRdevWLfR9KxoVdzR//73rzfTwYZg3z42RbExREnNJYfC8wfzwW3gHVGh+YnNGX5K77lezurP5ww8/ZOTIkbz//vuoKv3792fr1q0AjB49OqN31XTDhw+nfPnyDB06FIC3336b22+/nT///JOXX36Zc889l4MHDzJgwAASExNJSEjg2Wef5YILLshy/oEDB+jbty9r166lcePGOY7ANm/ePB544AHS0tKoVq0aCxYsYMeOHVx33XUkJyfTqlWrjM75+vfvz6ZNm+jSpQv9+vVjyJAhuTpeYSMU+ZLCxx+7oTCPO85dXWSD35iiKOaSQlE0a9Ysnn32WebOnUuVKlW4/vrrGTJkCO3bt2fr1q1cfPHFrFu3LtttpKam8u233zJ37lweeeQR5s+fn3EPw6pVq/jxxx/p1KkTGzZsyHL+Cy+8QNmyZVm5ciUrV66kRYsWWe4vOTmZW265hUWLFlG3bl127twJwCOPPEL79u15+OGH+fDDD5k0aRLgqrvmzZvHwoULqVatWjgOW56oFO2Swquvwk03wWmnwdy5cPLJkY7ImNB8TQoicgkwBogHJqvqE5mWVwJeB07xYnlaVacWZJ+5/UVfWBYuXEhiYiKffPIJFStWBGD+/PmsXbs2Y53du3ezZ8+ebLfTrVs3AM466yy2bNkCwOLFi7njjjsA1+tq7dq12bBhQ5bzFy1axJ133glAs2bNaNasWZb7++abb+jQoUNGVdBxxx0HuDEbZs6cCcBll11GlSpV8nQ8/FKU2xRGjYJ77oGLLoKZM8H7GBhTJPmWFEQkHhgP/B1IAr4TkTmqujZotYHAWlW9XESqA+tFZLqqHvYrrsJWr149Nm3axIYNG2jZ0vVFFQgEWLJkCWXKlMn1dkp5LZHx8fGkpqYCWY+rEI7xFlQ1y3WL5JgNQpHrOlsVRoyAYcNcF9evvgolS0Y6KmOy5+e3qBWwUVU3eSf5GUDXTOsoUEHcWaY8sBNI9TGmQle7dm1mzpxJr169WLNmDQCdOnVi3LhxGev8kM9BpTt06JAxvOaGDRvYunUrDRs2zNX81atXs3Llyiy33aZNG7744gs2b94MkFF9FLyNjz76iF27duUrdj/ESdFJCqrwwAMuIfTuDdOnW0Iw0cHPb1EN4Oeg50nevGDjgMbANmAVcJeqBjJvSERuFZFEEUlMTk72K17fNGzYkOnTp9O9e3f+97//MXbsWBITE2nWrBlNmjRh4sSJ+dru7bffTlpaGk2bNuXaa69l2rRplCpVKsv5AwYMYO/evTRr1oynnnqKVq1aZbnt6tWrM2nSJLp168YZZ5yRMTzosGHDWLRoES1atOCTTz7hlFNOyVfs4aaSdcmmsKm6oTKfeMKNgTBlCsRH1311phjzbTwFEekOXKyqN3vPbwRaqeodQetcDbQD/gn8DfgUOENVd2e13aI6noLJnt//o/h742moDVn71NqcV/aRKtx5J4wb5zq3e+456+HUFA25HU/Bz5JCElAr6HlNXIkgWF9gpjobgc1AIx9jMjFK0YhXH6nCvfe6hPDPf1pCMNHJz6uPvgPqi0hd4BegB3B9pnW2AhcBX4rICUBDYJOPMZkQWrduzaFDh46a99prr9G0adMIRZQPReA+hREj3JVGAwbA009bQjDRybekoKqpIjII+Bh3SeoUVV0jIv295ROBEcA0EVkFCHCvqv6Rz/0VmTrlaLN06VJft18YQ75Guu+jZ5450qg8bpwlBBO9fL1PQVXnAnMzzZsYNL0N6FTQ/ZQuXZodO3ZQtWpVSwxFjKqyY8cOSvs9ynwESwovvghDh0L37q4fo7iicxGUMXkWE3c016xZk6SkJKLxyqTioHTp0tSsWdPfnUQoKbz3Htx+u+vP6PXXISEmvlGmOIuJj3CJEiUi0gmbKUIicPPakiXQowe0bAlvvWX3IZjYEBNJwZjCLils2ACXXw41asD770O5coW2a2N8ZUnBxAaBeCmcO8R++82NkBYX57q/Pv74QtmtMYXCkoKJDUKh9JJ64AB06QK//w6ffw6nnur7Lo0pVJYUTNRTVYiDePW3pKDqur9OTIRZs+Dss33dnTERYUnBRL1AwHWX5fcdzU88AW++CY89Bl0zd+1oTIywK6pN1EtJSwH8TQrvved6Pb3+erjvPt92Y0zEWVIwUS81zfW27ldSWLkSevZ01UWTJ9vdyia2WVIwUS8l1ZUU/Lj6aNcuuOIKqFQJZs+GPIyLZExUsjYFE/XSSwrxceFNCoEA9OoFSUmwaJGNq2yKB0sKJur51abw5JPwwQcwdiycc05YN21MkWXVRybq+VF9tHAhPPig68Zi0KCwbdaYIs+Sgol6aWlpAMSFqXvSX35xyaBBA3jpJWtYNsWLVR+ZqJdefRSOkkJKClx7Lezb50oL5csXeJPGRBVLCibqhbOhefhw+OoreOMNaNKkwJszJupY9ZGJeuG6T2HhQnj8cdeVxXXXhSMyY6KPJQUT9cJRUtixA268EerXhzFjwhWZMdHHqo9M1MsoKeSzoVkVbrkFtm+Hb76xsRFM8WZJwUS9gjY0T5rkej19+mlo0SKckRkTfaz6yES99EtSE+Ly/htn3ToYMgQ6dXJ/jSnuLCmYqJcayF/1UUqK6+iufHl45RU3kpoxxZ1VH5mol987mh97DJYvd1VHJ57oR2TGRB/7bWSiXlrAVR/l5eqj77+HkSPhhhtcL6jGGMeSgol6eU0Khw653k+PP951dmeMOcKqj0zUy7j6KJdJYdgwWLMG5s6FKlX8jMyY6GMlBRP18nLz2jffwKhRcPPN8I9/+B2ZMdHHkoKJeunVRzldknrgAPTuDTVrwjPPFEZkxkQfqz4yUS+3JYURI2DDBvj0U6hYsTAiMyb6WEnBRL3cJIUVK+Cpp6BvX+jYsbAiMyb6+JoUROQSEVkvIhtF5L4s1jlfRH4QkTUi8oWf8ZjYlFF9FB+64Jua6trDSZnMAAAYp0lEQVQQqlZ1XVkYY7LmW/WRiMQD44G/A0nAdyIyR1XXBq1TGZgAXKKqW0XkeL/iMbErPSlkdUfz2LGQmAgzZsBxxxVmZMZEHz9LCq2Ajaq6SVUPAzOArpnWuR6YqapbAVR1u4/xmBiV3s1FqOqjzZvhoYegc2e45prCjsyY6ONnUqgB/Bz0PMmbF6wBUEVEPheRZSLSK9SGRORWEUkUkcTk5GSfwjXRKr1NIfPVR6rQv7/r02jCBBtr2Zjc8PPqo1BfQQ2x/7OAi4AywBIR+UZVNxz1ItVJwCSAli1bZt6GKeayuiT19dfhk09g3DioVSsSkRkTffxMCklA8FexJrAtxDp/qOo+YJ+ILALOADZgTC6F6uZi50745z+hTRsYMCBSkRkTffysPvoOqC8idUWkJNADmJNpnfeAc0UkQUTKAq2BdT7GZGJQRvVR0NVH//437NoFEydal9jG5IVvJQVVTRWRQcDHQDwwRVXXiEh/b/lEVV0nIvOAlUAAmKyqq/2KycSmzJekfvcdvPgi3HUXNGsWyciMiT6+3tGsqnOBuZnmTcz0fBQwys84TGwLrj5KS3PVRSeeCI88EuHAjIlC1s2FiXrBSWHSJFi2DN54w7qyMCY/rLbVRL30+xT27S3FAw/AhRdCjx4RDsqYKGVJwUS9QCAAwIw3T2TfPhg/3u5JMCa/LCmYqJdeUvhiYRXuvhsaNYpwQMZEMUsKJuqlpLqSQrWqAR58MMLBGBPlLCmYqLdiZXUAbrvtT8qVi3AwxkQ5Swomqu3YAUu/dV1qndvuUISjMSb6ZZsURGRa0HRv36MxJo+GDYNDKW46ISHnMZqNMdnLqaRwRtD0XX4GYkxerV7turFo3MR1qVUirkSEIzIm+uWUFKxHUlMkqcLgwe4GteZnbgYgIcHuxTSmoHL6FtUUkbG4brDTpzOo6p2+RWZMNubMgQUL4Pnn4bu0Q7D/2K6zjTF5l9O36F9B04l+BmJMbh065LrFbtLEDaLzzfPZj9FsjMm9bL9FqvpKYQViTG6NHg2bNrkBdBISIE1dUiiRYG0KxhRUjpekikhvEfleRPZ5j8Sshs00xm+//QYjR0KXLvD3v7t56d1chBqj2RiTN9mWFLyT/2Dgn8D3uLaFFsAoEUFVX/U/RGOOeOABV330zDNH5qWXFKz6yJiCy6mkcDtwpaouVNW/VPVPVf0MuMpbZkyhWbYMpk6FIUPg1FOPzE/vOrtEvFUfGVNQOSWFiqq6JfNMb571Vm8KjSoMHQrVq7uhNoMF1FUf2SWpxhRcTt+iA/lcZkxYffABfP656xY78+A5Vn1kTPjk9C1qLCIrQ8wXoJ4P8RhzjJQU+Ne/oGFDuOWWY5enNzRb9ZExBZdTUjgDOAH4OdP82sA2XyIyJpPJk2H9enjvPSgR4ryfXn1kl6QaU3A5tSk8B+xW1Z+CH8B+b5kxvtq9G4YPhw4d4PLLQ69j1UfGhE9O36I6qnpM9ZGqJopIHV8iMibIU0/B9u2uTSGrITbTrz6ypGBMweVUUiidzbIy4QzEmMySktz9CNdfD2efnfV6GdVH1qZgTIHllBS+E5FjmvZE5CZgmT8hGeM8+KC7FPX//i/79dKTQlycjRllTEHlVN4eDMwSkZ4cSQItgZLAlX4GZoq3H36AV191Vx3VqZP9ugENQKBQwjIm5uXUId7vQFsRuQA43Zv9oXdXszG+SL9R7bjj4P77c14/vaRgjCm4XLXMqepCYKHPsRgDwLx5bqyEMWOgcuWc17ekYEz42OUapkhJTXVVRqee6sZKyI30S1KNMQVnScEUKdOmwZo18O67ULJk7l5jJQVjwseSgiky9u6Fhx6Cdu3gyjxcxmBJwZjwsaRgioynn3aD6MyalfWNaqGoKkIeXmCMyZKvF3aLyCUisl5ENorIfdmsd7aIpInI1X7GY4qubdtg1Ci45ho455y8vTZAANSfuIwpbnxLCiISD4wH/gE0Aa4TkSZZrPck8LFfsZiib9gw1xvq44/n/bVpmmZJwZgw8bOk0ArYqKqbVPUwMAPoGmK9O4B3ge0+xmKKsNWrYcoUGDQI6uWjQ3arPjImfPxMCjU4usvtJG9eBhGpgbszemJ2GxKRW0UkUUQSk5OTwx6oiax77nED5zz4YP5eH1CrPjImXPxMCqF+umX+6o4G7lXN/kJzVZ2kqi1VtWX16tXDFqCJvAUL4KOP3BCbxx2Xv20ECFhJwZgw8fPqoySgVtDzmhw7ME9LYIa4S02qAZeKSKqqzvYxLlNEBALuRrXatV3VUb63Y5ekGhM2fiaF74D6IlIX+AXoAVwfvIKq1k2fFpFpwAeWEIqPN96A5cth+nQonV0n7TlQFFErKRgTDr4lBVVNFZFBuKuK4oEpqrpGRPp7y7NtRzCx7eBBV2V01lnQo0fBthXQQOjKSmNMnvl685qqzgXmZpoXMhmoah8/YzFFy9ixsHWr69aioMMgBNTaFIwJFxuVxBS6HTvgscfgssvgggsKvj3FLkk1JlwsKZhCN2IE7Nnjxl8OB7uj2ZjwsaRgCtX//gcTJsBNN0GTY+5vzx+7ec2Y8LGkYArVAw9AiRLwyCPh26bdp2BM+FhSMIVm6VJ46y031OZJJ4Vvu3ZJqjHhY0nBFIr0cZdPOMHdsBbWbVtDszFhY+MpmELx3nuweDFMnAjly4d32wENIHkZgMEYkyUrKRjfpaTAvfdC48augTncrKRgTPhYScH4bvJk2LAB5syBBB8+cZYUjAkfKykYX+3ZA8OHw3nnQefO/uwjQMAamo0JEyspGF89/jhs3w4ffJC3cZfzQlFrUzAmTKykYHyzZQs8+yzceCOcfbZ/+7HqI2PCx5KC8c2997rO7h57zN/9WFIwJnwsKRhfLF7sblS7916oWdPffdnNa8aEjyUFE3aBAAwZAjVquBvW/KYocWIfZWPCwRqaTdi9/jokJsJrr0G5cv7vT8Wqj4wJF/t5ZcJq3z64/37XsHz99TmvHw7WpmBM+FhJwYTVU0/Btm3w9tsFH1EttxQlzn7fGBMW9k0yYfPzzzBqFFx7LbRtW7j7tvsUjAkPSwombO6/3zUyP/lk4e7X2hSMCR9LCiYsvvoKpk+Hu++G2rULd99WfWRM+Fibgimw1FQYOBBq1XIjqxU2FUsKxoSLJQVTYC++CCtWuMblwrgENTMrKRgTPpYUTIEkJ8ODD8JFF8FVV0UmBispGBM+lhRMgdx3H+zdC88/718vqDkSLCkYEyb2TTL5tnQpTJniurRo3DhycahYNxfGhIt9k0y+pKW5xuWTT4aHHopsLCpKPPGRDcKYGGHVRyZfJk+GZcvgjTegQoUIByMQL5YUjAkHSwomz37/3bUlnH8+9OgR6WispGBMOFlSMHk2eDDs3w8TJ0awcTmYNTQbEza+fpNE5BIRWS8iG0XkvhDLe4rISu/xtYic4Wc8puA++ghmzIB//xsaNox0NB6rPjImbHxLCiISD4wH/gE0Aa4TkSaZVtsMnKeqzYARwCS/4jEFt28fDBjgrjS6995IRxMkzpKCMeHiZ/VRK2Cjqm4CEJEZQFdgbfoKqvp10PrfAD4P3GgKYtgw+Okn+PJLKFUq0tEEsZKCMWHjZ/VRDeDnoOdJ3rys3AR8FGqBiNwqIokikpicnBzGEE1uff89PPcc3HYbtG8f6WgyicPuUzAmTPz8JoVqgtSQK4pcgEsKISslVHWSqrZU1ZbVq1cPY4gmN1JS4Oab4fjj4YknIh3N0dLS0lxJIc5KCsaEg5/VR0lAraDnNYFtmVcSkWbAZOAfqrrDx3hMPj32GCxfDrNmQeXKkY7maIdTDwNWfWRMuPhZUvgOqC8idUWkJNADmBO8goicAswEblTVDT7GYvLp++9h5Ei44Qa44opIR3OswykuKSTE2dXVxoSDb98kVU0VkUHAx0A8MEVV14hIf2/5ROBhoCowwRtOMVVVW/oVk8mbQ4egd2+oXh3Gjo10NKFZScGY8PL155WqzgXmZpo3MWj6ZuBmP2Mw+ffII7B6NXzwAVSpEuloQstICtamYExY2CUbJqSlS91Yy337wmWXRTqarKWkpgCWFIwJF0sK5hi7d0PPnlCjhrsMtShLLylYm4Ix4WHfJHOMgQNh82b44guoVCnS0WTPqo+MCS8rKZijvPYavP66u3u5yN2kFoJdfWRMeFlSMBn++1+4/Xbo0MF1eBcN0tsULCkYEx6WFAwABw/CdddBiRKupBAfJbUxKWnW0GxMONnPK4Oqa0dYtgzeew9q1cr5NUWFNTQbE15WUjC89BJMmQIPPghdukQ6mrzJaFOIt6RgTDhYUijmli6FO+6Aiy+G4cMjHU3epVcflYgvEeFIjIkNlhSKsd9/h6uucvcjvPFG9LQjBLOGZmPCy75JxdT+/XD55bBrFyxeDMcdF+mI8ie9TaFEgpUUjAkHSwrFUFqa6/U0MdF1h33mmZGOKP8OHj4IQOkSpSMciTGxwZJCMXTPPS4ZjB4NXbtGOpqC2X94PwBlSpaJcCTGxAZrUyhmnn3WPe64A+66K9LRFNyBwwcAKJVQlAaNNiZ6WVIoRiZNgrvvhquvLvod3eXWoZRDAJQtVTbCkRgTGywpFBPTp0P//nDppW46Gq80CuVgimtTsOojY8LDkkIx8PbbbgS1886Dd96BkiUjHVH4pJcUSpe0hmZjwsGSQoybMgV69IA2bWDOHCgTYz+o00sKVn1kTHhYUohho0fDTTdBx47w8cdQoUKkIwo/SwrGhJclhRgUCMADD8CQIdCtmyshlI3Rc+bhNHfzmrUpGBMedp9CjNm3D3r1gpkz4ZZbYMIESIjh//LBVFdSKFe6XIQjMSY2xPDpovjZvNldbrp8ubsXYfBgEIl0VP5K7+bCqo+MCQ9LCjFi1izo29eNjTBnDnTuHOmICkd69VHZ0pYUjAkHa1OIcvv2waBBru2gfn1XSiguCQGOlBTKly4f4UiMiQ2WFKLYggXQtCmMH++qihYvhnr1Ih1V4TocsJKCMeFkSSEK/fKLuxmtY0d3Z/KiRa7bilLFsPufQ6mHIM3GaDYmXCwpRJE9e+DRR6FBA5gxA+67D1asgHPPjXRkkbM3ZS9yOMZb040pRNbQHAV27oTnn4cxY9ygOFddBU89VfyqikLZl7qPeKyUYEy4WFIowpYvdz2bvv467N0LXbq4m9Jat450ZEXH/rT9JIh9jI0JF/s2FTGbN7vLS994A5Ytg9Kl4ZprYOhQ16hsjnZQD1KSGOrhz5gIs6QQYXv3wldfweefu/6Jli9385s3h7Fj3bCZVapENMQi7bAcphx2N7Mx4eJrUhCRS4AxQDwwWVWfyLRcvOWXAvuBPqr6vZ8xRYoqJCfDhg3www+ugXj5cvc3NdV1RXHOOTBqlLvnwNoLcudw/GGqa/VIh2FMzPAtKYhIPDAe+DuQBHwnInNUdW3Qav8A6nuP1sAL3t8iSxUOH4b9++HAgaP/7toFf/wBO3a4v8nJkJTkqoS2bHHrpKtaFc48E/71L7jgAmjbFsrZD948+W3nb6SWT6VenGVQY8LFz5JCK2Cjqm4CEJEZQFcgOCl0BV5VVQW+EZHKInKSqv4a7mD6jHieV7aPAbzLF9VbcNTVjDld2pjNctFj55UMwN9SkFNTEVKJk1REUhE5zF+k8jnw+SEYNQ+Yl4vdZ0MJsf9jYsz/9nPcT7iuCs1iO6H2qyUVEuCihheFaefGGD+TQg3g56DnSRxbCgi1Tg3gqKQgIrcCtwKccsop+QrmhCpVKLOxNunZwHUUp95D3Lko48SuR52bjqwLIoGMR1z6dFz68zTi41OIi0slLi7Ve50AJbyHt71QZz5vlqiEnp+Ls26W6wTNztj+MbvJfvvZxVyQbR+zzPsXSG568jsMDSs05P7u9+e8rjEmV/xMCqG+1Zl/7uVmHVR1EjAJoGXLlrn4SXysJwfdwJPckJ+XGmNMseHnHc1JQK2g5zWBbflYxxhjTCHxMyl8B9QXkboiUhLoAczJtM4coJc45wB/+dGeYIwxJnd8qz5S1VQRGQR8jLskdYqqrhGR/t7yicBc3OWoG3GXpPb1Kx5jjDE58/U+BVWdizvxB8+bGDStwEA/YzDGGJN71kuqMcaYDJYUjDHGZLCkYIwxJoMlBWOMMRnEtfVGDxFJBn7K58urAX+EMZxoZ8fjaHY8jrBjcbRYOB61VXPuPTLqkkJBiEiiqraMdBxFhR2Po9nxOMKOxdGK0/Gw6iNjjDEZLCkYY4zJUNySwqRIB1DE2PE4mh2PI+xYHK3YHI9i1aZgjDEme8WtpGCMMSYblhSMMcZkiOmkICLxIrJcRD7wnh8nIp+KyH+9v1UiHWNh8YY6fUdEfhSRdSLSppgfjyEiskZEVovImyJSujgdDxGZIiLbRWR10Lws37+I3C8iG0VkvYhcHJmo/ZPF8RjlfV9WisgsEakctCxmj0dMJwXgLmBd0PP7gAWqWh9Y4D0vLsYA81S1EXAG7rgUy+MhIjWAO4GWqno6rmv3HhSv4zENuCTTvJDvX0Sa4I7Pad5rJohIfOGFWiimcezx+BQ4XVWbARuA+yH2j0fMJgURqQlcBkwOmt0VeMWbfgW4orDjigQRqQh0AF4GUNXDqvonxfR4eBKAMiKSAJTFjfhXbI6Hqi4CdmaandX77wrMUNVDqroZN/5Jq0IJtJCEOh6q+omqpnpPv8GNDAkxfjxiNikAo4F7gEDQvBPSR3bz/h4ficAioB6QDEz1qtMmi0g5iunxUNVfgKeBrcCvuBH/PqGYHo8gWb3/GsDPQeslefOKk37AR950TB+PmEwKItIZ2K6qyyIdSxGRALQAXlDVM4F9xHbVSLa8uvKuQF3gZKCciNwQ2aiKNAkxr9hcyy4i/wZSgenps0KsFjPHIyaTAtAO6CIiW4AZwIUi8jrwu4icBOD93R65EAtVEpCkqku95+/gkkRxPR4dgc2qmqyqKcBMoC3F93iky+r9JwG1gtariatui3ki0hvoDPTUIzd1xfTxiMmkoKr3q2pNVa2DaxD6TFVvAOYAvb3VegPvRSjEQqWqvwE/i0hDb9ZFwFqK6fHAVRudIyJlRURwx2Mdxfd4pMvq/c8BeohIKRGpC9QHvo1AfIVKRC4B7gW6qOr+oEUxfTx8HaO5CHoCeEtEbsKdGLpHOJ7CdAcwXURKApuAvrgfBcXueKjqUhF5B/geVy2wHNeNQXmKyfEQkTeB84FqIpIEDCOL74eqrhGRt3A/JFKBgaqaFpHAfZLF8bgfKAV86n478I2q9o/142HdXBhjjMkQk9VHxhhj8seSgjHGmAyWFIwxxmSwpGCMMSaDJQVjjDEZLCmYmCMiz4nI4KDnH4vI5KDnz4jIP3PYxte52M8WEakWYv75ItI2m9ddISIPe9PDRWSoN/25iORqcHgRqS4i83KzrjF5YUnBxKKvcXcoIyJxQDVcj5bp2gJfZbcBVc3ypJ4L56fvPwv3ABMKsH1UNRn4VUTaFWQ7xmRmScHEoq84clI+DVgN7BGRKiJSCmiMu2ENEfmXiHzn9Zn/SPoGRGSv9zdORCZ4Yy98ICJzReTqoH3dISLfi8gqEWkkInWA/sAQEflBRM4NDkxEGgCHVPWPLGK/QUS+9sZ5aOW95jxvWz94HRpW8NadDfTM91EyJoTidkezKQZUdZuIpIrIKbjksATXi2Ub4C9gpaoeFpFOuC4KWuE6OZsjIh28bpTTdQPqAE1xvYauA6YELf9DVVuIyO3AUFW9WUQmAntV9ekQ4bXD3UmdlXKq2lZEOnj7OR0Yirtr9isRKQ8c9NZNBEbm9rgYkxtWUjCxKr20kJ4UlgQ9T28v6OQ9luNO1I1wSSJYe+BtVQ14fUgtzLR8pvd3GS555OQkXDfmWXkTMvr3r+iN9vUV8KyI3AlUDurjfzuul1djwsaSgolV6e0KTXHVR9/gSgrB7QkCPK6qzb3Hqar6cqbthOomOdgh728auSt5HwBKZ7M8c78zqqpPADcDZYBvRKSRt6y0tz1jwsaSgolVX+G6PN6pqmmquhOojEsMS7x1Pgb6eVUyiEgNEck8sM5i4CqvbeEEXCNyTvYAFbJYtg44NZvXXuvF0h43+M9fIvI3VV2lqk/iqozSk0IDXMIzJmwsKZhYtQp31dE3meb9ld7I64229gawRERW4caZyHwyfxfXf/5q4EVgKa5dIjvvA1eGamgGFgFnel12h7LLuxx2InCTN2+w1/C8AlcySB8B7ALgwxxiMSZPrJdUY3IgIuVVda+IVMX1m9/Oa1/I7/bGAO+r6vwCxrUI6KqquwqyHWOC2dVHxuTsA6/BtyQwoiAJwfMY0LogGxCR6sCzlhBMuFlJwRhjTAZrUzDGGJPBkoIxxpgMlhSMMcZksKRgjDEmgyUFY4wxGf4fRK9Tk7fEh2sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(year_10_pdf['lbs'], year_10_pdf['posterior_cdf'], label = 'posterior_cdf',\n",
    "        linestyle = '-', color = 'k')\n",
    "plt.plot(year_10_pdf['lbs'], year_10_pdf['prior_cdf'], label = 'prior_cdf',\n",
    "        linestyle = '-', color = 'b')\n",
    "plt.plot(year_10_pdf['lbs'], year_10_pdf['likelihood_cdf'], label = 'likelihood_cdf',\n",
    "        linestyle = '-', color = 'g')\n",
    "plt.xlabel('Weight (lbs)')\n",
    "plt.ylabel('CDF')\n",
    "plt.title('Bayesian Baby Weight @ 10 years')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Posterior mode at 10 years: 75.24 lbs\n"
     ]
    }
   ],
   "source": [
    "# Location of max posterior\n",
    "max_posterior_lbs_10_yrs = (year_10_pdf\n",
    "                            .loc[year_10_pdf['posterior_pdf'] == np.max(year_10_pdf['posterior_pdf']), 'lbs']\n",
    "                            .iloc[0])\n",
    "print(\"Posterior mode at 10 years:\", str(np.round(max_posterior_lbs_10_yrs, decimals = 2)), \"lbs\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}