{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "This page has been moved to [https://econ.pages.code.wm.edu/414/syllabus/docs/index.html](https://econ.pages.code.wm.edu/407/notes/docs/index.html) and is no longer being maintained here." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "from IPython.display import YouTubeVideo\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from scipy.stats import norm\n", "from scipy.stats import gamma\n", "from scipy.stats import uniform\n", "import statsmodels.formula.api as smf \n", "\n", "np.random.seed(12578)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Frequentists Statistics\n", "\n", "The statsmodels package has some of the functionality of stata. It can estimate most models you would see as a student taking econometrics, time series, and cross section here at William and Mary. Here is an example OLS regression:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ideducln_wagepexptimeabilitymeducfeducbroken_homesiblingspexp2
04122.14240.261210144
16151.91440.4412160216
28132.32840.5112151264
311141.64141.821617121
412132.1664-1.3013120536
\n", "
" ], "text/plain": [ " id educ ln_wage pexp time ability meduc feduc broken_home \\\n", "0 4 12 2.14 2 4 0.26 12 10 1 \n", "1 6 15 1.91 4 4 0.44 12 16 0 \n", "2 8 13 2.32 8 4 0.51 12 15 1 \n", "3 11 14 1.64 1 4 1.82 16 17 1 \n", "4 12 13 2.16 6 4 -1.30 13 12 0 \n", "\n", " siblings pexp2 \n", "0 4 4 \n", "1 2 16 \n", "2 2 64 \n", "3 2 1 \n", "4 5 36 " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# load data and create dataframe\n", "tobias_koop=pd.read_csv('https://rlhick.people.wm.edu/econ407/data/tobias_koop_t_4.csv')\n", "\n", "tobias_koop.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that statsmodels works seemlessly with Pandas dataframes (tobias_koop)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: ln_wage R-squared: 0.166\n", "Model: OLS Adj. R-squared: 0.163\n", "Method: Least Squares F-statistic: 51.36\n", "Date: Mon, 12 Mar 2018 Prob (F-statistic): 1.83e-39\n", "Time: 15:39:05 Log-Likelihood: -583.66\n", "No. Observations: 1034 AIC: 1177.\n", "Df Residuals: 1029 BIC: 1202.\n", "Df Model: 4 \n", "Covariance Type: nonrobust \n", "===============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 0.4603 0.137 3.353 0.001 0.191 0.730\n", "educ 0.0853 0.009 9.179 0.000 0.067 0.104\n", "pexp 0.2035 0.024 8.629 0.000 0.157 0.250\n", "pexp2 -0.0124 0.002 -5.438 0.000 -0.017 -0.008\n", "broken_home -0.0087 0.036 -0.244 0.807 -0.079 0.061\n", "==============================================================================\n", "Omnibus: 55.892 Durbin-Watson: 1.761\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 112.050\n", "Skew: -0.355 Prob(JB): 4.66e-25\n", "Kurtosis: 4.448 Cond. No. 391.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "formula = \"ln_wage ~ educ + pexp + pexp2 + broken_home\"\n", "results = smf.ols(formula,tobias_koop).fit()\n", "print(results.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the results object contains the other types of information (including robust standard errors). " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Intercept 0.131223\n", "educ 0.009905\n", "pexp 0.022629\n", "pexp2 0.002239\n", "broken_home 0.031162\n", "dtype: float64" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results.HC0_se" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to specify which variance covariance matrix to use, use the `cov_type` argument." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: ln_wage R-squared: 0.166\n", "Model: OLS Adj. R-squared: 0.163\n", "Method: Least Squares F-statistic: 65.13\n", "Date: Mon, 12 Mar 2018 Prob (F-statistic): 3.89e-49\n", "Time: 15:39:05 Log-Likelihood: -583.66\n", "No. Observations: 1034 AIC: 1177.\n", "Df Residuals: 1029 BIC: 1202.\n", "Df Model: 4 \n", "Covariance Type: HC0 \n", "===============================================================================\n", " coef std err z P>|z| [0.025 0.975]\n", "-------------------------------------------------------------------------------\n", "Intercept 0.4603 0.131 3.508 0.000 0.203 0.718\n", "educ 0.0853 0.010 8.609 0.000 0.066 0.105\n", "pexp 0.2035 0.023 8.994 0.000 0.159 0.248\n", "pexp2 -0.0124 0.002 -5.543 0.000 -0.017 -0.008\n", "broken_home -0.0087 0.031 -0.280 0.779 -0.070 0.052\n", "==============================================================================\n", "Omnibus: 55.892 Durbin-Watson: 1.761\n", "Prob(Omnibus): 0.000 Jarque-Bera (JB): 112.050\n", "Skew: -0.355 Prob(JB): 4.66e-25\n", "Kurtosis: 4.448 Cond. No. 391.\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors are heteroscedasticity robust (HC0)\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/robhicks/anaconda/envs/python36/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n", " from pandas.core import datetools\n" ] } ], "source": [ "results = smf.ols(formula,tobias_koop).fit(cov_type='HC0')\n", "print(results.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Statistical Distributions\n", "\n", "In the course, we will be evaluating cdf and pdf's, building likelihood and log-likelihood functions, taking random samples and perhaps other statistical operations. We will primarily use the `scipy` family of distributions. At times, I may \"slip-up\" and use some `numpy` functions, but in my opinion the `scipy` functions are more powerful by saving time coding. Here are some examples:\n", "\n", "### Evaluating pdf's and cdf's:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PDF value\n", "0.12951759566589174\n", "CDF value\n", "0.9331927987311419\n" ] } ], "source": [ "mean = 0\n", "y = 1.5\n", "std = 1\n", "print(\"PDF value\") \n", "print(norm(mean,std).pdf(y))\n", "print(\"CDF value\")\n", "print(norm(mean,std).cdf(y))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAHVCAYAAADLvzPyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd4FVXixvHvSS8ktNATSKiBEGooKipYECwgCiouKjbW3t21txV17WWxYFkbYBdBURQXBKSXhNBCL6GGACFA+p3fHxf4AQZImWRueT/Pw5Pcm8nJG014mblzzjGWZSEiIiKeI8DpACIiInIslbOIiIiHUTmLiIh4GJWziIiIh1E5i4iIeBiVs4iIiIdROYuIiHgYlbOIiIiHUTmLiIh4mCCnvnBMTIwVHx/v1JcXERGpVgsXLtxlWVa9shzrWDnHx8ezYMECp768iIhItTLGbCzrsbqsLSIi4mFUziIiIh5G5SwiIuJhHHvNuTRFRUVkZmaSn5/vdBTxQmFhYcTGxhIcHOx0FBGRSvGocs7MzCQqKor4+HiMMU7HES9iWRbZ2dlkZmaSkJDgdBwRkUrxqMva+fn51K1bV8Us5WaMoW7durrqIiI+waPKGVAxS4XpZ0dEfIXHlbOIiIi/UzkfZ+TIkSQlJdGhQwc6derE3LlzAXj99dc5ePCgbV8nPj6eXbt2Vfjzp02bxsUXX1zq88YYJk6ceOS5iy++mGnTplX4a1XEib6/+Ph4kpOT6dixI3379mX79u3HPJ+cnEy7du147LHHKCgoAGDDhg2Eh4fTqVOnI38KCwur9fsREalOKuejzJ49mx9//JFFixaxZMkSpkyZQlxcHGB/OZdXSUlJmY+NjY1l5MiR1fK1KmLq1KmkpaWRkpLCc889d8zz6enpzJs3j3Xr1jFixIgjH2vRogWpqalH/oSEhFRpRhERJ3nU3dpHe3riMpZv3WfrmO0aR/PkJUkn/Pi2bduIiYkhNDQUgJiYGADefPNNtm7dSp8+fYiJiWHq1KnceuutzJ8/n7y8PAYPHszTTz8NuM8Ar7vuOiZOnEhRURFff/01iYmJZGdnM3ToULKysujevTuWZR35updeeimbN28mPz+fu++++0gp1ahRg/vuu4/JkyfzyiuvsH//fu655x5iYmLo0qXLCb+Pjh07UlRUxG+//cb5559/zMd+//13HnjgAYqLi+nWrRvvvPMOoaGhxMfHc8MNN/Drr79yxx138O6779K5c2cWLlxIVlYWn376Kc8//zzp6elceeWVPPvssyfNXhZnnXUWb7755l+er1GjBu+++y5xcXHs3r27zOOJiPgKnTkfpW/fvmzevJnWrVtz22238ccffwBw11130bhxY6ZOncrUqVMB9+XvBQsWsGTJEv744w+WLFlyZJyYmBgWLVrErbfeyssvvwzA008/Ta9evVi8eDEDBgxg06ZNR47/6KOPWLhwIQsWLODNN98kOzsbgAMHDtC+fXvmzp1LSkoKN998MxMnTmTGjBlHLgefyGOPPXakQA/Lz89n+PDhfPnll6Snp1NcXMw777xz5ONhYWHMnDmTq666CoCQkBCmT5/OLbfcwsCBAxk1ahRLly7l448/PpLxRNnL4scffyQ5ObnUj0VHR5OQkMDq1asBWLt27ZFL2rfffnuZv4aIiDfy2DPnk53hVpUaNWqwcOFCZsyYwdSpU7nyyit54YUXGD58+F+O/eqrrxg9ejTFxcVs27aN5cuX06FDBwAuu+wyALp27cp3330HwPTp04+8f9FFF1G7du0jY7355pt8//33AGzevJnVq1dTt25dAgMDufzyywFYuXIlCQkJtGrVCoBhw4YxevToE34vZ555JgAzZsw48lxGRgYJCQm0bt0agOuuu45Ro0Zxzz33AHDllVceM8aAAQMASE5OJikpiUaNGgHQvHlzNm/eTN26dU+Y/WT69OlDYGAgHTp0+Ms/II529NWFw5e1RUT8gceWs1MCAwPp3bs3vXv3Jjk5mU8++eQv5bx+/Xpefvll5s+fT+3atRk+fPgx82sPXxYPDAykuLj4yPOlTfWZNm0aU6ZMYfbs2URERNC7d+8jY4WFhREYGHjSzz+ZRx99lJEjRxIU5P7ffHTZlSYyMvKYx4e/j4CAgCPvH35cXFx80uwnM3Xq1CMvGZxIbm4uGzZsoHXr1uTk5JxyTBERX3LKy9rGmI+MMTuNMUtP8HFjjHnTGLPGGLPEGHPiF0M9XEZGxpHLqACpqak0a9YMgKioKHJzcwHYt28fkZGR1KxZkx07dvDzzz+fcuyzzjqLMWPGAPDzzz+zZ88eAHJycqhduzYRERGsXLmSOXPmlPr5iYmJrF+/nrVr1wIwbty4U37Nvn37smfPHtLS0o6MsWHDBtasWQPAZ599xtlnn33KcU6krNnLa//+/dx2221ceumlx1xhEBHxF2U5c/4Y+A/w6Qk+3h9odehPD+CdQ2+9zv79+7nzzjvZu3cvQUFBtGzZ8sil4xEjRtC/f38aNWrE1KlT6dy5M0lJSTRv3pwzzjjjlGM/+eSTDB06lC5dunD22WfTtGlTAPr168e7775Lhw4daNOmDT179iz188PCwhg9ejQXXXQRMTEx9OrVi6VLS/330jEeffRRBg4ceGSM//73vwwZMuTIDWG33HJLWf/z/EVZs5dVnz59sCwLl8vFoEGDePzxxys1noiItzKnutQJYIyJB360LKt9KR97D5hmWda4Q48zgN6WZW072ZgpKSnWggULjnluxYoVtG3btszhRY6nnyERKTfLAlcJuIrBKgHLdeiPBeG1bPsyxpiFlmWllOVYO15zbgJsPupx5qHnTlrOIiLix1wlULgfCvYf+/bI+7lQlAfF+VBccOht4bGPS4577Cr+/5I9+q11+P3Dzx9VxK5idxGXJqwmPLSp9I9VMTvKubS7lEo9HTfGjABGAEcu64qIiA8oLoT922HfNsjdBgez4eBuyNtdyvt7oKA8N3oaCA6HoFAICoPAEPfbw4+DQiG8NgQEQ0DgoT9BYA69Lc9zJuD/3waFnjpaFbGjnDOBuKMexwJbSzvQsqzRwGhwX9a24WuLiEh1KMqHvRthzwbYvd79ds8G2JfpLuSDJ1iOOKQGRNSB8DoQURfqNHe/DasFoVEQWsN9TGiU+21I5FHvR7jLNyAI/GxjGzvKeQJwhzHmC9w3guWc6vVmERHxUEV5kJUBO1fAzuWQtdL9fs7mY48LjoTa8VAzFpp0hahG7j/RjSGqIUTEuEvZwbNPb3bKcjbGjAN6AzHGmEzgSSAYwLKsd4FJwIXAGuAgcH1VhRURERuVFMPOZZC5ALYsdL/dtYojr0wGhkBMG2jaE2KudZdx7QT328gYvzubrU6nLGfLsoae4uMWoPUURUQ8XXEBbJ4H66bBxlmwdTEU57k/FhEDsSmQNAgatIN6bd2XoAO1VpUTtLb2cbRlpD1O9P3t37+fv//977Ro0YKkpCTOOuusI/+NAwMD6dSpE0lJSXTs2JFXX30Vl8t9F+W0adOoWbPmkfW1zzvvvGr9fkS8kmXBjmUwexR8Phj+HQ+fXAwzX3Pf6dx1OFz+IdydBg+ugau/hD4PQ7uBUK+1itlB+i9/lKO3jAwNDWXXrl1H9g1+/fXXGTZsGBEREY5kKykpOWYpz5M5vGXkJZdcUuVfq7xuuummIxtaBAQEsG7dOlasWAFAeHj4kfWzd+7cydVXX01OTs6RHb/OPPNMfvzxxyrJJeIzLMt9Rrx8PCz/wX3TFkDdVtB5GDTvA/FnuKcJicfy3HL++SHYnm7vmA2Tof8LJ/ywtoys2i0j165dy9y5cxkzZgwBAe6LNs2bN6d58+Z/ObZ+/fqMHj2abt268dRTT51wTBHBXchbFsHy792FvHeT+w7n5r2h173Q4lyoFXeqUcSD6LL2UbRlZNVuGbls2TI6depU5rPy5s2b43K52LlzJ+DeYevwZe2RI0eWaQwRn5a3F+aOhndOhw/OgTnvum/gGjgKHlgNw751X7pWMXsdzz1zPskZblXRlpHVs2VkeRx9hUGXtUX4/7PkhR9B+rfuG7oadYJL3nC/VhyuzWJ8geeWs0O0ZeT/s3vLyKSkJNLS0nC5XEcua5/MunXrCAwMpH79+kdelxbxWyXF7teRZ70J29IgOAI6DIGu10MTr90MUE5Al7WPoi0jy6e8W0a2aNGClJQUnnzyySP/UFi9ejU//PDDX47Nysrilltu4Y477ij3P0pEfEpxAcz/AN7qAt/e6F4k5MKX4f6VMOAtFbOP0pnzUbRlZPlUZMvIDz74gPvvv5+WLVsSERFB3bp1eemllwDIy8ujU6dOFBUVERQUxDXXXMN9991X4XwiXq2kCFLHwB8vuZfIjO0GFzwHbS6EMlx5Eu9Wpi0jq4K2jJSqoJ8h8XqWBSsmwG9Pwp717lLu86j7zmtdRfJq1b1lpIiI2GHLQvjlEdg8B+q3g6u/glZ9Vcp+SOUsIuK0g7vh92dg4cfuNasveQM6DdMKXX7M4/7PW5alG4CkQpx6iUakwiwLlnwJkx9xz1nueRv0fgjCop1OJg7zqHIOCwsjOzubunXrqqClXCzLIjs7m7CwMKejiJTNvq0w8R5YPRliu8PFr7pXMRTBw8o5NjaWzMxMsrKynI4iXigsLIzY2FinY4ic3OGz5Un/cG8+0e8F6P533YEtx/Cocg4ODiYhIcHpGCIiVSN/H/x0H6R/DXE94dK3oW4Lp1OJB/KochYR8VmZC+Gb6yEn0z016sz7IaBqdn8T76dyFhGpSpYFCz6CXx6CGg3h+p+haQ+nU4mHUzmLiFSVojz48T5IGwstz4fLRkNEHadTiRdQOYuIVIXc7TBuKGxdBGc/BGf/Uzd9SZmpnEVE7LY9HcZeBXm74aqxkHiR04nEy6icRUTstPo3+Oo6CKsJN/wCjTo6nUi8kMpZRMQuaV/A+NugQZJ7XezoRk4nEi+lF0BEROww6y34/u8QfwYM/0nFLJWiM2cRkcqwLJj6HEx/EdoNhMveh6BQp1OJl1M5i4hUlGXB70/DzNeg8zC45E0tLCK2UDmLiFSEZcGvj8Hs/0DKDXDhK5oqJbbRT5KISHlZFvz2hLuYu/8dLnpVxSy20k+TiEh5TX8JZr0J3W6C/v8GbXErNlM5i4iUx+xRMHUkdLwa+r+kYpYqoXIWESmr1LEw+RFodykMeEuXsqXK6CdLRKQsVk+BH+6A5r3d06UCdT+tVB2Vs4jIqWxdDF9dCw3awRWfQVCI04nEx6mcRUROZu8mGDMEIurC376BsGinE4kf0HUZEZETKch17y5VXOhekjOqodOJxE+onEVESuMqge9GQNZKGPYN1GvjdCLxIypnEZHS/P40ZExyT5dqcY7TacTP6DVnEZHjLf0W/nzDvSxn95udTiN+SOUsInK0nSvghzshrgf00+pf4gyVs4jIYfk58MXfICQShnyiKVPiGL3mLCIC7s0sxt8GezbAdRMhupHTicSPqZxFRADmvgsrf4S+IyH+DKfTiJ/TZW0Rka2L4dfHoXV/OO12p9OIqJxFxM/l74Ovr4ca9eHSt3UDmHgEXdYWEf/2033uJTqH/wQRdZxOIwLozFlE/Fn6N5D+NfR+CJqd5nQakSNUziLin3Iy4cf7ILY79LrP6TQix1A5i4j/cblg/K3gKobL3tPezOJx9BMpIv5n7ruwfjpc8ibUae50GpG/0JmziPiXXWvcm1q07gddrnU6jUipVM4i4j9cLphwBwSFwiVvaNqUeCxd1hYR/zFvNGyaDZe+A1ENnU4jckI6cxYR/7B7vftydsvzoeNQp9OInJTKWUR8n2XBxLshIEiXs8Ur6LK2iPi+JV/C+j/golegZhOn04icks6cRcS3HdwNkx9xLzbS9Qan04iUicpZRHzbr49Dfg5c8joE6K888Q76SRUR37VhJqR+DqffCQ2SnE4jUmYqZxHxTSVF8NP9UKspnPUPp9OIlItuCBMR3zT3PchaCUO/gJAIp9OIlIvOnEXE9+Ruh2kvQKu+7mU6RbyMyllEfM9vT0BJAfR7QXOaxSupnEXEt2yc7Z7XfPpdULeF02lEKkTlLCK+w+WCX/4J0U3gzPucTiNSYbohTER8R9o42JYGl30AIZFOpxGpMJ05i4hvKNgPvz8DTVIgebDTaUQqRWfOIuIb/nwd9m+HKz/XTWDi9XTmLCLeb+9mmPUWtB8Mcd2cTiNSaSpnEfF+/3vW/fa8p5xMIWIblbOIeLft6e6pUz1ugVpxTqcRsYXKWUS8229PQngt6HWv00lEbFOmcjbG9DPGZBhj1hhjHirl402NMVONMYuNMUuMMRfaH1VE5DjrpsHa3+HMB9wFLeIjTlnOxphAYBTQH2gHDDXGtDvusMeAryzL6gxcBbxtd1ARkWO4XO5lOmvGQbebnE4jYquynDl3B9ZYlrXOsqxC4Atg4HHHWED0ofdrAlvtiygiUorl37sXHDnnMQgOczqNiK3KMs+5CbD5qMeZQI/jjnkK+NUYcycQCZxnSzoRkdK4SmDq81C/HSRf4XQaEduV5cy5tNn81nGPhwIfW5YVC1wIfGaM+cvYxpgRxpgFxpgFWVlZ5U8rIgKQ/jVkr4beD0OA7msV31OWn+pM4Oj5CbH89bL1jcBXAJZlzQbCgJjjB7Isa7RlWSmWZaXUq1evYolFxL+VFLv3am6YDIkXO51GpEqUpZznA62MMQnGmBDcN3xNOO6YTcC5AMaYtrjLWafGImK/tHGwZz30fkRnzeKzTvmTbVlWMXAHMBlYgfuu7GXGmGeMMQMOHXY/cLMxJg0YBwy3LOv4S98iIpVTXAjTX4TGnaFNf6fTiFSZMm18YVnWJGDScc89cdT7y4Ez7I0mInKc1M9h7ya46FVtbiE+TdeERMQ7FBfA9Jchtju01IQQ8W3aMlJEvMPCT2DfFhg4SmfN4vN05iwinq8oD2a8Ak1Ph+a9nU4jUuV05iwinm/BR7B/Owz+UGfN4hd05iwinq3wAMx8DRLOgvheTqcRqRY6cxYRz7bwEziQBb0/czqJSLXRmbOIeK7iQpj9H2h2BjQ7zek0ItVG5SwinmvJl+47tHvd53QSkWqlchYRz+QqgT9fh4YdoOW5TqcRqVYqZxHxTCsmQvYa6HWv7tAWv6NyFhHPY1kw81Wo0wLaDXQ6jUi1UzmLiOdZ+z/Ylga97oGAQKfTiFQ7lbOIeJ4Zr0JUY+hwldNJRByhchYRz7J5HmycCaffAUEhTqcRcYTKWUQ8y4xXIbw2dLnO6SQijlE5i4jn2LEcVv0MPW6B0BpOpxFxjMpZRDzHzNcgOBK6j3A6iYijVM4i4hn2boKl30LK9RBRx+k0Io5SOYuIZ5j7nvttz1udzSHiAVTOIuK8glxY9CkkXQo1Y51OI+I4lbOIOG/x51CwD3re7nQSEY+gchYRZ7lKYM47ENcTYrs6nUbEI6icRcRZK3+CvRvhtNucTiLiMVTOIuKs2aOgVjNIvNjpJCIeQ+UsIs7JXAib57gXHdEGFyJHqJxFxDlzRkFoNHQe5nQSEY+ichYRZ+RkwrLx0OVaCIt2Oo2IR1E5i4gz5r4HWFqqU6QUKmcRqX4F+2HhJ9B2ANRu5nQaEY+jchaR6pc6Fgpy4DQtOiJSGpWziFQvlwvmvgOx3SCuu9NpRDySyllEqtfa32H3Ovf0KREplcpZRKrXvPchsr779WYRKZXKWUSqz+71sPpX6DocgkKcTiPisVTOIlJ9FnwIJgBSrnc6iYhHUzmLSPUoPAiLPoO2F0N0Y6fTiHg0lbOIVI+l30L+Xuh2s9NJRDyeyllEqp5lwfz3oV5biO/ldBoRj6dyFpGql7kAtqVB95vAGKfTiHg8lbOIVL15o927T3W4yukkIl5B5SwiVWt/FiwfDx2HQmgNp9OIeAWVs4hUrUWfQEkhdLvJ6SQiXkPlLCJVp6QYFnwEzXtDvdZOpxHxGipnEak6q36GfVu0Z7NIOamcRaTqzP8AasZB635OJxHxKipnEaka2Wth3TToeh0EBDqdRsSrqJxFpGos+gRMIHQa5nQSEa+jchYR+xUXwuIx0KY/RDdyOo2I11E5i4j9Vk6Eg7ugq3afEqkIlbOI2G/hx1CrKbQ4x+kkIl5J5Swi9speC+unQ5frIEB/xYhUhH5zRMReC/8LAUHQWTeCiVSUyllE7FNcAKlj3TeCRTV0Oo2I11I5i4h9VkyEg9nQdbjTSUS8mspZROxz+Eaw5roRTKQyVM4iYo9da2DDDN0IJmID/QaJiD2O3Ah2jdNJRLyeyllEKq8o/9CNYBdCVAOn04h4PZWziFTeyh8hb7duBBOxicpZRCpv4cdQqxk07+N0EhGfoHIWkcrZvc59I1jna3QjmIhN9JskIpWzeAyYAOh0tdNJRHyGyllEKq6k2H0jWItzoWYTp9OI+AyVs4hU3NrfIXcrdNH0KRE7qZxFpOIWfQoRMdC6v9NJRHyKyllEKmb/Tlj1C3S8CoJCnE4j4lNUziJSMWlfgKtYK4KJVAGVs4iUn2XB4s8gtjvUT3Q6jYjPUTmLSPltnge7VulGMJEqonIWkfJb/CkER0LSIKeTiPgklbOIlE9BLiz9HtoPgtAop9OI+KQylbMxpp8xJsMYs8YY89AJjrnCGLPcGLPMGDPW3pgi4jGWfQ9FB6DztU4nEfFZQac6wBgTCIwCzgcygfnGmAmWZS0/6phWwMPAGZZl7THG1K+qwCLisEWfQUxriOvudBIRn1WWM+fuwBrLstZZllUIfAEMPO6Ym4FRlmXtAbAsa6e9MUXEI+xcCZnzoMu1YIzTaUR8VlnKuQmw+ajHmYeeO1proLUx5k9jzBxjTL/SBjLGjDDGLDDGLMjKyqpYYhFxTurnEBAEHa5yOomITytLOZf2z2PruMdBQCugNzAU+MAYU+svn2RZoy3LSrEsK6VevXrlzSoiTiophiVfQasLoIZ+f0WqUlnKOROIO+pxLLC1lGN+sCyryLKs9UAG7rIWEV+x9n+wfwd0Gup0EhGfV5Zyng+0MsYkGGNCgKuACccdMx7oA2CMicF9mXudnUFFxGFpYyG8jvvMWUSq1CnL2bKsYuAOYDKwAvjKsqxlxphnjDEDDh02Gcg2xiwHpgIPWpaVXVWhRaSa5e2BlT9B8hBtciFSDU45lQrAsqxJwKTjnnviqPct4L5Df0TE1yz9DkoKodPVTicR8QtaIUxETi11LNRvB406Op1ExC+onEXk5LJWwZYF7rNmzW0WqRYqZxE5ubSxYAIh+Qqnk4j4DZWziJyYqwTSvoSW50FUA6fTiPgNlbOInNi6aZC7VXObRaqZyllETixtHITVgtb9nU4i4ldUziJSuvwcWDER2l8OwWFOpxHxKypnESndsvFQnK+5zSIOUDmLSOlSx7r3bW7S1ekkIn5H5Swif5W9FjbPgY5DNbdZxAEqZxH5q7RxYAKgo/ZtFnGCyllEjuVyQdoX0Lw3RDd2Oo2IX1I5i8ixNsyAnM3Q6W9OJxHxWypnETlW2jgIjYbEi5xOIuK3VM4i8v8KcmH5D5A0CILDnU4j4rdUziLy/5ZPgKKDmtss4jCVs4j8v9SxUKc5xPVwOomIX1M5i4jbng2wcSZ01L7NIk5TOYuIW9oXgNHcZhEPoHIWEffc5tSxkHAm1IpzOo2I31M5iwhsmg17N2pus4iHUDmLCKSNhZAa0PYSp5OICCpnESk84N4est2lEBLpdBoRQeUsIit+hML90Gmo00lE5BCVs4i/Sx0DtZpB09OdTiIih6icRfxZTiasn+7etzlAfx2IeAr9Nor4s7QvAEtzm0U8jMpZxF9Zlntuc7MzoE6C02lE5CgqZxF/lTkfdq/VJhciHkjlLOKvUsdAcAS0G+h0EhE5jspZxB8V5cHS76HtAAiNcjqNiBxH5Szij1b+BAU5mtss4qFUziL+KG0cRMdC/FlOJxGRUqicRfzNvm2w9n/u6VOa2yzikfSbKeJvlnwJlkt3aYt4MJWziD85PLc5rgfUbeF0GhE5AZWziD/Zugh2ZeisWcTDqZxF/EnqWAgKg6RBTicRkZNQOYv4i+ICSP8GEi+GsJpOpxGRk1A5i/iLjJ8hf6/mNot4AZWziL9IGwdRjaB5H6eTiMgpqJxF/MH+nbD6N+hwJQQEOp1GRE5B5SziD5Z8BVaJ7tIW8RIqZxFfd3huc5OuUK+N02lEpAxUziK+bvsS2LkMOupGMBFvoXIW8XWpYyEwBNpf7nQSESkjlbOILysuhPSvoU1/iKjjdBoRKSOVs4gvW/0rHMyGTn9zOomIlIPKWcSXpY2DyPrQ4lynk4hIOaicRXzVgV2w6hfocAUEBjmdRkTKQeUs4qvSvwZXsS5pi3ghlbOIr1o8Bhp3hgbtnE4iIuWkchbxRduWwI50nTWLeCmVs4gv0txmEa+mchbxNcWFkP4VtLlQc5tFvJTKWcTXrJ6suc0iXk7lLOJrUsdCjYbQ4hynk4hIBamcRXzJ/p2wajJ0vFJzm0W8mMpZxJekf+3et7mj9m0W8WYqZxFfYVnuuc1NukL9RKfTiEglqJxFfMW2NPe+zZ101izi7VTOIr4idSwEhmpus4gPUDmL+ILiAvfc5sSLILy202lEpJJUziK+YNUvkLdHc5tFfITKWcQXpI6FqEbQoo/TSUTEBipnEW+XuwNW/wYdr4KAQKfTiIgNVM4i3i79K81tFvExKmcRb3Z4bnNsN6jX2uk0ImITlbOIN9u6GLJW6EYwER+jchbxZqljISgMkgY5nUREbFSmcjbG9DPGZBhj1hhjHjrJcYONMZYxJsW+iCJSquIC91raiRdDeC2n04iIjU5ZzsaYQGAU0B9oBww1xrQr5bgo4C5grt0hRaQUGT9D/l7orEvaIr6mLHvKdQfWWJa1DsAY8wUwEFh+3HH/Al4EHrA1oYiUbvHnEN0EEs4+5aFFJS7W7NzPjn355OQVsS+/mNCgAGqGB1M7IoSW9WtQJzKkGkKLSFmUpZybAJuPepwJ9Dj6AGNMZyDOsqwfjTEnLGdjzAhgBEDTpk3Ln1ZE3HK2wNrf4cz7S53bbFk+xZZmAAAgAElEQVQWqZv3Mil9G/M37GHFtn0UFLtOOmSTWuF0jKvJOYkNOL9dA2qGB1dVehE5hbKUsynlOevIB40JAF4Dhp9qIMuyRgOjAVJSUqxTHC4iJ5I6FizXX+7Szsot4NPZG/hu0Ra27M0jJDCATk1rcU3PZiTH1iS2djg1w0OIDg+ioMhFTl4R2QcKWbltH+lbcli4cQ+T0rcTHGg4q1U9hp8RT6+WMRhT2l8DIlJVylLOmUDcUY9jga1HPY4C2gPTDv0CNwQmGGMGWJa1wK6gInKIywWLP4OEs6BOAgCbdx/k3T/W8vXCTIpKXPRuXY/7zm/Neac4Az78i31263qA+4w7LTOHSenbGL94C9d8OI+kxtHc2rsFF7ZvRECASlqkOpSlnOcDrYwxCcAW4CrgyFJElmXlADGHHxtjpgEPqJhFqsiG6bB3I5zzOHmFJbzzx1re/WMtWHB51ybcfGZzmterUaGhjTF0iqtFp7ha3N+3NeMXb+G96eu4Y+xiOjddzzMD2pMcW9Pmb0hEjnfKcrYsq9gYcwcwGQgEPrIsa5kx5hlggWVZE6o6pIgcZdFnEFaLPwJ78sirf7Blbx4DOjbm4QsTaVQz3LYvExoUyJXdmjKkaxzfLsrk37+sZMComVzdvSkPX9iWGqFl+be9iFSEsSxnXvpNSUmxFizQybVIueTtwXq5DXNrX8xVmZfTukENnhnYnp7N61b5l96XX8Rrv63ik1kbiKsTwetXdqJzU+0dLVJWxpiFlmWVaR0QrRAm4kV2/vkZpqSAZ7Z05cZeCUy8s1e1FDNAdFgwT16SxBcjTqO4xGLwu7MZNXUNLpfu7RSxm86cRbzE/1buoNG4vmAC2Hn1r0du4nJCTl4Rj36fzo9LtnFhckNeGdKJ8BBtVylyMjpzFvEhlmUxevpaXv30a9qaDTQ+52ZHixmgZngwbw3tzCMXJvLz0u0MeW8W23LyHM0k4ktUziIezOWyeHLCMp6btJIH683DCgqjZjfP2LfZGMOIs1rwwbUprM86wKWj/mT1jlynY4n4BJWziIcqLnHxwNdpfDp7I7f1asxZ+VMxbQd43CYX57ZtwDe3nk6JC64cPYelW3KcjiTi9VTOIh6ooLiE28cu4rvFW3igb2sejFuFKdgHXa5xOlqp2jaK5utbTiM8OJCho+cwf8NupyOJeDWVs4iHKSpxccfYxUxetoMnL2nHHee0wiz+DGrHQ7NeTsc7oYSYSL6+5TTqRYVy3UfzWLhxj9ORRLyWylnEg5S4LO7/Ko3flu/g6QFJXH9GAuxeBxtmQOdhEODZv7KNa4XzxYie1I8KZfh/5+kSt0gFefZvuogfcbksHvkunQlpW/lnv0SuOz3e/YFFn4EJ+MsmF56qfnQYY27uSXRYMNd+NE83iYlUgMpZxEO8/GsGXy7YzF3ntOTW3i3cT5YUufdtbnUBRDd2NmA5NKkVzuc39SAwwHDNh/PYnpPvdCQRr6JyFvEAY+du4u1paxnavSn3nt/6/z+QMQkO7ISU650LV0EJMZF8cn13cvOLuP7j+eTmFzkdScRrqJxFHDY1YyeP/7CUPm3q8a+BScfunbzgvxAdCy3Pcy5gJbRrHM3bw7qyakcut49dTFGJy+lIIl5B5SzioBXb9nH7mEUkNoziP1d3ISjwqF/J3eth3VToci0EeO/SmGe3rsdzg9ozfVUWT05Y5nQcEa+gPd9EHLLnQCEjPltAVFgQHw3vRuTxWzAu+gRMoMfObS6PK7s1ZUP2Qd6ZtpZ2jaIZ1rOZ05FEPJrOnEUcUFzi4vaxi9iRU8B716TQIDrsuAMK3TeCte7nVTeCncwDfdvQp009npqwjLnrsp2OI+LRVM4iDnhu0kpmrc1m5KD2dIorZTnOjElwIAu6Dq/2bFUlMMDw+lWdaVongtvGLGLLXm2UIXIiKmeRajYhbSsf/bme4afHMyQlrvSDFv4XasZBy3OrN1wVqxkezOhrUygsdnHb5wspKC5xOpKIR1I5i1SjNTv389C3S0hpVptHL2pb+kHZa2HdNOhynVffCHYiLevX4KUhHUjLzOH5SSudjiPikVTOItUkr7CE28YsJCw4kLeu7kxw4Al+/RZ96r4RrPOw6g1Yjfq1b8QNZyTw8awN/LRkm9NxRDyOylmkmjw2fimrd+7n9Ss70ahmeOkHFRdC6hho0x+iG1VvwGr2UP9EOjetxT+/XcL6XQecjiPiUVTOItXgh9QtfLsokzv7tOSs1vVOfODKHw/dCOZ9K4KVV0hQwKG53Ya7xmmBEpGjqZxFqtiWvXk8Nn4pXZvV5q5zW5384IUfQ82m0KJPtWRzWpNa4bxwWTLpW3J4fcoqp+OIeAyVs0gVKnFZ3PdlKi6XxWtXdDp2BbDjZa+F9X9AV+9eEay8+rVvxBUpsbw9bS3z1u92Oo6IR1A5i1Sh92esY+763Tw1IImmdSNOfvCCjyAgCDp7/4pg5fXEJUnE1Y7g3i9T2acNMkRUziJVZemWHF75NYP+7RsyuGvsyQ8uPAiLP4O2AyCqYfUE9CA1QoN47cpObMvJ4ymtvy2ichapCvlFJdzzZSp1IkN4blDysTtNlSb9a8jPge43V09AD9S1WW3uOKcV3y3awo9LtjodR8RRKmeRKvDCzytZs3M/Lw/pSO3IkJMfbFkw731o0B6anlY9AT3Unee0pFNcLR75Lp1tOVreU/yXylnEZtMydvLxrA3ccEYCZ7Y6ybSpwzbPhR3p0O0mONUZto8LDgzgtSs7UeyyuP+rNFwuy+lIIo5QOYvYKOdgEf/4ZgmtG9TgH/3alO2T5r0PoTWhwxVVG85LJMRE8sTF7Zi1NpuPZ21wOo6II1TOIjZ69qflZB8o5NUrOhEWXIbpULk7YPkP0PlvEBJZ9QG9xJXd4ujdph4vTc5gU/ZBp+OIVDuVs4hN/liVxdcLM7nl7Oa0b1KzbJ+06BNwFbkvacsRxhieG5RMYIDhoe+WYFm6vC3+ReUsYoP9BcU88l06LepFcuc5p1gF7LCSIvfc5hbnQN0WVRvQCzWuFc7DFyYya202X8zf7HQckWqlchaxwb9/XsnWnDxeHNyxbJezAVb+BLnboPuIqg3nxYZ2a8ppzevy3E8rdPe2+BWVs0glzV2XzWdzNnL96Ql0bVa77J84/wP3Otqt+lZdOC8XEGB44fJkilwuHv1+qS5vi99QOYtUQl5hCf/8dglN60TwwAWty/6JO5bDhhnQ7Ua/Wke7IprVjeTBCxL538qd/JCqxUnEP6icRSrhtSmr2JB9kBcuTyYiJKjsnzj/AwgM9ct1tCti+OnxdGlai6cmLiMrt8DpOCJVTuUsUkGpm/fywYx1XN2jKae3iCn7J+bnQNoX0P5yiKxbdQF9SGCA4cXBHThYUKK1t8UvqJxFKqCoxMU/v1lCg+gwHu6fWL5PXvQZFB2AHroRrDxa1o/i7vNa8VP6Nn5bvsPpOCJVSuUsUgEfzlxPxo5cnhnYnqiw4LJ/YkkxzH0Xmp0BjTtXXUAfNeKs5iQ2jOLJH5ZyoKDY6TgiVUblLFJOm3cf5PUpq+jbrgHnt2tQvk9eORFyNkPP26omnI8LDgxg5KD2bM3J57XfVjkdR6TKqJxFysGyLJ74YSmBxvDUgKTyDzD7baidAG362x/OT3RtVoeh3Zvy31kbWLY1x+k4IlVC5SxSDj8v3c7UjCzuPb81jWuFl++TN8+HzHnQ81ZNn6qkh/olUjsimEe+X0qJdq4SH6RyFimj3Pwinp64jHaNohl+enz5B5gzyr37VKe/2Z7N39SMCOaxi9qRtnkvY+dudDqOiO1UziJl9Mqvq9iZW8BzlyUTFFjOX529m9y7T3W9DkJrVE1APzOwU2N6tYzhxV8y2Lkv3+k4IrZSOYuUwZLMvXwyewPX9GxGp7ha5R9g7nuAgR5/tzua3zLG8K9L21NQ4uKZH5c7HUfEVipnkVMoLnHxyPfp1KsRygMXtCn/AAW5sOhTSLoUasbaH9CPJcREcnvvlvy4ZBt/rMpyOo6IbVTOIqfw6eyNLN2yjycuaUd0eeY0H7b4cyjYBz1vtz+ccEvv5jSvF8nj45eSX1TidBwRW6icRU5iW04er/yawdmt63FRcqPyD+AqgTnvQFxPiO1qf0AhNCiQZy9tz6bdB3nrf6udjiNiC5WzyEk8PWE5xS6Lfw1sjzGm/AOs/An2boTTtOhIVTq9RQyXdWnC6OnrWL0j1+k4IpWmchY5gSnLd/DLsu3cdW4rmtaNqNggc96GWs0g8WJ7w8lfPHphWyJDg3jk+3RcmvssXk7lLFKKg4XFPDlhGa3q1+DmM5tXbJBNc2HTbC06Uk3q1gjl4f6JzN+wh68XbnY6jkilqJxFSvHGlNVs2ZvHc5clExJUwV+Tma9CeB3ocq294eSEhnSNo1t8bZ7/eSXZ+7Xvs3gvlbPIcTK25/LhzPVckRJLt/g6FRtkxzJY9Yv7rDkk0t6AckIBAYaRg5LZn1/MCz+vdDqOSIWpnEWO4nJZPDY+nRphQTzUv23FB5r5GoTUgO432xdOyqR1gyhuPDOBrxdmMn/DbqfjiFSIylnkKN8symT+hj083D+ROpEhFRtk93pY+i2kXA/hte0NKGVy97mtaFIrnMe+X0pRicvpOCLlpnIWOWTPgUKen7SClGa1GdI1ruIDzXoTAoK06IiDIkKCeGpAEhk7cvlo5nqn44iUm8pZ5JB//7KSffnFPDuoPQEBFZjTDJC7AxaPgU5XQ3QFFi0R25zfrgHntW3A64du7hPxJipnEWDhxt18MX8zN/ZKILFhdMUHmjMKXEVwxt32hZMKe2pAOwCenrDM4SQi5aNyFr9XVOLi0e+X0rhmGHef26riA+XthfkfQdIgqFPBudFiq9jaEdx1bit+Xb6DKct3OB1HpMxUzuL3Pv5zAyu35/LEJUlEhgZVfKD570NhLvS6175wUmk39kqgVf0aPDlhGQcLi52OI1ImKmfxa1v35vHalFWck1ifC5IaVHygwoPuDS5a9YWGyfYFlEoLCQrg2Uvbs2VvHm/9b43TcUTKROUsfu2ZictxWRZPD0iq2MYWhy3+DA5mQ6/77AsntunRvC6Du8byvjbGEC+hcha/NXXlTn5Ztp07z2lFXJ0KbmwBUFwAf74BTU+DZqfZF1Bs9XD/RCJDg3hs/FIsSxtjiGdTOYtfyiss4YkJS2lRL7LiG1sctvAT2LcFej9kTzipEnVrhPJQ/0Tmrt/Nd4u2OB1H5KRUzuKXRk1dw+bdefzr0vYV39gCoCgPZrwCzc6AhLPtCyhV4sqUOLo0rcVzk1aw92Ch03FETkjlLH5nzc79vDd9LYM6N+H0FjGVG2zBR7B/O/R5BCrzmrVUi4AAw7OXJrM3r4gXJ2c4HUfkhFTO4lcsy+Lx8UsJDw7kkQsrsbEFQOEB9wYXCWdDfC97AkqVa9c4muGnxzNu3iYWb9rjdByRUqmcxa/8kLqV2euyebBfIvWiQis32Lz34UAW9HnUnnBSbe49vzUNosJ49PulFGtjDPFAKmfxGzl5RTz703I6xtXi6u5NKzdYQa77Du2W50HTHvYElGpTIzSIJy5px/Jt+/h09kan44j8hcpZ/MbLkzPYfaCQkZe2J7CiG1scNvddyNvtfq1ZvFL/9g3p3aYer/62ih378p2OI3KMMpWzMaafMSbDGLPGGPOX+SLGmPuMMcuNMUuMMb8bY5rZH1Wk4lI37+XzuRu59rR42jepWbnB8nNg1lvQuj806WpPQKl2xhieGdCeohIXT0/UxhjiWU5ZzsaYQGAU0B9oBww1xrQ77rDFQIplWR2Ab4AX7Q4qUlFFJS4e+nYJDaLCuL9v68oPOPttd0H3ebjyY4mjmtZ1b4wxKX27NsYQj1KWM+fuwBrLstZZllUIfAEMPPoAy7KmWpZ18NDDOUCsvTFFKu79GetYuT2XZwYmERUWXLnBDu6GOW9D20ugUUd7AoqjRpzVnDYNonj8h6XsL9DGGOIZylLOTYDNRz3OPPTcidwI/FzaB4wxI4wxC4wxC7KyssqeUqSCNuw6wBtTVtMvqSF9kxpWfsBZb7lvBuuts2ZfERwYwPOXJ7N9Xz4va+6zeIiylHNpd86UujCtMWYYkAK8VNrHLcsabVlWimVZKfXq1St7SpEKsCyLR8enExIYwNMDkyo/4L5t7p2n2l8GDWwYTzxGl6a1ubZnMz6ZvYHUzXudjiNSpnLOBOKOehwLbD3+IGPMecCjwADLsgrsiSdScd8u2sKfa7L5Z/9EGkSHVX7Aac+DqxjOebzyY4nHeeCCNjSICuOhb5dQpLnP4rCylPN8oJUxJsEYEwJcBUw4+gBjTGfgPdzFvNP+mCLlk72/gGd/Wk7XZrUrP6cZICvDvS1ktxuhTkLlxxOPExUWzDMDk1i5PZf3Z6xzOo74uVOWs2VZxcAdwGRgBfCVZVnLjDHPGGMGHDrsJaAG8LUxJtUYM+EEw4lUi2d/WsGBgmJeuCyZgMrOaQaY8jQER8JZD1Z+LPFYfZMa0r99Q96YspoNuw44HUf8WFBZDrIsaxIw6bjnnjjq/fNsziVSYdNXZfH94i3cdW4rWjWIqvyAm+ZAxk9wzmMQWcmNMsTjPTUgiZmrd/Ho+HQ+v7EHRhuaiAO0Qpj4lIOFxTw6Pp3m9SK5rXeLyg/ocsHkR6FGQ+h5W+XHE4/XIDqMf/ZP5M812XyzMNPpOOKnVM7iU178JYPNu/N44bIOhAUHVn7Apd/ClgVw7hMQEln58cQrXN29Kd3j6/CvH5draU9xhMpZfMbcddl8PGsDw0+Pp3tCncoPWHgQpjzpXmyk49DKjydeIyDA8O/BHSgscfHId+lYVqmzR0WqjMpZfEJeYQn/+HYJTetE8I9+bewZdPZ/YN8WuOB5CNCvir9JiInkgb5t+H3lTsanbnE6jvgZ/Y0jPuGlyRlszD7Ivy/vQERIme5zPLl9W2Hma9B2AMSfUfnxxCtdf0YCXZvV5qkJy9mpy9tSjVTO4vUWbNjNf2et55qezTitRV17Bp3ylHvBkfOfsWc88UqBAYYXB3cgv6iER75fqsvbUm1UzuLV8otKePCbJTSpFc5D/RPtGXTjLFjyJZx+lxYcEVrUq8EDfdswZcUOJqT9ZXFEkSqhchav9sqvGazfdYAXL+9AZKgNl7NLimHSgxAdC2feV/nxxCfc0CuBzk1r8eSEZWTlanViqXoqZ/FaCzfu4cOZ67m6R1NOb2nT4iALPoIdS+GCkZo6JUcEBhheGtyBg4UlPDZed29L1VM5i1fKKyzhwW/SaFQznIftupx9YBdMfRaa94Z2A091tPiZlvWjuO/81kxetoMfUnV5W6qWylm80shJy1mXdYAXB3cgKizYnkF/fRwKD0D/F0FLNkopburlvnv78R+WkrnnoNNxxIepnMXr/L5iB5/P2cTNZyZwhl2Xs9f9AWlj4Yy7oZ5N86TF5wQFBvDaFZ1wuSzu/yqNEpcub0vVUDmLV8nKLeAf3ywhsWEUD1xgU4kW5cOP90LtBO06JafUtG4ETw1IYu763Yyerq0lpWqonMVrWJbFQ98uIbegmDeu6kxokA1rZwPMfBV2r4WLX4XgcHvGFJ82uGssFyY35NXfMli6JcfpOOKDVM7iNcbM3cTvK3fyUL9E2jS0YStIgKwMmPEqJF8BLc6xZ0zxecYYRl6aTJ3IEO7+YjF5hSVORxIfo3IWr7A2az/P/rScM1vFMPz0eHsGdZXAD7dDaA331CmRcqgdGcLLQzqyNusAL/y8wuk44mNUzuLxCotd3PNFKuHBgbw8pCMBATbdST3nHcic7747u0Z9e8YUv3Jmq3rc2CuBT2ZvZGrGTqfjiA9ROYvHe23KKtK35PD8Zck0iA6zZ9DstfC/f0Hr/pA8xJ4xxS89eEEb2jSI4sGvl2j1MLGNylk82tSMnbwzbS1XdYujX/tG9gzqcsGEOyEwFC5+TXOapVLCggN5Y2gncvOLuOfLxZpeJbZQOYvH2rI3j3u/TCWxYRRPDUiyb+C578LGP6HfcxBtU+GLX0tsGM2/BrbnzzXZvPH7aqfjiA9QOYtHKix2cfuYRRSXWLwzrCthwTZNm9qx3L0dZJsLodPf7BlTBBiSEsvlXWJ563+rmb4qy+k44uVUzuKRXvh5Jamb9/Li4A4kxNi0AUVxAXw3AkKj4JI3dTlbbGWM4dlL29O6fhT3fJnKtpw8pyOJF1M5i8f5OX0bH/25nuGnx3Nhso2Xnac9DzvSYcBbUKOefeOKHBIeEsjbw7pQUFTCHWMXU1TicjqSeCmVs3iUDbsO8I9vltAxrhaPXNjWvoHXz4CZr0OXayHxQvvGFTlOi3o1eP7yDizcuIeXJmc4HUe8lMpZPEZ+UQm3jVlEQIBh1NWdCQmy6cfzwC747mao2wIueN6eMUVOYkDHxlzTsxmjp6/j12XbnY4jXkjlLB7Bsiwe/X4py7ft47UrOxJbO8KegV0uGH8rHNwNQz52rwYmUg0eu7gtyU1qcv/XaazZud/pOOJlVM7iEUZPX8e3izK5+9xWnJPYwL6B54yC1b+6l+dsmGzfuCKnEBoUyDvDuhASGMBNn8xn78FCpyOJF1E5i+N+X7GDF35ZyUXJjbj73Fb2DbxprnvaVNsB0O0m+8YVKaPY2hG8d01Xtu7N5/axi3SDmJSZylkclbE9l7vGLaZ945r2rpuduwO+uhZqxrnvzta0KXFISnwdnrssmT/XZPPMxOVOxxEvEeR0APFf2fsLuPGT+USGBvH+tSmEh9i00EhJEXxzPeTnwLBvIbyWPeOKVNDgrrGs3pHLe9PX0bpBDa45Ld7pSOLhdOYsjigsdnHr54vIyi3g/WtTaFjTpg0tAH570r0854A3oWF7+8YVqYR/9Evk3MT6PDVxOTNX73I6jng4lbNUO8uyeGx8OvM27OalIR3pGGfjme3iMe6bwLr/HTpcYd+4IpUUGGB4Y2hnWtarwW1jFrIuS3dwy4mpnKXavfPHWr5akMld57RkQMfG9g28cTZMvBsSznbfnS3iYWqEBvHBdSkEBQZww8fzyd6vLSaldCpnqVbfLMzkxV8yGNCxMfec19q+gfdshC//BrWbwRWfQGCwfWOL2CiuTgTvX9uVbTn53PDxfA4UFDsdSTyQylmqzdSMnfzz2yX0ahlj753ZeXth3FXgKoahX0J4bXvGFakiXZvV4T9XdyF9S46mWEmpVM5SLRZt2sNtny8isWGUe2EGu5bmLC6AL4fBrtVwxacQ09KecUWq2PntGjByUDLTMrL45zdLcLkspyOJB9FUKqlyy7bmMPyjedSPDuW/13cjKsymS84uF3x/C2yYAZe9D8172zOuSDUZ2r0pu3ILeOW3VUSGBvHMwCSM5uQLKmepYmt27ufaD+cRGRrEmJt6UD/KpilTlgW/PgbLvoPzn9Gd2eK17jinJbkFxYyevo6I0EAe6peoghaVs1SdjdkHGPbBXIyBMTf1sG8zC4A//u2eMtXjFjj9LvvGFalmxhge7p/IgYJi3vtjHRHBQdx9no3L2IpXUjlLldiw6wBD359DfnEJ427uSfN6Nu4GNestmPY8dBrm3gJSZxni5Ywx/Gtge/KKSnhtyioAFbSfUzmL7dbvOsDQ0XMoKC5h7E09adso2r7B53/ovpyddJl7BbAA3dMoviEgwPDS4I4YDK9NWYXLsrj3fBunG4pXUTmLrdbszOVvH8ylqMRi7M02F/O892HSA9C6P1w2GgJsWotbxEMEBhheHNyBAANv/L6aYpeLB/q20WvQfkjlLLZZkrmX6z6aR2BAAGNv7kFiQxuLefYomPwItLkIhnysRUbEZwUGGP59eQeCAg2jpq4lN7+Ypy5Jsm9dAPEKKmexxey12dz86QJqRQTz+Y09iI+JtGdgy4KZr8HvT7v3ZR78kYpZfF5AgOG5QclEhwXz3vR17Msr4qUhHQkO1Ms4/kLlLJX205Jt3PtVKk3rRPD5jT3s22HK5YLfHofZ/4H2g2HQuypm8RvGGB7qn0h0eDAvTc5gz8EiRv2tCzVC9de2P9A/w6TCLMvi/enruH3sIpKb1OTrv59mXzGXFMH4W9zF3P3v7kVGVMziZ4wx3N6nJc9flszMNbu44t3Z7NiX73QsqQYqZ6mQ4hIXT01YxshJK7gwuSFjbupB7cgQewbP2wtjr4AlX8I5j0H/f+uubPFrQ7s35YPrUtiYfYBBo/5k5fZ9TkeSKqa/8aTc9hwo5Lr/zuOT2Ru5+cwE/jO0C2HBNt05vXs9fNgX1k+HAf+Bsx7UPGYRoE+b+nx1y2mUWBaXvT2LX5ZudzqSVCGVs5RLxvZcBo76k/nr9/Di4A48elE7++4i3TgLPjgX9u+Aa8ZDl2vsGVfERyQ1rsmEO3rRukEUt3y+kNenrNKGGT5K5SxlNn7xFi4d9Sd5RSWMG9GTK1Li7BnYsmD22/Dxxe7tHm/6HRLOtGdsER/TIDqML0b05PIusbw+ZTU3fbqAvQcLnY4lNlM5yynlF5XwyPfp3PNlKu2bRDPxjl50bWbTnskF++GbG2Dyw9CmP9z8P237KHIKYcGBvDykA88MTGLG6iwuenMmqZv3Oh1LbKRylpNavSOXQW/PYuzcTdxydgvG3dzTvjuyty6G986C5ePhvKfgys8hrKY9Y4v4OGMM154Wzzf/196dB8dZ33ccf3/30Oq0ZF22bMlYvpDPQPBBYpcEc9hQJ4YWappkwgxkUqZDCem0TYA2AZoQSDuZBMhMQ4EUKAmFEIgpJkBtThOMbfB9X7Jk2bLk9VrHSnt++8cjG2GwLVu7fta739fMzu6z+8yz35+O/ezze57n97vliwBc/x/v8au3dlo3d5awcDafSVX5r+W7WfDQuxzs6OWxG6fz/asa8KViEIRkEpY/CI9eAfFeuPElmPNdOw/f+nkAAA6ZSURBVPHLmDPwuboyXr5tDpc1DOMnr2zha4++T0uox+2yzCCJqjvfsqZPn66rVq1y5b3NyTUFw9z5wnre2d7O3IZqHvjLaVSVBFKz8eAu+MPfQeO7MPEr8JUHobA8Nds2JoepKs+tbuaexRvxeIR/WTCJ6y+qtXG5M4iIrFbV6QNZ14aaMcckk8qTf9rDT1/digA/umYKX581KjX/3MmEM3HF0nvA43cuk7rwG7a3bEyKiAh/Nb2OWfXl/ONz6/in363jpbUt3HftVOrKUziXujkrbM/ZALBh3xH++cUNrGkKccmEKu67dgq1Q1P0D92yBl7+e9i3GsZfCQt+DqUjU7NtY8ynJJPK0ysauf+VLSQVbrtsPDfPqSfPZ0cy3XQ6e84WzjnuSDjGv7+2lf9e0UhFUR53Xj2Ray8cmZq95XAQ3vwJrHwUCith3o9h6vW2t2zMWdJ8OMw9L23i9U2tjK0q4t6FU5g9rtLtsnKWhbM5pUg8wVN/auShZTvo7I3xzS+M5rtXTKC0IAXjV8ejTiC/9QBEOmDGt+DSu6CgbPDbNsactmVbWrl78Sb2BsPMbajmjqsaGD+sxO2yco6FszmhRFJZvHYfP3t9G03BHi6ZUMUdVzUwsSYFcy8nE7DheXjjPji8G8ZcClf+CIZPGfy2jTGD0htL8MR7e3j4jR10R+IsmlHHrXPHM7KswO3ScoaFs/mURFL533Ut/GLpdna1dTOxZgh3XNXAJROqBr/xZAI2/QHevB/at8KwKc51y+Muty5sYzJMsDvKg0u38/SKRgAWzajjb788jhEW0mln4WyO6Y0leP7DZv7z7V3sORTm/GEl3H75eOZNHj74MbFjvbDuGeea5eBOqDwfLr0DJi60WaSMyXD7Qj08vGwHz61qAmDBtBpunjOGqbU2EFC6WDgbmoJh/mdlE8+s3Et7V5RptaXc8qWxzE9FKB/eA6ufgI+egu42GHEhzL7duW7Zk6LZqYwxZ0VTMMzjy3fz7MomuqMJZo4u56Y59VwxaRjeVE1qYwAL55wVTyRZuuUgv1mxl7e3twHONHPf+rN6vjCmYnBnYCfisO2PsPrXsGOp0109YT7M+huo/5J1XxtzjuvojfHsyiZ+vXwP+0I9jCov5BsXj2LhBSMZNiRFQ/bmOAvnHLO9tZPFa1t4dlUTrR0Rhg0JsGh6HYtmjhrcyR6qcGC9M/b1mt9CZwuUjIDPf9OZzrG0NnWNMMZkhHgiyasbW3l8+W5WNx7GIzB7XCV/8fmRzJs8nMI8G7vqTFk4ZzlVZWtrJ0vWH2DJ+v3sONiFCFwyvoqvzRrFZQ3VZz4GtirsX+uc4LXpRWe4TfHA2Lkw/SYYPw+89s9pTC7Y2dbFix/t4/cf7mNfqIfCPC/zpwznz6fWMHtcJfl+O4x1Oiycs1AknuDDxhDvbG/jjxsPsKutG4/AzPpyrp5aw/zJw6k+066nSBc0vge73oCtS5xjyuKF+ktg8jXQsACKbOACY3JVMqms3BPkhY/28fK6/XRG4uT7PcweW8ncidXMbaimptTO9j4VC+cskEwqWw50snxHO+/uaOeD3UF6Ygm8HuHiMeVcNaWGeZOHn9mEFImYs3e88w0nkJs+gGQMvAEYPQcmLewL5IrUN8wYc06LxBN8sDvI0s0HWbqllaagMwPWxJohfHFsBTPry5kxupzyojyXK808Fs7noMPdUdY0h1izN8SaphBrm0OEwjEAxlUXM2dcJXPGVTJrTDkl+acxipcqHGmC5lXO2NbNq2D/GmeqRoDh02Dspc6AIaMuBr99+zXGDIyqsrOti6WbD7Jsy0E+agoRjScBGF9dzIz6cmaOLmfKyFLqK4ty/uxvC+cMFk8kaQyG2Xagk22tXWxr7WRDyxEaD4UB8AhMGFbCBXVlTB9dzpxxlQwvHWB3dbQb2rbAwc19t01wYAN0H3Re9wag5nNQOx1qZzjd1tZdbYxJkUg8wfrmI6zYHWTlniCr9xymMxIHoMDvpaGmhMkjhjB5RCmTaoYwpqro9HY2znEWzi6LJZK0hHrYGwyzNximKdhDUzDMrvZudrZ1HftmKQKjygtpGF7CBXVDuaCujGm1pRQFTnDClSr0hiC0F4K7nWPDh/c4Q2UGdznPH+UNQNX5UD0JRl4EtRfBsKngs64mY8zZkUgq21o72djSwcaWI2xs6WBzS8exwAaoKglQX1nE2Koi6iuLGF1RxIiyAkaWFVBW6M+q+ahTHs4iMh/4BeAFHlXV+497PQA8CVwEHAIWqeqek23zXAvneCJJqCdGKBwjFI5yOBzjcHeUg529HOjopbUjQmtHL60dvbR1Rkj2+7HmeT3UDi3gvIpCJgwrOXYbV11Mgd/jTA7RE3KCt7sdulqhcz90tkLXAejsu3W1ftwdfVRBOZTXw9DRUNUA1ROdQB462gYEMcZknGRSaTocZsuBTna3d7OrrYtdbd3sbu/mUHf0E+sW+L2MKMtnRFkBw4bkU1GcR0VRHhVFASqK86gsDlBelEd5Ud45ceb46YTzKa+JEREv8EvgCqAZWCkii1V1U7/VbgYOq+o4EbkBeABYdPqlD04yqUQTSSLxJJF4gmg8STTuLEfjSee1WJJoIkFPNEl3NE53JE44mqArEiccidMVSRCOxumKxOnocQK4s6eHaG8PfuLkESdPYs49cfKIUZ4PwwuFSQUJqqsTVNXFqQokqAjEGeqLUSQRPLEwRLugPQTNoY/DuPcIaPKzGxQYAsXDoGQ41M38+HHZKCd8h46GfBtqzxhz7vB4hPMqijivouhTrx0Jx2gMdtMS6qUl1OPcjvSw73APOw920d4dPdbzeLzigI/SAj8l+T6KAz6K832U5PspDvg+fi7gI9/vJd/vId/vJeA78X2g37IbBnLB6kxgh6ruAhCRZ4CFQP9wXgjc3ff4d8DDIiJ6lvrM31ryDJ73fwmaxEsSryTx4Dw+eu9HCZCkrN9zXpJ4JIkHxUsSnxx9XvFKEp/G8eOclMXJDvsq0N13O554Ia8I/IUQKIb8MiisgPKxzhSK+WWfvC+scAK4eDjkFab+h2WMMRmqtNDPtMIypp1gfCNVpSsS51BXlEPdUQ51RY7dt3dF6eyN0xWJ0dkbJ9gdZe+hMJ2ROJ29MXpjJ9gJOomSfB/r7543yFadmYGE80igqd9yMzDrROuoalxEjgAVQHv/lUTk28C3AUaNGnWGJX/aqFIfgSEJJwg9PkS84PEi/W9eH+Lx4vF4EY8Pr9eL1+fD5/Ph8zrLeLx92+i79/rBm+ccp/UGTvLYD76AE8B5RR+HcV6Rs14WHTMxxhi3iAgl+X5K8v2Mrvz0nvfJxBJJuiNxemNOz2r/+95Ygkj8k/e9sQQeFz+7BxLOn1Xd8XvEA1kHVX0EeAScY84DeO8BqZ99Hcy+LlWbM8YYk2X8Xg9lhefOCbEDGeOxGajrt1wLtJxoHRHxAaVAMBUFGmOMMblmIOG8EhgvIvUikgfcACw+bp3FwI19j68Dlp2t483GGGNMtjllt3bfMeRbgVdxLqV6XFU3isi9wCpVXQw8BjwlIjtw9phvSGfRxhhjTDYb0PRCqroEWHLccz/o97gXuD61pRljjDG56QznFTTGGGNMulg4G2OMMRnGwtkYY4zJMBbOxhhjTIaxcDbGGGMyjIWzMcYYk2EsnI0xxpgMY+FsjDHGZBgLZ2OMMSbDWDgbY4wxGcbC2RhjjMkwFs7GGGNMhrFwNsYYYzKMhbMxxhiTYURV3XljkTag0ZU3T61KoN3tIs4Ca2f2yZW25ko7IXfaeq628zxVrRrIiq6Fc7YQkVWqOt3tOtLN2pl9cqWtudJOyJ225kI7rVvbGGOMyTAWzsYYY0yGsXAevEfcLuAssXZmn1xpa660E3KnrVnfTjvmbIwxxmQY23M2xhhjMoyFszHGGJNhLJwHSUT+VUTWicgaEXlNREa4XVO6iMi/iciWvva+ICJlbteUDiJyvYhsFJGkiGTd5RoiMl9EtorIDhH5vtv1pIuIPC4iB0Vkg9u1pJOI1InIGyKyue/v9jtu15QuIpIvIh+IyNq+tt7jdk3pYsecB0lEhqhqR9/j24BJqnqLy2WlhYhcCSxT1biIPACgqt9zuayUE5GJQBL4FfAPqrrK5ZJSRkS8wDbgCqAZWAn8tapucrWwNBCRS4Au4ElVneJ2PekiIjVAjap+KCIlwGrgmiz9nQpQpKpdIuIH3gW+o6rvu1xaytme8yAdDeY+RUDWfttR1ddUNd63+D5Q62Y96aKqm1V1q9t1pMlMYIeq7lLVKPAMsNDlmtJCVd8Ggm7XkW6qul9VP+x73AlsBka6W1V6qKOrb9Hfd8vKz1wL5xQQkR+LSBPwdeAHbtdzltwEvOJ2Eea0jQSa+i03k6Uf5LlIREYDFwIr3K0kfUTEKyJrgIPA66qalW21cB4AEfk/EdnwGbeFAKp6l6rWAU8Dt7pb7eCcqq1969wFxHHae04aSDuzlHzGc1m555FrRKQYeB64/bgevayiqglVvQCn526miGTlIQuf2wWcC1T18gGu+hvgZeCHaSwnrU7VVhG5EVgAXKbn8AkLp/E7zTbNQF2/5VqgxaVaTIr0HX99HnhaVX/vdj1ng6qGRORNYD6QdSf92Z7zIInI+H6LXwW2uFVLuonIfOB7wFdVNex2PeaMrATGi0i9iOQBNwCLXa7JDELfSVKPAZtV9Wdu15NOIlJ19CoRESkALidLP3PtbO1BEpHngfNxzu5tBG5R1X3uVpUeIrIDCACH+p56PxvPTBeRa4GHgCogBKxR1XnuVpU6InI18HPACzyuqj92uaS0EJHfAl/GmV6wFfihqj7malFpICJzgHeA9TifQwB3quoS96pKDxGZBjyB87frAZ5V1XvdrSo9LJyNMcaYDGPd2sYYY0yGsXA2xhhjMoyFszHGGJNhLJyNMcaYDGPhbIwxxmQYC2djjDEmw1g4G2OMMRnm/wGndb6kZGQulgAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# use these for some plots\n", "y = np.arange(-3.5,3.5,.05)\n", "\n", "pdf_cdf_values = np.zeros((y.shape[0],3))\n", "\n", "row_index=0\n", "for i in y:\n", " pdf_cdf_values[row_index,0] = i\n", " pdf_cdf_values[row_index,1] = norm(mean,std).pdf(i)\n", " pdf_cdf_values[row_index,2] = norm(mean,std).cdf(i)\n", " row_index+=1\n", " \n", "plt.figure(figsize=(8,8))\n", "plt.plot(pdf_cdf_values[:,0],pdf_cdf_values[:,1],label=\"Standard Normal PDF\")\n", "plt.plot(pdf_cdf_values[:,0],pdf_cdf_values[:,2],label=\"Standard Normal CDF\")\n", "plt.legend(loc='upper left')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Random Variates\n", "\n", "Let's draw some random numbers from the normal distribution. This time, instead of using mean 0 and standard deviation of 1, let's assume the distribution is centered at 10 with a standard deviation of 2. Draw 100 random variates:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[10.41686472 6.72611306 9.20348461 11.67079106 14.15377413 9.78287809\n", " 9.83408784 7.83980489 6.21018425 12.23916019 9.45461662 10.63889929\n", " 6.58849774 9.96389407 7.96884601 10.42344446 9.45020002 12.93908039\n", " 10.11263851 7.98512175 11.87855743 8.42314993 9.09156988 10.90302077\n", " 7.3771751 7.03648814 10.54189927 14.09773549 9.19448907 12.93931697\n", " 10.49027688 9.22235715 5.73973794 11.20407959 9.82559841 8.78135374\n", " 5.69799647 13.62382972 8.26666461 11.31502973 9.24674708 9.09477411\n", " 12.31188227 11.20571449 12.89555143 10.90910426 12.65193488 10.73500682\n", " 9.1105223 12.11264092 9.56369209 11.47744907 9.81428745 8.91313828\n", " 11.91558726 8.0820239 7.91563515 11.06442387 5.56267745 6.09232808\n", " 10.34458283 6.64327682 7.85285525 7.73673878 9.12125366 7.71302735\n", " 9.56898416 8.09889655 9.40693208 12.77956292 13.20303638 10.8201574\n", " 10.30760335 11.8707502 12.27990302 9.86902746 6.13998553 9.44624862\n", " 10.44409674 9.52436141 9.05273226 12.713509 10.0567424 6.28136942\n", " 8.25894703 10.41825466 12.67761413 12.86492231 7.79962787 11.20064669\n", " 10.09174688 9.51563211 8.94945377 8.95720001 7.96935112 11.64306391\n", " 10.29524047 6.46765218 8.75384597 10.67989401]\n" ] } ], "source": [ "mean = 10\n", "std = 2\n", "N = 100\n", "y = norm(mean,std).rvs(N)\n", "print(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To convince ourselves, we have decent random numbers, let's view a histogram:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9.777465558583781\n", "2.037619083453264\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAETdJREFUeJzt3X+MZWV9x/H3p6BtQVpABuTXutYSWjSCZIJaUqNScAUCtrEtxNptxawabLWxqVgTaDRpaKyathjpVrZgQ9EURUlYlA01QRNFBrrAIuJSusqwW3YVBa02du23f8zZOs7euzPec3fuuM/7lUzuOc95znm+ewKfeebce+5JVSFJasfPTLoASdLyMvglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTl40gUMctRRR9Xq1asnXYYk/dS4++67v1FVU0vpuyKDf/Xq1czMzEy6DEn6qZHka0vt66UeSWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqzIq8c1dazOrLbpnIuNuuPG8i40rj5Ixfkhpj8EtSYwx+SWqMwS9JjTH4JakxiwZ/khOTfDbJg0keSPLWrv3IJJuSbO1ejxiy/9quz9Yka8f9D5Ak/WSWMuPfDby9qn4VeDFwaZJTgMuA26vqJOD2bv3HJDkSuAJ4EXAGcMWwXxCSpOWxaPBX1Y6quqdb/g7wIHA8cCFwXdftOuDVA3Z/JbCpqp6oqm8Bm4A14yhckjSan+gaf5LVwAuBO4FjqmoHzP1yAI4esMvxwKPz1me7NknShCw5+JM8A/g48Laqemqpuw1oqyHHX5dkJsnMrl27llqWJOkntKTgT/I05kL/+qr6RNf8eJJju+3HAjsH7DoLnDhv/QRg+6Axqmp9VU1X1fTU1JIeFC9JGsFSPtUT4Brgwap6/7xNNwN7PqWzFvjUgN0/A5yT5IjuTd1zujZJ0oQsZcZ/JvA64BVJNnc/5wJXAmcn2Qqc3a2TZDrJhwGq6gngPcBd3c+7uzZJ0oQs+u2cVfV5Bl+rBzhrQP8Z4A3z1jcAG0YtUJI0Xt65K0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqzKIPYkmyATgf2FlVz+/aPgac3HU5HPh2VZ02YN9twHeAHwK7q2p6THVLkka0aPAD1wJXAR/Z01BVv7tnOcn7gCf3sf/Lq+oboxYoSRqvpTx68Y4kqwdt6x7E/jvAK8ZbliRpf+l7jf/XgcerauuQ7QXcluTuJOv2daAk65LMJJnZtWtXz7IkScP0Df6LgRv2sf3MqjodeBVwaZKXDutYVeurarqqpqempnqWJUkaZuTgT3Iw8FvAx4b1qart3etO4CbgjFHHkySNR58Z/28AX6mq2UEbkxya5LA9y8A5wJYe40mSxmDR4E9yA/AF4OQks0ku6TZdxILLPEmOS7KxWz0G+HySe4EvAbdU1afHV7okaRRL+VTPxUPa/2BA23bg3G75EeDUnvVJksbMO3clqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxS3kQi6TO6stumci42648byLj6sDkjF+SGmPwS1JjDH5JaozBL0mNMfglqTFLeRDLhiQ7k2yZ1/YXSR5Lsrn7OXfIvmuSPJTk4SSXjbNwSdJoljLjvxZYM6D9A1V1WvezceHGJAcBH2TuQeunABcnOaVPsZKk/hYN/qq6A3hihGOfATxcVY9U1Q+AjwIXjnAcSdIY9bnG/5Yk93WXgo4YsP144NF567NdmyRpgkYN/g8BzwVOA3YA7xvQJwPaatgBk6xLMpNkZteuXSOWJUlazEjBX1WPV9UPq+p/gX9g7rLOQrPAifPWTwC27+OY66tquqqmp6amRilLkrQEIwV/kmPnrf4msGVAt7uAk5I8J8nTgYuAm0cZT5I0Pot+SVuSG4CXAUclmQWuAF6W5DTmLt1sA97Y9T0O+HBVnVtVu5O8BfgMcBCwoaoe2C//CknSki0a/FV18YDma4b03Q6cO299I7DXRz0lSZPjnbuS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYsGvxJNiTZmWTLvLb3JvlKkvuS3JTk8CH7bktyf5LNSWbGWbgkaTRLmfFfC6xZ0LYJeH5VvQD4KvDOfez/8qo6raqmRytRkjROiwZ/Vd0BPLGg7baq2t2tfhE4YT/UJknaD8Zxjf/1wK1DthVwW5K7k6wbw1iSpJ4Wfdj6viR5F7AbuH5IlzOranuSo4FNSb7S/QUx6FjrgHUAq1at6lOWJGkfRp7xJ1kLnA+8tqpqUJ+q2t697gRuAs4YdryqWl9V01U1PTU1NWpZkqRFjBT8SdYA7wAuqKrvDelzaJLD9iwD5wBbBvWVJC2fpXyc8wbgC8DJSWaTXAJcBRzG3OWbzUmu7voel2Rjt+sxwOeT3At8Cbilqj69X/4VkqQlW/Qaf1VdPKD5miF9twPndsuPAKf2qk5LsvqyWyYy7rYrz5vIuJL68c5dSWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjlhT8STYk2Zlky7y2I5NsSrK1ez1iyL5ruz5bu+f0SpImaKkz/muBNQvaLgNur6qTgNu79R+T5EjgCuBFzD1o/YphvyAkSctjScFfVXcATyxovhC4rlu+Dnj1gF1fCWyqqieq6lvAJvb+BSJJWkaLPnN3H46pqh0AVbUjydED+hwPPDpvfbZr20uSdcA6gFWrVvUoSzrwTOq5yuCzlQ9E+/vN3Qxoq0Edq2p9VU1X1fTU1NR+LkuS2tUn+B9PcixA97pzQJ9Z4MR56ycA23uMKUnqqU/w3wzs+ZTOWuBTA/p8BjgnyRHdm7rndG2SpAlZ6sc5bwC+AJycZDbJJcCVwNlJtgJnd+skmU7yYYCqegJ4D3BX9/Purk2SNCFLenO3qi4esumsAX1ngDfMW98AbBipOknS2HnnriQ1xuCXpMYY/JLUGINfkhrT585dSQ2Y5F3Dk3Kg363sjF+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhrjnbsaWYt3dEoHgpFn/ElOTrJ53s9TSd62oM/Lkjw5r8/l/UuWJPUx8oy/qh4CTgNIchDwGHDTgK6fq6rzRx1HkjRe47rGfxbw71X1tTEdT5K0n4wr+C8Cbhiy7SVJ7k1ya5LnjWk8SdKIegd/kqcDFwD/MmDzPcCzq+pU4O+AT+7jOOuSzCSZ2bVrV9+yJElDjGPG/yrgnqp6fOGGqnqqqr7bLW8EnpbkqEEHqar1VTVdVdNTU1NjKEuSNMg4gv9ihlzmSfKsJOmWz+jG++YYxpQkjajX5/iTHAKcDbxxXtubAKrqauA1wJuT7Aa+D1xUVdVnTElSP72Cv6q+BzxzQdvV85avAq7qM4YkabwOuDt3J3U36YH+jE5JBw6/q0eSGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGjONh69uS3J9kc5KZAduT5G+TPJzkviSn9x1TkjS6cT2I5eVV9Y0h214FnNT9vAj4UPcqSZqA5bjUcyHwkZrzReDwJMcuw7iSpAHGMeMv4LYkBfx9Va1fsP144NF567Nd2475nZKsA9YBrFq1agxlLa9JPfJRkn5S45jxn1lVpzN3SefSJC9dsD0D9qm9GqrWV9V0VU1PTU2NoSxJ0iC9g7+qtnevO4GbgDMWdJkFTpy3fgKwve+4kqTR9Ar+JIcmOWzPMnAOsGVBt5uB3+8+3fNi4Mmq2oEkaSL6XuM/BrgpyZ5j/XNVfTrJmwCq6mpgI3Au8DDwPeAPe44pSeqhV/BX1SPAqQPar563XMClfcaRJI2Pd+5KUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMaM45m7knRAmdQztLdded6yjDPyjD/JiUk+m+TBJA8keeuAPi9L8mSSzd3P5f3KlST11WfGvxt4e1Xd0z1+8e4km6rqywv6fa6qzu8xjiRpjEae8VfVjqq6p1v+DvAgcPy4CpMk7R9jeXM3yWrghcCdAza/JMm9SW5N8rxxjCdJGl3vN3eTPAP4OPC2qnpqweZ7gGdX1XeTnAt8EjhpyHHWAesAVq1a1bcsSdIQvWb8SZ7GXOhfX1WfWLi9qp6qqu92yxuBpyU5atCxqmp9VU1X1fTU1FSfsiRJ+9DnUz0BrgEerKr3D+nzrK4fSc7oxvvmqGNKkvrrc6nnTOB1wP1JNndtfw6sAqiqq4HXAG9Oshv4PnBRVVWPMSVJPY0c/FX1eSCL9LkKuGrUMSRJ4+dXNkhSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGtP3mbtrkjyU5OEklw3Y/rNJPtZtvzPJ6j7jSZL66/PM3YOADwKvAk4BLk5yyoJulwDfqqpfBj4A/NWo40mSxqPPjP8M4OGqeqSqfgB8FLhwQZ8Lgeu65RuBs/Y8fF2SNBl9gv944NF567Nd28A+VbUbeBJ4Zo8xJUk9jfywdQY/aL1G6DPXMVkHrOtWv5vkoR61jeoo4BsTGPengedmOM/NcJ6b4fY6N+l3MfzZS+3YJ/hngRPnrZ8AbB/SZzbJwcAvAk8MOlhVrQfW96intyQzVTU9yRpWKs/NcJ6b4Tw3w03y3PS51HMXcFKS5yR5OnARcPOCPjcDa7vl1wD/WlUDZ/ySpOUx8oy/qnYneQvwGeAgYENVPZDk3cBMVd0MXAP8U5KHmZvpXzSOoiVJo+tzqYeq2ghsXNB2+bzl/wZ+u88Yy2yil5pWOM/NcJ6b4Tw3w03s3MQrL5LUFr+yQZIaY/B3khye5MYkX0nyYJKXTLqmlSLJnyR5IMmWJDck+blJ1zQpSTYk2Zlky7y2I5NsSrK1ez1ikjVOypBz897u/6n7ktyU5PBJ1jgpg87NvG1/mqSSHLVc9Rj8P/I3wKer6leAU4EHJ1zPipDkeOCPgemqej5zb+S3/Cb9tcCaBW2XAbdX1UnA7d16i65l73OzCXh+Vb0A+CrwzuUuaoW4lr3PDUlOBM4Gvr6cxRj8QJJfAF7K3KeQqKofVNW3J1vVinIw8PPdvRiHsPf9Gs2oqjvY+16U+V9Nch3w6mUtaoUYdG6q6rburn2ALzJ3v09zhvx3A3PfYfZnDLmxdX8x+Of8ErAL+Mck/5bkw0kOnXRRK0FVPQb8NXMzkh3Ak1V122SrWnGOqaodAN3r0ROuZ6V6PXDrpItYKZJcADxWVfcu99gG/5yDgdOBD1XVC4H/ot0/139Md736QuA5wHHAoUl+b7JV6adNkncBu4HrJ13LSpDkEOBdwOWL9d0fDP45s8BsVd3Zrd/I3C8CwW8A/1FVu6rqf4BPAL824ZpWmseTHAvQve6ccD0rSpK1wPnAa71z//89l7nJ1L1JtjF3CeyeJM9ajsENfqCq/hN4NMnJXdNZwJcnWNJK8nXgxUkO6b5S+yx843uh+V9Nshb41ARrWVGSrAHeAVxQVd+bdD0rRVXdX1VHV9XqqlrN3OTz9C6L9juD/0f+CLg+yX3AacBfTrieFaH7K+hG4B7gfub+m2n2bswkNwBfAE5OMpvkEuBK4OwkW5n7hMaVk6xxUoacm6uAw4BNSTYnuXqiRU7IkHMzuXr8y0uS2uKMX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktSY/wOJBpixd9tGuAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print(np.mean(y))\n", "print(np.std(y))\n", "plt.hist(y,bins=10)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Jury is still out on that distribution. But with $N=100$, it is a small sample.\n", "\n", "### Likelihood functions\n", "\n", "Suppose the values y are observed and we want to calculate the likelihood of y given a mean for each of the 100 observations in y. Even though we know the mean above should be 10, let's calculate the likelihood of y given a mean of 10.3 and a standard deviation of 2 assuming y is distributed normal. Note these are simply pdf values for each datapoint:\n", "\n", "$$\n", "\\mathcal{L}_i(y_i|\\mu,\\sigma) = \\frac{1}{\\sigma \\sqrt{2\\pi}}e^{\\frac{(y_i - \\mu)^2}{2\\sigma^2}}\n", "$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.1991309 , 0.04041036, 0.17163606, 0.15771444, 0.03116273,\n", " 0.19291365, 0.19413141, 0.09360648, 0.02465169, 0.12466458,\n", " 0.18242433, 0.19662787, 0.03565051, 0.19667422, 0.10112776,\n", " 0.19909155, 0.18225368, 0.08351919, 0.19859777, 0.10208818,\n", " 0.14608507, 0.12842548, 0.16618999, 0.1906073 , 0.06856666,\n", " 0.05268632, 0.19801745, 0.03287854, 0.17121161, 0.08350615,\n", " 0.19857044, 0.17251864, 0.01482286, 0.18009772, 0.19393779,\n", " 0.14951306, 0.01413092, 0.05013348, 0.118969 , 0.17536752,\n", " 0.17364307, 0.16635073, 0.12026658, 0.18003112, 0.0859322 ,\n", " 0.19043169, 0.09990501, 0.19480824, 0.16713676, 0.13228533,\n", " 0.18640112, 0.16773264, 0.19367471, 0.15684317, 0.14394111,\n", " 0.1078501 , 0.09800516, 0.18542057, 0.012066 , 0.02181519,\n", " 0.19942159, 0.03749533, 0.09435883, 0.08774028, 0.16766856,\n", " 0.0864111 , 0.18658214, 0.10885999, 0.18054378, 0.0924937 ,\n", " 0.0695619 , 0.19283774, 0.1994697 , 0.14653474, 0.12220104,\n", " 0.19489333, 0.02293019, 0.18210039, 0.19895409, 0.18502069,\n", " 0.16422049, 0.09630703, 0.19800114, 0.02649605, 0.11850231,\n", " 0.19912277, 0.09839977, 0.0876469 , 0.0913033 , 0.18023724,\n", " 0.1983927 , 0.18470601, 0.15880431, 0.159219 , 0.10115753,\n", " 0.15920489, 0.19947058, 0.03181088, 0.14794573, 0.19590496])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(10.3,2).pdf(y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also easily calculate the joint likelihood. This is\n", "\n", "$$\n", "\\mathcal{L}(\\mathbf{y}|\\mu) = \\prod_{i \\in N}\\mathcal{L}_i(y_i|\\mu,\\sigma)\n", "$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9.256632406264856e-95" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(10.3,2).pdf(y).prod()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For computational advantages, we usually, work with log likelihoods, noting that\n", "\n", "$$\n", "ln(\\mathcal{L}(\\mathbf{y}|\\mu)) = \\sum_{i \\in N}ln(\\mathcal{L}_i(y_i|\\mu,\\sigma))\n", "$$\n", "\n", "Implemented in code for our vector $\\mathbf{y}$:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-216.52024352295027" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(10.3,2).logpdf(y).sum()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's visualize why the log-likelihood is better to work with than the likelihood:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "mu_candidate = np.arange(7,13,.2) # plot likelihood and log-likelihood in this range" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAE/CAYAAACEto0QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl8XOV1//HP0b5YlmxLsmzL+4ptbIMNmC3BrCZha9lDAwnNj9KEJmmTlqRp0jZL2zRpktKkSUhCSBoIkBXCbggBg9kM2MbG+y4vkrxII1kaaUZ6fn/MHXswkjWSZuaOZr7v12teaO6dmXvGaK7OPM+5zzHnHCIiIiIyeDl+ByAiIiKSKZRYiYiIiCSIEisRERGRBFFiJSIiIpIgSqxEREREEkSJlYiIiEiCKLEaIszsXDPbGHN/h5ldOIDX+Rcz+4X38wQzazWzXO/+n8zsY4mLutcYPmJmLyb7OINhZuvM7Dy/4xDJVmZ2k5k9HXPfmdm0AbzOvWb2Ve/nhJxHBxDD0fNuuvL+FkzxO45MoMQqzfT2QXfOLXfOzUzksZxzu5xzw5xzXYl83cEws0neCbQ15rY6ycc8euKNcs7Ncc79KZnHFUlnKUw6ekyYnHP3OecuTuSxknEeHSwzO8/Muo875/0hycd8z5do72/BtmQeN1vk+R2ASC8qnHNhv4MQEUmBvc65Wr+DkMTQiNUQ4X2rqetl3ywz225mN3j3x5rZb8ys0dv+yV6eFx0dik2wJ5rZS2bWYmZPm1llzOOv8KbImrxvPCfF7DvJ29bkPeaKmH2jzOwRMwuY2WvA1AH+G7xrOP34+L3jf+UE8Z9jZiu8GHd7U5K3ATcB/xD7TTH227qZFZrZd8xsr3f7jpkVevvOM7M6M/uMmTWY2T4z++hA3p/IUGFm/8/MtpjZIe+zPTZm38VmttHMms3sf83s+YGUGJyoZMD7LO82syXe/VlmtsyLZ6OZXdfL83o6jy4wszVevA+aWVGc7/MsM3vde97rZnZWzL7J3vtuMbNlQCUDcPxo+vHxe+epz54g/ivNbJV37t1qZkvN7GvAucB3vXPed73HHh05NLNyM/u59zdkp5n9k5nlePs+YmYvmtk3zeywRf7GXDqQ95eplFgNcWZ2KvA08DfOuQe8X/4/AKuBccAFwKfN7JI4X/JDwEeBaqAA+Kx3nBnAL4FPA1XA48AfzKzAzPK9Yz7tPe9vgPvMLDrk/j0gCIwBbvVuydJb/BOAJ4D/8eJfAKxyzt0N3Af8pzcUfnkPr/kFYLH3nPnA6cA/xeyvAcqJ/Hv/JfA9MxuR+Lcm4j8zOx/4d+A6Ip/pncAD3r5K4NfA54FRwEbgrJ5facDHv4TIuehq59xzZlYKLAPuJ/K5vxH4XzObE+dLXgcsBSYD84CPeMc50fscCTwG3EXkfX4LeMzMRnmveT/wBpGE6ivALQN/xwOO/3Tg58DfAxXA+4AdzrkvAMuBO7xz3h09vOb/EDmnTQHeD9xM5LwadQaR/7eVwH8CPzEzS/g7G6J8TazM7B7vW/7aBL3e181srXe7Pmb7T8xstZfV/9rMhiXieGngXOAR4Bbn3KPettOAKufcl51znd6c+Y+AG+J8zZ865zY559qBh4gkEwDXA48555Y550LAN4FiIifNxcAw4D+8Y/4ReBS40SKF8VcDX3LOHXHOrQV+FkccBywystRkZp+NM/YTxX8T8Ixz7pfOuZBz7qBzblWcr3kT8GXnXINzrhH4V+DDMftD3v6Qc+5xoBVIqzoOkQS6CbjHOfemc66DSBJ1pplNAj4ArHPO/dabyr8L2J/AY18L3A18wDn3mrftMiIJw0+dc2Hn3JvAb4Br4nzNu5xze51zh4h8QYw9Z/T2Pj8IbHbO/Z93zF8CG4DLvS9xpwFfdM51OOde8F73RMbGnO+aehtx62f8f+nFv8w51+2c2+Oc29DXi3nn7OuBzzvnWpxzO4D/4t3nvJ3OuR959bk/I5J4ju5HzBnN7xGre4lk2oNmZh8ETiXyS3UG8PdmNtzb/bfOufnOuXnALqCnDH0ouh1Y4Zx7LmbbRI77kAL/SPy/9LEnwTYiCRPAWCLf2ABwznUDu4mM0owFdnvbonZ6+6qI1PLtPm5fXyqdcxXe7Ztxxn6i+McDW/vxOrHe9d69n8fG3D94XD1Y7HFFMs3x54JW4CAx54KYfQ6InbpaZ8cKtM8dwLE/DTzknHs7ZttE4Izjznk3ERlJjke857zj3+fx57GdMfsOO+eOHLfvRPbGnO8qnHMPxRn7ieIf6Dmvksho//HnvHE9HdM51+b9qHOex9fEysvkD8VuM7OpZvakmb1hZsvNbFacLzcbeN779nCEyFTYUu84Ae+1jcgoi0vYm/DX7cAEM/t2zLbdwPbjPqRlzrkPDPJYe4mcwICj/5bjgT3evvHROXjPBG9fIxD2Hhu7byCOACUx9+M9cULk36W32q6+fh/e9d6JxL+3H8cWySTHnwtKiUyH7QH2AbUx+yz2vne17TDvtnwAx74WuMrMPh2zbTeRc3/sOW+Yc+6vB/D6sU70Po8/J8Cxc94+YIT3+Nh9A+HHOe8AkVH44895e/px7Kzm94hVT+4mUi+0kEh9zP/G+bzVwKVmVuLN8y8h5o+5mf2USJY9i8j8cTrLN7OimFtvV2+2EEke32dm/+Ftew0ImNmdZlZsZrlmNtfMThtkTA8BHzSzC7yaqs8AHcAK4FUiJ4B/MLN8i6z/dDnwgDdU/FvgX7z/N7MZeL3BKiLvdYKZlRMZmo/XfcCFZnadmeVZpKA+OmReT6SWoDe/BP7JzKq8360vAWm9Jo1IgvR0Lrof+KiZLbDIRRz/BrzqTRk9BpxsZld5j/0E8SUDBccdJ7eXx+0lUjf6STP7uLftUWCGmX3YO//km9lpFnNxzQCd6H0+7h3zQ9755HoiX+4fdc7tBFYC/+rVoJ5D5Hw4EKuAD5jZSDOrITJiF6+fePFfYGY5ZjYuZqCi13Oed85+CPiamZWZ2UTg79A5L25plVh5tU9nAb8ys1XAD4nM3WJmfx5TPxV7ewrAOfc0kV/2FUT+EL5MZKQEb/9HiQzRricyf5zOHgfaY27/0tsDnXNNwEVEksqveB+Ky4lMiW4n8u3jx0QKEQfMObcR+AsiSekB7xiXezVVncAVwKXevv8Fbo6Zz7+DyDDxfiLTvz8dYAzLgAeBNUQKQx898TPe9dxdROo/PkNklHQVkUJ0iJyAZnvTCL/v4elfJXKiXAO8DbzpbRPJdO85FznnngW+SKSOaR+RUZEbAJxzB4iMKv0nkWmz2UQ+Ox19HGfdccfp9cpa77N8AXCnmX3MOdcCXOzFsJfIeebrQGH/3+67jnOi93mQSG3XZ4i8z38ALvPeP0QuojmDyLnmn4kUkQ/E/xEZNNhB5OKgB/sR/2tE/h2/DTQDz3NsFOq/gWssclXfXT08/W+IfFneBrxIJMm8Z2BvIftYZArcxwAihYCPOufmejVRG51zYxLwuvcDv/CKiWO3vx/4e+fcZYM9hoiI9M4rD6gDbjquFlQkY6XViJVXC7XdzK6FyPy8mc3v42l4j80171JXM5tH5LLTp73XiK7NYURGWvq8MkJERPrPzC4xswpv+uwfAQNe8TkskZTxdeV1M/slcB5QaZFFz/6ZyNUc3zezfwLyiawbEk9Lk3xgeSR3IgD8hXMu7H1j+pk3Gmbeaw22qFFERHp2JpGpowLgHeAqb/kTkazg+1SgiIiISKZIq6lAERERkaFMiZWIiIhIgvhWY1VZWekmTZrk1+FFxAdvvPHGAedcld9xJILOYSLZJd7zl2+J1aRJk1i5cqVfhxcRH5hZPO2MhgSdw0SyS7znL00FioiIiCSIEisRERGRBFFiJSIiIpIgSqxEREREEkSJlYiIiEiCKLESERERSRAlViIiIiIJosRKREREJEGUWImIiIgkiG8rr0t22H2ojY5wN9Oqh/kdiohIn7q6HQdbO9gfCLK/OUhLMMykyhKmVZdRXpzvd3gyBCixkqT60sNraWoP8buPn+13KCIiR3V3O17dfohl79Szt6md/YEg9YEgDS0ddHW7Hp8zenghM0aXMb26jOmjhzFj9DBmjC6jrEgJlxyjxEqSak9TO22dXX6HISICREbRf/NmHb95s47dh9opys+hdkQJNcOLmDq1kpryQmqGFzF6eBE15UUMK8xj+4EjbKpvZXNDC5vrW7n/tZ0EQ90A5OUYl80bw63nTGZebYXP707SgRIrSar6QAfdvXz7ExFJhbbOME+u3c+vVtbx8raDmMHZUyv5zEUzuWRODcUFuSd8/pSqYVxw0uij97u7HXWH29nc0MLyzQf41crd/H7VXk6bNIJbz57MxXNqyM2xZL8tSVNKrCRpgqEumttDAIS6usnP1bUSIpI6+5rb+e9nNvPomn20doSZOKqEz1w0gz9fWMu4iuIBv25OjjFhVAkTRpVwwUmj+buLZ/DQ67u5d8UO/vq+N6kdUcxHzprEdaeNZ7imCbOOEitJmsaWjqM/B9pDjBpW6GM0IpItursdD7y+m39/fD2h7m4unzeWaxbWcvrkkZglfiRpeFE+Hzt3Ch89ezLL3qnnnpe289XH1vPtZZu48fQJ/O1FMygt1J/bbKH/05I09YHg0Z+blFiJSArsOHCEz/12Da9sO8TZ00bx7382jwmjSlJy7NwcY+ncGpbOrWHtnmbueXE797y0nWfW1/OdG05hwXjVYGUDzc1I0jTEjFg1tYV8jEREMl24q5u7X9jKJd95gXV7A3z96pP5xV+ekbKk6nhzx5XzresX8MBtZxLqclz9/RX8z7Obe73iUDKHRqwkaWJHrJrbO32MREQy2Yb9Ae789RpW1zVz0ezRfPWquYweXuR3WACcPnkkj3/qXL74+7X817JNvLC5kW9fv4DaEf4kfJJ8GrGSpKkPaMRKRJIn3NXNt5Zt4rK7XqTucDvf/dAp3P3hhWmTVEWVF+dz142n8J3rF7B+XwuXfmc5D6/a43dYkiQasZKkaWgJMqwwj9aOsBIrEUmoYKiLTz3wFk+tq+eqBWP50uVzGFla4HdYJ3TVKeNYOHEEn35wFZ96YBXPbWjgy1fN1ZWDGUYjVpI0DYEOplaVYhYpXhcRSYSWYIiP/vR1nlpXzz9fPpvv3HBK2idVUeNHlvDgbYv5u4tm8Ic1+7jsrhfZfajN77AkgZRYSdI0tAQZU17M8KJ8mttUYyXJZ2bfMLMNZrbGzH5nZhXe9ovM7A0ze9v77/kxz1nobd9iZndZMq7Hl4Q52NrBh370Kq/vOMR3rl/AR8+e7HdI/ZaXm8MnL5jOQ3+1mOb2ENf/8GV2Hjzid1iSIEqsJGnqAx1UDy+koiRfI1aSKsuAuc65ecAm4PPe9gPA5c65k4FbgP+Lec73gduA6d5taerClf7Y09TOtT94mc0NLfzo5kVcdco4v0MalIUTR3L//zuD9lAX1//wFbYfUHKVCZRYSVJEV10fPbyIiuJ81VhJSjjnnnbOhb27rwC13va3nHN7ve3rgCIzKzSzMcBw59zLzjkH/By4KuWBS5+2NLRwzfdXcKC1g1/85RksmVXtd0gJMWdsOb+8bTGhrm6u/+HLbGlo9TskGSQlVpIU0VXXq8sKKS8p0IiV+OFW4Iketl8NvOWc6wDGAXUx++q8bZJGVu1u4tofvEy42/HgX53Jokkj/Q4poWbVDOeB2xbT7eCGu19hc32L3yHJICixkqSIrmFV7Y1YqcZKEsXMnjGztT3crox5zBeAMHDfcc+dA3wd+Kvoph4O0esKjmZ2m5mtNLOVjY2Ng38z0qcXNx/gQz96hbKifH59+5mcNGa43yElxfTRZTxw22JyLJJcbdyv5GqoUmIlSRFdw2q0aqwkwZxzFzrn5vZwexjAzG4BLgNu8qb38LbXAr8DbnbObfU21+FNF3pqgb30wjl3t3NukXNuUVVVVaLfmhznlW0HufXe15kwsoRf334mE0eV+h1SUk2rHsYDty0mL9e44e6XeWdvwO+QZACUWElSNLR4I1Zl3ohVe4hutXKQJDOzpcCdwBXOubaY7RXAY8DnnXMvRbc75/YBLWa22Lsa8Gbg4RSHLT3YefAIt//iDcaPLObB286kOs0W/UyWKVXDePC2MynOz+VDP36FtXua/Q5J+kmJlSRFfaCD/FxjREk+5SUFOActwXDfTxQZnO8CZcAyM1tlZj/wtt8BTAO+6G1fZWbR6ue/Bn4MbAG20nNdlqRQIBji1ntfB+Ant5xGeUl2LaA5qbKUB//qTEoL8viLn7zKroNa52ooUWIlSdEQCFJdVoSZUVEcOSk2qV+gJJlzbppzbrxzboF3u93b/lXnXGnM9gXOuQZv30pvKnGqc+6O2OlDSb1wVzd33P8WOw+28YO/WMikysye/uvN+JEl3PexM+judtz2fytp69QX06FCiZUkRUNLZA0rgArv26aWXBCRvnz1sfW8sKmRr/3ZXBZPGeV3OL6aVFnKXTeewsb6Fj73m7dRzj80KLGSpKgPBBldFqmJOJpYqYBdRE7g/17Zyb0rdvCxcyZz/WkT/A4nLZw3s5rPXjyTR1bv5Scvbvc7HIlDn4mVmY03s+fMbL2ZrTOzT/XwmPPMrDmmduFLyQlXhor6QJDR3ohVeXGkh1eTllwQkV68uPkA//LIOs6fVc3nP3CS3+GklY+fN5VL59bwb4+vZ8WWA36HI32IZ8QqDHzGOXcSsBj4hJnN7uFxy2NqF76c0ChlSAmGuggEw0ev4omOWDVrxEpEerC1sZWP3/cG06qG8d83LCA3R+0aY5kZ37h2PlOrhvGJ+9+k7rCK2dNZn4mVc26fc+5N7+cWYD1amVhOoCFwbNV1gPJi1ViJSM+a2jr52M9Wkp+bw49vWURZUXZdARivYYV53H3zIsLdjtt/8QbBUJffIUkv+lVjZWaTgFOAV3vYfaaZrTazJ7zVjSVL1XtrWI32Rqzyc3MYVpinxEpE3iXc1c3H73uTPYfb+eGHFzJ+ZInfIaW1yZWl/PcNC1i3N8A//k7F7Okq7sTKzIYBvwE+7Zw7fjnYN4GJzrn5wP8Av+/lNdQOIgscHbHyaqwgMmql5RZEJNYPX9jGiq0H+dqfzc24/n/Jcv6s0fzthTP47Zt7+NmKHX6HIz2IK7Eys3wiSdV9zrnfHr/fORdwzrV6Pz8O5JtZZQ+PUzuILBDtExi9KhAidVbNGrESEc87ewN855lNfHDeGK5dNN7vcIaUO5ZM46LZo/nKY+t5ZdtBv8OR48RzVaABPwHWO+e+1ctjarzHYWane6+r/9tZqr4lSEFuztGidUD9AkXkqI5wF3/30CrKiwv46pVz/Q5nyMnJMb513Xwmjirhjvvf4vARzQakk3hGrM4GPgycH7OcwgfM7HYzu917zDXAWjNbDdwF3KDVi7NXY6CDqrJCvFwbgIriAi23ICIA3PXsZjbsb+HrV5/MiNICv8MZksqK8vnujafS1NbJVx59x+9wJEZeXw9wzr0InPDaV+fcd4n06BKhvuXYGlZR5SX5Wm5BRHhz12G+/6etXLeolgtOGu13OEPa7LHD+fh5U7nrj1u4fP5Ylsyq7vtJknRaeV0Srj7QcfSKwKiK4nya2kK6ikUki7V3dvHZh1YzpryYL17W03KI0l+fOH8aM0YP4x9/9zYtQX15TQdKrCThIg2Y3z1iVVGST7jbcaRTa6+IZKuvP7mBbQeO8I1r52m9qgQpzMvlP6+ZT30gyL8/scHvcAQlVpJg7Z3vXnU9qkJtbUSy2ootB7h3xQ4+ctYkzpr6novGZRAWjK/gY+dO4f5Xd6nlTRpQYiUJ1eAtDnr8iFV5iVZfF8lWLcEQf//rNUypLOXOpbP8Dicj/e2FM5g0qoQ7f7uGts6w3+FkNSVWklANLZHFQXuqsQL1CxTJRl959B32NbfzzevmU1yQ63c4Gam4IJevXz2P3Yfa+eZTm/wOJ6spsZKEOro46PGJVUl0KlCJlUg2eXZ9PQ+trOOvz5vKqRNG+B1ORjtjyig+vHgiP12xnTd2HvI7nKylxEoSqv64BsxR0cVC1dZGJHsEgiE+99u3mVVTxicvmO53OFnhzktnMba8mH/49Ro1avaJEitJqIYeVl2HSK9A0IiVSDb5n2c3c6C1g29cM5/CPE0BpsKwwjz+/c9PZmvjEe56drPf4WQlJVaSUA2BDqqHv3vVdYCi/FyK8nNUYyWSJbYfOMK9K3Zw3cLxnFxb7nc4WeV9M6q4dmEtP3xhG2v3NPsdTtZRYiUJVd/DGlZRamsjkj2+9tg7FObl8tlLZvodSlb6pw/OZmRpAf/w6zV0dWth5lRSYiUJ1dDy3lXXoypK8jUVKJIFXtjUyDPrG7jj/GlU9fJFS5KrvCSfL142m3f2Bfjtm3V+h5NVlFhJQtUHgr0mVuXF+TRpKlAko4W7uvnKo+8wcVQJHz17kt/hZLXL541hwfgKvvn0RtrV9SJllFhJwrR3dtESDPf6DbWiJJ9mjViJZLT7Xt3F5oZWvvCBk1Sw7jMz4wsfPIn6QAc/Wr7N73CyhhIrSZjoquu9TgUWF2i5BZEMdvhIJ99atomzp43iotmj/Q5HgNMmjWTpnBp+8PzWo+doSS4lVpIw0TWsRg/vfcRKNVYimes7z2yiJRjiS5fNec+VweKfOy+dRWe4m28v0/ILqaDEShLmWJ/AXmqsSvLpCHdr0TqRDLSpvoVfvLqLm86YyMyaMr/DkRiTK0v58JkTefD1XWzc3+J3OBlPiZUkTJ8jVsWRtjZay0okszjn+Mqj71BakMvfXjTD73CkB588fzqlhXn8+xPr/Q4l4ymxkoRpCAQpyMs5usr68Y62tdF0oEhG+eOGBpZvPsCnL5zByNICv8ORHowoLeBvzp/GnzY2snxzo9/hZDQlVpIwDS0dVJe9d9X1qIqjbW1UwC6SKTrD3Xz1sfVMrYpMN0n6uuWsSdSOKOZrj63XoqFJpMRKEuZEa1hBpMYK0FpWIhnkZyt2sP3AEb542Wzyc/UnJZ0V5uVy59JZbNjfwm+0aGjS6FMgCXOidjYAFSVejZWmAkUyQlNbJ3f9cTNLZlZx3sxqv8OROFzmLRr6X09vpK0z7Hc4GUmJlSTMidrZQMxUoNayEskI97y4nZZgmDsvneV3KBInM+OfvEVDf7x8u9/hZCQlVpIQbZ1hWoJhqnu5IhCgpCCX/FxT8bpIBmhuC/HTl3Zw6dwaZtUM9zsc6YdFk0Zy6Vxv0dCAFg1NNCVWkhAN3lILva1hBZFvSuoXKJIZ7nlpOy0dYT55wXS/Q5EBuHOpt2joM5v8DiXjKLGShGhoOfEaVlHlxeoXKMljZt8wsw1mtsbMfmdmFcftn2BmrWb22ZhtS81so5ltMbPPpT7qoae5PcQ9L21n6ZwaThqj0aqhaFJlKTedMYFfrayj7nCb3+FkFCVWkhD1gRP3CYyqKFG/QEmqZcBc59w8YBPw+eP2fxt4InrHzHKB7wGXArOBG81sdopiHbJ++lKktkqjVUPb7edNxQx++LwaNCeSEitJiGhidaKrAiFSwK4aK0kW59zTzrnopU6vALXRfWZ2FbANWBfzlNOBLc65bc65TuAB4MpUxTsUNbeH+MmL27lkzmhmj9Vo1VA2pryYaxbW8uDK3aq1SiAlVpIQjS0dJ1x1PapcjZgldW7FG50ys1LgTuBfj3vMOGB3zP06b5v04t6Xdmi0KoP89fun0dXt+NFyjVolihIrSYjI4qC9r7oeVVFcoF6BMihm9oyZre3hdmXMY74AhIH7vE3/CnzbOdd6/Mv1cIhel6Q2s9vMbKWZrWxszL62IIFgiJ+8uI2LZo9mzthyv8ORBJgwqoQr5o/lF6/s4tARlWkkQp7fAUhmqA90nPCKwKiKknxaO8KEurq1SrMMiHPuwhPtN7NbgMuAC5xz0STpDOAaM/tPoALoNrMg8AYwPubptcDeExz7buBugEWLFmVdT5B7X9pBIBjmUxqtyigfP28qv1+1h5++tJ3PXDzT73CGPP1lk4RoaAn2eUUgHGvErFErSQYzW0pkyu8K59zRS52cc+c65yY55yYB3wH+zTn3XeB1YLqZTTazAuAG4BEfQk97kdGq7Vx40mjmjtNoVSaZPrqMpXNquHfFDgJBnZsHS4mVJERDnCNW5UcbMevDK0nxXaAMWGZmq8zsByd6sFfofgfwFLAeeMg5t+5Ez8lWP3tpB83tIY1WZahPLJlGSzDM/7280+9QhjxNBcqgtXWGaek48arrUUf7BWrJBUkC59y0OB7zL8fdfxx4PFkxZYKWYIgfv7idC0+q5uRajVZlornjylkys4ofL9/GR8+eREmB0oOB0oiVDFp01fXR8dRYacRKZMj5+cs7vdGqGX6HIkl0x/nTOdwW4v5Xd/kdypCmxEoG7egaVv2osVJiJTI0tHaE+dHybZw/S6NVmW7hxBGcNXUUd7+wjWCoy+9whiwlVjJo9Ufb2cQzYhWZClS/QJGh4WcrdtDUptqqbHHHkmk0tHTw6zfq/A5lyOozsTKz8Wb2nJmtN7N1ZvapHh5jZnaX12trjZmdmpxwJR1FV+yNZyqwrCgPM2huU42VSLpr7+zix8u3sWRmFfPHV/T9BBnyzpw6ilMnVPD9P20l1NXtdzhDUjwjVmHgM865k4DFwCd66KV1KTDdu90GfD+hUUpaa/BWXR9e3HexY06OUV6crxErkSHgd2/t4XBbiNvfP9XvUCRFzIw7zp/GnqZ2fv/WHr/DGZL6TKycc/ucc296P7cQuST5+JYPVwI/dxGvABVmNibh0UpainfV9Sj1CxRJf8457nlpO3PHDef0ySP9DkdSaMnMamaPGc73/7SVru6sWwd30PpVY2Vmk4BTgFeP2xVXv61sbweRqRoCHXFNA0aVlxRoxEokzS3ffIAtDa3cevbkuL80SWaIjlptO3CEx9/e53c4Q07ciZWZDQN+A3zaORc4fncPT3lPmuucu9s5t8g5t6iqqqp/kUraqm8JxnVFYFRFcb5qrETS3D0vbaeqrJAPztPkQzZaOqeGadXD+MHzWznWGUriEVcHOHZvAAAgAElEQVRiZWb5RJKq+5xzv+3hIXX0o9+WZJbGOFddj6ooUY2VSDrb0tDKnzY28uHFEynMy/U7HPFBTo5x69mTWbc3wMqdh/0OZ0iJ56pAA34CrHfOfauXhz0C3OxdHbgYaHbOafwwCxzpiKy6Hs9SC1GqsRJJb/eu2E5BXg4fOmOC36GIj646ZSzDi/K4d8UOv0MZUuJZs/5s4MPA22a2ytv2j8AEAOfcD4i0g/gAsAVoAz6a+FAlHTV4a1hVl8U/FVheUkAgGKKr25Gbo9oNkXTS1NbJb97Yw1ULxlI5LP7PtWSekoI8bjh9Aj95cTv7m4PUlMf/BTqb9ZlYOedepOcaqtjHOOATiQpKho6ja1j1c8TKuUj/sWjvQBFJD798bTftoS5uPWey36FIGvjw4on8aPk27nt1J5+5eKbf4QwJWnldBuXYquv9KF5XWxuRtBTq6ubnL+/g7GmjmFUz3O9wJA2MH1nCBbNGc/+ru9TmJk5KrGRQoiNW/S1eB7W1EUk3T67dz77mILeerdEqOeYjZ03i4JFOHluj0ul4KLGSQWlo6aAwzlXXo8qj/QK15IJIWrnnpe1MGlXCkpnVfociaeTsaaOYVj2Mn728Q0svxEGJlQxKQyCyhlV/FhCMjlg1a8RKJG28ueswb+1q4qNnTyZHF5VIDDPjlrMmsaaumbd2N/kdTtpTYiWD0tQeYkQ/C9ArilVjJZJufvrSDsqK8rhmYa3foUga+vNTxlFWlMe9L+3wO5S0p8RKBqW5PcTwovx+PadciZVIWtnX3M7jb+/jxtMnUFoY/7S+ZI/SwjyuXTiex9/ed7S2VnqmxEoGJdAeOpooxSsvN4eywjya2lVjJZIOfv7yTpxz3HzmRL9DkTR285kT6XKO+17d5XcoaU2JlQxKIBjuV+F6VHlJPs0asRLxXXtnF/e/uoulc2uoHVHidziSxiZVlnLejCrue3UXneFuv8NJW0qsZFACA5gKBPULFEkXv3mzjub2kJZYkLh85OzJHGjt4PG3tfRCb5RYyYAFQ110hLsZ3s+pQICK4gIttyDiM+cc967YwbzachZOHOF3ODIEnDutkimVpeofeAJKrGTAAsHIiNPwooFNBWrESsRfK3ceZktDK3+xeGK/lkyR7JWTY9x85kRW7W5itZZe6JESKxmwQHsYYIAjVqqxEvHbA6/tpqwwj8vmjfE7FBlCrl5YS2lBLj/TqFWPlFjJgB0dsRpIYuWNWGkVXxF/BIIhHnt7L1csGEtJgZZYkPiVFeVz7aLx/GHNXhq9frFyjBIrGbBAe3QqcGA1Vl3djtaOcKLDEpE4PLxqL8FQNzecNsHvUGQIuvnMiYS6HL98TUsvHE+JlQxYIBhJisoHuNwCaJFQEb88+PouZo8Zztxxw/0ORYagKVXDOHd6JQ+8touubs08xFJiJQM2uBEr9QsU8cvaPc2s3RPghtPHq2hdBuz608aztznIiq0H/A4lrSixkgGLJkUDq7GK9BfUiJVI6j34+m4K83K4cv44v0ORIeyi2aOpKMnnoZV1foeSVpRYyYAFgiEK8nIoys/t93MrolOBamsjklLtnV38ftUePnjymKNT8iIDUZiXy5Xzx/LUuv26yjuGEisZsEB7eEDTgHBsKlAjViKp9fjb+2gJhrn+tPF+hyIZ4NpF4+kMd/PI6j1+h5I2lFjJgAWCoQH1CYRj04eqsRJJrQdf383kylJOnzzS71AkA8wdV87sMcM1HRhDiZUM2ED7BAIU5edSnJ+rtjYiKbS1sZXXdhzi+tNUtC6Jc92iWt7e08w7ewN+h5IWlFjJgAWC4QEVrkdVlORrKlAkhR56fTd5OcbVp9b6HYpkkCsXjKMgN4dfvbHb71DSghIrGbDIiNXAV2wuL1a/QEksM/uGmW0wszVm9jszq4jZN8/MXjazdWb2tpkVedsXeve3mNldlqFDOZ3hbn7zZh0XnjSaqrJCv8ORDDKitICL5ozm92/toSPc5Xc4vlNiJQMWaA9RPsgRK11JIgm2DJjrnJsHbAI+D2BmecAvgNudc3OA84DoL9/3gduA6d5taYpjToln19dzoLWT609X0bok3rULazncFuLZ9Q1+h+I7JVYyIM45r3h9EIlVcYGWW5CEcs497ZyL9kl6BYjOeV0MrHHOrfYed9A512VmY4DhzrmXXaRx5c+Bq1IeeAo88PpuxpQX8b7pVX6HIhno3OlV1Awv4qGVmg5UYiUDEgx1E+pyAy5eB9VYSdLdCjzh/TwDcGb2lJm9aWb/4G0fB8RezlTnbcsoe5raeWFzI9cuGk9uTkbOdIrPcnOMaxbW8sKmRvY3B/0Ox1dKrGRAAsHoquuDqLEqidRYRQYKROJjZs+Y2doeblfGPOYLQBi4z9uUB5wD3OT998/M7AKgpyyj119IM7vNzFaa2crGxsaEvadk+5U3inDdIhWtS/Jcs7CWbge/eTO7l15QYiUD0jyIPoFRFcUFdIa7CYa6ExWWZAHn3IXOubk93B4GMLNbgMuAm9yxrL0OeN45d8A51wY8DpzqbY/NNmqBvSc49t3OuUXOuUVVVUNjSq2r2/GrlXWcO72K2hElfocjGWxSZSlnTB7Jr1buzuovzEqsZEACg+gTGKW2NpJoZrYUuBO4wkugop4C5plZiVfI/n7gHefcPqDFzBZ7VwPeDDyc8sCTaPnmRvY0tXODVlqXFLhu0Xh2HGzj9R2H/Q7FN0qsZECiU4GDuipQbW0k8b4LlAHLzGyVmf0AwDl3GPgW8DqwCnjTOfeY95y/Bn4MbAG2cqwuKyM8+PpuRpYWcOFJo/0ORbLApSfXMKwwL6uL2AdeICNZLdAeufBqUOtYlSixksRyzk07wb5fEFly4fjtK4G5yYzLL83e5e8fPnMiBXn6Hi3JV1KQx2XzxvDwqr38yxVzGFaYfWmGPmkyIMeK1wdXYwXQrKlAkaR4Yu0+Oru6uWpBxl3oKGns2kXjaQ918fiafX6H4gslVjIg0RqrskGMWFVoxEokqX6/ag9TqkqZO26436FIFjl1QgVTq0qzdjpQiZUMSHN7iKL8HArzcgf8GseK15VYiSTavuZ2Xt1+iCvnj1PDZUkpM+O6ReNZufMwWxtb/Q4n5ZRYyYAE2sODWmoBoDg/l4LcHI1YiSTBH1bvxTm4csFYv0ORLPRnp44jN8f41crsW9Oqz8TKzO4xswYzW9vL/vPMrNm7AmeVmX0p8WFKugkEB9cnECLfaspL8lVjJZIEv39rLwvGVzCpstTvUCQLVZcV8f4ZVTyyag/d3dm1plU8I1b30ndT0uXOuQXe7cuDD0vS3WD7BEZVFKutjUiiba5v4Z19AY1Wia+umD+Wvc1B3tiVXWta9ZlYOedeAA6lIBYZQiJTgYO/jLaiJJ/DbRqxEkmkh1ftJcfgsnlKrMQ/F84eTWFeDn9Y3Wszg4yUqBqrM81stZk9YWZzEvSaksYSNWJVXlygESuRBHLO8fDqPZw9rZKqskK/w5EsNqwwjwtPGs3jb+8j3JU9rcsSkVi9CUx0zs0H/gf4fW8PHKoNTOW9mttDgy5eBxhRoqlAkUR6c1cTuw+1a+0qSQuXzx/LgdZOVmw96HcoKTPoxMo5F3DOtXo/Pw7km1llL48dcg1M5b2ccwTaQwwvHvxU4IjSAk0FiiTQw6v2UJiXw8Vz1MJG/HfezCrKCvOyajpw0ImVmdV4zUsxs9O918ye1DQLHensotsNrk9gVEVJPh3hbto7uxIQmUh2C3V189iafVw4ezRlCRhRFhmsovxcLp5Tw5Pr9tMRzo7zfDzLLfwSeBmYaWZ1ZvaXZna7md3uPeQaYK2ZrQbuAm5wzmXXtZVZJrrqemKmAiNtbTRqJTJ4L245wMEjnVw5X0Xrkj6uWDCWlmCY5zdmRwlQn3M5zrkb+9j/XSId5SVLJKJPYNQIb/X1w22djK0oHvTriWSzR1btpbw4n/NmVvsdishRZ00dxcjSAh5ZvZeL59T4HU7SaeV16bdAexhIzIhVhTdipQJ2kcFp6wzz1Lr9fODkGgrydGqX9JGfm8MHTq7hmfX1HOkI+x1O0unTJ/3WHJ0KTETxuqYCRRJi2Tv1tHV2caWuBpQ0dPm8sQRD3Tyzvt7vUJJOiZX0W7TGKhHF68emAjViJTIYj6zay5jyIk6fNNLvUETe47RJIxlTXpQVVwcqsZJ+O1pjlcipwCMasRIZqENHOnl+UyNXzB9LTo75HY7Ie+TkGJfNG8PzmxppyvAZCiVW0m/RGquyBLS0KcjLobQgVyNWIoPw+Nv7CHc7rlBvQEljl88fS6jL8dS6/X6HklRKrKTfAsEQpQW55OUm5tenoqQg47/BiCTTw6v2ML16GLPHDPc7FJFenTyunEmjSngkw6cDlVhJvzW3J6ZPYNSIUjViFhmousNtvL7jMFcuGIu3VrNIWjIzrpg/lpe3HqShJeh3OEmjxEr6LZCgPoFRI0oKNBUoMkDRb/+6GlCGgsvnj6XbweNr9vkdStIosZJ+CwRDCbkiMEpTgSID98iqvZw6oYLxI0v8DkWkT9NHlzGrpiyjpwOVWEm/BdrDCVnDKmpESb5GrEQGYMeBI2zY38IH56loXYaOKxaM5c1dTew+1OZ3KEmhxEr6LRBM7FRgRUkBgWCIrm61mBTpjyfWRq6uumTOaJ8jEYnf5d4XgUczdDpQiZX0WyDRxesl+Th3bEV3EYnPk+v2M6+2nNoRmgaUoWP8yBJOmVCRsdOBSqykX7q7HS0dYYYnYA2rKLW1Eem/vU3trN7dxNK5md/UVjLPFfPHsn5fgC0NLX6HknBKrKRfWjrCOEdCR6wqvLY2KmAXiV90kcWlc5RYydDzwZPHkGPwyOrMmw5UYiX9EjjagDmxyy0AHD6iqUCReD2xdj8zR5cxpWqY36GI9Fv18CJOmzSSJ9cqsZIsl8g+gVGaChTpn8aWDl7fcYhLNA0oQ9glc2rYVN/K9gNH/A4loZRYSb9E+wQmcrmFitLoVKBGrETiseydepyDS5VYyRB2sXc1a6b1DlRiJf2SjBGrssI88nJMI1YicXpi7T4mjSphVk2Z36GIDFjtiBLmjhuuxEqyW3RJhESuvG5mVGiRUJG4NLeFeHnrQS6ZW6PegDLkXTK7hrd2NVEfyJzegUqspF+OFq8ncMQK1NZGJF7PrK8n3O24dO4Yv0MRGbRoneDT79T7HEniKLGSfgkEw5hBWQLXsYJoWxslViJ9eWLtfsaUFzFvXLnfoYgM2vTqYUyuLOXpDJoOVGIl/RJoDzGsMI+cnMROQURGrDQVKINjZt8wsw1mtsbMfmdmFd72fDP7mZm9bWbrzezzMc9ZamYbzWyLmX3Ov+j7dqQjzAubG7lkTk3CP4MifjAzLp4zmpe3HqQ5Q/4GKLGSfkl0n8AojVhJgiwD5jrn5gGbgGgCdS1Q6Jw7GVgI/JWZTTKzXOB7wKXAbOBGM5vtQ9xxeW5jA53hbl0NKBnlkjk1hLsdf9yYGdOBSqykXxLdJzBqREkBh9tCOKdGzDJwzrmnnXNh7+4rQG10F1BqZnlAMdAJBIDTgS3OuW3OuU7gAeDKFIcdtyfW7qdyWAGLJo30OxSRhFlQW0F1WSFPrVViJVko0J7YPoFRFSUFdIa7aQ91Jfy1JWvdCjzh/fxr4AiwD9gFfNM5dwgYB+yOeU6dty3tBENdPLehgYtm15CraUDJIDk5kenA5zc1EsyAvwFKrKRfAsFkjVhFXlNLLkhfzOwZM1vbw+3KmMd8AQgD93mbTge6gLHAZOAzZjYF6ClD6XXY1MxuM7OVZraysbExYe8pHss3H6Cts0tNlyUjXTKnhvZQFy9sSu3nKhkSP/QgGS3QHkroGlZR0UbMh490Mq6iOOGvL5nDOXfhifab2S3AZcAF7tjc8oeAJ51zIaDBzF4CFhEZrRof8/RaYO8Jjn03cDfAokWLUjpv/eTa/QwvyuPMKaNSeViRlFg8ZRTDi/J4al09Fw/xxuIasZJ+CQTDSSler/D6BerKQBkMM1sK3Alc4Zxri9m1CzjfIkqBxcAG4HVguplNNrMC4AbgkVTH3ZdQVzfPrK/nwtmjKcjTaVsyT35uDhecNJpnN9QT7ur2O5xB0SdU4hbu6qa1I5zQPoFRasQsCfJdoAxYZmarzOwH3vbvAcOAtUSSqZ8659Z4he53AE8B64GHnHPrfIj7hF7eepDm9hBLh/g3eZETuWTOaJraQry2/ZDfoQyKpgIlbi1BrwFzkpZbALT6ugyKc25aL9tbiSy50NO+x4HHkxnXYD25bj8lBbm8b0aV36GIJM37ZlRRmJfDU+v2c9a0Sr/DGTCNWEncjjZgTkqNVXTESlOBIrG6uh1Pr9vPklnVFOXn+h2OSNKUFOTxvhlVPP1O/ZBeekeJlcQt0B4dsUr8QGdBXg6lBbmaChQ5zsodhzjQ2qlpQMkKl8ypYV9zkDV1zX6HMmBKrCRu0RGrZFwVCGprI9KTJ9ftpyAvhyWzqv0ORSTpLjypmtwc46kh3DtQiZXELdCevKlAgBGlamsjEss5x7J36jl3WiXDClUSK5mvoqSAMyaPVGIl2SGZNVZwrK2NiERsbmil7nA7F5w02u9QRFLmkjk1bG08wpaGVr9DGZA+Eyszu8fMGsxsbS/7zczu8jrDrzGzUxMfpqSD5uiIVRJqrCA6FagRK5GoZ9c3AHC+pgEli1w8J/JFYqiOWsUzYnUvsPQE+y8Fpnu324DvDz4sSUeB9jA5BqUFyUmsRpTkc/iIEiuRqD9uqGfO2OHUlBf5HYpIyowpL2Z+bTlPZ2pi5Zx7ATjRal1XAj93Ea8AFWY2JlEBSvqI9gnMSVID2IqSAgLB8JBfdVckEQ4f6eSNnYe5QKNVkoUunlPD6rpm9jW3+x1KvyWixmrIdIeXwQm0h5KyOGhUdJHQ6JSjSDZ7flMj3Q7OV32VZKFLvOVFnl5X73Mk/ZeIxCru7vB+doaXwQsEk9POJmqEFgkVOerZDQ1UDitg3rhyv0MRSblp1cOYUlXKM+uzM7GqI87u8M65u51zi5xzi6qq1JphqEn2iFWF2tqIAJGmy89vbGDJzOqkTb2LpLvzZ1bz6rZDtHWG/Q6lXxKRWD0C3OxdHbgYaHbO7UvA60qaaU76VKBGrEQA3th5mEAwzAUnqb5KsteSWdV0dnXz0paDfofSL33O65jZL4HzgEozqwP+GcgHcM79gEjz0g8AW4A24KPJClb8FSleT8VUoEasJLs9t6GB/FzjnOka2ZfsddqkkQwrzOO5jQ1cNHvo1Br2+VfSOXdjH/sd8ImERSRpK9AeTlo7G4CKUk0FikCkvmrxlFFabV2yWkFeDudMq+S5DQ045zAbGtPiWnld4tIZ7qY91JXUqcCywjzyckxTgZLVdh6MrDitRUFFYMmsKvY1B9mwv8XvUOKmxEri0pLkdjYAZkZFSb5GrCSr/XGDVlsXiTpvZuRz8NzGBp8jiZ8SK4nL0XY2SayxgsgioYePaMRKstcfNzQwtaqUiaNK/Q5FxHejhxcxZ+xwntugxEoyTCAYudw1mVOB4LW10YiVZKnWjjCvbDuopssiMc6fVc0bOw/TPETKRJRYSVwC7cmfCoRoI+ah8eERSbQXNzcS6nKaBhSJcd7MarodPL95aCwsrsRK4hLwaqySeVUgaMRKstuz6xsYXpTHwokj/A5FJG0sGF/ByNKCITMdqMRK4hJoT9VUYAFN7SEiq3iIZI/ubsdzGxt4/8xq8nN1ahaJys0x3j+jiuc3NdLVnf5/G/TplbgEgqkrXo8u7SCSTdbsaeZAaycXaBpQ5D3Om1nFoSOdrK5r8juUPimxkrg0t4fIyzGK83OTepwRXr9ArWUl2eaP6+vJMXj/DK22LnK898+oIsfgT0NgOlCJlcQl0B5ieHF+0le+rYi2tTmiOivJLn/c2MDCiSMYUVrgdygiaaeipIBTJ4zgj0NgPSslVhKXQDDM8KLkt9eIjljpykDJJvubg6zdE+D8WVpmQaQ3S2ZVs3ZPgIZA0O9QTkiJlcQl0B5K+hWBwNFv67oyULJJdFXpC05SfZVIb5Z4q7D/aWN6L7ugxEriEgiGkr6GFUBFiRoxS/Z5dn0DtSOKmV49zO9QRNLWSWPKGFNedLTtU7pSYiVxCbSHkr7UAkBFcXTESlOBkh2CoS5e2nKAC2ZVJ72GUWQoMzPOm1nNi1sO0Bnu9jucXimxkrg0t4eTvtQCQEFeDsMK8zQVKFnj5W0HaQ91sUTLLIj0acnMKlo7wqzcccjvUHqlxEriEgimZsQKItOBKl6XbPHH9Q0U5+eyeMoov0MRSXtnT6ukIDfnaF1iOlJiJX0KhrroDHenpMYKIquva8RKssXzmxo5e9ooipK8RpxIJigtzOOMKSPTus5KiZX06diq66kbsVKNlWSDnQePsOtQG+/ToqAicVsys5qtjUfYdbDN71B6pMRK+nSsT2Dya6zA6xeoESvJAi9silw2fu50JVYi8YrWI6brdKASK+lTqkesRpTka+V1yQovbD7A+JHFTBpV4ncoIkPG5MpSJleWpu10oBIr6VNzu5dYpax4vYBAMEy4K30vpxUZrFBXNy9vPci506u0zIJIPy2ZWR25orazy+9Q3kOJlfQp4CVW5SlYbgGOtbWJJnQi/WFmXzGzNWa2ysyeNrOx3nYzs7vMbIu3/9SY59xiZpu92y2piHPV7iZaO8K8b3plKg4nklGWzKqiM9zNiq0H/A7lPZRYSZ8CQa/GKlVTgaVaJFQG5RvOuXnOuQXAo8CXvO2XAtO9223A9wHMbCTwz8AZwOnAP5vZiGQH+cKmRnJzjDOnKrES6a/TJ4+kpCA3LdvbKLGSPgV8mAoEtbWRgXHOBWLulgLO+/lK4Ocu4hWgwszGAJcAy5xzh5xzh4FlwNJkx/nC5gMsGF+Rkh6cIpmmMC+XMyaP5MUtGrGSISgQDFGQl5OydXaiU4EasZKBMrOvmdlu4CaOjViNA3bHPKzO29bb9p5e9zYzW2lmKxsbB/5NuamtkzV1TZyraUCRATtnehXbDxyh7nB6LbugxEr6lKo+gVEjSqJTgRqxkp6Z2TNmtraH25UAzrkvOOfGA/cBd0Sf1sNLuRNsf+9G5+52zi1yzi2qqhr4EgkvbjmAc1pmQWQwol9MXtycXqNWSqykT4EU9QmMqvBGrDQVKL1xzl3onJvbw+3h4x56P3C193MdMD5mXy2w9wTbk2b5pgOUFeUxv7Y8mYcRyWjTq4cxenghy9NsOlCJlfQplX0CAYYV5pGXY5oKlAExs+kxd68ANng/PwLc7F0duBhods7tA54CLjazEV7R+sXetqRwzrF8cyPnTKskL1enYJGBMjPOmVbFii0H6O7ucZDZF/pUS58C7aGUFtiaGRVafV0G7j+8acE1RJKkT3nbHwe2AVuAHwEfB3DOHQK+Arzu3b7sbUuKrY1H2Nsc1DSgSAKcO72Sw20h1u0N9P3gFEnd/I4MWYFgmAmjSlN6zMjq6xqxkv5zzl3dy3YHfKKXffcA9yQzrqhjbWxUuC4yWGdPi3yOXtjcyMlpMrWuESvpU6R4PbU5+IiSAhWvS0ZavrmRKZWljB+pNjYig1VVVsismrK0KmBXYiUn5JyjuT2UssVBoypK8mlSjZVkmI5wF69sO6TRKpEEOnd6JW/sPJw27W2UWMkJtYe6CHe7lBavg0asJDO9sfMw7aEu1VeJJNA506vo7Orm1e0H/Q4FUGIlfQi0R9vZpHYqsKI0MmIVKYsRyQwvbDpAXo6xeOoov0MRyRinTxpJQW5O2kwHKrGSEwoEow2YUz9i1dnVTVuaDO2KJMLyzY0snDiCYYW6bkgkUYoLclk0aUTatLeJK7Eys6VmttHrCv+5HvZ/xMwavW7yq8zsY4kPVfyQ6j6BUcfa2mg6UDLDgdYO1u0N8L4ZmgYUSbRzpleyYX8LDS1Bv0PpO7Eys1zge0Q6w88GbjSz2T089EHn3ALv9uMExyk+iY5Ypb54PdqIWQXskhmi0xQqXBdJvPd5dYsvpcGoVTwjVqcDW5xz25xzncADRLrESxZoPjpileIaq2KNWElmeWFzIyNK8pk7Nj3W2hHJJLPHDGdkaQHL06DOKp7EKt7O71eb2Roz+7WZje9hvwxBx4rXUzwVWBptxKwRKxn6Im1sDnDO9Cpycnrq+Swig5GTY5w1dRQvbj7g+0VP8SRW8XR+/wMwyTk3D3gG+FmPL2R2m5mtNLOVjY2N/YtUfOFXjZUaMUsm2VjfQmNLh6YBRZLo3OmVNLR0sLmh1dc44kms+uz87pw76Jzr8O7+CFjY0ws55+52zi1yzi2qqlIB51AQCIYozs+lIC+1F5BWFHsjVmprIxlAbWxEku8cr84q+nnzSzx/LV8HppvZZDMrAG4g0iX+KDMbE3P3CmB94kIUPx1s7Tw6epRKBXk5DCvMU42VZITlmw8wY/QwxpQX+x2KSMYaV1HMlMpS35dd6DOxcs6FgTuAp4gkTA8559aZ2ZfN7ArvYZ80s3Vmthr4JPCRZAUsqbV2bzOzasp8OXakrY0SKxnagqEuXt1+SKuti6TAOdMreXXbITrC/q2BGNf8jnPucefcDOfcVOfc17xtX3LOPeL9/Hnn3Bzn3Hzn3BLn3IZkBi2p0doRZnNDK/NqK3w5fqStjaYCZWh7bfshOsPdmgYUSYFzplXSHurizZ1NvsWgldelV2v3NOMcLBjvT2KlESvJBC9saqQgL4czJquNjUiynTl1FLk5xotb/KuzUmIlvVpTF8n459X6s+6ORqwkE5QV5XPZyWMoLsj1OxSRjFdWlM8p4yt87RuohlXSq9W7m6kdUcyoYYW+HH9ESb6K12XI+9SF0/0OQSSrnDO9kv9+djNNbZ1Hu3ikkkaspFer65qY71N9FUTa2rQEw4S7ukKta5gAAA45SURBVH2LQUREhpZzp1fiHKzYetCX4yuxkh4dbO2g7nA788f7134j2oi5qV3TgSIiEp/5tRWUFeb51t5GiZX0aE1dM4BvVwTCsbY2KmAXEZF45eXmsHjqKN8K2JVYSY9W7W4ix+Dkcf6NWEXnxlXALiIi/XHu9Ep2H2pn58EjKT+2Eivp0Zq6JqZVD6O00L/rG6JTgYePaMRKRETid860yLpxL/gwHajESt7DOcfqumZfC9chstwCQJNGrEREpB8mV5YytryIV7alvoBdiZW8R93hdg4d6WSeTwuDRkV7FGrJBRER6Q8zY/GUUby67SDOuZQeW4mVvMdqb2HQBT6PWA0rzCMvx1RjJSIi/bZ4yigOtHaypaE1pcdVYiXvsaaumYLcHGb61Hw5ysyoKCnQVYEiItJvZ06NtJFK9XSgEit5j1W7m5g9djgFef7/eowszaexpcPvMEREZIipHVHMuIpiXlZiJX7q6nas3dPMfJ/6Ax7v1AkjeGXbQYKhLr9DERGRIcTMOGPKSF7ZdiildVZKrORdtjS00tbZxXyfC9ejLplbw5HOLlZs9a+hpoiIDE1nThnFoSOdbKpPXZ2VEit5l2jhup8rrsc6e2olZYV5PLl2v9+hiIjIELN4SurrrJRYybus3t1EWWEeUypL/Q4FgIK8HC44qZpl79SrGbOIiPTL+JEljKsoVmIl/llT18zJteXk5JjfoRy1dG4Nh9tCvLb9kN+hiIjIEHPm1FG8su0g3d2pqbNSYiVHBUNdrN8XSJv6qqj3zaiiKD+HJ9dpOlBERPpn8ZRRHG4LsamhJSXHU2IlR63fFyDc7dLmisCokoI8zptRzZNr96fsG4eIiGSGxVNGAvDy1tRMByqxkqNW706vwvVYS+fW0NDSwVtejCI9MbOvmNkaM1tlZk+b2Vhv+03e9jVmtsLM5sc8Z6mZbTSzLWb2Of+iF5FkqB1RwviRqauzUmIlR62pa6aqrJAx5UV+h/IeS2ZVk59rPKXpQDmxbzjn5jnnFgCPAl/ytm8H3u+cmwd8BbgbwMxyge8BlwKzgRvNbHbqwxaRZFo8eRSvbj+UklkPJVZy1Kq6JubXlmOWPoXrUeXF+Zw1tZIn1+5PeUNNGTqcc4GYu6WA87avcM4d9ra/AtR6P58ObHHObXPOdQIPAFemKl4RSY3FU0bR1BZiw/7k11kpsRIAAsEQ2xqPMD8NpwGjLp1bw65DbbyzL9D3gyVrmdnXzGw3cBPHRqxi/SXwhPfzOGB3zL46b5uIZJDFKewbqMRKAHi7rhmAeWl2RWCsC2ePJsfgKS0WmtXM7BkzW9vD7UoA59wXnHPjgfuAO4577hIiidWd0U09HKLXIVEzu83MVprZysbGxsS8IRFJunEVxUwYWaLESlInuuJ6ul0RGKtyWCGnTRqpZReynHPuQufc3B5uDx/30PuBq6N3zGwe8GPgSudc9OxaB4yPeU4tsPcEx77bObfIObeoqqoqMW9IRFJi8ZSRKamzUmIlQOSKwImjSqgoKfA7lBO6dG4Nm+pb+f/t3X1sVfUdx/H397YrpR1PLULlsa2iLhBhrmIxU7e4KVsM+DAzFl2I8zGZm9tfc3Ex7sFsxu2PJduyGHU4M53T6QbKBgaT6R8iJQt0QEGwqK0gMCqgqyCl3/1xD3q5tHh777n9nWM/r+Tm3Ht6evnklH7z7e/87u+8tm/47vsk6WFms3JeLgK2RvtnAE8D33T3V3OOaQNmmVmTmVUBS4Dlw5VXRIbPgjPqOfj+UTreLu90EjVWAmQ/EZjk+VXHXTa7AUD3DpTB/CK6LNgOXAbcEe2/G6gHfhctxbAewN37yF4uXAV0AH9x980BcotImV3QlJ1nVe71rCrL+u6SCnsPHWb3wcOcm+DLgMdNGT+audPHs2rz23z7i2eGjiMJ4+7XDLL/JuCmQb62ElhZzlwiEt6U8aOZWV/D2s4ebrqouWz/jkashI3RxPV5CZ64nmvh7Abauw/y1oH3Q0cREZEUWdBcz7qd+zlWxnlWaqyEjV0HqMgYs6ckf8QKsquwgy4HiojI0LQ213PocB8dZVy2R42VsLH7AGdNHsPoqorQUQrSNLGWcxrGaNkFEREZktbm8q9npcZqhHN3NnYdSPQyCwO5fHYDbW/0sO/dI6GjiIhISjSMq6ZpYm1ZJ7CrsRrhXt/fy6HDfcxNyfyq4xbOacAdVm/RqJWIiBSutbmOdTt7yjbPSo3VCNceLQyahk8E5jqnYQyN9TWaZyUiIkPS2lzPu0f62LKrPPOs1FiNUO8d6ePBlzq597kOPj2qkrMmjwkdaUjMjMvnNPDya/s52Hs0dBwREUmJ4/OsXu78b1nev6DGyswWmtk2M9thZncO8PVRZvZE9PVXzKwx7qASj/3vHeFXq7dx4c/X8LPnOmg+rZY/3HA+n6pIX4+9cHYDff3OivZduJf3FgUiIvLJMHlsNc0Ta1nb2VOW9//YBULNrAL4LfBlsvfVajOz5e6+JeewG4F33P1MM1sC3Ad8vRyBpThdPb08+FInT6zv4vDRfi6fPZnbLjmDz86YEDpa0eZOG8/M+hp+9LdN/HrNduY31jG/qY7zG+s4p2EMmcxA99cVEZGRrvWMelZs2EXfsX4qYx5YKGTl9fnADnfvBDCzPwOLgdzGajFwT/T8KeA3ZmYewzBCf7+zZuveUt/mE63fnWP9Tl+/c6y/n6PHcl4f62dD1wFWtO8mY3DlvKncekkzZ05K16W/gWQyxpO3LmD1lj20vd5D284envvPbgDGVlfS0phtsk4fV40ZZMw+3GZ7ruzWTA1YKb5w9mmpHPEUkZGrtbmex155k827DsX+4a1CGqupQFfO627ggsGOcfc+MztI9r5cJ1zANLNbgFsAZsyYUVDAY+7c/Mf1BR0rA6upquCGCxu58aImTh83OnScWE0aW831rTO5vnUmAN3v9LJuZw9tr/ewbmcPL6gpL7v2ey5TYyUiqdLaVAfAy537gzRWA/05nz8SVcgxuPsDwAMALS0tBY1mVWaMZ7/z+UIOHdEqK4zKTIbKjFGRMSorom0mQ+2oCkZVpmPxz1JNm1DDtAk1XH3eNAB6/vcBB3o/wMmu2dXv4J4d5et3R1OzSldbpVuOiki6TBpbzVO3LSjLHUcKqYjdwPSc19OAXYMc021mlcA4IJZZYWbGnKnpWgpAkqOutoq62qrQMUREJGFaGuvK8r6FjN+3AbPMrMnMqoAlwPK8Y5YDS6PnXwNeiGN+lYiIiEiafOyIVTRn6nZgFVABPOzum83sJ8B6d18OPAQ8amY7yI5ULSlnaBEREZEkKmhyhLuvBFbm7bs75/lh4Np4o4mIiIikiz7KIyIiIhITNVYiIiIiMVFjJSIiIhITNVYiIiIiMVFjJSIiIhITNVYiIiIiMVFjJSIiIhITC7VAupntA94YwrdMJO+mzgmS5GygfKVSvuLlZ5vp7qeFChOnIdawJP+MQPlKkeRsoHylys1XUP0K1lgNlZmtd/eW0DkGkuRsoHylUr7iJTnbcEr6eVC+4iU5GyhfqYrJp0uBIiIiIjFRYyUiIiISkzQ1Vg+EDnAKSc4Gylcq5StekrMNp6SfB+UrXpKzgfKVasj5UjPHSkRERCTp0jRiJSIiIpJoiW6szOxsM9uQ8zhkZt8LnSuXmX3fzDab2SYze9zMqkNnymVmd0TZNifh3JnZw2a218w25eyrM7PnzWx7tJ2QoGzXRueu38yCfnJlkHz3m9lWM2s3s2fMbHzC8v00yrbBzFab2ZRQ+UJIeg1T/RpynsTWr1PkUw0rPltR9SvRjZW7b3P3ee4+D/gc0As8EzjWh8xsKvBdoMXd5wAVwJKwqT5iZnOAm4H5wFzgCjObFTYVy4CFefvuBNa4+yxgTfQ6hGWcnG0TcDXw4rCnOdkyTs73PDDH3c8FXgV+ONyhcizj5Hz3u/u50e/ws8Ddw54qoCTXMNWvoiwjufULVMNKsYyY6leiG6s8lwKvuftQFhUdDpXAaDOrBGqAXYHz5PoMsNbde929D/gXcFXIQO7+ItCTt3sx8Ej0/BHgymENFRkom7t3uPu2EHnyDZJvdfSzBVgLTBv2YB9lGSjfoZyXtcBIntSZxBqm+jUESa5foBpWijjrV5oaqyXA46FD5HL3t4BfAm8Cu4GD7r46bKoTbAIuNrN6M6sBvgpMD5xpIJPdfTdAtJ0UOE9afQv4R+gQ+czsXjPrAq5jhI1Y5UlUDVP9io3qV3wSV8OKqV+paKzMrApYBDwZOkuu6Fr6YqAJmALUmtn1YVN9xN07gPvIDrX+E9gI9J3ymySVzOwusj/bP4XOks/d73L36WSz3R46TwhJrGGqX5IkSa1hxdSvVDRWwFeAf7v7ntBB8nwJ2Onu+9z9KPA0cGHgTCdw94fc/Tx3v5jsMOf20JkGsMfMTgeItnsD50kVM1sKXAFc58leP+Ux4JrQIQJJYg1T/YqH6leJUlLDCq5faWmsvkGChtBzvAm0mlmNmRnZORQdgTOdwMwmRdsZZCcwJvE8LgeWRs+XAn8PmCVVzGwh8ANgkbv3hs6TL2+y8SJga6gsgSWxhql+xUP1qwRJrmFF1y93T/SD7ITK/cC40FkGyffj6GRvAh4FRoXOlJfvJWAL2WH0SxOQ53Gy8zmOAt3AjUA92U/TbI+2dQnKdlX0/AiwB1iVsHO3A+gCNkSP3ycs31+j3412YAUwNfT/wQDnJbE1TPVryHkSW79OkU81rPhsRdUvrbwuIiIiEpO0XAoUERERSTw1ViIiIiIxUWMlIiIiEhM1ViIiIiIxUWMlIiIiEhM1ViIiIiIxUWMlIiIiEhM1ViIiIiIx+T+AVZrF8HnpdQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "like_loglike_values = np.zeros((mu_candidate.shape[0],3))\n", "\n", "row_index=0\n", "for i in mu_candidate:\n", " like_loglike_values[row_index,0] = i\n", " like_loglike_values[row_index,1] = norm(i,2).pdf(y).prod()\n", " like_loglike_values[row_index,2] = norm(i,2).logpdf(y).sum()\n", " row_index+=1\n", "\n", "plt.figure(figsize=(10,5))\n", "plt.subplot(121)\n", "plt.title('Likelihood Function')\n", "plt.plot(like_loglike_values[:,0],like_loglike_values[:,1],label=\"Likelihood Function\")\n", "plt.subplot(122)\n", "plt.plot(like_loglike_values[:,0],like_loglike_values[:,2],label=\"Log-Likelihood Function\")\n", "plt.title('Log-Likelihood Function')\n", "plt.show()" ] } ], "metadata": { "hide_input": false, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" } }, "nbformat": 4, "nbformat_minor": 4 }