{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Density Estimation"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Preliminaries\n",
"\n",
"- Goal \n",
" - Simple maximum likelihood estimates for Gaussian and categorical distributions\n",
"- Materials \n",
" - Mandatory\n",
" - These lecture notes\n",
" - Optional\n",
" - Bishop pp. 67-70, 74-76, 93-94 "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Why Density Estimation?\n",
"\n",
"Density estimation relates to building a model $p(x|\\theta)$ from observations $D=\\{x_1,\\dotsc,x_N\\}$. \n",
"\n",
"Why is this interesting? Some examples:\n",
"\n",
"- **Outlier detection**. Suppose $D=\\{x_n\\}$ are benign mammogram images. Build $p(x | \\theta)$ from $D$. Then low value for $p(x^\\prime | \\theta)$ indicates that $x^\\prime$ is a risky mammogram."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- **Compression**. Code a new data item based on **entropy**, which is a functional of $p(x|\\theta)$: \n",
"$$\n",
"H[p] = -\\sum_x p(x | \\theta)\\log p(x |\\theta)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- **Classification**. Let $p(x | \\theta_1)$ be a model of attributes $x$ for credit-card holders that paid on time and $p(x | \\theta_2)$ for clients that defaulted on payments. Then, assign a potential new client $x^\\prime$ to either class based on the relative probability of $p(x^\\prime | \\theta_1)$ vs. $p(x^\\prime|\\theta_2)$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
"### Example Problem\n",
"\n",
"\n",
"Consider a set of observations $D=\\{x_1,…,x_N\\}$ in the 2-dimensional plane (see Figure). All observations were generated by the same process. We now draw an extra observation $x_\\bullet = (a,b)$ from the same data generating process. What is the probability that $x_\\bullet$ lies within the shaded rectangle $S$?\n",
"\n",
"\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG2CAYAAAB20iz+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X2QVFV+//HPTEO3IoICq4sFCD4wuqHcCpBkxyAMsD4QAsgvcePulms2bqpUQBGTUowz03RjmFIJvy3XEHwozB9RkvxcFGpYdJQnzYYsmKV8nqy6ZljRXcC1R0C7l577+wO77enpnumevt33nnPfr6qptXu6e07fmeV++tzv+Z46x3EcAQAAWKje6wEAAABUC0EHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFjL2KCzZs0a1dXVafny5V4PBQAA+JSRQWffvn165JFHdNlll3k9FAAA4GPGBZ1jx47pu9/9rh599FGdffbZXg8HAAD42BCvB1CuJUuWaP78+frmN7+p1atX9/vYZDKpZDKZvd3T06OPP/5Yo0ePVl1dXbWHCgAAXOA4jj799FOdd955qq8vb47GqKCzadMm/fd//7f27dtX0uPXrFmjVatWVXlUAACgFg4ePKhx48aV9Rxjgs7Bgwd1++236/nnn9dpp51W0nNWrlypFStWZG8nEglNmDBBBw8e1IgRI6o1VAAA4KLu7m6NHz9eZ555ZtnPrXMcx6nCmFz3zDPPaPHixQqFQtn70um06urqVF9fr2Qy2et7hXR3d2vkyJFKJBIEHQAADFHJ+duYGZ25c+fqtdde63Xf97//fV1yySW66667Bgw5AAAgeIwJOmeeeaamTJnS674zzjhDo0eP7nM/AACAZODycgAAgFIZM6NTyK5du7weAgAA8DFmdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAeBL0WhU8Xi84Pfi8bii0WhtBwTASAQdAL4UCoXU0tLSJ+zE43G1tLSw7QuAkhjdMBCAvZqbmyVJLS0t2duZkBOLxbLfB4D+GLN7uRvYvRwwTybchMNhpVIpQg4QQJWcvwk6AHwvEokolUopHA4rmUx6PRwANVbJ+ZsaHQC+Fo/HsyEnlUoVLVAGgEIIOgB8K7cmJ5lMKhaLFSxQBoBiKEYG4EuFCo8LFSgDQH8IOgB8KZ1OFyw8ztxOp9NeDAuAYShGBgAAvkYxMgAAQAEEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AABGiEajisfjBb8Xj8cVjUZrOyAYgaADADBCKBRSS0tLn7ATj8fV0tKiUCjk0cjgZ0O8HgAAAKVobm6WJLW0tGRvZ0JOLBbLfh/IVec4juP1IGqlu7tbI0eOVCKR0IgRI7weDgBgEDLhJhwOK5VKEXICoJLztzGXrtavX6/LLrtMI0aM0IgRI9TY2Kif/OQnXg8LAFBjzc3N2ZATDocJOeiXMUFn3Lhxamtr0/79+7V//37NmTNHixYt0htvvOH10AAANRSPx7MhJ5VKFS1QBiSDgs6CBQv0J3/yJ5o8ebImT56s++67T8OHD9fevXu9HhoAoEZya3KSyaRisVjBAmUgw8hi5HQ6rX//93/X8ePH1djYWPRxyWRSyWQye7u7u7sWwwMAVEGhwuNCBcpALqOCzmuvvabGxkZ9/vnnGj58uDZv3qyvfe1rRR+/Zs0arVq1qoYjBACzRKNRhUKhggEhHo8rnU77pj9NOp0uWHicuZ1Op70YFnzOqFVXqVRKXV1d+uSTT/T000/rscce0+7du4uGnUIzOuPHj2fVFQB8odjybJZtw08qWXVl1IxOOBzWRRddJEmaPn269u3bpx/+8IfasGFDwcdHIhFFIpFaDhEAjEJvGtjOqKCTz3GcXjM2AIDy5Yad1atX05sGVjHm0tU999yjefPmafz48fr000+1adMmtbW1afv27bryyitLeg0aBgJAcZFIJLtsmw+R8JNANAz89a9/rRtuuEENDQ2aO3eu/uu//quskAOgNth40Uz0poGtjAk6jz/+uN5//30lk0n95je/0QsvvEDIAXyIjRfNQ28a2MzoGh0A/kNxq1lK7U1j0jJ0IBdBB4DrKG41R6m9aTIzdbnfk3oHJcCPjClGdgPFyEBtUdxql/zZH2bqUCuBKEYGYBaKW+3T3Nycrd+JRCL9hhyK0uEXBB0ArqO41V7Nzc3Z8BoOh4vO5FCUDt9wAiSRSDiSnEQi4fVQAGvFYjFHkhOLxUq6H2bJ/B7D4fCAv8/83zl/AxisSs7fFCMDcBUbL9qrWI2OVHjXcIrS4QcUIwMABlTJ5p8UpaNSgdnUEwDgjcHO1BUqSmdGB7XEjA4AoCpYjg63MKMDAPCVUjsuA9VG0AEAuI6idPgFl64AAICv0RkZgLHooIta4W8tmAg6ADxFB13UCn9rwUSNDgBPFSpQZXUOqoG/tWCiRgeAL2ROOJl+K5x4UC38rZmnkvM3QQeAb9BBF7XC35pZKEYGYLxCHXSBauBvLVgIOgA8l1snkUwmFYvFChaNApXiby14KEYG4Ck66KJW+FsLJoIOAE/RQbe2otGoQqFQwRN6PB5XOp22tp8Mf2vBRDEyAARIseXULLOGn7GpJwCgJPSSQdAwowMAAUQvGZiEPjolIugAwJfoJQNT0EcHAFAWeskgKAg6ABAw9JJBkFCMDAABQi8ZBA1BBwAChF4yCBqKkQEAgK9RjAwA6CUajRatuYnH49Z2PwbyEXQAwAfcDiahUKhggXGmRicUCg12qK4hjKEWCDoA4ANuB5Pm5uY+q6n81gHZhDAGCzgBkkgkHElOIpHweigA0EcsFnMkObFYrODtSl4zHA5X/FrVUI33DPtUcv6mGBkAfKQaWzP4vQMy21FgIGwBUSKCDgATuBlMTAkRfg9j8BarrgDAEm5uzWBKB2S2o0BVuXwZzdeo0QHgZ27WqxR7rt9qYKjRQSkqOX/TGRkAfMDtrRlM6IDMdhSoBYIOAPiAm8EkGo0qFAoVDAnxeFzpdNoXPWpMCGMwH8XIAGCZYv1y/NZHByhVJedvZnQAuM7LGQVTZjOqqdDlH5tDDr9z9IdVVwBc51bH28FsEUC33VNyOyNHIhFrQ47E7xwDcL002sdYdQXUjhuraQa7cijz/aamJicWixV8fCwWc1pbW8t/Y4bJdEQOh8OuvWZra2u/x96L48rqLbtVcv4m6ACoGje2HxjsCSzzuMxXfsgJwkmwWts/+HXput+3u8DgEXRKRNABas+NGYXBnsAyjw/iJ/1qz3D4dQalGjNY8B5Bp0QEHaC23PyEXe4JLP9nB+mTfq1mXPw2g+K38cA9BJ0SEXSA2qlGl99ST2DFfnZQPunXsobGLzMofpth8mMdk8kIOiUi6AC14eaMQrknsGKFx8XqdTB4fplB8WPNkB/HZLJABJ2///u/d6ZPn+4MHz7c+cpXvuIsWrTIefvtt8t6DYKO9/iUEwxu/Z4Hc7LI/9m5j43FYk5TUxMnGhf4aQbFr/+u+OkYmS4QQefqq692Nm7c6Lz++uvOgQMHnPnz5zsTJkxwjh07VvJrEHS8x6cclKPSExh/b9XBcS2dX2a9TFfJ+dvYLSAOHz6sc845R7t379bMmTNLeg5bQPhDfodWmzu2wlu17JgbpO68QXqvbohEIkqlUgqHw0omk14Px0iB3AIikUhIkkaNGlX0MclkstcfVXd3d9XHhYHltqdfvXq1UqkUIQeuSCQS+t3vfpe9vXTpUknSkSNH+jz2lltuKfq9wUgmk2pra9OJEyd05513Zu9fu3at2tradPfdd7v2s7xWy+NqunXr1mVDTiqVUjwe59+6GjNyRsdxHC1atEi//e1v9dJLLxV9XDQa1apVq/rcz4yOP/ApB25KJBLasGFD9kOQF15++WXt2bNHM2fO1IwZM/rcRrBkfv/33HOP7rvvPmavKxC4GZ2lS5fq1Vdf1csvv9zv41auXKkVK1Zkb3d3d2v8+PHVHh5KEI/H+ZQDV/3ud79TIpHQ6aefrmHDhnkyhkWLFmnYsGHavn27/uM//kPpdFrXXHONrrzyypKe/9xzz6m+vr7g4zs6OtTT06Orr77a7WGjCjo6OrIh94477pDUezb74MGDOnz4sF555RX9+te/1llnnaULLrhAl19+udauXevl0K1jXNBZtmyZtmzZoj179mjcuHH9PjYSiSgSidRoZChVsRodSYQdVGzYsGEaPny4Zz9/8eLFeuGFF3Ty5EkNGTJEixcvLvm5p512mrZs2aJwOKz58+dn729vb9f27du1cOFCT98bSjdkyBBdc801mjp1aq/7m5ub1dnZqccee0yzZ8/W/fffr7Fjx+rDDz/U/v37tWnTJoKOy4wJOo7jaNmyZdq8ebN27dqlSZMmeT0kDEKhqdvcTzm5twETtbe3Z0POyZMn1d7e3iu09CfzuC1btmRvt7e3a8uWLVq4cGHJrwPvLViwQMeOHdPRo0f7fO/gwYO64IIL9Nxzz2nIkC9Pw9dff73uv//+Wg4zEIwJOkuWLNGTTz6pZ599VmeeeaY++ugjSdLIkSN1+umnezw6lCqdThe8Pp25nU6nvRgW4Ir8UJK5LWlQYWfbtm06efIkIccyR48e1ZgxY3qFnIz6+noPRmQ3Y4LO+vXrJUlNTU297t+4caP+8i//svYDwqD0t+SUmZxTWLpbG1u3blV9fX3BANHe3q6enh4tWLCg5NcrNPNSaIamFPPnz8+GnCFDhhByLNPY2KjHHntMt912m7773e9q6tSpGjp0qNfDspYx0dE51dywzxchB7YJhUJqaWlRPB7vdX/msl8oFPJoZHapr6/Xli1b1N7e3uv+TGAp95N1T09PwZmX+fPna+HCherp6Sn5tQpd/oI92traNGPGDD300EP6xje+oTPOOEN//Md/rLa2Nh07dszr4VnHmBkdICgK1SyxLNV9/dXDNDQ0lD3T09/sTzkzMvkzQ2vXri06IzSYmSd4b/To0XrppZe0f/9+vfjii9q/f7927dqllStXasOGDdq3b5/GjBnj9TCtQdABfIimirVRqB6moaFBnZ2dfYqIcwNItRS6/HXJJZfof/7nf/qEnVqMB9U1ffp0TZ8+XdKp9gh33XWX1q1bp/vvv5+iZBcZc+kKCJrm5uZsn6FwOEzI+UI0Gu1zWS/j5Zdf1nPPPVfW682fPz97iWjIkCFasWKFFi5c2OuyVq1WPhW6/JW59CVJb7/9dk3HUytbt24tenmuvb1dW7durfGIam/o0KFqbW2VJL3++usej8YuzOgAPkVTxcIyNUxS7wL2tWvXas+ePbrmmmvKer3+loPXeuVTsUtQueNZsmSJdSuxMvVSkmo+i+aFDz/8UGPHju1z/1tvvSVJOu+882o9JKsRdAAfoqliccVqmNra2jRz5sySuxBLAy8H99PKJ7+Nx01B6x909dVXa9y4cVqwYIEuueQS9fT06MCBA1q7dq2GDx+u22+/3eshWoWgA/gMTRUHVqiG6e677y5rpdRAy8E7OzsH3fivkEqXs1fSiNAE5c6iud0eoJbuvfdePfvss1q3bp0+/PBDJZNJjR07Vt/85je1cuVKXXrppV4P0SoEHcBnaKpYmubm5mzICYfDuvPOO7Vu3bqSn9/fcvDOzk51dnZW1PgvXyWXZ9xoROgHpYST3Hqp/t6byZe7vvWtb+lb3/qW18MIDIIO4DO2N1V0qyFifg3T2rVry5rRKfZpv729vVfIkQbf+C/XYC/PuNmI0GsDhZOGhoaSZ62CdrkLg0fQAVBTxYqJcy/ZDaRYDdPMmTO1aNGiisbX30xP5vuDNZgi52qOp9YG6l1U7iwa22WgFAQdADVVbkPE/Bmg/EAUjUYVjUZ14sQJtbW16YMPPtDv/d7vDbp2w63Gf/29RjlFxZWOx2+1LP31LhrMrJXNRdpwB310AAP110smHo/7fi+s5uZmxWIxtbS0KBKJ9Nv1OX9LjEwNk6ReW2LceeedmjlzpiS5urWD22q9vYPbW124Ib930cUXXzzo7TPYLgMDYUYHMJAbl3+8ll9MXKz+KH8GKBPyCoWjGTNmaPTo0dq9e7cvaze8KCr2Yy1LfjgpNuMk9X9cbCnSRnURdAAD2bAfVjkNEcvdEsOPtRteFhX76Xi4FU5sKtJGdRF0AEP5cT+sUldUFSsm3rFjh+bMmVP0PYRCoZK3xPBb7YbXRcV+OB5uhhOvjyfMQY0OYDC/7YeVX0+TkQkyoVCoaEPEWCymXbt29fv8dDrdawaoP36r3ViwYEG/l2eqXRDsh+PRXzgZqBYnn9fHE+YIRf1eteiiZDKptrY2rVy5UpFIxOvhABWLx+N64YUXsif/UCikWbNmeTaeWbNmZcNOZiz5wWbnzp2aO3dun1CWea7jONq4cWOf50tSLBZTR0dHn58hSSdOnNDevXs1bNgwdXR0ZGcOli9fni3Ira+v1+TJk2t+XLyWO5Pi5fFoaGgo+vMmT56shoaGmo2lFlKplD777DM1NjZq2LBhXg/HaJWcv7l0BRjKr/thDXRJbaCGiLnvJfN8SSVvidHR0aHt27dTu/EFalkQdAQdwEB+3w+r1BVVpTy/vr5e0Wi05C0xqN3ojeOBoCPoGMitFvowl9/3wypnRVUpzy+m0GteffXVGj58eMHHB3HmotoNEAG/oxjZQKUUfMJuhWY4Mpqbmz0NurmzTclkMtsYcKDiYbee76atW7cWLdptb2/X1q1bazyiYOC4w03M6BjIhh4q8KdKZwsrvaTmt0tyJu+QbTKOO9xE0DGUH3uowHyVdlyu9JKa3y7J+bGrcBBw3OGmOsdxHK8HUSvd3d0aOXKkEomERowY4fVwXBGJRLK1DMlk0uvhwALFVnP5PUgfOXJE69at0+jRo4vW6AxW5iSb6UHDybY2TD/ux44d09GjR3XHHXdozJgxXg/HaJWcv6nRMVihgk+gUuVsuBkU+ZtQmnSyNRnHHW4g6BjKTwWbsI/fOi57zQ9dhYOI4w43UKNjIL8VbMI+lS4Ptwk7ZHuD4w63EHQM5LeCTdjFrx2XvUBXYW9w3OEmgo6BBmqhDwxWfsjJLDfPXBqVvvwbC0JzSroKe4PjDjcRdABk5c8WZpabx2IxxWKx7GxhqcvNTWdTV+GtW7eqvr6+4Ljb29vV09Pjmx2/bTru8B7FyIBhotFo0aLzeDyupqamfr/f3wxMfsfl3BVYuT+blVjmyTThyy/ozVwmqq/ndAA7MaMDGGagpn5z5sypqOlfPluaU5o0o1ENNOFDUBF0AMOUsgVIfgFxpbMwle5G7gdsK9A77Gzbts3IJnxAuQg68Ay7sA/eQLMsbs/CmLjcPH8GJ/ck39nZqYsvvjgbfoJ0sp8/f3425NCED0HARVl4hl3YKzNQUz+3mv75pTnlQLVJ999/f6/7CtWkzJ8/Xw0NDers7NS2bdsCF3IkmvB5ZaC/Xz7UVQ9BB57JLXTN/ANAoWvpBtoCxI0tQoo1p/Qi7JQbjOfPn6+FCxf2Cjvt7e3q7OxUfX29enp6AjejkXuZ7uGHH+5zfFA9fLDzTigaoBiZTCbV1tamlStXKhKJeD0cSJo1a1b2H4A1a9bohRdeIOSUIDeAdHR0ZI9hKBTSrFmzBvx+qXbu3Km5c+f2+X1kfm/pdFpNTU0uv7vCcv9WCr3P2267TXv37tWwYcMUDoclSZMnT87O7Gzfvl1vvfWWGhoadOTIkeyMRn19vSZPnlyT9+ClQoXHuccnKMehllKplD777DM1NjZq3rx5/f798m9e/yo5f7N7OXyBXdhLV+wfx9xVVzt27Cj6fdP/Uc28j8xMVeb99Ld7+ZIlS7KhJrcZXZBWHQV91ZkXCu1eXuzvF/2r5PxN0IHn+D9+eQYq4n7xxRcLzsJkvm9DkXehYFws6OT2ienp6VFDQ4NWrFjR5/tBCDuorUJBR+KD3WBUcv5m1RU8xb5K5RtoC5D+jptpx7RQqMvUHoVCoQFXgOWGmJ6eHv3iF79QZ2en2tvb+6zGYlsB1IKJKxhNR9CBZ9iFHQPJb46Yf3kutzniLbfc0uu5xWZqCu2CzUwOaoEPdt4g6MAz7MKOgTQ3N2vHjh1qaWnRrl27suEmU4MknZqJaWlp0YkTJ3ptY8DGkPATPth5hxodAL6W+6k3s9IrE3Jyt7T49NNPNXTo0ILFyIAXcmt0fvSjH9EgtQLU6ACwVu6n3nQ6nV06nv/pOFOMDPjRQLV1qB4aBgIwSiqVsmKZPIDaIOgAFrC5vXxubUNmNgcASkXQASxQTnt5k0JRbsiRlF2WK8mT/bYAmIegA1ignH3D8kNRJvgUCkVeB59ChceZjUUlaceOHZ6NDYAZKEYGLJFbtLt69eqiXabzl7Tm9qrJfXz+bIoXMiGsv2W5XjdcY2sFwN8IOoBFmpubsyEnHA4XDQC5QaFQ3Yuf9sXye7+lzKaYkgo2Jly4cKFXQwMggg5glXLay+eHonvvvXfA2SAv+H1Zbibc5IYd9s8C/IOgA1ii3Pby+aFIUva/+5sNQl+5YWfbtm06efIkIQfwCaOKkffs2aMFCxbovPPOU11dnZ555hmvhwT4QrE6lvwC5UKPzxT3trS09JkNQunmz5+vIUOG6OTJkxoyZAghB/AJo2Z0jh8/rq9//ev6/ve/rz/7sz/zejiAb5RTxzJQ/c29994rif13ytXe3p4NOSdPnuy1QzoA71QcdDJbZdXV1VU8mIHMmzdP8+bNq/rPAUxTTh1LfijKX12Vu+cOYac0+TU5hXZIB+CNQQedxx9/XOvWrdMvfvELSdLFF1+s5cuX6wc/+IFrg6tUMplUMpnM3u7u7vZwNIA/5Iciv69qKteJEydq+vM6Ojq0fft2XXPNNZo1a5aOHTumWbNmKZVKacuWLUqlUrryyitrOiYbPffcc6qvry94LDs6OtTT06Orr77ag5EVV+u/RRQ2qKDT3NysdevWadmyZWpsbJQk/ed//qfuuOMOvf/++1q9erWrgxysNWvWaNWqVV4PA/A1v69qKtXQoUOzuxt/9tlnNfu5x44d08yZMzV16lQdPXo0e//UqVN14sSJ7A7WqMznn3+uPXv26MSJE5oxY0b2/pdffll79uzRzJkzfXmcR44cqaFDh3o9jECrczLXnsowZswYPfTQQ/r2t7/d6/6nnnpKy5Yt05EjR1wbYDF1dXXavHmzrr322qKPKTSjM378+EFt8w74WTQaVSgUKrq6KvdylMkGep/Hjx/X3/zN33gwMtTC2rVr1dbWprvvvlt33nlnn9t+lAngqEx3d3f2g0y55+9Bzeik02lNnz69z/3Tpk3TyZMnB/OSVRGJRBSJRLweBlB1ud2Nc0OAH7obu6mU9zlmzBivhocqW7NmjYYNG6aWlhb9wz/8g6/6PcHHnEFYunSpc8cdd/S5/84773RuvfXWwbxk2SQ5mzdvLus5iUTCkeQkEokqjQrwTiwWcyQ5sVis4G1bBOV9orhwOOxIcsLhsNdDQY1Ucv4ueUZnxYoV2f+uq6vTY489pueff17f+MY3JEl79+7VwYMH9b3vfc+9FJbn2LFjeuedd7K3f/nLX+rAgQMaNWqUJkyYULWfC5ig1L2uTBeU94nCyun+DUhl1OjMnj27tBesq6vajsK7du0qOI4bb7xRTzzxxIDPr+QaH2CKSCSSPRHk1qjZJijvE18q1v2boGu/mtTo7Ny5s+yBua2pqUkl5jLATum09NJL0ocfSmPHSldcIYVC2W8H5dNuUN4nvjTQLva5t4FcRm0BAQTaj38sTZwozZ4tfec7p/534sRT96v4tg7V3MohGo0Wff14PF6VlV5evE+/8eK4e62/fk+xWMy4fk+oIbcLhvyMYmQY6+mnHaeuznGk3l91dY5TV+fErr++YEFutQt1i71+tX5urX+eX3EcEDQ1KUYG4JF0Wrr99lPRJp/jSHV1Sm/bplg0WvPuxoUuHVSzbsK2Ls6DVevjDphsUA0DTUUxMoy0a9epy1QD2blTamqq9mgKypxkMzUznGxrg+OOoKjk/E3QAfzuqadO1eQM5Mknpbxu5bXEKihvcNwRBJWcvylGBnL4sshz7Fh3H1cFhVZBofo47sDACDpAjswWA/knjMwlglDOUu6aueIKadw4qa6u8Pfr6qTx4089zgOsgvIGxx0okcuF0b7GqiuUwpdbDGRWXeWvvMrc9/TTngyL1T/e4LgjaFh1BbjIl1sM/J//I/2//3dq9dWvfvXl/ePGSf/3/576fgnc3uWcVVDe4LgDpaMYGSjCl0WeA3RGHkixJcjlLk12OzCZiGMA1A7FyIDLfFvkGQqdWkL+7W+f+t8ya4YyXWRzazkG03/Fl7VMNcYxAAzh+oU0H6NGB6XwZY2OyzLvKRwOD/q9BeE4DYRjANRGJedvgg6QI0hFnpmQEw6HB/0abgQm03EMvtTa2lr0/cdiMae1tbW2A4I1Kjl/c+kKyBGUjQPdujTX3NycfY1wOBzIrrymHYNq9orich58qQrBy7eY0QHcvdzCbIZ5x6Das5ZczkM1cOmqRAQdBJ2bJzlOaOYeg2qP27TwB/8j6JSIoIOgc6uGIki1TMWYfgyqHUbcqAEDMio5f9NHB0DZatlDxq/9asoZl1/fQ7V6RbGrOtxW0fnb9djlY8zoAOYxfebEcfz5Hqo1o2Pq5Tz4G5euSkTQAcxkw8nTT++hWmPxY6CDHQg6JSLoAOayocDVD++hmmGEPjqoFmp0SkSNDmA2X+4/Viav34Nf64WA/lRy/mb3cgBGKNTk0LQCVz+8h/5CjGnHEygFnZEBl1Wz82xQ5W48mkwm+2xMagIb3gO5vezjAAAUlklEQVRgJNcvpPkYNTqoBQoy3WXD8bThPQBequT8zaUrwGWZ6f+Wlpbs7dxP8yZfHvCivqO//ccy3/c7G94DYCqCDlAFuWFn9erV1jRNy2zaKPWu58gNcoPRX4AKhUJFg4Apx5O6GMA71OgAVWLartalyOzinltb4sZsFbteA6gWgg5QJYVW2JguMzORCTuRSKTXTM5gL1tVK0ABAEEHqAJbV9jkXrrKna2SVPHMS27YyQ1QhBwAFalCcbRvseoKtWD7CpvM+1BOh1833xe7XgPIV8n5mxkdwGX9rbCJxWKssPlCoX5Dmct9oVDImst9ADxWheDlW8zoAJXJnZXKnXkZzGxVsY0l58yZ0+t/TZoBY68noDroowMEkJc9bST1KrSWVPZsVe4S/F27dmnHjh2aM2eOduzYkZ0Ry9Q65T7ez6q1/B5ABaoQvHyLGR27BP3Ts1e1QMVmYgb78zLPD4VCRd+PSb9Lt48PgMrO3wQdGMv2ot9S1PqkWq1jblsBcuZ4ZN5XEP4WgWoi6JSIoGMfPj3X9qRajVk0W0OBbeEN8BJBp0QEHTvZeqIsh6knVVuDKn+TgLsIOiUi6NjL1BO9G0w9qdp66dHW8AZ4iVVXCLRCWy2YsELHDfnbJJi0SsnGHb0LbVtRaDd7ADVUheDlW8zo2Ce390ru7dxP0yat2CmHrTMiflJuTVLQVwIC1cKMDgIjt3dM5tNzpvfK3LlzdcUVV2T3S8r0ZrG1d4mNMyKSN/2Biim3L05/42ImB/BIFYKXbzGjY77c2YrcT8+53XRzbzOrYR6/zVRRcwN4j2LkEhF07DDQ1gG1KsrlMkX1+C1cmFrwDdiCoFMigo45BgoRTU1NBU88tVx95beZB9v4LVwEeWUf4DWCTokIOuYoJUTkn3i8ODH6bebBNn4JF34LXUDQEHRKRNAxS38hIv/Ek1+TU8vAMdiTIJe++ueXcEGYBbxH0CkRQcc8hU52+Sea/ELk/OfW4oQ0mJkHLn0V55dwwe8I8AeCTokIOmbKDRGFTjCtra1FV1nVYmakkpkHv5zQ/cTtcFHJzBmzboA/EHRKRNAxT36IaGpq8tWJx42g4pdLNH7hVrjIvE6h30luQbubx5tgBFRHoILOww8/7EycONGJRCLO1KlTnT179pT8XIKOWfw+2+HmzINfim5tUqimK/e/q/G3xKUuoDoCE3Q2bdrkDB061Hn00UedN99807n99tudM844w/nf//3fkp5P0DGHCScMtz69M6NTPcUCTjWPs98DOmCiwASdP/zDP3RuvvnmXvddcsklzt13313S8wk65gjKJYBanBSDciyLyQ+StZg5I7wC7gpE0Ekmk04oFHJ+/OMf97r/tttuc2bOnFnSaxB04Ce1mrUyYXas2vJDTi3eN5cjAfcEYlPPI0eOKJ1O69xzz+11/7nnnquPPvqo4HOSyaSSyWT2dnd3d1XHCJSjVptyZl4vd3PK3E0pbd9sMh6PK5VKZW/fe++9klRws063f2Y4HFYqlVI8Hrf+OAO+VYXgVRUffPCBI8n56U9/2uv+1atXOw0NDQWf09ra2uuafOaLGR0EURAvp+QXHhcrUK7Gz6RGB3APl66KXLr6/PPPnUQikf06ePAgQQeBFqTLKZmAkd+SID/suFmjxGVCoDoCcekqHA5r2rRp6ujo0OLFi7P3d3R0aNGiRQWfE4lEFIlEajVEwNeCdjmllEuD0Wi05j8TQI1VIXhVTWZ5+eOPP+68+eabzvLly50zzjjDef/990t6PsXICCoupwAwWSBmdCTpL/7iL3T06FHFYjF9+OGHmjJlirZt26bzzz/f66EBvlWo8LhQgTIA2MiooCNJt956q2699VavhwEYg8spAIKsznEcx+tB1Ep3d7dGjhypRCKhESNGeD0cAABQgkrO3/VVGhMAAIDnCDoAAMBaBB0LRaNRxePxgt+Lx+OuL6kFAMCvCDoWCoVCamlp6RN2MqtvQqGQRyMDAKC2jFt1hYEFfW8jmCsajSoUChX8G43H41Vp8gfAbgQdS+WGndWrVyuVShFy4HuZ2Uipd2+f3KAOAOVgebnlIpFItu1/7k7ugF/lzz4yGwmgkvM3MzoWC9reRrADs5EA3EQxsqVyPwUnk0nFYrGCBcqAHzU3N2cDejgcJuQAGDRmdCzE3kYwHbORANxC0LEQexvBZMVqdCQCOoDyUYwMwDeKFR5TkAwEG8XIAKzAbCQAtzGjAwAAfI3dywEAAAog6AAAAGsRdAAAgLUIOpaIRqNFmwHG43E2QgQABBJBxxKZzRDzw05mWW4oFPJoZAAAeIfl5ZYo1PmY3iMAgKBjebllMuEm0zqfkBNs0WhUoVCo4N9APB5XOp3msiYA32N5ObLYDBG5uKQJIOi4dGUZNkNELi5pAgg6go5F2AwRheSGndWrV3NJE0CgUKNjCTZDxEAikUh2ti+ZTHo9HAAoGZt6gs0Q0S8uaQIIKmZ0AMsVu6TJLB8AUzCjA6CgQqGmUIEyANiKoANYjEuaAIKOS1cAAMDXaBgIAABQAEEHAABYi6ADAACsRdABAADWIugAAABrEXQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0LFYNBpVPB4v+L14PK5oNFrbAQEAUGMEHYuFQiG1tLT0CTvxeFwtLS0KhUIejQwAgNoY4vUAUD3Nzc2SpJaWluztTMiJxWLZ7wMAYKs6x3EcrwdRK5Vs826yTLgJh8NKpVKEHACAUSo5fxN0AiISiSiVSikcDiuZTHo9HAAASlbJ+ZsanQCIx+PZkJNKpYoWKAMAYBuCjuVya3KSyaRisVjBAmUAAGxEMbLFChUeFypQBgDAVgQdi6XT6YKFx5nb6XTai2EBAFAzxhQj33fffWpvb9eBAwcUDof1ySeflP0aQS5GBgDAVIEoRk6lUrruuut0yy23eD0UAABgCGMuXa1atUqS9MQTT3g7kICKRqMKhUIFa3ri8bjS6TRbSgAAfMeYGZ3BSCaT6u7u7vWFwWE7CQCAiYyZ0RmMNWvWZGeCUBm2kwAAmMjTGZ1oNKq6urp+v/bv3z/o11+5cqUSiUT26+DBgy6OPniam5uzfXgikQghBwDge56uujpy5IiOHDnS72MmTpyo0047LXv7iSee0PLly1l15SG2kwAA1FIl529PL12NGTNGY8aM8XIIKFOh7SSY0QEA+JUxxchdXV06cOCAurq6lE6ndeDAAR04cEDHjh3zemiBwXYSAADTGFOM3NLSon/+53/O3v793/99SdLOnTvV1NTk0aiCg+0kAAAmMiboPPHEE/TQ8RDbSQAATGTMFhBuoBgZAADzBGILCAAAgHIRdAAAgLUIOgAAwFoEHQAAYC2Czhei0WjRfjDxeJyduQEAMBBB5wvszg0AgH2M6aNTbezODQCAfeijkycTbjJ7ORFyAADwViV9dAg6BbA7NwAA/kHDQBcV2p0bAACYiaCTg9250R9W5gGAeShG/gK7c2MgmZV5Uu+/hdy/HQCAvxB0vsDu3BgIK/MAwDwUIwNlYmUeANQWq65KRNCBW1iZBwC1w6oroIZYmQcA5iDoAGVgZR4AmIViZKBErMwDAPMQdIASsTIPAMxDMTIAAPA1ipHhO3QRBgD4AUEHVZHpIpwfdjJ1LqFQyKORAQCChBodVAVdhAEAfkCNDqqKLsIAgErRGblEBB1v0EUYAFAJipHhW3QRBgB4iaCDqqGLMADAaxQjoyroIgwA8AOCDqqCLsIAAD+gGBkAAPgaxcgAAAAFEHQAAIC1CDoAAMBaBB0AAGAtgg4AALAWQQcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLSOCzvvvv6+bbrpJkyZN0umnn64LL7xQra2tSqVSXg8NAAD42BCvB1CKt99+Wz09PdqwYYMuuugivf766/rrv/5rHT9+XA8++KDXwwMAAD5V5ziO4/UgBuOBBx7Q+vXr9d5775X8nO7ubo0cOVKJREIjRoyo4ugAAIBbKjl/GzGjU0gikdCoUaP6fUwymVQymez1HOnUAQMAAGbInLcHNTfjGOidd95xRowY4Tz66KP9Pq61tdWRxBdffPHFF198WfD17rvvlp0ZPL10FY1GtWrVqn4fs2/fPk2fPj17+9ChQ5o1a5ZmzZqlxx57rN/n5s/ofPLJJzr//PPV1dWlkSNHVjb4AOvu7tb48eN18OBBLgFWiGPpHo6lOziO7uFYuieRSGjChAn67W9/q7POOqus53p66Wrp0qW6/vrr+33MxIkTs/996NAhzZ49W42NjXrkkUcGfP1IJKJIJNLn/pEjR/JH54IRI0ZwHF3CsXQPx9IdHEf3cCzdU19f/mJxT4POmDFjNGbMmJIe+8EHH2j27NmaNm2aNm7cOKg3CwAAgsWIYuRDhw6pqalJEyZM0IMPPqjDhw9nv/fVr37Vw5EBAAA/MyLoPP/883rnnXf0zjvvaNy4cb2+V06JUSQSUWtra8HLWSgdx9E9HEv3cCzdwXF0D8fSPZUcS2P76AAAAAyEQhcAAGAtgg4AALAWQQcAAFiLoAMAAKwV2KCzcOFCTZgwQaeddprGjh2rG264QYcOHfJ6WMZ5//33ddNNN2nSpEk6/fTTdeGFF6q1tVWpVMrroRnnvvvu0+WXX65hw4aV3fkz6P7xH/9RkyZN0mmnnaZp06bppZde8npIRtqzZ48WLFig8847T3V1dXrmmWe8HpKR1qxZoz/4gz/QmWeeqXPOOUfXXnutOjs7vR6WcdavX6/LLrss23CxsbFRP/nJT8p+ncAGndmzZ+vf/u3f1NnZqaefflrvvvuu/vzP/9zrYRnn7bffVk9PjzZs2KA33nhD69at0z/90z/pnnvu8XpoxkmlUrruuut0yy23eD0Uo/zrv/6rli9frr/7u7/Tz3/+c11xxRWaN2+eurq6vB6acY4fP66vf/3r+tGPfuT1UIy2e/duLVmyRHv37lVHR4dOnjypq666SsePH/d6aEYZN26c2tratH//fu3fv19z5szRokWL9MYbb5T1Oiwv/8KWLVt07bXXKplMaujQoV4Px2gPPPCA1q9fr/fee8/roRjpiSee0PLly/XJJ594PRQj/NEf/ZGmTp2q9evXZ++79NJLde2112rNmjUejsxsdXV12rx5s6699lqvh2K8w4cP65xzztHu3bs1c+ZMr4djtFGjRumBBx7QTTfdVPJzAjujk+vjjz/Wv/zLv+jyyy8n5LggkUho1KhRXg8DAZBKpfTKK6/oqquu6nX/VVddpZ/+9KcejQroLZFISBL/LlYgnU5r06ZNOn78uBobG8t6bqCDzl133aUzzjhDo0ePVldXl5599lmvh2S8d999Vw899JBuvvlmr4eCADhy5IjS6bTOPffcXvefe+65+uijjzwaFfAlx3G0YsUKzZgxQ1OmTPF6OMZ57bXXNHz4cEUiEd18883avHmzvva1r5X1GlYFnWg0qrq6un6/9u/fn3383/7t3+rnP/+5nn/+eYVCIX3ve98ra0sJm5V7LKVTe5Jdc801uu666/SDH/zAo5H7y2COI8pXV1fX67bjOH3uA7ywdOlSvfrqq3rqqae8HoqRGhoadODAAe3du1e33HKLbrzxRr355ptlvYYRe12VaunSpbr++uv7fczEiROz/53ZPX3y5Mm69NJLNX78eO3du7fsaTEblXssDx06pNmzZ6uxsVGPPPJIlUdnjnKPI8ozZswYhUKhPrM3v/nNb/rM8gC1tmzZMm3ZskV79uzps08jShMOh3XRRRdJkqZPn659+/bphz/8oTZs2FDya1gVdDLBZTAyMznJZNLNIRmrnGP5wQcfaPbs2Zo2bZo2btyo+nqrJgorUsnfJAYWDoc1bdo0dXR0aPHixdn7Ozo6tGjRIg9HhiBzHEfLli3T5s2btWvXLk2aNMnrIVnDcZyyz9NWBZ1S/exnP9PPfvYzzZgxQ2effbbee+89tbS06MILL2Q2p0yHDh1SU1OTJkyYoAcffFCHDx/Ofu+rX/2qhyMzT1dXlz7++GN1dXUpnU7rwIEDkqSLLrpIw4cP93h0/rVixQrdcMMNmj59enZGsaurizqxQTh27Jjeeeed7O1f/vKXOnDggEaNGqUJEyZ4ODKzLFmyRE8++aSeffZZnXnmmdkZx5EjR+r000/3eHTmuOeeezRv3jyNHz9en376qTZt2qRdu3Zp+/bt5b2QE0CvvvqqM3v2bGfUqFFOJBJxJk6c6Nx8883Or371K6+HZpyNGzc6kgp+oTw33nhjweO4c+dOr4fmew8//LBz/vnnO+Fw2Jk6daqze/dur4dkpJ07dxb8G7zxxhu9HppRiv2buHHjRq+HZpS/+qu/yv7/+itf+Yozd+5c5/nnny/7deijAwAArEUxBQAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAaxF0AACAtQg6AADAWgQdAABgLYIOAACwFkEHAABYi6ADwGjbt2/XjBkzdNZZZ2n06NH60z/9U7377rteDwuATxB0ABjt+PHjWrFihfbt26cXX3xR9fX1Wrx4sXp6erweGgAfYFNPAFY5fPiwzjnnHL322muaMmWK18MB4DFmdAAY7d1339V3vvMdXXDBBRoxYoQmTZokSerq6vJ4ZAD8YIjXAwCASixYsEDjx4/Xo48+qvPOO089PT2aMmWKUqmU10MD4AMEHQDGOnr0qN566y1t2LBBV1xxhSTp5Zdf9nhUAPyEoAPAWGeffbZGjx6tRx55RGPHjlVXV5fuvvtur4cFwEeo0QFgrPr6em3atEmvvPKKpkyZojvuuEMPPPCA18MC4COsugIAANZiRgcAAFiLoAMAAKxF0AEAANYi6AAAAGsRdAAAgLUIOgAAwFoEHQAAYC2CDgAAsBZBBwAAWIugAwAArEXQAQAA1iLoAAAAa/1/kcp2fHeFWvMAAAAASUVORK5CYII=",
"text/plain": [
"PyPlot.Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"using Distributions, PyPlot\n",
"N = 100\n",
"generative_dist = MvNormal([0,1.], [0.8 0.5; 0.5 1.0])\n",
"function plotObservations(obs::Matrix)\n",
" plot(obs[1,:]', obs[2,:]', \"kx\", zorder=3)\n",
" fill_between([0., 2.], 1., 2., color=\"k\", alpha=0.4, zorder=2) # Shaded area\n",
" text(2.05, 1.8, \"S\", fontsize=12)\n",
" xlim([-3,3]); ylim([-2,4]); xlabel(\"a\"); ylabel(\"b\")\n",
"end\n",
"D = rand(generative_dist, N) # Generate observations from generative_dist\n",
"plotObservations(D)\n",
"x_dot = rand(generative_dist) # Generate x∙\n",
"plot(x_dot[1], x_dot[2], \"ro\");"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Log-Likelihood for a Multivariate Gaussian (MVG)\n",
"\n",
"- Assume we are given a set of IID data points $D=\\{x_1,\\ldots,x_N\\}$, where $x_n \\in \\mathbb{R}^D$. We want to build a model for these data."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- **Model specification**: Let's assume a MVG model $x_n=\\mu+\\epsilon_n$ with $\\epsilon_n \\sim \\mathcal{N}(0,\\Sigma)$, or equivalently,\n",
"\n",
"$$\\begin{align*}\n",
"p(x_n|\\mu,\\Sigma) &= \\mathcal{N}(x_n|\\mu,\\Sigma) \n",
"= |2 \\pi \\Sigma|^{-1/2} \\mathrm{exp} \\left\\{-\\frac{1}{2}(x_n-\\mu)^T\n",
"\\Sigma^{-1} (x_n-\\mu) \\right\\}\n",
"\\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Since the data are IID, $p(D|\\theta)$ factorizes as\n",
"\n",
"$$ \n",
"p(D|\\theta) = p(x_1,\\ldots,x_N|\\theta) \\stackrel{\\text{IID}}{=} \\prod_n p(x_n|\\theta)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true,
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- This choice of model yields the following log-likelihood (use (B-C.9) and (B-C.4)),\n",
"\n",
"$$\\begin{align*}\n",
" \\log &p(D|\\theta) = \\log \\prod_n p(x_n|\\theta) = \\sum_n \\log \\mathcal{N}(x_n|\\mu,\\Sigma) \\tag{1}\\\\\n",
" &= N \\cdot \\log | 2\\pi\\Sigma |^{-1/2} - \\frac{1}{2} \\sum\\nolimits_{n} (x_n-\\mu)^T \\Sigma^{-1} (x_n-\\mu) \n",
" \\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Maximum Likelihood estimation of mean of MVG\n",
"\n",
"- We want to maximize $\\log p(D|\\theta)$ wrt the parameters $\\theta=\\{\\mu,\\Sigma\\}$. Let's take derivatives; first to mean $\\mu$, (making use of (B-C.25) and (B-C.27)),\n",
"\n",
"$$\\begin{align*}\n",
"\\nabla_\\mu \\log p(D|\\theta) &= -\\frac{1}{2}\\sum_n \\nabla_\\mu \\left[ (x_n-\\mu)^T \\Sigma^{-1} (x_n-\\mu) \\right] \\\\\n",
"&= -\\frac{1}{2}\\sum_n \\nabla_\\mu \\mathrm{Tr} \\left[ -2\\mu^T\\Sigma^{-1}x_n + \\mu^T\\Sigma^{-1}\\mu \\right] \\\\\n",
"&= -\\frac{1}{2}\\sum_n \\left( -2\\Sigma^{-1}x_n + 2\\Sigma^{-1}\\mu \\right) \\\\\n",
"&= \\Sigma^{-1}\\,\\sum_n \\left( x_n-\\mu \\right)\n",
"\\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Setting the derivative to zero yields the **sample mean**\n",
"\n",
"$$\\begin{equation*}\n",
"\\boxed{\n",
"\\hat \\mu = \\frac{1}{N} \\sum_n x_n\n",
"}\n",
"\\end{equation*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Maximum Likelihood estimation of variance of MVG\n",
"\n",
"- Now we take the gradient of the log-likelihood wrt the **precision matrix** $\\Sigma^{-1}$ (making use of B-C.28 and B-C.24)\n",
"\n",
"$$\\begin{align*}\n",
"\\nabla_{\\Sigma^{-1}} &\\log p(D|\\theta) \\\\\n",
"&= \\nabla_{\\Sigma^{-1}} \\left[ \\frac{N}{2} \\log |2\\pi\\Sigma|^{-1} - \\frac{1}{2} \\sum_{n=1}^N (x_n-\\mu)^T \\Sigma^{-1} (x_n-\\mu)\\right] \\\\\n",
"&= \\nabla_{\\Sigma^{-1}} \\left[ \\frac{N}{2} \\log |\\Sigma^{-1}| - \\frac{1}{2} \\sum_{n=1}^N \\mathrm{Tr} \\left[ (x_n-\\mu) (x_n-\\mu)^T \\Sigma^{-1}\\right] \\right]\\\\\n",
"&= \\frac{N}{2}\\Sigma -\\frac{1}{2}\\sum_n (x_n-\\mu)(x_n-\\mu)^T\n",
"\\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"Get optimum by setting the gradient to zero,\n",
"$$\\begin{equation*}\n",
"\\boxed{\n",
"\\hat \\Sigma = \\frac{1}{N} \\sum_n (x_n-\\hat\\mu)(x_n - \\hat\\mu)^T}\n",
"\\end{equation*}$$\n",
"which is also known as the **sample variance**."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Sufficient Statistics\n",
"\n",
"- Note that the ML estimates can also be written as\n",
"$$\\begin{equation*}\n",
"\\hat \\Sigma = \\sum_n x_n x_n^T - \\left( \\sum_n x_n\\right)\\left( \\sum_n x_n\\right)^T, \\quad \\hat \\mu = \\frac{1}{N} \\sum_n x_n\n",
"\\end{equation*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- I.o.w., the two statistics (a 'statistic' is a function of the data) $\\sum_n x_n$ and $\\sum_n x_n x_n^T$ are sufficient to estimate the parameters $\\mu$ and $\\Sigma$ from $N$ observations. In the literature, $\\sum_n x_n$ and $\\sum_n x_n x_n^T$ are called **sufficient statistics** for the Gaussian PDF."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The actual parametrization of a PDF is always a re-parameteriation of the sufficient statistics. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Sufficient statistics are useful because they summarize all there is to learn about the data set in a minimal set of variables. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Solution to Example Problem\n",
"\n",
"\n",
"We apply maximum likelihood estimation to fit a 2-dimensional Gaussian model ($m$) to data set $D$. Next, we evaluate $p(x_\\bullet \\in S | m)$ by (numerical) integration of the Gaussian pdf over $S$: $p(x_\\bullet \\in S | m) = \\int_S p(x|m) \\mathrm{d}x$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG2CAYAAAB20iz+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd4W+X1xz+SbHnEI8vZOyGDEVYYARogQAk1SYCW0QEtLWVTVimEIseRAqS0KVBKISWUX0sLFEqBpIZAGCEECMSEhOyQvYfjxNuSJb2/P44Uy47s2Bq2ZZ/P89xHvtfSva+kK92vzvs951iMMQZFURRFUZR2iLW1B6AoiqIoihIvVOgoiqIoitJuUaGjKIqiKEq7RYWOoiiKoijtFhU6iqIoiqK0W1ToKIqiKIrSblGhoyiKoihKu0WFjqIoiqIo7RYVOoqiKIqitFtU6CiKoiiK0m5JWKHz6KOPYrFYuOuuu1p7KIqiKIqitFESUugsWbKEv/71r4wePbq1h6IoiqIoShsm4YROeXk5P/7xj3nuuefo0qVLaw9HURRFUZQ2TFJrD6C53HbbbeTm5nLhhRcyffr0Ru/rdrtxu92H1/1+P8XFxXTr1g2LxRLvoSqKoiiKEgOMMZSVldGnTx+s1ubFaBJK6LzyyissXbqUJUuWNOn+jz76KNOmTYvzqBRFURRFaQm2b99Ov379mvWYhBE627dv58477+S9994jNTW1SY+ZMmUK99xzz+H1kpISBgwYwPbt28nKyorXUBVFURRFiSGlpaX079+fzMzMZj/WYowxcRhTzHnzzTe5/PLLsdlsh7f5fD4sFgtWqxW3213nf+EoLS0lOzubkpISFTqKoiiKkiBEc/1OmIjOBRdcwIoVK+psu/766xk5ciT333//UUWOoiiKoigdj4QROpmZmRx//PF1tnXq1Ilu3bodsV1RFEVRFAUSML1cURRFURSlqSRMRCccCxYsaO0hKIqiKIrShtGIjqIoiqIo7RYVOoqiKIqitFtU6CiKoiiK0m5RoaMoiqIoSrtFhY6iKIqiKO0WFTqKoiiKorRbVOgoitImyc/Px+Vyhf2fy+UiPz+/ZQekKEpCokJHUZQ2ic1mIy8v7wix43K5yMvL07YviqI0iYQuGKgoSvvF4XAAkJeXd3g9KHKcTufh/yuKojRGwnQvjwXavVxREo+guLHb7Xg8HhU5itIBieb6rUJHUZQ2T0pKCh6PB7vdjtvtbu3hKIrSwkRz/VaPjqIobRqXy3VY5Hg8ngYNyoqiKOFQoaMoSpsl1JPjdrtxOp1hDcqKoigNoWZkRVHaJOGMx+EMyoqiKI2hQkdRlDaJz+cLazwOrvt8vtYYlqIoCYaakRVFURRFadOoGVlRFEVRFCUMKnQURVEURWm3qNBRFEVRFKXdokJHURRFUZR2iwodRVEURVHaLSp0FEVRFEVpt6jQURRFURSl3aJCR1EURVGUdosKHUVRFEVR2i0qdBRFURRFabeo0FEURVEUpd2iQkdRFEVRlHaLCh1FURRFUdotKnQURVEURWm3qNBRFEVRFKXdokJHURRFUZR2iwodRVEURVHaLSp0FEVRFEVpt6jQURRFURSl3aJCR1EURVGUdosKHUVRFEVR2i0qdBRFURRFabeo0FEURVEUpd2iQkdRFEVRlHaLCh1FURRFUdotKnQURVEURWm3qNBRFEVRFKXdokJHURRFUZR2iwodRVEURVHaLSp0FEVRFEVpt6jQURRFURSl3aJCR1EURUkI8vPzcblcYf/ncrnIz89v2QEpCYEKHUVRFCUhsNls5OXlHSF2XC4XeXl52Gy2VhqZ0pZJau0BKIqiKEpTcDgcAOTl5R1eD4ocp9N5+P+KEorFGGNaexAtRWlpKdnZ2ZSUlJCVldXaw1EURVEiIChu7HY7Ho9HRU4HIJrrd8JMXT3zzDOMHj2arKwssrKyGDt2LO+8805rD0tRFEVpYRwOx2GRY7fbVeQojZIwQqdfv37MmDGDwsJCCgsLGT9+PJMnT2bVqlWtPTRFURSlBXG5XIdFjsfjadCgrCiQQEJn4sSJfO9732P48OEMHz6chx9+mIyMDBYvXtzaQ1MURVFaiFBPjtvtxul0hjUoK0qQhDQj+3w+XnvtNSoqKhg7dmyD93O73bjd7sPrpaWlLTE8RVEUJQ6EMx6HMygrSigJJXRWrFjB2LFjqa6uJiMjgzfeeINjjz22wfs/+uijTJs2rQVHqCiKkljk5+djs9nCCgSXy4XP52sz9Wl8Pl9Y43Fw3efztcawlDZOQmVdeTwetm3bxqFDh3j99deZPXs2H3/8cYNiJ1xEp3///pp1pSiKEqCh9GxN21baEtFkXSVURMdutzNs2DAAxowZw5IlS3jyySeZNWtW2PunpKSQkpLSkkNUFCUBMMBXQBEwoZXH0tpobZrmUQn4gYzWHojSZBJK6NTHGFMnYqMoitIYO4B/Av8A1gDDgPWApTUH1QYIFTvTp0/X2jQNcBCYCKQB/wP0Z3RikDBTVw8++CCXXHIJ/fv3p6ysjFdeeYUZM2Ywb948LrrooibtQwsGKkrHowaYC8wC5iPRHIBU4HLgWUC/DYSUlJTDadv6I/JIlgLjgArgauBfgDadaBk6RMHAvXv3cu211zJixAguuOACvvjii2aJHEVRWoa20nhxK/AQMAD4PvAeInLOBZ4H9gIvoSIniNamOTqnAP8FkoF/A3dSK5yVtkvCCJ3nn3+eLVu24Ha72bdvH++//76KHEVpg7Rm40Uv8BbwPWAw8DCwB+gJTAE2AguAn6MCJxStTdN0vgv8HZnufBr4XesOR2kKpgNRUlJiAFNSUtLaQ1GUdo3T6TSAcTqdYddjzTZjTJ4xpq8xhpDlAmPMa8YYd1yO2j5o6L2pv33q1KkNvn9Op9NMnTo13kNtU/zJyDlmNcYsauWxdASiuX4ntBlZUZS2SUuYWw3wEfAkYgz1B7Z3B64HfgkcE7OjtV+aWpsmGKkL/R/UjQZ1JO4AlgAvAj8ClgOdW3VESkMkjBk5FqgZWVFalniYW6uAl4EngBUh288DbkIMxpoNEx/qp5139DT0MuAkYBNwJeLb6egZfPGiw9TRURQlcQhnbo3mYrgd8UQ8BxQHtqUDP0N+XY+McrztCa+BEi9U+qDSD1V++bvKD34j0S+D/G0AuwXsVkixQmpgybZB5yRIC7FUNSdSl0gVlyMlExHdZwOvAVcA17TqiJRwqNBRFCXmNPTLH5rfi+hL4HHkQhIs8D8AuB24AegSs1G3ffwG9tXAHo8sewO3+zxw0CvLIS+UxbATQopFBE+3ZOhth37XO0gKiBy73c5DD4V/PzvKVNfpSHZfPvAbYDJSZ0dpO6jQURQlpsSi8aIPeBOYCXwesv084C7gUtp3/ZJKH2yqhg1VsK0atrnldqcb3E00G1iAdKtEZNKtkBaI1NgsslgAa2CexeMHjwG3X5Zqvwgmj5Hj7a2RZXUl7J7twuvxYEmWSN3Q21xMuMvBqHQYmQ7HdoI+9o5Vcfk3SMmC7cAfgd+27nCUeqjQURQlpkTTeLEaqVz8e6RiMUjNkh8hAuek2A+3VTEBEbGmAtZVwbeVIm52ehp+jA3oZYeegaWXHXokS8SlcxJ0SYIuyZBpE0ETzdgq/bVRoqIaePpRF189m8eZdzrpc4ODL592sfmZPN4y8OUNte93Lzucngln3ubgQX/7r7ichqSZ/wh4FDHD92nVESmhqBlZUZRW5wDwF+DPwL7Ats7AbcgUVa9WGlesKa6BbypgZYWIm7WVUNKA7uuWBMPSYFAq9E+FgSly28sOSa3geA0XjfH44ddTXTw1PY8L7hHxs75KPEJBrMBXY1Pw17TvissG8ep8jpyzT7XucNodakZWFCUh2YNMTz2DlNUH6I9UnL0RMXsmKsbAlmr4qhyWl8M35eEjNTZgaBqMSJdlWJosndvYt3O4SJ3dCn9yOcixy//zR0GVD74uhy/L4LMS+PTPLvw1tVNdZ9/tYnqeg7Oz5fHtBQvgBC4C/gZMA7q26oiUIBrRURSlxdmGTE89BwR/358E3Iek6Sa30riiwRjY7oYlZVBYBkvL4ID3yPsNSYXRGXBcwNMyNK19XfBDCUaBfvAbJ/afOljwlItdz+bR52YnI29ycGUPuCoHuibiGx4GA5yM1NR5FHigdYfTrtCIjqIoCcFGYAZSQr8msO1MwAFcQuLVIDnkFVHzZSl8UXpkxCbFIqLmpAwY3QmO7wSZHeRbN9xU14Y/OrizE7w9UwzKs29w8OIeuKw7XN8buie44LEA9wA/Raau7gHsrToiBVToKIrSAqwGHkFqjgQrGJ+PpOWeT+IIHL8RX80nJfBpCayprNvUMckigmZMpizHdZLaNB2RcFNdw9Kg4A8OpmXBxgofpEsm17/3w5tFcFUPuKE3dErglLprkMjkLuB9pO+a0rro1JWiKHFjOeBCOj4Hv2i+h6TfntVag2om1X5YXAoLD4m4qT8dNSQVzsiSLKNTMyE9gS/SLY0xMtX3zC5YETBp9UyG+wbAeQncT+EOxFh/HRK9VKInmuu3Ch1FUWLOGqSA2qsh265ABM4p9e7bFivoVvhgUQl8eBA+LRWxEyTdKsLmnGw4KwtydG4iaoyR1/sP22un/87vDPf1hx4xfH1b6lz7FDgHyAL2AqlR71GJ5vrdQYOqiqLEg7XAT4DjqRU5VwMrgdc5UuRAbQVdl8tVZ3vQ42GztUyIpNQLc4vgzm/hwuXw283wwSEROb3t8MMe8PQx8P6J8PuhMLm7ipxYYbHAdzrDv4+Dn/WSTLSPDsFVq+H9g7E7Tkuda2OBfkAp8F5M9qhERcx6qCcA0bR5VxSlYTYbY64zxliNMQSWycaY5U18vNPpNIBxOp1h1+NFaY0xc4uMufNbY874yphTC2uXy1cY8+cdxqwuN8bvj+swlHp8W2nMdatr34tnd8buPWipc+1WI5+DX8V0rx2XaK7fOnWlKErE7AMeRurgBLOoJgFTCR+9aYzgr+pgE9B4VdD1Gvi8BAqKxXfjCfkGHJoKF3WF8Z1hcKpEGpTWwWvg6Z3w4l5Zn9AV8gbGJhW/Jc61V5Fo5onAspjuuWOiHp0mokJHUWJDKdLTZyZQHth2AVI75LQo9puSknK4WWQsK+iaQLbU28Uwr1jaGgQZnArf7QIXdoHB2o2xzfFmETy6VfqfnZQBfxgam2KK8TrXguwBeiMZhQfoWM1n44F6dBRFaRHcwJPAUKTyazlwKjAfSaWNRuS4XK7DFx6Px3OEjyISimrgn3vhmtVw7Vp4eZ+InK5J4rn55yh49Vj4ZR8VOW2Vy7rDU8dAhg2WlcNt66EsTCHG5hCPc60+vYDhSLbh4pjvXWkOKnQURTkqPuBFYCTSXLMI+RJ/FVgCXBjl/kOLy7ndbpxOZ1jTaFPwGvjoINy9AXK/gSd2wMZqKd53URd4Yhi8PRru7S+ViXV6qu1zehY8P0IE6roquHND3Uy45hDLc+1oBJvQror5npXmoAUDFUVplHeRAmgrAuu9kdTxnxObL5BwFXSDt3l5eXXWG2OvR6Y53iyC/TW120d3gku7icjpKFWJ2yND0yTr7ab10hh16mZ4dAhYmyFUY3WuNZVRgds1MdujEgn6sVcUJSwrEYEzL7CejfTu+RWQHsPjhKugC7UXHJ+vgfbeSKXixaXw+n6pwxK8Z5ckmNgNJnWX7t9KLW2xblFTOSYdZg6FW76V1P9Zu+CWvk1/fDTnWiSMDNyq0Gld1IysKEoddgN5SAdmP9Jg83akXUNb6cZc4oU5RfCf/XX7S52SAd/PkWJzidYo0+eH0koorQoslVBWDVUeqPZAdY0s7hoReEGC3+D2JEhNhlQ7pCTJbXY6dM2Azp2gSyfolALTpx8Z1YDw0Y62SsEBmLpFjL7PjRCTclvkayT7sCdiTlYiR5t6KooSNVXAH5Cmm5WBbd8PrA9rrUHVY20l/HsfvFcM7sAFPsMm0ZvLu8OQNm4o9nhh7yHYfQj2l8pyoExuD1XWipZ4YU+C3sc7uOLnMlVTVAaPuhzMfCxxRA5Abjdppjr3AEzbAi8fC6ltUNj2CtzuR340tMEhdgg0oqMoHRyDVC3+NbA1sG0sInraQj8qv4HPSiV7qrCsdvvwNGkCeXEXSGtj/aWMgYMVsHU/bC2CXcWw+yDsL2tczFgtkJkGWWkSjclMg3Q7pCTXjdZYrSGNUC2AAY9Poj3VHnB7JRJ0qELGcbBcokOhLC1wUTgnD2uSHb/Xw9W/FJEzrJccr61T5pVsur01kkF3b//WHtGR1FDbvXwfkNOKY0l0tI5OE1Ghoyh1WYZkUX0cWO8PPIYUOmvtZKQqnxT1e3kvbA2UObEBF3SBq3uIybitZExVVMOmfbBxjwib7UVHCosgnVKgdxfomQ3dMqF7JnTPktvMtOaZa5tDjRcOlMPOYthxQJY7J6bg83qwJtm54Wl5kZNsMKIPnDEMThkMyW047v95CdyxQf5+YQSc0AansLoBxYiZ//hWHksio1NXiqI0iwOI52YWEtFJBX4D3E9sjcaRcLAGXt0Pr+6DkoA3tJMVrsgRgdOrDfSXKi6HtTthwx7YtFemoupjtUDfrjCgO/TvJuKmTxcRM60h0JKToFdnWU4dEjAee2tryRR94WL4BQ6Ky2HVdlle+UwEz3dGyXNpa4zNloy6/x2Ap3fBs8Nbe0RHko0InfKj3VGJGyp0FKUD4QeeB6YgYgckevM7YGBrDSrADjf8a6+YjIP+m34pcE0P8eB0asXpqbIqWLtLxM26XeKpqU/PbBjWCwbliLjp01U8MW2R+sbj4Pq0gTDtVw6WbITP1omg+2iVLINy4JKT4cSBbSeSBnBzH3i3WKY1vyqDUzNbe0R1Cc4C1jR6LyWetNGPoaIosWYJcFvgFiSM/jQwLk7HM6ZpF8SNVfDCHjEYB2vAHZsO1/WS7ClbK1xU/Qa27YcV22HlNvHahM7xWy0wqAcc0wuG9oKhPSEjQdLYj1ZLxmKR9dyTYc1OWLQWlm2BLfvhmfdExE0+DY5vI56YXnaY1A1eL4K/7VahoxyJCh1FaecUAw8Cf0Uu1pmAExE98fCcbq2G/ilN85qsrYSfhBQZGZsFP+0Fp2a0fNTAXQOrdsDyLbByO5TX89j06wYj+8DIvhK5SWsDU2iR0NRaMlYrHNdfltJK+HClLNuK4Kl3YPQAuPos8Re1Nj/tJYUivyiDbyul3k5bIfgZi7JrhRIFakZWlHaKAf6OFP0rCmz7CWI27h2H4y0phd9tl2hIJ5tUIv5RT0hqRLAYAz9dCz3t8Ive0pKhJSmrgm+2SsRizU6oCakXl5osouaEARK96NypZcfWFimrgnnLRPD4DSTb4MqxMG5U609n3bcRPjoE1/aEO/u17lhCORYpGLgAOLd1h5LQqBlZUZQ6rABuBRYF1o8D/kL8pqlWVcAj26S+yUkZ0nzxL7ugxsCVOZDVwDeNxQKzR7Rscb9KN3y9BQo3iucmtPhe90w4ebBEK4b2AlsCFz6JRwXkzDQRNmePhJcXwfrd8NIiWL8LfjKudaNcl3QVofNuMdzRN37Za80lWJOqDQWZOhwqdBSlHVGBTEvNRNohpCN9qe4iPtNUPiMemtUVko5+TQ8p4Dcm4JNYcAiyk+AHjRQQaQmR4/FK1GbJBpme8oU0hBzQHU4aJEufLq0XmYi1MLHZbGH7N4V6dCKlTxe4+1J4/xt440so3CRp9bd+V0zYTSWWz/nsbDn39tXA0vLac7C1UaHT+qjQUZR2wjtIFGdLYP1y4EmkNk6s+boMTs6sNQpvqBJTaIZNuocnWeC6nrCuEj4tgdMyYWALm3X9fli3G774FpZuFg9OkL5dYcxQOG0o5LSRWexYC5NwzSpj2ebBaoHvniiRr9kfSCba7+fCry6BwT2ato9YPucUK4zvDHMOwMeH2obQMUAwQU9nPlsR04EoKSkxgCkpKWntoShKzNhljLnaGENg6W+MeStOx1pSaszkFcb8cJUxxR5jfH7Z/kGxMacXGlPmlfWawPaPDhpz9Spj3i6K04DCsK/EmP9+Ycxv/mnMjbNqlwdfMubNL43ZeaDlxtJcnE6nAYzT6Qy7Hs0+7XZ71PtqiLIqYx59Q17nO5435tvdzR9fLJ7zOweMObXQmJ+sbvZD48JBU/u5rGrlsSQ60Vy/1YysKAmKAWYjZuMSpI/OXcA0INYFYg95IW+z1Cm5qseRHojdbrj1WzgrC+4bUBvVAfjRavl1fU9/8cPEwztR44Nlm+GTtVLnJki6XSI3ZxwjKeCtbZhtCsFoRrCQXyyiLykpKXg8UhzQ7XbHaKR1qa6BZ98TU3e6HX49qelFBmP1nPd6IHeFfBY+Oql1ay8BrAVGAZ2Bg607lIRHW0A0ERU6SnthA3Aj8FFgfQySPn5yHI613wP3b4JN1fCf46B7wOxT7pOpKoBqP7y+H/68E14aBYPTakXNQ5uh1At/Oib2Y9tbAgtXw+froSJw/bYAx/aHs0fA6IGSGZRoxFKYxEM4NYTHC08UwMa90qvrt1fIbVOI1XOetAJ2eeDPx8CZrfw1/xEwHhiJZF4pkRPN9TuBcwoUpePhRZptjka+RNOAPwKLiY/IAcixS0RmVDrsdMOiErh1Pfx6I9zxLayskM7RE7rCyRkwZbPcz2qRxovbqmF8l9iNx+cXY/ETBZD3b3h/hYicLp0g9xR4+IfiEzl1SGKKHJfLdfiC7/F4cLlcUe0r6HVxu904nU7y8vKi2mdj2JPgtgliVi6phBcXNq0jeyyf83EBM8yGqoh3ETO2BG7bSG3FDouakRUlQVgJXA8UBtbHA88BQ2J8nM1V0mOql722r9QlXaW43z0bJIpzeQ7U+OGrcpiyCR4YIFkvDw+GX66Hm9dLdeNvqyDNGhtjaHm1RG8WrpGO3CDRm+MHwLnHwnH9pMhdW8IYqC6DqjKorgB3OVSXg7sSvB7w+0IWP/zjTRfPvZrHLdc6ufVnDp57WYRKVSlMczlIboah+2gVkEPXY0mnFLjhAnjkv7Bim0wnjhvV9HEG1yMdX78Uud0Rnxm6ZrE+cNsGW3B1KFToKEobpwYp8jct8Hc2EsW5nth2GK/2w/StkrHSPwX218Ave8Nl3WUqalxnKex3VU5t5dlfAj9bC+8Uyy/pLsnw5DBYUSG1dU7OlJTzaNhzCD5YIdNTwYJ+malSy+U7o6T2TWti/FBWDAd3yVJaBOUHoOwAlBeDr4m1/wuWuphTmMekMU5OSnPw2b/hOKuDSWPg0T/ksWoBXPEdBxldoUsfyBkAOYOg+0BIClM7oKkVkONB365w2enwn8Xw2udSUbpH9pH3i4cYU6Gj1EeFjqK0YVYCPwWWBtYnIh3H41HZeG6RTDPNHgGdk+CVfeK72e2RSrMTukprhsFpcv9gDZ1Lu8GsXZAVmCbqmyLLhCi6XRsjxfzeXyHtGIIM6A4XnNB601K+GjiwE/Zvgf1b4eBOOLgHvI1dVC2Q2glSOkFqRuC2EyTZwWqrXQpLfdw41Mm1kx143RLxqamGa/s5sCVDjduHuwLcFXBgO2z4QnZvTRLRM+AEGHoaZOXEp1hgc7ngBInorNsF//wE7rn0yPvEQ4wFhc7ONiB01gVu42BPU5qBCh1FaYN4kShOPhLF6QL8CfgxsY3ihPJ2MQxNgxGBaM1NfaBbMvxlp1Q8HpYGGWm197ciYmdjFfS2g9svXp1oMpv8fql58+5y6akE8nxHD4QLR0sTzZbMnKoshT3fwu71sG+ziBx/mKZF1iTo3FMiLdk9ILMbZHST205dwNaEb9rTL89v8H8/cDjwVEmEqLRIhM7+rbB/s0yL7d0ky5K3IGcg7N9s4y//iE+xwKZitcB158LUf4vYWb8Lhvepe5/GxFak02pdAq91aSs3l6pCsq4ATmjNgSgqdBSlrbEOuA74MrAezyhOkGq/ZFOlh3hcUq1wYRcp+DdzOzwTEn8PRnMWl8Dyckk5TwuJsDQ3olDjhc/Ww/xvpPAcSMTm7JFw4QnNK+oXTTSjuhx2rIFd60TclOw98j4pnURM9BgM3fqJuMnKkahMPLGnQde+sgw6UbYZA2VFsHMtbPoKdq0VAXRiqkx55eXlUXEQZvwxtsUCm0r3THkPP14N85YfKXRiReh7nh54H6oC1a9bKoJVn+VIdfIeQN8WPbJSHxU6itJG8CP9qH6D/BrMRqI41xK/KE6QVKss292wx1NrQu6ZDFd0h8d3SP+qkzJkeuu/RdKl/MtS+GUf8fGE0tSKtx6vGIzfXQ6lgSyZTilw/nFw/vHwhxn5bPiweaKlOdV2/T6Zhtq+Gnasgn1bkAJFIXTtB72HQc9h0GMQZHZvO/V4LBYRWVk5MOo7UFUKm7+GjUsgF3nuv3s8j5l/mo7XF9/U8oa4aLQYyFdth+0HoH+32B8j9D2/c4o8P4+BaU4X+VNbJoJVn+B086nE//OrNI4KHUVpA+xAzMXvB9YvBP5GbNNSl5XXzaQKEqx3c3UPuGGdtHf4bleJ2FgsYjzungy73CJ0+qVAhU+MyQWjxc9Tn6O1H7jvfgfzlkmvpLJqeUyXTtJS4OwRkBIw10bSIiD02B9++CHjx48/vO50Oplyv4Mty8DpdFG828f3RufXeXzXvtB3FPQeDr2GiZ8mUUjLgmPPleXQHhg+z8HbvxCRk2S1M26gg+py8QpFQ3OiZjlZcOpg6Yf1yRr40TnRHTscoe+51wC5DnbPdpH/bMtGsEIJWKg4tcWPrBxBjKs0t2m0BYTSFnnFGNPZSJn4NGPMU8YYXwz3/2WJMZO+MebSb4y5eLkxD20yZk2F/M/vr3t73wZjfrzamFXltY+v8BpzzlJjCkLaOFQ3cYD12w/kTXWagq+Mufv/6rZmWLjamBpv4/tobouA4P2Cy69+6TTv/sWY2bcZM2mM/G/SGKf5v7uMmT/LmDWfGFNe3LTnlSgEX4PkJPvh5/vP+43Zuyk2+63/HjS0/Zut8l7f92Jt25B4EDy+JTlwvk2LfbuLpjLQyGf63VYbQfsimuu3Ch1FaSUOGWN+Ymp74ZxujFkWbVlYAAAgAElEQVQb42PsdRtz/Rpjnt5hzG63MZ8cMuZHq4y5ZZ0xWwPNd7x+WYyR/lWXrzDm1xuM+bpMts8vFvGzOcJmPUGRk5xsN/f8vVbgOF4x5rN1xnibIJoi6ddU4zbGHrjgBS/ys26sFTk3XOU0O9cZ42tAYCU69UXHb+6qFXfP3WLMig9qBW4s9t+YAPV4jbnjb/K+b9ob+TGbQvAcsSTb4yqqGmOzkc90kjGmrHWG0O5QodNEVOgobYVFxphBRj6AVmNMnjHGE4fjLC4x5qyvjNkSIlI+LzHm5nXGPLCxdpsvROwsOmTM3d/K465fY8zYr4yZvSuy4+fny8XPFogojJnkNA+9YswX3xrja2bYKngBs9vtDd7H7zNm51pjPnrBmMvOlGMnWWvFTjCyMa0Vf+m3BA2JjqkOZx3R98Hz0Qm95gjQWfNF6Ly+OPLjNXU8wYhOPBqYNoX/M/LZHtsqR2+fqNBpIip0lNamxoiosRr58A02xnwax+PNL5ZO45sqa7f5/Ma8ts+Yid8Y89mh8I+r8hnzRYkxb+43Zn8ECsznN+bGO52Hxc2Ns4z5zg9kPT8CkXG0C2rpfmMK5xjz0oOmTtTm++c4zef/MeaBe2qnsRoTSu2FqVOnNniRnzbNaW65bqr5683yWr3/XHRipykC1BiJ3t04y5jfvRn5sRojeI7c5XCaUwuNGXxL9J3fI+UaYwxTp5pxDRzb6XSaqVOntuiYEh0VOk1EhY7SmmwxxpxlaqeqrjPGxPtM3FgpEZkP63lPtlQZc+8GY6ZulvX9HmNe2G3MyvIjdtFsVm03ZvzVtSLnvheN+WilTF801V8TSkNTJNPynWbjV8b873G5YAeXYCTnvjudxu+ru4/g0lq/9NsSW5YZ89wt8pp98q/I9tGciM6OAyJ0fvW32Pt0Qs+Rz0uMObXQmKtWRna+RYvHGJNtjKGZPialcTqE0HnkkUfMmDFjTEZGhsnJyTGTJ082a9c2z9GgQqf1aeyXZnv+lfOaqTUcZxljXmrBY9/5rTE3rDWmvN6v9j9sM+b29fL3klJjLlhmzLM7Iz/Otv3GPP4/uZidculUc+ZlTlOw1JjqehGh5rzP4S4KZcXG3PSjulMws24y5n9/NGb9YmMcD9U9x0L34XQ6zXnnnacXmgCblsprN+tGY1Z93LzHNtck7vUZc+tzcn7sj/FXcOj3ytwiETq3ra8dZ0t+r3xk5HPe3RiTH6GRXjmSDiF0Lr74YvPCCy+YlStXmmXLlpnc3FwzYMAAU17e9J+gKnRan+ZmayQ6lcaYm0xtFOcMY8zGRh8Re9ZXGHN6oUxXuUN8Mc/vMuZ7y2u9OR8fjMycerDcmBc+MuamgMn4lueM+fenxpRFaF4OJXgB8/uN2bnOmHf/Ysxfb6qdnrps7FTzxRsydRWOjna+RcLSt+X1/OvNxhzY0bTHRPq6uv4j58jXm6McdCM8v0uETjBa2dLcbWojtsZEZqRXjiSa63fC1NGZN29enfUXXniBHj168NVXXzFu3LhWGpXSXI5WX6U16l3Ei7XAVcAKpGDY/YATCNN/Ma4ckw4/6wXP7ZK2DZd0lYqtKyvgkm5SLwekaWdz8HilkvG7y8AdKLc/Zghcfjp0b0Yl48Zw/DafjYXw34el7UGQ3sPhjzc4GHRy4+0VWrKxZVvoLxUJJ02Q9hZbl8Nn/4bcu49eEDHS17VbphQNLK2MydDDsi3Q4yrY86ol8QOvBv6+PHDrcDiYPn06Ho8Hu93err7jEoWEETr1KSkpAaBr14Y7B7rdbtzu2s5upaWlcR+XcnRCxU7wC6C9iZwXgVuACqQE/D+Bi1pxPLf0hVKfNN98swgOecFuhVsjqE1vDCzZCP/9Ag5WyLbBPeDKsTC0Z2zG666AVQtkqQp8bG3JMHwsHH++tF0IR0lJCTU1te3Cb7/9dgCKioqOuO8tt9zS4P8iGrPbzYwZM6isrOTee+89vH3mzJnMmDGDBx54IGbHijUjLoA1hbB+OfRZIA1CGyPS19VXJW02du2Boii72jfE6l3grYKu2dDSL/fnwE4gEzgNKAIef/zxwyLH4/Hgcrna1XddImAxxpij361tYYxh8uTJHDx4kE8++aTB++Xn5zNt2rQjtpeUlJCVFaOfnErEpKSkHP4CCBWkiUwlcAdS1RhgPPAvoFcLHNsYqDEiYMLh8cOmalhfCUkW+F4Epfh3HIBXPoNvd8t61wy44nQYMzQ2bREqDsGK92HNQqgJnBLpneG486TFQWMVfUtKSpg1a9bhH0GxxvjBVxFYysFfAX43mJrA4oUl2xdRuHshp/UZx5j+51C4axFLti/kjCHjOPOEc0jKAmsmJGWApY39zNy1HvZuAHs6jBoH1gbOo2hYsxM27xNhPCoODaCMgRf3gtfAFTnhq3bHk/eAQuB4YBKwaNEiFi5cyIMPPsjDDz/cbqPXLUFpaSnZ2dkRXb/b2Eetadx+++188803LFq0qNH7TZkyhXvuuefwemlpKf37x7KovhIpLper3f3KWQ/8gNqpqnzgt0Ccez0CUFQDD2+FDBu4Boe/j90KI9NlaS6VbnirUJozGiMNNy85WfoY2WPwLVKyD5a/C+sX13YH79oXTrwYho5pWsPMmpoaSkpKSEtLIz09gicZwO8D7wHwHABvEdQUybq3jDp9sKyBJZRzrJNJzkzn813zWLrrU3z4GJs5gTPcF0GhTBn6kI70yd0hZQCk9IPUfmBNhXfffRer1cpFFx0Z/5s/fz5+v5+LL7444ufWGJ3HgLsIPNXgL4OcIbE/Rtcq2F0OqenQLQ49rw7VgKkCuwUG9aidlm0J/Mh3AMAYYOn8+SxcuJBx48Zx9913A3Wj2du3b2f//v189dVX7N27l86dOzNkyBDOOussZs6c2XID7wAknNC54447mDNnDgsXLqRfv36N3jclJYWUlFaYqFUapf6vmuA6kLBi5z/Az4EyZKrqZSSa0xK8fxAe3QolPki2wI29oX9qbPZtDHy+Xqapgj2pTh0CPzhTojnRUrIPvn4bvv1CIiYg/aVOmgD9j48sSpSenk5GRtMH5y0F905ZPDvBsxtMiMXECtgRoWixgS0bkrLAlgW2dLDYwWoP3CbDAC7ny6ffx+f3YrMmcfn3L5co0CGoOQTeQ2DcwEHwH4Sq5VBlhZTegDuVeUvnYLfbyc3NPTyGgoIC5s2bx6RJk5r13JrL0BNg8zIo3wtDR8d+/ympkJwK6RkQj6exuVzek74pkJ0Z+/03xhqgFEhDhM68pCQmTJjAKaecUud+DoeDdevWMXv2bM4//3wee+wxevfuze7duyksLOSVV15RoRNjEkboGGO44447eOONN1iwYAGDBzfws1Vp04QL3YYzKCcKNYjJ+PHA+neAV4AGLCQxpcwLj22Hd4plfXgaOAfHTuTsPgj/WlQ7TdW7C1xzFoyMwZRDOIHT/3g4+RIROvHEVwnVW6B6s9x6Dx15H2sqJOeAvQck5YA9B5K7gjUdLEeZ0ikoKMDn95KUlITX62VRcUEd0WKMTHtVbw+MYSt4i0VonZaUi/0SmDNnDgC5ubkUFBQwZ84cJk2aVGc/8SBnsAidkr3EpPlnffyBiJg1TpGWHUEjsr3x+8WDzwO3pyEJBxMnTqS8vJwDBw4ccd/t27czZMgQ3n33XZKSai/D11xzDY899liLjLcjkTBC57bbbuOll17irbfeIjMzkz179gCQnZ1NWlpaK49OaSotmQUTb3YjWVXBCdT7gEdomQ/V0jLI2wJ7PBJxuL4X3NAbkmPgq/B44e2v4b3l4PPL1NSlp8CFo8EW5f7Li6FwLny7uFbgDDgBTrkUegyKeuhhMUaiNJXroGoT1OypdwcLJPeQ6SN7X0jpC0ldIosm1RclwXXgsEixWMCWAZ1GyQJQUwJVG2TaZeI5uViTRey8/fbbeL3eFhE5IJ3as3uK0Nm/RcRnLPEF3vN4+H+gVuj0beFAfjXwdeDvM5tw/wMHDtC9e/c6IieINV4vTgcmYYTOM888A8B5551XZ/sLL7zAz372s5YfkBIRjaXXJlIkZxFwJbAHyAL+DlwWo303lqY8zenii0M+9v84HwP0tYsnZ3SMfnmv2Qn/+gT2BzKdThgAPzxb0oKjoaoMvn4HVn9c68GJp8AxPnBvF3FTuR589RIuk3NgUeVc7F2tTLomF2u9CEBBQQF+v5+JEyc2+ZjhIi/B2/pipz7J2ZB8au16bm7uYZGTlJTUsMgxkFSdSnJlOvbyDKw1yRibD2+Km6puB/CleJo8/iA5A0ToHNobe6FTE3jvk+NgXKvwwf5Awl2soppN5UvAg0xbN8XaNHbsWGbPns2vfvUrfvzjH3PKKaeQnNzShSc6DgkjdBIwOUxphxjgaeBuwItkV/wXOCaGx7DZbGGn8X491cVMZx59bnbSG5jYDX7dHzrF4KJR4YbXF8On62S9cye4+iw4eVB02VSealgxH76ZX5tF1Xs4nH459Iyx2dUY8OyDg8vh4Ebwh9RqsdghbQikDYe0wYGISoGVOXPmkNyNI/wwQcHSHPx+f9jIS3Dd7/c3eV8FBQWHRY7X66WgoCCs2Om2fjg9Vh5P6qHO+G0+jM0HxoIvxU1V12J2nrYEd+fmZaFlBCp2VISZ0ouWioDuSo9DxGVTlXw+eyaLKb+lMMCCwN/jkESEozFjxgzWrl3LU089xVNPPUVycjKnnXYaEydO5Pbbb4+rD6sjkjBCR1FamyrgZuAfgfWrgdlArL+SwnmWfjLFxb9miMgZcZODBwfChV1ic7yvN8NLi6C0StbPOxYuOx3SovA5+H2w7lMonCPRHIDuA0Tg9B0Vm1T0INVb4MDbsHEOHFgO2TZIt4E1DdKOgfQRkDpYjMKhhIu2BEXOiBEjwgqLxiI9jUV/mjPtVD8yNHPmzLARoW7rhtNn4RkUps+jy0TwZJRjrH5sNcmkF3Wn91en0vfL09ly/kf4k71NPn6nQOFIT5UIVXsMoyOVAUN7pzh4aDYGzt+hLexk2ITUzkkCzmriY7p168Ynn3xCYWEhH3zwAYWFhSxYsIApU6Ywa9YslixZQvfu3eM36A6GCh1FaQLbkUqnXyGemMeAe2jar7dICBU706ZPx+fx0OdmJ7l3OXAOhl4xuFCUVsHLi2DpZlnvmQ3XnQvDoiz6s2M1fP4aHNwl61k94LTJMOSUoxt5m4rfDQc/hKI3oHypbKvxSm2a9JHQ4xRIHSRZUo0RKnaCU0UjRoxg3bp1R0RRIo30NIdw018jR45k/fr1R4idtE8H8adtD1B55hZye9SO05fqxpNZTkWPfRz72pUkVaXhSS5r8hhsyZCWAVXlEtWxx6gIlDEhEZ0YTy35jNSIAhjWwkLno8Dt6UCnZj52zJgxjBkzBpDyCPfffz+PP/44jz32mJqSY4i6nhTlKHyGZFJ8BXQD5gP3Ej+RE+Tqex1Yk+34PB4syXam5Tl4ZnhsRM7STTDtNRE5VgtMOAkc349O5BzaA+88BW8/KSInJR3GXgVXTpVaOLEQOdVb4O6L87l9lIstjoDIsUL2OTDgQfh20CIKre+SNvToIidIbm7u4SmipKQk7rnnHiZNmsScOXMoKCgAwguQeBBu+is3N/ewuFq7du3h8ZhKK4NG921wPBa/FSwGq6/58zgpgSt2TVWzH9ogFdXg84kRuVO9qau5c+cefq3rU1BQwNy5cxvd99ZqqPJDuhX6tKARuQgpEAhwfpT7Sk5OZurUqQCsXLkyyr0poWhER1EaYQdSD8cNjAbeAgbF+ZjGwH+L4O48F/4aETmmxsPu2S5sURq2K6qlsvGXG2S9X1f46XkwIIoouacalhZIRWPjF0Fz3PlwSq5k8USL8UPpYtj3CpR+BtW7bTy7Kw9LEvz2PgfdJ4G9J0yZMpNPFi1kwoQJzdp/Y36Yls58amj6K3Q8t912G16vl8nn3Ml303/EjgOLcWeWYTEWMBas3iSSq9Lo++XplPfcize1utnjsAauDP4YJkEWB9qFZKcfmb1ntVrDTs81NYq2KrDvUZ1atkjgfMSjcywwoBmP2717N7179z5i+5o1awDo06clClR0HFToKEoj9AMeBL5BMqticN1ulAofPLIV/u/3LnY9m8eZdzp55w8Onno0+qKKK7bBiwuhpFI8MhNOhNxTI8+AMQY2LoHF/4HKgN91wAlw5pXQOQY9r3xVcOB/InDcWwMbLXDvVQ66lsEjz+fRaw84ekrRyRkzZjBu3LiwVYUb4mjp4E3KfGoh6o+H7+/E8v4ohv/vUip67KMmrQqL34qtJpmUkmz8yTXsPPtTvGnNFzq2gJ/JV9P4/ZpDcbnchis02Zhf6mgC0+uHdYHI07GRF8RuNmXAp4G/m1ur+uKLL6Zfv35MnDiRkSNH4vf7WbZsGTNnziQjI4M777wzxqPt2KjQUTokxUDD7WDr8lDgNt7zvBur4DcbYfHTInKuuM/Jf37nwGKJrqiixwuvfQ4L5cciPbPh+vOl31CkHNwFi16G3YGa91k5cNbVR28G2RS8pbDv3yJwfAEBZe0E3SdBj6ulZcLDOEgdWLcx7AMPPNCsGiRHSwdft25dkzKfmsrcuXOxWq3NNjmH3id0PG999B9yJ1WRvWUg2dsGkFwl5pSatCpKBmyjeNiGZpmQQ7EFxK8vsoeH5UDAJtRQRe3mRtGCr+cx5+fi9kumVbBjeSTlAZrLx0jB0AHAiGY+9qGHHuKtt97i8ccfZ/fu3bjdbnr37s2FF17IlClTGDVqVOwH3IFRoaN0OH4BDAUeoGnipSWMbPOKYfpWqPZDOj5uecjJX1zRF1XcVgSzP4C9AcFwwfGSURVpfypvDXxdAMvfk2kNWzKc/D0YfREkRVkGpKYI9r4E+/9Tmxqe0g96XAPdJoKtXjjN4XAcFjl2u517772Xxx9//MgdN0Bj6eDr1q1j3bp1jRb+ay7RTM+EizwVzCnA7k/joonjKRm0NfwDDRGZyYKRHFuMSrv4/LCvFLZ+M5fUA1ZGXN6w2Av1SzX2Wgdfz6EVwBm5HJcufrOWMI1XAx8E/v4uzX+Jr7rqKq666qrYDkppEBU6Sofie8AypIppqICJ8HoQNTV+eGIH/Hu/rJ+eCdP/lE/XBi4wTY3k+A3MXy6NOH1+6JwOPzsPRjXeHq5Rdq6BT16C0n2yPmA0nH0NZDazOWP9gog1RbD7Bcmgem6rC5/xcde5+fS6Hrpc0LCpuH5j2JkzZzYrotPQr/2CgoI6IgeaXvivMSKdnmko8vSd0qsYtOkUPvrvC5x4xTDwh5zBFiMndIQntTeQHZUUozTwojIxIicnWZk/bw5p9vBib8SIEU2OouXm5uLxw7z/zSG7BkZf2XLtMj4BKoEc4NSj3FdpfVToKB2G8UiWxFqkmvEhpAFfDVILp6XFTlEN3L8RlgeMlNf3gpv7RG+mPFQBf/sI1gXSu08eBD8ZBxkRpvRWl0u6+LeLZT29M5x9NQw6ObJ6OMGCiL4quKGHg32vSZPL2btdPLsrjwd/4WTUc43vu6HGsOPGjWPy5MmRPdEAsSz8V59ITM4NjWfgab04uHgHVQTML9bYFVWtCQid5BhlMO05KLdnnZvLcf0brl3U3Cha73NzyS6FkoVzmPZZy5jGaxATMsAENHU5EVCho3QIrkWqly5FRM5s4H/At4jB+PfAuS04nhXlcN8mETudrNKM89zO0e935TZ4YQGUV0NKklQ3PmtE5AX6Nn8Ni/4VKPpngWPPhdMvA3sUtUoevNtB6Zcw7dE89vaBG3o7+LtfRM60aU7y8upGrepHgEJFTvD/+fn5VFZWMmPGDHbu3Mlxxx0XsRcmVoX/GttHc0zODY3n4LCNMAzO5KRGHx+JN8gTmDpMipHQ2RFoPNurC5zTSO2i5kTR/Aa+LoPsc3Ip+7TlTOMLgRKgM03ra6W0PipGlQ7BVYhp8FXgdmAqUuDrJ8Aw5JfZe4H7xrvZyJtFcON6ETlDUuHFUc0XOfn5+bhcrsPrPj+8/gU8NQ8Wvu7i2/fz+e0VcPbIyEROdTl8MBvmPysip0tvmPwbOOeHkYscf40YjFdMhmt2O7i5j5Nnd+Vx1jcpPPW1CJf6IgdqI0DB5xtsDAtiRrYFnLP33nsv48aNA6hTAydIMELQ2k0Tw6WzxwW/BYvXhtVia9brUeOWkgEA6dnRD6O0Cg6Wy3nYLzDNWb920THHHNNgFG3SpElho2ibqqHYC+WfFuD3tcDriXhz3g6ODY0UJAr6PintGj+i5iciLRxuAHoC/wLOC7nftcBdiHcnXvXGvAZmbofXAn6c8ztD/qDIelWF9sO69S4Hsz+ETXthaYGLwjl5TM130jPCCNHmpbDoJRE4FiuceDGcmhu5MdUYOPQB7HxaGm2CtGR45DEHf7uk1kzckP+ofsZZUOSFTl0FOeecc+jWrRsff/xxw16YCbmwF+nIugeZw6wMWQJTiaTUW7oA/QNLTyL6mdiU7uZNxVadgs1jx5MVSGcKKnRjAashtSSbrO39ufKMLmAxTfYGVQZ6XKV2it5gDrA1cL736gypgf3VF3sNRZyC4w3Hl6VQsqiAkoWxeT2bwgdAOeLNOTvme1fihQodpV1SjoSXS4FgouZVyAm/ConmhHI2MrV1CLmGxZpDXvHjfBWwU9zcB37eS7JEIiH04v/O13DCBAffzBORU//i31TclfDpy7DhS1nv0gfO+ynkDIpsjAAVq2H7TKhYLutJ3aDPTZIqPv3RumZil8vVJLETzLRq7HnW8cIUvI3X52VSv0nkfpULBUC0hfBSgOHAGORkaoLoiaa7eThyVh9Ll82DWfP9/8qG4LlkEcXjTXGTuasPtprkZnmDKgIZerGI5oBk/gEMzJHbWIi9fR5Y9r6InItyY2sab4gy4N3A35OBFuwbqkSJCh2l3bEeSSEvAbYAPwP+FPjfFcBFQLCuWDDiUwmcBGTGYTxbquGuDbDDLSXqXTHw4/j9cOIlDsYUwudv5PHl3On4vI1f/Btj13pY8AKUF8sUw4kTIoviBP00U+5wsPNpKHoTMGBJgZe7uUg73ofziiMjMsH1Dz/8kPHjxzf4HGw221EjQJQCyyB3Vy5vIyIniSRyd4Rc9JIRRdsL6evRCTkpgosVmafwIGWxq4F9wDakg6MbWBFY5iFXvpNo1M0ea5OzN7Uai99KUmUaVp8NqzdJlppk2V6dSnJlOqkHuxw+TlO8QcGsuoymFppqhOJymbayWmXaKlZi77MSMH4/Qy+YxA8mxd40Ho63kLe9P5pplWio0FHaFSuRnjM/QwTNQeBHwEjg1sB9QsWMFalu+ntgCrUCKFYsLoUHNkG5D/rY4Y/Dom86WF4Nz38Iq3fAKbkOlr0zHW/NUS7+DeDzSofx5e8BBjK7w/ifQ8+hkY3NapEptd0vwC+6yli6XgJ/r3Hxh9/l4TzNGXbaKTRis2DBgjrboNaADBwZATKI4HgF+D9gtzymgAK8iMjx4qXgmAJyL86F3ki1yEitOn6ky+s3wIeB4z2L9Aa5HDnZwhBrk3NNeiWpB7sw8s3LpP1DAGOVC7w/yYu1Jpny3vKCNNbq4vBjjfQsA+gcg2aeGwP76tdVpq1iIfaKamB1JXQeN5HrjuyiUGd/sWInsCjw99WouTXRUKGjtBuKgd8ifpvfh2xfgFQxvZXaCA7AGuAfSAbWbcCvYjyeV/eJJ8cHnJQBvx8CXaL0PGzZD3+dDwfKpXVDzQrXYZFztOmf+hzaI4bjAwHfzIizpQmnPcI09IpVcOUGB7v7wDOb87B1gkf+4eCJ/7mYNr1W2OTn54eNPAXXP/zwwzoVoENFTv0IEB+BY4dD0udCKMguYE7JHCadNYncH+ZSMD8wPTIKck+I8iJoBQYGlgsQF/sHSPjwceBkRGnHuEN3fYzVj7H52DWmEH+S9/C6sfoxVj9+m8zPeVOrmzxdVFkiRmSrDTJzohufzwebA/6cYLPYWIi9z0pE2w5Pi02D26bwBnLMk4FjWuaQSgxRoaO0G2xIjYsx9bYPQooEhmKQMHQF8BhwfQzH4TNSBPDlwBTApd3gwQFgj/Jn4Ofr4Z+fgNcHPbKg+hsXTz955PQPHL2w4PrPxXDs9Uin6nHXwuCTIxuXrwp2PQv7Xgb8cOMwB9lnw+9ey+OvZx7pp8nPz29wXw6Ho85zCfpxICBy7nHA38HxYSAC9FHg+aY6YBxgg/lV85m3YF7MC/6FJR24DAkjvo1Ukvsa8CLqOY6FmYzVj8/uoXjYhkZr6DRnuuhQIBqWnVPbBiJSthaBpwbSU8SIHAuKa2obeJ4dIw/R0ViLBAwtSMBOSTxU6Cjthmzgr0gjTpBrTRJSFDA4XRXUGhbEUjGD2E5XVfvBsRk+CmSu3N4Xftoz8jo2UJs6/sEKWT9xIOz53MUfHm14+id0PZQatxiO138u631GwPk/h04RXohKv4StD4Nnp6x3nQD97oaTuzl4POXoGVUNEdrewWq1kn9HPo6DDuiLmK8Ah8UBg8E32gcvIj6ax8FfGL+Cfw2SDfwQMSbPRK6MC4lrcSZPRjllfXZh9dnwW73h6yIYC8Znmvx67N8mt137Rjc2Y2BN4Jw4prd4dGLBJyUSlR2aCn3ilR4Zgg+ZEQV5K+ORqKDEHxU6CUj9AmqhuFwufD5fo7+a2xPrEO9pb6SAV1Dk+KjNivAh0ZsgzyGp5rcRW5FTXAP3bISVFZBskdTxi6M0dFa44bn3ay8auafApaeC83Nfo9M/4fphFe+CD/4KB3eL8Dp1Ipx0SWQXIW8Z7HgCDrwl68k9YeCDkB3Iua3fnqE5U2p1Hp9kx+P1wJMh/xyKhOCuBceAkH0Gar9cfPHFZGSE7xwZ9w7kQxHH+2uBZQRieI4D7s4lbL7gw9oN4cS0xXDppEsb3Efo61FdAaWBqabuA6Mb20KIlVMAACAASURBVN4SMSHbrCJ0YsF+T200JxbFNZvCx4gFqxMQv85ZSrxRT1UCUr+AWpBguN8Wbcw5QZgJXIhcV84Gfo4InyDB7/3qkG2/B24CBhPb9NBt1XD9WhE5WTb4yzHRi5zdB2HGGyJy7Elw44UwaYykpOfn5zeail1f6K7/HN54REROejbk3g2n5EYmcko+g9VX1YqcnCvhuFfripyg2djtdgcKAR55vjaEa1rg8TlO3F43TpzkkYfrGJdMD61HzFgDmj/25jJ37twGC9AVFBQwd+7cI/8xHjEk1yAGsGhT2RvB5rZj8YV/Ey0+K9lbB2AvbVou4f5AX9DsnpAS5S+ANTvkdmiv2to5zSHc674w4M1J/qKAwvfCvO4xpgyYE/j7MkTsKImJRnQSkHBTFA0VUGuvvAm4gBeQKaj3gf8g3tAFSLVjN1LupAoRNk8g18cPqVssMFpWVcCdG6RWTl87PHkMDIrSiLpmB8x6H6o80C0Tbv1ubVXZ5uDzSp+q1Qtkvd+xcP71kJYV/v6NRQudU10c/NTHTw7lA5AyAAblQUZIB4KjZVSFrh+BF1zXuMh7PQ8nThz7HZAFjusDj38yDwrBcUnLnd8RdRy3IhGnaUh21krgxNiPzVadwon/uI6iUWvYftZnGJu/TsM2i99KzqrjcGeVsv2cTxvdlzGwb6P8nRNlNKe4HHYGWj6MjHAKrP7rvr0a1lbWFgg8Lo6dyYO8inx39AfOifvRlHiiQidBaW4BtfZGMI08aA78JWKP+E1geyG18+l+JLtqLuIVPSOG4/isBH6zSbw5o9LhiWHQLcrMqk/WwEuLpJfPsF5w80WQGUFKemUJzJ8FewMXsFMvlSiOpZEoTmjF5dBzyXGHi+l/zuPmPk7oDTlXQ787wFpP0AXbMzRnSg0/YoSYBr71PhE5XR1wD9KvIxscOKBbA4+PI5F2HKczckIuIG5Cx5/kBSBrRz8Gfnwu277zCf5k72Gx40/2su/4lfT/7Kyj7qt0nxQKtNmiKxAJsHyL3A7qAVkRllIIfd2NgX2n5h5RBTmerAS+RDTjj9Gpj0RHhU4CE2rYjMTwmcjUIF9G1dRm8Z6ImIvvQoTPa0hEZxBSvX8JMCSGYyg4AM4tMjNxZhY8NgTSo5gP8/vhv1/C/G9k/YxhcO25kkbeXPZshPdnidixp4nheODooz+ufvTloQcd3P99F79/S0TOLaMdDJoKWQ10MzxaRtURvAfcz+G0uPxu+XAvInDqzbi01vkdScdxAI5DhM4q6kRaYoWx+sFi2HDxuwz54AIGfziereMW4k2rnayt7nyI5Kqjq41d6+U2ZzAkR5GyXVQaiOZY4IQopxbriMy33waflwmXxl/kVCMtYkBmIQfH9WhKS6BCNYEJZ/jsKIxBBM6bSHZVkNFICZNtwObAtl8AW4mtyPnnXpi6RUTOJV3h8aHRiRyPV6aqgiJn0hi4/vzIRM7aRfC/mSJyuvSBy6c0TeQEcTgch301KfaUwyLnvuscHPtKwyKnWawAvgtcjIicLGA68qZNIT4lqqOgfhPKJl1sRyBGsANIb61YE0gp96W4WTv5LZKqUxn84XhSi7scvkun/TnUpFU12qnWUwVFgWyrPsOjG9I3gf0M6QHZMXD6X3hJLhZbEvi8WG1JXD4xzmZypENIMVIwe3Lcj6a0BCp0EpRoDZ+JzkTEi/oIMk0V/B63Ad9Dal+sDmyzAOFzcJqPMfD0TqmTA/DjHjBtECRH8Ukqr4YnCmDZFkiywi/GS3ZVJCnp5cXw6Svg98HgU+Cy+8Vc2lzuvNhBssVOjd9DssXO9FkOBj8CSdHWLikGbkGMVfORVgx3ARsRA1UbEzhBIuo4noKkAwIUxW9sVm8SfnsN6ybNwW/zMeT9C+n19UnkrDyOfovP5MDw9Y0+fsdqOa+zcqJr+7DnIOwKtBA5PkZG8WdfL8D4vFhsSfh98e1MDvIDaX7g72uIX4NfpWXRqasEJCrDZ4ISGvkPpo6/wf+zd97xTdXrH3+n6QAKFErZlI0IIqjgQgEZThQ3zuv6KW6viAOVlJK6F17XBdfV67puBFERGU5A2XtvKKPQSUea5Pn98SQkLU0zTlIKnvfr1VdOTk6+OUlOcz7n+T7P59HpqvuAJ9GWD6Ch52OBZlHeB7fA837dx+9uBTcbLJ3dVwivfA+78qBeItx5LhxTxZjiUXLBxE/9VBhwI+Tv9uTjhCmWxAXZ72jlU7k4SIhTsfPGoixsFxo4ptzA+2gSlfekfyU61xjNUFsMiErH8RgZB7oSyrG4fSp7w3nTaDm/N6nrOxPnjKeg9Q72Hrci4OuXl0K2x1W6bY/I98MtsMATQu3SMvLcHH8+nTyVlT9NJqX/MG6/bCirZ8a2M7kLPUQFjRiHEQQ1qeWYQucIJKKEzyMUB3qQen/KBRU5TiARmIvOgDyK+uP0QyuxGhLdSgmnQOZm+GG/njNGt4XLDVrkb9sHr34P+cXQOBnuOx9aBbiiFgm9FLzzyZHtjzMPNj4O4ydlMWFnBg8MsvP8VBtPPh+643KVLEX7b3gLf7oDrxPd0rcYYbgJpTfUGCOhs/ms2TiSD/heQyC7zwKyey/AWpaEq05Ztc/fsVor8+qn6jRnxPuxR31zEqzQMwrRnG+/ncrMqSpyTj1nKJ3qQqcYdiYH7Uy+HfXWuiqqI5scbkyhcwQSdsLnEUgJ2oyzEJ3t+D9gKJpY7EYPXCdqSPs9eiU2Gf2x6o2KnmjhcMOjG+HnfBVZWR3gHIMeOeuy4bUfoLQcWjVWkdO4ivm1+ZMhf49GaI7pq07GqQZOSIEoXgsbRsEbC1XkjL7OztMf2g6Wm3unRsF3jAU1pywD7MCz6OVyMjAWnaoyWJlWUxhuQul9OEZCJ7/9loorLL7bYCKnvBR2eIyn2vaI3L273KXTrgDHtYU6Ueg/lV3mJqX/MNL6DWWIL+UoZu7WO4BvPctXoRdKJkcPptAxqXUUASejouYa1En/XeBLYCLaVK8c37kyFT133o9GgKI5r17mhoc3wO8FkGSBZzvBmQbzVJZvhQnT9QTRpaV65NSrYqdnvw/Za7UHVZ1kWDgVtq+AzqfoX7TI/Qk2Z4K7FCwNXNj+acf+sooZb7m53W7HbrcfjBb6T59WyRLgBjSaA+rq+DJqSnIEYagJZTng6XdGBB5IoZBQlExCSV1KUverjw7qr5NYVB+rI5E4l/XgFNY3074iLi7u4H5vXQ6uckhuDHOXTkWWuKt9v4FYtgWKyyC5DhwbBRFe5IKCUy4ixQ39UiCl0lkq2pGccuAdVIsfT3TtJ0xqB6bQMal1zECThz9Fr6z+D62uegP1tJgEeH9PV6P5q15fsmiKnFI3PLAe/ixUkTO+M5xi8FLvrw3w7kzNaTi+rbodJ1bxX7hns4qcc+/y9R3qcpqKne8/2ELpV0t46JlDTdOysrKYMWMGgwcPDtoiRNywcyLsekcfa3gavPJUJvF+79E/98tut5OZmVm9OaUT7ZKaiZ5B0oB/A1eE9TEdHWzHF8kyOM0ZiNZ/nQxiYVvfP3BZNYJTLyeN1n+djNWRiLNOKcm7m7P6kkkVTPgG9R9KtidHeVPBVKbPDmB+GITcA772JCd3gniDduMiOj1c6oaWicb/30LhazSiUx/V5jHsw2pymDCrrkxqHQXAOnTaysslwEh0/vwBNOqzFfXLecRzP5oUu+C+dSpy6sbBK12M/+j+ugremaEi5+ROcOc5VYscgPgE7SzuKNH74oZmHaDv1SBxJexZ2ZinbW/qY548EP8WIMFahLiKYePDPpHT/B/Q+WUqiBwvFcrNk5ICi5x1aGLU46jIuRg1O6olIieidg5G8PobtCdmZ8/EovoUp+VUmKYqTssBsZDXbgvbT52HO6GchOJ6DB06lGHDhjF58mQ+/e9U3G5Yv98ncsKNlIjAn+v1Nj0tMufuyqwqhjXFemIa2gSsMVYda9ELK4AbMaesjlZMoWNS62iLRmwW4ktxALVcGY76r63xbHcq0J/olY8DlLjg/vWwsAiS4+C1LtDbYNnzrBXw4a+am9q/G9wyUBseBsISp395u/S+oCXjDdPgVvuxdEjvztrfksjKysJiqTiVNGPGjEPsBvwfH32njbW3Q95ssCRC+3HQ5p9gqSa+a7PZDvo1VWlO+QWaHDUPPVu8j14q16J2z96IRmWx4004jotWi20vqz23Ma4qs4hHDXgEryvRARahqGU2B1rspqxhAVaHJs4MHTqUswcNY87SyXzy293MWxG50/C6bNjrcVPuE4X3WOTSaA7AGSnQPAq5PtVRjBYugOpzs8rq6MWcujI5bATqqzQAKN20idvq1mVuixa096yPQ4t37KipV2/ghSjvU6lbO5B7Rc7rx0APg938floGn8/R5bN7wuWnBk/8bNwSjj0DfvtYK2LadNcrZ7dLm3IOuLIJ8sV53D+uXZUtQAK1CHnoHzZW3wKOnRDfCDq9BPVD+IUP2I28HA2pjfdseCbwMbUiF2fKlCkVclL8q6XWrFlDly5dDoqfqLcVKEKjWQAnRm/YypQ2yqNOXiOspb4KqwY7W4FYEIvPJdArdNwuaF9vKHGW73BLGOaHlThQCos8EasT2ml+jhFE4NscKHFD8wQ4I8ahFUEP0/3o7OqVsX05k8OMGdExOWxUN8WyqVcvyhMTuQXfDADoD1RPfDk60aTMDQ9ugL8KoV4cvNrFuMiZvtQncs47IbDIKS3S7tHbV0JRrv7w9xkGnU6GaW/A+j/1eXF+ORBt0ptTr079gFGWylGYkRf5RE5SG+j6n9BFTpXmlA9maYm4V+Q8hHZMjZHI8eYHBdrH5557rsK6qiI4Q4cOpWvXrqxZs4bvvvsuNiIHtN+ICw07RtjYMhRyO24keXdzWs3vQ539jWmwoxWt/zwFR/0incICyhr4JoF3rIJ5y6aqyLGGYX7ohwjMW6/J9GkNoWsU/hkXFcGGUoi3wMVp1Uc7o8E89CuyALfiayMTS4Idv9VV05oYw4zomBw2gnVhvzY1lYFo5dUI1MRrLjAHTQOJJg43PLQB5hb4cnJ6GpwPm7YEvpqny0NPgot6BxY53zyneTn7dkDTtpqPc8ql0O86SKoHs95VY7eWXcAaD4u+hz2uXyk4sO/QKIuHylGYh4ZlcWszG8k9oNN4SGh86L5UJqA55WbIeNHzvaXY4D00kSqGBGo46t3H0aNHV9g+UEPONWvWEBcXh9vtjjiiERSPuOX06A99EIHC1jvZc/wyWs3vQ9qarohFONB8N9m9F+Cspwlee45fhrNOKSVFMGXKVJZumcygM4dx1T8iMz/ctEcdkOPi4PQuofs7BSK3HGbk6vKARtA0xlNWe9FoDqjDek31sgp2/AasYDQxjCl0TA4rwbqwLwRuQa1YitED9m10eitaOAUe3wR/FECdOO1AfoJBkTNjmU/kXNhbRU5ViMC01yGtrQobazys/FkjOz+8BoP+D/pepQLnr290fWJd2Oucx6iX+h/8rLw/lnCoYLyrt40xI7KYsC2DpHR4aYLtkK7jgajSnPK/YPvAY07Z3KVGgJ0i/KDCIJgwvvPOOxk/fnyF51TVkNMb0fFv5xBVsbMGba5mRX0SYoVHNOd0W82+Luuom9sYV6KDspQCfcBjJ17Uchfiho9fncqSzZM5pfswhl8fgfkh2q7kr/W6fHxbSDEY8XS5YVIOOATaJsEpMW4B4gAmoBZPnYHzY/tyFQh2/B4tHmi1EYuISPDNjg4KCgpISUkhPz+fhg3N/PraRFJS0sHoQ1lZRaOzMtSOJBftQh7NmRG3gH0LfLsPEiwqck41eGj8vBI+/k2Xh56kDToDUZgD0/4NZ14DLTrrOlc5bFsByzzlIANu0D5ETgeUHYDx/3oR27gHD/lx9P5oDho0iJkzZ2K327njOBtbnwYEPmqYxfhZBn5UBS0b9154XokmHUfB7j8cvO/TG6nyvp+cnBzGjx9PkyZNqF+/olK9++67cTqdByM5lds5RHX66iVU7AxAXS9rmspuzKLGgF9/OQVrXBw33TOUupUExdSpU3G7q/fRcbth+jJNQE5rCOf0NB7NmZGrUdQ6cXBry0M9c6LNf1FdXh8Yg/6exJKioiL27dvHyJEjSUtLAwIfvybVY+T8bebomBx2gnVhT0LFTU+iK3JE4IVtKnKswNMdjYuc31f7RM65vQJHcrwk1NES8t0bfeusCdD+BOg5BNxOWDdX18cnqrmbk8KALUC8pn52u50RHW1sfQoQaDocXvzJVsH0LyzKgH/gEzmjgf9R4yIHQqgAq4S3IadX5HTt2rVCgrK35DoqDSPXoSLHCpxnfLiIsFChnL2kEDYvgp7tLuLiyw4VOaCfQzCzwBXbVeQkWOGMrsZFzoYSFTkAF6bGXuTMQ0WONy8n1iInEOEevybGMYWOyWHFP3S7qYa7sL+ZDZ95GnSObQ9nNTI23l8b4INfdHlwD7j0lODVVfFJ0KKTOh4XVupw3a4XtOwKq38HR6lvfWZmZsAfR5vNxuzZs7kt3cZ2zyxOi5sg/SEtV7fZbOEnPRaiMf6P0BP4m8DT1MivR1UJnF5hbLVaqxTG/vhHbC644IKD01aVE5SHDRtmvK2AAN94lvuilt2HGbcb1swBlwsaNdcp0EjYWwBLPd0m+nSGBgYFboETJnuO994NoKvBKbBg7EIPX9BWMt1i+3LVEuzCziT6mDk6JocNf5FzrM1GR+Admw07se/C/ukeeCtblx9JhwsMmp0t3wb/maXnugHd4crTg4scEU1A7jMMJj0Dv30C/a/XqI2Xtj1gyxJwlkFiiHk1uz+EHa/pcqs7oOWtEb0lJReNTPyJxvu/RLuo1hCVEzgrT88NGjTo4ON33nlnhecGmpaqKgE3KtNWf6ARnQQMJ3+IG8q2Qx2DDTK3LoWCvRolPOZ0FbvhUlquZpci0L4pdGxmbJ9cAl/nQLGnlHywwQuMYJSi5txlwDGo0DlcVM7JqZxbZxIbTKFjctjwTrEMtNkYgv4QzQVejXEX9h/265QVwO0t4UqDP9wbdsGEHzWx8uROcPUZgUVOwV5IStZKKovFYwLYFC4cBd+/Aj9OgFMvg6btdFpr5xrdxhrif+qeT2H7y7psWOTsA84GFqHRiWlo6VsNYrPZmDlzJhkZGcyePfuguPHmIIE2eMzIyKC4uLiC6Z/hhpzhkA987lm+CMO9rQrmqKFjw9Og8eDIxsjN1n5WAMecBnUiSLAXgTlrtJdVg7pwapfIm396mZUL28sgKQ4ubwoJMYwMCppGtgttAHwbh28aI2AFI7G/sPu7Ywodk8NGZmYmy4F+qMi5GO37CLH7h/8jH8Zu0h/A4U01AdII2/dpF/JyF/RIh5sHQlyAE8HmxfD7/+D4wdq3qm4D9cVxu7Uj+aWPwvSJ8MsHEBcP9RvD7g1w9h0qjoKx92vY9rwut/g/gyJnDzAEWIb2aZqBdjw8DAwaNOigyLFarRVEjvfEMWjQIAoLCysIHUMNOcPlE6AE9c0ZYmyosmzI80yBJqRFNoajBNb8rsstOqtwjoRVO2CHp5S8XzdIMHjGWFMM8zy2Phc1gcYx7mI/Ha3ctAK3c3hbPFRZwYjvty5WF3YmptAxOYzsRCP8eajdyMfoD1KsWFMMj2xUH7dzG8OD6cauTvcXwavfQ7EDOjWH288ObHSWmw3zvoSEJFj5CzjLdSqhfmM9ibjd6nh88cNqDliQo9NaJ12oOTzByP0JTTxG+1a1uiPy98V+fCKnJSpyDmNSg/9Vr8vlIjEx8eB9/xOHt+qqxpmLRr0saFdIAwexqxRyvgbcUO9YSI6gL4HbDat/07yu5BToFGEUbleez/24T0dINWi5kOOAKZ68nFMbQtd6xsYLxgrgK8/ycGrEAaFaqsuNMyM5scUUOiaHhSLgQrTBc1fgW7RhZ6zY7dD+VSVu9erIbB848hIKJQ4VOXnF0Kox3HNe4Aad4oZd66FRSzUAXP07rPoZykvg2H6Q0swjdlwa4el8Snj7UjgfNtkAgbTLofV9BgRcIXABKnJaALPRxIZahMPhqD3eI7vxuc9diKGyQHHDvsngzIX4FEi9ILLvcfNiyNut053H9gt92tOfolK/vJxm0MVg5LPEBZ/vhTKPX87AGOfl7AbeQiO3fYmu75bJkYdZdWVS47hQe5FFaJ+ZqcS2QKXIBSPXw95y6FgHnu1oLC/A6YKJ02FnLqTUg3vPh3pJgbe3xEHrY+G4szRqc9IFcML5GrlZOh3279Dt4qxQUhDevhSvhfWjYMKWTD6sk0Xbhw89OYZsL1+GuhvPQ7+Q6dQKkeOf2+CN5tQKytGzaRnQBRWIBiiYCyXrwGKFtMvAGkFl097NaioJGjFMjkBQOF3qBVVWrlGc0wzm5XiTj/c7tYT8sqax7UruTT4uQfupXkvMmsebHCGYQsekxnkAmIL640wmtiFlp8DojbC2BFLj1RCwgYE4poh2IV+1A5LiNZITSki/YVNtzOnluLPg9OGwdRks/kH7XBXshR//DVuWhrYvjl2w/j5wH4C6ray8/HsGTzx1aCl2RkYGVqtvPqXKnjtu4AbImplFZkKmJh73CG0/Yklle3xvWS5QYzYEAfkS2AYkA/+HoV/Tkk2afAzQ+BxIiqB/1IFcWOvxXGrTPbK8HBGYuxZyi6BOglYQxhucT56VC5tK1ZDzyqaQHMP5aTfajSQbTT6+Ay2CM/l7YwodkxrlbeAVz/IHxLYVEMBL23zOqy93hlbVRF5CYdoSmLNWp71uGwJtI0gW9XqRd+wNA2/Wq/D538CUFzUJuV0IeRmuElj/AJTnQJ1O8OIs2yEeRIHs5Ss3U83MzCRrUBZZn2WRQQbWa60Hq6sOd7NBbwIn+HJyvI1FAWbOnHl4dmw2MMuzfBOG3OfKcyDnK0A0J6d+BN3OHSWwYja4nNC4hRpORsLybbB5r0Zw+nUz3pV8cVHF5OPmMQ7IfYNGir3JxymxfTmTIwQzR8ekxpgJeJ1O7GgHgVjyxV6fIWBWe+hu0JRsyWaY9KcuD++rvX6qQ6TqkL/F4hM7rbrCgJtg8vOadHzRqOD7IW7YPBZK1kJ8KnT5F8Q3CN43zEvlklbrKisZP+uy/VI7tvf08drQbNAbfaquLLdyM9OYswR1hQYYBlO2TCFuW1yVlVzBWiu4DqglgLtUO8qnnhf+NJHLqf3RSg9oJd+xZ0bmWrx5rx7jACd3huYG82i2lKqVA0C/FOgWY1PA34AfPMv/4PAnH5vUHkyhY1IjrAEuB5xoN/IxMX69+YXw/FZdvqsVDDTo975jP7zjZwg48Ljgz/GesKoSPN775WWw9EdokAoX3B/avmS/DXkzwRIPnZ6DxBa+x2w220GRU529vL9QSMTvMtsTTahNzQZrVVnuZnxZrmcAF0Dcd3FVNsX0NyysCnc57P0CnHkQ3wiaXgFxYc6ziBvW/qFVegmJ0GOg+i+FS06B+uUAHNsajjGYfJxbDl/u1fyc7vVU6MSSlcCHnuWhxD5SbHJkYQodk5iThxakeMvI3yW2yYE7yuCRDZr0fF4q3Nwi6FOqpagUXp+myZldW8FVfQNvu3Yu5O7QE0+7XtC8A6Q018eqEjzlpdqo88JRWk4ejLzZkP2mLrd9DOpXmqKoyl4+oNi5y8YTGU/gwEGiJZExY8cEjQYdDmpNWe4u4DU0Cbk7cB1gqboDeLBmoeKGfVPU/TiuDjS7CqwRRDw2L4a9WzWC020A1I3AKKaoBGavVMPL1qlwUofwx/CnxAWf7tUKx5aJcGET4yaD1bET7UoiwKmoX6OJiT+m0DGJKW40jLweaAdMAgxO+1dLqRse3gD5Lr2SHNPO2I+sW7S1w75CaNoQbh8S2CtnxSz4c5Lm3pQUaKQmIQlOvADaHl/1ftRLgaEjteIqGGU7YPM4XW52DaRVChSEZS8vkDUg66DIcYgD4jCbDQYiBxiPlt+3QRNA/L4zf7Hz3Xff4XQ6A4scgf3ToHiVp8Lq8siMAbevgm2eCqsup2ovq3ApLYcZy6HUAY2T4cxjjTXrLHfrdPG+cmhoheHNYut8nA+8ilZYdUZ/a8wKK5PKHFHJyL/88gsXXXQRrVq1wmKxMGnSpMO9SyZBeAL1yKmDmncZ7LZQLSLw5BZYUwKN4uG5TpqEbIRvF2iCZoIV7jg7cHJm0X5YPgv6XQ8DboBhD2krh4bNtIeVtxoGYPMSrbRyluv9UESOuxw2PgquQkg+Htr8s+LjgfJYAjVJzRqeRcaKDOxxdsr+9DVTNZsNVsF+4CU0JNkC+CdVqvWhQ4cSHx+P0+kkPj4+oPty3mwoWghYoMkwqNs+/F3avQk2LtDl9idA8wgSUspdMGs5FJbocT2whzHnY7fAN/u0vUOdOLi6GdSPYYVVGfA6+vU0Q/P/zAork6o4oiI6Bw4coFevXtx8881cfvnlh3t3TILwPZDpWZ4AnBTj1/toD3y/Xy+0n+0ILQxWeCzZDFMX6vL1/aFNNf2LXE5wFEMdv+mH9B6QnKpRnWU/6WPpx8HmRZC9TqurGodYRrzjVSheCdaG0OEpzc+p8Pph5LFkPZhFxhcZ2LFjG2vTCqvvfc8ZM0YzqMz+O2hT0xfRvl9NgZEE7CMwderUgyLH6XQyderUQ8RO/hwo+EOXU8+D5O5VDBRsl3ZqXg5A6656TIWL2w2/rdJIZWICDOpRvRdUKPyUq+7jVk8ZedMYVli5gXeALWh1/71oz1kTk6owLHTEUz5iieUkrIfzzz+f88832BbYpEbYgBp1CXqldWOMX29RIby6XZdHpkPvBsbGyymE/8zW5YHHqWladSQlQ4M0LRVv1dUXpUltBd36+Rostj0ezrgacraFLnIK5sIej/tu+wxIqiJRNJw8FtcUl4qck23w6KHVVS6X6+B4f2uxsxttvrYfdbZ8UgVqRQAAIABJREFUAAhQiVQ5J6eqDukFf2kSOUCjQdAgAuWfv0crrESgaXudJg33p9ctMGedJthb42BgdzW+NMIf+fCXp4x8WBNoG8P5aUFbiy1BT2B3EdtIscmRT8RC55133mH8+PGsW7cOgC5dunD//fdz661GOglGl7KyMsrKyg7eLygI03bWJCKKgcvwJR+/XP3mhslzwuObNPn4/FS4qqmx8VxueHuGtnno0AyuOC34c+okQ7OOOn3VvJOKHS9N0rXf0J9fwwnnaX+rlkGE08F9KYLNnhmkpldCo7PCfjsV+QEy12bqf/57QEItq2qKAsXFxcYH2QFMBA6gncjvABLR3iWVmD59Oj/88APnnXceAwYMoKioiAEDBuBwOJg8eTIOh4PT0s4mz+O70+BUsB4PRVWMVR2F+7RPmtcrp/XxcCDMtyqi/as27FKBdPIxUNca/r74s6QIZuTqcv9G0FaMjReIadOmERcXh+vss/H0POVadDbx6+nTcbvdnHvuudF/YQNE5Vg0MUxEQsdmszF+/HjuvfdeTj9dC/nmzJnDyJEj2bx5M0888URUdzJSnn76acaNG3e4d+NvhaB5mkuB5sAXQCw9wtwCmZthTzm0S4JH2xqv8Jg8HzbtgbqJcOvg0J1h+w5Xd+Of3oKzblKx462kqt9Eq6/C7Tu07UUo360eK63vC++5h+BA80sA7kMrh6hFVU0GSUhIICUlhfz8fEpKSiIfaBvwKfp5NUP9ENzo9FUVFBUV0b9/f0466ST27fNtdNJJJ1FcXEzupiK2LNb19XuCsyvsCzBWIIoLYP08cJXrsZTWFXLzwn9ra3bCht2asNurnaYahbsv/mwogZ89+9GrPnRwGBuvOkpLS/nll1+guBjOPJPz0Lzwb377jV9++YX+/ftX+PxrCykpKSQkmNlDhxOLeOeewiAtLY1XX32Va665psL6Tz75hHvvvZecnJyo7WAgLBYLX3/9NZdccknAbaqK6KSnp5Ofn0/DhhHUYZoE5XXgHjRPZgaxb6b34W54eTskWuD9Y6GLwRD8qh3wr6kq2EYMgd4dw3u+2w2z39NE0RPO1emF5Mbw20daYXXOncFG8JH/G6y/H7BA17cOLSX3JzMzE6vVWqUwycrK0umoepnwCKpA13BE2sYGe58HDhzgwQcfjPwFpqBzrWVoOPJDAubkhMKezyF7oi43uwZa3BK+EM/NhmlvaA5Y0/YwZITmfYXLT0vh+8W6fPlp0NdgH7O5BZCxSTXgsCZwb+vYlpF/B9z44ovwzDOcMXo0k0aN4sUXX+SZZ55h9OjRjBoVgtvmYcArwE2MUVBQcPBCJtzzd0QRHZfLRZ8+fQ5Z37t3b5xOZyRDxoSkpCSSkgxm2JmEzCI0jQHgeWIvclYc8OXlPJBuXOQUl8H7s1XknHls+CIHtDR30C2Qlq5NO5fP0saK4YocdylsfU6Xm11bvcgBX1sHqBiFOZh/M8YOz3hWPs0RKXIghPdpt5OWFkGtNugX/xEqci5CozoRNNYEnSLaOQFK3tEKwBY3Q6u7whcCOdvgj/9Aolsbww69HxIj2Kcfl8Cs9VCnPlx+KpzTK/wx/PmrAJ7Kg7hGMDQVMttrW5RYMR24DeDpp+ldrx6/Z2TQ+qWXapXfk0ktRiLgnnvukZEjRx6yftSoUXLXXXdFMmTYAPL111+H9Zz8/HwBJD8/P0Z79felSES6iggiMkxE3LF+PafIxctEes8XeXi9iDsKL/j2DJERE0XGfCJS6ojCPu4Xyc3WP5czvOfumCgyv7fIkgtEnMWhPcdutwsgdrv90Pv/Ev1yOotImPtS26j2fRolR0TsIlIe+RBup8jmp/T7m99bZOfbkR2fezaJvHe/yMQRIl89JVJaFNn+TF+ix/WIiSJT5kc2hj9LCkXOXKj/eyPXiZTH+J/9VxGpJ3r4Xib61SQmJgogiYmJsX1xk1qDkfN3yBGdBx544OCyxWLh7bff5scff+S00zRTc+7cuWzbto0bbrgheiqsEkVFRaxfv/7g/U2bNrF48WJSU1Np2zZI4yGTmHIfOhvSCi37jHUN3nNb1a+jRSI8btAUEGD+BvhzvY5z80BIisKUenKEbSfKdsKu93W5zf1gDfEKPmCvq8ds4E1+HkkFo7sjkVB7eh1CLsGbbzYBDAQH3OWwOQNypwMWaDsamkbghLF7A3z3ijpnN+8I598XWSRnxjL43OPhNPQkuLB3+GP4s7oY7luvrsenNoCnO0J8DP/ZF6AtHYqB84CPgafDcP82MYEwcnQGDhwY2oAWS8w6Cs+ePbvK/bjxxht57733gj7fyByfSWA+AG5Axc0MILQjJXJ+2A9jNqnb5Ztd4QSDBhr5xZD5uU5dXXAiXHxyVHYzYjY8AnkzoEEf6PLv8EVcUlLSwRNBWVmZ2lFfCqSiibYGp/hqC4e8z0DsQo1WFqOC7xw0KTvKJ2hXkX53hfPU56jDE9B4SPjj7FwD017XPmgtu8C590BiBOXaM5fDpx6/nfNP0OPayAXB+hK4Y61WOZ5QH17trBVbsWIFOv29z3P7HfBiAPdvc/rq6KdGcnRmzZoV9o5Fm7POOosQdZlJDbEB9bEANQeMtcjJKddoDsD/tTQuckBPBsVlkN7E+BWvUYpXq8jBAukPVnFicrng118hOxtatoR+/cDqO9tU2etqpecEcAtHjcgJuadXNlo11Rh4FvgdGI22dHgEMOi35KVspyaOl26EuLrQ6XloGIItQWU2LYQZ74DbqTk559wVWeLxD4vh6z91+dxexkXOumK4c52KnO714OUYi5x1wBBU5JyM5ohXFjlQMbrnf9/ExJ8jqgWESe2iHO1rWAT0Ax6P8euJwNNboMAFx9ZToWOUJVtgwUZNpLxhQOA+VsHYs9nX0sEIOz0VOqnnQd3OlR786ito3x4GDoRrr9Xb9u11PRUTcsvKfG0dsr72GPFcanz/KpOZmRmwVURWVla1peuREvB9+u+H1/pnBeos9wxq7vQiMA4NPb4bnf05sBxW36QiJ6GpVshFInJW/wY/vakip/2JGskJJHICfe4icN0dWTxuywR0uurSU4xHcvxFzmtdYtvaYTMwCA3E9QR+QPVodX5Pdrv9iPN7MqlBop0wVJsxk5GjS4ZogmCKiGypgdf7fp8mQJ66QGTtAePjlZSJPPKhJml+MSfycYrzRf7zT5FPHhcp2Bv5OEXLPQmsJ4uUbK704JdfilgsInou8/1ZLCIWi9ivvrrKhFz7zZ5E3Xr2mCQhB0oEjmqCcIivF0+8TDt3WsUnjBeRE0Vkj9+6fSJyk4gMEJGdxvZn/3SRBX31e1txjUjZrvDHcLtFFn2vSccTR4jMfl/E5ar+OVV9Dm63yFUjdH2fYXb5bmH4+1KZdcUigxfr/90/VooUGEjSDoVtItJB9HflWBHZHduXMzmCqJFkZBMTf+ahDTtBDWRjnQq+rxye90xZ3drSeCk5wJQFkHsA0hrARYe6JYTM/CngKNEWEMmpkY+T/Y7eNrkA6rTze8Dlgn/+U6VNZUTAYsH13XfYMzMPvdpt6nE37uKKSRJyVVMHscybCHhVX2xjdPxo9i7dq9NSjdH32xNYiba5boqavqQCF6CNOr9BXY/DRNyQ/Zb+AaT0gw5PgjXM49LthjmfworZer/XuXDKpcEjMJU/98cft3HViCy+eCeDPsPsPP2EjSHHh7cvlVldDHevhXyXRnJe7wINYnjG2IFGcjYBndCgm9nawSQqxEB41VrMiE50KBWRbqJXXdfW0GuO3qBXldesiE456879Ine8qdGcZQbCUXm7Rd68Q6/Ed66NfJzSbSLz+2hkoGRTpQdnzTo0klPV36xZhw58kegX9Ubk+xYK3giDt+w32pGcgLwnIokicoqITK7i8RIR6SgiD3nue6NauSLSV0QelrAjXeUFIuvu95WPb31BS8rDxVEq8v1rnkjO7SJLpoc/hvdzj09IPBjJ+XlF+ONUZnGhyIBFvkhOfowjOdtFnQ8QkXYiUjmgaWJi5Pxt5uiYhM06tM9hc+CVGni9X/Jgeq5eoGe0N17OKqIJyG6Bnm2hh4Fw1IIpenWf3iP0/lVVsfcLQKDh6VCnfaUHs7NDG6Sq7TZ5biMwPwwHm812MDE4MTGxZpJCdwCfAa3REONFaPLxDqDUs00icDdq2b0dPYgEbc7pdYj2rguBkvWw+gbI/xUsidBuLKSPAkuY0bLifJjyAmxdCtYEdTvuGUGF1kOjbVjjE3GWO4iLT+TfL9noH0FHdH/mF8I966DIpcn+bxwDDWMYydmJRnLWA+2Bn4F21T3BxCRMTKFjEjY90BzPb1DbkVhS5IJnPFNW1zXXJGSjLNqkrR7irTC8b8XHwkmu3b8D1v+lyydfHPn+uEog5xtdbnZVFRu0DDHruvJ2Amz0LMdY6FRVBRVzWgPXo6LlRdTM6Sx0WqoPMBdf87UOaEn5TrSsvAj11fF2EA9BPO//UZOOy7ZBYkvo+g6kXRT+bu/fAZOehZyt6lR84UjoGEEn86JSuOiGLFxOFTlup4PvPzb2uf+eD/9c5/PJebVzbBOPs1GRsxYVN7MwRY5J9DGFjklENAFOrYHXeXOnNuxskwQjWhkfz+mCL+fp8rm9oGklOwZvi4HKJ2pv3onVr5R74VRAoMNJkGYgKpQ3G1yFkNgaGvatYoN+/aBNm8CJGxYLpKfrdv4Uok5roE6OMSKkKiijfAO8iubagK+qqj/QDTX5241Gbsag3R7vBr4FktEKq/moh87zaNvrjUAIAtVdBlufhU2PaWuOBqdAtw8guVv4b2PzEvjmWSjaBynN4OJHtNt9uOzJhwtvyOKn/2Vw+qV21u8w/rn/sB8eWA9lAv1S4KUYl5BvRb++NWiO3yw0omNiEnViMJVWazFzdI4sNhSLnDJf8wT+yIvOmDOWaV7Og/8N3OYhlBYD+XtE3rxd8ytythrbp/UPaa7H9ter2chbdVW58sq77ssvD32OS0SsookP243tYyBiXnW1TkQuFRGLaK7NrCq2+UZEJopWU3nZJSL9ROQu0f4kIiJzROQRETlfRK6TkCquijeKrLjal4+z/bXI8nHcbpGFUzUXZ+IIkSkvipQUhj+OiMiGXSJ9L9PPt98VdsnO9T0W6ef+yW79P+s9X+SxDbFv67BeRNqKHprtRWRDbF/O5CjArLoyOeoQgRe26YX7gBQ4PQpNKEsd8N1CXb6wd+A2D6G0GFj2k+5jm+7QJD3yfXKXQoHHvbZRdW6Ll10GX3yh1Vfbt/vWt2kDL7+sj1cmDs1D2YmakrQOsct5GN431XmbeB+PGBcahXGhjnEPohGaE9DpKheaX+PNbfFOa7rR952OJpQle9af5vkrBYI4DYvAvsmw7Xn9juIbQ/tMSDkj/LfhdMDs92HjfL3f/SzoOxziIoiWLN4Mb8+AUoeLwVfb+fIdGyl+07nhfu4i8GY2vOVJ7xreFB5Mj22DzjXodNVO4BjgJ/SrMjGJFabQMamV/JwPfxZCogVGRulXcPoyKCyFZinanbw6bDbbQZFTObm2tAjWeMRJr3OM7VPBPD2RJraAekH2icsug4svrtYZ+RC8QsdzIgul+3coeAVTVaIoEsFUJVbUFrc9KlBWoY3UzgXOxlcuXzlvKw7IQ53njkPzdLwnbiGoyHEWwNZnIPdHvd/gFOhgh4RKTdFDEY2j7slk+kTNx7HEwRnXQPf+wd74oYho36ov5upbuPnuTG4dDHWqEOuhJoK7RF3Gv8zR+7e3VOsGo33jqmMZ+tXtBrqjJeQtYvdyJiaAmaNjUgspd8MrnqDFdc01P8coJQ6YuUyXL+4T3AG5uuTatXP1Kr1JOrQKJk6CUOi5yk/pF+IJxmqFs86Ca67R2+pEDuglM8B0vfG6yPrnckTiexNOLpMhOqEiBzTZ2IJGd/Z51vlXSwlQBuwBnkQTju+mYqJxkM+4YC6svMYjcqzQ+h7o8tqhIgeCfwZFOVa+eqpi0nEkIsfpgg9/1eacAvTvBneeU7XICZVSNzyyUUWOBXgkHW5rFVuRMxftWbUbtTeajSlyTGqIGEyl1VrMHJ0jA2++wNmLRQqj5Ob7wyLNzcn4VMQVJP+guhwdt1vks0zNs1gx2/h+rbpZcz9yvjU+VpV8L5oI0UhEin2ro+F7E0ouU1A2isgCEckJsp3Xx+VfItJFRKZWsc27InKJiDQVkZNF5M/Qd8N5QGTL075cnGUXixQtC/68QJ/B7dfaD+bjfPWUSEGw9xeAwhKRFybrsXv7myLTl2i+jxFyy0VuXqX/Y6cvEPlpv7HxQuFHEakneiieLiI18JImRxlGzt+m0DGpVeSViwz0GJV9uSf49qFQVq7JxyMmivyxpvptgyXXPjLSLhNHiLx9t0ipwTYU7nJf+4BDTAKjhUs02xNRcz0/vCInMTEx4uEjFkz7ReRKEUkWkZ4iki4i00TEmyBe+WTuf7+niNwuvrYO3n/nDSIySkRmhvceChaILB3mEzlbnhFxFgd/npfKn8H159sPtnP49SMRZ4Ck92Bk54qM+USP2/veFVkahT4rO0pFLlum/19nLRJZFGFCdCDGjh17yDHwpainI3a7dBw79mBuuIlJOJiGgSZHDe/t0qadHevAsCqmCyLh99VQUAJN6sMplRtlViJY48DdmzXJs8OJkGTQ06dkI0gZxCVDUqx6aMQBt3mWX+CgkV60fG8iNgp8DZ3D+B34FOiLdhP3NtqsbOBnwVdSfg9ai/wDOkV1NVp23hF9j9UldfvhKoKtz8Ha28GxQ/OkurwBbR8Ba93QxoCKn0G8NZF+6TasCXDWzXDmtWoIGC5Lt8DTX8OeAj1uH74Yjjd4jKw8ADevhjmvZFLwbhbvdFVDQH+MNmKtajpvAeDIyoKMDK63Wg/mhpuY1BgxEF61FjOiU7vZ6xDpu0CvNn+NUjm5yy0y5n96VTxrubGx3C6R/z6oV+pbQ5jWCEbuLxpBWHm98bGqZY+INBGN6lwkYh8ThSknDxFFdEpFozL3+63LFZG7ReQ40dJwEY1GVUWe6FScRURaicjHYe+25M4SWXK+L4qzyS7ijDC6MW6cpw1DnH4GVw6wS862yMZyuUWmzNfjdcREkee+EckPI7oUiJn7ff9bPe6JriVA5eBb5XHGee5n1lRbEJOjEnPqKkRMoVO7eWGr/hDftMp4HoKXZVv0hPHP/4iURDiF4GXXBhU5794X+XSEP3sn6Ul27X3GxwrKdBFJErHjOQk9bPwkF3GOzk7RHJpXK62fLSK9RftPiVTdg+obEWkg2hDpo5B39SBlu0TWj6qYi5NvoHP9mEf1PQ/ro9NVN18SuWgsLhN5Y5pP5Hz0q0i5wRw1t1vkg10ifTweOfeu1by3aORXTRLVnFVx2HqfmRy1mEInREyhU3vZU6aJkb3ni8yJ4tfzr+/0pPHZH8bHmveVCp3pbxofS0Qk+z+eaEJmdMYLym8iY+uOFTt2kQ4isrjiw3a7XcaOHRvSUIaNAk8T7Qhb6reuSNTQ71QR2e1Zt1JElvpt87uIPBvSLlbAXS6y62ORhf08IucUNf9zlYQ/lpcH7vKJnHfuFVnjEUyRiIbsXJGxn+qxetdbIr+tiny/vJS7RZ7a7DMCfHpLRSPASMVIqYicICJnSWChIxKdHDATEy+mYaDJEc8Hu8EhmjNwaoPojJlTACu2aXrHWccZH2/rcr1t38v4WABOTzsDa8Pqt4saZ0Dm4kw4H21/cCJwE3An0Dt0/xWIglHgSOBGNC+np2ddMtAL+A71/RHgFlhRsIJJV0zi8XGPay6PX5uMUDx7CubCtpeg1NP3K7kntHsM6gbJ1wqEoxTmfg5r57oY1sfOjZfYGHwrpDTXx6v6DKrz3LljZBaLNro4YWgmjZLhjrOhQ7PI9s1LvhMe3aheVBbg/jZwbbOK5ePVeUUFYjkwGPVsnAQESmWqKgesRhq9mphURQyEV63FjOjUTgrKRfot1KvO36KUmyPiy3UYH4XS7bJin33/gSjt444JGl3Y/FR0xguZPSJyhWg81/vXSrSSaaqIRBLlcIrIVtHKqSdEZJiI/Lea7V0i0l20FYNfCwNZKJp7s95zf4zIpEsnRRQ9Kt0msu4B3zTVokEie77UXKtI2blW5OPH9DiYeLvIH5+FNo1Z1b46nCJX3Krr+wyzywuTRfIMVvKJiGwsFrnEU1l15kKRWblVbxduROdHEYkTkVv81v0lIoukYtAtKrYDJiaVMKeuQsQUOrWT97L1R3n48ujl5rjdIo97SnPnrDU+3rYVenL7+FHjY3nZ9ZGegDdEccywmCMqeOpLRdFTR0SOFZEhInKTiDwuIq+IyEsi8pSIZIjm0dwhIueKyDHiqR+u9HdDkNefLiJ1ReQZ8fnoPCMig+UQX51wTp7lBSLbXhFZcJpvmmrrCyLlBv7tnQ6ROV/4xO5Hj4rsWB3eGP77nFMgMuRqn8j5ap6I04AA8/Jrnkh/z0XDhUtF1gYQTpGIkf+JHhrvimrT60VTqpqIWheNE5HMrKzY9j4z+dtiTl2ZHLE43PDJHl2+rnn0nFk37oG9BZAUDye2Nz7ebs+0RySdpgMR75mychVFb8ywOA34HHUTngVM9vztAFZ7/sIhHugA9EFbNwwIsv0Q4FG0rcOHQDPUPvcloEnFTUPpP+Z2wN7PIftdcHmmBRueBm0egLodw3wvfuzZDD+/D7k79f4xfbVXVWIYJeiV38M4+xO4nA5Ov9TOhPE2eraLfP9AW0R8sBte3aEzfifWh+c6QuMqSturcsL23zf/+6Ctw+KAq9CpqwzgKaAd8C/P46uB24GhvXpht9sZE4veZyYmEWIRkcqOFUctBQUFpKSkkJ+fT8OGNZUYYVId3+6DzM3QLAG+6QEJUXJ2+myO9gY6tTPcMsj4eD9OgM2L4PThcPxg4+MB5M+B9fdCUjr0+Do6YxpG0EaY21DBs91zuxdIQPtEef+S0T5UHVGB0wZf/6lQcaJnyV9QX527UcETgKSkpIO5H2VlZbrLLtj/PeycAI5dul2dDtD6Pkg5M3Lx7HTA/CmwbLoKiToNoP/10P6EyMYrd8JXf8K1A5JwOx1Y4xPZtb+MNIM5aSUuyNoCP+bq/UvT4OH0wP9LoTZ2zQPK0a+opd82/wBygbcqrR+NWiKtJHDujolJpBg5f5sRHZPDymeeaM4VTaMnckTUcA3ghA7RGbNwr96mGEwS9af+8UAclG0Dxx5IjOLYEWNB+2MdE2zDKBEP9PD8BaFygqvdnsW9fW3sfNOXaJzQFFqOgLSLwGLg1y17Lfz8ARR4js9OJ0Pfq6BuhKJkV552Hf/mv1m4nQ7iExJxljv498vGknS3l8GDG2B9iWrMUelwZdPqxV11idvefZkM2FGhk4P6MF7j2eYd4E8qihyAxmjk5wCm0DGpXZjOyCaHjeUHYGUxJFj0KjRaZOfptFV8HHRvY3w8ESjwdHhu2DT49pmZmQGdhv2dZ631oV5XXV+00Ph+Hs34T7eUlpTx6I12xo7N4LEbsyjdCNYG2oCzx9fQ9NLIRY6jBH79CKa8qCKnXiM49y4YfGtkIkcEfl8DT36lImf+5AzufsBOuaPskOaq4fJbPvxjlYqcJvEw4RgY3sz49O8rwHXAzUAmGsG5Dm0GD5AInFnF87LR2cpUYy9vYhJ1zIiOyWHDG805p3HVuQSRsswTzena2liHZy9lB6Dc0zqhfpPqtwWfDT5UzHXwP1l7adAbildBwTxIPc/4vh6NeD+3cZl27j7Nxqrr4PJ1NnJbwYSdGdTvDU9/ZCPewBSQCGycD3M+h2JPfs+xZ8JpV4Sfi+OlqBQ+/g0WbISFU1XkPDbGzpNZwfNiqsMl8Ha2/glwfDI82xGaJUa2n/78AXwA/Ae4wrNuGDAd+Ax4uNL2ZUABGvH5yLONefVsUtswhY7JYeGAC2Z4cgquCCFKEg5rs/X2uChEc0CFDkBCHYgPQThVdQKrKgEUIKU/7P4Q9v+oEYmEEITU4STU/I5o4ixz8chwO5cutLFpiq6LS9bPtcV2kHiXIZGTvxt++wR2rNL7Kc2g3/XQqmvkY67cDu/NhvxiiLPAMS1cXDTOTkaGsSTdPCfYNsGcAr1/eRo8WE0+Trg0A5qiueRerGgUp7KOKkG9dMajOTsz8FkimZjUJkyhY3JYmJkLZQLtkqBHFLv8ud2wYbcud24RnTFdTr0NpzljKFVCAPVPhOQecGA57PkftL47OvscK8KJVoVDVQLKmQd7PoOcT6yUl7hwtAJrCjQbDs2u0aq1sUSe3+J0wJJpsPgH/Y6t8XDC+dDr3NAErUjgaaLfV6vIaZ6iyfDtb8sMOE6okZwlRfDYRthdDkkWeKwdDI2CMBY0NQugM/Ax0Mhz34EKnDigYaXtk9DUqn+gFVdRCCiZmMQEU+iYHBam7tfboU2iV1IOsGM/lDh0yqpNlKIjB4VOmP8toTjPWizQ4ibY8CDs+QSaXgGJzY3vc6wIJ1oVDv4C6oErbOz9DPZPg7e2ZDFhZwZ3d7WT/jCkDYO4Osbeg4hW0M35Aor26bo23eGMa4Inm2/LgR25cFqX6o/ba8+EtIYw9CRINPgrKwIf7YFXt2sD97ZJ8Fwn6ByFjN830SqpUjTZ+HRU5JSjRXZe8bIf8AbNLMAaIB043vNnYlKbMadTTWqcPQ6YX6jL50c5c9EbzenYHKxROrrdnlmFuDDHq8oGvypSBkByL3CXwsbH9La24p2W8ibSJiUlVYjkRDpt9fjDNkZfp2M+OCSLfVN8IufRG+y8utxGs+HGRc7+HTB1PEyfqCInuTEMGQHn31e9yDlQCv/+UZOKS8pUfFRHch249BTjIifPCQ9sgJc9IufsxvBBt+iInNuAx9AIzS/AfcCtnscSUH8cgCLPX3fP/ceB3qjrgInJkYApdEykwQFPAAAgAElEQVRqnF88yZ49k6FlUnTH3unJ+2kbxSqupHp6W3og9Of4RznKyqqvsLFYoN3jWoV1YAlsfBTEGaWdjzL+kRevgEtM1Ov+jIwMrNbwjHSKV8PW52Dp+XDFaht3tLIzYWcGfRcnMWFnBuPG2XnqfZuhUnGAkgL47WP4Mgt2rtHo3IkXwPBx0LF39dGZbxfAQx/qtOi44TCwR3SjkIFYXATXrYRf87Uy8ZF0eKoDJIfrVVQFK9AS8W9R078VqNCZDYzwbBOHiqvtaDSnAXAlGgX6HuhifDdMTGoEc+rKpMaZnae3AxpVv10kZHuETsvG0RuzXorelpdCeRkkBBFn4TrPgjr3dhoP6+6B/F9hy5PQLqNmTqjhUPl9eMVOOFNX5bmQOw1yJkPJWt/6hOaQcaeNd+/yTfdVTt4NF6cDlv0Ei6f5Kuc6nASnXQ4NQhDDm/eq0BnQHa7x1FTv3A9JCdCgrkZsqsvViQSXwLvZ8Fa2RlXaJsFTHeHYetF7jc2ouPHm61vQqSsrcAeajHyb536RZ/veQHNgPZASvV0xMYk5ptAxqVEKnfCXp2LkrBgKnVZRFDoJdSA+UU+axfnB8zgi7ezd4ETo+BRseAj2TdEIT5uRYDkK4q7OAsibBbnToeAvwAUTd2oC8kPX20i7GBqcDE88pdN9VqvVUNdrtxvWzYH5k+GAR1intdVy8WDVVDv3a2TwpA7Qvimc3Bn2F2mZ+O9roKAYih3QOBkuOxU6RTGnardDq6oWetqCXJAKj7SNThTHnzboVNRcfGXk9YCL0Zydl9Au5R3R6qpS4P+AV6O7GyYmNYIpdExqlL8KNRzeLgnaGcy3qEyJAwo9V+3NonjJabFA/VTI26WlyMGETijOs4FoNECnsbZkaXLygZWQPgqSu1f7tBrDP1rln2g9ZsyYQ6JV5TmQ/xvkzoLCeRWn4+p1g9SOVp79PINW9cB2mi+pedCgQcycOZNBgwaF7TEjApsXq8Dx9qaq3wROuQQ69aleNBaUwEe/wpLNcHxb6OPpazasDzz1NWzNgZ7tNMF4T766b0/4ER65hINtHIyU38/IhSe3QIEL6sbB6LbRqaqqio5oB48v0ZZn3shOY+AidGpqi2e7fqggOiU2u2JiEnNMoWNSo/zlSUI+NQatxgpL9DYpPjpGgf4066hCJ3sttI1xmUnaxWBJgK1Pac7O6hsgdaiWnvu3iTgcnjbeaBVQIdEaYNw4O6W7Xex8S6ffildWfG7dLtB4CDQ+B+qkwzPYSM7SabDZs2cfFDczZ848GBHzih+oXuyIwPaVMP8b2OsxjEysByeeD8cNDF4u/sVc+GkZ9GoHg3rA8m3gcmtCe9OGcHEfKHPC4B4Qb9UoTu+O8PgnMG+dih+RyMrvD7jghW0wxVMB1r0ePNEB2kbxQsC/hNyN5ts8BwxEIzt343M0PgXYCWzwPO5dZ2JypGIKHZMa5U/PtFUfg40M4dATvTeaU79u9E/0rY6BtX9A9rqoDBeUJhdAg5Ngx+vasHL/VMibAc1vgObXgTU5dp421eFtb5GRoYnCD99gI/MRvX9XRzu3NM4ke65v+3rHQaN+0Ggw1K2i75h/zo/Vaq0gcvwfr85QL3st/PUN7Fqv9+OT4PhB0PNsSAri0bS3AGz/g/Q0GHUhdGkJv6yEVTv0sRae6dX+3cFRriLHS2I8dGgGmzwO3xZL+OX3i4tg7CbY4VAhcmMLuL1ldAwAPV1LaIivTNzbidyNRmpeBMZ61o1AzQLXoH2sjjW+CyYmtQP5G5Gfny+A5OfnH+5d+VuSXSbSe77IyfNFCsqNj2e32wUQu90uIiKLN4mMmCgy5OqK66NBQY7IxBEib94hUlYctWFDomiZyKqbReb31r+FZ4ise0Bk7ySRsaMrvtfKn0k0Kc8Veew2Hf/+M+yy5FzfPt3RStff2dYu6x/UfXPsDX3sxMREASQxMTGk7d1ukW0rRSY/r9/LxBEib98l8senIsVh/nsv3VLx/vpskdsniuQU+F6rKvYVimR8KjJz+aHbeb8H7/uq/H04XCKvbdf/hd7zRS5cKrKgILz9ro53RaSDiBwnIj1F5E0RKfI85hQR/7eUISInikgLEblBRJqIyLXR2xUTk6hg5PxtEQnmCHH0YKTNu4lxftwPj22CbvXUCyQa+F8tn3eNjbtGaU8hI+Z1gfgsE/Kyoe/V0GNg0M2jigjk/gQ7J0DZFr8HLPC+K4tXF2WQGJ+Iw+kgc4ydsVmRvXdxQvk+KN2qr1O6RZdLN4Ij25NAbLFya0sd35II9Xtpz643FmUR18TFOHtmWK/p/Q59XckDf3cisG05LJwKezbpurh46NpXy8XrB0lCLyqFnEJNJE4JUMW0rxCenwwX9oYzqwhrFJfp35fzNK/nxgFV54QlJSUdnN4rKys7uH59CWRsgrWeqdYLUuHhtlA/SgnHHwN3Ac+iJeDvA8uBXsC7nm1c6HSWN6S/GO1Y7gLaoonHJia1CSPnb3PqyqTGWFmst9Fs+eA/VZD1xBOUOxycc230RQ7AcQPg9//Bilm6XJPVUBYLpJ6tOS4layHvF8j/RRuC3hhnY4LlCRxOBwmWRC783sayRRDfBOIbaGdva339E7eKGSnXW3cZOHNV3Dj3a9sFqrn0ue+0TOp1g+RuUK87JB8HcZ5y+6wI2jFUntYJlJPjdsHGBbDkR9i3TddZE6BbP+h5TnCBAypMfl6hoiTvAAw7GU7trKXibtGeVKDTUwlWrdzyp9wFq3fAr6tg3S7o3BxuHwINqxBMVZlFPjrGxvu7tGzcKZBihUfbwZAoVggC/AxchrZlAM2z+TfalXw08AxaNg5QjFZbneD5MzE5GjGFjkmNscJjuNc9in4gULHVQlx8IkOuir7IAehyOvw5SSuvtq+E9B4xeZlqsVigXlf9a3UbOPbA2PuzKF/gICEukXK3g7ezs7gVG45dEb6IFZJaQ512kNRWb+u009e01o/eewnFb2j0QzZW/w7LZvjaNcQnQvcBmoNTL8TqugUbYdlWuP1sFTqzVqjo2ZULw/v6thPRSE+Dupp707+7TwQlWKFZQy05P/9Ezc8J5X1573+5F+Jv1PfXLwUebwdpUU6aB/W9Kfa7bwGuBfKA/6JVVpcAC4APgZsxm3GaHN2YQsekRnALrPb8+naPYkQHfFfPCYmJlDsc/Pi/LB4aFn2xk1gHup4By2fAwu+0P9Lh9rh5dmIWz3xy6Ek17WIYdbUNVyG4CsFZCO4DQJxWdFniIS5Bl+MbQ0IqJKRBfCrEp4Alyr4tVVGd35CjBLYud/Hxo1DmOW7qNtAKqu4DoE6YgmveOmhSH45L1/uXn6pJx1/OU5+cDs00ghMXp9VWbdPUS8fhrNjGoXkjuOCkwK9TWeQ4BdJvs9FmJyx5NYP2Tngty8YFqbEzg0wHpqNVU54KeRoBl6Nl4lNRobMWnebqhSl0TI5uTKFjUiPsLYcSt4bMo+mf4++9MvGjGVxzexY/fZJBVjdfxUs0q696DoHVv8LuDbD6d506OVwEi4gktgzdf+ZwUPk7EdGqthWzoO1+G+mNVOSkNNPoTZfTNJoTLi431Euq2J/KGgcndtAy8k//gNGXqMgR0ceSEjQP58knxmG1xpERYgm/v3hbUwxZW1TgN7/VRpskOLOBK2beOF4eBiagrR1ewFdx1RVt2vmm5/41QHvPOhOToxlT6JjUCNs8uZitkiDewJWsf0l5ZYO5W64ajLthP0652F7BmyWaZdb1U6HPxTD3c5j7BbQ+Fho2jdrwYRGpA3Ntw1EC6/+ClT/D/u0wZX4mcXFWbrvGRo9B0Lanr6FqJMLVGqdRmexcNfrzJg43SoZ+3eCT32DVdujWRiOPVgu0aQLTFkPL+HqMtT2MhdBK+DMzMylzw+s74L+7NLm3oRUeSIeh42010tIjFXgHuAo4Bp2a8gZR09ES8lzUHNAUOSZ/B0yhY1IjbPN43KQbbOLp7x3jf6IfPHgwM2fOpFVXuPCBGdTJm32IJ0u06DEQNi3UqM6Mt+CiByOLNBjFiAPz4UYE9m6GVb/Chvng9AhhawI062Dl3S8z6DMMLjwhNH+g9btUuKQm+0SRiOZVx1ngrO7aeXxtNjRp4Ots3zpV7+/KV6HjXZ/eBDKHQ4tGDxEnpSH74iwq1CjOVs/7GdIYHkyPTS5OdVwOPAGMAvai7RzSgTeAzqhhoInJ3wVT6JjUCDvVPJfWBoWO/9SMf16Kv6vu23cn4XZWX6ZshDgrDL4Vrh+aiSy0klDHxrl3aU8sL7FyJj7SKSmEDX/Bmt9h33bf+pTmOg14TF/4v2Qb7bNCM91bvwv++7Pm0gC08yQKt2+qOTAWNErTKlWdjGcsgzap0N6TSNykgUZ5vJVm3sTjFn592PyPOW/Se+X9yHfCqztgkselLy1BWzjEop9bqIz23E5Fp7DaoGLnQ9Qg0MTk74Lpo2NSIzy5Bb7OUdfX21oF3z5Ye4OZM2cye/bsQ7xXEhKTcJY7SEhIxOEoq2Lk6PHw/Vk8/68MhvWxc8sVNs67F5LqVZ0783fGWQ5bl8LaueqBI56ybWs8dOgN3c6EFl0OTc4N5q9TUAJv/aSJw4OPh9158NWf2v7j6jM0WuN2AxYVL//f3p3HR12dexz/JEM2AiQkEHYIu1CKsigFlUWvKK0oWqleW7QVa7GIIretRZsQEhEt9treqshSwXqrqG2tXreCiqClKCARAUEBMUBYEpaEBLLNnPvHmclGCAmZ5Dcz+b5fr/NK5pfMzGEImYdzzvM8p0tg/qv2+vgL7QHkL7Ph5X/DDy8/exaVT011cYyBt4/BE/vhuDfYmtQO7usCrQPkv5E52L5VZdiMK5FgpDo6EvBOeN8E4ur4E1eX9gbr1q0rf+Px/a+/rNSmmJeWnn/n67r67e9TCA+Hx57wbqWVpfBZUQbpDyvIcZfBgR229s3ezfYcjk/7HvZgcZ8REF1LBl7lsgG+v+PKcvJhz2GYPBISWtkxORzezoTXNsDPr7bbWB6PXamJibTfu2Y7PP5/0K+T3cq6tL8NlmpTU12c23+ZwqNZ8Im3f1uvaHiwB1zkxxT8c1kHPI/dkjrb8Z/23iHSXCnQkSaRV89A51w9g6BqU0nfGZ2Zv0rnVO8Utv6zbs0gG+rR/04hLBwe/V0qb336MGWeEn5xb92CHCeacjYmdxlk77Rnbr7JrEgLB4htC31H2ACnbae6PV5NwUXl16qgyKZ7V36D79cZjuTD25shcy9clFxxZgdsV/K+nWzhv9x8uOGSugU5lVfo5qTb20uyIenOFKLC4M5O8KMO/ulRVRelQAYwD9u3ajiqZixyVn5uRxHQ1OvKOVO2254+a4/X73419Qyq3s/piiuuMIC54oorTGmZMTOetT2v7vtl4/V9qs43vxbhkWbRXca8u8SYgnP8Wc/Wl6ox+1X52+mTxnz5b2NWLTLm2Xsr+k4tusuYP//CmA//YsyBHcZ43PV73OqvQU2vSW6+MdOXGvPxV1XveyTPmEWrjFnyrr1dUGTMe58b8/Xh+v/5qj/v2uPGXLfFmM7T7PUR96WbfUX1f9yG2GGMGW6MwTt+ZIw50bRTEGlyDXn/1oqONKn6ptdW374Azjj/cvnltpjN+++/z/xHMhg+OoV/7YQh300hvXXjp1lXX3l489MMvkcKWVtg2EQYdIU9wFzTn8335/HdDvTzPR4P5H4DB76Afdts5lnlU34xraHnUOg1zJ67CT+PFY66VExOSUkhsTVc2MMeMB7UzdbKAWjfxh4y/ibH1tCZ97tnOZ5wM1cNiS0/hFz5uWpbOfNl9t3xqxTu3wUf5tnrF92dwvi20D3CTdcGHrCvKw/wFPAAcBqbHr4Qm0YuImenQEcCWvUgoqaUcd+blO9Na1R/+NdOW/b/tw+kEN2Iqd9nK/cfGw9je6Ww/q+w7QP49pXQbyRExlS9f10yepxkPHD8EBz6CvZ/YbemSk5V/Z6ELrbWTY/BkJTc8GrR9akPdN3FkP5XWPcljBkAEd7faHEx8NW+AubN+wMxeFix5Md0mjSUG0fMBqoeaK+tztIDKWk8dxh+sA2KjS14+aMOMLUTtHzkzL+jxtqOzMLWw3nfe/s/gOVAl3o/kkgz1AgrTI3qqaeeMsnJySYqKsoMHTrUrF27ts731daVc3xbV2vqsXVVl+2Lmng8xqSssNtX729tyKzrN7/q12dMTTfL76/Yylk63Zj3lhqzb5sx7mpbOb6tr8jIyMabcB2UnDYme6cxn75lzNt/NGb5zKrbUYvuMmbZfcb882ljtq42Jj/H0ekaY4x5Y5Mx//VnY9ZsM+Z0sTHFpcY8s9KYWf/9b0NY2Bnbnb7Pa/tZ8niMefuoMd/9zP7cDttozLSdxuw5Vftc/L0d6THGPGuMaWPsNlWMMeZJY0w9dwJFgl6z2bp66aWXmDlzJk8//TSXXnopixYtYsKECWzfvp3u3bs7PT2pRax366awjrtIdd2+qElYGIwbBCv+BW9ugu/0tRk3/laXlYdbZ8NX62Hr+3DiEOz6xI7IltCpr62s/MLbtR+6bQzGA4Un7JxysyB3n+0KnneEM7qXt4iEpJ7Qub/t79Wue81bcY0hK9d2C//BKNtUsybfGwonT8Mbn8JHO+F0sc2ymvGj7xBfMLf85yg9Pb38Zwc468rZjlOwIAs+8zah7RwJM7vCuPhzb736czsyG9uB/A3v7ZHAc0DfOj+CiECQ1dEZMWIEQ4cOZeHCheXXBgwYwKRJk5g/f/457686Os6ZvQdWHYdZXeHWDuf+/oZuAbg9MPcVOJwHEy6CSZc0YPJ+4KsE/OW/bcE8X0bSm59m8PrGVG4clc7UH6Tw6r8yWPpKKr+YkU7a3BRi2pxfUGEMFBdC4XF4+JE03KUupkxMIe+wDWbyDkNZiX1+j8fNxOFp5fdtGQ8dekHH3tCxDyR2a7rABuw5oC1Z8P5W2Jltr902Gi694Oz3KXXbTuRZuUCYTRn3qV6PB6hSC8fnUAksPABvHbOxXnQ4/KSj3aqKqud23LlqANXGYAOa+7EdxyOxGVb/hd06E2mOGvL+HTSBTklJCS1btuSVV17hhhtuKL9+3333kZmZyZo1a875GAp0nPPbLHg5x75xTG+igwWZe2HhSrsSkH6zrbMSCDxuu4KSlprBwudTmTQinQkXVrwJ+oKf64an872hKUS2tIUIo2PtSpCrhnXY0mJbq6a0CEqK7Oeespofr/x5Nmfw+oZUbr8unZl3p5DYDdp1gxiH/mkUFsO6nbB6Gxz11qYJD7MVjScMsYX+zpev2B9wRvBR4IbnDsELh+05HICr28K9XaFDA1YCayoweC77sKs4b3tvDweWAYPOfxoiIaFZFAzMzc3F7XbToUPV5YAOHTpw6NChGu9TXFxc5RdMfn5+o85Rzi7B2+snt7TpnvPCHrZmylcHbZuAeyecXxaQv4W77KHdpF526+vB2SnkfgPHD9qtpO7fTiEqFk4X2H2+klN2nMyt/3NFt4bbrk8hrgM8/2YqXQfCr2alsPgFG+QEwsHnfUdh7Xb4+Cso9gZnsVFw2QUwZqDNoGoI34F2n9/85jeA3V7aVgjH/jOFY97nHdIK7usKg2opZFif56zrdqQBlmJXbU4CUcBc7+2g+SUtEqCC7t9QWLVNcmPMGdd85s+fz9y5c5tiWnIOyd4+ULtP1/59/hQWBj+8DB55Fb44YCvmfm9o0z3/uVTeeuvQ2w6fa+5JweO220/Fp6Co0Pt5YUULhcpLsRGRNqMrIgYio+3nMa1tk0yA75NCX2//qKUrnM/uKimzWXFrtsPXRyqud25r2zlc0sd2HG8o3xYSUJ5dlZqaym2z07lgejovPZZK5zwYMT2Fe7vCmLj6l0A423NWz8SDms+VfQXcBXzgvf0d4FlgQMOmISI+/j0X3XiKi4uNy+Uyf//736tcv/fee83o0aNrvE9RUZHJy8srH/v27VPWlUO+Pm0zV0ZtMqbM07TPvW6nzcD62SJjMr9u2ucONE5nd32TY8yLHxkzc7n9O7lrkTHTFtsCfzsO2Gwnf/FlOo0dO9akp6cbj8eYf+cZM3iGvd55Wrrp9fN08/1fzDElfkpjqk/WVYkx5hFjTJSpyKj6nTGmzD9TEQkpzSLrKjIykmHDhrFq1aoqZ3RWrVrF9ddfX+N9oqKiiIpqompeUqtuURAVZs9AHCiG7tHnvo+/jOxneyKt/QKWvg/3fRf6dGy65w8U9d1O8ZeCIvhkl61ttP9oxfXEVnD5ABjVH+Ja+v95K2fFbSmAu7+CjSch4vYUepTB4Bg3L/42rTwj0N/PWVn1GkAbgDuBLd6vXwU8A/Ty31RExCtoDiODTS+fMmUKzzzzDCNHjmTx4sUsWbKEbdu20aNHj3PeX4eRnfXjHbC1EOb0gInn6C/kb24PPPVP2LbPHk6+80rbB6m5ONt2SmNtX5W6YWuWDXC2fANl3u22FuFwYTKM6gcDuzb+mamdp2BhNnzkrWgcEQaT28MdnSDegf/m5QEp2ArHHiAB+D3wI87elFNEmslhZICbb76Zo0ePkp6ezsGDBxk0aBBvvfVWnYIccd532thAZ11+0wc6rnCYdhUsedemLj+zCm69FEYPbNp5OKEhNYnqw2Ng10H4eBd8ugdOVZz/pVuiTfm+pA/ENsFq3t4iWJRtSxqATcu+NhF+2hk6NmKl7LMxwF+B+4CD3mu3Ak8ASWe7k4j4RVCt6DSUVnSc9VkBTN0JrV2w6kJo4cB/Yd0eeOEj+GiHvf3dITBxuE1jDlWN2SXd44GvDtnAZvNeyKvUHiI+Fi7uDSP6QLcmCmx3n4Y/HbQBju8X29Vt4a7O0KMJt0urzAm4B3jHe7sP8DR2u0pE6qZZ1NHxBwU6znIbuOozyHfD4n4wtIFpw+fLGPi/TfDmp/Z2n442O6tzA+q0NCdlbpuyv3kvbP4a8itl0sVEwtCeMKIv9O3YdOn8u07D0oPw7vGKa2PiYFpn6NsI53/qohhYAMwDirCF/2YDvwYcirlEglaz2bqS4OYKg9Hx8MZReP2oc4FOWBhcN9zWZ1nxL9h1CDL+BlcNtunnURHOzCuQFRbb802f7YWt+6CoUj2k2Ch73mlIT7igy9lbNTSGLwph2SF4/0TFtSvi4c5O0M+hAAdgFTAdmzoOcCV2FaefYzMSab60oiNNaksB3LHTZmC9NRjiHA61jxXAS+tsFWWwmUCTR9o37obWUwlmHg98kwvb99ux57A9g+PTOgYGd4ehvWBAF3sGqqkYY7Onlh2CT05WXL8y3p7B6RNz9vs2tmxgFvCS93ZH7Dmcm9FhY5GG0IqOVNGYZzIa6tux0C8GvjxtV3Z+WIe+V40poRXcPd6uVKxYB0cL7EHlrolwzYX2jbwp38SdYgwcOgFfHrT9pb44AKeqdS3o3BYG97AVp5OTmv5ck9vABydsu4bt3rNALuDqBLito7MBTjHwB2xPqgIgHJiBrW4c59y0RAQFOiHJ5XLVmE1TOfvGKWFhcGN7eDQLXjwCN7Wvf8PExnBhst12eXszvL/N1ntZ+j7ErYfLL7ANJQOlV5Y/uD2Qfdyu1OzMtgHOyWpVq2Mi4YLOMLAbDOwC7RxaBC3y2KD4L4dhnzf4igqD69vZhpudHS6V9Q42m+pL7+3vAAuBixybkYhUpkAnBNWUOtzYdVPq49pEePag7Rb94hH4cYAU74uKsF3Orxpsg521220W0Ruf2pHc3gZEQ5KhY3zwbG0ZY7fo9h217Rb2HIZvcir6SvlEuKBXB+jXyW5HJSc1/WpW5dXIvDJ4JQdeOgLHy+Dg0gwijZtfpaRxc1JF/zSnfI3tMP6a93YH4DFgCnZFR0QCgwKdEFU52Hn4Yed7G1UWHQ4/7wJpe2HZQbg+EdoG0AHg2GiYOAwmXGSzitZ8YevD7M2x47UN0CEOvtUNeiVBzyR7sDkQAp+CIjh8wm5D7T9mV6b2H61a08YnOsIGM/062ZGc1LQHiWviW41ccwJO/zClvJv4yWUZZD+TSkpaOnd3cXaO+cB87NmbYuwv0XuBVLRNJRKIdBg5xEVFRZWX/a/cyd1pHgO37YAdp+D77WB2gNd8zDtlK/xm7oUdByoq/fq0joGe7W2Kevs2FSM+1r9nWUrK4HihXaHxjaMn4UieDW4Kz/JXHB5mz9gkewOzXkl2VSoQurmDPX/zYR68eBje+L0NajpPS2fMjBRa/G8Gz893fjXSjW22+RvA14f0SuB/gGZQd1LEUaqjU0fNLdDxbVf5ehs5/UZR3aaT8DPvwYYn+9rKycGgqAS27bdp6XsO2y0ht6fm723hgjYxNgU7Ntr7McqunISH262h8DA73B4ocds6NSVlUFpmV2IKiuwoLDpzu6kmbWOhQzx0SbAVibsmQqd4O5dAc6IMXsuFv+fAAe+qkwuI/N8MPvp94PzsfgDMBD7z3u4HPA5ci7KpRJqCAp06ak6BTlP3Njpfj2bBX3MgsQW8OND5cxfno7QMsnLtttbhPMjJh9yTdqXlbAFQQ0S2sGnwbVvZA9IJrSCpjV2hSYoL/DpAxsC2U/DKEVvBuMT7G6iNC25oB5OTbJuGQFiN3IYt8PeG93Y8MAdbIyfAX2aRkKL0cqmiqXob+cPMrvDpSdhTBLP3wFP9nGkN0RARLaB3Rzsq83jsNlP+abuldKrYuzJTbFdtPB67ZePx2IDIFW4fK9JlP0a4oGUUtIquGLFRNhsqEM4D1VdeGbx9DF7PteUFfAa0tI02xyfY81vgXKd1nwPYgGYZtvmmC7gLSAeauE2biDSQAp0Q5Ha7a1y58d12u91OTKtG0eHwWC+4fQdsKoDHsmB299DoPRUebg8pJzpUAToQGAOZBfD3XHiv0upNZBhc1dau3nyrZdXA7WyrkdD4AXo+tm3D7wBfLHYj8JPQJFEAABehSURBVAjQv1GfWUQai7auJCCsPg6/2mMbMd7UHn7VLTSCnebqUAm8ddTWv8mqtOvUN8bWv5mQUHNV7LNtsTb21msxtvbNPCDXe+1S4LfAKL8/m4jUl7auJOiNawtzkmHuXntmBxTsBJvTblu5+P+OwoaTFd3Do8NtB/Eb28PAlrVvuzX1aqQbeAFIAb7xXuuPXcG5AR00FgkFWtGRgPLGURvsGGza+QMhso0Vqko9sC4f/nkM1ubZKsY+Q1vBxES4oi3EBljGl8EW+ksFPvde6wykAT9B/wMUCTRa0ZGQcW2i/Th3L/wtF/YXw9ye0E4pLgHD7T13s/IYvHsc8iotsnSJhO8l2tHF4dYMZ/Mu8CCwwXs7DpiN7U3lYMNzEWkkCnQk4FybaDOvMvbCxyfh1u0wNxlGquysY9wGNhfYA8XvH4ejler5tIuA8W1tc81zbU05aT3wEPC+93YstkfVL4C2Tk1KRBqdAh0JSNck2C7nD34Nu07DjF0wpQP8vDNEBEg131BX5IFP8m3F4rUnqgY3rV0wJt7+PV3cGlwBGtwArMN2EV/pvR0JTMOu6nRwalIi0mQU6ISIys0Qq8vIyMDtdpOWltb0E2uAXjGw/AL4w37b3PH5w7aa8gPd4VuxTs8uNB0rhY/y7KHij/Mp7zUFtqDfmHj4j7ZwSevADzjXYWvhvOu97QJu814L8I4jIuJHCnRChK8ZInDWtNxgFB1uA5tLWkP6N7D9lK25c2W8bQzaI9rpGQY3j4EvTsG/8myA88WpimwpsBWKR8fB6HgY3jo4ijl+CDxMxQpOC+DH2HM4vRyak4g4R4FOiKip8nGgtn04H+PawsBYWJgNbx6F907YVYfr2sFdnaB9pNMzDB7ZxfDJSdiQb9PAj1Xrn3VBSxgTZ1dv+sYE7pmbygzwHpABrPVea4HNoHoQSHZmWiISAJReHmICvZGnP+w6DU8dsGdHwFbZHZ8AN7eHAdrSqmLOnDRO4WLsvSlkFtgzN77mmQeXZmDcbvrencaINnBpHIxqE1xBowfbh+oR4GPvtUhsgPMA0NOheYmIf6mpZx01h0AHAqMZYlPYfBKePACfFVZcGxQL1yXawKdVgNVuaQpuA18XwZYC+LQAXn48g6+eTqXztHQ63WkDXhdQ+lwGm/+YyrSH0vmf9JSAP29TXSmwAngU2O69Fg38DJtF1dWheYlI41AdHSnndDPEpjSkNSztD1sL4aUcW9Nla6Edj++z2UCj4+HyOEgKolWKujIGDpfCjlOwrRA+L4TthXCqUtG+Nnek0MUDB55J5cJYeDAlhX/+PoOH/xicW5pFwHLgMWCv91ocNovqfpRFJSJn0opOCDlbM8RgfEM7H0dLbX+l147C3qKqXxvY0m7NDI61qz6tgyzEP+22f6Y9RfDVadh5Cr48VbVYn0/LcHueaUgrW534261gwbzg3tLMwfaiego44r3WHpgF/BwIvX/NIlKZtq7qKJQDHaeaIQYi492+WXPCtiXYWlg1kygM6BkNg1vZAKhXDCRHQ7zDwU+RBw4UV4zsEhvcfF1km2TWxAX0jLEdwAfF2qCmZ3TNdW2CcUtzO/AE8Dy28SbYbalfAneiSsYizYW2rqTJmyEGsrAwG7z0ioGfdIJcb22YTSft9s7+YrsysqcI/lHpfnEum67eNQraR9hDue0i7IhvAbHh0NIFMeF1679VZqDYA4VuyHdDXhmc9H7MLYWcUjhSYj8/Umo/1qZtCxvE9ImBfi1tdlTPaIiqw/maYNrSNNjqxU8Ab1a6Phz4L+D7gDqCiEhdaUVHmp1jpbClED4vsNtAe2pZMalJGDbYiQizAU+Y91p4WEVwU+yxnbHrq5XL9ovqEmVHj2i72tSQFadg2dIsBl7EBjhbvNfCgEnYLapLUTdxkeZKKzoi9ZAQAWPj7fA57YasYrtVdNC7ypJbCjnez/PccMpt05kNVQ/8nosLaNPCVhaOa2HbJyREQJJ31ah9hB2do+yqkj/r1tQU1NRUc8lJB4AlwDPAYe+1ltgU8fuAvg7NS0RCgwIdESDGBf1b2nE2xtiWCIXeoKfU2KDHYCsMe7CVg6PC7HZSZLj9PDrcuaJ7gbqlaYAPgCeB16hY/eoC3Av8FDXaFBH/0NaViDSZk8AL2OypzytdHw3cjc7fiEjNtHUlIgHtU2ARNsgp8F6LxTbZnA58y6F5iUjoU6AjIo0iD3gJWApsqHS9P7aC8U+A+BruJyLiTwp0RMRvDLAOWAy8Apz2Xo/Abkv9DBiDsqdEpOko0BGRBjuI3ZZ6loreUwADsSs3twFJDsxLRESBjoicl9PYjKk/A//EZp2BTQ2/BZs5NQKt3oiIsxToiEidlWGrFq8A/o49h+MzCrtycwu20aaISCBQoCMitTLYSsV/Bv5CRVE/gB7AFGyAo8J+IhKIFOiISI2+wGZNvez93CcRmAz8J3AZUIdWWyIijlGgIyKAXbnZjt2SehnYWulrkcB12NWba7y3RUSCgQIdkWbMA3yMDW7+Aeyq9LUI4GrgB9ggR+duRCQYKdARaWZKgDXA69isqX2VvhYJ/AdwE7ZruPpNiUiwU6ATwtLS0nC5XDV2p87IyMDtdpOWlnZej50PnAI6NmiG0lSOAiuxwc1b2L8/n1bAROBG7ApO6yafnYhI49E5whDmcrlITU0lIyOjyvWMjAxSU1NxuVzn/diPAb2BVGyjRgksZdgKxXOwtWzaA7di08LzgQ7AndjAJwdb7O8mFOSISOjRik4I863kpKamlt/2BTnp6ek1rvTUhQH+jV3RyQCeAX4NTMMWixNn7Mau2qwEVlO1xg3AIOBa4HrgEvS/HBFpHsKMMcbpSTSVhrR5D2a+4CYyMpKSkpIGBTk+BnuAdTbwlfdaEjbYuQvo0qBHl7o4gj1rsxob3Oyu9vW2wFXY7ajxQNcmnZ2IiP805P1bgU4zERUVRUlJCZGRkRQXF/vtcUuB54B5wF7vNRdwAzAdNXD0p33Y7ag1wAdUrW0Ddnl2FDaouQoYhv27EBEJdg15/9bqdTOQkZFRHuSUlJSccWanISKwZz2+xJ7/uBxwA38FxgHfAh7hzNUGqd0pbFDz39jifF2B7tj2CgupCHIGAzOw2VNHsUHQQ9itKQU5IiI6oxPyqp/J8d0GGrx9VVkEcLN3bMG+GT+PfUN+yDuGYA+83gj0Rys9Pm5soPixd6wHPvder8wFXIgNJsd6PyY22SxFRIKTAp0QVtPB45oOKPvbYGyg8yh2Zecl4D1gs3c8hF2duNI7rgA6+X0WgekktvpwpndsxgaGp2v43g7YjKmR3jEciG2aaYqIhAwFOiHM7XbXePDYd9vtrr5m4F9xwFTvyMFur/wV2/06C1jmHQAXABdj38yHARcRvG/qbmA/drtuF3ZVa7t37D/LfVoCQ4HvYIOb76DDwyIi/hA0h5HnzZvHm2++SWZmJpGRkZw4caLej9GcDyMHkkLgI+wqj2+lp/oPYTjQDxsA9cNudfUD+mCzu5w8XFYEHMIGa994P+7zjj3eUVLL/TtiV72GYAO6Idg/l87UiIjUrCHv30GzolNSUsLkyZMZOXIkf/rTn5yejjRALDbl+Wrv7aPYujybvGMjcBDY4R3VRWC3urpgVz06YlOp2wLx3o+tgWggqtJoge3t5MEGVh5s1tipSuM0tqDeMeB4pXHEO6eDQF1C7AigF7aoYn9goHcMQG0VRESaUtAEOnPnzgVg+fLlzk6kmWrMdhKJ2EJ211a6dhB7duVL79jp/ZiFDU6yvMMpUdhAKxnoAXTzjmTs6kw3tEIjIhIIgibQOR/FxcVVasbk5+fX8t1SG187Cah6gLnygWd/6uQdV1e7XordNjpQaRymYuXlhPfjSaC42ijDBh9h2K2vMOzKS0vviPF+bI1ddUmo9LFdpTl1xK4cKWtMRCTwhXSgM3/+/PKVIGmYxmonUV8RVKyeiIiInIujBQPT0tIICwurdWzcuPG8H3/27Nnk5eWVj3379vlx9s1PSkoK6enppKamEhUV1eRBjoiISH05mnWVm5tLbm5urd+TnJxMdHR0+e3ly5czc+ZMZV05qLHaSYiIiNQkaLOu2rVrR7t27ZycgtRTTe0ktKIjIiKBKmh6XWVlZZGZmUlWVhZut5vMzEwyMzMpKChwemrNRuUzOcXFxeXbWP7snSUiIuJPQXMYOTU1leeee6789pAhQwBYvXo1Y8eOdWhWzYdT7SREREQaImgCneXLl6uGjoOcbichIiJyPoKmBYQ/6DCyiIhI8GnI+3fQnNERERERqS8FOiIiIhKyFOiIiIhIyFKgIyIiIiFLgY5XWlraWevBZGRknHdnbhEREXGOAh0vX3fu6sGOr36My+VyaGYiIiJyvoKmjk5jC5Tu3CIiIuI/qqNTjS+48fVyUpAjIiLirIbU0VGgUwN15xYREQkcKhjoRzV15xYREZHgpECnEnXnltooM09EJPjoMLKXunPLufgy86Dqz0Llnx0REQksCnS81J1bzkWZeSIiwUeHkUXqSZl5IiJNS1lXdaRAR/xFmXkiIk1HWVciTUiZeSIiwUOBjkg9KDNPRCS46DCySB0pM09EJPgo0BGpI2XmiYgEHx1GFhERkYCmw8gScFRFWEREAoECHWkUvirC1YMd3zkXl8vl0MxERKQ50RkdaRSqIiwiIoFAZ3SkUamKsIiINJQqI9eRAh1nqIqwiIg0hA4jS8BSFWEREXGSAh1pNKoiLCIiTtNhZGkUqiIsIiKBQIGONApVERYRkUCgw8giIiIS0HQYWURERKQGCnREREQkZCnQERERkZClQEdERERClgIdERERCVkKdERERCRkKdARERGRkKVAR0REREKWAh0REREJWQp0REREJGQp0BEREZGQpUBHREREQpYCHREREQlZCnREREQkZCnQERERkZClQEdERERClgIdERERCVkKdERERCRkKdARERGRkKVAR0REREJWUAQ6e/fuZerUqfTs2ZOYmBh69+7NnDlzKCkpcXpqIiIiEsBaOD2ButixYwcej4dFixbRp08ftm7dyk9/+lMKCwt5/PHHnZ6eiIiIBKgwY4xxehLnY8GCBSxcuJA9e/bU+T75+fnExcWRl5dHmzZtGnF2IiIi4i8Nef8OihWdmuTl5ZGQkFDr9xQXF1NcXFzlPmBfMBEREQkOvvft81qbMUFo165dpk2bNmbJkiW1ft+cOXMMoKGhoaGhoRECY/fu3fWOGRzdukpLS2Pu3Lm1fs+GDRsYPnx4+e3s7GzGjBnDmDFjWLp0aa33rb6ic+LECXr06EFWVhZxcXENm3wzlp+fT7du3di3b5+2ABtIr6X/6LX0D72O/qPX0n/y8vLo3r07x48fJz4+vl73dXTr6p577uGWW26p9XuSk5PLP8/OzmbcuHGMHDmSxYsXn/Pxo6KiiIqKOuN6XFycfuj8oE2bNnod/USvpf/otfQPvY7+o9fSf8LD658s7mig065dO9q1a1en7z1w4ADjxo1j2LBhLFu27Lz+sCIiItK8BMVh5OzsbMaOHUv37t15/PHHycnJKf9ax44dHZyZiIiIBLKgCHRWrlzJrl272LVrF127dq3ytfocMYqKimLOnDk1bmdJ3el19B+9lv6j19I/9Dr6j15L/2nIaxm0dXREREREzkUHXURERCRkKdARERGRkKVAR0REREKWAh0REREJWc020Lnuuuvo3r070dHRdOrUiSlTppCdne30tILO3r17mTp1Kj179iQmJobevXszZ84cSkpKnJ5a0Jk3bx6jRo2iZcuW9a782dw9/fTT9OzZk+joaIYNG8aHH37o9JSC0tq1a5k4cSKdO3cmLCyMf/zjH05PKSjNnz+fiy++mNatW5OUlMSkSZPYuXOn09MKOgsXLmTw4MHlBRdHjhzJ22+/Xe/HabaBzrhx43j55ZfZuXMnf/vb39i9ezc33XST09MKOjt27MDj8bBo0SK2bdvGE088wTPPPMODDz7o9NSCTklJCZMnT+buu+92eipB5aWXXmLmzJk89NBDbN68mcsvv5wJEyaQlZXl9NSCTmFhIRdeeCFPPvmk01MJamvWrGH69OmsX7+eVatWUVZWxvjx4yksLHR6akGla9euPProo2zcuJGNGzdyxRVXcP3117Nt27Z6PY7Sy71ef/11Jk2aRHFxMREREU5PJ6gtWLCAhQsXsmfPHqenEpSWL1/OzJkzOXHihNNTCQojRoxg6NChLFy4sPzagAEDmDRpEvPnz3dwZsEtLCyMV199lUmTJjk9laCXk5NDUlISa9asYfTo0U5PJ6glJCSwYMECpk6dWuf7NNsVncqOHTvGX/7yF0aNGqUgxw/y8vJISEhwehrSDJSUlLBp0ybGjx9f5fr48eNZt26dQ7MSqSovLw9AvxcbwO12s2LFCgoLCxk5cmS97tusA50HHniA2NhYEhMTycrK4rXXXnN6SkFv9+7d/PGPf2TatGlOT0WagdzcXNxuNx06dKhyvUOHDhw6dMihWYlUMMYwa9YsLrvsMgYNGuT0dILO559/TqtWrYiKimLatGm8+uqrDBw4sF6PEVKBTlpaGmFhYbWOjRs3ln//L3/5SzZv3szKlStxuVzcdttt9WopEcrq+1qC7Ul2zTXXMHnyZO68806HZh5Yzud1lPoLCwurctsYc8Y1ESfcc889bNmyhRdffNHpqQSl/v37k5mZyfr167n77ru5/fbb2b59e70eIyh6XdXVPffcwy233FLr9yQnJ5d/7uue3q9fPwYMGEC3bt1Yv359vZfFQlF9X8vs7GzGjRvHyJEjWbx4cSPPLnjU93WU+mnXrh0ul+uM1ZsjR46cscoj0tRmzJjB66+/ztq1a8/o0yh1ExkZSZ8+fQAYPnw4GzZs4A9/+AOLFi2q82OEVKDjC1zOh28lp7i42J9TClr1eS0PHDjAuHHjGDZsGMuWLSM8PKQWChukIT+Tcm6RkZEMGzaMVatWccMNN5RfX7VqFddff72DM5PmzBjDjBkzePXVV/nggw/o2bOn01MKGcaYer9Ph1SgU1effPIJn3zyCZdddhlt27Zlz549pKam0rt3b63m1FN2djZjx46le/fuPP744+Tk5JR/rWPHjg7OLPhkZWVx7NgxsrKycLvdZGZmAtCnTx9atWrl8OwC16xZs5gyZQrDhw8vX1HMysrSObHzUFBQwK5du8pvf/3112RmZpKQkED37t0dnFlwmT59Oi+88AKvvfYarVu3Ll9xjIuLIyYmxuHZBY8HH3yQCRMm0K1bN06ePMmKFSv44IMPeOedd+r3QKYZ2rJlixk3bpxJSEgwUVFRJjk52UybNs3s37/f6akFnWXLlhmgxiH1c/vtt9f4Oq5evdrpqQW8p556yvTo0cNERkaaoUOHmjVr1jg9paC0evXqGn8Gb7/9dqenFlTO9jtx2bJlTk8tqNxxxx3l/67bt29vrrzySrNy5cp6P47q6IiIiEjI0mEKERERCVkKdERERCRkKdARERGRkKVAR0REREKWAh0REREJWQp0REREJGQp0BEREZGQpUBHREREQpYCHREREQlZCnREREQkZCnQEZGg9s4773DZZZcRHx9PYmIi1157Lbt373Z6WiISIBToiEhQKywsZNasWWzYsIH33nuP8PBwbrjhBjwej9NTE5EAoKaeIhJScnJySEpK4vPPP2fQoEFOT0dEHKYVHREJart37+bWW2+lV69etGnThp49ewKQlZXl8MxEJBC0cHoCIiINMXHiRLp168aSJUvo3LkzHo+HQYMGUVJS4vTURCQAKNARkaB19OhRvvjiCxYtWsTll18OwEcffeTwrEQkkCjQEZGg1bZtWxITE1m8eDGdOnUiKyuLX//6105PS0QCiM7oiEjQCg8PZ8WKFWzatIlBgwZx//33s2DBAqenJSIBRFlXIiIiErK0oiMiIiIhS4GOiIiIhCwFOiIiIhKyFOiIiIhIyFKgIyIiIiFLgY6IiIiELAU6IiIiErIU6IiIiEjIUqAjIiIiIUuBjoiIiIQsBToiIiISshToiIiISMj6f2HDROHMBJlHAAAAAElFTkSuQmCC",
"text/plain": [
"PyPlot.Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"p(x⋅∈S|m) ≈ 0.1788220649957509\n"
]
}
],
"source": [
"using Cubature # Numerical integration package\n",
"\n",
"# Maximum likelihood estimation of 2D Gaussian\n",
"μ = 1/N * sum(D,2)[:,1]\n",
"D_min_μ = D - repmat(μ, 1, N)\n",
"Σ = Hermitian(1/N * D_min_μ*D_min_μ')\n",
"m = MvNormal(μ, convert(Matrix, Σ))\n",
"\n",
"# Contour plot of estimated Gaussian density\n",
"A = Matrix{Float64}(100,100); B = Matrix{Float64}(100,100)\n",
"density = Matrix{Float64}(100,100)\n",
"for i=1:100\n",
" for j=1:100\n",
" A[i,j] = a = (i-1)*6/100.-2\n",
" B[i,j] = b = (j-1)*6/100.-3\n",
" density[i,j] = pdf(m, [a,b])\n",
" end\n",
"end\n",
"c = contour(A, B, density, 6, zorder=1)\n",
"PyPlot.set_cmap(\"cool\")\n",
"clabel(c, inline=1, fontsize=10)\n",
"\n",
"# Plot observations, x∙, and the countours of the estimated Gausian density\n",
"plotObservations(D)\n",
"plot(x_dot[1], x_dot[2], \"ro\")\n",
"\n",
"# Numerical integration of p(x|m) over S:\n",
"(val,err) = hcubature((x)->pdf(m,x), [0., 1.], [2., 2.])\n",
"println(\"p(x⋅∈S|m) ≈ $(val)\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Discrete Data: the 1-of-K Coding Scheme\n",
"\n",
"- Consider a coin-tossing experiment with outcomes $x \\in\\{0,1\\}$ (tail and head) and let $0\\leq \\mu \\leq 1$ represent the probability of heads. This model can written as a **Bernoulli distribution**:\n",
"$$ \n",
"p(x|\\mu) = \\mu^{x}(1-\\mu)^{1-x}\n",
"$$\n",
" - Note that the variable $x$ acts as a (binary) **selector** for the tail or head probabilities. Think of this as an 'if'-statement in programming. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- **1-of-K coding scheme**. Now consider a $K$-sided coin (a _die_ (pl.: dice)). It is convenient to code the outcomes by $x=(x_1,\\ldots,x_K)^T$ with **binary selection variables**\n",
"$$\n",
"x_k = \\begin{cases} 1 & \\text{if die landed on $k$th face}\\\\\n",
"0 & \\text{otherwise} \\end{cases}\n",
"$$\n",
" - E.g., for $K=6$, if the die lands on the 3rd face $\\,\\Rightarrow x=(0,0,1,0,0,0)^T$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Assume the probabilities $p(x_k=1) = \\mu_k$ with $\\sum_k \\mu_k = 1$. The data generating distribution is then (note the similarity to the Bernoulli distribution)\n",
"\n",
"$$\n",
"p(x|\\mu) = \\mu_1^{x_1} \\mu_2^{x_2} \\cdots \\mu_k^{x_k}=\\prod_k \\mu_k^{x_k}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- This generalized Bernoulli distribution is called the **categorical distribution** (or sometimes the 'multi-noulli' distribution).\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Categorical vs. Multinomial Distribution\n",
"\n",
"- Observe a data set $D=\\{x_1,\\ldots,x_N\\}$ of $N$ IID rolls of a $K$-sided die, with generating PDF\n",
"$$\n",
"p(D|\\mu) = \\prod_n \\prod_k \\mu_k^{x_{nk}} = \\prod_k \\mu_k^{\\sum_n x_{nk}} = \\prod_k \\mu_k^{m_k}\n",
"$$\n",
"where $m_k= \\sum_n x_{nk}$ is the total number of occurrences that we 'threw' $k$ eyes."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- This distribution depends on the observations **only** through the quantities $\\{m_k\\}$, with generally $K \\ll N$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- A related distribution is the distribution over $D_m=\\{m_1,\\ldots,m_K\\}$, which is called the **multinomial distribution**,\n",
"$$\n",
"p(D_m|\\mu) =\\frac{N!}{m_1! m_2!\\ldots m_K!} \\,\\prod_k \\mu_k^{m_k}\\,.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The catagorical distribution $p(D|\\mu) = p(\\,x_1,\\ldots,x_N\\,|\\,\\mu\\,)$ is a distribution over the **observations** $\\{x_1,\\ldots,x_N\\}$, whereas the multinomial distribution $p(D_m|\\mu) = p(\\,m_1,\\ldots,m_K\\,|\\,\\mu\\,)$ is a distribution over the **data frequencies** $\\{m_1,\\ldots,m_K\\}$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Maximum Likelihood Estimation for the Multinomial\n",
"\n",
"- Now let's find the ML estimate for $\\mu$, based on $N$ throws of a $K$-sided die. Again we use the shorthand $m_k \\triangleq \\sum_n x_{nk}$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- The log-likelihood for the multinomial distribution is given by\n",
"\n",
"$$\\begin{align*}\n",
"\\mathrm{L}(\\mu) &\\triangleq \\log p(D_m|\\mu) \\propto \\log \\prod_k \\mu_k^{m_k} = \\sum_k m_k \\log \\mu_k \\tag{2}\n",
"\\end{align*}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- When doing ML estimation, we must obey the constraint $\\sum_k \\mu_k = 1$, which can be accomplished by a Lagrange multiplier. The **augmented log-likelihood** with Lagrange multiplier is then\n",
"\n",
"$$\n",
"\\mathrm{L}^\\prime(\\mu) = \\sum_k m_k \\log \\mu_k + \\lambda \\cdot (1 - \\sum_k \\mu_k )\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- Set derivative to zero yields the **sample proportion** for $\\mu_k$ \n",
"\n",
"$$\\begin{equation*}\n",
"\\nabla_{\\mu_k} \\mathrm{L}^\\prime = \\frac{m_k }\n",
"{\\hat\\mu_k } - \\lambda \\overset{!}{=} 0 \\; \\Rightarrow \\; \\boxed{\\hat\\mu_k = \\frac{m_k }{N}}\n",
"\\end{equation*}$$\n",
"\n",
"where we get $\\lambda$ from the constraint \n",
"\n",
"$$\\begin{equation*}\n",
"\\sum_k \\hat \\mu_k = \\sum_k \\frac{m_k}\n",
"{\\lambda} = \\frac{N}{\\lambda} \\overset{!}{=} 1\n",
"\\end{equation*}$$\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Recap ML for Density Estimation\n",
"\n",
"Given $N$ IID observations $D=\\{x_1,\\dotsc,x_N\\}$.\n",
"\n",
"- For a **multivariate Gaussian** model $p(x_n|\\theta) = \\mathcal{N}(x_n|\\mu,\\Sigma)$, we obtain ML estimates\n",
"\n",
"$$\\begin{align}\n",
"\\hat \\mu &= \\frac{1}{N} \\sum_n x_n \\tag{sample mean} \\\\\n",
"\\hat \\Sigma &= \\frac{1}{N} \\sum_n (x_n-\\hat\\mu)(x_n - \\hat \\mu)^T \\tag{sample variance}\n",
"\\end{align}$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- For discrete outcomes modeled by a 1-of-K **categorical distribution** we find\n",
"\n",
"$$\\begin{align}\n",
"\\hat\\mu_k = \\frac{1}{N} \\sum_n x_{nk} \\quad \\left(= \\frac{m_k}{N} \\right) \\tag{sample proportion}\n",
"\\end{align}$$\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
" \n",
"- Note the similarity for the means between discrete and continuous data. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"- We didn't use a co-variance matrix for discrete data. Why?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"---\n",
"The cell below loads the style file."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"open(\"../../styles/aipstyle.html\") do f\n",
" display(\"text/html\", readstring(f))\n",
"end"
]
}
],
"metadata": {
"anaconda-cloud": {},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 0.6.1",
"language": "julia",
"name": "julia-0.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.6.1"
}
},
"nbformat": 4,
"nbformat_minor": 1
}