bayesian_logistic_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Methods - Bayesian Logistic Regression\n",
    "\n",
    "Using six features (seven including constant) we are to classify whether or not a woman is working or not. This is done by logistic regression, in which the weighting parameters follow some prior distribution, in the Bayesian approach. Recall, in logistic regression we are trying to find parameters $\\mathbf{w}$ for the features which are then parsed into the Sigmoid function.\n",
    "\n",
    "$$\n",
    "\\sigma(z) \\triangleq \\frac{1}{1-e^{-z}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(Y_i = 1 | \\mathbf{x}, \\mathbf{w}) = \\sigma(\\mathbf{w}^{\\text{T}}\\mathbf{x})\n",
    "$$\n",
    "\n",
    "The prior is set to:\n",
    "$$\n",
    "\\mathbf{w} \\sim N(\\mathbf{0}, \\tau^2\\mathbf{I}_7)\n",
    "$$\n",
    "\n",
    "With the likelihood function being:\n",
    "\n",
    "$$\n",
    "p(\\mathcal{D}|\\mathbf{w}) = \\prod_{(\\mathbf{x}, y) \\in \\mathcal{D}} [\\sigma(\\mathbf{w}^{\\text{T}}\\mathbf{x})]^y [1-\\sigma(\\mathbf{w}^{\\text{T}}\\mathbf{x})]^{1-y}\n",
    "$$\n",
    "\n",
    "And the posterior distribution having the form derived from Bayes' rule:\n",
    "$$\n",
    "p(\\mathbf{w}|\\mathcal{D}) = \\frac{p(\\mathcal{D}|\\mathbf{w})p(\\mathbf{w})} {p(\\mathbf{\\mathcal{D}})}\n",
    "$$\n",
    "\n",
    "Where the evidence is calculated as:\n",
    "$$\n",
    "p(\\mathbf{\\mathcal{D}}) = \\underset{}{\\int_{\\mathbb{R}} \\cdots \\int_{\\mathbb{R}}} \\prod_{(\\mathbf{x}, y) \\in \\mathcal{D}} [\\sigma(\\mathbf{w}^{\\text{T}}\\mathbf{x})]^y [1-\\sigma(\\mathbf{w}^{\\text{T}}\\mathbf{x})]^{1-y}p(\\mathbf{w})~dw_1\\cdots dw_7\n",
    "$$\n",
    "\n",
    "\n",
    "The prior distribution is non-conjugate. In addition, there exists no closed-form expression of the posterior distribution, which is easy to see, since Bayes' factor is computationally intractable to calculated due to the curse of high dimensionality.\n",
    "##### Laplace Approximation\n",
    "However, for large samples sizes, we can normally approximate the posterior distribution. Then it is simple to find samples from the posterior. This is also called **Laplace Approximation**.\n",
    "\n",
    "The approximation has the form:\n",
    "$$\n",
    "\\mathbf{w}|\\mathcal{D} \\sim N(\\mathbf{w}_{\\text{MAP}}, -\\nabla^2\\ln(p(\\mathbf{w}_{\\text{MAP}}|\\mathcal{D})))\n",
    "$$\n",
    "\n",
    "\n",
    "We can arrive at the posterior normal approximation by Taylor expanding the posterior distribution around $\\mathbf{w}_{\\text{MAP}}$. There are two paths forward. (i) is to calculate the Taylor approximation of the log posterior by hand (requiring the usage of the likelihood function and prior). (ii) is to find the $\\mathbf{w}_{\\text{MAP}}$ numerically with some optimzation algorithm, by maximizing the likelihood times the prior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 488,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np \n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(['science', 'notebook', 'grid'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 441,
   "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>Work</th>\n",
       "      <th>Constant</th>\n",
       "      <th>HusbandInc</th>\n",
       "      <th>EducYears</th>\n",
       "      <th>ExpYears</th>\n",
       "      <th>Age</th>\n",
       "      <th>NSmallChild</th>\n",
       "      <th>NBigChild</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>22.394940</td>\n",
       "      <td>12.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>7.232000</td>\n",
       "      <td>8.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>34.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>18.271990</td>\n",
       "      <td>12.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>41.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>28.069000</td>\n",
       "      <td>14.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>7.799889</td>\n",
       "      <td>12.0</td>\n",
       "      <td>10.0</td>\n",
       "      <td>31.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>163</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>25.075040</td>\n",
       "      <td>12.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>43.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>164</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>27.799960</td>\n",
       "      <td>17.0</td>\n",
       "      <td>8.0</td>\n",
       "      <td>54.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>165</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>20.200000</td>\n",
       "      <td>13.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>53.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>166</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>42.250000</td>\n",
       "      <td>12.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>30.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>167</th>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>13.390000</td>\n",
       "      <td>12.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>44.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>168 rows × 8 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     Work          Constant      HusbandInc    EducYears      ExpYears  \\\n",
       "0     1.0               1.0       22.394940         12.0           7.0   \n",
       "1     0.0               1.0        7.232000          8.0          10.0   \n",
       "2     1.0               1.0       18.271990         12.0           4.0   \n",
       "3     0.0               1.0       28.069000         14.0           2.0   \n",
       "4     1.0               1.0        7.799889         12.0          10.0   \n",
       "..    ...               ...             ...          ...           ...   \n",
       "163   1.0               1.0       25.075040         12.0           9.0   \n",
       "164   1.0               1.0       27.799960         17.0           8.0   \n",
       "165   0.0               1.0       20.200000         13.0           7.0   \n",
       "166   0.0               1.0       42.250000         12.0           3.0   \n",
       "167   0.0               1.0       13.390000         12.0           3.0   \n",
       "\n",
       "         Age            NSmallChild    NBigChild  \n",
       "0       43.0                    0.0          3.0  \n",
       "1       34.0                    0.0          7.0  \n",
       "2       41.0                    1.0          5.0  \n",
       "3       43.0                    0.0          2.0  \n",
       "4       31.0                    0.0          1.0  \n",
       "..       ...                    ...          ...  \n",
       "163     43.0                    0.0          3.0  \n",
       "164     54.0                    0.0          0.0  \n",
       "165     53.0                    0.0          0.0  \n",
       "166     30.0                    1.0          1.0  \n",
       "167     44.0                    0.0          1.0  \n",
       "\n",
       "[168 rows x 8 columns]"
      ]
     },
     "execution_count": 441,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Data import. \n",
    "df = pd.read_csv('data/WomenAtWork.dat', delimiter='\\t', engine='python')\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 463,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x: np.array) -> np.array:\n",
    "    return np.exp(x) / (1 + np.exp(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 443,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(X: np.array, t: np.array, w: np.array) -> float:\n",
    "    epsilon = 1e-3\n",
    "    prob = sigmoid(X @ w)\n",
    "    likelihood_vector = t @ np.log(prob + epsilon) + (1 - t) @ np.log(1 - prob + epsilon)\n",
    "    return np.sum(likelihood_vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 444,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import multivariate_normal\n",
    "\n",
    "\n",
    "def log_prior(w: np.array, init_sigma: float) -> float:\n",
    "    epsilon = 1e-3\n",
    "    mean = np.zeros(len(w))\n",
    "    covariance_matrix = init_sigma * np.identity(len(w))\n",
    "    return np.log(multivariate_normal.pdf(w, mean, covariance_matrix) + epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 445,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(w: np.array, X: np.array, t: np.array, init_sigma: float) -> float:\n",
    "    return -(log_likelihood(X, t, w) + log_prior(w, init_sigma=init_sigma))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 459,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "# Stable initial solution.\n",
    "x0 = np.array(\n",
    "[-1.7012604510755014,\n",
    " -1.1196676407493864,\n",
    " 2.639602500390759,\n",
    " 1.088785147860549,\n",
    " 0.3254163608102913,\n",
    " 3.4308174226161103,\n",
    " 1.8217426564308559]\n",
    " )\n",
    "\n",
    "# Parameters.\n",
    "X = df[df.columns[1:]]\n",
    "X = np.array(X)\n",
    "\n",
    "t = df[df.columns[0]]\n",
    "t = np.array(t)\n",
    "init_sigma = 5\n",
    "\n",
    "\n",
    "# Optimization results.\n",
    "solution  = minimize(objective, x0=x0, args=(X, t, init_sigma))\n",
    "\n",
    "# MAP weights\n",
    "w_map = solution.x\n",
    "\n",
    "# Inverse Hessian\n",
    "inv_hessian = solution.hess_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 475,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_posterior(x: np.array, mean: np.array, covariance_matrix: np.array) -> np.array:\n",
    "    w = multivariate_normal.rvs(mean, covariance_matrix)\n",
    "    return sigmoid(x.T @ w) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 489,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAINCAYAAACDPSPaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAABYlAAAWJQFJUiTwAABEXElEQVR4nO3deZwkdX3/8fcHuQQGOUTYBZEoGYN44ajxwozGI2hGjfcuInjtYtRfEkcf0ShmMMZ7jEmMZtcLRXY9QR1FxWu8zxFPlBUVBXcxHhzDrezn90dV79QMfXy7p7/9rf726/l49GOru79V9eneT1V95tv1rTJ3FwAAAOprt9QBAAAAoD0KNgAAgJqjYAMAAKg5CjYAAICao2ADAACoOQo2AACAmqNgAwAAqDkKNgAAgJqjYBsRZubl4+LUsaxWp89iZkdV2swPNrp4cvo/TMXM7mdmZ5nZz8zsWr7TQt23mVjxhfz/B+xvyKEO2O76Y/fUAaBQJu7tmrzlkq6RdJWkKyVdKOl8Sd+Q9Fl3/9OgYgxhZjPl5BXu/saEodQe39Vgmdmpkt4syfq0vIvVfJuVpJtUbK+XSPq2pPeq2F65tcwIGuVtvd/b3SijYKs/k7Rf+Vgr6RhJjynfu8zM3inp1e5+VZrwbuZfy39/KemNCeMYBnxXA2Jmh0h6g5YOGudJ+pakK8rnV/Z5lbeQdFD5uJukZ0j6kpk9xd1/1ed1of5GcltPsN1ljYKtnv5X0s8qz/eSdICkwyTdU9IdVWwAh0l6saQTzewUd//8gOMEhsVDJd2ynN7q7uv7vPyV2+wekg6XdIKk25evHS/pc2Z2b3f/Q5/XD9RR7O1upFCw1dP73H2+1ZtmdqSkv5f0jyqKuSMlnWtmD3H3rzSbx92z6Y7O6bN0Y1Q/d58cWZn+UoTlN91mzWw3SS+R9PLypTtIepmKbReJ9WObYrtsK/Z2N1IYdDCE3P1X7v4iSfeV9Ovy5b0lfcTMDk4XGVBbe1emrxvUSt19p7v/m6T3VV5+SlnIAblLst3lip3GEHP38yU9VtKN5UsHS3pBs7aho3LM7Dgze5OZnW9mV5jZn8pRPZeY2efN7F/LET+3qMyzawRXZVG3q6yz+phvNl/jdTPbw8yeZmafNLNfmdmN5ftHdftZVnyuvczsVDP7opldZmbXm9nFZnaGmd07YP7JynrPCGg/U2l/SrPPXGne8bvq9nOXn3eDmX3czC4tP+/lZvZDM/svMzsuYBnN/n/MzNab2afMbLuZ3WBmO8zsHDN7WKdldms1n6P6f6alc4gk6Z1Nvuuj+h37Cu+sTB+spZ9JG7H2tC1U5l/1/3eLZfa8zZTLMDO7r5mdbmaftmI/cl25rO1lHv2jme2XKL5Vj1ZstYxetvXy+/pF5fVjAmM4oTLPF3v9LJXl1Xa7M7NXVeb/Qoe2u5vZuZX255nZHt2uszbcnUcNHpIuVjEi1CVNdjnv2yrzLkraq0mbxvsXt1iGSXqdpJ2Vtu0ed67Me1TgPC5pvsV88yq6z7/VYr6juvgsK5d7uIqReq1iuknSv3f4jicr7c8I+D+ZqbQ/ZbXfVcjnrrS7t6RfdFj2TklvkbR7m+Ws/B4PUnHScLvlzvZxm1jV51jxf9bpcVQP8V1cmb/tNitpfMX67tvhuw7aFvr1/91k/aveZsrlfinw+/+NpAd2kYv9iq/jNtWpTav31ft+8aWV118XmIsfrMzz1My3uwMk/b6yjIe2aXtGpd03Je23mu8m9YNz2PLwPypGoUnFaNK/lNTtX1nP0/Leua+ouHTI71Sc63iIpLuo+Bl2rxXz/kHSC8vp15X/Xi7plU3Wc0mL9e8t6cOSjpP0f5I+pmKnsZ+k+6vY4Hqxh4qd2YSKz/JhFQfagyQ9UsUAjt0k/YuZ/dHdZ3pcT6h+fFctmdlfSvqspH3LlxYlfUTF5WD2lfRgFTtkk3SqpENV9NJ2sruK7/FBki6TNKdixNuYihPr71q2e76ZLbj7lm5jj/A5fqal7/phKk6AloqfJ7+9om3sQQAre5CubdM2eFuI9P/dz23mkPLfX6vYn/xUxcjAPVScz/c3km5TPj5hZvdy9wsGGF9MvW7r71DxB98tJJ1kZi/2NpdvMrNbS3pU+fRKFd9NT4Zhu3P3K8zsVVr6Tv9N0qebfJZXSzq5fHqhpEe4+9Xdrq9WUleMPIqHVtfDZiqGSTfmf2mTNp3+SvxV+f6Nkh7SZl37SjpF0uEt3u/4F2ul7VG6+V9cZ0jat8N8nT5LdbmNHsM5Sbda0W43FV32jbZ/knTPFsucrMYY8NlmKu1PWe13Ffi5b6liZ9lo90VJhzZpd2L5/9xo95wu/n/eIumWTfLvjZU2P1nlttDXzxH6/9FDnBdXltl2m5X0/BV5dsBqt4V+fk8xtpmy/emS7tXm/T0lvbayvM8G5GI/4+u4DXZq049lNGn/kco8j+nQ9p8qbd+yinweiu2uXO7eWjpmuaSpFe8/r/LeryXdrl/rTvlIHgCP8j9iFQVbOf8XK/O/tcn7LXcYKs6pabw/t8rPEbxj0s0PUl+UZKtdR5Pl/kRNfiautD+j0vbsFm0mK23OCIix446qh514p8+9sdLmUkn7t1nW/6u03S5pz4Dv8WNtlrfHih3osavIob5+jtD/jx7iDNpmJd1WRa9ko+2X+rEt9PN7irHNdPldVguUowO+n77EF7INBmx3q15Gk/Z/W5mn7T5Z0g8qbVsWpwHrHIrtrrLsp1WWfX5je5H0RBU/ibuKHs0793O9KR8MOshHtWv5oC7nvUVlep8+xNKrl3u5xUVY7g1t3j9NxcYtSY8yswMjxDAIJ1emX+HtL6b8JhUFliStUfHTRScvbfWGu/9R0kcrL3V9kntF7M8RXXmy8+3M7NkqzkU7tPL2TMAiQraFmN/ToLeZrZXp4wPa575Nf0JF0SRJJ5jZmmaNysEVdy6ffs/dV/7k2I1h2+7eJelH5fTdJT3ezB4k6UwVPa3Xqeh5+2GC2KKgYMtH9bf5/buc97cqegAk6UFm9g8JRtJcJ+nzEZa7U8sLiZtx90tUHFSlongNGmFWJ2a2l4pzehransfi7jslfajy0gM6rOIyd/9uhzbbKtO36dC2qQF8jlg+Xx39JumPKnrg3qylYm2npH9y9892WFbHbSHy9xRlmzGzw8vRjM8oR4a+oPGQ9FeVpnfssKjst2l3v0nFuWxSEf/JLZo+ozL99l7XN4zbXRnDv1Reeo2K8xn3VNHD9iR3//Kg44qJgi0fY5Xprm5TVf4l/5ryaeN8pB1m9h4z22hmd7f41436abmT6reLPOxE0+9XpjsdMOrodip2VJJ0qbv/LmCe71Sm/7xD218GLG+xMt31ZRpKsT9HCjtVnBT9QA+7j2TIthDze+rrNmNmTzKzb6voMTpXxaj2/1Bx0njjcWpllgM6rHdUtul3qMgdqfj5bxkz20fSk8un10t6zyrWNZTbnbt/VFKjKPszLXVWPMvd51LEFBOjRPNR/Rn0993O7O5vLK+F9BIVJ3QerOLk0hPLJleY2Zyk/3X3r6422CZi3VMu9Luo7qCG7ecTaXnMITvble06/Yx+fcDyqj/h9Vrgx/4csay8NdWfVPzh9CtJ3/HubkUVsi3E/J76ss2YmUl6q5b3AoXYu8P7I7FNu/svzew8FSNpx83seHev3i3g8VoqUM5x98tXsbph3e4k6d1a3sP3Mnd/Z6pgYqJgy0DZ+3XXykshvSE34+6vMLO3SVov6SEqLiHQ2CEcIOkkFcPMz1DxF0zLoeY92Nm5CVBbbW8n16VctoVnaKlY+6OKHqM5ST9WcbmS6xo9ieW5R59LEWTNbVZRsEnFd1kt2KqF8NsGFlGNmNldVIwyrkp5HnZU/CSah7tr+Xlrba/+3I67X+bub3D3R6j4q+seKq7P9s1Ks1NUDNcfBqG36rp1ZbrZX6rV3qOQewcOeqdRjbmXz1yXm5Hn8jlii/k99Wub+X+V6Se4+6nu/nF3/7m7X73iZ99bBa6zn/ENgzktnV/8eDMbkyQzO1pLgzN+rtWf/zt0252Z3U7SJ7X0E3ojn/7BzNYOOp5BoGDLw99Xpq/S8uKqZ17cB/F8d59197/U0sUQJenUAZzX1g9Hm9m+nZst66G8sMn71XNmQs7POrJzk776pZZuUXZbC7unbHUk57aWrQYrl88RW8zvadXbTHl+1V3Kpz939490WNadO7zf1/iGRfkrRuPnvX21dM7a07X0h+M7+jC6fqi2u/JiwedJahRmr1BxjUipuJ7czCDjGZRhOOCijXJY90mVl97k7je2ar9Kb5B0TTl9kJauYl7V+Jn0Fk3eS2E3LV0FvCkzO0LSPcunN6l5wfubyvR4h+XtqeKOAJ307bsqL3GwUHnp8e3al8X24yovfWW1MfRDLp8jtsjfUz+2meo5USG9L4/r3GSXfm3Tg7Sabb1x60FJeroV93FujBq9ScU151ZlmLa78lzrc7W0H367u58m6d+1dAeRp5vZMA40aYuCbYiZ2d0lna2l0T2/kzQbcZW30PLza65r0qZxwnSdTvI9rSyiWnm5lraFuWYn77r7pSoufyJJdzazdkXbcxV2WYt+f1fvqkz/S+PnkxaerWJkmCTtkPSpPsXQD7l8jthifk+r3WYu11KRcYyZ3bLVgszsFBWndXRj1dv0gPW8rbv7z7V0ft99VNw1o9Gz9El3//Xqw5M0BNtdebmpD0m6V/nSnIoL/srdL1Nxm0apOFY1uwXYUKNgG0JmdmR5L7Wvq7gJslSM4ntUl6PRGst7oJmda2aPbLUTNLPdVWwAjY34ey0urNj46WHfsvcvNZd0jKT3m9my69OZ2W5m9lItDZnfKelVbZZ1dmX6beVfesuY2dMlvTowtn5/V+/W0kjFIyV91MxuVjia2ZNV9JY2/Ht54du6yOVzxBbre1r1NuPu12rpvpH7StpUXutrZWynSNqk5eeIdtLPbXpQVrutb65MVwuRnq+91kStt7ty1PEZWrpI79dUXGutei7ka7R0eaHHlvdGzQajROvpSWZ2z8rzPVWclHuYir8s/kLLT3z/pYrbfnytx/XtpuIG3idIutLMvqHiCtK/U3Gj98NV3FT5sLK9a/kFC6s+Lul+5fScmZ1VxtfYqH7t7uf0GGcvvqbir61HS7rIzM4p42ncKPovKm1f6e7tfjp5rYqfn/dRccLvhWb2QRUnBd9axY7kzuXzj0l6ZofY+vpduft1Znaiir/G91FxO62fmtnKmzdXd2IfVnFx19rI5XPEFvF76tc282otXVz1JEkPNLNzVVyPrbG9HKsi31+p4pJCIfq5TQ/Karf1D6vYH99aS8ft36joYeqLIdju3qDiCgZSMdL4b9192a887v57M3uDivvJSkUOhpyeMhwGcf8rHp0fWn5fwtDHdhW/248FLL/lvexUXL4jdJ2XS3pym/WMqdiYWs0/X2l7VLPXV/NZmi1X0hEqLvLYKqadkl4duO7HqOjNbLWsn6v4eWem8topq/2uQj53pd29A/Jpp4prh+3eZjld/f+oGD3caD/Th22iL5+jXFbH/49VbrOTq1xWT9tCv76nWNuMpJdp6WbtzR5Xq7jW42TltTMGGF/HbapTm8BldLWtt1jG61fM89p+5PEwbHeSXlRZziWSbtum7f4qrtXXaH9CjO8pxYOfRIfDNSrOE/ixipskz0h6uIqkfYm7L7aZtyN3/4qK8xE2qugW/66KE4X/JOkGFYXhp1Vc3uNod39vm2UtqvgL7KUq/hJuLCcZL84/u5+k56k4Ofb/VIyIukTFfefu6+4vClzWh1WMftusoji7XtIVKg4eL5Z0nHe+hVNjWVG+Ky96FO6o4urxn1Dx/3ejivNoLlBxnseEF5dZSPp/004unyO2GN9Tv7YZd3+5it7o90v6tYrrsf1BxV0IXi3pru5+VvCH7XN8g9Knbf2tK5738+fQXeq23ZU/mzd+1r5CRQF2Sav2Xpyq85rKS68akisadNS4uz0AAKgpMzte0hfLp1929+PbtUd+sqg6AQDIXPV+oiN5Z4NRRw8bAAA1Vl7I9lcqBgNcLulwX3HCPfJHDxsAAPX2ci3d7u6tFGujiR42AABqpLxW271VXErjQSoGmUnFSffj7v7bFrMiY1yHDQCAenmElq4l1uCSnk2xNrr4SRQAgPq6TNInJR3f7pJKyB8/iQIAANQcPWwAAAA1R8EGAABQcxRsAAAANZftKFEz4+Q8AAAwVNzdmr1ODxsAAEDNZdvD1tBsFKyZtXwvhpzXl/NnG/T6cv5sg15fzp9t0OvL+bPlvr6cP9ug1zeIdTXW0Qo9bAAAADVHwQYAAFBzFGwAAAA1R8EGAABQc9nemqpxWY9cP187gz7xE8OJPEEI8gShyJXVqXx/TUcfZD9KtNmoC5IJAACk0mlEaDP8JAoAAFBz2few0ZsGAADqpN01Yluhhw0AAKDmKNgAAABqLvufREcRPwMjBHmCEOQJQpErcdHDlqnZ2dnUIWAIkCcIQZ4gFLkSDwVbphYXF1OHgCFAniAEeYJQ5Eo8FGwAAAA1R8GWqTVr1qQOAUOAPEEI8gShyJV4uDUVAABAYp1uTUUPW6bm5uZSh4AhQJ4gBHmCUORKPBRsmVpYWEgdAoYAeYIQ5AlCkSvxcB02AOhCtzdt5rQMAP1ADxsAAEDNMeggU4uLixobG0sdBmqOPOleo4dtv5POatvu6jNPlJTHPog8QShypXcMOhhR27dvTx0ChgB5ghDkCUKRK/FQsGVq69atqUPAECBPEII8QShyJR4KNgAAgJqjYAMAAKg5CrZMTU1NpQ4BQ4A8QQjyBKHIlXgYJQoAXRjFUaIA4mOU6IiamZlJHQKGAHmCEOQJQpEr8VCwAQAA1BwFGwAAQM1RsGVqfHw8dQgYAuQJQpAnCEWuxMOgAwDoAoMOAMTAoIMRtWXLltQhYAiQJwhBniAUuRIPBVumtm3bljoEDAHyBCHIE4QiV+KhYAMAAKg5CjYAAICay37QQTO5fmYA8THoAMBqNfYjzTDoYMQsLCykDgFDgDxBCPIEociVeLIv2Nz9Zo9RMDc3lzoEDAHyBCHIE4QiV8L0UptkX7ABAAAMOwo2AACAmqNgy9S6detSh4AhQJ4gBHmCUORKPBRsmVq7dm3qEDAEyBOEIE8QilyJh4ItU7Ozs6lDwBAgTxCCPEEociUeCjYAAICao2ADAACoOQq2TE1MTKQOAUOAPEEI8gShyJV4sr81Va6fD0Aa3JoKQAyNfQu3phoxmzZtSh0ChgB5ghDkCUKRK/FQsGVqx44dqUPAECBPEII8QShyJR4KNgAAgJqjYMvU2NhY6hAwBMgThCBPEIpciYdBBwDQBQYdAIiBQQcjan5+PnUIGALkCUKQJwhFrsRDwZYpNhqEIE8QgjxBKHIlHgo2AACAmqNgAwAAqLlVDTows/0kPVHSoyXdWdKhknZK+o2kiyV9QdJn3P3rHZazl6SnSnqCpDtJOkTSHyT9VNLZkt7h7ld1GdtIDzrYvn271q5dmzoM1Bx50r1RHHRAniAUudK7ToMOdl/Fgp8k6Q2Smv3PjEk6WtJDJD1T0lFtlnNnSe+VdOyKtw4rH8dL+iczO9Hdv9xrvAAAAMOqp59EzezFKoqsRrH2E0mbJc1IOr2cPk/StR2Wc5SkT2upWLtC0tslnSbpzZIuK18/UtK5ZnZcL/GOos2bN6cOAUOAPEEI8gShyJV4uu5hM7NTJL2yfPprSc9090+2aLuXpHZF1ttV9KJJ0lckPdrdf1+Z/wWStkh6jIpeu/eY2V3d/aZu4wYAABhWXfWwmdkRkt5YPv2tpPu2KtYkyd1vaHX+mpk9VNKDy6d/kPS4arFWzn+dpBMlXVS+dCdJJ3cTMwAAwLDr9ifR50u6VTn9Ane/ZBXrfmplepO7/6ZZI3e/VtLrW8yHFiYnJ1OHgCFAniAEeYJQ5Eo8waNEzeyWkrZLOkDFuWaHuvuNPa3UbDcVvWqN4u84d/9um/YHq+jRMxWjUA919991WMdIjxIFEMcojhIFEF8/b011HxXFmiR93d1vNLPbmdlrzewCM7vGzK4sp99sZndts6zba6lYW5T0/XYrLn8q/XEl5nbLhqTZ2dnUIWAIkCcIQZ4gFLkST7cFW8NPzezJkn4g6YWSjpG0j6T9y+lnS/qumc2WvWkrVS/h8TN33xmw/osq0ysvAYIVFhcXU4eAIUCeIAR5glDkSjzdjBL988r0PVUUZburuPTGRyT9StLBkk5QUbSZinPe9pO0ccWyDqtMXxq4/ur5coe1bLVCo4sxBD9dAKOrm30FAITo536lm4LtwMr0fct/3y/p6e5+TeMNM3uhpBdLekX50gYz+6i7f7wy/1hl+hqFqV7Tbaxlq1WYmZmRJK1Zs0YbN27U3NycFhYWdr0/PT2t7du3a+vWrbtem5qa0sTExK55JWl8fFzr16/Xli1btG3btmXLX1hY0Nzc3K7X1q1bp7Vr1y7rRp6YmNDU1JQ2bdqkHTt2SJLGxsY0PT2t+fn5ZTfX3bBhg6Tl176ZnJzUmjVrNDs7u+uvnRw+0+TkJJ+pz5/pwAMPXDZ/Dp9pNf9PMSwsLAx97pnZsvhT/z/lmHu5fKZvf/vbu9aXy2dazf9TP3Uz6OAzkv668tIPJd3D3f/Yov27tDSi8yvu/oDKe6dJenn59Ex37zjy08xeruKCupL0dnd/Zof2DDoAEKzbwQQMOgDQT/0cdHDdiueva1WslV5Zmb6fmVV76KrL2jNw/XtXptveQQFa9hcF0Ap5ghDkCUKRK/F0U7CtPJPwc+0au/uFKi4DIhXns1XveFBd1r6B69+nTSxYodqtDLRCniAEeYJQ5Eo83RRsOyrTrqVirJ3qgIJbV6Yvq0wfEbj+arumF9kFAADIUTcF248q014+OmnV5oLK9O1bXPpjpaNbxAIAAJC1bgYd3FXS9yovHeHuv+4wz6WSDi+fPsjd58vXV97p4O7u/r2bL2HXcg6S9Dtxp4Ngi4uLGhuLMpgWGSFPljDooDXyBKHIld71bdCBu39f0s8qL/11q7blise1VKzdJOk7lWXtlPSxSvMndFj941QUa5L05U7FGqTt20N+scaoI08QgjxBKHIlnm5v/v7OyvQLzGyPNm1fUpn+jLtfteL9MyvTp5rZbZotpLyH6QsrL707KNIRV71uDdAKeYIQ5AlCkSvxdFuw/Yekxs+gd5F0ppktG+VpZruZ2b9o6RpsO7V0zbVd3P1TkubLpwdL+lD502d1WbeU9B4t3WXhJ5Le1WXMAAAAQ62bOx3I3a81s3WSzlNxXbQnSforMztHxa2jDpb0CBW3pmr4V3f/aotFPkPS1yTdRtIDJF1kZh+U9EsVt596nKQ1ZdurJT3F3f/UTcwAAADDrquCTZLc/Utm9kgVPV1HqCisnt2k6fWSXuTu/9lmWT83s4dKeq+KIu9ASc9q0vRSSSe6Oxd4CTQ1NZU6BAwB8gQhyBOEIlfi6bpgkyR3/5yZ3UnSKZL+TtIdJR2i4oK2P1fRA/fmTqNIy2V938yOk3SyisEHd1JxzbbLJV0k6WxJb2tyDhzamJiYSB0ChgB5ghDkCUKRK/F0ew7bLu6+6O7/7e4PdvfD3X1Pdz/Y3e/l7i8JKdYqy7rB3Te7+0PLZe3l7oe5+wPc/Q0Ua92r3hgXaIU8QQjyBKHIlXh6LtgAAAAwGBRsAAAANUfBlqnx8fHUIWAIkCcIQZ4gFLkST/CtqYbNqN+aCkB3uDUVgJT6dmsqDJctW7akDgFDgDxBCPIEociVeCjYMrVt27bUIWAIkCcIQZ4gFLkSDwUbAABAzVGwAQAA1ByDDgBADDoAkBaDDkbUwgK3XUVn5AlCkCcIRa7EQ8GWqbm5udQhYAiQJwhBniAUuRIPBRsAAEDNUbABAADUHAVbptatW5c6BAwB8gQhyBOEIlfioWDL1Nq1a1OHgCFAniAEeYJQ5Eo8FGyZmp2dTR0ChgB5ghDkCUKRK/HsnjqA2BrXNaniukjA8Gu2bbfDdg+gLrrdf0n0sAEAANRe9j1so/pX9cTEROoQMARyyJPQOw6gdznkCQaDXAnTrDbp1OtGD1umpqamUoeAIUCeIAR5glDkSjwUbJnatGlT6hAwBMgThCBPEIpciYeCLVM7duxIHQKGAHmCEOQJQpEr8VCwAQAA1BwFW6bGxsZSh4AhQJ4gBHmCUORKPJbrKEozc2l0R4kCuWuMqAodJdppX9Dt8vq1XgCQlvZB7t50uCg9bJman59PHQKGAHmCEOQJQpEr8VCwZYqNBiHIE4QgTxCKXImHgg0AAKDmKNgAAABqjoItUxs2bEgdAoYAeYIQ5AlCkSvxULABAADUHAVbpjZv3pw6BAwB8gQhyBOEIlfioWADAACoOQo2AACAmqNgy9Tk5GTqEDAEyBOEIE8QilyJh4ItU2w0CEGeIAR5glDkSjwUbJmanZ1NHQKGAHmCEOQJQpEr8VCwZWpxcTF1CBgC5AlCkCcIRa7EQ8EGAABQcxRsmVqzZk3qEDAEyBOEIE8QilyJx9w9dQxRmJlLUq6fDxh1ZiZJ2u+ks9q2u/rMEyV13hd0u7x+rRcApKV9kLtbs/fpYcvU3Nxc6hAwBMgThCBPEIpciYeCLVMLCwupQ8AQIE8QgjxBKHIlHgo2AACAmqNgAwAAqDkKtkxNT0+nDgFDgDxBCPIEociVeLIv2MzsZo9RsH379tQhYAgMKk+abYftHqgX9icIRa6E6WW/l33BNqq2bt2aOgQMAfIEIcgThCJX4tk9dQCxcQ0koD5Cr10GADlrVpt06mWjhw0AAKDmKNgyNTU1lToEDAHyBCHIE4QiV+KhYMvUxMRE6hAwBMgThCBPEIpciYeCLVMzMzOpQ8AQIE8QgjxBKHIlHgo2AACAmqNgAwAAqDkKtkyNj4+nDgFDgDxBCPIEociVeCjYMrV+/frUIWAIkCcIQZ4gFLkSDwVbprZs2ZI6BAwB8gQhyBOEIlfioWDL1LZt21KHgCFAniAEeYJQ5Eo8FGwAAAA1R8EGAABQcxRsmeLihQhBniAEeYJQ5Eo8FGyZWlhYSB0ChgB5ghDkCUKRK/FQsGVqbm4udQgYAuQJQpAnCEWuxEPBBgAAUHO7pw4AQP2YWeoQAAAV9LBlat26dalDwBAgTxCCPEEociUeetgytXbt2tQhYAh0ypP9Tjqr7ftXn3liV+0wnNifIBS5Eg89bJmanZ1NHQKGAHmCEOQJQpEr8VCwAQAA1BwFGwAAQM1RsGVqYmIidQgYAuQJQpAnCEWuxEPBlqmpqanUIWAIkCcIQZ4gFLkSDwVbpjZt2pQ6BAwB8gQhyBOEIlfioWDL1I4dO1KHgCFAniAEeYJQ5Eo8FGwAAAA1R8GWqbGxsdQhYAiQJwhBniAUuRJP9nc6aHZPRHdPEMlgTU9Ppw4BQ6CueZLTvUxDP0ud90t1zRPUD7kSppd9HD1smZqfn08dAoYAeYIQ5AlCkSvxZN/DVue/WmOan5/X5ORk6jBQc3XNk073JpWG5/6kOdxnta55gvohV8I0q0069brRwwYAAFBzFGwAAAA1R8GWqQ0bNqQOAUOAPEEI8gShyJV4sj+HDcCSnEZfAsAooYctU5s3b04dAoBMsD9BKHIlHnrYgBGUw8hFABgl9LABAADUHAVbprgODoB+YX+CUORKPBRsmWKjAdAv7E8QilyJh4ItU7Ozs6lDwACZWdBjlPHd9I79CUKRK/FQsGVqcXExdQgAMsH+BKHIlXgYJQpkhNGfrfHdABhm9LBlas2aNalDAJAJ9icIRa7EQ8GWqY0bN6YOAUAm2J8gFLkSDwVbpubm5lKHACAT7E8QilyJh4ItUwsLC6lDAJAJ9icIRa7EQ8EGAABQc30r2MzsrmZ2o5l55XFU4Lz7m9k/mtkXzWyHmd1gZpea2Xlm9kwz27NfcQIAAAybvlzWw8x2l3SGpD16mPeBkt4j6bYr3jq8fDxU0j+Y2ZPc/YJVhjoypqenU4cAIBPsTxCKXImnXz1sL5Z0nKQbu5nJzO4p6VwtFWuXSfofSadJerukK8rX7yzp02Z2ZD+CHQXbt29PHQKATLA/QShyJZ5VF2xmdhdJLy2f/nsX8+2uomdt3/KlcyTd3t2f6+6vcPdnSjpa0lfL99dKeutq4x0VW7duTR0CgEywP0EociWeVRVsZdH1Tkl7SvqhpFd1MfvTJN2xnL5I0lPc/bpqA3f/vaTHSbq8fOlhZvbg1cQMAAAwbFbbw/bPkiYk7ZT0DHf/YxfzPrUy/Tp3v7ZZI3e/TNKmFvMBAABkr+eCzcyOlfSy8ul/uvs3u5j3EEn3K5+6pA92mOUDlelHmxmXI+lgamoqdQgAMsH+BKHIlXh6KnzM7BZa+in0F1o6hy3U3SrrvsDd/9Ch/XclLZbTB0g6qsv1jZyJiYnUIQDIBPsThCJX4um1p+qFku5VTj+r1c+ZbRxbmf5pp8buvlNFYdhsfjQxMzOTOgQAmWB/glDkSjxdF2xmdoykmfLpO9z9sz2s97DK9KWB81zSYn4AAICsdXXh3MpPoXupuGZar1fIG6tMXxM4T7UXb6xlqxXMLLSp3D24LQD0E/sqID/dbNeddHung2lJf1lOP8fdr+hxvbesTIdebPf6yvQ+Pa63rUZX7po1a7Rx40bNzc0tu5Ht9PS0tm/fvuw6M1NTU5qYmFjWDTw+Pq7169dry5Yt2rZt27LlLywsaG5ubtdr69at09q1azU7O7vrtYmJCU1NTWnTpk3asWOHJGlsbEzT09Oan5/X/Pz8rrYbNmyQJG3evHnXa5OTkxofH9fs7KwWFxez+UyTk5N8pg6fCaNhbm5uoLm33377LcuzUdme+Ezdf6Yf/OAHu9aXy2dazf9TP1noX2pmdkcVJ//vLelsd39ci3bVBf6Zu1/cpM2bJD2nfPoad39RwPrfL+kJ5dNpd39Dh/Yu8ZcoRkPjr7j9Tjqrbburzzwxi3bDEGPMz8x+DchPYz/u7k275YLOYSsvo/FOFcXaFVoqtnq1WJnet2Wr5aq9aostW0GStGXLltQhAMgE+xOEIlfiCR108DRJ9y2np8uL2a5Gdf4jAueptvvNKtefvWq3MACsBvsThCJX4gkt2G5bmX67mXmrx4r5flF5b6by+gWV6T/vtPKyh+/2lZd+FBg3AADA0Ot20EG/fFfFHQ5M0jFmdqC7X96m/d20NDL0Si2/JhuQrX6OMAIADK/Qgu27kt4V2PbkyvSHJF1dWYYkyd1/a2ZflXR/Fb18j5P0tjbLfHxl+iPlhXTRBiMHAfQL+xOEIlfiCSrY3P3Dkj4c0tbMqgXbC5qNEi2dqaJgk6QXmtlZ7n5dk+UdKqk6NvbdIXGMuoWFBW4RkpHQkYZADOxPEIpciSflTdTfoaXbUo1LOtPMqtdnk5kdpOLG8AeXL32mxzsrjJzqdWcAYDXYnyAUuRJPqnPY5O5/NLOnSPq8ikt2PE7S/czsQypGkR5VvnZgOctlkp6VIFQAAICkkhVskuTu3zSzR0p6j6TDJa2R9NwmTS+Q9KQ2P68CAABkK2nBJknuPm9mx0p6hqTHSjpaRa/a71QUah+Q9C53vyFdlMNn3bp1qUNAC4z8xLBhf4JQ5Eo8fS/YWt1SocM8V0p6Q/lAH6xduzZ1CAAywf4EociVeJL3sCGO2dlZhlfXXDf3jgRSYn+CUORKPClHiQIAACAABRsAAEDNUbBligsXAugX9icIRa7EwzlsmZqamkodAoBIQkcau3tf1sf+BKHIlXjoYcvUpk2bUocAIBPsTxCKXImHHrZM7dixI3UIACIZ9L1l2Z8gFLkSDz1sAAAANUfBlqmxsbHUIQDIBPsThCJX4qFgy9T09HTqEABkgv0JQpEr8VCwZWp+fj51CAAywf4EociVeCjYMsVGA6Bf2J8gFLkSDwUbAABAzVGwAQAA1BwFW6Y2bNiQOgQAmWB/glDkSjwUbAAAADVHwZapzZs3pw4BQCbYnyAUuRJP9remanaT5H7dEBkA6mzQN4kHECZ026yihw0AAKDmsu9hG9W/HCcnJ1OHACCxft0knv0JQpErYZrVJp163ehhyxQbDYB+YX+CUORKPBRsmZqdnU0dAoBMsD9BKHIlHgq2TC0uLqYOAUAm2J8gFLkSDwUbAABAzVGwZWrNmjWpQwCQCfYnCEWuxEPBlqmNGzemDgFAJtifIBS5Eg8FW6bm5uZShwAgE+xPEIpciYeCLVMLCwupQwCQCfYnCEWuxEPBBgAAUHMUbAAAADVHwZap6enp1CEAyAT7E4QiV+KhYMvU9u3bU4cAIBPsTxCKXImHgi1TW7duTR0CgEywP0EociUeCjYAAICao2ADAACoOQq2TE1NTaUOAUAm2J8gFLkSDwVbpiYmJlKHACAT7E8QilyJh4ItUzMzM6lDAJAJ9icIRa7EQ8EGAABQcxRsAAAANUfBlqnx8fHUIQDIBPsThCJX4qFgy9T69etThwAgE+xPEIpciYeCLVNbtmxJHQKATLA/QShyJR4Ktkxt27YtdQgAMsH+BKHIlXgo2AAAAGqOgg0AAKDmKNgyxcULAfQL+xOEIlfiyb5gM7ObPUbBwsJC6hAAZIL9CUKRK2F6qU2yL9hG1dzcXOoQAGSC/QlCkSvx7J46gNjcPXUIAAAAuzSrTTr1stHDBgAAUHMUbJlat25d6hAAZIL9CUKRK/FQsGVq7dq1qUMAkAn2JwhFrsRDwZap2dnZ1CEAyAT7E4QiV+KhYAMAAKg5CjYAAICao2DL1MTEROoQAGSC/QlCkSvxULBlampqKnUIADLB/gShyJV4KNgytWnTptQhAMgE+xOEIlfioWDL1I4dO1KHACAT7E8QilyJh4INAACg5ijYMjU2NpY6BACZYH+CUORKPBRsmZqenk4dAoBMsD9BKHIlHgq2TM3Pz6cOAUAm2J8gFLkSDwVbpthoAPQL+xOEIlfioWADAACoOQo2AACAmqNgy9SGDRtShwAgE+xPEIpciYeCDQAAoOYo2DK1efPm1CEAyAT7E4QiV+KhYAMAAKg5CjYAAICao2DL1OTkZOoQAGSC/QlCkSvxULBlio0GQL+wP0EociUeCrZMzc7Opg4BQCbYnyAUuRIPBVumFhcXU4cAIBPsTxCKXImHgg0AAKDmKNgytWbNmtQhAMgE+xOEIlfi2T11ALGZ2c1ec/cEkQzWxo0bU4eQjWY51M4o5BdGC/sThCJXwnR7XJHoYcvW3Nxc6hAAZIL9CUKRK/Fk38M2qr0dCwsLmpqaSh1GVvY76ay271995okDigQYLPYnCEWuhGlWm3TqdaOHDQAAoOYo2AAAAGqOgi1T09PTqUMAkAn2JwhFrsRDwZap7du3pw4BQCb2339/mVnwA6OLY088FGyZ2rp1a+oQAAAjhmNPPNmPEgUA9AcjpYF06GEDAACoOQq2THEdHADAoHHsiYeCLVMTExOpQwAAjBiOPfFQsGVqZmYmdQgAgBHDsSceCjYAAICao2ADAACoOQq2TI2Pj6cOAQAwYjj2xEPBlqn169enDgEAMGI49sRDwZapLVu2pA4BADBiOPbEw50OMrVt27bUIdQe9zwECmwL6BeOPfHQwwYAAFBzXfWwWfFn2H0kPUTSfSXdSdJtVBR+l0v6kaTPSXqHu1/WxXJ3k/RESesk3V3SoZKukvQLSR+V9DZ3/003sQKhuD8iRh3bAFB/wQWbmT1K0lskrW3R5LDy8deSTjOz09z99QHLPULSeyXdf8Vbh5SPe0t6vpk9093PCY131HHxQgDAoHHsiaebn0TvpOXF2s8kvVvSKyS9RNLbJDV6wfaW9Dozm223QDM7UNJ5WirWrpN0lqSXSfoPST8vXz9I0vvN7G+6iHekLSwspA4BADBiOPbE0+05bIuS3iDpLu5+tLuf7O6nufsr3f1Zko6StKnS/vlm9qA2y3u9pGPK6Z9IupO7P8Xd/83dny9pXNKby/d3l/RuM9u/y5hH0tzcXOoQAAAjhmNPPN0UbOdKOsrdp939h80auPv17n6qpE9UXn5us7Zm9heSTimf/lHS49394hXLu0nS8yR9pXzpEEkv6CJmAEDNmFlXDwBdFGzu/n13/0Ng8/+qTN+3RZsTK+v/oLv/qMV6d0r6t8pLJwXGAAAAkIVY12H7RWX64BZtHl2Zfn+H5X1a0hWSDpB0lJkd5+7n9xzdCFi3bl3qEACgLUan5odjTzyxrsN2eGX6tyvfNLO9VAxiaPhyu4WVvWxfrbx03KqiGwFr17YazAsAQBwce+KJVbA9sTL9pSbv/4WkW5TTV7j77wKWeVFl+theAxsVs7NtB+gCANB3HHvi6XvBZmbHSnp65aVNTZodVpm+NHDRl7SYHwAAIGt9PYfNzPaV9B5Je5QvnePu802ajlWmrwlc/LUt5u8UU2hTuXtwWwAAgHb6Ocq5bwVbeXupM1XcWkoqes6e1aL5LSvTNwau4vrK9D5dBReocYXmNWvWaOPGjZqbm1t2EcDp6Wlt375dW7du3fXa1NSUJiYmll3deXx8XOvXr9eWLVuW3Qh3ZmZGCwsLy65Ts27dOq1du3ZZN/LExISmpqa0adMm7dixQ5I0Njam6elpzc/Pa35+flfbDRs2SJI2b96867XJyUlNTExodnZWi4uL2XymycnJvn6mmEI+E5Cr6nbWbB/RrVT7iBz3e7E/04UXXrhrfbl8ptX8P/WT9aNXqbzH6GZJzyxfukrSg9296SWPzezxkj5QPv2Gu98nYB1/L+l/yqdz7v6oDu1dotcMrTX+8gkdqRbarlPOha63l3WPSrthiHGUP3O/toHQ5QE5aGwX7t60W65f57C9SUvF2qKkv2lVrFXaNOwbuI5qr9piy1aQJG3a1OzUQQAA4uHYE8+qCzYz+29Jf18+vVrSCe7+tQ6zXVaZPiJwVdV2v2nZCpK0q1sXAIBB4dgTz6oKtrJYa9x66hpJj3D3r7SZpeFCSTvL6QPM7NYB8xxdmW56VwQAAIAc9VywmdmbtFSsXSvpke7e7JprN+Pu10u6oPLS/TusazdJ96u8xF0OOhgbCx5ICwBAX3Dsiaengs3M/kfSc8qnjWLtC10u5iOV6Sd0aPvXkg4sp3/p7t/pcl0jZ3p6OnUIAIARw7Ennq4LtrJYa5yzdq2kv21xrbVOztLSz6JPMLNjWqxvN0mnVV46s4d1jZzqMGQAAAaBY088XRVsLYq1z/eyYnf/saR3l0/3lPQhM7vdivXdQtJ/Sjq+fOl3kl7fy/pGDRsNAGDQOPbEE3zhXDP7Fy0Va5L0WUkTZjYRMPv73P2SJq+/QNJ9Jd1R0jGSLjCzs1UMSjhQ0qMl3aFs+ydJp7j7laExAwAA5KCbOx2Mr3g+VT5CfFvL7wUqSXL335vZwyS9V0Xhto+kpzSZ/3JJz3L3j4eHCwAAkIe+3ku0F+7+KzN7gKQnSVqn4tZWh6q4OO4vJH1U0lvd/bKWC8HNNG6fgcHr573jgGHENjC6OPbEE1ywufspkk6JEYS775S0tXwAAACgInkPG+LYvHnzwG50juVC748I5IptYHRx7ImnX/cSBQAAQCQUbAAAADVHwZapycnJ1CEAAEYMx554KNgyxUYDABg0jj3xULBlanZ2NnUIAIARw7EnHgq2TC0uLqYOAQAwYjj2xEPBBgAAUHMUbJlas2ZN6hAAACOGY088FGyZ2rhxY+oQAAAjhmNPPBRsmZqbm0sdAgBgxHDsiYeCLVMLCwupQwAAjBiOPfFQsAEAANQcBRsAAEDNUbBlanp6OnUIAIARw7EnHgq2TG3fvj11CACAEcOxJ57sCzYzu9ljFGzdujV1CACAEcOxJ0wvtUn2BRsAAMCw2z11ALG5e+oQAAAAdmlWm3TqZaOHLVNTU1OpQwAAjBiOPfFQsGVqYmIidQgAgBHDsSceCrZMzczMpA4BADBiOPbEQ8EGAABQc9kPOkA+ur0kCwNOAAC5oIctU+Pj46lDAACMGI498dDDlqn169enDiGa/U46q+37V5954oAiAQBU5XzsSY0etkxt2bIldQgAgBHDsSceCrZMbdu2LXUIAIARw7EnHgo2AACAmuMcNgBArXU7QjwEo8gxbOhhyxQXLwQADBrHnnjoYcvUwsICtwgBkIVOI8OlpdHhjCJPi2NPPPSwZWpubi51CACAEcOxJx4KNgAAgJqjYAMAAKg5CrZMrVu3LnUIAIARw7EnHgq2TK1duzZ1CACAEcOxJx4KtkzNzs6mDgEAMGI49sRDwQYAAFBzFGwAAAA1R8GWKS5cCAAYNI498XCng0xNTU0lXX8v9/7r9739Ytx/EEAeQvcP3HO0O6mPPTmjhy1TmzZtSh0CAGDEcOyJhx62TO3YsSN1CJK6uwfgoNfNPQWB0cX+IY66HHtyRA8bAABAzVGwZWpsbCx1CACAEcOxJx4KtkxNT0+nDgEAMGI49sSTfcFmZjd7jIL5+fnUIQAARgzHnjC91CbZF2yjio0GADBoHHviyX6UKNfQAQAAddKsNunUy0YPGwAAQM1RsGVqw4YNqUMAAIwYjj3xULABAADUHAVbpjZv3pw6hK41GzUzaqN7AWCYDeOxZ1hQsAEAANRc9qNEMTy4tx8AAM3Rw5apycnJ1CEAAEYMx554KNgyxUYDABg0jj3xULBlanZ2NnUIAIARw7EnHgq2TC0uLqYOAQAwYjj2xEPBBgAAUHMUbJlas2ZN6hAAACOGY088FGyZ2rhxY+oQAAAjhmNPPBRsmZqbm0sdAgBgxHDsiYeCLVMLCwupQwAAjBiOPfFQsAEAANQcBRsAAEDNUbBlanp6OnUIAIARw7EnHgq2TG3fvj11CACAEcOxJx4Ktkxt3bo1dQgAgBHDsSceCjYAAICao2ADAACoOQq2TE1NTaUOAQAwYjj2xEPBlqmJiYnUIQAARgzHnngo2DI1MzOTOgQAwIjh2BMPBRsAAEDNUbABAADU3O6pA4jNzG72mrsniGSwxsfHJTX//O2MwncDAIijcexBe90emyV62LK1fv361CEAAEYMx554su9hG9Ueoy1btizbcPY76ay27a8+88TYIQEAMrfy2IPmmtUmnXrd6GHL1LZt21KHAAAYMRx74qFgAwAAqDkKNgAAgJrL/hy2XIWMMDn99NMHEAkAoJNRGbHPhXPjoYcNAAD0xcLCQuoQskUP25Bj9CcADI/c99lzc3PcTzQSetgAAABqjoINAACg5ijYAABAX6xbty51CNmiYAMAAH2xdu3a1CFki4INy5hZ0AMARkGqfeKw7otnZ2dTh5AtCjYAAICa47IeWCb3IecA0I1U+0T2xViJHjYAANAXXIMtHgo2AADQF1NTU6lDyBYFGwAA6ItNmzalDiFbnMO2Cr2M0BnWG/oCAFqr44jNFHbs2JE6hGzRwwYAAFbNzHT66aenDiNbtelhM7PbS3qmpEdIOlLS3pJ2SPqGpPe4+7kJw2ur02geiRE9AJAzRnUitloUbGb2HEmvk3TLFW/dvnysM7MPSzrZ3a8acHgAAABJJS/YzOxUSW+qvPR9SZ+UdK2ku0maUhHnYySdY2YnuPuNg44TAAAglaQFm5ndQdJ/Vl56sbu/ekWb4yR9QtKhkh4saVrSqwYWJAAAQGKpBx28XNKe5fSWlcWaJLn7+ZKeWnnpRWZ2wABiS2JY7h9X9/gAYBSEHjO6eaSKLYVhiLEhWcFmZvtJemz51FUUb025+3mSvl4+3V/Fz6OrWXfyLx4AAPRu1I7lKX8SfZiKkaCS9H13v7BD+w9Iuk85/XeSzogUVy3UfcRR3eMDgFHQzVUKuC9qa8MQY8qfRO9Rmf5yQPsvVaaP63MsAAAAtZWyYDu2Mv3TgPYXVaZva2b79zkeAACAWkpZsB1Wmb60U2N3v1zSNZWXDu17RAAAADVkqe5taWY/1FIv2wnu/smAef5P0iHl0wl3/06btty0EwAADBV3bzqSImUPW/WuBqEXwr2+Mr1PH2MBAACorZSjRK+rTO/ZstVye1emr23XsFWFCgAAMGxS9rAtVqb3DZyn2qu22LIVAABARlIWbJdVpo/o1Li8u0G1sPtNvwMCAACoo5QF2wWV6T8PaH90ZfpSd7+qz/EAAADUUsqCrTrC8/4B7Y+vTJ/f51gAAABqK2XB9iktjfq8m5mNd2j/+Mr0OXFCAgAAqJ9kBZu7Xy3pw+VTk3Raq7Zm9hBJ9yufLlbmAwAAyF7KHjZJepmkP5bTTzGzF65sYGZ3k/TuykuvKe96kBUzu72ZvdLMvmtmfzCza83sZ2a2xcweMcA43mdmXnmcMah1I8wgc2VFLoQ8ntvP9aN3KfcpZjZhZq81s++Y2W/M7AYz22Fm55vZ283sKWZ2UMwYEG5QuWJmF/ewT2k8jupXHMMq2Z0OdgVg9hxJb6q89D1Jn1RxnbW7SZqStEf53hckPdzdbxhokJGV38HrtPxiwit9WNLJMQdbmNljJX1oxcvvcvdTYq0T3Rl0rvRwx5DnufubOjdDTKn2KWZ2G0lvlLQuoPnT3P2Mfq0bvRlkrpjZxZJu18Os10o6tPxlbmQlL9gkycyeJ+m1Wn5h3JXmJJ3k7lcOJqrBMLNTJb2l8tL3dfOCtXGB48+puI1X6J0huonjYBUjd2+j4s4TjYsZU7DVRIpcqRRsl0t6ZcAsn3H3765mnVidVPsUMztS0me0NOr/GhV/ZP9A0pWSDpa0RtI9JY2Lgi25QeeKmT1L0q0Cmz9P0pHl9Lvd/eRe15sNd6/FQ9IdJL1aRQ/b5SruhPALSe+V9MjU8UX8zDdI8vLxoiZtjlNxzbpGmxdHimVrufz/k/TflfWdkfp74pEuVyrLujj1d8Cj1nmyt6QfV5a5SdLBbdr/maQjU39fo/yo0/GnyXoPU3G6VGO9D0z9fdXhkTyAUX5IOquSkGe1afewSrsrJR3Q5zgeU1n+OkkzFGz1eqTKFQq24XokzJPZyvJenfp74FHfXAmM7Z8r67ww9XdVl0fqQQcjy8z2k/TY8qlLenmrtu5+nqSvl0/3V1Fg9SuOgyT9b/n0Y+6+tV/LRn/UJVdQb6nyxMwOkfSc8ulPJb2012VhMIZgn/KMyvTbB7C+oUDBls7DtHTO3vfd/cIO7T9Qmf67PsbxX5IOlXSVpGf3cbnon7rkCuotVZ48Q9Je5fRb3f1Pq1gWBqO2+xQze6CWzoP8k6R3xVzfMKFgS+celekvB7T/UmX6uH4EYGZTkk4sn/6zu1/aj+Wi75LniqSDzOzj5aUZbjSzy83sx+UlGk7o0zqwOqnypPr//0UrPMnMPlnmS+OSHp8ys2eb2V4tl4RBqcM+pZVq79rH3J37hpco2NI5tjL904D2F1Wmb2tm+69m5WZ2oIoTgyXpi5Vp1E/SXCmNSXqEipOB95B0gKS/kPR0Seea2bfMLOSewIhn4HliZrdQMeqz4beSzlMxWOzhKvJlz/Lfh0l6s6SfWHF9TaRTh33KzZjZrbT8rkb8HFqxe+cmiOSwynTHni13v9zMrpG0b/lS42fMXv2niiH210t6lpdneqKWUufKTknfUnH/3x3l88MlTUo6pmxzT0nfMrMHuPsPV7Eu9C5FnhwuaZ/K8zO0dN/nz0r6qorRfneR9CgVP50epaIn7l7uvq3L9aE/Uu9TWlmnpXz6taRPRFjH0KJgS2esMn1N4DzXammDGWvXsB0z+1tJJ5VPT2enWXvJckXSiyS9093/r9mbZvYoSW+TdIiK6yt91Mzu5O7XN2uPqFLkyYErnh+v4kD+WHf/bPUNM7uDpHNVXINtfxWjFO/Vwzqxein3Ke1Ufw49w91virSeocRPoulUryodeiHC6kFwn5at2jCzA7T08+f5kl7fy3IwUElyRZLc/TWtirXy/Y9KeoiWdvp/Jmljr+vDqqTIk2YH7lNXFmuS5O4/k/ToSmz3tOI+0Ri8ZPuUVszsrlr6ed0lvaPf6xh2FGzpXFeZ3rNlq+Wqd4K4tsf1/oektSpG3zyDEV1DIVWuBHH376v4ib3hiTHXh5ZS5Ml1K57/XMX5a025+08knV156W97WCdWr477lGdWpj/v7j+PsI6hRsGWzmJlet+WrZar/lWz2LJVC+VovlPKp6939/O7XQaSGHiu9OD9lel7mxn7lsFLkScr5/l8wPmwn6tM36NlK8RUq31KOXL4xMpLDDZogp1qOpdVpo/o1Lj8KbO6YfUy1Hlz+e9PJZ3ew/xII0WudKs6imx3FfeNxGClyJPLVjwPuTRQtc2te1gnVq9u+5S/k3RQOX25lvfCokTBls4FlemQyyEcXZm+1N17GaHT2DD/XNJ1ZubNHpL+tTLPySvew+ClyJVukRvpDTxPynkuqb4UMlu360Hf1W2fUv059D0MWmqOgi2d71Sm7x/Q/vjKND9ljpZhyJXqDv0mSb8f0HqxJFWefLsy3bG3ZkWblgNaEFVt9ilm9meSHlx5iZ9DW+CyHul8SsWom70l3c3MxjtcXqN6McFzelxn6C0+7i6pcWHLnynsStiIJ0WudOsJlelvufvOAa0XS1Llydlaul3Rg83MOpzHVj04f3MV60Xv6rRPeZokK6e/7e7f6/Py85H67vOj/JC0VcXPAy7pzDbtHlJpd5WkAyPHNVNZ3xmpvyce9c2Vcp3HqjgJubHe56f+vkb1kSJPVFza4w+V5T25Tds7Srqh0vY+qb+zUX3UYZ+i4le+X1WWf2rq76XOD34STetlKq4CLklPMbMXrmxQ3sLl3ZWXXuPulzdpN8m5ZlkbeK6Y2bfN7AlmtnebNo9UMepvv/KlX0l6S+ePg0gGnifuvqjlg5g2mdmDV7YrL5z7ES1dRuJT7v71jp8IsdTh+PNwSbctp6+VtKWLeUeOlVUuEjGz50h6U+Wl70n6pIrkvZukKRX3bpSkL0h6uLvf0GQ5k5I+33ju7rayTRcxzWhp4MG73P2UXpeF/hl0rlR2vFdL+rqkH6o4N22nimv5PUjSnSqzXCnpr5yfNJJKsU8p7yl6torbTzV8RsWtqf4k6c7le43i/9eS7unuK0eZYoBSH3/M7IOSHlc+5VjTAeewJebu/1Nes+q1Ks8n0NL5Y1Vzkk5qtrFgNCTMlf1U/CzS7qr0X5V0srtf1KYNBiBFnrj7TWb2JBU3d39a+XKrnPmmpMdTrKWX8vhjZodoeYHPYIMOKNhqwN3/28zOlfQsSSdIOlLFxnOZpG+oOL/g4wlDRE0MOFeOk3Sf8nEnFfcLvbWKG3hfIemXKnre3u/uX+rTOtEHKfYpXlyK4elm9lZJJ0uaVNETu4eK0aDfkPQ+SWc7P+3URsLjz0la6r27kH1IZ/wkCgAAUHMMOgAAAKg5CjYAAICao2ADAACoOQo2AACAmqNgAwAAqDkKNgAAgJqjYAMAAKg5CjYAAICao2ADAACoOQo2AACAmqNgAwAAqDkKNgAAgJqjYAMAAKi5/w88V+hP2ptvJAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 262,
       "width": 310
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Testing posterior simulation.\n",
    "NUMBER_OF_SIMULATIONS = 1000\n",
    "prob_list = []\n",
    "x = np.array([1, 20, 12, 8, 43, 0, 2])\n",
    "for _ in range(NUMBER_OF_SIMULATIONS):\n",
    "    prob_list.append(simulate_posterior(x, w_map, -inv_hessian))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist(prob_list, bins=50, edgecolor='black')\n",
    "ax.set_title('Distribution of Probability of $x$')\n",
    "fig.set_size_inches(5, 4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given we can simulate the probability for a woman working with a feture vector $x$, we can use this for other purposes. For instance, we can examine the probability of $k$ women working out of $n$ given $K\\sim \\text{Bin}(n, p)$ in which $p$ is the previously ascertained probability. "
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit",
   "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.9.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}