{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy as sp\n", "import scipy.stats\n", "from mpl_toolkits.mplot3d import Axes3D\n", "import sympy as sy\n", "sy.init_printing() \n", "import matplotlib as mpl\n", "mpl.rcParams['text.latex.preamble'] = r'\\usepackage{amsmath}'" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Perhaps this one of the most important application of linear algebra. We will build up intuition gradually before turning to multivariate normal distribution." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Univariate Normal Distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PDF of univariate normal distribution is given by\n", "\n", "$$\n", "p(x; \\mu, \\sigma^2) = \\frac{1}{\\sqrt{2\\pi}\\sigma}\\exp{\\left(-\\frac{1}{2\\sigma^2}(x-\\mu)^2\\right)}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $p(x; \\mu, \\sigma^2)$ mean random variable is $x$, parameters are $\\mu$ and $\\sigma^2$, it is not a conditional sign which commonly looks like $p(x|y)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that $-\\frac{1}{2\\sigma^2}(x-\\mu)^2$ is a quadratic function, which is the univariate version of quadratic form we have seen in earlier chapters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we define $\\sigma = 2$, $\\mu = 1$, let's plot the quadratic function and its exponential." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAD4CAYAAADb0iMYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3dd3xc5ZX4/8+RbbnL3ca23HtvcsMUAwZMCSYsLIZkAwFCSDbJN9kluyzkF5ZsNuVLXr/N/pKwCRAWElgTIAS8GLApJjQ3uXdbuMrGvRdZlnR+fzwzvndkyRppyp07c96vl16+d3Q1c2yPnjn3KecRVcUYY4wxxqRGXtABGGOMMcZkM0u2jDHGGGNSyJItY4wxxpgUsmTLGGOMMSaFLNkyxhhjjEmhxkEHcCEdO3bU3r17Bx2GMSZNli5dekBVOwUdRzJY+2VM7qmtDcvoZKt3794UFxcHHYYxJk1EZHvQMSSLtV/G5J7a2jAbRjTGGGOMSSFLtowxxhhjUsiSLWOMMcaYFLJkyxhjjDEmhSzZMsYYY4xJoaQkWyIyXUQ2ikiJiDxUw/ebisifIt9fJCK9k/G6xhiTKBF5RkT2iciaWr4vIvL/RdqvVSIyNt0xGmPCLeHSDyLSCPgNcDVQCiwRkdmqus532b3AYVXtLyIzgZ8Dtyf62iZLlJXBoUNw5sz532vZEtq1gyZN0h+XyRXPAr8G/lDL968DBkS+JgL/FfnTGOfMGdizB6qqYh9v1Qo6ZUXZOJOgZNTZmgCUqOoWABF5EZgB+JOtGcC/Ro5fAX4tIqKqmoTXN2Gh6pKqkhLYtcsdHzoEp07V/bMFBdC+PXToAD17Qr9+riEzJkGq+mEdve0zgD9E2quFItJWRLqq6udpCdBkrmPHYOFCWLq05ptFgB494OKLYdAgyLOZO7kqGclWd2Cn77yU8+/6zl2jqhUichToAByo/mQicj9wP0DPnj2TEJ4JVEUFfPaZS7BKSuDw4YY9z7Fj7mvbNtewAXTtCv37w4ABrkETSVrYxvjU1MZ1B85Ltqz9yhF798Knn8Lq1ef3ZlW3cyf86U/uZnHyZBg92nrqc1Aykq2aPuGq91jFc417UPVJ4EmAoqIi6/kKq6NHobgYli2DkycvfG3jxm6osGXL2ISpqgqOH4cjR2pu0D7/3H199BF07Ajjx7uGrGnT5P5dTK6z9ss4qvDOOy7Rqq59e2jTxjuvqoLSUqisdOeHDsGcOfDxx3DnndClS3piNhkhGclWKdDDd14I7K7lmlIRaQy0AQ4l4bVNJlGF7dth0SLYsMGdV5efD336QN++0Lmza6AKCi7cK1VZ6ZK3Q4dcclVS4u4W/QnYgQPw1lvw3nswahRMmGBzJUyyxNPGmWxXUQGvvQZrqq2j6NXLDRMOHHh+O3b8OCxeDEuWuLmp4NqyZ56BO+4A2zszZyQj2VoCDBCRPsAuYCZwZ7VrZgN3AQuAW4H3bb5Wltm7F+bOhS1bzv9eQQGMGOGG/Hr2hEaN6vfcjRq5pKx9e/ccl17q5kds3QqbNsHatd58ifJy17AVF8OYMXDFFdC6deJ/P5PLZgPfisxHnQgctflaOaaszA0Fbt3qPdavn2tfCgtr/7nWreGqq+CSS2D5cpg/37VVZ87AH/8It9wCw4alPn4TuISTrcgcrG8Bc4FGwDOqulZEfgQUq+ps4PfAH0WkBNejNTPR1zUZ4vhx14AsX35+T1afPq6HKRUTQ5s2hcGD3de118KqVe4Ocv9+931VN4S5Zg1MmeLuPG2ehKmBiMwCpgIdRaQUeBRoAqCqvwXeBK4HSoBTwFeDidQE4vhxeP55d0MZNWECTJ8ef7vWtClMmuTaxOefd89ZWQmvvOKOJ01KTewmY0gmdzAVFRVpcXFx0GGYmlRWunkLH33kepOiRFyP0uTJ6R/Giw5jfvyxG2r0KyhwSZndRWY0EVmqqkVBx5EM1n5lgSNH4Nln3Z9R0Z6qhi7IOXLEJVwHfOvDLr/c9ZKZ0KutDUvGMKLJNQcPwquvuvINfgMGwNVXu7lYQRBxcyB693bJ1rx5sG+f+96xY/Dyy24u2Q03QLNmwcRojAmHigp46SUv0crLg5tucotwEtG2LdxzD8ya5eaeAvz1r67dtJvBrGXJlolfdGju7bfh7Fnv8c6dXa9Rv37BxVZd//5uEv7y5fD++96KyNWrYccOuPlm16VvjDE1mTsXdkfWQTRqBDNnuhvKZGjRAr7yFTcPLNoLP3s2XHSRqyVoso5VWDPxOXnS3Yn97/96iVajRjBtGjzwQGYlWlF5eTBuHHznO25oM+roUfjDH1zPV0VFcPEZYzLT6tVuoU3UNdckL9GKatIEbr3Vlb0BN2n+pZdib2RN1rBky9Tt88/hd79zK/+iOnWC++5zcxcyvSpy06YwYwbcfru7owTXS/fpp/Dcc3XXATPG5I79+91NZdSwYW5CfCo0awZ/+7eu1iC4Sfhz5qTmtUygMvxT0gRuwwZXE+bYMe+xSZPg/vtdBfcwGTIEvvENN8QYtXMnPPVU7EojY0xuKi93vUvRRT8dOrh5WqncnaJrV7juOu98xQo3/cFkFUu2TM1U3aq+P/3J69Zu1gy+9CW35DmsZRRat3Z/h2uv9RrQI0fg97+P7bkzxuQWVXjjDa98TJMmrtcpHTtSjB3rijFHzZnjNrY2WcOSLXO+igp4/XV4912vdlb79m7YMNnzFoIg4kpT3Hmn15CWl7s5aQsW1Fz53hiT3TZtcvX6om64IX1b6oi414uu5I62wXXtu2hCw5ItE+vsWXjxRdeVHdWrl0u0OnYMLq5UGDAA7r3XLcUGl2TNneu2/LGEy5jccfasW2UdNXp04iUe6is/P3b+1uefu9XfJitYsmU8Z8+63h1/QdAxY9wS5ejE8mzTuTN87WvQw7f13ccfu5WKlnAZkxs+/RQOH3bHzZu71YdB6NjRbUcW9d57cOpUMLGYpLJkyzjl5fDCC7F7G15+uZscWt+9DMOmZUu46y639U/UggXuTtcSLmOy25EjbieMqKuuCvbm8uKLvXIQp0+7OoEm9CzZMq6+y/PPw7Zt3mNXXum2j0jlKpxM0rgx3HabW7EYtWiRm6hqCZcx2WvuXK/eXteubrJ6kJo0cYuQopYu9YqrmtCyZCvXRROtHTu8x66+Gi67LLiYgtKokSsy6N8yo7jY1dyxhMuY7FNSAuvXe+fXX58ZdQMHDvQWI6nCm29aGxRyGfCuMoGprHQ1ZaL7c4EriTBlSnAxBa1RI/ibv4ERI7zHli1zcyeMMdmjshLeess7Hz06du5mkERc71Z0CkdpKaxcGWxMJiGWbOUqVbcX12efeY9dd50riZDr8vLgi1+MrXvz8ceweHFwMRljkmvhQjh40B03beq2HsskHTq4+VtR77wDZWXBxWMSYslWrnrvvdg7palTYeLEwMLJOHl5boufgQO9x956K3bIwRgTTqdPw4cfeudXXAGtWgUXT20uvRQKCtzxyZNu1aQJJUu2ctHixa6nJmrcOLfy0MTKy3NzuAoL3bkq/PnPsH17sHEZYxKzaJGbrwqu3ML48cHGU5v8fDeHNmrRIpcomtCxZCvXrFsXO09h4EBXuThXVh3WV34+3HGH69IHt2pp1izYty/YuIwxDVNW5oYQoy67LLPL2wwb5hWUPnMmNnYTGpZs5ZLPP4dXX/VWtRQWup6bTFh9k8latoQvf9kbZigrcwmX3WEaEz4LF3pznzp0gOHDg42nLnl5sSMPCxda2xNC9imbK06dcptKR+vJdOjg9gbMzw82rrBo185tYB399zp8GF55xfYuMyZMqvdqXX55OG42rXcr9ELwLjMJq6pyicGRI+68aVOXaGXrFjyp0rUr3HKLd/7ZZ1bd2ZgwWbQoXL1aUdV7t/x/DxMKlmzlgvfei92G55ZbvDlIpn4GD44t+Prxx24enDEms5WVuW24oi67LBy9WlH+3q3qPXQm44XonWYaZO1a+OQT7/zyy2HQoODiyQZTp0L//t75a6/B/v2BhWOMiYO/N6h9+9jCxWGQlxd7o+efe2YyniVb2WzfPnj9de984ECXKJjE5OW5KvPRzWLLy+HFF63hMyZThXWuVnXDh3ujEmVlLoE0oZDQu01E2ovIOyKyOfJnu1quqxSRFZGv2Ym8ponT2bNuK57ycnfevr0bPrQSD8nRvDnMnOk2jQVXidr2UDQmMy1Z4q3gC2OvVlT13q0FC7x6YSajJZraPwS8p6oDgPci5zU5raqjI183JfiaJh5z58KBA+64SROXGDRrFmxM2aZLF1dlPmrtWtu/zJhMU1kZ2wMUtrla1Y0YEdu7tWJFsPGYuCT6jpsBPBc5fg64OcHnM8mwYQMUF3vn118PnTsHF082Gz7cVeCPevNNOHQouHiMMbHWrIETJ9xx69bh7dWKysuL3cN24UIrQRMCiSZbXVT1c4DIn7V9ojcTkWIRWSgiF0zIROT+yLXF+23Scf0dP+42mI4aOtTtZm9S59prvTvN8nK3pU9lZbAxGWPcsL5/rtaECZldLT5eI0e6qQzgav5t2hRsPKZOdSZbIvKuiKyp4WtGXT/r01NVi4A7gV+KSL/aLlTVJ1W1SFWLOnXqVI+XMKi6lXGnTrnzggL4whdsnlaq5ee7CfPRRnzXLvjrX4ONyRgDO3a4nTMAGjeO7YUOs/z82L+LlYHIeHUmW6o6TVWH1/D1OrBXRLoCRP6sccM4Vd0d+XML8AEwJml/A+NZuNAV2gSXYN1yi3f3Y1KrWze48krv/KOPbMNqY4LmT0JGjcquQs4TJnhzz7Zt85JKk5ESHUacDdwVOb4LeL36BSLSTkSaRo47AlMAqwKZbHv2wLvveudTpkDv3oGFk5Muvhj69HHHqm4fSisHYUwwDh9281ejJk0KLpZUKChw00SirHcroyWabP0MuFpENgNXR84RkSIReTpyzRCgWERWAvOBn6mqJVvJVFXl6mlF5wl16wZXXBFsTLlIBL74Ra838ehReOedYGMyJlctWuSVYunfH7JxWop/orx/IYDJOAklW6p6UFWvUtUBkT8PRR4vVtX7IsefquoIVR0V+fP3yQjc+Hz6aey8hFtuyY5JoGFUUAA33uidL10KW7cGF48xuejMGVi+3DvPtl6tqO7doUcPd1xZ6eqJmYwU4mIjBnC1tD74wDu/4gpv/ywTjKFDYcgQ73z2bFdk1mQsEZkuIhtFpEREzqsXKCI9RWS+iCwXkVUicn0QcZo4LV/uFfvs1An61bomK/z8ieSSJVBREVwsplaWbIWZqvsgj/5ydesW261sgiHiaptFi8gePgzz5wcbk6mViDQCfgNcBwwF7hCRodUu+wHwkqqOAWYCT6Q3ShO3qqrYIqYTJ2b3iuwhQ6BNG3d86hSsWhVsPKZGlmyFWXGxW9oMblXKTTeFuzJyNmnd2tXfilqwwJWEMJloAlCiqltUtRx4EVew2U+BgshxG2B3GuMz9bF5s7vBATd/ctSoYONJtbw8l1BG+eeqmYxhn8xhdeRI7OTrSy6Biy4KLh5zvtGjoW9fd6wau4jBZJLuwE7feWnkMb9/Bb4sIqXAm8C3a3oiK8qcAfy7Z4wb5+1fms3GjvX+nnv32o1dBrJkK4xU4Y03vE2mO3aM3ZzUZAYRV1Q22gju2+fqb5lMU9MYU/WugTuAZ1W1ELge+KOInNd+WlHmgB05AiUl7lgke4qY1qVZM7d1WJQ/4TQZwZKtMFq3LrZBmTHDrUI0maddO7jqKu/8o4/g4MHg4jE1KQV6+M4LOX+Y8F7gJQBVXQA0A2wlSqZZtswbQuvXz/3+5YqiIu94zRo4fTq4WMx5LNkKm/JymDvXOx8/3lv6azLThAlQWOiOKyvh7bdtTkVmWQIMEJE+IpKPmwA/u9o1O4CrAERkCC7ZsnHCTFJZ6ZKtqFzp1Yrq1s2bSlJRYRPlM4wlW2Hz0Udw7Jg7btkydosYk5ny8tzqxOiKqM2bbePYDKKqFcC3gLnAetyqw7Ui8iMRuSly2T8CX4sUZ54F3K1qGXNG2bTJK+rZujUMHBhsPOkmEtu7VVxsN3UZxMaewuTgQVfANOrqq73yAiazdevm7rSjcynefttNns+FybshoKpv4ia++x/7oe94HW6rMZOp/POUxozJzcLOI0bAvHluBGT/fti5E3r2DDoqg/VshYcqvPWWt5qtR4/sX9Kcba680tvK5/Bh+OSTYOMxJlscOgSffeaORdzqvFzUtKlLuKJsonzGsGQrLDZujJ0U7x+WMuHQokXsZPmPP/bqARljGs4/V6t/f2jbNrhYguafq7ZunSt0agJnyVYYnD3rhp2ixo2Drl2Di8c03Nix3v9dRUXsYgdjTP1VVsbug+ift5SLunVzX+DamJUrg43HAJZshcMnn7j6MXB+74gJl7w8uOEG73zDBq/H0hhTfxs2wMmT7rigAAYMCDaeTOBPOJcutYnyGcCSrUx37Fjs3J6rrvLm/ZhwKix0E3ij5s51+7kZY+rPPy9p7FjbsgxcgdOmTd3xgQOwfXuw8RhLtjLe/PluGBHc8JP/Q9qE17RpXmO4f3/sMIgxJj6HDsHWre44lyfGV5efDyNHeuf+OW0mEJZsZbK9e2HFCu/86qvtri1btGwJU3yVBObP97ZfMsbEx98+DhjghhGN4088162DsrLgYjGWbGW0d97xxtoHDPA2NTbZYfJk78PhxInYGmrGmAurqopNtqzXP1bXrrEV5desCTaeHGfJVqbasiW21MPVVwcbj0m+Jk3giiu8808/9SpgG2MubOtWbzeNFi1yr2J8PPwJqE1VCJQlW5lI1VUBjho9Gjp3Di4ekzqjRkGXLu64vBw++CDQcIwJDX/yMGpUblaMr8uIEd6/y65dsG9fsPHkMEu2MtGqVbBnjzuu3vthskteXmyv5bJlbsK8MaZ2p0+7kg9Ro0cHF0sma9ECBg/2zv3DriatLNnKNGfPwvvve+f+eT0mO/Xr583Hq6qCd98NNh5jMt3q1W4eErgCntHeYXM+fyK6cqW35ZtJK0u2Ms3ixXD0qDuuvmLNZKfonLzo9ksbN1pdHGMuxCbGx69fP++G/eRJ2Lw52HhyVELJlojcJiJrRaRKRGrdI0FEpovIRhEpEZGHEnnNrHbmTGwB06lTvVpMJrt17RpbF+f9963qszE12bMHdu92x40buwKepnZ5eW5OW5RNlA9Eoj1ba4BbgA9ru0BEGgG/Aa4DhgJ3iMjQBF83Oy1a5G0a2q6dFejLNVOnenXUtm/3ijUaYzz+Xq3Bg21HjXj4hxI3b7ZVzwFIKNlS1fWqurGOyyYAJaq6RVXLgReBGYm8blYqK4uts3TZZba6Jte0axc7JDJ/vvVuGeNXWekWEEXZEGJ8OnSAXr3ccVVV7L+hSYt0zNnqDuz0nZdGHquRiNwvIsUiUrw/l1ZlLVzoVfht3z6229fkDn+SvXOnbVJtjN/GjV7vf5s20KdPsPGEib93a/lyu5FLszqTLRF5V0TW1PAVb++U1PBYrf/LqvqkqhapalGnTp3ifImQO30aFizwzv3DSSa3tGkD48Z559a7ZYzHP4Q4erS1k/UxbJjbMxFceZldu4KNJ8fU+U5V1WmqOryGr9fjfI1SoIfvvBDY3ZBgs9ann7rJ8QAdO9qEz1x36aVu4i+4icCbNgUbjzGZ4MSJ2J5e6/2vn/x8l3BFrVwZXCw5KB23BUuAASLSR0TygZnA7DS8bjicPOkmxkdZr5Zp3RrGj/fOrXfLGFdbq6rKHffq5aZbmPrxJ6hr1ni1ykzKJVr64YsiUgpMBuaIyNzI491E5E0AVa0AvgXMBdYDL6nq2sTCziKffOK2aQG3JY//zsPkrilT3O4B4Ja6r18fbDzGBM3fE2O9Wg3Tqxe0beuOT5+2XvM0SnQ14l9UtVBVm6pqF1W9NvL4blW93nfdm6o6UFX7qeq/Jxp01jhxApYs8c6vuMIrbGlyW6tWMGGCd/7BB9a7ZXLXnj3eFmaNG8NQqx7UICKxiaoNJaaNjVcFacECtz0PuKKW/j2sjJkyxZvQum9f7F5wxuQSf1IwZAg0axZcLGHnT7Y2b3ZTWUzKWbIVlNOnY3u1LrvMerVMrBYtYnu3PvrIerdM7qleW8uGEBPTvj307OmOq6rcXDiTcpZsBWXxYm+uVqdO1qtlajZpUuzKxC1bgo3HmHT77DOv96V1a2/TdtNwNpSYdpZsBaG83BUxjbr0UuvVMjVr1Sp226aPPgouFmOC4E8GRo601drJMGyYdxP3+eewd2+w8eQAe9cGYelSN4wIbmWI1dUyF3Lxxd4HzLZtsGNHoOEYkzanT8fOVbQhxORo1ix2NMV6t1LOkq10q6iI3QPxkkvsTs1cWNu27o4+6uOPg4vFmHRau9bN2QLo1s2VxzHJ4U9cV63yapiZlLBP+XRbtQqOH3fHrVrF7ldlTG2mTPGGmjdt8pbBG5PN/NvzWK9WcvXr5z6DwJUh+uyzYOPJcpZspVNVVWyvxMUXe+PmxlxIp05uyXuU9W6ZbHfwIJSWuuO8PJtukWx5ebE95jaUmFKWbKXTunVw6JA7bt48dsNhY+pyySXe8dq17sPImGzl//AfOBBatgwulmzl7y3csAHKyoKLJctZspUuqrErySZOhKZNg4vHhE+3btC/vztWdVs9GZONVG17nnTo0sUV1AY3n3jdumDjyWKWbKXLli3e8tr8/NhilcbE69JLveOVK91cC2OyzfbtcPSoO27eHAYMCDaebGY1t9LCkq10WbDAOx4zxlUHN6a+evaE7t3dcWVl7C4EpsFEZLqIbBSREhF5qJZr/lZE1onIWhH5n3THmFP8H/rDh9vc1lQaPtxbEb99Oxw+HGw8WcqSrXTYtw9KStyxiKsKbkxDiMDkyd75kiXe/pqmQUSkEfAb4DpgKHCHiAytds0A4F+AKao6DPhu2gPNFWfPujmJUTaEmFqtWnnTEyB2aySTNJZspYO/V2vIEGjXLrhYTPgNHQpt2rjjU6escUzcBKBEVbeoajnwIjCj2jVfA36jqocBVHVfmmPMHRs2eFuZdejg9eSa1Kk+lGh7sCadJVupduJE7Iehv1fCmIbIy4vtHV2wwBrHxHQHdvrOSyOP+Q0EBorIJyKyUESm1/REInK/iBSLSPH+/ftTFG6W8w8hjh5tW5mlw6BBrqo8uBXzO3de+HpTb5ZspdqSJV4F5MJC92VMosaM8VazHjjgDVObhqjp07x69toYGABMBe4AnhaRtuf9kOqTqlqkqkWdOnVKeqBZ7/jx2OKa/jpQJnUaN3b7JUbZRPmks2Qrlc6ejZ3APHmy3aWZ5GjWLHaDav9QtamvUqCH77wQ2F3DNa+r6llV3QpsxCVfJplWr/Z6afv08YbLTer5hxLXrnWlIEzSWLKVSqtWuTk14Pa381cANyZREyd6yfuWLbaFT8MtAQaISB8RyQdmArOrXfMacAWAiHTEDStuSWuU2U7VtucJUo8e3nzisjLYuDHYeLKMJVupohrb2zBxom04bZKrbVs3WT7KercaRFUrgG8Bc4H1wEuqulZEfiQiN0UumwscFJF1wHzg+6pqJfyTae9et3IboEkTuzlNNxGruZVC9umfKps3u7k04ObW+Id8jEkW/4KLNWu8Tc5Nvajqm6o6UFX7qeq/Rx77oarOjhyrqv6Dqg5V1RGq+mKwEWchf6/WkCG2w0YQ/MlWSQmcPBlcLFnGkq1UWbjQOx471hoOkxqFha77H6zIqQmvyko3XyvKhhCD0a6dK5wMUFUV+39iEmLJVirs3+/m0IDrmp04Mdh4THbz924tXWoTW034fPaZ14vSurWbHG+C4U90/b2NJiGWbKWCv3dh0CA3t8aYVBk8GAoK3PHJk7aZrAkf/4f6yJE2vzVIw4Z52yPt2WMLb5IkoXe0iNwW2SesSkSKLnDdNhFZLSIrRKQ4kdfMeGfOxE4stA2nTarl5UGR79fPhhJNmJw+HbvybfTo4GIxrqyMf3GCTZRPikRvH9YAtwAfxnHtFao6WlVrTcqywqpVLuEC6NjRusNNeowdC40aueOdO2F39TJRxmSoNWu8ws/du4MVgw2efyhx1Srv/8c0WELJlqquV1UrxhGlCosXe+fjx1sRU5MerVrFVoC23i0TFv4hROvVygx9+7q5c+CmJvir+psGSdfAuALzRGSpiNx/oQtDvbfYtm1ucjxAfr41HCa9xo/3jlev9grqGpOp9u+HXbvccaNGMHx4sPEYJy/PJsonWZ3Jloi8KyJraviaUY/XmaKqY4HrgL8XkctquzDUe4v5exNGjbJyDya9Cguha1d3XFEBy5cHG48xdfF/iA8aBM2bBxeLieVPtjZutJu3BNWZbKnqNFUdXsPX6/G+iKrujvy5D/gLkH2zxo8dgw0bvHN/L4Mx6SASuyCjuNjVyjEmE1VVuflAUTYSkFk6dXJz6MDN2VqzJth4Qi7lw4gi0lJEWkePgWtwE+uzi/+DrXdv6Nw50HBMjho+3OsdOHzYVYE2JhNt2eLteNCqFfTrF2w85nz+BNhWJSYk0dIPXxSRUmAyMEdE5kYe7yYib0Yu6wJ8LCIrgcXAHFV9O5HXzTgVFa6YZJSVezBBadIkdmso/4INYzKJfwhxxAhvNa3JHMOHe/8vu3Z5c5JNvSW6GvEvqlqoqk1VtYuqXht5fLeqXh853qKqoyJfw6L7jmWV9eu96scFBa7IpDFBKSryVsGWlMChQ8HGY0x1ZWWx0y5sCDEzNW/u5tJF2UT5BrMyvcng79UaN86qH5tgtWsHAwZ458uWBReLMTVZu9bbVqprV+jSJdh4TO38ifCqVTYPtIEsK0jUwYOu5AO4JMs/hGNMUMaN846XL7eihCazWG2t8OjXz82pAzfHzmpuNYglW4ny92oNHOgVgjMmSAMGxO6XuNFqD5sMsX+/2+UArLZWGDRq5ParjLKhxAaxZCsRFRWxbzzr1TKZIi8Pxozxzv03BcYEyV//bdAgaNkyuFhMfPy9jxs2WM2tBrBkKxH+Qm8FBdC/f7DxGOM3Zow3Uf6zz1wpCGOCVFkZW0LAf0NgMlfnzm3ilHcAACAASURBVK5oMrj/Q399NBMXS7YS4e8tGDvWJsabzNK2bewNgFWUN0HbvDl25bbV1goPf2K8fLnbC9jEzbKDhjp0yBXlA9d7YHdoJhP5h7aXL7eVRCZY/pWxo0bZDWqYDB/u6vgB7N0Lu3cHG0/I2Du9ofyNxoAB0KZNcLEYU5uBA2NXEm3aFGw8JncdP+56tqLsBjVcmjaFYcO8c+sprxdLthqistImxptwaNTIJsqbzLBypTf01Ls3tG8faDimAfxtyerVcPZscLGEjCVbDbFpE5w44Y5bt3a9B8ZkKv/NQEkJHD0aXCwmN6nG9oRYr1Y49ewJHTq44zNnYN26YOMJEUu2GsLfOzBmjM07MJmtXTvo29cdV//QMyYdduxwBaDBDUcNHRpsPKZhqs9PtrYkbpYl1NeRI14FXZsYb8LCX1F+2TKbKG/Sy/+hPGKEN9HahI9/YcO2bbb3apws2aov/7yDvn1dr4ExmW7wYK945LFjsHVrsPGY3HHmjNsLMcpuUMOtdWsrKdMAlmzVh2rsxHhrNExYNGrkehSibMsNky5r1ngTqTt3hm7dgo3HJM4/D3TFCuspj4MlW/WxY4dXhbtZM7fVhDFh4d9yY/16KCsLLhaTO6rPcY3uamDCa8AAr6e8ekkPUyNLturD3xvgL/BmTBhcdJH7Arevp39ox5hU2L3bK37ZuLGb72PCr3pJmeLi4GIJCUu24lVeHvvh5O8lMCYs/O9bG0o0qeb/EB42DFq0CC4Wk1z+RTclJW7xmKmVJVvxWr/eJVwAHTtC9+7BxmNMQ4wY4a0k2rkTDhwINh6TvcrKXOHLKP+Hswm/du28ifKqVjC5DpZsxcvfCzB6tM07MOHUsmXsXMOVK4OLxWQ3f4Xxzp2hR49g4zHJV1TkHS9f7nZXMTWyZCseR454S+VFbN6BCTf/UOLKlbaSyCSfauwQYlGR3aBmo4EDXSkIcLuqbNwYbDwZzJKtePjv/vv3995cxoRR//5Wc6saEZkuIhtFpEREHrrAdbeKiIpIUW3XGKC0FPbudcdNmsDIkcHGY1IjLy+2DIRNlK+VJVt1qV5byybGm7Br1Cj2wy/HixKKSCPgN8B1wFDgDhE5bz8ZEWkNfAdYlN4IQ8j/oTtihCuVY7LT2LFer+WWLd62TCaGJVt12b7damuZ7OO/adiwIddrbk0ASlR1i6qWAy8CM2q47t+A/wvk9D9WnU6fjl25XWSdgFmtTRs3nBhlE+VrlFCyJSKPi8gGEVklIn8Rkba1XBdXF31G8vdqjRjhasUYE3ZdukDXru64osJV+c5d3YGdvvPSyGPniMgYoIeqvnGhJxKR+0WkWESK9+/fn/xIw2DlSveeAlct3irGZz9/Qr1ihff/b85JtGfrHWC4qo4ENgH/Uv2CeLvoM9LZs7BunXduQ4gmm1SfKJ+7apq5ree+KZIH/Afwj3U9kao+qapFqlrUqVOnJIYYEtUnxlu5h9zQr5/r4QI4dSr2c9MACSZbqjpPVaMp7EKgsIbL4u2izzwbN8bW1rI7NJNNhg+PrbkVHS7PPaWAvy5BIbDbd94aGA58ICLbgEnAbJskX4Nt27zabU2bxu7HabJXXl5sYr1kSXCxZKhkztm6B3irhsfr7KL3y6hu+FWrvOORI23psskuLVt6RQkhtgBlblkCDBCRPiKSD8wEZke/qapHVbWjqvZW1d64G8ubVNWWXlW3yLd2YORIyM8PLhaTXmPHusU34G7edu++8PU5ps5kS0TeFZE1NXzN8F3zCFABvFDTU9TwmNbwmPtGpnTDnzrltiCIsjs0k4387+tVq9wwUI6J9M5/C5gLrAdeUtW1IvIjEbkp2OhC5PDh2DpLEycGF4tJv1at3JZMUYts0a5fnbO9VXXahb4vIncBNwJXqdbYUtfVRZ+Z1q71ij326OG2JjAm2wwa5Hofysvd8M+ePd7E+Ryiqm8Cb1Z77Ie1XDs1HTGFzuLFXrLev7+bemFyy8SJ3ojQmjVw9dUuCTMJr0acDvwzrkv9VC2XXbCLPmNVH0I0Jhvl58OQId65/31vTLzOnIFly7xz69XKTd27Q2Fk6nZlpRU59Ul0ztavcZNH3xGRFSLyWwAR6SYib0LtXfQJvm5qHT7sxpzBTfzzd40ak238Q4mrV9v2Pab+Vq50CRdAhw6xcwFNbpk0yTsuLrYyEBEJFY1S1Rp/o1R1N3C97/y8LvqM5p8o3L8/tGgRXCzGpFrfvm6y/MmTbn+zrVvdUm5j4qEaOz9n4kRbTJTLhgyBggK3FdiJE25Kju0nbBXkz6NqQ4gmt+Tlnd+7ZUy8Skq8LVqaNbN6hLmuUSMYP947X7QoJxfeVGfJVnWff+7VicnPt+15TG7wJ1vr1rmCvsbEw9+rNWaMlXswruZWdLeV3bvdxuQ5zpKt6vx39UOGuB3rjcl23bq5uTbgVib6l/AbU5sDB7wSOSIwYUKw8ZjM0KJF7A3cwoXBxZIhLNnyq6qKTbZsCNHkCpHY97utSjTx8PdqDRpkJXKMxz9Rfv16OHo0uFgygCVbflu3ugl94GqD9OkTbDzGpJP/TrSkxBX2NaY2p065TYejrNyD8evSBXr3dsdVVTlf5NSSLT9/r5Z/3zhjckH79l6NnKoq20zWXNjixd7cvosu8j5YjYmaPNk7Li6G06eDiyVglk1EVVTAhg3euW3PY3KR/32/Zk1wcZjMVl4e21NxySVW7sGcb+BAiG67V16e0xtUW7IV9dlnUFbmjtu1cxOGjck1Q4d6H5rbt8Px48HGYzLT0qVeL0W7du59Y0x1Ii4Rj1q0KGdXOluyFbXWV9R+2DC7SzO5qXVr6NXLHavaUKI5X2UlLFjgnU+ZYlMuTO2GD4c2bdzxyZOwfHmw8QTEfkPAZdr+IUTbnsfkMv/7f21m76xlArBqlasODm4hkRUxNRfSqBFcfLF3/umnLmHPMZZsgVt5VV7ujjt0cJM9jclV/qHEHTtyfsm28VGFTz7xzidN8opXGlObsWO9be+OHMnJmzhLtiB2IrANIZpc17JlbNkTG0o0URs2eDtsNG0KRUXBxmPCoUmT2NIgH3+cc1v4WLJVXg6bNnnnw4cHF4sxmcL/e2CrEg24D8ePP/bOx493eyEaE48JE7ytnPbtg82bg40nzSzZ2rTJWx3RubP7MibXDRniTXretQsOHw42HhO8bdvcewHc0KG/QrgxdWne3O2ZGOVP3HOAJVvVhxCNMa5h7NfPO8/BORbGRxX++lfvfPRoNznemPqYPNlNmAc3H3Tr1mDjSaPcTrbOnPE2UQUbQjTGz//7YMlWbtu61fVsgevxnDIl0HBMSBUUxK5eff/9nJm7ldvJ1oYNrnI8uBWIHToEG48xmWTQIO8u9PPP4eDBYOMxwVB1H4pRY8bYhtOm4S67zGtXdu6M7fDIYrmdbPnv1q1Xy5hYzZrBgAHeufVu5aZNm6C01B03bgyXXx5sPCbc2rSJXcWaI71buZtsnT7ttuiJsvlaxpzP/3thqxJzjyrMn++dFxW5oSBjEnHppa4cBLhe8/Xrg40nDXI32dq40ati2727dYsbU5NBg7xGcd8+r8aSyQ3r1sGePe64SZPYfe6MaahWrVwpiKj586GqKrh40iB3ky1/oUbbRNWYmuXnQ//+3rkVOM0dVVWxvVoTJ9oKRJM8U6a4wrgA+/dnfc95biZbZ87EDiFasmVM7fy/HznQ3W8iVq2KrRZvKxBNMrVo4UpBRM2fn9V7JiaUbInI4yKyQURWichfRKRtLddtE5HVIrJCRIoTec2k2LTJ+0+96CIbQjTmQgYOjF2VaAVOs19lZWxdrYsvdrXXjEmmSZO899Xhw7BiRbDxpFCiPVvvAMNVdSSwCfiXC1x7haqOVtXgN9OyIURj4te0aWyBUxtKzH7FxV5S3aKFVYs3qdGsWWyP6QcfuC30slBCyZaqzlPVSKEqFgKFiYeUYuXlsXU9hgwJLhZjwsL/e2JDidnt1KnYuVqXXOLNrTEm2SZM8OYCHj+etdv4JHPO1j3AW7V8T4F5IrJURO6/0JOIyP0iUiwixfv3709ieBElJd5eiJ06uS9jzIUNGuTtlVhaCseOBRuPSZ3586GszB23bx+7asyYZMvPh2nTvPNPP4UjR4KLJ0XqTLZE5F0RWVPD1wzfNY8AFcALtTzNFFUdC1wH/L2IXFbb66nqk6papKpFnVKRCPnvyq1Xy5j4tGgBvXt759a7lZ327nVDiFHXXusKmRqTSqNGuRJM4HZ1mTcv2HhSoM5kS1WnqerwGr5eBxCRu4AbgS+p1lwGVlV3R/7cB/wFCOZWqaLCTY6PsvlaxsTP//ti87ayjyq8/bZXzbtfP7c4wphUE4Hp073zdeu8vTizRKKrEacD/wzcpKqnarmmpYi0jh4D1wDBFNTYssWVfQDXPd6lSyBhGBNKgwe7RhFgxw44cSLYeExybdjgNpwGN2R87bXe/7cxqdajB4wc6Z2/9VZWFTpNdM7Wr4HWwDuRsg6/BRCRbiLyZuSaLsDHIrISWAzMUdW3E3zdhvHfjQ8ZYg2JMfXRqhX07OmOVd2Hc5YQkekislFESkTkoRq+/w8isi5S5uY9EekVRJwpU33oZvx46Nw5uHhMbpo2zduxYu9eWLYs2HiSKNHViP1VtUekpMNoVX0g8vhuVb0+crxFVUdFvoap6r8nI/B6q6x0W/RE2RCiMfWXhQVORaQR8BvcnNKhwB0iUr2BWA4URcrcvAL83/RGmWILFnilHpo3h6lTAw3H5KiCArdvYtT777t9jLNA7lSQ37bN+09r0wa6dQs0HGNCyb+oZOvWbGkIJwAlkRvDcuBFYIb/AlWd75sqEY4yN/E6cgQ++sg7v/JKK2BqgjN5MrSN1Ec/dcolXFkgd5Kt6qsQbQjRmPorKIDCSJ5RVRXbWxxe3YGdvvPSyGO1uZfay9yEiyrMnu0VkuzcGcaNCzYmk9uaNIFrrvHOlyzJisnyuZFsVVVZyQdjksX/+5MdqxJruvOqcWW1iHwZKAIer+X7qa0TmGzLl7uFQ+BuQG+6yaunZkxQhgxxtf2iZs/26mOGVG4UUCkthZMn3XHLlm7VQxpVVVVRWlrKyWgMxiRRy5YtKSwsJC9dH5JDh8I777jjLVtcr0h+fnpeOzVKAX+jUAjsrn6RiEwDHgEuV9UzNT2Rqj4JPAlQVFRUY8IWhGPHjrFv3z7O+j+wqqrc/92117rzpk1dBe8smYtnQm7kSOjTxytFsnp1jcPbaW//Gig3ki3/UIe/EnaaHDhwABFh0KBBGf+GMOFSVVXFrl27OHDgAJ3TtXqsXTtXNmXvXreK7bPPwt5bvAQYICJ9gF3ATOBO/wUiMgb4HTA9Ui8wNI4dO8bevXvp3r07zZs3R0TcB9ihQ+7/EtxG4506Wa+WySynTsVWk+/YMebGLpD2r4Gy/zer+hL1wYPTHsKRI0fo0qWLJVom6fLy8ujSpQtHjx5N7wv7u/hDXgIisr/rt4C5wHrgJVVdKyI/EpGbIpc9DrQCXo6UuZkdULj1tm/fPrp3706LFi1cogVuYcMZX+dc27aWaJnM07x57L6cR454PV0E2P41QPb3bB04AAcPuuP8fOjbN+0hVFZW0iRaO8SYJGvSpAkVFRV1X5hMgwfDhx+6402b3JBUiD+sVfVN4M1qj/3QdzztvB8KibNnz9LcP/xSWQn+D6cWLWyjaZOZRFz1gP37XZJVUeGGugsKzl0SSPvXAOFtHePlv+vu3z+wfb7EVj+aFAnkvdW1q9fgnT4N27enPwYTt3PvEVWXaEV7Bxo1ivngMibjNG4c+x49ccJbPUt4PluzP9mqPl/LGJM4kdgh+ewoAZH9Tp6EsjLvvE2bUPdImhzRokXsIpzDh0O3lU92/5YdP+5WIoJrUGxTVWOSp/q8rZr3oTeZorwcjh3zzlu2hGbNgovHmHiJuHmF0V6sykqXcIWozcnuZMt/t92rl1VFzkLTpk2jY8eO/PjHPw46lNzTu7f3YX3kiFudaDJTVZW3HQ+4wpE2fBgqOd/WNW7sVZYHt8AjROWUsjvZCngVokm9Z599ll/84hdBh5ESS5cuZcqUKVx22WVceeWVbIkWn8wUjRrBgAHeechXJWa1I0dcbwC43oF27WwXjRTbsWMHO3bsSNrzZWtbd6F27uOPP469uHlz1yMbdexYzPytTJa9ydaZM27vtiibr5WVCguzZ4u66rp168bbb7/Nhx9+yIMPPsijjz4adEjns3lbme/Mmdh5Wm3bBrZQKFecOHGC5557jp49eybtObO1rbtQO9eqVStmzZoV+wMFBa5nNiok87eyN9kqKfHu5C66KLb70eS8uXPncql/d/mA3X777fz+97+Peaxr1660bt0agPz8fBpn4gdk//6uhwvg889jCxCa4O3YEbtZeMuWNp0iDX7xi1/wpS99KegwgMxv6y7Uzo0ePZqFCxdy2D8EXr1ntrLSvcczPOHK3mTLhhBNLVSV733vezz22GNBh3LOY489xsMPP8xp/wdjxMmTJ3nooYd48MEHA4isDk2bui01oqx3K3McPAj+XgGbp5UW+/fvZ+nSpfQNoKZjdWFq62pr526++WZ+9rOfxT5J48be7gfg9k18662MnjCfnclWZSVs3uydW7JlfObNm0d5eTlXXHFF0KGcM3jwYPr3739el3l5eTm33XYbP/jBDxg2bFhA0dXBhhIzz8mT8PzzXq9WXp7N00qTV155hUsuuSToMIDwtHUXaucuv/xyXnnlFbR6ItWsGbRq5Z0vWQKffprK0BOSncnW9u3eHIW2bd0+bqZOVVVVtGnThnnz5sU8/sUvfjEze1WAe+65h8cff5xnn32WL3zhC3H9zGuvvca0adNiiuH96le/onfv3pyJbGGyfv16LrroIl5++eWkxPnHP/6RUaNGxTz27W9/m29+85vnzq+++mpee+21c+eVlZXceeed3HLLLdx8881JiSMl/PMht22LHbYy6Xf2LPzP/3irD0WgfftQzNNKdRt06tQpHnzwQfr06UP79u2ZPn06JSUlAOzevZsuXbrw/PPPn7v+3nvv5YorrqAyMiWld+/e/OhHP+KSSy6hVatWFBUVsWTJkpjXmDNnDmPGjEk41urC0NbF085BbFtXVzuXl5dH586dKS4uPv8FW7eOLV/yzjuwZk3Cf49UyM5kyz+EOGiQ3c3FKS8vj4kTJ7Jo0aJzj7377rssWLCAH/7whzHXfvOb36Rt27a1fp3X7ZsizzzzDGvXrqWkpIT//d//jetnli1bxtChQ2Me+8Y3vkGzZs144okn2Lp1K9dccw0//elPue2225IS54oVKxg/fnzMY8XFxYwdO/bc+YgRI1i2bNm585dffpm3336b559/nqlTp/Ltb387KbEkXevWEJ28W1UV26ts0quqCv78Z9i1y52LnF8QMoOlug2677772LBhAwsXLmTPnj1MnDiRG2+8kbNnz9KtWzdeeOEFvvnNb7J+/Xr+8Ic/MGfOHGbNmkWj6LxE4Le//S3/+Z//yaFDh7j11lu5/vrrOearX7Zw4UIGpqCmYxjaunjaOYht6+Jp5wYNGsTChQvPf8Ho/C3/jcRf/pKRO1pk/q1Ofam6vdqiMnEV4r/+a8a+5uTJk1m8eDEAFRUVfPe73+UnP/kJBdXmejzxxBM88cQTDQ7n7rvv5rnnnqv1+4888kjK6skcPnz4vL9P48aNefzxx/nqV7/Kr3/9a77//e/z1a9+tcafb0jsy5cv5/bbbz93XlFRwcqVKxk3bty5xwoKCjh06NC585kzZzJz5sx6/d0CM2iQV0B40yYYOTLYeHKRKrz9duzN5nXXxa7cgoxufyB1bdCBAweYNWsW27dvp0tktOPRRx/ll7/8JYsWLeKSSy5h2rRp/MM//AMzZsxgz549vPbaa1x00UUxz3Pvvfee+73953/+Z5544gneeOMN7rzzTs6ePcuhQ4do06bNea8fRJuXSFuXqnYOYtu6eNq5goIC9tZWx0/ELfzo1MntoVhZ6eYq3nuveyxDZF+ytX+/tyKqaVNXzNTE7eKLL+a//uu/ANeYNW/evNakIxG//vWvL1gzpkWLFkl/zah27drF3IlGDR8+nJMnTzJu3Di+853v1PrzDYl9xYoVPP744+fO161bR2VlJcOHDz/32LFjx2jfvn28f43MMmgQvPeeO46uBPb1BpgUU4W5cyGSpAAwZQpMmADr1wcXVwOkqg3aGikFNLLajcDZs2fZuXPnufMHHniAn/70p0yaNIkrr7zyvOfp3bv3uWMRoWfPnpRGbjQOHDiAqp5bXecXRJuXSFuXqnYO6t/WFRQUcODAgdovEIEvfQmeftrtnVhWBs89B1/5CnTuHPfrpFL2JVv+Xi3/snQTl0mTJnHgwAGKi4t57LHHmDNnTo0bfT7wwAMxcxuqe/jhh3n44Ydr/X6rVq1o5Z/cWIdkbjb69a9/nXXr1sU8tnv3bqZNm8bXv/51fve737FhwwYG17Kwor6xb9u2jcOHDzNixIhzj82dO5cRI0bQxNfrsGbNmpTM9UiLTp3c/MgjR1xDt3OnqzBvUk8V5swB/5yW4cNh2rTgYkpAqtqgXpEb782bN9Oplh6Pqqoq7rrrLm688UYWLFjAM888wz333BNzzbZt284dqyo7duw4VwMrmoCUl5efV6ol3nYjWW2dqjJmzJgGt3Wpaueg/m3dmTNn6k5G27aFO++EZ591hU5PnHDHf/d30LVr3K+VMqqasV/jxo3Tenv6adVHH3VfK1bU/+dTYN26dUGHUC/Dhg3Tvn376le+8pWgQ2mwP//5z3rffffpjBkzdO/evTHfe+utt7R///7nzvft26dDhgzRxx57TFVV7733Xr3hhhuSFsurr76qgM6bN08rKir0gw8+0Hbt2umXvvQlPXPmzLnrpkyZok8//XSDXiMj3mNz5ni/e2+/3aCnAIo1A9qeZHw1qP2qr8pK1b/8xft3f/RR1ZdeUq2oOHdJRrw36ilVbdCdd96pt956q5aWlqqq6uHDh/XVV1/V48ePq6rqY489pkOGDNETJ07o/PnztXXr1rp69epzP9+rVy/t1q2bLl26VMvLy/XnP/+5dujQQY8cOXLumqZNm+ru3buTGveFZEpbF287p1r/tu5rX/ua/uxnP6v1+zHv8e3bVX/yE+/34ac/VY38f6dDbW1Ydk2QP3XKmzciEruViInb5MmT2b9/f9omucejvlvX3HLLLTz11FN89atfPW9i5bXXXkvjxo354IMPOHr0KNdeey3XX3/9uQm4jz32GO+//z7vvvtuUmJfsWIF1113HT/4wQ/o1KkTTzzxBF//+teZM2fOuZVMGzduZPPmzdx5551Jec1A+CcF+3uYTWpUVbnJwCtWeI+NHAl/8zeh79FPVRv01FNPMWjQIKZOnUrr1q0ZMWIEL7/8MiLC/Pnz+cUvfsHLL79My5YtmTp1Kv/0T//EbbfdxknfHnz3338/3/nOd2jXrh1/+tOfmDNnTswcrb59+9Y+vygOYW3r4mnnoGFt3d69e+nXr198F/fs6YYPo6sUy8rgD39wBX6DVFMGVp8v4N+AVcAKYB7QrZbr7gI2R77uiue5631nuGKFl802sIcgFcJ2Z3nVVVfpz3/+86DDiLF79249duyYqqrOmTNHv/zlL9f5M2fOnNF77rlHT548ed733nrrLb300kuTHmdNvvCFL+gvf/nLC14zc+ZMfeqppxr8GhnxHjt7VvXHP/Z+Bw8cqPdTYD1b8SkrU33hhdgerddecz1d1WTEe6OeMrENUnU9W3/84x8veM33vvc9/e///u8Gv0ZY27p42jnVhrV1ffv21cOHD9f6/Rrf47t3q/78597vx49/rLp+fb1etyFqa8OS0bP1uKqOVNXRwBvAD6tfICLtgUeBicAE4FERaVf9uoRl+irEEHjyySfZs2cP3/ve94IOJcaFtnSYNWsWkydPZtiwYeTl5TFkyBAqKyv57ne/y8MPP1zjWP/06dP58MMP0xL78uXLY+Yx1GTWrFncd999aYknZRo3Bv/dp/VupcbBg24isP/fd/x4uOkmV7w05DK1DYrXzTffzIIFCxr882Ft6+Jp56D+bd3nn39Onz59aFvfLfe6doW77/YKn549Cy++CH/9ayCV5hP+zVRV/1KHlkBNf4trgXdU9ZCqHgbeAaYn+toxKivdKqioFNQ5yWaLFy+mTZs2PPHEE7zyyivnTWjMFNW3dNi0aRO/+tWvmD9/PmvWrGHo0KF88skn/Md//Adr1qzh8ccfT6jhS9TBgwcpLS09bzVO1rKhxNQqKYGnnnKrrqMuuQSuvz709QTD0gbV5bLLLmPfvn3nioY2VJjaulS2c88//zyPPPJIw364c2eXcPm39pk/H156yU2iT6OkrEYUkX8HvgIcBWraF6A7sNN3Xhp5rKbnuh+4H6jfjuk7drjd7cGtSsig+hphMGHCBI4ePRp0GBdU05YOr7/+OnfffTfNIuPzTZo0oXHjxjz44IMZUfW+Q4cO0WH03OCfJxndycFf4dk0jKrbiuTdd7278saNXW9WltQ0C0Mb5F+JeCEPPfQQs2bN4u67727Q64StrUtVO1dWVsauXbv4/ve/3/An6dgRvvY1ePlliJT/YP1610N8xx2xiVgKxdWzJSLvisiaGr5mAKjqI6raA3gB+FZNT1HDYzX+z6jqk6papKpFtS3PrZH/LnrgwNDf5ZlYtW3pcPr0acoiWzPNmzeP3r17n1fEz6RR69bQrZs7rqqCzz4LNp5scOwYvPCC24ok+oFWUAD33JM1iVa2mThxIs2aNWO/vwcyTtbWeV555RUeffTRxJ+oRQtXAmLSJO+xffvgt7+F5cvTMqwYV8+WqsZbsOV/gDm4+Vl+pcBU33kh8EGcHXhVHwAACHtJREFUzxmf6smWySrRLR0OHDjA888/z4gRI/jVr37FXXfdxe23386rr75Kly5dePrpp4MO1QwcCLt3u+NNmyBTN9DOdKqwahW89Za31ytAjx5w++2xm/CajDNz5sxzyVF9WFvnufXWW8/15CUsLw+mT3d7Jb/xhpt6dOYMvP46rFsHX/iCu4lJkYSHEUVkgKpGN0O7CdhQw2VzgZ/4JsVfA/xLoq99zsGD7gvcHmBWTDHr1LalQ69evWreM8sEZ+BA+OADd7x5s+vhyoKJ22l1/Lj7QNi40XtMBCZOdMVKQ7CptKFBiYK1dZ6kJVp+Y8a4aUavvgrR7dE2b4YnnnDbW40cmZKRsWT8xv5MRAYBVcB24AEAESkCHlDV+1T1kIj8GxAttvEjVT1U89M1gL9Xq29fa4iMCVLXrm448fhxV/tu1y7XG2Pqdvas23Lno49ie7PatYObb7btx4xJhsJC+MY33BZjixa5XuSyMle3bvlyuOYabzpEkiSclajq39TyeDFwn+/8GeCZRF+vRv67PxtCNCZY0YLCy5a5840bLdmqS3TI8P33ofok8QkTXG9Wfn4wsRmTjZo0ccOKQ4bAa6/B4cPu8W3b4MknYcQIuOoqt+AuCcLfBVRWFlsZ1qrGGxO8QYO8ZGvTptDu05cWW7bAvHmwZ0/s4x06wI03Qp8+wcRlTC7o1cvr5VqyxE17AFi92s3lmjgRLr0UmjdP6GXCn2yVlHj/ON27u+ELY0yw+vRxw/kVFW7Vz5EjSbtDzCqrV8Of/xz7WMuWMHUqjB0b+m13jAmF/Hw3X2v8eFdeZUNk6nllpSu5smEDfOtbCc09Df+s1c2bvWMbQjQmM+Tnx/bI+H9PjWfwYG8FVJMmcNll8J3vuEbfEi1j0qtjR5g505VVKSz0Hp8wIeFFPuHv2brhBhg61A1VDB4cdDS1UlXEan+ZFMjYoqmjR7sKzoMGxTZcxtOkiZsXsn07XHFFynrmq6qqyLMVoSYLpaT969kT7r3XFT9duhSKihJ+yvAnW/n5rjHP4L0QmzVrxsGDB+nQoYMlXCapVJWDBw+mZol0ooYNsxpb8Rg1yn2lSMuWLdm1axddunShSZMm1gaZrJHS9k/EdeQMHZqUpwt/shUChYWFlJaWNqiasDF1adasGYXWc2RqUVhYyIEDB9i+fTsVFRVBh2NMUoWl/bNkKw2aNGlCH1tRZIwJQF5eHp07d6Zz585Bh2JMzrJBfGOMMcaYFLJkyxiT80RkuohsFJESEXmohu83FZE/Rb6/SER6pz9KY0xYWbJljMlpItII+A1wHTAUuENEqs+KvRc4rKr9gf8Afp7eKI0xYWbJljEm100ASlR1i6qWAy8CM6pdMwN4LnL8CnCV2LI+Y0ycLNkyxuS67sBO33lp5LEar1HVCuAo0KH6E4nI/SJSLCLFtvrYGBOV0asRly5dekBEtsd5eUfgQCrjSRGLO73CGjeEN/b6xN0rlYHUoqYequqVEuO5BlV9EngSQET250D7BeGN3eJOr1yJu8Y2LKOTLVXtFO+1IlKsqomXeU0zizu9who3hDf2EMRdCvTwnRcCu2u5plREGgNtgEMXetJcaL8gvLFb3OmV63HbMKIxJtctAQaISB8RyQdmArOrXTMbuCtyfCvwvmbsPknGmEyT0T1bxhiTaqpaISLfAuYCjYBnVHWtiPwIKFbV2cDvgT+KSAmuR2tmcBEbY8Imm5KtJ4MOoIEs7vQKa9wQ3tgzPm5VfRN4s9pjP/QdlwG3pTCEjP83uoCwxm5xp1dOxy3WE26MMcYYkzo2Z8sYY4wxJoUs2TLGGGOMSaGsTLZE5EERURHpGHQs8RCRx0Vkg4isEpG/iEjboGO6kLr2kctEItJDROaLyHoRWSsi/yfomOpDRBqJyHIReSPoWOIlIm1F5JXIe3u9iEwOOqYwsPYrtaz9Sr8wtl+Q3DYs65ItEekBXA3sCDqWengHGK6qI4FNwL8EHE+t4txHLhNVAP+oqkOAScDfhyTuqP8DrA86iHr6T+BtVR0MjCJ88aedtV+pZe1XYMLYfkES27CsS7Zwm8T+EzVUd85UqjovsgUIwEJcUcVMFc8+chlHVT9X1WWR4+O4X5rqW7JkJBEpBG4Ang46lniJSAFwGa5kAqparqpHgo0qFKz9Si1rv9IsjO0XJL8Ny6pkS0RuAnap6sqgY0nAPcBbQQdxAfHsI5fRRKQ3MAZYFGwkcfsl7gO4KuhA6qEvsB/478jwwdMi0jLooDKZtV9pYe1X+oWx/YIkt2Ghq7MlIu8CF9XwrUeAh4Fr0htRfC4Ut6q+HrnmEVx38QvpjK2e4tojLlOJSCvgz8B3VfVY0PHURURuBPap6lIRmRp0PPXQGBgLfFtVF4nIfwIPAf9PsGEFy9qvwFn7lUYhbr8gyW1Y6JItVZ1W0+MiMgLoA6wUEXBd2ctEZIKq7kljiDWqLe4oEbkLuBG4KsO3AYlnH7mMJCJNcA3VC6r6atDxxGkKcJOIXA80AwpE5HlV/XLAcdWlFChV1ejd9yu4hiqnWfsVOGu/0ius7RckuQ3L2qKmIrINKFLVjN9lXESmA/8vcLmq7g86nguJbMK7CbgK2IXbV+5OVV0baGB1EPcJ9hxwSFW/G3Q8DRG5M3xQVW8MOpZ4iMhHwH2qulFE/hVoqarfDzisULD2KzWs/QpO2NovSG4bFrqerSz1a6Ap8E7krnahqj4QbEg1q20fuYDDiscU4O+A1SKyIvLYw5FtWkxqfBt4IbK58xbgqwHHY1LD2q/Us/YrGElrw7K2Z8sYY4wxJhNk1WpEY4wxxphMY8mWMcYYY0wKWbJljDHGGJNClmwZY4wxxqSQJVvGGGOMMSlkyZYxxhhjTApZsmWMMcYYk0L/P1G9jkCCmdQxAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sigma = 2\n", "mu = 1\n", "x = np.linspace(-4, 6)\n", "y = -1/(2*sigma**2)*(x-mu)**2\n", "\n", "fig, axs = plt.subplots(1, 2, figsize = (10, 4))\n", "axs[0].plot(x, y, lw = 3, color = 'r', alpha = .5,\n", " label = r'$y=-\\frac{1}{2\\sigma^2}(x-\\mu)^2$')\n", "axs[0].legend(loc ='best', fontsize = 13)\n", "\n", "axs[1].plot(x,np.exp(y), lw = 3, color = 'r', alpha = .5, \n", " label = r'$y=\\exp{\\left(-\\frac{1}{2\\sigma^2}(x-\\mu)^2\\right)}$')\n", "axs[1].legend(loc ='best', fontsize = 13)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The constant in front, $\\frac{1}{\\sqrt{2\\pi}\\sigma}$, is a normalizing factor which ensures the integral of the whole function equals to $1$. \n", "\n", "$$\n", " \\int_{-\\infty}^{\\infty} \\frac{1}{\\sqrt{2 \\pi} \\sigma}\\exp \\left(-\\frac{1}{2 \\sigma^{2}}(x-\\mu)^{2}\\right)dx=1\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Easiest way to plot a univariate normal PDF is to use Scipy's normal distribution function, ```sp.stats.norm.pdf()```, we can specify the $\\mu$ and $\\sigma$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEvCAYAAAC6xJMcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deZhU5Z33//eX7qZZZN/3JYDa7NjsiAqiaIgYNQnqJGbijL/MjM/MJJn5TfLML5vJk+vJzFyPmZkkk/jEGDMuaDQaEnFBQUWQpdkaENEWEZpFGhoQ6Iamu+/fH3cVp6qppqvXU8vndV11ee5Tp6q+VTb1rfuc+/7e5pxDRERE0ku7sAMQERGRxlMCFxERSUNK4CIiImlICVxERCQNKYGLiIikISVwERGRNJQbdgCN0bt3bzd8+PCwwxAREWkTmzZtOuqc65PovrRK4MOHD6eoqCjsMERERNqEmX1U3306hS4iIpKGlMBFRETSkBK4iIhIGlICFxERSUNK4CIiImkoqQRuZgvNbLeZlZjZNxPc/3Uze8fMis3sNTMbFnPfPWb2fuR2T8z+q8xse+Q5/8PMrGXekoiISOZrMIGbWQ7wM+AmoAC408wK6hy2BSh0zk0AngH+JfLYnsB3genANOC7ZtYj8pj/Au4DRkduC5v9bkRERLJEMj3waUCJc26Pc64KWAosjj3AObfKOVcRaa4DBke2bwRWOOfKnXPHgRXAQjMbAHR1zr3t/ILkvwVubYH3IyIikhWSKeQyCNgf0y7F96jrcy/w4iUeOyhyK02wX0TCVlkJFRWXPqZjR+jUqW3iEZGEkkngia5Nu4QHmv0ZUAhc08BjG/Oc9+FPtTN06NCGYhWRptq/H9auhXffBZfwn2O8yy+HmTNh2DDQEBaRNpdMAi8FhsS0BwMH6x5kZtcD/wxc45w7F/PYa+s89vXI/sF19l/0nADOuYeAhwAKCwuT+FYRkaTV1vqEvXYtlJY2fHys3bv9beBAn8gLCiAnp3XiFJGLJJPANwKjzWwEcABYAtwVe4CZTQZ+CSx0zh2Juetl4EcxA9duAL7lnCs3s1NmNgNYD3wJ+M/mvRURSVpVFWzZAuvWwfHjF9/fo0f9vWrn4h9z8CA8+yy8+ipMnw5TpkCHDq0Tt4hc0GACd85Vm9n9+GScA/zaObfTzB4Aipxzy4B/BS4DfheZDbbPOXdLJFH/AP8jAOAB51x5ZPuvgN8AHfHXzF9ERFpXRYXvbRcVwdmz8ffl5MD48b433a/fpZ/n6FF4+23Ytg2qq/2+kyfhlVfgjTd8Ep89Gy67rHXeh4hgLplrXSmisLDQaTUykSY6ehQeewxOnIjf37EjTJ3qb126NO45z5yBjRv97cyZ+Pu6doW77274x4CI1MvMNjnnChPepwQukgX27YMnn/QjzKN69vS97YkToX375j3/+fNQXOx75UePBvvz8+ELX4CRI5v3/CJZ6lIJPK3WAxeRJti5E557LjjVnZcHt9wCY8dCuxaqppyXB1dd5U+dv/suPP88nDvnb4895l9v0qSWeS0RAZTARTKXc36Q2ssvB/s6d/antQcObJ3XNIMrr/S9+8cfh08+8SPdn3/eXyOfO1dTzkRaiBYzEclEtbXw0kvxybt3b/iLv2i95B2rXz//WrHXv1etgj/+EWpqWv/1RbKAErhIpjl/Hn73O1i/Ptg3dCh85St+elhb6doV/vzP469/b97sr8WfO1f/40QkKUrgIpnkzBl49FHYtSvYN3YsfOlL4ZQ+7dDBn7KfODHYV1ICv/kNnDrV9vGIZBAlcJFMUVkJjzwSX1Ft1iy44w7IDXG4S04O3HorXHNNsO/QIXj44YunnolI0pTARTKBc/D73wdTuMzgppvghhtSY9CYGVx3nR+NHh35fuKEP9VfWxtubCJpSglcJBO8/jq8/37Qvu02X9Y01UyZAp//fPCjYu9eX4JVRBpNCVwk3e3e7cuXRs2Z40uipqorrvC98ai1a2HHjvDiEUlTSuAi6ezYMV+kJWrkSJg3L7x4knX11X450qg//AGOHKn/eBG5iBK4SLqqqoKnngoWJenWzQ9Ya6nqaq3JDD77WejVy7fPn4elSy9eYEVE6pUG/9JF5CLOwbJlQa81N9fXHA9jqlhTdejgY47WYS8v9wPx0mh9BpEwKYGLpKN16+KvGy9a1DYV1lpa376weHHQfu89ePPN8OIRSSNK4CLp5sMPYcWKoF1YmN4LhYwd6+erR9UdUS8iCSmBi6STkyfhmWeCudODB8PCheHG1BKuvx5GjPDbzsGzz/pT6iJSLyVwkXRRXQ1PPx1UL+vc2c+pDrPKWktp184PwOva1bfPnvUD9Kqqwo1LJIUpgYukizffhAMH/Ha7dj55RxNeJujc2Q9qy8nx7Y8/hpUrw41JJIUpgYukg7IyWLMmaC9YAMOGhRdPaxk0CG6+OWivX+/rpovIRZTARVKdc/CnPwXraA8ZAjNmhBtTa5oyJViC1Dm/hrjqpYtcRAlcJNVt3QoffeS327XzU8ZSYYGS1mLm32P02v7Bg7BxY7gxiaQgJXCRVHbmDLzyStCeNQv69QsvnrbSsyfMnRu0V66ETz4JLx6RFKQELpLKXnnFr/MN0L17/JramW7WLOjd22+fOwcvvhhuPCIpJqkEbmYLzWy3mZWY2TcT3D/XzDabWbWZ3RGz/zoz2xpzO2tmt0bu+42ZfRhzXxpXohBpBR9+CNu2Be1Pfxry8sKLp63l5sJnPhO0d+3yK6+JCJBEAjezHOBnwE1AAXCnmRXUOWwf8GXgididzrlVzrlJzrlJwDygAog5H8g/Ru93zm1t+tsQyTDV1X7gWtTYsTB6dHjxhGXYMJg8OWgvX6654SIRyfTApwElzrk9zrkqYCmwOPYA59xe51wxcKmhoncALzrnKpocrUi2eOstv1QoQH5+ZlRba6oFC4JFWk6e9KVWRSSpBD4I2B/TLo3sa6wlwJN19v0vMys2swfNLL8JzymSeY4ehdWrg/b8+dClS3jxhK1TJ7jhhqC9bh0cPhxePCIpIpkEnmi+SqPW+zOzAcB44OWY3d8CrgCmAj2Bf6rnsfeZWZGZFZWVlTXmZUXSj3PwwgvBnO9Bg/xiJdlu4kQYPtxv19b6ywuaGy5ZLpkEXgoMiWkPBg428nU+DzznnDsf3eGcO+S8c8Aj+FP1F3HOPeScK3TOFfbp06eRLyuSZoqL/eA18HO+P/MZ/99sF50bHi2zWloKmzaFG5NIyJL5ZtgIjDazEWbWHn8qfFkjX+dO6pw+j/TKMTMDbgV2JHicSPaoqICXY05SzZgB/fuHF0+q6d0brr46aL/6Kpw6FV48IiFrMIE756qB+/Gnv3cBTzvndprZA2Z2C4CZTTWzUuBzwC/NbGf08WY2HN+Df6POUz9uZtuB7UBv4IfNfzsiaWzlSp/EAbp1g2uvDTWclDRnDvTq5bfPnYsvciOSZZJah9A5txxYXmffd2K2N+JPrSd67F4SDHpzzs1rTKAiGe3oUdi8OWjffDO0bx9ePKkqN9efSn/0Ud/evt0XfBkwINy4REKgi2siqWDlymBQ1siRcPnl4caTykaMgCuuCNqvvRZeLCIhUgIXCduBA/DOO0F7/vzwYkkX8+YFC7qUlMDevaGGIxIGJXCRsMX2IAsK/NQxubS+ff3UsqhXX/VT8ESyiBK4SJj27PE38NPF5mloSNKuuy5+WpnqpEuWUQIXCYtzvucYNXlysPqWNKxbN5gWUz7itddU3EWyihK4SFjeeQcORmoi5eZm11KhLeXqq32teICysvjV20QynBK4SBhqavzI86jp06Fr1/DiSVedOvlpZFGvv+5XchPJAkrgImHYujVYbaxDB1+gRJpm5kzo3NlvnzwJGzeGG49IG1ECF2lr58/HL4k5Zw507BhaOGmvffv4yw9vvglnz4YXj0gbUQIXaWvr1wc1vLt08afPpXmuugp69PDblZWwdm248Yi0ASVwkbZUWQlvvRW0r7kG8vLCiydT5OT4aWVRb78Np0+HF49IG1ACF2lLa9YEp3d79fJTx6RljB8P/fr57fPn/al0kQymBC7SVj75BNatC9rz5gWFSKT5zOD664N2URGUl4cXj0grUwIXaStvvBFMcRo40JdNlZY1ahQMG+a3a2th1apw4xFpRUrgIm3hxAnYsiVoz58fLMYhLaduL3zHDl/gRSQDKYGLtIW33grKfA4f7pcMldYxZAiMGeO3nYPVq8ONR6SVKIGLtLZTp+J733Pnqvfd2ubODbZ37NC1cMlISuAirW3tWl86FWDwYBgxItx4ssHgwcFZjtpaP/pfJMMogYu0pjNn/GjoKPW+205sL3zrVj8LQCSDKIGLtKZ16/ycZID+/WH06HDjySbDhsHQoX67pkbV2STjKIGLtJazZ2HDhqB99dXqfbclM/+ZR23a5M+IiGQIJXCR1rJhA5w757d794Yrrww3nmw0ahQMGOC3z5/3JVZFMkRSCdzMFprZbjMrMbNvJrh/rpltNrNqM7ujzn01ZrY1clsWs3+Ema03s/fN7Ckza9/8tyOSIqqq4quuXX01tNPv5TZnFn8tfONGX49eJAM0+I1iZjnAz4CbgALgTjOrW0JqH/Bl4IkET1HpnJsUud0Ss//HwIPOudHAceDeJsQvkpqKiqCiwm/36AHjxoUbTza74gro08dvnzsXf1lDJI0l0yWYBpQ45/Y456qApcDi2AOcc3udc8VAbTIvamYGzAOeiex6FLg16ahFUll1dfyAqdmzVfM8THWvha9bF1zaEEljySTwQcD+mHZpZF+yOphZkZmtM7Noku4FnHDOVTfxOUVS15YtwVKWXbrApEnhxiP+DEjPnn67sjJ+ap9ImkomgScaNusa8RpDnXOFwF3AT8zsU415TjO7L/IDoKhMNY0l1dXUxK/3PXs25OaGF4947drBnDlB++23g+l9ImkqmQReCgyJaQ8GDib7As65g5H/7gFeByYDR4HuZhb9Zqv3OZ1zDznnCp1zhX2i17FEUlVxMZw86bc7d4YpU8KNRwITJ0LXrn779On48rYiaSiZBL4RGB0ZNd4eWAIsa+AxAJhZDzPLj2z3BmYD7zjnHLAKiI5Yvwf4Q2ODF0kptbXxve8ZM6C9JlekjJwcf0Yk6q23ghK3ImmowQQeuU59P/AysAt42jm308weMLNbAMxsqpmVAp8DfmlmOyMPvxIoMrNt+IT9v51z70Tu+yfg62ZWgr8m/nBLvjGRNvfOO3DsmN/u0AGmTQs3HrnYlCn+zAj40qrbtoUbj0gzJHVxzjm3HFheZ993YrY34k+D133cWmB8Pc+5Bz/CXST9OQdvvhm0p0+H/Pzw4pHE8vJg1ixYscK333rLDzLUHH1JQ/qrFWkJ770HR4747fbtfQKX1FRYCB07+u3ycti589LHi6QoJXCRlhA777uwEDp1Ci8WubT8/PgfWGvX+jMoImlGCVykuUpL4aOP/Ha7dn7wmqS2adOC6X2HDsHevaGGI9IUSuAizRXb+x4/PpiqJKmrUyeYPDlor1kTXiwiTaQELtIc5eWwa1fQnjUrvFikcWbODJZ3LSkJxjCIpAklcJHmePvt4PrpqFHQr1+48UjyevaMX+I19kyKSBpQAhdpqooK2Lo1aKv3nX5i/59t3+7nhoukCSVwkabauDGop92/P4wYEW480niDB8PQoX67pgbWrw83HpFGUAIXaYrz5+PXlZ49O7ieKukltrxqUZGWGpW0oQQu0hTbtsGZM367WzcoKAg3Hmm6MWOgd2+/fe4cbNoUbjwiSVICF2ms2lo/eC1qxgy/UIakJzM/Ij1q3TotciJpQQlcpLF2745ftERLhqa/iRPjFzlReVVJA0rgIo1Vt2yqFi1Jf7m58eVV16xReVVJeUrgIo2xf7+/gT9trkVLMkdhoV+tDODjj2HPnnDjEWmAErhIY8SW3JwwAbp0CS8WaVmdOsVfDlFhF0lxSuAiyTp2zF//jood+CSZYcaMYDrgBx/A4cPhxiNyCUrgIsmKLZs6ejT07RtuPNLyevSInxKoXrikMCVwkWScORNfNjW2+Idkltj/tzt2wMmT4cUicglK4CLJ2LABqqv99sCBMGxYuPFI6xk4EIYP99u1tX5euEgKUgIXacj5877uedSsWSqbmuliFznZvFnlVSUlKYGLNKS42K88Biqbmi1Gj44vr7p5c7jxiCSgBC5yKc7Fn0KdMQPa6Z9NxqtbXnX9en86XSSF6JtI5FJKSqCszG/n58PkyeHGI21nwgQ/NxzgxAnYtSvceETqSCqBm9lCM9ttZiVm9s0E9881s81mVm1md8Tsn2Rmb5vZTjMrNrMvxNz3GzP70My2Rm6TWuYtibSg2EVLpkzxtc8lO+TlwdSpQTt2GqFICmgwgZtZDvAz4CagALjTzOpeBNwHfBl4os7+CuBLzrmxwELgJ2bWPeb+f3TOTYrctiKSSg4fDsppmqlsajaaOjVYaa601N9EUkQyPfBpQIlzbo9zrgpYCiyOPcA5t9c5VwzU1tn/nnPu/cj2QeAI0KdFIhdpbbHXvgsKoHv3+o+VzHTZZf5UelTsGRmRkCWTwAcB+2PapZF9jWJm04D2wAcxu/9X5NT6g2aWcEknM7vPzIrMrKgsei1SpLWdOgXbtwdtlU3NXrH/73ftguPHw4tFJEYyCTzRhNdGXQgyswHAfwN/7pyL9tK/BVwBTAV6Av+U6LHOuYecc4XOucI+fdR5lzayYQPU1PjtoUNh8OBw45Hw9O0Lo0b5bef8iHSRFJBMAi8FhsS0BwMHk30BM+sKvAD8f865C+cknXOHnHcOeAR/ql4kfFVVUFQUtNX7lti/gc2b4ezZ8GIRiUgmgW8ERpvZCDNrDywBliXz5JHjnwN+65z7XZ37BkT+a8CtwI7GBC7SarZtg8pKv92jB1x+ebjxSPhGjgwWr6mqgk2bwo1HhCQSuHOuGrgfeBnYBTztnNtpZg+Y2S0AZjbVzEqBzwG/NLOdkYd/HpgLfDnBdLHHzWw7sB3oDfywRd+ZSFOocIskkqiwS/QSi0hIcpM5yDm3HFheZ993YrY34k+t133cY8Bj9TznvEZFKtIW3nvPr/sNfs63CrdI1Pjx8NprcPo0fPIJvPOO3ycSEnUtRGLFThO66ipo3z68WCS15OaqsIukFCVwkaiDB2HvXr/drp0Kt8jFCgt9Igf/97JvX7jxSFZTAheJiu19jx0LXbuGF4ukps6dYeLEoK3CLhIiJXARgJMnYefOoK2pY1Kf2L+N3buDMRMibUwJXAR84ZbocpHDh8PAgaGGIymsd28YM8Zv1521INKGlMBFzp2Ln9er3rc0JPZvZOvWoG6ASBtSAhfZsiWorNWrV9C7EqnP8OHQv7/fPn8+vnKfSBtRApfsVlsbfwp05kxftEPkUsxg1qygvX49VFeHF49kJSVwyW67dsGJE367U6f4EcYilxI7U+H06fjV60TagBK4ZC/nYO3aoD11KuTlhRePpJecnPhaASrsIm1MCVyy1759cOCA365bZUskGbHV+o4cgQ8+CDceySpK4JK9YotwTJgAl10WXiySnjp0gClTgnbsGR2RVqYELtnp2DFfhCNKU8ekqWbMCAY+7tkDhw+HG49kDSVwyU7r1gXXK0ePhj59wo1H0lf37lBQELRVXlXaiBK4ZJ+KCj/3Oyp2OpBIU8T+DW3f7pcbFWllSuCSfTZuDObsDhjgi3KINMegQTBsmN+urfXzwkVamRK4ZJfqal/3PEqFW6SlxI6j2LTJl+gVaUVK4JJdiovhzBm/3bWrL8Yh0hIuv9yX4gVfmjf2Mo1IK1ACl+zhXPwAoxkzfDEOkZZgFt8LX7cuWOFOpBUogUv2KCmBsjK/nZ8fP39XpCVMnOhL8oIv0btrV7jxSEZTApfsEVtkY8oUX4RDpCXl5cVX9Fu7VuVVpdUogUt2OHQIPvzQb7drF1/DWqQlTZ3qS/OCL9W7b1+48UjGSiqBm9lCM9ttZiVm9s0E9881s81mVm1md9S57x4zez9yuydm/1Vmtj3ynP9hpqHA0opir30XFPjiGyKt4bLLfGneKBV2kVbSYAI3sxzgZ8BNQAFwp5kV1DlsH/Bl4Ik6j+0JfBeYDkwDvmtmPSJ3/xdwHzA6clvY5HchciknT8KOHUFbhVuktcUOZtu9G44eDS8WyVjJ9MCnASXOuT3OuSpgKbA49gDn3F7nXDFQd8jljcAK51y5c+44sAJYaGYDgK7Oubedcw74LXBrc9+MSEJvvx2MBh42DAYODDceyXx9+sCYMX677rK1Ii0kmQQ+CNgf0y6N7EtGfY8dFNluynOKJK+iwhfViJozJ7xYJLvMnh1sb9um8qrS4pJJ4ImuTSc7rLK+xyb9nGZ2n5kVmVlRWXQKkEiyNmyA8+f9dr9+MGpUuPFI9hg6FIYM8ds1NX5euEgLSiaBlwJDYtqDgYNJPn99jy2NbDf4nM65h5xzhc65wj5aMUoao6oqvib1nDkqmyptxyz+jE9REVRWhhePZJxkEvhGYLSZjTCz9sASYFmSz/8ycIOZ9YgMXrsBeNk5dwg4ZWYzIqPPvwT8oQnxi9Rvy5bgC7N7d5VNlbY3ZkywVG1VlU/iIi2kwQTunKsG7scn413A0865nWb2gJndAmBmU82sFPgc8Esz2xl5bDnwA/yPgI3AA5F9AH8F/AooAT4AXmzRdybZraYmfuDQ7Nl+/rdIW6rbC1+3LrikI9JMuckc5JxbDiyvs+87MdsbiT8lHnvcr4FfJ9hfBIxrTLAiSduxw08fA+jcGSZNCjceyV7jxsHKlf7v8cwZ2Lo1vlqbSBOpSyKZxzlYsyZoT5/uS1yKhCEnJ35e+Nq1WuREWoQSuGSe996DI0f8dvv26u1I+KZMgY4d/fbx47BzZ7jxSEZQApfM89ZbwXZhYfDFKRKW9u3j6++vWaNFTqTZlMAls+zbB/sjtYNycvya3yKpYNq04FLO4cPwwQfhxiNpTwlcMkts73viROjaNbxYRGJ16gRXXRW0Y/9WRZpACVwyx8cf++vf4KfvaNESSTUzZwbTGffuhdLSSx4ucilK4JI5YkeeX3EF9O4dXiwiiXTrBuPHB231wqUZlMAlM5w4Eb9kqBYtkVQVu8jJu++C1niQJlICl8wQO7d2xAgYpMXtJEX17QuXXx60tdSoNJESuKS/M2d83fMo9b4l1cX+jRYXa6lRaRIlcEl/sfWlBwyAkSPDjUekIUOGwLBhfrumJn78hkiSlMAlvVVUaMlQSU9XXx1sb9oEp06FF4ukJSVwSW9vv+2XaQR/bbGgINx4RJL1qU8FYzWqq9ULl0ZTApf0VVkJGzYE7WuuUe9b0ocZXHtt0C4qUi9cGkUJXNLX22/DuXN+u08f9b4l/YwaBQMH+u3qao1Il0ZRApf0VFkZf+1bvW9JR3V74Rs3wunToYUj6UUJXNLTunVB77t3b/W+JX2NHh3fC9e1cEmSErikn8pKn8CjrrkmqC8tkm7M/N9wVFGReuGSFH3rSfpZvz6+9z12bLjxiDTXmDG+hgH4mga6Fi5JUAKX9HL2bHzve+5c9b4l/dXthW/c6CsMilyCvvkkvaxb55M4QK9eMG5cuPGItJTLL4f+/f22euGSBCVwSR91e9+69i2ZpO6I9A0b1AuXS0rq28/MFprZbjMrMbNvJrg/38yeity/3syGR/bfbWZbY261ZjYpct/rkeeM3te3Jd+YZKD169X7lsymXrg0QoMJ3MxygJ8BNwEFwJ1mVnfOzr3AcefcKOBB4McAzrnHnXOTnHOTgC8Ce51zW2Med3f0fufckRZ4P5Kpzp71hVuidO1bMlGia+EVFeHFIyktmW/AaUCJc26Pc64KWAosrnPMYuDRyPYzwHyzi6pq3Ak82ZxgJYtt2BD0vnv2hPHjw41HpLVccQX06+e3q6rUC5d6JZPABwH7Y9qlkX0Jj3HOVQMngV51jvkCFyfwRyKnz7+dIOGLeOfOqfct2aNuL3zDBvXCJaFkvgUTJVbXmGPMbDpQ4ZzbEXP/3c658cDVkdsXE7642X1mVmRmRWVlZUmEKxln/XpfvAV873vChHDjEWltV17pV9cD9cKlXskk8FJgSEx7MHCwvmPMLBfoBpTH3L+EOr1v59yByH9PAU/gT9VfxDn3kHOu0DlX2KdPnyTClYxSURFfWlK9b8kGdUekr1+vlcrkIsl8E24ERpvZCDNrj0/Gy+ocswy4J7J9B7DSOecAzKwd8Dn8tXMi+3LNrHdkOw9YBOxApK7Vq+Orrqn3Ldniyivjq7O9/nqo4UjqaTCBR65p3w+8DOwCnnbO7TSzB8zslshhDwO9zKwE+DoQO9VsLlDqnNsTsy8feNnMioGtwAHg/zb73UhmOXEifr3v+fPV+5bsYQbXXx+0t2yBo0fDi0dSTm4yBznnlgPL6+z7Tsz2WXwvO9FjXwdm1Nl3BriqkbFKtlm5Empq/PaQIX50rkg2+dSnYORI2LMHamvh1VdhyZKwo5IUoe6MpKbDh2H79qC9YIHW+5bstGBBsP3uu7B/f/3HSlZRApfUtGIFuMhEhssvh6FDw41HJCwDBsTXPYj9tyFZTQlcUs+ePfDBB3677nVAkWw0bx7k5Pjtfftg9+5w45GUoAQuqcU538OImjwZNH1Qsl2PHlBYGLRfe81fE5espgQuqWXnTjh0yG/n5sbPhRXJZnPnQn6+3y4rg61bL328ZDwlcEkdNTW+ZxE1YwZ07RpePCKppHNnmD07aK9a5eeHS9ZSApfUUVQEx4/77Y4dYc6ccOMRSTUzZkCXLn771ClYty7ceCRUSuCSGs6dgzfeCNpz50KHDuHFI5KK2rePv6z01lta6CSLKYFLalizJvgi6t4dpk4NNx6RVDV5si8rDP6H7+rV4cYjoVECl/CdOhW/XOi8eX4Am4hcrF07X1Y4asMGX3ZYso4SuITvjTeCwTj9+8cXrRCRi11xhS8vDH7w58qV4cYjoVACl3AdPgybNgXt669XyVSRhtQtcFRcrBKrWUgJXMLjHLzwQlAW8lOf8jcRadiwYfEL/Lzwgm/9i70AABx4SURBVIq7ZBklcAnPtm1BryEnB266Sb1vkcZYuBDy8vz24cN+KqZkDSVwCUdlZXzJ1FmzgpG1IpKc7t3h6quD9sqVcPp0ePFIm1ICl3CsXAlnzvjtbt3iv4REJHmzZkGvXn777Nn4H8aS0ZTApe0dPBh/qm/hQl+gQkQaLzfXX36K2rYNPvoovHikzSiBS9tyDpYvDwaujRoVPxBHRBpv1CgoKAjay5drQFsWUAKXtrVlC5SW+m0NXBNpOTfeGAxo+/hjX+BFMpoSuLSdigp49dWgPWdOcO1ORJqnWze45pqgvWqVr3IoGUsJXNrOypXx9c612phIy5o5M75Ouga0ZTQlcGkbBw7EV1y76abgdJ+ItIycHLj55qBdXAx794YWjrQuJXBpfbW18RXXxozxNxFpeSNHwrhxQfuFF3y9dMk4SSVwM1toZrvNrMTMvpng/nwzeypy/3ozGx7ZP9zMKs1sa+T2i5jHXGVm2yOP+Q8zjWTKWJs3+6ljEEx50f9ukdZzww3B1MyyMli/Ptx4pFU0mMDNLAf4GXATUADcaWYFdQ67FzjunBsFPAj8OOa+D5xzkyK3r8bs/y/gPmB05Law6W9DUtaZM/Daa0F7zhzo0SO8eESyQdeucO21Qfv11+HkybCikVaSTA98GlDinNvjnKsClgKL6xyzGHg0sv0MMP9SPWozGwB0dc697ZxzwG+BWxsdvaQ25+CPf/RlU8En7tmzw41JJFtMnw59+/rtqir4wx+Cy1iSEZJJ4IOA2HXqSiP7Eh7jnKsGTgLR+UEjzGyLmb1hZlfHHF/awHMCYGb3mVmRmRWVlZUlEa6kjG3b4N13g/aiRRq4JtJWcnLgM58JLlft2QMbN4Ybk7SoZBJ4op503Z9x9R1zCBjqnJsMfB14wsy6JvmcfqdzDznnCp1zhX369EkiXEkJJ0/Ciy8G7alTtVSoSFsbMsTXSo9asQKOHQsvHmlRySTwUmBITHswcLC+Y8wsF+gGlDvnzjnnjgE45zYBHwBjIscPbuA5JV05B88/7+ehAvTsCQsWhBuTSLa67rrgVPr58/DccyqzmiGSSeAbgdFmNsLM2gNLgGV1jlkG3BPZvgNY6ZxzZtYnMggOMxuJH6y2xzl3CDhlZjMi18q/BPyhBd6PpIING+DDD/22GXz2s1qsRCQsublw223+lDr4UsZr1oQbk7SIBhN45Jr2/cDLwC7gaefcTjN7wMxuiRz2MNDLzErwp8qjU83mAsVmtg0/uO2rzrnyyH1/BfwKKMH3zGPOt0raOno0vvrT7Nn+NJ6IhKd///gyq6+/DocPhxaOtAxzaTQqsbCw0BXFLkMpqaW2Fh5+2FddA+jXD/7yL30PQETCVVsLv/51sJiQ/n2mBTPb5JwrTHSfKrFJy3nrrSB55+T403b6chBJDe3a+ctZsSuWvf56qCFJ8yiBS8s4dCj+y+Daa/0vfBFJHb16wfXXB+01a2D//vqPl5SmBC7NV10dP7J1yBAVbBFJVdOm+Xrp4GeMPPecL/QiaUcJXJpv1So4csRv5+XBrbf603UiknrMYPFiyM/37fJyLTuapvQtK82zdy+sXRu0Fyzwp+lEJHV16+YXFYrauBHefz+8eKRJlMCl6T75BH73u6C+8siRvuKaiKS+iRPhiiuC9u9/D8ePhxePNJoSuDRNdTU8/bRfbQygUyd/6lzLhIqkBzNfK71LF9+urISnnvLV2iQtKIFL07z4YjCftF07+Nzn/BKGIpI+OneGL3whqNJ2+DAsW6ZVy9KEErg03qZN/ha1YAGMGBFePCLSdIMHw803B+3t22H9+vDikaQpgUvjlJbC8uVBe/x4mDEjvHhEpPmuusrfol55JVjPQFKWErgk7/Rpf927psa3+/WDW27RdW+RTHDTTb43Dr6mwzPP+GWBJWUpgUtyamr8iPNPPvHtjh1hyZKgLKOIpLfcXPj85+Gyy3z7zBk/qK26Oty4pF5K4JKcV16Bjz7y22Zw++3Qo0e4MYlIy+ra1Q9IjRZiOngQXnhBg9pSlBK4NGzbtvhBLfPnw6hR4cUjIq1n2DBYuDBob9kCWgUyJSmBy6UdOgR//GPQLihQnXORTDd1qi/0EvXSS7BvX3jxSEJK4FK/8nJ4/PHgGlifPr6GsgatiWQ2M1i0CAYM8O2aGli6FMrKwo1L4iiBS2KffAK//a0feQ5+4YMlS4IFEEQks+Xl+SIvnTr5dkWF/05QudWUoQQuFztzxv9DPXHCt3Nz4a67tEiJSLbp3h3uvhvat/ftU6f8d8OpU+HGJYASuNR19iw89hgcPerbOTn+V/iwYeHGJSLhGDTI/4DPzfXt48d9Eq+oCDcuUQKXGOfPwxNP+IFr4K+D3XYbjB4dblwiEq7hw/0c8ej0srIy/0P/3LlQw8p2SuDi1dT4og2xI00/8xkYOza8mEQkdYwZ43/QRwexHjwITz6p1ctCpAQuvmzi738PJSXBvhtvhClTwotJRFLPuHF+dHrU3r2+QmO0vLK0qaQSuJktNLPdZlZiZt9McH++mT0VuX+9mQ2P7F9gZpvMbHvkv/NiHvN65Dm3Rm59W+pNSSM45+d579wZ7LvmGpg5M7yYRCR1XXWVX4Ew6r334LnnfEdA2lRuQweYWQ7wM2ABUApsNLNlzrl3Yg67FzjunBtlZkuAHwNfAI4Cn3HOHTSzccDLwKCYx93tnFOJn7A4By+/7CstRc2YAddeG1pIIpIGZs/217/ffNO3d+zwI9UXLQquk0urazCBA9OAEufcHgAzWwosBmIT+GLge5HtZ4Cfmpk552IyAzuBDmaW75zTyIew1dTAH/4AxcXBvkmT/KnzFCzUcvbsWcrKyjh79izVWlxBskheXh59+/ala9euYYcS77rr/KyVDRt8e/NmqKz018m1yFGbSCaBDwL2x7RLgen1HeOcqzazk0AvfA886nZgS53k/YiZ1QDPAj90ThXz28TZs37AWux6vwUFKbs06MmTJ/n444/p06cP/fv3Jzc3F0vBOEVamnOOyspKDhw4AJBaSdzML0FaVQVbt/p9u3b5KWZ33hkUgJFWk8y5jkTflHUT7SWPMbOx+NPq/0/M/Xc758YDV0duX0z44mb3mVmRmRWVqYxf833yCTzySHzyLiyEO+5I2VNfR48eZfDgwfTo0YO8vDwlb8kaZkanTp0YNGgQR44cCTuci5n58sqxY2b274eHH1bFtjaQzDd2KTAkpj0YOFjfMWaWC3QDyiPtwcBzwJeccx9EH+CcOxD57yngCfyp+os45x5yzhU65wr79OmTzHuS+hw5Ar/6FXz8cbBv/nz49KdTNnkDVFVV0bFjx7DDEAlNx44dOZ+q07XM/KW32Mtvx475JH6wbqqQlpTMt/ZGYLSZjTCz9sASYFmdY5YB90S27wBWOuecmXUHXgC+5ZxbEz3YzHLNrHdkOw9YBOxo3luRS/rwQ/j1r30PHHzC/uxn4eqrU/K0eV3qdUs2S4u//5kz/Zm8aMW206fhN7+B998PNaxM1mACd85VA/fjR5DvAp52zu00swfM7JbIYQ8DvcysBPg6EJ1qdj8wCvh2neli+cDLZlYMbAUOAP+3Jd+YxNi+3VdNOnvWt/PzfX3j2OUCRUSaa+xY+OIXIXrGrKrKF3vZvDncuDKUpdO4scLCQlekheWT5xysXQsrVgT7unTxybt///DiaqRdu3Zx5ZVXhh2GSKjS6t9BWZlfiji6IBL4+hLXXpsWZ/xSiZltcs4VJrovdS98SvNUVMDTT8cn7z594C/+Iq2St4ikoT594N57g/XEAd54w6+1EF2iWJpNCTwTffgh/OIXfkpH1PDh8JWvQLduoYUlIlmkSxf48pdh1Khg3/vvw3/9V3zZZmkyJfBMUlMDr73m52FGB6sBTJsGf/ZnwXUpkRinT5/me9/7HjfffDN9+vTBzPje977X6Oepqqri29/+NkOHDqVDhw5MmDCBJ598suUDDsHGjRv527/9W8aPH89ll13GoEGDWLRoEY29pJfJn1FC+fl+TvisWcG+M2f8mJyXXgIVZWoWJfBMUV7uR5mvXu2vfYMvpHDXXXDzzcHIUJE6jh49yve//32Ki4uZ0owFbO69915+9KMfsXjxYv7zP/+TQYMGcdddd/H444+3YLTh+PGPf8xTTz3F3LlzefDBB/m7v/s7du3axfTp01m+fHnSz5PJn1G9cnLghhv84LbLLgv2r1vnp7WqvkeTaRBbunPOl0N94QU/4jNq5Eg/TaxLl/BiayFpNXgnDZ07d45jx44xcOBASktLGTJkCN/97ncb1QvftGkThYWFcY9zzjF37lxKSkrYt28feWlcXnPt2rUUFhbSvn37C/uOHTtGQUEBgwYNYnMSo6yb+xllxL+DM2d8Cef33gv25eXBwoV+9UMNcLuIBrFlqspKvwzoc88FyTv2124GJO9sMnnyZG666aaL9v/93/89nTt3praVVnvKz89n4MCBzXqOp59+GjPjb/7mby7sMzP++q//msOHD/NmdNGLFvD4448zc+ZMunTpgplddCstLW2x14qaNWtWXPIG6NWrF9deey3vvPNOPY+K15afUcrq3NmfUo89K3j+vF8R8emnfYKXpOm8ajqqrYWiIli1yifxqF694PbboZlfxtL2qqur2bVrFzfeeONF923bto1x48bRLkG1vPPnz3Py5MmkXqNLly7k5+c3O9ZENm/ezPDhw6lbLXHatGkX7p8/f36zX+e73/0uDzzwADNnzuSBBx7g1KlT/OQnP+H48ePMmTOHgQMHMnjw4LjHtOZndPDgQXr16pXUsW31GaU8Mz8uZ9gwePZZXyES/KDbPXv8dLPp031nRC5JCTzdfPCBH/xR97rR5Ml+YYE6vQRJD7t37+bcuXNMmjTpovuKi4u5/fbbEz5uzZo1XHfddUm9xiOPPMKXv/zl5oRZr0OHDjEgdspQRHTfwRYoqVlUVMQPfvADPve5z/Hkk0+SE/mCHzduHLfffjtLliyJ691GtdZntHr1atasWcPXvva1pI5vi88orfTrB3/5l36qa3RFs3Pn4JVXYNMmX5p19GidVr8EJfB0ceyY/8PevTt+f48e/g/9iivCiStMTRgp3eqaGNP27dsBmFinOt7+/fspLy+/aH/UxIkTWRE71/8Sxo4d26TYklFZWUnfvn0v2t+hQ4cL9zfXgw8+SH5+Pj//+c8vJG+AayPr17/77rsJH9can9GhQ4e48847GTp0KN/5zneSekxbfEZpJy/Pn04fMwZefNF/z4H/7xNP+CloN97o55XLRZTAU93Zs/Dmm7B+vZ8mFtW+PcydCzNmaIR5BiguLqZjx46MGTMmbv+2bduAixN7VI8ePbj++utbPb6GdOzYkXPnzl20/2ykfG9zF6NxzvHSSy8xb948evfuHXdf9HXrW2qzpT+jkydPcvPNN3P69GlWr15NtyRrK7T2Z5TWRo2Cv/5r3xN/442g7HNJiT+tPnWqr+KWzZ9RAvrmT1VVVb5+8OrVFw/smDwZ5s3TILUMUlxczLhx4+J6lgBbt27FzJgwYULCx1VVVVFeXp7Ua3Tr1q3VksSAAQMoSVCc49ChQwDNHiR34MABysvLGT9+/EX3bdq0CYCCgoKEj23Jz6iiooJFixaxe/duXnnllYTx1Ke1P6O0l5PjF0SZMMGP79m0yc+yqa31HZjiYpgzB666CiJnLbKdEniqOXXK/7EWFQW/QqOGDPHXubP9H3pUKp5Cb6Li4uILp4JjvfTSSwwfPrze3uXatWtT4hr4lClTeO211zhy5EjcaeL169dfuL85PokUJqo7Ehxg6dKldOjQIeEIfmi5z6iqqorbbruNdevW8fzzzzNnzpzkgo9o7c8oY3TuDIsWQWGhH++zd6/fX1npr5e/+aZP4tOnZ31lSSXwVHHkiF94ZPv2+FPlAF27+qlhY8dqQEcGOnnyJPv37+fj2HXagRdeeIE1a9awePHieh/b1tfAz58/zwcffEC3bt3iBmTdcccd/Mu//As///nP4+Y4/+IXv6Bfv37MnTu3Wa87ePBg2rVrx2uvvcb3v//9C8trLl++nCeeeIJvfOMb9OzZM+FjW+Izqqmp4a677mLFihX893//N5/+9KfrfY6wPqOM078/3HMPvPuuH/9z/Ljff+6c/65ctw7GjfO99gSDA7OBEniYnPN1y9euTVwbuGdP/8c5aZIf7CEZqbi4GPCjpb/yla8wefJkduzYwfPPPw/Anj17WLZsGbfccstFj22p67s//elPOXHixIWe7ptvvskPf/hDAL74xS8ybNgwwJ/KvvLKK7nnnnv4zW9+c+HxU6dO5c477+QHP/gB5eXlTJgwgd///vesXr2aRx99NK5AiZlxzTXX8PrrrycdX9euXfn85z/P0qVLue2221iwYAHFxcX86le/4rrrrrsQayIt8Rn9wz/8A88++ywLFiygtraWxx57LO7+z372s3Tu3Blomc9IIszgyiv9aPTiYv9defSov6+21u8rLvaFq2bO9NfSs6iTowQehuPHfU97+/bEZQSHDPG1gy+/HBLM/ZXMEk3gv/vd7/jGN77BE088weTJk3nppZf4+te/zrZt2zjdyis4/du//RsfffTRhfaqVatYtWoVAHPmzLmQwC/lkUceYcSIEfz2t7/ll7/8JWPGjOGxxx7j7rvvvnBM9H0kmk7VkIceeoju3bvz7LPPsnz5cj71qU/xox/9iK997Wutnvy2bNkCwIoVKxL25j/88MMLCfxSkvmMJIHcXF+pbfJkvyDK2rXBqXXwA9327PG1MMaP97ck5+enM5VSbSunT8POnT5pJ6oUZeangs2a5RO4XJARJSQv4atf/Sp/+tOfWqWCWKpZvnw5ixYtYtu2bY0aACaZ/++g0Q4e9In8nXd8b7yugQN9Ih871l+GTFOXKqWqHnhrOnvWVxfavt2fKk/0Yykvz58inznTnzKXrFNcXFzvCOpMs2rVKpYsWaLkLc03cCDccQecOOGvh2/eHL8exMGD/vbKK77q2/jxUFCQUVPRlMBbUk0NHDgQnM4pLU38y7BdO3+tZvx4f5pc1dOylnOOHTt2cO+994YdSpv413/917BDkEzTvbtfDGX+fL9Iyvbt/jR7dDCwc/50+9698Kc/waBB/pr5yJEweHBa19FI38hTgXN+9Hg0YX/0UfwvwLqGD/ejJgsK/FKfkvX27t3LqVOnsqYHLtJq8vL86fKxY/2Us127YMeO+LOfzvmOVWmpn46Wl+d759GE3q9fWg2CUwJPlnN+jvahQ/60zKFDvrfd0Oo5Awf6pD1uXFpfh5HWMWLECNJpHIpIWujY0Q96mzLFf2/v3OmT+YED8Zcyz5/3M4Cis4A6dfI99AED/Hf3gAH+eztFk7oSeCI1Nf66SllZkKwPHfID0RrSvXvwa27ECF+UQEREwtGliy85PWOG75l/+GFw1rRuhb6KCn/6/f33g32dO/tEHk3qffr4NShSYLW07E3gtbVw8qQvml9eHv/fEycSX7tOpGNHn6ijSbtHj5T9tSYiktU6dvSXMKOXrI4fj0/oFRUXP+bMmfheOvjv+O7d/VS1nj39f6Pb3bu32fTfpBK4mS0E/h3IAX7lnPvfde7PB34LXAUcA77gnNsbue9bwL1ADfC3zrmXk3nOVrVmDaxceXHFs4a0bx/8Eov+GuvVS3O1RUTSUY8e/jZlij+1fuxY/GXSQ4d85be6nPPJP1odLlZOjq/ZnmT53uZoMIGbWQ7wM2ABUApsNLNlzrl3Yg67FzjunBtlZkuAHwNfMLMCYAkwFhgIvGpm0eWWGnrO1tOhQ8PJu2tXn5z79w+uhfTqpd61iEgmMoPevf0tOs3ROX9mNprMDx3ySf6TTxJPCwafW9posZVkeuDTgBLn3B4AM1sKLAZik+1i4HuR7WeAn5ovVrwYWOqcOwd8aGYlkecjiedsPdEKPZddlvgUSM+eKl2aYpxzF+pfi2QbDXQMiVmQG8aNC/afP+9734kuwZ461WZV4JJJ4IOA/THtUmB6fcc456rN7CTQK7J/XZ3HDopsN/ScrWfIEPjWtyA/v81eUpquffv2VFZW0klT7yRLVVZWqlZ6KsnLg759/a2uqqo2u6yazKsk6vbU/TlY3zGN3X/xi5vdZ2ZFZlZUlqhueFPk5Ch5p5HevXtTWlpKeXk558+fV29EsoZzjoqKCg4cOBC3BKmksPbt26w4TDKvUgrEFuceDBys55hSM8sFugHlDTy2oecEwDn3EPAQ+FroScQrGaZbt27k5+dTVlbGsWPHqK6uDjskkTaTl5dHv3796l0TXrJXMgl8IzDazEYAB/CD0u6qc8wy4B7gbeAOYKVzzpnZMuAJM/s/+EFso4EN+B54Q88pckGHDh0YokVeREQuaDCBR65p3w+8jJ/y9Wvn3E4zewAocs4tAx4G/jsySK0cn5CJHPc0fnBaNfA3zrkagETP2fJvT0REJDNpOVEREZEUdanlRFWBREREJA0pgYuIiKQhJXAREZE0pAQuIiKShpTARURE0lBajUI3szLgo7DjaAW9gaNhB5FC9HkE9FnE0+cR0GcRL1M/j2HOuT6J7kirBJ6pzKyovmkC2UifR0CfRTx9HgF9FvGy8fPQKXQREZE0pAQuIiKShpTAU8NDYQeQYvR5BPRZxNPnEdBnES/rPg9dAxcREUlD6oGLiIikISXwFGNm/2Bmzsx6hx1LWMzsX83sXTMrNrPnzKx72DGFwcwWmtluMysxs2+GHU9YzGyIma0ys11mttPM/i7smFKBmeWY2RYz+1PYsYTNzLqb2TOR741dZjYz7JjaghJ4CjGzIcACYF/YsYRsBTDOOTcBeA/4VsjxtDkzywF+BtwEFAB3mllBuFGFphr4hnPuSmAG8DdZ/FnE+jtgV9hBpIh/B15yzl0BTCRLPhcl8NTyIPD/Alk9MME594pzrjrSXAcMDjOekEwDSpxze5xzVcBSYHHIMYXCOXfIObc5sn0K/+U8KNyowmVmg4FPA78KO5awmVlXYC7wMIBzrso5dyLcqNqGEniKMLNbgAPOuW1hx5JivgK8GHYQIRgE7I9pl5LlSQvAzIYDk4H14UYSup/gf+zXhh1IChgJlAGPRC4p/MrMOocdVFvIDTuAbGJmrwL9E9z1z8D/BG5o24jCc6nPwjn3h8gx/4w/ffp4W8aWIizBvqw+M2NmlwHPAn/vnPsk7HjCYmaLgCPOuU1mdm3Y8aSAXGAK8D+cc+vN7N+BbwLfDjes1qcE3oacc9cn2m9m44ERwDYzA3/KeLOZTXPOHW7DENtMfZ9FlJndAywC5rvsnOtYCgyJaQ8GDoYUS+jMLA+fvB93zv0+7HhCNhu4xcxuBjoAXc3sMefcn4UcV1hKgVLnXPSszDP4BJ7xNA88BZnZXqDQOZeJhfkbZGYLgf8DXOOcKws7njCYWS5+AN984ACwEbjLObcz1MBCYP5X7aNAuXPu78OOJ5VEeuD/4JxbFHYsYTKz1cBfOOd2m9n3gM7OuX8MOaxWpx64pKKfAvnAisgZiXXOua+GG1Lbcs5Vm9n9wMtADvDrbEzeEbOBLwLbzWxrZN//dM4tDzEmSS3/A3jczNoDe4A/DzmeNqEeuIiISBrSKHQREZE0pAQuIiKShpTARURE0pASuIiISBpSAhcREUlDSuAiIiJpSAlcREQkDSmBi4iIpKH/H/U/aH3JJiOiAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(-5, 7)\n", "mu = 1\n", "sigma = 2\n", "\n", "y = sp.stats.norm.pdf(x, loc = mu, scale = sigma)\n", "\n", "fig, ax = plt.subplots(figsize = (8, 5))\n", "ax.plot(x, y, lw = 3, color = 'r', alpha = .5, \n", " label = r'$\\mu = %.1f,\\ \\sigma = %.1f$'%(mu, sigma))\n", "ax.legend(loc ='best', fontsize= 17)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's all we need to know about univariate distribution! \n", "\n", "For more details, check out my notebook of basic statistics." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Multivariate Normal Distribution (MND)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The PDF of a multivariate normal distribution is \n", "\n", "$$\n", "p(\\mathbf{x} ; \\mu, \\Sigma)=\\frac{1}{(2 \\pi)^{n / 2}|\\Sigma|^{1 / 2}} \\exp \\left(-\\frac{1}{2}(x-\\mu)^{T} \\Sigma^{-1}(x-\\mu)\\right) \\tag{1}\\label{1}\n", "$$\n", "\n", "Before we analyze the PDF, we should know some basics of random vectors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Expectation and Covariance Matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a random $n$-vector $\\mathbf{x}$, its expection is defined as\n", "\n", "$$\n", "E(\\mathbf{x}) = \n", "\\left[\n", "\\begin{matrix}\n", "E(x_1)\\\\E(x_2)\\\\ \\vdots \\\\E(x_n)\n", "\\end{matrix}\n", "\\right]\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The variance of $\\mathbf{x}$ is a covariance matrix, denoted as\n", "\n", "$$\n", "\\begin{aligned}\n", "\\Sigma_{\\mathbf{x} \\mathbf{x}}=\\operatorname{Var}(\\mathbf{x}) &=E\\left[(\\mathbf{x}-\\mu)(\\mathbf{x}-\\mu)^T\\right]\\\\\n", "&=\\left[\\begin{array}{cccc}\n", "\\operatorname{Var}\\left(x_{1}\\right) & \\operatorname{Cov}\\left(x_{1}, x_{2}\\right) & \\dots & \\operatorname{Cov}\\left(x_{1}, x_{n}\\right) \\\\\n", "\\operatorname{Cov}\\left(x_{2}, x_{1}\\right) & \\operatorname{Var}\\left(x_{2}\\right) & \\dots & \\operatorname{Cov}\\left(x_{2}, x_{n}\\right) \\\\\n", "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n", "\\operatorname{Cov}\\left(x_{n}, x_{1}\\right) & \\operatorname{Cov}\\left(x_{n}, x_{2}\\right) & \\dots & \\operatorname{Var}\\left(x_{n}\\right)\n", "\\end{array}\\right]\n", "\\end{aligned}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Linear Combination of Normal Distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you have guessed, the linear combination of a sequence of normal distribution will be a normal distribution as well. Let's say we have another random vector $\\pmb{z}$\n", "\n", "$$\n", "\\mathbf{z}=\n", "\\left[\n", "\\begin{matrix}\n", "z_1\\\\z_2\\\\ \\vdots \\\\z_n\n", "\\end{matrix}\n", "\\right]\n", "$$\n", "\n", "where $z_i\\sim iid(0,\\ \\sigma^2)$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For any matrix $A_{n\\times n}$ of full rank, random normal vector $\\mathbf{x}$ can be write as\n", "\n", "$$\n", "\\mathbf{x} = A\\mathbf{z}\n", "$$\n", "\n", "It simply states that each $x_i, i=(1,2,...,n) $ is a linear combination of $\\mathbf{z}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $\\mathbf{\\mu} = \\mathbf{0}$, the variance of $\\mathbf{x}$ is \n", "\n", "$$\n", "\\operatorname{Var}(\\mathbf{x})=E\\left(\\mathbf{x} \\mathbf{x}^{T}\\right)={A }E\\left(\\mathbf{z} \\mathbf{z}^{T}\\right) {A}^{T}={A} \\mathbf{I} {A}^{T}={A} {A}^{T}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Covariance matrix ${AA}^T$ is positive semi-definite matrix, let $\\mathbf{x}$ be an $n$-vector, then\n", "\n", "$$\n", "\\mathbf{x}^T{AA}^T\\mathbf{x}= ({A}^T\\mathbf{x})^T({A}^T\\mathbf{x}) = \\|{A}^T\\mathbf{x}\\|^2 \\geq 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It also means that all eigenvalues of ${AA}^T$ are non-negative as well, refer to Chapter 17, the section of positive definite matrix." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inverse and Positive Definite" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $ {A}_{n\\times n}$ is positive definite and symmetric, then all the eigenvalues are larger than $0$. But what if there is $0$ eigenvalue, then we get\n", "\n", "$$\n", " {A\\mathbf{x}} = {0}\n", "$$\n", "\n", "since $ {\\mathbf{x}}$ is a nontrivial solution(eigenvector is always nonzero), therefore $ {A}$ must be non-invertible.\n", "\n", "Thus, if $ {A}$ is positive definite, $ {A}$ does not have $ {0}$ eigenvector which means it is invertible." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Inverse and Symmetry" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The same $ {A}$ as in last section, we have $ {A} {A}^{-1}=\\mathbf{I}$. Taking the transpose, we get\n", "\n", "$$\n", "( {A}^{-1})^T {A} = ( {A}^{T})^{-1} {A}= {A}^{-1} {A}=\\mathbf{I}\n", "$$\n", "\n", "We can see that $( {A}^{-1})^T= {A}^{-1}$, $ {A}^{-1}$ is also a symmetric matrix." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can show that $ A^{-1}$ is also positive definite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "A \\mathbf{v}=\\lambda \\mathbf{v} \\Longrightarrow A^{-1} A \\mathbf{v}=\\lambda A^{-1} \\mathbf{v} \\Longrightarrow\\frac{1}{\\lambda} \\mathbf{v} = A^{-1}\\mathbf{v}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have proved that if $ {A}$ has eigenvalue $\\lambda$, then $ {A}^{-1}$ has $\\frac{1}{\\lambda}$ as its eigenvalue.\n", "\n", "If $\\lambda>0$, certainly $\\frac{1}{\\lambda}>0$, thus $ {A}^{-1}$ is also positive definite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bivariate Normal Distribution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the PDF of MND, the argument of exponential function is $-\\frac{1}{2}(\\mathbf{x}-\\mathbf{\\mu})^{T} \\Sigma^{-1}(\\mathbf{x}-\\mathbf\\mu)$, which is a quadratic form. \n", "\n", "$\\Sigma$ is symmetric positive semi-definite, so is $\\Sigma^{-1}$ for any vector $\\mathbf{x} = \\mathbf{\\mu}$. \n", "\n", "With a minus sign in front, we get negative semi-definite quadratic form\n", "\n", "$$\n", "-\\frac{1}{2}(\\mathbf{x}-\\mathbf{\\mu})^T\\Sigma^{-1}( \\mathbf{x}-{\\mu})\\leq 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we define a simple bivariate case, the quadratic form is \n", "\n", "\\begin{align}\n", "-\\frac{1}{2}\\left[\\begin{array}{l}\n", "x_{1}-\\mu_{1} \\\\\n", "x_{2}-\\mu_{2}\n", "\\end{array}\\right]^{T}\\left[\\begin{array}{cc}\n", "\\sigma_{1}^{2} & 0 \\\\\n", "0 & \\sigma_{2}^{2}\n", "\\end{array}\\right]^{-1}\\left[\\begin{array}{l}\n", "x_{1}-\\mu_{1} \\\\\n", "x_{2}-\\mu_{2}\n", "\\end{array}\\right]&=\n", "-\\frac{1}{2}\\left[\\begin{array}{l}\n", "x_{1}-\\mu_{1} \\\\\n", "x_{2}-\\mu_{2}\n", "\\end{array}\\right]^{T}\\left[\\begin{array}{cc}\n", "\\frac{1}{\\sigma_{1}^{2}} & 0 \\\\\n", "0 & \\frac{1}{\\sigma_{2}^{2}}\n", "\\end{array}\\right]\\left[\\begin{array}{l}\n", "x_{1}-\\mu_{1} \\\\\n", "x_{2}-\\mu_{2}\n", "\\end{array}\\right]\\\\\n", "&=-\\frac{1}{2}\\left[\\begin{array}{l}\n", "x_{1}-\\mu_{1} \\\\\n", "x_{2}-\\mu_{2}\n", "\\end{array}\\right]^{T}\\left[\\begin{array}{l}\n", "\\frac{1}{\\sigma_{1}^{2}}\\left(x_{1}-\\mu_{1}\\right) \\\\\n", "\\frac{1}{\\sigma_{2}^{2}}\\left(x_{2}-\\mu_{2}\\right)\n", "\\end{array}\\right]\\\\\n", "& = -\\frac{1}{2 \\sigma_{1}^{2}}\\left(x_{1}-\\mu_{1}\\right)^{2}-\\frac{1}{2 \\sigma_{2}^{2}}\\left(x_{2}-\\mu_{2}\\right)^{2}\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Further, we give value to $\\sigma$ and $\\mu$\n", "\n", "$$\n", "\\sigma_1 = 2,\\ \\sigma_2 = 3,\\ \\mu_1 = 0,\\ \\mu_2 = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can visualize quadratic form and its exponential" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " // select the cell after this one\n", " var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n", " IPython.notebook.select(index + 1);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib notebook\n", "mu_1 = 0\n", "sigma_1 = 2\n", "\n", "mu_2 = 0\n", "sigma_2 = 3\n", "\n", "#Create grid and multivariate normal\n", "x = np.linspace(-5,5,30)\n", "y = np.linspace(-5,5,30)\n", "X, Y = np.meshgrid(x,y)\n", "pos = np.empty(X.shape + (2,))\n", "pos[:, :, 0] = X; pos[:, :, 1] = Y \n", "norm = sp.stats.multivariate_normal([mu_1, mu_2], [[sigma_1, 0], [0, sigma_2]]) # frozen \n", "\n", "#Make a 3D plot\n", "fig = plt.figure(figsize = (8, 5))\n", "ax = fig.gca(projection='3d')\n", "ax.plot_surface(X, Y, norm.pdf(pos),cmap='viridis',linewidth=0)\n", "ax.set_xlabel('X axis')\n", "ax.set_ylabel('Y axis')\n", "ax.set_zlabel('Z axis')\n", "\n", "ax.set_title('Bivariate Normal Distribution, $\\sigma_1 = 2$, $\\sigma_2 = 3$, $\\mu_1 = 0$, $\\mu_2 = 0$')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we have expanded the bivariate quadratic form, back to the MVN PDF $\\eqref{1}$, we have\n", "\n", "\n", "\\begin{aligned}\n", " p(\\mathbf{x};\\mathbf{\\mu}, \\Sigma)\n", "&=\\frac{1}{2 \\pi \\sigma_{1} \\sigma_{2}} \\exp \\left(-\\frac{1}{2 \\sigma_{1}^{2}}\\left(x_{1}-\\mu_{1}\\right)^{2}-\\frac{1}{2 \\sigma_{2}^{2}}\\left(x_{2}-\\mu_{2}\\right)^{2}\\right)\\\\\n", "&=\\frac{1}{\\sqrt{2 \\pi} \\sigma_{1}} \\exp \\left(-\\frac{1}{2 \\sigma_{1}^{2}}\\left(x_{1}-\\mu_{1}\\right)^{2}\\right) \\cdot \\frac{1}{\\sqrt{2 \\pi} \\sigma_{2}} \\exp \\left(-\\frac{1}{2 \\sigma_{2}^{2}}\\left(x_{2}-\\mu_{2}\\right)^{2}\\right)\n", "\\end{aligned}\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We find that bivariate normal distribution can be decomposed into a product of two single variate normal distribution!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Covariance Matrix" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Covariance matrix is the most important factor that shapes the distribution. Use Scipy multivariate normal random generator, we can learn some intuition of the covariance matrix." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support. ' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
');\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('