Assignment2_KonstantinosGeorgiou.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"KonstantinosGeorgiou_Assignment2_MATH525.ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyPjXx/RKXSdEYJqFTrmydg6"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","metadata":{"id":"x1wTmc4rSmRU"},"source":["# MATH 525 - Statistics I, Fall 2021\n","## Assignment 2 - Konstantinos Georgiou\n","Due: Monday, October 10, by 8:00PM"]},{"cell_type":"markdown","metadata":{"id":"uY7BDsyhTodf"},"source":["## Initializations"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"OyNhoKzgTgu0","executionInfo":{"status":"ok","timestamp":1634026866488,"user_tz":240,"elapsed":26983,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"7eff71f3-87f0-45e9-f2d1-6099a1ddc59c"},"source":["# (For Google Collab) Mount Google Drive and set working directory\n","from google.colab import drive\n","import os\n","# Mount Google Drive\n","drive.mount('/content/drive')\n","# Set working directory\n","project_path = '/content/drive/MyDrive/Documents/Work Documents/UTK/Courses/MATH525/Assignments/Assignment 2'\n","os.chdir(project_path)\n","print(os.listdir(project_path))"],"execution_count":3,"outputs":[{"output_type":"stream","name":"stdout","text":["Mounted at /content/drive\n","['assignment_2.pdf', 'data', 'p1t3.png', 'p2t3.png', 'KonstantinosGeorgiou_Assignment2_MATH525.ipynb']\n"]}]},{"cell_type":"code","metadata":{"id":"N7_bR5v6Tr-Q","executionInfo":{"status":"ok","timestamp":1634026867373,"user_tz":240,"elapsed":890,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}}},"source":["# Imports\n","import numpy as np\n","from scipy.io import loadmat\n","from scipy.stats import beta as beta_dist, norm\n","from scipy.optimize import fminbound, fmin\n","from scipy.special import erfc, beta as beta_func\n","import matplotlib as mpl\n","from matplotlib import pyplot as plt\n","from typing import *"],"execution_count":4,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"8OXR7C57TzfN"},"source":["## Load Data & Exploration"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"R-cgVX4TT4C8","executionInfo":{"status":"ok","timestamp":1634026868251,"user_tz":240,"elapsed":888,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"0bde8390-12c2-4c85-9145-9b96a8950bec"},"source":["# Load datasets\n","print(os.listdir('/content/drive/MyDrive/Documents/Work Documents/UTK/Courses/MATH525/Assignments/Assignment 2/data'))\n","geo_beta = loadmat('data/geo_beta.mat')\n","normal_normal = loadmat('data/normal_normal.mat')"],"execution_count":5,"outputs":[{"output_type":"stream","name":"stdout","text":["['geo_beta.mat', 'normal_normal.mat']\n"]}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"q7yiTtVoUMYr","executionInfo":{"status":"ok","timestamp":1634027231955,"user_tz":240,"elapsed":207,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"8785b8cc-57d7-4f57-b3c0-75b83001d003"},"source":["# Data exploration\n","# Geo Beta\n","print(\"--------\\nComponents of Geo Beta: \", geo_beta.keys())\n","geo_beta_w = np.fix(geo_beta['w']).astype(\"int64\").astype(\"float64\")\n","print(f\"Geo Beta_w W:\\n\\tType: {geo_beta_w.dtype}, Shape: {geo_beta_w.shape}\")\n","print(f\"\\tContents: {geo_beta_w[:, :].T}\")\n","# Normal Normal\n","print(\"--------\\nComponents of normal_normal: \", normal_normal.keys())\n","normal_normal_w = normal_normal['w']\n","print(f\"Normal Normal w:\\n\\tType: {normal_normal_w.dtype}, Shape: {normal_normal_w.shape}\")\n","print(f\"\\tContents: {normal_normal_w[:, :].T}\")"],"execution_count":8,"outputs":[{"output_type":"stream","name":"stdout","text":["--------\n","Components of Geo Beta:  dict_keys(['__header__', '__version__', '__globals__', 'w'])\n","Geo Beta_w W:\n","\tType: float64, Shape: (14, 1)\n","\tContents: [[ 6.  2. 67.  8. 14. 17. 12.  7.  6.  4.  6.  2. 11.  0.]]\n","--------\n","Components of normal_normal:  dict_keys(['__header__', '__version__', '__globals__', 'w'])\n","Normal Normal w:\n","\tType: float64, Shape: (25, 1)\n","\tContents: [[ 0.85098623  2.68116604  0.7415467   0.39038696  0.65444876  0.92733513\n","   0.94131924  1.67838023  1.30313855  2.08644262  0.64811303  2.30032071\n","  -0.00940472  2.37587415  1.25721046  1.66681344 -0.46541871  0.22992861\n","   2.67517127  3.52916018  1.22484276  2.10365845  3.28125189  3.27365833\n","  -0.36512257]]\n"]}]},{"cell_type":"markdown","metadata":{"id":"12d6DcnVVUda"},"source":["## Problem 1(iii) - compare graphically the prior and posterior of µ."]},{"cell_type":"code","metadata":{"id":"RKglXOclV1jC","executionInfo":{"status":"ok","timestamp":1634007801072,"user_tz":240,"elapsed":32,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}}},"source":["# from scipy.interpolate import make_interp_spline, BSpline\n","# def interpolate(x, y, k=3):\n","#   x, y = x[:,0], y[:,0]\n","#   xnew = np.linspace(x.min(), x.max(), 300)\n","#   spl = make_interp_spline(x, y, k=k)  # type: BSpline\n","#   ynew = spl(xnew)\n","#   return xnew, ynew"],"execution_count":151,"outputs":[]},{"cell_type":"code","metadata":{"id":"6Ve15n73VehZ","colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"status":"ok","timestamp":1634014608219,"user_tz":240,"elapsed":195,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"a8dc5e92-75c0-4d60-ca96-1eea1d657c9e"},"source":["def get_normal_pdf(x: np.ndarray, mean: float, variance: float) -> np.ndarray:\n","  pdf = 1/np.sqrt(2*np.pi*variance)*np.exp(-(x-mean)**2/(2*variance))\n","  return pdf\n","\n","def get_posterior_mean_and_var(data: np.ndarray, \n","                               mean_prior: float, variance_prior: float,\n","                               variance_likelihood: float) -> Tuple[np.ndarray, np.ndarray]:\n","  data_sum = np.sum(data)\n","  N = data.shape[0] \n","  mean_post = (variance_likelihood*mean_prior+variance_prior*data_sum)/(variance_likelihood+N*variance_prior)\n","  variance_post = (variance_likelihood*variance_likelihood)/(variance_likelihood+N*variance_prior)\n","  return mean_post, variance_post\n","\n","# Set the values we know\n","mean_prior = 0\n","variance_prior = 10\n","variance_likelihood = 1\n","print(f\"mean_prior={mean_prior}, variance_prior={variance_prior}\")\n","print(f\"variance_likelihood={variance_likelihood}\")\n","# Generate x values\n","x = np.linspace(-100, 100, 1000000)\n","# Calculate prior pdf\n","prior_pdf = get_normal_pdf(x=x, \n","                           mean=mean_prior, variance=variance_prior)\n","# Calculate posterior mean and variance\n","mean_post, variance_post = get_posterior_mean_and_var(data=normal_normal_w, \n","                                                      mean_prior=mean_prior,\n","                                                      variance_prior=variance_prior,\n","                                                      variance_likelihood=variance_likelihood)\n","print(f\"mean_post={mean_post}, variance_post={variance_post}\")\n","# Calculate posterior pdf\n","post_pdf = get_normal_pdf(x=x, \n","                          mean=mean_post, variance=variance_post)"],"execution_count":209,"outputs":[{"output_type":"stream","name":"stdout","text":["mean_prior=0, variance_prior=10\n","variance_likelihood=1\n","mean_post=1.433514252429064, variance_post=0.00398406374501992\n"]}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":459},"id":"GVQVY89xb7tF","executionInfo":{"status":"ok","timestamp":1634014610366,"user_tz":240,"elapsed":2149,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"763aba47-8505-4e3e-ab98-db409fd4b5bf"},"source":["print(\"Plot\")\n","# Create figure\n","fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n","# Plots\n","ax.plot(x, prior_pdf, label='Prior')\n","# ax.plot(*interpolate(data_sorted, prior_sorted), label='Prior Interpolated')\n","ax.plot(x, post_pdf, label='Posterior')\n","# ax.plot(*interpolate(data_sorted, post_sorted, k=2), label='Posterior Interpolated')\n","# Ranges\n","ax.set_xlim(-0.5, 3.5)\n","# Labels\n","ax.legend()\n","ax.set_title(\"Prior and Posterior of μ\")\n","ax.set_xlabel('μ')\n","ax.set_ylabel('Density')\n","fig.tight_layout()"],"execution_count":210,"outputs":[{"output_type":"stream","name":"stdout","text":["Plot\n"]},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5hddX3v8c937klmkkASEiBAELwQAgQM0UMMWK2CRSl6apWqFKWl1nqqtVKpT5+iPb2gh6o91iPi5RGr1niUoiLUy5GICGICDRggyEU04ZZJQiYzmcye2/f8sfbK7NnZM7P23uu31t6T9+t55tkze6/1W9+9GDKf9dvftZa5uwAAAABEWvIuAAAAAGgkBGQAAACgBAEZAAAAKEFABgAAAEoQkAEAAIASBGQAAACgBAEZwKxjZgNm9ry86yhnZivMzM2sLe9akjKzt5jZ9zPa1hwz+46Z9ZnZ/81imwBQCQEZQMMzsyfM7EAx+D5rZl80s+6plnf3bnd/PMsa01Dt+0ww3ofM7Mv11OTuX3H3V9czRhV+T9JSSYvc/Y0ZbRMADkFABtAsXufu3ZLOkrRG0t+UL1DvzGyDzOzO+D6zUs/+sEi1f2NOkPRLdx+tdbsAkAYCMoCm4u5PSrpV0ipJKrYs/JmZPSLpkZLnTi5+v8DMvmRmvWb2azP7mzi4mdllZvZTM/u4me2W9KHy7ZnZWjO7y8z2mtnTZvavZtZR8rqb2TvN7JHiMp8yMyu+1mpm15rZLjN7XNKFdbzPi8zsgeI2NprZKSU1fMDMnjSzfjN72MxeaWYXSPqgpDcVZ6TvK9kfny++lyfN7O/NrHWq/VF87o6SbZ1jZpuKbRCbzOycktc2mtk/mNlPJQ1KOqTNxcxOKS63t/h+Lio+/2FJf1tS7+UV1v2imf19yc8vN7MdSfcpACTVCLMlAJCYmR0n6Xck3Vjy9MWSXiLpQIVVPilpgaKwtkjS9yU9LenzxddfIulrij7ab6+w/pikv5C0WdJyRaH1XZI+UbLMayWdLWm+pHskfUfSf0r64+JrZ0raL+mbtbxPM3uBpH8vvs+NxXq+Y2YrJZ0o6d2Sznb3p8xshaRWd3/MzP5R0snu/taSob8oaaekkyXNk3SzpO2SPjPF/nhTSU1HSvqupD8v1vNGSd81s5PdfXdxsbdJeo2khyVZ2XtqL+6bL0h6taSXSfqWma1x96vNzCvUCwCZYwYZQLO4ycz2SrpD0o8l/WPJa//k7nvcfVJALs6MvlnSX7t7v7s/IemfFYW42FPu/kl3Hy1fX5Lc/R53/1nx9ScUBcnzyha7xt33uvtvJN0maXXx+d+X9Al33+7ueyT9U43v802SvuvuP3D3EUnXSpoj6RxFAb5T0koza3f3J9z9sUoDm9lSRaH7ve6+3913Svp4cR8l2R8XSnrE3f+t+Pq/S9om6XUly3zR3R8ovj5Stv5LJXUX99ewu/9IUUC/JMF+AYDMMIMMoFlc7O4/nOK17VM8v1jRLOivS577taRjE6wrSSrO3n5MUT/wXEX/bt5TttgzJd8PKgqBknRM2fildUzlkPdpZseUruvu42a2XdKx7r7RzN6rqD3kVDP7nqT3uftTFcY+QdH+eLrYBSJFEyWlNU63PybVUVTN/jxG0nZ3H59mfQDIHTPIAGYDn+L5XZJGFAXD2PGSnkywbuzTimZJn+/u8xX19dr0qxz0tKTjyrZdi6dU8h6KPc7Hqfg+3P2r7v6y4jIu6SPFRcvf23ZJBUmL3X1h8Wu+u59assx0+2NSHUXV7M+nJB1XdvJe+fozmVPy/YIq1gOAxAjIAGYtdx+T9HVJ/2BmPWZ2gqT3Sarm0mc9kvZJGjCzF0n60yrW/bqkPzez5WZ2hKSrqli3fJwLiyfftUv6S0VB904ze6GZvcLMOiUNKerDjmdon5W0Ig6k7v60oh7sfzaz+WbWYmYnmVl5y8hUbpH0AjP7AzNrM7M3SVqpqE0iibsVzbD/lZm1m9nLFbVnfC3h+pL0WjNbatHl794paa6ZdVWxPgDMiIAMYLb7H4pOkHtcUV/vVxWdJJbU+yX9gaR+SZ+VtKGKdT8r6XuS7pN0ryafWJiYuz8s6a2KTjjcpShUvs7dhxX1H19TfP4ZSUdJ+uviqvHNNnab2b3F7y+V1CHpQUnPSfqGpKMT1rFb0UmHfylpt6S/kvRad9+VcP3hYu2vKdb7fyRd6u7bkqxf9ICiExV/qeiEwT2SrqtifQCYkbnP9OkiAAD5M7MvStrh7rldGxrA4YEZZAAAAKAEARkAAAAoQYsFAAAAUIIZZAAAAKBEQ90oZPHixb5ixYq8ywAAAMAscM899+xy9yXVrtdQAXnFihXavHlz3mUAAABgFjCzJHcwPQQtFgAAAEAJAjIAAABQgoAMAAAAlGioHmQAAIDD1cjIiHbs2KGhoaG8S2k6XV1dWr58udrb21MZj4AMAADQAHbs2KGenh6tWLFCZpZ3OU3D3bV7927t2LFDJ554Yipj0mIBAADQAIaGhrRo0SLCcZXMTIsWLUp15p2ADAAA0CAIx7VJe78RkAEAAIASBGQAAABIklpbW7V69WqtWrVKb3zjGzU4OFhxuXPOOSfjyrJFQAYAAIAkac6cOdqyZYu2bt2qjo4OXXfddZNeHx0dlSTdeeediceM12kmBGQAAAAcYv369Xr00Ue1ceNGrV+/XhdddJFWrlwpSeru7pYUXUHiyiuv1KpVq3Taaadpw4YNklRxnWbCZd4AAAAazIe/84AefGpfqmOuPGa+rn7dqYmWHR0d1a233qoLLrhAknTvvfdq69ath1xG7cYbb9SWLVt03333adeuXTr77LN17rnnTrtOM2AGGQAAAJKkAwcOaPXq1VqzZo2OP/54XX755ZKktWvXVgy6d9xxhy655BK1trZq6dKlOu+887Rp06Zp12kGzCADAAA0mKQzvWmLe5DLzZs3r+qxalmnUTCDDAB5Gh+Xbr1K2vA2aagv72oAoCrr16/Xhg0bNDY2pt7eXt1+++1au3Zt3mXVjRlkAMjTw9+V7v509P2ik6XfvjrfegCgCq9//et111136YwzzpCZ6aMf/aiWLVumbdu25V1aXczd867hoDVr1vjmzZvzLgMAsvPVN0vP/EJaulJ69gHpvVulFj7cAw5HDz30kE455ZS8y2halfafmd3j7muqHYt/hQEgL2Oj0q9/Kj3/VdIpF0n7npR2PZx3VQBw2CMgA0BenrlfKuyTVrxMOv6l0XPbf55vTQAAAjIA5Obp4pniy8+O+o/nHCE9SZsZAOSNgAwAedm5TerolhYeL5lJR50q9dJiAQB5IyADQF56t0mLXxCFY0la8oLouQY6eRoADkcEZADIS+826aiSM66XvCi6FvLAs/nVBAAgIANALuIgvPgFE88tfn70uOuX+dQE4LDX2tqq1atXa9WqVXrjG9+owcHBqtZ/4okn9NWvfrWmbZ9zzjk1rRcCARkA8rD3N9HjESdMPLew+P3e7dnXAwCauNX01q1b1dHRoeuuu66q9WsJyKOjo5KkO++8s+p1QiEgA0Ae4hC84PiJ5xYsl2QT4RkAcrR+/Xo9+uij2rNnjy6++GKdfvrpeulLX6r7779fkvTjH/9Yq1ev1urVq3XmmWeqv79fV111lX7yk59o9erV+vjHP66xsTFdeeWVOvvss3X66afrM5/5jCRp48aNWr9+vS666CKtXLlSktTd3S1JcnddeeWVWrVqlU477TRt2LBhynVC4VbTAJCHvmJAXnjcxHNtnVLP0QRkANKtV0V32UzTstOk11yTaNHR0VHdeuutuuCCC3T11VfrzDPP1E033aQf/ehHuvTSS7VlyxZde+21+tSnPqV169ZpYGBAXV1duuaaa3Tttdfq5ptvliRdf/31WrBggTZt2qRCoaB169bp1a9+tSTp3nvv1datW3XiiSdO2vaNN96oLVu26L777tOuXbt09tln69xzz512nbQRkAEgD3t/I7V1SfOWTH5+4fEEZAC5OXDggFavXi0pmkG+/PLL9ZKXvETf/OY3JUmveMUrtHv3bu3bt0/r1q3T+973Pr3lLW/RG97wBi1fvvyQ8b7//e/r/vvv1ze+8Q1JUl9fnx555BF1dHRo7dq1FYPuHXfcoUsuuUStra1aunSpzjvvPG3atEnz58+fcp20BQ3IZrZQ0uckrZLkkt7h7neF3CYANIW+7VFLRXyJt9jC46XtP8unJgCNI+FMb9riHuQkrrrqKl144YW65ZZbtG7dOn3ve987ZBl31yc/+Umdf/75k57fuHGj5s2bV3V9taxTi9A9yP8i6T/d/UWSzpD0UODtAUBz2LtdWnDcoc8vPE7qe1IaH8u+JgCoYP369frKV74iKQq2ixcv1vz58/XYY4/ptNNO0wc+8AGdffbZ2rZtm3p6etTf339w3fPPP1+f/vSnNTIyIkn65S9/qf3798+4vQ0bNmhsbEy9vb26/fbbtXbt2nBvsIJgM8hmtkDSuZIukyR3H5Y0HGp7ANBU+rZLy1Yd+nzP0ZKPSft3ST1Ls68LAMp86EMf0jve8Q6dfvrpmjt3rm644QZJ0ic+8Qnddtttamlp0amnnqrXvOY1amlpUWtrq8444wxddtlles973qMnnnhCZ511ltxdS5Ys0U033TTt9l7/+tfrrrvu0hlnnCEz00c/+lEtW7ZM27Zty+LtSpLMA92xycxWS7pe0oOKZo/vkfQed5/ysGHNmjW+efPmIPUAQMMYG5H+5xLpvA9Iv/XXk1978NvS198m/cnt0tFn5FMfgFw89NBDOuWUU2ZeEBVV2n9mdo+7r6l2rJAtFm2SzpL0aXc/U9J+SVeVL2RmV5jZZjPb3NvbG7AcAGgQ+3dJcqn7qENf61kWPfZzNz0AyEvIgLxD0g53v7v48zcUBeZJ3P16d1/j7muWLFlS/jIAzD7xraS7K7RQxM8NPJNdPQCASYIFZHd/RtJ2M3th8alXKmq3AIDD28DO6LFSQGYGGTishWp9ne3S3m+hr4P8PyR9xcw6JD0u6e2BtwcAjW9/HJArfGrW1inNOYIZZOAw1NXVpd27d2vRokWy8ktAYkrurt27d6urqyu1MYMGZHffIqnqxmgAmNXiFot5FXqQJal7mdRPQAYON8uXL9eOHTvEOVnV6+rqqnijklpxJz0AyNrATqlzvtQxt/LrPUsnQjSAw0Z7e3smd4nDzELfKAQAUG7g2cpXsIh1L6MHGQByREAGgKwN7Jy6vUKS5i2WBndlVw8AYBICMgBkbWDn9DPIcxdJI4PS8GB2NQEADiIgA0DWBnZWvsRbbN7i6HFwdzb1AAAmISADQJZGC1KhT5o3zY2R5i6KHgnIAJALAjIAZGlwT/Q494iplzkYkOlDBoA8EJABIEsH4oC8aOpl5sYtFnvC1wMAOAQBGQCyFIfeOUdOvczc4mv7mUEGgDwQkAEgSwdnkKcJyF0LJWulBxkAckJABoAsJZlBbmmJAjQBGQByQUAGgCzFoXe6GWQp6lHmJD0AyAUBGQCydOA5qW2O1D5n+uXmLuYkPQDICQEZALI0uGfm2WMpWoaT9AAgFwRkAMjSgT3T9x/H5i2mxQIAckJABoAsDe6Z/iYhsTlHSAf2Su7hawIATEJABoAsHdgz/U1CYl0LJR+TCv3hawIATEJABoAsDSZssZhTnGUe2hu2HgDAIQjIAJCV8bEo8CY5SW/OwujxwHNhawIAHIKADABZGeqTfDzZDHJXHJCZQQaArBGQASArgwluMx2LZ5BpsQCAzBGQASArBxLcZjrGDDIA5IaADABZicNuPDs8nfgkPXqQASBzBGQAyMpQX/TYlSAgd8yTWtposQCAHBCQASArcdjtWjDzsmZRkKbFAgAyR0AGgKwcnEFOEJClqBWDGWQAyBwBGQCyMtQntXZK7V3Jlu9aSA8yAOSAgAwAWRnqSz57LEUn6tFiAQCZIyADQFaqDsi0WABAHgjIAJCVagMyJ+kBQC4IyACQlZpmkPuk8fFwNQEADkFABoCsDPUlu0lIrGuhJJcKfcFKAgAcioAMAFmp5SQ9iTYLAMgYARkAsuBeW4uFxIl6AJAxAjIAZGHkgDQ+Ul1A7pwfPQ7RYgEAWSIgA0AWqr2LniR1xQF5X/r1AACmREAGgCzUEpDjGeQCARkAskRABoAs1DSDXFyWGWQAyBQBGQCycDAgV3GZt86e6JEZZADIFAEZALJQywxya7vUPpcZZADIGAEZALIQX6qtmoAsRX3I3CgEADJFQAaALMQzyPGJd0l1zWcGGQAyRkAGgCwM9UltXVJ7V3Xrdc6nBxkAMkZABoAsVHsXvRgzyACQubaQg5vZE5L6JY1JGnX3NSG3BwANq9aA3Dlf2vub9OsBAEwpaEAu+i1335XBdgCgcTGDDABNgxYLAMhCPTPI9CADQKZCB2SX9H0zu8fMrqi0gJldYWabzWxzb29v4HIAICc1zyAvkEaHpNHh9GsCAFQUOiC/zN3PkvQaSX9mZueWL+Du17v7Gndfs2TJksDlAEBOCvsm7oxXjfiycMwiA0BmggZkd3+y+LhT0n9IWhtyewDQsAr91V8DWYp6kKWJ6ygDAIILFpDNbJ6Z9cTfS3q1pK2htgcADWt0OGqTYAYZAJpCyKtYLJX0H2YWb+er7v6fAbcHAI1peCB6rCUgH5xBJiADQFaCBWR3f1zSGaHGB4CmUeiPHplBBoCmwGXeACC0egIyM8gAkDkCMgCEVtcMcvHScMwgA0BmCMgAENrBgFzPVSwIyACQFQIyAIQWz/7WMoPc2i61zWEGGQAyREAGgNDiGeSO7trW75rPdZABIEMEZAAIrZ7LvElRawYzyACQGQIyAISWygwyARkAskJABoDQCv1SR4/UUuM/ucwgA0CmCMgAEFphX+3tFRIzyACQMQIyAIRW6K8vIHf2TPQxAwCCIyADQGiFAamzxv5jqdhiQUAGgKwQkAEgtFRmkPul8fH0agIATImADACh1RuQ46tf0GYBAJkgIANAaIX+2m4zHYvDNQEZADJBQAaA0NJosYjHAQAER0AGgJDco/5hAjIANA0CMgCENDIo+TgBGQCaCAEZAEKq9zbTpesSkAEgEwRkAAgpDrWcpAcATYOADAAhFYq3iK6rxaIYrplBBoBMEJABIKSDM8j1BOS4xWJf/fUAAGZEQAaAkNIIyG2dUmsHt5sGgIwQkAEgpDjU1hOQpehEPVosACATBGQACCmNGeR4fQIyAGSCgAwAIaVxkp4UnajHVSwAIBMEZAAIqdAf9Q+3ddY3TictFgCQFQIyAIRUqPM20zFaLAAgMwRkAAiJgAwATYeADAAhDQ+kE5C5igUAZIaADAAhFfrru810rLOHk/QAICMEZAAIqbAvmv2tV2ePNDIojY3WPxYAYFoEZAAIKc0eZEkaps0CAEIjIANASGkHZG43DQDBEZABIKS0AnLcpsGJegAQHAEZAEIZG5FGh1I6Sa84BifqAUBwBGQACCWe7U21xWJf/WMBAKZFQAaAUA4G5DSuYkGLBQBkhYAMAKEEmUGmxQIAQiMgA0AoaQZkTtIDgMwQkAEglIMBOaU76UmcpAcAGSAgA0AowynOILe2S21zOEkPADJAQAaAUNJssZCiE/VosQCA4AjIABBK6gG5h5P0ACADwQOymbWa2X+Z2c2htwUADSUOyO3z0hmvgxlkAMhCFjPI75H0UAbbAYDGUuiXOnqklpT+qe2cT0AGgAwEDchmtlzShZI+F3I7ANCQCvvSa6+QorGGCcgAEFroGeRPSPorSeNTLWBmV5jZZjPb3NvbG7gcAMhQYSDlgEyLBQBkIVhANrPXStrp7vdMt5y7X+/ua9x9zZIlS0KVAwDZK/SnP4PMSXoAEFzIGeR1ki4ysyckfU3SK8zsywG3BwCNJUhAZgYZAEILFpDd/a/dfbm7r5D0Zkk/cve3htoeADScQn/UFpGWjh5prCCNDqc3JgDgEFwHGQBCKfSnc5vpGLebBoBMZBKQ3X2ju782i20BQMMYTrvFojgbze2mASAoZpABIAT3MD3IEn3IABAYARkAQhgZlHw8uvtdWg4GZFosACAkAjIAhBDP8nal2IPcwQwyAGSBgAwAIcQhtiNAiwV30wOAoAjIABBCHJBD9CAPcZIeAIREQAaAEA4G5DR7kItjcZk3AAiKgAwAIcQhNs0Z5PiEP3qQASAoAjIAhBCixaKlNQrJXMUCAIIiIANACCFO0pOiwM2NQgAgKAIyAIQQYgY5Ho8WCwAIioAMACEU+qWWdqmtM91xCcgAEBwBGQBCKPRHV50wS3fcjm4CMgAERkAGgBCGB9Jvr5CiMbnMGwAERUAGgBAK/VJnireZjnXOZwYZAAIjIANACIV9E9ctThNXsQCA4AjIABBCIWCLRaFfck9/bACAJAIyAIQRn6SXts4eycelkcH0xwYASCIgA0AYhf5AM8jcbhoAQiMgA0AIwwPhTtKTuN00AAREQAaAtI2PRwE51El6EifqAUBABGQASFt8neJQJ+lJtFgAQEAEZABIWxxeCcgA0JQSBWQzu9HMLjQzAjUAzORgQA7ZYkFABoBQkgbe/yPpDyQ9YmbXmNkLA9YEAM3tYItFgJP0OnombwMAkLpEAdndf+jub5F0lqQnJP3QzO40s7ebWXvIAgGg6cQn0AVtseAkPQAIJXHLhJktknSZpD+S9F+S/kVRYP5BkMoAoFnF7Q8hrmLR1im1tNNiAQABtSVZyMz+Q9ILJf2bpNe5+9PFlzaY2eZQxQFAUyoEvIqF2cTtpgEAQSQKyJI+6+63lD5hZp3uXnD3NQHqAoDmFfIqFvG4BGQACCZpi8XfV3jurjQLAYBZI2SLhRSd/EdABoBgpp1BNrNlko6VNMfMzpRkxZfmS5obuDYAaE7D/VJbl9TWEWb8zm4CMgAENFOLxfmKTsxbLuljJc/3S/pgoJoAoLkV+sPNHktRi8XAznDjA8BhbtqA7O43SLrBzP67u38zo5oAoLkVBsL1H0vR2LsfCzc+ABzmZmqxeKu7f1nSCjN7X/nr7v6xCqsBwOGt0B8+INNiAQDBzNRiMa/4GPCzQgCYZQjIANDUZmqx+Ezx8cPZlAMAs8Bwv9RzdLjxO+dLoweksVGpNenVOgEASSW6zJuZfdTM5ptZu5n9PzPrNbO3hi4OAJpS6Bnk+ATAYWaRASCEpNdBfrW775P0WklPSDpZ0pWhigKAppbFVSzi7QAAUpc0IMef4V0o6f+6e1+gegCg+WVxFQuJgAwAgSRtXrvZzLZJOiDpT81siaShcGUBQJMaG4n6gwnIANC0Es0gu/tVks6RtMbdRyTtl/S7IQsDgKYUh9agAXn+5G0BAFJVzenPL1J0PeTSdb6Ucj0A0NyGB6JHZpABoGklCshm9m+STpK0RdJY8WkXARkAJotDa9CT9LonbwsAkKqkM8hrJK10dw9ZDAA0vQIzyADQ7JJexWKrpGXVDGxmXWb2czO7z8weMDNuNgJg9jvYgzw/3DY6mEEGgJCSziAvlvSgmf1cUiF+0t0vmmadgqRXuPuAmbVLusPMbnX3n9VeLgA0uMK+6LEzYItFS2sUkgnIABBE0oD8oWoHLrZjFD9rVHvxixYNALNbFifpxeNzJz0ACCLpZd5+rOgOeu3F7zdJunem9cys1cy2SNop6QfufneFZa4ws81mtrm3t7eq4gGg4WRxmbd4fGaQASCIRAHZzP5Y0jckfab41LGSbpppPXcfc/fVkpZLWmtmqyosc727r3H3NUuWLEleOQA0oiyuYhGPT0AGgCCSnqT3Z5LWSdonSe7+iKSjkm7E3fdKuk3SBdUWCABNpdAvtc+L+oRDYgYZAIJJGpAL7j4c/1C8Wci0/cRmtsTMFha/nyPpVZK21VooADSFQn/YE/RiBGQACCbpSXo/NrMPSppjZq+S9C5J35lhnaMl3WBmrYqC+Nfd/ebaSwWAJlDoD99/LEWXkSMgA0AQSQPyVZIul/QLSX8i6RZJn5tuBXe/X9KZdVUHAM1meCCjgMwMMgCEkiggu/u4md0k6SZ351ITADCVQn/4E/SkiYDsLpmF3x4AHEam7UG2yIfMbJekhyU9bGa9Zva32ZQHAE1maF/Yu+jFOrslH5NGDoTfFgAcZmY6Se8vFF294mx3P9Ldj5T0EknrzOwvglcHAM2m0C91ZRGQeya2BwBI1UwB+W2SLnH3X8VPuPvjkt4q6dKQhQFAUyr0ZTSDXNwGARkAUjdTQG53913lTxb7kNvDlAQATco9+xlkbjcNAKmbKSAP1/gaABx+hgckH89oBpkWCwAIZaarWJxhZvsqPG+SugLUAwDNa6j4zyU9yADQ1KYNyO4e+F6pADCLFIoBOYsZ5PhScgRkAEhd0ltNAwBmkukMMifpAUAoBGQASEuWM8gHWywqdcEBAOpBQAaAtAz1RY9ZBOS2TqmlXSoMhN8WABxmCMgAkJZChi0WZhO3mwYApIqADABpGcqwxUIiIANAIARkAEhLYZ9krVLHvGy2R0AGgCAIyACQlqF9UWg1y2Z7nT2cpAcAARCQASAthX3Z9B/HuhZMnBgIAEgNARkA0lLolzoXZLe9rgXMIANAAARkAEjLUMYzyJ3zmUEGgAAIyACQlkJfdlewkIotFvsk9+y2CQCHAQIyAKQlPkkvK10LJB+Thvdnt00AOAwQkAEgLZmfpFfcFm0WAJAqAjIApMG9OIOccYuFREAGgJQRkAEgDSODUbtD1ifpSVzJAgBSRkAGgDRkfZtpSepaWNw2M8gAkCYCMgCkIZ7F7cr4OsgSARkAUkZABoA05DKDzEl6ABACARkA0lAohtQ8epAJyACQKgIyAKQhjxnk9i6prYuADAApIyADQBriHuQsbxQiRYGcq1gAQKoIyACQhngGOcsWC6l4u2lmkAEgTQRkAEhDoV+SSR0ZzyB3LZgI5wCAVBCQASANhX1Re0VLxv+sds1nBhkAUkZABoA0ZH2b6RgtFgCQOgIyAKRhqC/bm4TEuhZwkh4ApIyADABpGNorzVmY/XY7abEAgLQRkAEgDQf2Sl05BOSuBdLokDRayAV37Y4AABcFSURBVH7bADBLEZABIA15zSDHbR1cyQIAUkNABoA0HNibXw+yRJsFAKSIgAwA9RobkUb259diIRGQASBFBGQAqFccTvNssSgQkAEgLQRkAKjXgb3RYx4zyPG1l5lBBoDUEJABoF5DxYCc60l6BGQASAsBGQDqFQfkXE7Si2eQuYoFAKSFgAwA9cqzxaKjW7IWZpABIEXBArKZHWdmt5nZg2b2gJm9J9S2ACBXebZYmEUz1wRkAEhNW8CxRyX9pbvfa2Y9ku4xsx+4+4MBtwkA2ctzBlkqBuS9+WwbAGahYDPI7v60u99b/L5f0kOSjg21PQDIzVCf1NYltXfls/05R0yEdABA3TLpQTazFZLOlHR3hdeuMLPNZra5t7c3i3IAIF1DOd1FLzbnCOnAc/ltHwBmmeAB2cy6JX1T0nvd/ZDTrN39endf4+5rlixZErocAEjfgb35tVdIBGQASFnQgGxm7YrC8Vfc/caQ2wKA3AztzecEvRgBGQBSFfIqFibp85IecvePhdoOAOSuEWaQh/ZK4+P51QAAs0jIGeR1kt4m6RVmtqX49TsBtwcA+Rjqy38G2celAjcLAYA0BLvMm7vfIclCjQ8ADaMRTtKTojaLPIM6AMwS3EkPAOoxPh7d5jnvFguJPmQASAkBGQDqUeiT5Pm3WEgEZABICQEZAOoR3+KZGWQAmDUIyABQj/gOdswgA8CsQUAGgHrEoTTPGeR429xuGgBSQUAGgHoc2BM9zj0yvxraOqSObmaQASAlBGQAqMdgHJAX5VsHd9MDgNQQkAGgHnFAjvuA8zJnIQEZAFJCQAaAehzYI3XOl1rb862DGWQASA0BGQDqMbgn/9ljiYAMACkiIANAPQ7syfcEvRgBGQBSQ0AGgHoM7s7/BD1pIiC7510JADQ9AjIA1GNwjzSnQWaQx0ek4f15VwIATY+ADAD1OPBc47RYSLRZAEAKCMgAUKuxEamwr3FmkCUCMgCkgIAMALUabIC76MUOBuQ9+dYBALMAARkAatUIt5mOxScKDu7Otw4AmAUIyABQq4N30WuEgLw4etxPQAaAehGQAaBWDTWDfKQkkwZ35V0JADQ9AjIA1KqRZpBbWqM+5P0EZACoFwEZAGoV9/s2wo1CJGneYmaQASAFBGQAqNWBPVJbl9QxN+9KInMX04MMACkgIANArQafa4z2iti8RcwgA0AKCMgAUKvBXY3TXiEVZ5AJyABQLwIyANRqf6/UvSTvKibMWxy1fYyP510JADQ1AjIA1Gp/rzSvgQLy3MWSj3O7aQCoEwEZAGq1f1djBeR5xZuF0IcMAHUhIANALQoD0shgYwXkuB+aPmQAqAsBGQBqsb83emykgMwMMgCkgoAMALWIZ2m7j8q3jlJxWGcGGQDqQkAGgFrs3xk9xrO2jSBusRjkZiEAUA8CMgDUohFbLFrbpa4FzCADQJ0IyABQi0YMyFJ0qTd6kAGgLgRkAKjFQK/UuUBq68y7ksnmLZ4I7wCAmhCQAaAW+3sbq/841n2UNLAz7yoAoKkRkAGgFvt7G+sKFrHuZVL/M3lXAQBNjYAMALXYv6sxZ5B7lkpDe6WRobwrAYCmRUAGgFrs39l4J+hJ0QyyJA08m28dANDECMgAUK2xEWlwjzSvAVssegjIAFAvAjIAVGtgpySfCKONJO6Lpg8ZAGpGQAaAasXhc/4x+dZRCS0WAFA3AjIAVKv/qeixEWeQ5y2WrIUZZACoAwEZAKoVh8+eo/Oto5KW1qg3eoCADAC1ChaQzewLZrbTzLaG2gYA5KL/aclao9s6N6KepVI/LRYAUKuQM8hflHRBwPEBIB/9z0TtFS0N+iFc9zJmkAGgDsH+dXf32yXtCTU+AOSm/+nGbK+IMYMMAHXJffrDzK4ws81mtrm3tzfvcgBgZvEMcqPqXhbdCntsNO9KAKAp5R6Q3f16d1/j7muWLGnAu1IBQLl9TzX+DLI8utsfAKBquQdkAGgqIwekob2NPYM8f3n0uO+pfOsAgCZFQAaAajTyTUJiC4oBuW97vnUAQJMKeZm3f5d0l6QXmtkOM7s81LYAIDP9T0eP3UvzrWM6BwPyjnzrAIAm1RZqYHe/JNTYAJCbOHQuOC7fOqbTtUDq6CEgA0CNaLEAgGrs/XX0uLCBA7JZNItMQAaAmhCQAaAae7dL85ZI7XPyrmR6C5bTgwwANSIgA0A1+rY3dntFbMFyqe/JvKsAgKZEQAaAauzd3tjtFbEFx0qDu6LL0gEAqkJABoCk3KO+3qaYQS7WyCwyAFSNgAwASe3fJY0ekBYen3clM+NayABQMwIyACTV95vosSlmkAnIAFArAjIAJLW3GDaboQd5/nKppU3a86u8KwGApkNABoCk9jbRDHJrm7TwBGnP43lXAgBNh4AMAEnteUyau0iaszDvSpI58nkEZACoAQEZAJLa/Zh05El5V5Hckc+LWizc864EAJoKARkAktr9mLSoyQLycH909Q0AQGIEZABIYni/1P9U880gS7RZAECVCMgAkEQcMhc9L986qkFABoCaEJABIIndj0WPi07Ot45qLDxeshYCMgBUiYAMAEnsKQbkI5toBrmtIwrJux/JuxIAaCoEZABIYtejUvdSqbMn70qqs+QUaee2vKsAgKZCQAaAJHY+IB11St5VVO+oU6IZ5NHhvCsBgKZBQAaAmYyPSb0PS0edmncl1TtqpTQ+OtEiAgCYEQEZAGay51fS6JC0dGXelVQvnvXe+WC+dQBAEyEgA8BMdj4QPR7VhAF58fMla5V2PpR3JQDQNAjIADCTZx+UZNKSF+VdSfXaOqO7/xGQASAxAjIAzOTZrdHl3Trm5l1JbZaeKj1zf95VAEDTICADwEye+i/pmDPzrqJ2x5wl7f2NtH9X3pUAQFMgIAPAdPY9Le17Ujr2xXlXUru49ifvybcOAGgSBGQAmM5T90aPzRyQjz4juuU0ARkAEiEgA8B0dmyWWtqko0/Pu5LadXZHd9QjIANAIgRkAJjOjk3RSW7tc/KupD7HnhUF5PHxvCsBgIZHQAaAqYwMRQH5hHV5V1K/E86RDjw3cU1nAMCUCMgAMJUdP4/uoHfiuXlXUr8Tz4seH/9xvnUAQBMgIAPAVH51e3QXutkwg7zgWGnRydKvCMgAMBMCMgBM5fGN0fWPu+bnXUk6nvdy6YmfSqPDeVcCAA2NgAwAlfQ/E13B4gXn511Jek56pTSyX3riJ3lXAgANjYAMAJVs+64kl055Xd6VpOek35I6uqUHv5V3JQDQ0AjIAFDJQ9+JenaXvCjvStLTPieaEd92szQ2mnc1ANCwCMgAUG7v9qj/+NTXS2Z5V5OulRdLg7ulx36UdyUA0LAIyABQ7r++HD2e+bZ86wjhBRdI846SNn8+70oAoGERkAGg1MiQdM8XpZNeIR1xQt7VpK+tQ3rxH0q//J6051d5VwMADakt7wKy5O5lP1dYJsF6Uy9XabyZt1lJkrGmXq58mWT1l4s/WLaSj5hLP2yOn7aSZyt9Gl36XLzs5OdKl7WptzPbPupGY7r3S9LAM9Ibrs+7knDWvEP66b9It18rXfypvKsBgIbTUAH5F0/26aQP3nLw53qCKQ4fE0G99LlDQ32loK4pg3rZclVsZ+oxpz44KB8n+tFkNvFz6cGFlY0Xv176WnGIST+Xvq/SAw+bYhyVvla6zdKaymosH+eQbejQ5yq919Jxyt9H+X6YeK9Tj1O6DZW/Xhync2xA737gGu3pXq0vPbBY9uCDxfVtYlsmtRTXbbGJfdRS+prZwfFbyverFZ+T1NJSPm7pf4vyca2sltLlJ6/f0jJ5e6XjRLV0asUL3qpl931RD6y4VIUjnj/pPbbY5H3XkrA+lbyng9urVF/xhZYp67ODr8WPAJClhgrIS3o69afnnTTpuYozkpVWrrBgpeXKF7MKSyXdZuXZ0pn/Ia88fsp1lC1Z69+X0oOP0hns+HmvsGyl5SaPWfl1r/icT7lc6QuV6php/anemyq+t+nrqOm9T1GHe/Ra9HzxZ5/8Xlwlzx382Uuen/hZPjGuStaruJ0pxil/j6U1HVxuXHKNTzmOfGLvTFq3wna8Qj3x+9QhY0+MM+VrJe9z8j4sWdZdf6PrNdf36I9G3q+tm7ZPMU70/XjJus3oCK3VbZ0bNHzju/T7w1drTK15lzSt0rBdMcCrNHhPhPPyA4vyA45Kjy2lgd9MLS3x6yUHN2aT1yk5KCkdYyLkV1in5GCqfJ34YKt8nYMHVpXWUclzLYcelExap+TgbPI4E/t48jolB3tVrBO/j/J1yvdt6TKT6m0pGVeTt9NiJmuZ4r8HB1ioU0MF5GXzu/T+81+YdxkADkf3f1268YfSOe/WhvPfVdWq5aF53CcfRIx7SRgfn/zceFnInwjehz43Me5EOB8v2bY0+edJ40gaH58YRy4986t/0Fk/fa82nnm7Hlv9V5MOBkrr08F6K43rGh/XpHEPqbd8H1Wo95D3F9dQ3LbK9tXk/Ta5vsnjTrxWWsPksQ59HPfJ68U/Txp3PPp5bNw1MlayjkrWGT+0lvLtlP53ntju5For1Vi+DqZWfrAx3QFLix0axCWVhPlpDj5apgnzMxyQTVq+pcKBSPmYLVNvo+KYpQchUxy4HHIQNcX65csfug+nO8iaGL/ag6cZxyw5cIpfq1VDBWQAyMVD35G+9WfSCeukV15d9eoHZ/wqf77VuE5+u1T4hY7b/Fkdt2ypdO77a/+4CbmbNlSr+DheFsRVts54pTA/8f34eIV1isvEBy6lwb9yoJ98MFTVOsX3OT5e4YCldJ3xifc3cbA39bbLD1gqraMK9ZfXeMg6mthv8cFU/D7Gxn3Kg6ny/TH1Pqv83soPtkpfQzJBA7KZXSDpXyS1Svqcu18TcnsAUJVCv/Tjj0h3/qu0fI30pi9HV3k4nPzO/5IK+6Tb/l565n7pNR+R5h+Td1WoQdMeqCEz032KcjBElx9ETXEQU34ANX2wnzgASzZm8u1OqnP80AOFd36ktn0VLCCbWaukT0l6laQdkjaZ2bfd/cFQ2wSAGQ3vl57aIj18i7TlK9KB56QXv106/x+kjnl5V5e9llbpDZ+VjlopbbxGeuT70qlvkFZeJJ1wjtS1IO8KAaTkcDyIemeN64WcQV4r6VF3f1ySzOxrkn5X0tQBeeBZ6Sf/PP2oiT8fSLhc4o8bko6Xdn2MF368Rq6tCrPiv0WA8cYK0uCe6Kv/aem5J6JttHZIz3+19LK/iGaPD2dm0vr3SadeLN35yagf+76vRq/1HC0dsUKauygKy509UkvboV/1tGbU1dZx+PyhB5CdkAH5WEnbS37eIekl5QuZ2RWSrpCkFx/dIv2/vwtYEiIJ/6Ak/qN1OI3XyLXlOV7C4fKor7U9CndzjpSOWS2dcYm0bJW0Yr3UNT/hdg4TRz5Peu3HpfP/SdqxSdp+t7Tncem5X0cHFkN90tA+aXxEGh+TxkclH8u7agBIXe4n6bn79ZKul6Q1L36x62/uTLDWbAkVaW8XAFLQ3iWduD76mol7MSjX+olIHZ+kcMYRgJl8eE5Nq4UMyE9KOq7k5+XF56ZmJrV1BiwJAJAqs2iWHgBmkZaAY2+S9HwzO9HMOiS9WdK3A24PAAAAqFuwGWR3HzWzd0v6nqLLvH3B3R8ItT0AAAAgDUF7kN39Fkm3hNwGAAAAkKaQLRYAAABA0yEgAwAAACUIyAAAAEAJAjIAAABQgoAMAAAAlCAgAwAAACUIyAAAAEAJAjIAAABQgoAMAAAAlCAgAwAAACUIyAAAAEAJc/e8azjIzPolPZx3HbPAYkm78i5ilmBfpod9mR72ZXrYl+lhX6aHfZmeF7p7T7UrtYWopA4Pu/uavItodma2mf2YDvZletiX6WFfpod9mR72ZXrYl+kxs821rEeLBQAAAFCCgAwAAACUaLSAfH3eBcwS7Mf0sC/Tw75MD/syPezL9LAv08O+TE9N+7KhTtIDAAAA8tZoM8gAAABArgjIAAAAQIlcA7KZHWlmPzCzR4qPR0yx3JiZbSl+fTvrOhuVmV1gZg+b2aNmdlWF1zvNbEPx9bvNbEX2VTaHBPvyMjPrLfk9/KM86mx0ZvYFM9tpZluneN3M7H8X9/P9ZnZW1jU2iwT78uVm1lfyO/m3WdfYLMzsODO7zcweNLMHzOw9FZbhdzOBhPuS380EzKzLzH5uZvcV9+WHKyzD3/EZJNyP1f8Nd/fcviR9VNJVxe+vkvSRKZYbyLPORvyS1CrpMUnPk9Qh6T5JK8uWeZek64rfv1nShrzrbsSvhPvyMkn/mnetjf4l6VxJZ0naOsXrvyPpVkkm6aWS7s675kb9SrAvXy7p5rzrbIYvSUdLOqv4fY+kX1b4f5zfzfT2Jb+byfalSeouft8u6W5JLy1bhr/j6ezHqv+G591i8buSbih+f4Oki3OspdmslfSouz/u7sOSvqZof5Yq3b/fkPRKM7MMa2wWSfYlEnD32yXtmWaR35X0JY/8TNJCMzs6m+qaS4J9iYTc/Wl3v7f4fb+khyQdW7YYv5sJJNyXSKD4uzZQ/LG9+FV+5QT+js8g4X6sWt4Beam7P138/hlJS6dYrsvMNpvZz8yMEB05VtL2kp936NB/pA4u4+6jkvokLcqkuuaSZF9K0n8vfvT6DTM7LpvSZp2k+xrJ/Lfix4q3mtmpeRfTDIofUZ+paJapFL+bVZpmX0r8biZiZq1mtkXSTkk/cPcpfy/5Oz61BPtRqvJvePCAbGY/NLOtFb4mzdB5NAc+VeI/waNbLv6BpE+Y2Umh6wbKfEfSCnc/XdIPNHFED+TlXkX/Np4h6ZOSbsq5noZnZt2Svinpve6+L+96mtkM+5LfzYTcfczdV0taLmmtma3Ku6ZmlGA/Vv03PHhAdvffdvdVFb6+JenZ+COs4uPOKcZ4svj4uKSNio5YD3dPSio9AlpefK7iMmbWJmmBpN2ZVNdcZtyX7r7b3QvFHz8n6cUZ1TbbJPm9RQLuvi/+WNHdb5HUbmaLcy6rYZlZu6JA9xV3v7HCIvxuJjTTvuR3s3ruvlfSbZIuKHuJv+NVmGo/1vI3PO8Wi29L+sPi938o6VvlC5jZEWbWWfx+saR1kh7MrMLGtUnS883sRDPrUNS8X36Fj9L9+3uSflScqcdkM+7Lsl7EixT13aF635Z0afGKAS+V1FfSZoUqmNmyuBfRzNYq+vecP5wVFPfT5yU95O4fm2IxfjcTSLIv+d1MxsyWmNnC4vdzJL1K0rayxfg7PoMk+7GWv+FtaRZZg2skfd3MLpf0a0m/L0lmtkbSO939jySdIukzZjau6H+ya9z9sA/I7j5qZu+W9D1FV2H4grs/YGZ/J2mzu39b0T9i/2Zmjyo62efN+VXcuBLuyz83s4skjSral5flVnADM7N/V3QG+2Iz2yHpakUnTMjdr5N0i6KrBTwqaVDS2/OptPEl2Je/J+lPzWxU0gFJb+YP55TWSXqbpF8U+xQl6YOSjpf43axSkn3J72YyR0u6wcxaFeWbr7v7zfwdr1qS/Vj133BuNQ0AAACUyLvFAgAAAGgoBGQAAACgBAEZAAAAKEFABgAAAEoQkAEAAIASBGQAAACgBAEZAAAAKEFABoAmYWYrzGxryc/vN7MP5VgSAMxKBGQAAACgBAEZAAAAKEFABoDmYiXft+dWBQDMYgRkAGguJ5jZEjNrkXSupNa8CwKA2YaADADNZbekL0naLGmrpEvN7KR8SwKA2cXcPe8aAAAJmNkKSTe7+6qcSwGAWY0ZZAAAAKAEM8gAAABACWaQAQAAgBIEZAAAAKAEARkAAAAoQUAGAAAAShCQAQAAgBL/H5WayFbJNv5HAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 720x432 with 1 Axes>"]},"metadata":{"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"Y6h1oN8zubA3"},"source":["## Problem 1(iv) -  compute the prior and posterior probability that µ takes a negative value."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"r_5zBE1YucoU","executionInfo":{"status":"ok","timestamp":1634010252020,"user_tz":240,"elapsed":90,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"1155d9f8-2c6c-4f4c-a9d2-436c99dafeb4"},"source":["from scipy.integrate import quad\n","print(f\"mean_post={mean_post}, variance_post={variance_post}\")\n","# Set the values we know\n","mean_prior = 0\n","variance_prior = 10\n","variance_likelihood = 1\n","\n","# Calculate prior probability using integral\n","prior_I = quad(get_normal_pdf, -np.inf, 0, args=(mean_prior, variance_prior))\n","print(\"Integral for prior: \", prior_I[0])\n","\n","# Calculate posterior probability using integral\n","post_I = quad(get_normal_pdf, -np.inf, 0, args=(mean_post, variance_post))\n","print(\"Integral for posterio: \", post_I[0])"],"execution_count":183,"outputs":[{"output_type":"stream","name":"stdout","text":["mean_post=1.433514252429064, variance_post=0.00398406374501992\n","0.5\n","1.766059954858669e-114\n"]}]},{"cell_type":"markdown","metadata":{"id":"qyOeF1M6re5L"},"source":["## Problem 2(iii) - compare graphically the prior and posterior of µ."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"2fo84bF2re5M","executionInfo":{"status":"ok","timestamp":1634014639896,"user_tz":240,"elapsed":82,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"bfdf9093-f1d6-41ab-b923-ac81027aebda"},"source":["def get_beta_pdf(x: np.ndarray, alpha: float, beta: float, show: bool=False) -> np.ndarray:\n","  beta_func_out = beta_func(alpha, beta)\n","  left = x**(alpha-1)\n","  right = ((1-x)**(beta-1))/beta_func_out\n","  pdf = left * right\n","  if show:\n","    print(f\"alpha={alpha}, beta={beta}\")\n","    print(\"beta_func: \", beta_func_out)\n","    print(\"x:\\n\", x.T)\n","    print(\"left:\\n\", left.T)\n","    print(\"right:\\n\", right.T)\n","    print(\"pdf:\\n\", pdf.T)\n","  return pdf\n","\n","def get_posterior_alpha_and_beta(data: np.ndarray, \n","                                 alpha_prior: float,\n","                                 beta_prior: float) -> Tuple[np.ndarray, np.ndarray]:\n","  data_sum = np.sum(data)\n","  N = data.shape[0] \n","  alpha_post = N+alpha_prior\n","  beta_post = data_sum+beta_prior\n","  return alpha_post, beta_post\n","\n","# Set the values we know\n","alpha_prior = 1\n","beta_prior = 1\n","# Generate x values\n","x = np.linspace(0, 1, 100000)\n","print(f\"alpha_prior={alpha_prior}, beta_prior={beta_prior}\")\n","# Calculate prior pdf\n","prior_pdf = get_beta_pdf(x=x, \n","                         alpha=alpha_prior, beta=beta_prior)\n","# Calculate posterior alpha and beta\n","alpha_post, beta_post = get_posterior_alpha_and_beta(data=geo_beta_w, \n","                                                     alpha_prior=alpha_prior,\n","                                                     beta_prior=beta_prior)\n","print(f\"alpha_post={alpha_post}, beta_post={beta_post}\")\n","# Calculate posterior pdf\n","post_pdf = get_beta_pdf(x=x, \n","                        alpha=alpha_post, beta=beta_post, show=0)"],"execution_count":212,"outputs":[{"output_type":"stream","name":"stdout","text":["alpha_prior=1, beta_prior=1\n","alpha_post=15, beta_post=163.0\n"]}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":477},"id":"fpuTEZvWre5N","executionInfo":{"status":"ok","timestamp":1634014640638,"user_tz":240,"elapsed":637,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"02d72a5a-8cb2-4619-a9c4-80a918d4a96a"},"source":["print(\"Plot\")\n","# Create figure\n","fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n","# Plots\n","ax.plot(x, prior_pdf, label='Prior')\n","ax.plot(x, post_pdf, label='Posterior')\n","# Labels\n","ax.legend()\n","ax.set_title(\"Prior and Posterior of π\")\n","ax.set_xlabel('π')\n","ax.set_ylabel('Density')\n","fig.tight_layout()\n","print(np.sum(post_sorted, axis=0))"],"execution_count":213,"outputs":[{"output_type":"stream","name":"stdout","text":["Plot\n","[inf]\n"]},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de7xddX3n/9cnJ/f7lQAJEFS0cg0YQEW8Trl4Qe1oK623lpZpf3am1o5T2umj2namtY7VjpcR6cioHbXYFhlaoWhVRCpaAg0aBAQxSMIt95xDkpOcnM/vj7VP2BxOkn2Ss/faa+3X8/E4j7XPuuz1OWc/Au9881nfb2QmkiRJkgqTyi5AkiRJ6iYGZEmSJKmJAVmSJElqYkCWJEmSmhiQJUmSpCYGZEmSJKmJAVlSz4uIgYh4Vtl1jBYRKyIiI2Jy2bW0KiJ+KSK+2qF7zYiIf4iI7RHxt524p6TeYECWVDsRsS4idjWC7+MR8ZmImH2g8zNzdmY+2MkaJ8J4f84W3u/9EfF/j6SmzPx8Zl5wJO8xDm8ClgKLMvPNHbqnpB5gQJZUV6/LzNnAWcAq4A9Gn3CkI7NdMrJ7yJ+zU47k9xGF8f4/6QTgR5k5dLj3laSxGJAl1VpmbgBuBE4FaLQsvCsi7gfub9r3nMbreRHxuYjYGBEPRcQfjAS3iHhnRPxLRHwkIjYD7x99v4g4JyJui4htEfFoRHw8IqY2Hc+I+PWIuL9xziciIhrH+iLiQxGxKSIeBF5zBD/nJRFxd+MeN0fE85tq+N2I2BAR/RFxX0S8KiIuAn4f+IXGiPRdTb+PTzd+lg0R8d8iou9Av4/Gvlub7vXiiLi90QZxe0S8uOnYzRHx3yPiX4CdwDPaXCLi+Y3ztjV+nksa+/8I+MOmei8bdd35EbE7IvY1jg9ExHBE7Gr1dyqpdxmQJdVaRBwHvBr4t6bdbwDOBU4e45KPAfMowtrLgLcDv9x0/FzgQYp/2v/vY1y/D/htYDHwIuBVwP836pzXAmcDpwM/D1zY2P9rjWNnUowGv6mFHxF4+s8ZEc8Fvgi8G1gC3AD8Q0RMjYjnAb8JnJ2Zcxr3XpeZ/wT8KXBNo+XkjMZbfwYYAp7TqOsC4Fdb+X1ExELgK8BHgUXAh4GvRMSiptPeBlwOzAEeGnX9FOAfgK8CRwH/Efh8RDwvM983qt5PN1+bmd9u/JzfbhyfDTwCXNzK71NSbzMgS6qr6yJiG3Ar8C2KMDXizzJzS2Y+bTSxMTL6FuD3MrM/M9cBf0ER4kY8kpkfy8yh0dcDZOYdmfndxvF1wKcognazD2Tmtsz8KfBNYGVj/88Df5mZD2fmFuDPDvPn/AXgK5n5tczcC3wImAG8mCLATwNOjogpmbkuM3881htHxFKK0P3uzHwyM58APtL4HbXy+3gNcH9m/nXj+BeBe4HXNZ3zmcy8u3F876jrXwjMbvy+9mTmN4B/BC5t4fciSYetG/rnJKkd3pCZ/3yAYw8fYP9iYApPH8l8CFjWwrUANEZvP0wxAjyT4r+zd4w67bGm1zspQiDAsaPe/2kjqgfwjJ8zIo5tvjYzhyPiYWBZZt4cEe+maA85JSJuAt6TmY+M8d4nUPw+Hm10gUAxsNJc48F+H0+ro2E8v89jgYczc/gg10vShHMEWVIvygPs3wTspQiGI44HNrRw7YhPUoySnpSZcyn6euPgl+z3KHDcqHsfjkdo+hkaPc7H0fg5MvMLmfmSxjkJ/Hnj1NE/28PAILA4M+c3vuZm5ilN5xzs9/G0OhrG8/t8BDhu1MN7o6+XpAlnQJakhszcB3wJ+O8RMSciTgDeA4xn6rM5wA5gICJ+BviNcVz7JeA/RcTyiFgAXDGOa0e/z2saD99NAX6HIuh+JyKeFxGvjIhpwG5gFzAyQvs4sGIkkGbmoxT9v38REXMjYlJEPDsiRreMHMgNwHMj4hcjYnJE/AJF3/c/tnj99yhG2P9LREyJiJdTtGf8TYvXj7YbWHCY10rqIQZkSXq6/wg8SfHg2a3AF4Crx3H9fwZ+EegH/gq4ZhzX/hVwE3AXcCdw7Tiu3S8z7wPeSvHA4SaKUPm6zNxD0X/8gcb+xygefvu9xqUji21sjog7G6/fDkwFfghsBf4OOKbFOjZTPHT4O8Bm4L8Ar83MTS1ev6dR+8WNev8X8PbMvLeV68dwA/CliBjr4UxJ2i8yD/WvhZIkSVLvcARZkiRJamJAliRJkpoYkCVJkqQmBmRJkiSpSa0WClm8eHGuWLGi7DIkSZJUAXfcccemzFwyen+tAvKKFStYvXp12WVIkiSpAiJizBVLbbGQJEmSmhiQJUmSpCYGZEmSJKlJrXqQJUmS6mrv3r2sX7+e3bt3l11K5UyfPp3ly5czZcqUls5vW0COiOOAzwFLgQSuysz/GRELgWuAFcA64Oczc+sY178D+IPGt/8tMz/brlolSZK63fr165kzZw4rVqwgIsoupzIyk82bN7N+/XpOPPHElq5pZ4vFEPA7mXky8ELgXRFxMnAF8PXMPAn4euP7p2mE6PcB5wLnAO+LiAVtrFWSJKmr7d69m0WLFhmOxykiWLRo0bhG3tsWkDPz0cy8s/G6H7gHWAa8HhgZDf4s8IYxLr8Q+FpmbmmMLn8NuKhdtUqSJFWB4fjwjPf31pGH9CJiBXAm8D1gaWY+2jj0GEULxmjLgIebvl/f2DfWe18eEasjYvXGjRsnrGZJkiT1prYH5IiYDfw98O7M3NF8LDOToj/5sGXmVZm5KjNXLVnyjIVQJEmSNEH6+vpYuXIlp556Km9+85vZuXPnmOe9+MUv7nBlE6utATkiplCE489n5rWN3Y9HxDGN48cAT4xx6QbguKbvlzf2SZIkqSQzZsxgzZo1rF27lqlTp3LllVc+7fjQ0BAA3/nOd1p+z5FruknbAnIUzR6fBu7JzA83HboeeEfj9TuA/zfG5TcBF0TEgsbDeRc09kmSJKkLnH/++TzwwAPcfPPNnH/++VxyySWcfPLJAMyePRsoZpB473vfy6mnnsppp53GNddcAzDmNd2knfMgnwe8DfhBRKxp7Pt94APAlyLiMuAh4OcBImIV8OuZ+auZuSUi/gS4vXHdH2fmljbWKkmSVBl/9A9388NHdhz6xHE4+di5vO91p7R07tDQEDfeeCMXXVTMoXDnnXeydu3aZ0yjdu2117JmzRruuusuNm3axNlnn81LX/rSg17TDdoWkDPzVuBAjwy+aozzVwO/2vT91cDV7alOkiRJ47Vr1y5WrlwJFCPIl112Gd/5znc455xzxgy6t956K5deeil9fX0sXbqUl73sZdx+++3MnTv3gNd0A1fSkyRJqphWR3on2kgP8mizZs0a93sdzjWd0pFp3tSFHlsLn3ktfOt/lF2JJEmqqfPPP59rrrmGffv2sXHjRm655RbOOeecsss6JEeQe9HwMFx7OTxxN6z7Niw7C57zjK4XSZKkI/LGN76R2267jTPOOIOI4IMf/CBHH3009957b9mlHVQUUxHXw6pVq3L16tVll9H9fvo9uPoCeN3/hG/+KRx7Fvzi35RdlSRJOoh77rmH5z//+WWXUVlj/f4i4o7MXDX6XEeQe9EPr4O+aXDKz8Gm++F7n4I9T8LU7u0FkiRJ6hR7kHvRulvh+HNh+lx49itheC881PqE3pIkSXVmQO41e56Ex++G484tvj/hxTBpStGLLEmSJANyz9lwJ+S+pwLylBmw9GR45JlTtkiSJPUiA3KvefzuYnvMGU/tO2YlPHoX1OiBTUmSpMNlQO41m+6D6fNh1pKn9h27EnZvg63rSitLkiSpWxiQe83GH8GS50E0rQK+9LTGse6ek1CSJJWrr6+PlStXcuqpp/LmN7+ZnTt3juv6devW8YUvfOGw7v3iF7/4sK47HAbkXrPpPlj83KfvW/ycxrEfdb4eSZJUGSNLTa9du5apU6dy5ZVXjuv6wwnIQ0NDAHznO63PuDVyzeEyIPeSnVvgyY3PDMgzFhQtFwZkSZLUovPPP58HHniALVu28IY3vIHTTz+dF77whXz/+98H4Fvf+hYrV65k5cqVnHnmmfT393PFFVfw7W9/m5UrV/KRj3yEffv28d73vpezzz6b008/nU996lMA3HzzzZx//vlccsklnHzyyQDMnj0bgMzkve99L6eeeiqnnXYa11xzzQGvOVwuFNJLtjxYbBc955nHFj8XNj3Q2XokSdLhufEKeOwHE/ueR58GF3+gpVOHhoa48cYbueiii3jf+97HmWeeyXXXXcc3vvEN3v72t7NmzRo+9KEP8YlPfILzzjuPgYEBpk+fzgc+8AE+9KEP8Y//+I8AXHXVVcybN4/bb7+dwcFBzjvvPC644AIA7rzzTtauXcuJJ574tHtfe+21rFmzhrvuuotNmzZx9tln89KXvvSg14yXI8i9ZPvDxXb+8c88tug5jiBLkqSD2rVrFytXrmTVqlUcf/zxXHbZZdx666287W1vA+CVr3wlmzdvZseOHZx33nm85z3v4aMf/Sjbtm1j8uRnjst+9atf5XOf+xwrV67k3HPPZfPmzdx///0AnHPOOWMG3VtvvZVLL72Uvr4+li5dyste9jJuv/32g14zXo4g95JtIwH5uGceW3wS7NpStGHMXNjZuiRJ0vi0ONI70UZ6kFtxxRVX8JrXvIYbbriB8847j5tuuukZ52QmH/vYx7jwwguftv/mm29m1qxZ467vcK4ZiyPIvWT7wzBtHkyf98xj808ottt+2tmaJElSpZ1//vl8/vOfB4pgu3jxYubOncuPf/xjTjvtNH73d3+Xs88+m3vvvZc5c+bQ39+//9oLL7yQT37yk+zduxeAH/3oRzz55JOHvN8111zDvn372LhxI7fccgvnnHPOhP5MjiD3km0Pjz16DE+1XWx/uJgXWZIkqQXvf//7+ZVf+RVOP/10Zs6cyWc/+1kA/vIv/5JvfvObTJo0iVNOOYWLL76YSZMm0dfXxxlnnME73/lOfuu3fot169Zx1llnkZksWbKE66677qD3e+Mb38htt93GGWecQUTwwQ9+kKOPPpp775246Woja7R62qpVq3L16tVll9G9PnkezFsOv3jNM4/t3AIfPBEu/FN40bs6X5skSTqoe+65h+c///lll1FZY/3+IuKOzFw1+lxbLHrJtodh3gFGkGcsgKlzbLGQJEk9z4DcK3Zvh8HtB26xiCjaLAzIkiSpxxmQe8WOR4rt3GUHPseALElSV6tTa2wnjff3ZkDuFQOPF9s5Rx/4nJGA7B8+SZK6zvTp09m8ebMheZwyk82bNzN9+vSWr3EWi14x8ESxnb30wOfMPw4GdxTtGDPmd6YuSZLUkuXLl7N+/Xo2btxYdimVM336dJYvX97y+QbkXjEygjz7qAOfM+eYYtv/mAFZkqQuM2XKlAlZJU6HZotFrxh4HCZPh2lzD3zO/oD8aGdqkiRJ6kIG5F7R/3jRXhFx4HPmGpAlSZIMyL1i4PGD9x+DI8iSJEkYkHvHwBMH7z8GmDIDps+HHQZkSZLUuwzIvaKVEWQoRpEdQZYkST3MgNwLhvbAri2tBeS5BmRJktTbDMi94MnGfImHarGAxgjyY+2tR5IkqYsZkHvB/jmQW22xeAyG97W3JkmSpC7VtoVCIuJq4LXAE5l5amPfNcDzGqfMB7Zl5soxrl0H9AP7gKHMXNWuOnvCyAjyrCWHPnfO0ZD7imsOtiy1JElSTbVzJb3PAB8HPjeyIzN/YeR1RPwFsP0g178iMze1rbpesnNLsZ258NDnzj222PY/akCWJEk9qW0BOTNviYgVYx2LiAB+Hnhlu+6vJjs3F9uZiw597qxGn/KA67xLkqTeVFYP8vnA45l5/wGOJ/DViLgjIi7vYF31tGsLRB9Mn3foc2c32jCefKK9NUmSJHWpdrZYHMylwBcPcvwlmbkhIo4CvhYR92bmLWOd2AjQlwMcf/zxE19pHezcXLRXHGyZ6REjfcpPOoIsSZJ6U8dHkCNiMvBzwDUHOiczNzS2TwBfBs45yLlXZeaqzFy1ZEkLD6H1op1bWmuvAJg6C6bMssVCkiT1rDJaLP4dcG9mrh/rYETMiog5I6+BC4C1HayvfnZugRktPKA3YvYSWywkSVLPaltAjogvArcBz4uI9RFxWePQWxjVXhERx0bEDY1vlwK3RsRdwL8CX8nMf2pXnT1hpMWiVbOWwIABWZIk9aZ2zmJx6QH2v3OMfY8Ar268fhA4o1119aRdW2Dm2a2fP+so2LqubeVIkiR1M1fSq7vMxghyiz3IYIuFJEnqaQbkuhvsh+Gh8fUgzzqqCNUuNy1JknqQAbnuxrNIyIhZSyCHn7pWkiSphxiQ627XyDLT42yxAOdCliRJPcmAXHc7RwLyOFsswJksJElSTzIg193OwxhBdjU9SZLUwwzIdTfSRzxjQevXjLRYOIIsSZJ6kAG57nZtgZgE0+e3fs30+dA31RFkSZLUkwzIdbdzSxF4J43jo44o2iwMyJIkqQcZkOtu93aYMY7R4xEzFzrNmyRJ6kkG5LrbvR2mzxv/dTMXGZAlSVJPMiDX3e5t4+s/HmFAliRJPcqAXHeOIEuSJI2LAbnudm07zB7kRUW43rd34muSJEnqYgbkujuSEWSAXVsnth5JkqQuZ0Cus727Yd/gYQbkxtLUtllIkqQeY0Cus93biu3hPqQHBmRJktRzDMh1tnt7sT2SFgsDsiRJ6jEG5Drb1RhBPtyH9MCALEmSeo4Buc72jyAfRkCeYQ+yJEnqTQbkOtvfg3wYLRZTpsPU2bBzy8TWJEmS1OUMyHV2JCPIUMxk4QiyJEnqMQbkOjuSEWRwNT1JktSTDMh1tmsbTJkJk6ce3vUGZEmS1IMMyHV2uKvojTAgS5KkHmRArrPd2w6//xgaAdmH9CRJUm8xINfZEY8gL4Q9A8WS1ZIkST3CgFxnRxyQFxfbXY4iS5Kk3mFArrNd2w5vFb0RrqYnSZJ6kAG5zibiIT0wIEuSpJ5iQK6r4eFGQHYEWZIkaTwMyHW1px/II39ID5zJQpIk9RQDcl3tX2b6CALyjAXFdtfWI69HkiSpIgzIdTXYX2ynzTn89+ibAlPnGJAlSVJPaVtAjoirI+KJiFjbtO/9EbEhItY0vl59gGsvioj7IuKBiLiiXTXW2kQEZChGkW2xkCRJPaSdI8ifAS4aY/9HMnNl4+uG0Qcjog/4BHAxcDJwaUSc3MY662lwoNgecUCe7wiyJEnqKW0LyJl5C3A4Q4/nAA9k5oOZuQf4G+D1E1pcLxjcUWyPNCDPXGhAliRJPaWMHuTfjIjvN1owFoxxfBnwcNP36xv7xhQRl0fE6ohYvXHjxomutbr2NEaQp84+sveZscCALEmSekqnA/IngWcDK4FHgb840jfMzKsyc1VmrlqyZMmRvl19TGQPsktNS5KkHtLRgJyZj2fmvswcBv6Kop1itA3AcU3fL2/s03gMTtQIcqPFYnj4yGuSJEmqgI4G5Ig4punbNwJrxzjtduCkiDgxIqYCbwGu70R9tTK4A6bMhL7JR/Y+MxZADjcWHpEkSaq/I0xPBxYRXwReDiyOiPXA+4CXR8RKIIF1wH9onHss8L8z89WZORQRvwncBPQBV2fm3e2qs7b2DBz56DE8tVjIzi1HtuiIJElSRbQtIGfmpWPs/vQBzn0EeHXT9zcAz5gCTuMw2H/k/cfw1HLTu7YCJx75+0mSJHU5V9Krq8EBmDaBI8jOZCFJknqEAbmuBvth2twjfx8DsiRJ6jEG5Lra0z9BPcjNLRaSJEn1Z0Cuq4nqQZ4xv9gakCVJUo8wINfVRPUg902BqXOKWSwkSZJ6gAG5riZqBBlgpstNS5Kk3mFArqOhPbBvsBj5nQgzDMiSJKl3GJDraE9jmemJGkGesQB22WIhSZJ6gwG5jgYby0JPRA8yFDNZOIIsSZJ6hAG5jvYHZFssJEmSxsuAXEdtabHYCsPDE/N+kiRJXcyAXEcjI8gT9ZDezIWQwzC4Y2LeT5IkqYsZkOuoHS0WYJuFJEnqCQbkOprwh/QMyJIkqXcYkOuoHT3I4FRvkiSpJxiQ62h/D/IETvMGsGvbxLyfJElSFzMg19FgP0yZBZP6Jub9bLGQJEk9xIBcR4P9E9d/DDBjfrHdaYuFJEmqPwNyHQ32T1z/MUDfFJg21xFkSZLUEwzIdbRnYOL6j0fMmG9AliRJPcGAXEcTPYIMjdX0bLGQJEn1Z0Cuo8GBNgTkhY4gS5KknmBArqPBHW0aQTYgS5Kk+jMg11FbepAXOIuFJEnqCQbkOmpHD/LMhbB7GwwPT+z7SpIkdRkDct0MDcK+PRM7DzIUI8g5DIPbJ/Z9JUmSuowBuW4GB4rttLkT+76upidJknqEAblu9vQX2wnvQV5YbA3IkiSp5gzIdTPYCMjtmMUCDMiSJKn2DMh1s7/Fog09yAA7DciSJKneDMh1s38EeYJ7kGfaYiFJknqDAblu2tWDPH1+sTUgS5KkmjMg1027epD7Jhej0gZkSZJUc20LyBFxdUQ8ERFrm/b9j4i4NyK+HxFfjoj5B7h2XUT8ICLWRMTqdtVYS+3qQYbGctOupidJkuqtnSPInwEuGrXva8CpmXk68CPg9w5y/Ssyc2VmrmpTffU02KYWC2gEZEeQJUlSvbUtIGfmLcCWUfu+mplDjW+/Cyxv1/171p4BmDILJvVN/HsbkCVJUg8oswf5V4AbD3Asga9GxB0RcfnB3iQiLo+I1RGxeuPGjRNeZOUM7pj4/uMRMxbATlssJElSvZUSkCPivwJDwOcPcMpLMvMs4GLgXRHx0gO9V2ZelZmrMnPVkiVL2lBtxQwOtKf/GIqp3hxBliRJNdfxgBwR7wReC/xSZuZY52Tmhsb2CeDLwDkdK7DqBvvbO4K8exsMD7fn/SVJkrpARwNyRFwE/BfgkszceYBzZkXEnJHXwAXA2rHO1Rj2DLTnAT0oAnIOw+D29ry/JElSF2jnNG9fBG4DnhcR6yPiMuDjwBzga40p3K5snHtsRNzQuHQpcGtE3AX8K/CVzPyndtVZO4P9E7+K3ogZrqYnSZLqb3K73jgzLx1j96cPcO4jwKsbrx8EzmhXXbU32N++HuQZC4qtAVmSJNWYK+nVTbt7kAF2GpAlSVJ9GZDrpp09yDNtsZAkSfVnQK6ToUHYt6f9I8gGZEmSVGMG5DoZWWa6XQF5+vxiu8vFQiRJUn0ZkOuk3QG5bzJMm+cIsiRJqjUDcp2MBOR29SADzJhvQJYkSbVmQK6TPQPFtl0jyFD0Ie+0xUKSJNWXAblO2t1iAcVMFo4gS5KkGjMg10knAvKMBQZkSZJUawbkOulID/ICZ7GQJEm1ZkCuk470IC+EXdtgeLh995AkSSqRAblOOjWCTMLg9vbdQ5IkqUQG5DoZbCwzPamNH+vIanrOZCFJkmrKgFwngzvaO3oMxSwWULRZSJIk1ZABuU72DLS3/xieGkF2JgtJklRTBuQ6GezvYEC2xUKSJNWTAblOBgdgWptbLBxBliRJNWdArpPBfpg2t733mD6/2BqQJUlSTRmQ62RPf/sf0uubDNPmGZAlSVJtGZDrpBM9yAAz5jvNmyRJqi0Dcl1kdqYHGYqp3hxBliRJNWVArouhQRje26ER5AUGZEmSVFsG5LrYM1Bsp3YqINtiIUmS6smAXBeDO4ptR0aQbbGQJEn11VJAjohrI+I1EWGg7laDjRHkTvQgz1hQLDU9PNz+e0mSJHVYq4H3fwG/CNwfER+IiOe1sSYdjsH+YtupHmQSdm9r/70kSZI6rKWAnJn/nJm/BJwFrAP+OSK+ExG/HBFT2lmgWtTJHuSZC4utbRaSJKmGWm6ZiIhFwDuBXwX+DfifFIH5a22pTOPT8RFkijYLSZKkmpncykkR8WXgecBfA6/LzEcbh66JiNXtKk7jsD8gd6gHGZzJQpIk1VJLARn4q8y8oXlHREzLzMHMXNWGujReHR1BtsVCkiTVV6stFv9tjH23TWQhOkIjAXnKrPbfa/8IsgFZkiTVz0FHkCPiaGAZMCMizgSicWguMLPNtWk89gwUD+hN6sBMfNPnFdudtlhIkqT6OVSLxYUUD+YtBz7ctL8f+P021aTDMbijM/3HAH2Ti5DsCLIkSaqhgw43ZuZnM/MVwDsz8xVNX5dk5rWHevOIuDoinoiItU37FkbE1yLi/sZ2wQGufUfjnPsj4h3j/sl6zeBAZ/qPR8xYYECWJEm1dNCAHBFvbbxcERHvGf3Vwvt/Brho1L4rgK9n5knA1xvfj77vQuB9wLnAOcD7DhSk1TDYD1M7NIIMjYBsi4UkSaqfQzWsjjzxNRuYM8bXQWXmLcDoFPV64LON158F3jDGpRcCX8vMLZm5lWKu5dFBW832dHoEeaEjyJIkqZYO2oOcmZ9qbP9oAu+5tGke5ceApWOcswx4uOn79Y19zxARlwOXAxx//PETWGbFDPbDrCWdu9+MBbD1J527nyRJUoe0NOVBRHwwIuZGxJSI+HpEbGxqvzhsmZlAHuF7XJWZqzJz1ZIlHQyI3aaMHmRnsZAkSTXU6pxgF2TmDuC1wDrgOcB7D/Oej0fEMQCN7RNjnLMBOK7p++WNfTqQwR2d7UGeuRB2b4fhfZ27pyRJUge0GpBHWjFeA/xtZm4/gnteD4zMSvEO4P+Ncc5NwAURsaDxcN4FjX0aS2YJPcgLgCxCsiRJUo20GpD/MSLuBV4AfD0ilgC7D3VRRHyRYsW950XE+oi4DPgA8LMRcT/w7xrfExGrIuJ/A2TmFuBPgNsbX3/c2KexDO2G4aHOzYMMrqYnSZJq61ALhQCQmVdExAeB7Zm5LyKepJiN4lDXXXqAQ68a49zVwK82fX81cHUr9fW8wYFiO21u5+5pQJYkSTXVUkBu+BmK+ZCbr/ncBNejwzG4o9h2dB7khcXWgCxJkmqmpYAcEX8NPBtYA4w8lZUYkLvDnpER5E73IGNAliRJtdPqCPIq4OTGtGzqNoP9xbaMHmSnepMkSTXT6kN6a4Gj21mIjsBgGSPI84utI8iSJKlmWh1BXgz8MCL+FRgc2QYKi7AAABm0SURBVJmZl7SlKo3PyAjy1A4G5El9MH2eAVmSJNVOqwH5/e0sQkdoz0iLRQcDMhRtFrtssZAkSfXS6jRv34qIE4CTMvOfI2Im0Nfe0tSyMnqQoZjJwhFkSZJUMy31IEfErwF/B3yqsWsZcF27itI4DQ4AAVNmdfa+MxYYkCVJUu20+pDeu4DzgB0AmXk/cFS7itI4DfYXcyBPavXjnCAzF8LOzZ29pyRJUpu1mqgGM3PPyDeNxUKc8q1b7OnvfP8xwMxFTvMmSZJqp9WA/K2I+H1gRkT8LPC3wD+0ryyNy2B/5/uPoQjIgztgaM+hz5UkSaqIVgPyFcBG4AfAfwBuAP6gXUVpnAZLHEEG2ywkSVKttDqLxXBEXAdcl5kb21yTxmtwoPyAPPeYzt9fkiSpDQ46ghyF90fEJuA+4L6I2BgRf9iZ8tSSskaQZy0uto4gS5KkGjlUi8VvU8xecXZmLszMhcC5wHkR8dttr06tGeyHaXM7f9/9I8ibOn9vSZKkNjlUQH4bcGlm/mRkR2Y+CLwVeHs7C9M4jEzz1mn7A7IzWUiSpPo4VECekpnPGB5s9CFPaU9JGpfM8qZ5m7Gw2D7pCLIkSaqPQwXkg83f5dxe3WDvTsjhcgJy32SYPt8eZEmSVCuHmsXijIjYMcb+AKa3oR6N12B/sS0jIEPxoJ4BWZIk1chBA3Jm9nWqEB2msgPyzEU+pCdJkmql1YVC1K26IiD7kJ4kSaoPA3LVdUNA9iE9SZJUIwbkquuGgLxzczGbhiRJUg0YkKuu7IA8azEM732qDkmSpIozIFfdnoFiO7XEEWTwQT1JklQbBuSqG2zMwldmiwX4oJ4kSaoNA3LVDfbDpCkweVo595+5uNj6oJ4kSaoJA3LVDTaWmY4o5/4zG8tNu1iIJEmqCQNy1Q0OwLTZ5d1/VmME2R5kSZJUEwbkqhvsh2lzy7v/1NnQN9URZEmSVBsG5Kob3FHeA3pQtHbMXGxAliRJtWFArrqRHuQyzVwETxqQJUlSPRiQq27PQBcE5IWOIEuSpNroeECOiOdFxJqmrx0R8e5R57w8IrY3nfOHna6zMgb7iz7gMs1a7EN6kiSpNiZ3+oaZeR+wEiAi+oANwJfHOPXbmfnaTtZWSd3SYuEIsiRJqomyWyxeBfw4Mx8quY5q2jcEe3eWO4sFFA/p7d4O+/aWW4ckSdIEKDsgvwX44gGOvSgi7oqIGyPilE4WVRl7Bopt6SPII4uFuNy0JEmqvtICckRMBS4B/naMw3cCJ2TmGcDHgOsO8j6XR8TqiFi9cePG9hTbrQb7i23ZAXlksZAne+z3L0mSaqnMEeSLgTsz8/HRBzJzR2YONF7fAEyJiMVjvUlmXpWZqzJz1ZIlS9pbcbfZH5DLfkiv8Xs3IEuSpBooMyBfygHaKyLi6IiIxutzKOr0KbDRumYE+ahia0CWJEk10PFZLAAiYhbws8B/aNr36wCZeSXwJuA3ImII2AW8JTOzjFq72p6RgFzyQ3qzGyPIA0+UW4ckSdIEKCUgZ+aTwKJR+65sev1x4OOdrqtyumUEefp86JsKTxqQJUlS9ZU9i4WOxEhALnuhkIiiD3nAFgtJklR9BuQq65YRZCgCsiPIkiSpBgzIVdZNAXn2UT6kJ0mSasGAXGWD/TBlFkzqK7sSWywkSVJtGJCrbLC/O0aPodFisRGcbESSJFWcAbnKBvvLXyRkxOyjYHgv7NpadiWSJElHxIBcZV01guxiIZIkqR4MyFU22F/+FG8jXCxEkiTVhAG5ygZ3wPR5ZVdR2D+CbECWJEnVZkCust1dFJBnNwKyM1lIkqSKMyBX2e7t3ROQZyyAmGQPsiRJqjwDclUN74M9/TBtbtmVFCb1wczFtlhIkqTKMyBX1cgqetO7JCBD0WZhi4UkSao4A3JV7d5ebLulxQIai4U4gixJkqrNgFxVgzuKbbe0WIAjyJIkqRYMyFW1uxGQu6nFYmQE2eWmJUlShRmQq6obWyxmHwVDu5/qj5YkSaogA3JVdWOLhctNS5KkGjAgV9X+FosuGkGes7TY9j9Wbh2SJElHwIBcVSMtFt00gjznmGLb/2i5dUiSJB0BA3JVDW6HyTNg8tSyK3nKnKOLrSPIkiSpwgzIVbV7R3fNYAEwfT5Mnu4IsiRJqjQDclXt3t5d/ccAEcUosiPIkiSpwgzIVTW4o7v6j0fMOcaALEmSKs2AXFW7t3dfiwUUI8gDBmRJklRdBuSq2r2j+1oswBFkSZJUeQbkquraFoujYc+Aq+lJkqTKMiBXVTc+pAdNcyE7iixJkqrJgFxFQ3tgaHf39iCDU71JkqTKMiBX0WBjmelpjiBLkiRNNANyFY0sM92NLRazlxZbR5AlSVJFGZCraH9A7sIWi2lzYMosR5AlSVJlGZCraH+LRRcG5P2r6TmCLEmSqsmAXEXd3GIBzoUsSZIqrbSAHBHrIuIHEbEmIlaPcTwi4qMR8UBEfD8iziqjzq60uzGC3I0tFuAIsiRJqrTJJd//FZm56QDHLgZOanydC3yysdXubcV2+vxy6ziQOUcXI8iZRcuFJElShXRzi8Xrgc9l4bvA/Ig4puyiusKurRB9xQNx3Wje8mKe5p1byq5EkiRp3MoMyAl8NSLuiIjLxzi+DHi46fv1jX1PExGXR8TqiFi9cePGNpXaZXZthRnzu3d0dm7jY9qxvtw6JEmSDkOZAfklmXkWRSvFuyLipYfzJpl5VWauysxVS5YsmdgKu9WubTBjQdlVHNi85cV2uwFZkiRVT2kBOTM3NLZPAF8Gzhl1ygbguKbvlzf2adfWLg/IjY/NgCxJkiqolIAcEbMiYs7Ia+ACYO2o064H3t6YzeKFwPbMdGoEKB7S69YH9ABmLYa+aQZkSZJUSWXNYrEU+HIUPbSTgS9k5j9FxK8DZOaVwA3Aq4EHgJ3AL5dUa/fZtRUWnVR2FQcWAfOWGZAlSVIllRKQM/NB4Iwx9l/Z9DqBd3WyrsrYta14SK+bzVtuQJYkSZXUzdO8aSzD+4qV9Lq5Bxlg7nLYYcu4JEmqHgNy1ezeDmR39yBDMYLc/yjsGyq7EkmSpHExIFfNyCp63T6CPG855LBLTkuSpMoxIFfNrq3Ftut7kBuLhdiHLEmSKsaAXDW7qjKC7FzIkiSpmgzIVbN/BLnLA7LLTUuSpIoyIFfNSA9ytz+kN212UaMjyJIkqWIMyFVTlR5kKNostj1cdhWSJEnjYkCuml3bYMpMmDyt7EoObcEJsO2hsquQJEkaFwNy1eza1v39xyMWrICt6yCz7EokSZJaZkCuml1bu7//eMSCFTC0GwYeL7sSSZKklhmQq2Z3lUaQTyy2W9eVWoYkSdJ4GJCrZufmajygB8UIMhiQJUlSpRiQq2bnZpi1uOwqWjP/OCBgy0/KrkSSJKllBuQqGR6GnVtgZkUC8uRpxYIhjiBLkqQKMSBXye5tkPtg5qKyK2ndyEwWkiRJFWFArpKdm4ttVVosABauMCBLkqRKMSBXyUhAnrmw3DrGY8EKGHgM9uwsuxJJkqSWGJCr5MlNxbYqPcjw1FRv235abh2SJEktMiBXyf4R5Ir1IANsdSYLSZJUDQbkKtnZGEGuVA/ys4rt5gfKrUOSJKlFBuQqeXIzTJkFU2aUXUnrZi4sRrw33V92JZIkSS0xIFfJzs3Vaq8Ysfi5BmRJklQZBuQq2bkJZlUxIJ8Em35UdhWSJEktMSBXSVVHkBedVIT7nVvKrkSSJOmQDMhV8uTmak3xNmLxc4utD+pJkqQKMCBXSVVHkBefVGxts5AkSRVgQK6Kvbtg75PVWkVvxPwTYNIUA7IkSaoEA3JVDDxebGcvLbeOw9E3GRY9GzbZYiFJkrqfAbkqBp4otlUMyNCYyeK+squQJEk6JANyVYyMIM+paEBe8jOw5Sewd3fZlUiSJB2UAbkqqtxiAbD0FMh9jiJLkqSuZ0Cuiv7HgajmNG8AS08tto+tLbcOSZKkQ+h4QI6I4yLimxHxw4i4OyJ+a4xzXh4R2yNiTePrDztdZ9cZeBxmLSkeeKuihc+CydPh8bvLrkSSJOmgykhbQ8DvZOadETEHuCMivpaZPxx13rcz87Ul1NedBp6obnsFwKQ+OOr58LgjyJIkqbt1fAQ5Mx/NzDsbr/uBe4Blna6jcgYeg9lHlV3FkVl6ahGQM8uuRJIk6YBK7UGOiBXAmcD3xjj8ooi4KyJujIhTDvIel0fE6ohYvXHjxjZV2gUGnoA5R5ddxZFZemqxGuDIA4eSJEldqLSAHBGzgb8H3p2ZO0YdvhM4ITPPAD4GXHeg98nMqzJzVWauWrJkSfsKLtPwcKPFouojyI2/59hmIUmSulgpATkiplCE489n5rWjj2fmjswcaLy+AZgSERWdvmEC7N4Gw3ur3YMMTwXkx35Qbh2SJEkHUcYsFgF8GrgnMz98gHOObpxHRJxDUefmzlXZZfofK7ZVH0GeuRAWrIANd5ZdiSRJ0gGVMYvFecDbgB9ExJrGvt8HjgfIzCuBNwG/ERFDwC7gLZk9/GTXjkeK7Zxjy61jIix7Afx0rJZzSZKk7tDxgJyZtwJxiHM+Dny8MxVVwI71xXbe8nLrmAjLXgBr/74YFa/6Q4eSJKmWXEmvCrZvgJgEc44pu5Ijt+wFxdY2C0mS1KUMyFWwfT3MPrq6q+g1O/p0iD7YcEfZlUiSJI3JgFwFO9bXo70CYOrMYjYLA7IkSepSBuQq2L4B5tVoscFlLyhaLIb3lV2JJEnSMxiQu10m7NgAc2sUkE94MQxud8EQSZLUlQzI3W7nZhjaXZ8WC4ATziu2624ttw5JkqQxGJC73fYaTfE2Yt4yWHAirPuXsiuRJEl6BgNyt9v+cLGtU0AGWPESeOhfYHi47EokSZKexoDc7bY8WGwXnFhuHRNtxfmwe5t9yJIkqesYkLvdlp/AjIUwY37ZlUysFS8ptg/eXGoZkiRJoxmQu92WB2Hhs8quYuLNWwZHnQz3f7XsSiRJkp7GgNzttvwEFtasvWLEcy+En94Gu7aVXYkkSdJ+BuRuNjRYrKJXxxFkgOdeBMND8ONvlF2JJEnSfgbkbrbtp5DD9XtAb8Tys2HGAtssJElSVzEgd7PNPy62dR1BntQHJ10A990IQ3vKrkaSJAkwIHe3jfcU2yXPK7eOdjr13xfTvf3462VXIkmSBBiQu9sT98DcZfWb4q3Zs19ZTGP3g78ruxJJkiTAgNzdnvghHPX8sqtor74pcMob4L4bYHCg7GokSZKYXHYBVfbgxgHWPNyeKcpieIhLnriPH89axdo717flHt1i4bRX8PK9V3PHjVfz0PE/V3Y5kiSpw561ZDYrj+uefzE3IB+B2x7czH/9cnuWSn5WPMIbp+3hynumce3dd7XlHt1jEjdNXc60Oz/Ne777LCDKLkiSJHXQO150ggG5Ll53xrG85DmL2/LeM++7Fr4K//ntP8dvLTm1LffoJrN/8JssuvkKbnvrHPYcs6rsciRJUgfNntZdkbS7qqmYudOnMHf6lPa8+fa7YcpMjj3pBdDXAx/Ti94Gt/0px9zzf+DUl5VdjSRJ6mE+pNetNqyGY1b2RjgGmDYbzr4M7r6umL1DkiSpJAbkbjS0Bx79Pix/QdmVdNaL/xNMnQXf+vOyK5EkST3MgNyNHvk32DdYLMXcS2YuhHN/He7+Mjx8e9nVSJKkHmVA7kYPfhMIWHF+2ZV03kveDXOOga+8B/YNlV2NJEnqQQbkbvTjb8Cys4oR1V4zbQ5c9Gfw2Pfhu58ouxpJktSDDMjdZucWWL8anvWKsispz8lvgOe/Dr7+x7D+jrKrkSRJPcaA3G1++P8g98HJry+7kvJEwCUfg7nHwpfeDts3lF2RJEnqIQbkbvODv4PFz4WjTyu7knLNWAC/8H9h93b4vz8HT24uuyJJktQjDMjd5LG18NCtcMalxShqrzvmDLj0i7DlJ/Dpn4UtD5ZdkSRJ6gEG5G5y64dhyixY9ctlV9I9Tjwf3nE97NoCV728GGGXJElqIwNyt/jJt2Ht38MLf6NoL9BTjn8h/No3itaTv78M/vqN8MiasquSJEk1VUpAjoiLIuK+iHggIq4Y4/i0iLimcfx7EbGi81V20NaHiuC3YAWc/ztlV9OdFj4Lfvmf4MI/LRZSueplcPXFsOaL9idLkqQJNbnTN4yIPuATwM8C64HbI+L6zPxh02mXAVsz8zkR8Rbgz4Ff6HStbbdvCO65Hm78XRgahLddB1Nnll1V9+qbDC96F5z5Vlj9f+DOz8J1vw5EMW/0sWfBMafDkp+BuctgztEwqa/sqiVJUsVEZnb2hhEvAt6fmRc2vv89gMz8s6Zzbmqcc1tETAYeA5bkIYpdtWpVrl69un3Fj/boXfDAP0MOQ2ZjO+preN8zjw/vLUaNH10DOzfDUafAm66Go36mc7XXwfBwMZr8wNfgwW8Vi4vsGXjq+KTJMGsJTJsL0+fB9LkwdTZMngaTpkDfFOib2thOgWiE6QggntoecB9N+3hqvyRJGp9jV8KzX9nx20bEHZm5avT+jo8gA8uAh5u+Xw+ce6BzMnMoIrYDi4BNo98sIi4HLgc4/vjj21HvgW24o1jM4hlFTRr11df0OopRzbnL4KQL4GdeA8+9uBgd1fhMmgTLX1B8vfyKIjBv/Qls/jFsfxh2bICBx2H3DhjcUSzCsvWh4i8o+/bCvj2NbeN1DgNZ/GWGzv7FUZKknnbO5aUE5AOpfCrLzKuAq6AYQe7ozc98O5zxi0XwndT3VABWOSZNgkXPLr4mUjaF5v3/iDEqSHf4X2IkSaqVLmuJLCMgbwCOa/p+eWPfWOesb7RYzAO670msvsmO/PaCCP/iI0lSDyljFovbgZMi4sSImAq8Bbh+1DnXA+9ovH4T8I1D9R9LkiRJE6Hjw5+NnuLfBG4C+oCrM/PuiPhjYHVmXg98GvjriHgA2EIRoiVJkqS2K6U/IDNvAG4Yte8Pm17vBt7c6bokSZIkV9KTJEmSmhiQJUmSpCYGZEmSJKmJAVmSJElqYkCWJEmSmhiQJUmSpCYGZEmSJKmJAVmSJElqYkCWJEmSmhiQJUmSpCYGZEmSJKlJZGbZNUyYiNgIPNTh2y4GNnX4nuoMP9t68/OtLz/b+vKzra+yPtsTMnPJ6J21CshliIjVmbmq7Do08fxs683Pt778bOvLz7a+uu2ztcVCkiRJamJAliRJkpoYkI/cVWUXoLbxs603P9/68rOtLz/b+uqqz9YeZEmSJKmJI8iSJElSEwOyJEmS1MSA3KKIuCgi7ouIByLiijGOT4uIaxrHvxcRKzpfpQ5HC5/teyLihxHx/Yj4ekScUEadGr9DfbZN5/37iMiI6JophnRwrXy2EfHzjT+7d0fEFzpdow5PC/9NPj4ivhkR/9b47/Kry6hT4xcRV0fEExGx9gDHIyI+2vjsvx8RZ3W6xhEG5BZERB/wCeBi4GTg0og4edRplwFbM/M5wEeAP+9slTocLX62/wasyszTgb8DPtjZKnU4WvxsiYg5wG8B3+tshTpcrXy2EXES8HvAeZl5CvDujheqcWvxz+0fAF/KzDOBtwD/q7NV6gh8BrjoIMcvBk5qfF0OfLIDNY3JgNyac4AHMvPBzNwD/A3w+lHnvB74bOP13wGviojoYI06PIf8bDPzm5m5s/Htd4HlHa5Rh6eVP7cAf0LxF9rdnSxOR6SVz/bXgE9k5laAzHyiwzXq8LTy2SYwt/F6HvBIB+vTEcjMW4AtBznl9cDnsvBdYH5EHNOZ6p7OgNyaZcDDTd+vb+wb85zMHAK2A4s6Up2ORCufbbPLgBvbWpEmyiE/28Y/3x2XmV/pZGE6Yq38uX0u8NyI+JeI+G5EHGzUSt2jlc/2/cBbI2I9cAPwHztTmjpgvP9PbpvJZdxUqqKIeCuwCnhZ2bXoyEXEJODDwDtLLkXtMZnin2lfTvGvPrdExGmZua3UqjQRLgU+k5l/EREvAv46Ik7NzOGyC1N9OILcmg3AcU3fL2/sG/OciJhM8c8+mztSnY5EK58tEfHvgP8KXJKZgx2qTUfmUJ/tHOBU4OaIWAe8ELjeB/UqoZU/t+uB6zNzb2b+BPgRRWBWd2vls70M+BJAZt4GTAcWd6Q6tVtL/0/uBANya24HToqIEyNiKsVDAdePOud64B2N128CvpGuwlIFh/xsI+JM4FMU4dg+xuo46Gebmdszc3FmrsjMFRT95Zdk5upyytU4tPLf5OsoRo+JiMUULRcPdrJIHZZWPtufAq8CiIjnUwTkjR2tUu1yPfD2xmwWLwS2Z+ajZRRii0ULMnMoIn4TuAnoA67OzLsj4o+B1Zl5PfBpin/meYCiAf0t5VWsVrX42f4PYDbwt43nLn+amZeUVrRa0uJnqwpq8bO9CbggIn4I7APem5n+q16Xa/Gz/R3gryLityke2HunA1LVEBFfpPiL6+JGD/n7gCkAmXklRU/5q4EHgJ3AL5dTqUtNS5IkSU9ji4UkSZLUxIAsSZIkNTEgS5IkSU0MyJIkSVITA7IkSZLUxIAsSZIkNTEgS5IkSU0MyJJUQxHx0YjYFxFrIuIHEbEnIt5bdl2SVAUuFCJJNRURA5k5u7HU8urGktqSpENwBFmSJElqYkCWJEmSmhiQJan+hoCpZRchSVVhQJak+tsObIyIq8ouRJKqwIf0JEmSpCaOIEuSJElNDMiSJElSEwOyJEmS1MSALEmSJDUxIEuSJElNDMiSJElSEwOyJEmS1OT/B5iftTPN/xg1AAAAAElFTkSuQmCC\n","text/plain":["<Figure size 720x432 with 1 Axes>"]},"metadata":{"needs_background":"light"}}]},{"cell_type":"code","metadata":{"id":"jfIsDwv6WtBq","executionInfo":{"status":"ok","timestamp":1634010775480,"user_tz":240,"elapsed":5,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}}},"source":[""],"execution_count":196,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"pT9ie9lKodGF"},"source":["## Problem 2(iv) -  compute the prior and posterior probability that $\\pi$ takes a value above 0.15."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"HgamkFU7omPY","executionInfo":{"status":"ok","timestamp":1634012720768,"user_tz":240,"elapsed":127,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"867bb414-8662-47f7-bf41-964fe0001ad2"},"source":["# Find value for Beta(1, 1)\n","print(\"Beta(1, 1)=\", beta_func(1, 1))\n","# Values of a1 and b1\n","print(f\"alpha_prior={alpha_prior}, beta_prior={beta_prior}\")\n","print(f\"alpha_post={alpha_post}, beta_post={beta_post}\")\n","\n","# Calculate prior probability using integral\n","prior_I = quad(get_beta_pdf, 0.15, 1, args=(alpha_prior, beta_prior))\n","print(\"Integral for prior: \", prior_I[0])\n","# Calculate posterior probability using integral\n","post_I = quad(get_beta_pdf,  0.15, 1, args=(alpha_post, beta_post))\n","print(\"Integral for posterior: \", post_I[0])"],"execution_count":207,"outputs":[{"output_type":"stream","name":"stdout","text":["Beta(1, 1)= 1.0\n","alpha_prior=1, beta_prior=1\n","alpha_post=15, beta_post=163.0\n","Integral for prior:  0.85\n","Integral for posterior:  0.0033625668021733273\n"]}]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"-Htq4sjcwLJ6","executionInfo":{"status":"ok","timestamp":1634015028644,"user_tz":240,"elapsed":80,"user":{"displayName":"Kostas Georgiou","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GhzP5qrxeKXeuneTscU-PqTEiyBLJpOwQ-olOc2x5s=s64","userId":"10199933473437407446"}},"outputId":"d7c32dcc-0ba9-42a3-844d-10f12e739620"},"source":["beta_func(15, 163)"],"execution_count":214,"outputs":[{"output_type":"execute_result","data":{"text/plain":["3.0604098911703136e-23"]},"metadata":{},"execution_count":214}]}]}