{
"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+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X+UVPV9//HX7siMIuwqEA2WRYm2amLNqZL24AFcQI2eLYsmMcc0x9qE9EQCRMDTI9rO7jBDyh4RaE9LKP44pn80oU39BaGhckRYOK0NWDkxGjViyBLRRDDO8iPO1Nn7/YPvnczs3pmdX3fuvZ/7fJyzJ5nZ+fGZy7qf134+78/n02JZliUAAAADtXrdAAAAALcQdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQIbdNasWaOWlhYtW7bM66YAAACfCmTQ2b9/vx5++GFdffXVXjcFAAD4WOCCzsmTJ/XlL39ZjzzyiM4//3yvmwMAAHzsLK8bUK3Fixerq6tLN9xwg1avXl32sZlMRplMJn97aGhI77//viZOnKiWlha3mwoAABrAsiydOHFCF110kVpbqxujCVTQ2bJli/73f/9X+/fvr+jxa9as0apVq1xuFQAAaIYjR45oypQpVT0nMEHnyJEjuueee/Tss8/q7LPPrug5999/v1asWJG/nU6nNXXqVB05ckRtbW1uNRUAADTQ4OCgOjo6NH78+Kqf22JZluVCmxru6aef1m233aZIJJK/L5fLqaWlRa2trcpkMkXfczI4OKj29nal02mCDgAAAVFP/x2YEZ158+bp5ZdfLrrvK1/5iq644grdd999o4YcAAAQPoEJOuPHj9dVV11VdN+5556riRMnjrgfAABACuDycgAAgEoFZkTHye7du71uAgAA8DFGdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAgyQSCaVSKcfvpVIpJRKJ5jbIYwQdAAAMEolE1NPTMyLspFIp9fT0hG6D3UAvLwcAAMXi8bgkqaenJ3/bDjnJZDL//bAIzFlXjcBZVwCAsLDDTTQaVTabDXTIqaf/JugAAGCoWCymbDaraDSqTCbjdXNqVk//TY0OAAAGSqVS+ZCTzWZLFiibjqADAIBhCmtyMpmMksmkY4FyGFCMDACAQZwKj50KlMOCoAMAgEFyuZxj4bF9O5fLedEsz1CMDAAAfI1iZAAAAAcEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAAIhkUgolUo5fi+VSimRSDS3QQgEgg4AIBAikYh6enpGhJ1UKqWenh5FIhGPWgY/O8vrBgAAUIl4PC5J6unpyd+2Q04ymcx/HyjUYlmW5XUjmmVwcFDt7e1Kp9Nqa2vzujkAgBrY4SYajSqbzRJyQqCe/pugAwAInFgspmw2q2g0qkwm43Vz4LJ6+m9qdAAAgZJKpfIhJ5vNlixQBiSCDgAgQAprcjKZjJLJpGOBMmCjGBkAEAhOhcdOBcpAIYIOABgikUgoEok4dvapVEq5XC7Qe83kcjnHwmP7di6X86JZ8DmCDgAYwt5nRioe2SgcCQmyciGNkRyUQtABAEOwzwwwEsvLAcAw7DMD07CPToUIOgDCgn1mYBL20QEA5LmxzwwHaiKoCDoAYBC39pnhQE0EFcXIAGAIN/eZodAZQUXQAQBDuL3PTGHYWb16NYXOCASKkQEAVaHQGc1GMTIAoCk4UBNBQ9ABAFSEAzURRNToAABGxYGaCCqCDgBgVByoiaCiGBkAAPhaKIqRN23apKuvvlptbW1qa2vTjBkz9MMf/tDrZgEAAB8LTNCZMmWK+vr6dODAAR04cEBz587VggUL9Morr3jdNAAA4FOBnrqaMGGC1q5dq4ULF1b0eKauAAAInnr670AWI+dyOX3/+9/XqVOnNGPGjJKPy2QyRZtZDQ4ONqN5AADAJwIzdSVJL7/8ssaNG6dYLKa7775bTz31lD75yU+WfPyaNWvU3t6e/+ro6GhiawEAgNcCNXWVzWY1MDCgDz74QE888YQeffRR7dmzp2TYcRrR6ejoYOoKAIAAqWfqKlBBZ7gbbrhBl156qTZv3lzR46nRAQAgeEKxvNyJZVkcKAcAAEoKTDHyAw88oFtuuUUdHR06ceKEtmzZot27d2vHjh1eNw0AAPhUYILOr371K915551655131N7erquvvlo7duzQjTfe6HXTAACATwUm6Dz22GNeNwEA4BOJREKRSMTxINFUKqVcLqdEItH8hsF3Al2jAwAIp0gkop6eHqVSqaL77VPWI5GIRy2D3wRmRAcAAJs9ktPT05O/bYccp1PWEV6BXl5eLZaXA4BZ7HATjUaVzWYJOYYK7T461SLoAIB5YrGYstmsotEoW44YKrT76AAAwi2VSuVDTjabHVGzAxB0AACBVFiTk8lklEwmHQuUEW4UIwMAAsep8NipQBkg6AAAAieXyzkWHtu3c7mcF82CD1GMDAAAfI1iZAAAAAcEHQBAwyUSiZJFwalUiuMZ0DQEHQBAw3FEA/yCYmQAQMNxRAP8gmJkAIBrOKIBjcAREBUi6ABA83FEA+rFqisAgC9xRAO8RtABALiCIxrgBxQjAwAajiMaGiORSCgSiTheq1QqpVwux1L9URB0AAANxxENjWEv05eKg2FhkER5FCMDAOBjw0fHwrhMn1VXFSLoAPCzRk9TMO1hjrAv02fVFQAYoNG7CbM7sTni8Xg+5ESj0VCFnLpZIZJOpy1JVjqd9ropAOAomUxakqxkMul42+vXgzfsf7doNBrKf796+m+CDgD4TKM7tbB3kkFHWK2v/6ZGBwB8qNG7CbM7cTCVKjwOW0EyNToAYJBSuwknEomSm+2lUqmShcXsThxc5ZbpJ5NJlulXouHjSz7G1BUAvys3TVFqyqLcVAbTHjABNToVIugA8LPCENLb21s27HR2do74XrnXq+R+wK/q6b/ZGRkAfKJwmsKuwZCKdxO2l4Tv3r07X3djP2f4vjnDX8/eN4fdiREqLgQv32JEB0CQlJt2sldQRaPRko8f7X4gKBjRAQADFR6CuXr16vzojaQRxcXxeNzx0Mywrc4BhmPVFQD42PAdcSXlg0smk1EymSza/dhejdPT06NYLEbIQeixjw6AQArLOU7DzziSVNGeKuybA5Owjw6A0BntHKe9e/c6Pq/cfjN+UxhgMpmMOjs7HR83fE8V9s0BCjS8YsjHKEYGzFKqWHfu3LmBL8qttbCYfXNgIvbRqRBBBzBPqXOcgt7hF+6jM5y9z47T/UEPeIATzrqqEDU6gJlK1aMMr28xvSg3LHVLCJ96+m+CDoBAGy3MUJQLBB/FyABCaXix7vCl1hTlAmDDQACB5LSkunDDvN27d2vXrl1lj1QAYD6mrgD40mj1Js8995zmzZvn+P158+YVhZzC5/X09Gju3Ll67rnnHF+3EXUs1MoAjcXUFQDjjLZPTqmQI0mzZs1yLDyOx+OaO3eudu3aVfJ17UMz3Wz78PdIJBIjHmvf57TvT5D2AgI819D1Xz7H8nIgWNxaIt6MpefVvIfT9+z7St3PUnGECfvoVIigAwRPqX1y/Pq6tb5HqWBUSyCrZQ8ewM8IOhUi6ADBZAeFaDQaiNet9T2cglEtgYyNA2Eagk6FCDpA8IRlRMfmFIxqCWRe7AzNSBLcQtCpEEEHCJaw1OgMf069IzrlXs9NfhlJInCZh6BTIYIO0BjN6Ejc6jSb0RnX8h6NrNEp1IzpuUJ+OGPML4ELjVNP/82GgQCqZi+fluS4T00ymaz7PXK5XMkl4vb3/fS69byH0+aHTgo3RCy8XYrTztBub5ZY2MbVq1d7csaY03Wq9BrDQC4EL99iRAdoHD/85W4KpxEy+z6nEbJKRs28/vdp9kiSk2ZP3cE9oZi6+tu//Vtr+vTp1rhx46yPfexj1oIFC6zXXnutqtcg6ACNRUfiT15P3fjp58IPgQv1q6f/DszOyHv27NHixYv1wgsvaOfOnfroo49000036dSpU143DTCa0669hVpbW/PTI0wJ+EO5qbNkMtmQ6blSRjtotZk41BWSgjt19etf/9qSZO3Zs6fi5zCiA1RvtNEBVfGXO6thzOb1SFK592RqNdhCWYycTqclSRMmTCj5mEwmo0wmk789ODjoersA05Qr7JRU1engbhcxp9Np/d///V9dr4HanThxQitXrtSiRYt07Nix/P2LFi3S6dOndeLEiaL73bJu3Tr19fUVtcVuQ09Pj06fPq17773X9XZI0pgxY9Te3t6U94KzQJ5eblmWFixYoN/85jfau3dvycclEgmtWrVqxP2cXg5Uzw4j9jSApJKng5db2TL8MY1aDZNOp7V58+b8H0EIr/7+frW2tmrmzJkjvrdv3z4NDQ1p9uzZTWlLe3u7vv71rxN26lTP6eWBDDqLFy/W9u3btW/fPk2ZMqXk45xGdDo6Ogg6QI1isZiy2awikYh6e3sdg0kqlVIulyt7uvbw0NSIJb/Hjh3Thg0bdM4552js2LF1vRbQCKdPn9Zvf/tbLV++XJMmTSr63v/8z/+or69PL774on71q1/pvPPO0yc+8Qldd911WrdunUct9q96gk7gpq6WLl2qrVu3qr+/v2zIkc78Uo7FYk1qGWC24YWdpVQSWOLxeH6PlUYXMY8dO1bjxo1r2Os107Zt29Ta2qqurq4R39u+fbuGhoY0f/58D1qGWv32t78dcd/27dvV3d2tzs5OPfjgg5o8ebLeeecdHThwQFu2bPFt0EkkEopEIjX/geOVwKy6sixLS5Ys0ZNPPqldu3Zp2rRpXjcJCI1Gr6RhNYyz1tZWbd26Vdu3by+6f/v27dq6dataWwPzKxtlPPjgg5o2bZr+8z//U3fccYeuv/563XHHHXrooYc0MDDgdfNKsmvshv/3av9+iEQiHrWsvMCM6CxevFjf/e539cwzz2j8+PF69913JZ2Z/zznnHM8bh1gLqcammp36C33epUUMYeFPZKzdevW/G075HR3dzuO9CB4jh8/rkmTJumss0Z2wX4Os0HdcTowQWfTpk2SpM7OzqL7H3/8cf3FX/xF8xsEhEQjj0xodGgyUWHY+Y//+A999NFHhBzDzJgxQ48++qi++c1v6stf/rKuueYajRkzxutmVcQPR3xUK5DFyLWqp5gJQP3cnOO3i5EnTpzY0Bodr+pmFi9erI8++khnnXWWNm7c2PDXh/tOnjyp48ePjyhGPn78uG699Vbt27dP0pkl6J/5zGc0f/58LVmyJBA1ZvbChGg0WrToxy319N/+HSMDYJxEIlHyL794PO7PQkYP6ma2b9+eDzkfffTRiPdGsE2cOFF79+7V/v371dfXpwULFuiNN97Q/fffrz/8wz9syl5D9QhajR1BBwDK6OrqUnd3d1HYcbNupvC1N27cOOK9YY7p06frvvvu0/e//30dPXpUy5cv1+HDh/Xggw963bSS/HTER6UCU6MDAF5pVt2MU4ByKlCGecaMGaPe3l5t2LBBP/nJT7xujqOg1tgRdACgAl1dXfmQc9ZZZ7kSOIaGhhwDlH17aGio4e/pFvYEKu2dd97R5MmTR9z/05/+VJJ00UUXNbtJFWnkwoRmIugAQAWc6mYaHXbKdfxBG8mxa5uk4rYXjlqF1Wc/+1lNmTJF8+fP1xVXXKGhoSEdPHhQ69at07hx43TPPfd43URH5Wro/DiSYyPoAMAohk8p2bel4AWQZmFPoNL+5m/+Rs8884w2bNigd955R5lMRpMnT9YNN9yg+++/X1deeaXXTTQKQQcAyqBupnbsCeTsi1/8or74xS963YzQIOgAQBkm1c14oRm1TUA5BB0AKMOkuhkvNKO2CSiHoAMAcAW1TfADgg4AoOGobYJfEHQAAA1HbRP8gqADAGg4apvgF5x1BQAAjEXQAWCkbdu2lTwIc/v27dq2bVuTWwTACwQdAEayjyBYt25dUeCxi2RbW1uNDDwEPKAYNToAjFS4wueNN97I3194zpKJZy5xxhRQjKADwFiFYcfu/IeHnGYXxrp9qjdnTAHFCDoAjFZ4BIEkz89casaIC2dMAb9D0AFgtOFHEHh95lKzRlw4Ywo4g6ADwFjDR0nscOHWmUuVTkvVOuJSzbQXZ0wBZ7DqCoCRnEJOd3d30e1Sq5NqsW3bNv3sZz9zfN3169fnV3rZurq68iGk0hEXe9pr+OsXriQrvN3d3a2NGzequ7u74Z8XCApGdAD4WiKRUCQSUTweH/G9VCqlXC6nRCIx4nv2EQSSc+Hxa6+91tAzl1pbW/X666/r8ssvL3rd9evX5+8fXpNT7YhLJdNenDEFFCPoAPC1SCSinp4eScqHnUQiob1792rXrl1KJpNFj9+5c6fOOuus/BTOtm3bRoScrq6ufCho1JlLhWHCDjs/+MEPNDQ0pMsvv1wrVqzIP7aeU71Hm/bijCmgGEEHgK/Z4aYw7NghZ+7cuUUjPfv27VN/f3/RyqVmnrk0fORkaGhIra2tZUOO0/MqCTulCo05YwooRo0OAN+Lx+NKJpPq6elRLBbLh5xdu3YplUpJktatW6f+/n7dfPPNnnboXV1dRbU4Q0NDRbUx5UZcuru7KxpxcZr2CiJ2cUYzEHQA+EoikciHl0LxeFyRSETZbFbRaFTPPfdcUfjp6+vT7NmzdeONN3rQ6t9Zv359fiRHUn4ay+7QC1ddDdfV1TXqZoEmFRpXWlwN1IOfIgC+YtfkDA878+bNUy6Xy4edVCqleDyu1tbWfPiZOXNm0XOaPSpQWHi8adMmdXd3FxUo1xtGSk17BTXsOLWdXZzRaNToAPAVp5qcefPm5aernnvuOaVSKfX09Gj37t35qZ5sNqt9+/ZpwYIFkqrfabjeoxm2b9+eDzV2Tc7wAuV6C4FNLDRmF2e4jaADwHcKw86qVauUy+XyIcf+/u7du/Php7OzUz09Perv79fYsWMVjUarHhWo92iGSkJIPWdYSfUVGrt9xlY92MUZbiLoAPCleDyu1atXK5vNKhKJ5EOOdGb/HDvkzJo1S/F4XKdPn1ZfX5927NghSVWPCtR7NIPfVzv5+VRzdnGGmwg6AHwplUrla28Ka3IkKZfLKZlMFi0tv/fee7V27VrlcrmSIxejMXkaxa+nmtezpxBQCYIOAN+xa3DsMGPfls6M9DjthLxu3bp8sXIul6t5VKCaaRQ/Twc58VuQYxdnNANBB4CvDA85knOB8vDn2MvLFyxYoD179tTcUVYzjeLn6aBS/FQPY2JxNfyHoAPAV5ympSQVTVsVsoPRypUr8/uu1DoqUO00il+ng8rxUz2M3+uaYAaCDgBfcZqWsjkd7GkHo0WLFmnDhg35+6sdFah1GsVv00HlUA+DMCLoAAg0OxgdO3ZsxPeq6bzrmUbx03RQKdTDIKwIOgCg+qZR/DQdVAr1MAgrgg4A1CEo00HUwyCsCDoAUCOmgwD/I+gAaKhEIqFIJOJYOJxKpZTL5coWHDeS2/vcMB3kLGj7C8FsnF4OoKFKnT5uLwOPRCL5+xKJxIjHFT6+3kBk73Mz/FRveyTGXo5eq/nz55ccsenq6gptZ+72dQeqwYgO4CE/jX40itPmfk6bAEq/C0WFz5OKNw2sRxD3uTEB1x1+QtABPOR2R++VwrBjH8xZbhPASkJRrYK0z43bmjmlxHWHXzB+CHgoHo8rmUwWTfU0uqP3Sjwezx/IGY1GS36WwmsQi8Vc+exdXV35pd9+3eemGZo9pcR1hx8QdACPVdLRu13L4gan08dLqTQU1cppn5sw6urqUnd3d1HYcXNKiesOPyDoAD4wWkdfTYFvtdwIUYWjUplMZsSoldPjKw1FThKJhNatW+f4vfXr1+c78o0bN47o6E2wbdu2kp9n+/bt2rZtW/52YdhZvHixqyHH9OuOYCDoAD4wWkfv5hRXo0NUqdPHS4WdakNRqc/Q19enffv2Fd2/fv16vf7667r88suL9rkxrdOtdkrK7SmlUvsLmXbdEQwUIwMeGx4M7NtScYFypQW+1Wp0QXA1p4+XCkXD21PJZzh9+rT6+vo0duxY3Xbbbdq+fXs+5KxYsaLo8abtc1PtKie3j6xgfyH4SYtlWZbXjWiWwcFBtbe3K51Oq62tzevmACUDRbmgEYvF8qM/mUym4W2xR5WaUQzdyOX1x44d0+c//3n19/fnO/CwrfKxw025z1/qyIqwXatmOHnypI4fP67ly5dr0qRJXjdHUnC3tKin/2bqCvBQudGPZDJZNPoh1V/LUo7bBcFOEolE2dVY1f7CnTlzpiKRSGhX+Yw2JcWUEtys9/OrQAWd/v5+zZ8/XxdddJFaWlr09NNPe90koEi1hb3VdPSNqGUpx80Q1Sz79u1TLpcL7Sqf0VY5lZtS6u7uZkopBEze0qKUQNXonDp1Sp/+9Kf1la98RZ///Oe9bg4wglsbADaqlqXS1y9VJ+Rn69atU39/v26++eZ8jU6YDtas5BR1TjCH5F69n1/VHXTsEp+Wlpa6GzOaW265Rbfccovr7wPUyq2dfqsp8K2W2yGqGTUBqVRKfX19mj17tm688UZJ4TpFnFPUUa14PJ4POc2aqvZKzUHnscce04YNG/Szn/1MkvT7v//7WrZsmb72ta81rHH1ymQyRcWag4ODHrYGYeHGX0vlgkC9v6DcDFFSc465yOVyWrlypVpbW/WDH/xAsVhMN954o66//npls1l9+OGHOnnypHbu3KmhoSF99rOfrfs9/eTDDz/UpZdeqmw2q5MnT+bvtz//K6+8og8//NC4z+13p0+f9roJJTlNVZsadmoKOvF4XBs2bNDSpUs1Y8YMSdJ///d/a/ny5Tp8+LBWr17d0EbWas2aNVq1apXXzUAIBemvJTdDVOFruHmeVSKRUDqd1ubNm5XL5bRjxw6dPn1aM2fO1DXXXCNJeuaZZ9Tf36/Zs2fr+PHj+ef29/ertbVVM2fOHPG6+/bt09DQkGbPnl13G918n+nTp+vDDz8s+ty206dP69ChQ/q93/u9os+N5mhvb9eYMWO8bkYRE6aqq1FT0Nm0aZMeeeQRfelLX8rf193drauvvlpLly71TdC5//77i/bPGBwcVEdHh4ctQliE6a+lSjSjJqC9vV1f//rX9dWvflXr1q1TX1+frrvuOt177735+p2VK1fq3nvvzT/nwQcfVEtLi/r7+/OPtX3uc5/T3r17tXLlSi1fvrzu9g0NDRW1yVbYtnreZ/ny5RV/bjTPmDFj1N7e7nUz8tyeqvYlqwbnnXee9cYbb4y4//XXX7fa29trecmqSbKeeuqpqp6TTqctSVY6nXapVYBlJZNJS5KVTCYdb4dZNBq1JFnRaLTs43p7e0ter2QyafX29o76XvZ1t9/T6fXsx8ydO7foMfbtuXPnjv6hqtCMn41KPjfCqxH/bXmhnv67pqCzZMkSa/ny5SPuv/fee61vfOMbtbxk1Qg68KNSHZfXYccPv9yq6YAbdR0rCVbDw04kEnEl5Ax/PzeDSKWBEgiKpgSd5cuX57+WLl1qjR8/3vrUpz5lLVy40Fq4cKH1qU99ympra7OWLFlSdSMqdeLECeull16yXnrpJUuStX79euull16yfvGLX1T0fIKOWfzQeQ/nxzbZ7+1lAKtlJKPe0Y9agpX9FYlEqvh01XMziDCiAxM1Jeh0dnZW9DVnzpyqG1Gp559/vuiXkf111113VfR8go5ZvO68g8arKbV6/p1q7bRr+az2SI795dZ1cTOIMG0KUzV96iqoCDrm4Rd7dbz4a3/4KFfh7eGjXE6jXtWOftQSrJymrdy4Pm7+vBL8YTKCToUIOmZiqL46XtdvVNMh1/JvW+304fDC41IFyvVyO4j4ddoUaASCToUIOubyuvMOCr+EwkpGNpq5Qml44XHh/Y0KCAQRoHYEnQoRdMzkl87b7/w2zVfu361Z0zCEDyAYCDoVIuiYx2+dt1/5tX6j1EicHwKIH9oA4Ix6+u9AnV4OFB4Q6bSNuX1uk7E7fNbI7fOsalFu92i3j6WoRDPO6ALQBC4EL99iRCf4Ckcghq/e0bCRHf7i9q+gjMQFpZ2A6Zi6qhBBxwx0PsHm12m0UqgBA7xXT//dYlmW1fxxJG8MDg6qvb1d6XRabW1tXjcHdbCnD+xpD6dpmVoVTo85vW8ulys7tYLygnh9Y7FYfpotk8l43RwgdOrpv1tdahPgqng8ng850Wi0oXUbdm1GKpUqut8OV5FIpGHvFUaJRKLkv1c8HvddyHGqJQIQHAQdBJKbnU88Hs8XNNuvO7zwGeFQ+O+eyWRG/FwACICGT6T5GDU6ZmhWjQ61GeEWtFoiwGQsL0doOI2s2P/b6CXl8Xhcq1evdmV6DP7nxyX5AKpH0EGgNLPzKbfPC2rnh2LketvAzwEQHNToIFCaVchKbYZ7/FDs7Yc2AGgSF6bSfIsaHVSC2gz3+WEvJD+0AUBlqNEBGijotRl+mBoaTWFdlV0H1ewVbX5oAwD3sWEgYJhSS+H9uETeDxvx+aENAMpjw0AAeUHZB8gPG/H5oQ0AXNbwiTQfo0YHYeLnfYD8UB/jhzYAqAyHelaIoAMvFZ62Ppxbp63bIScajdb0fDfa7Idibz+0AUDl6um/mboCmqTZS5obMS3jRpvLFXsnk8mmFHu70oZcTtq9W/re9878r8+L1oHQcCF4+RYjOvBas4+vaMT7MMVTgSeesKwpUyxL+t3XlCln7gdQN6auKkTQgR+4XTvjxrSMn+t9PPfEE5bV0lIccqQz97W0EHaABqin/2Z5OeABN5c0u7WPDsuwHeRy0iWXSL/8pfP3W1qkKVOkn/9cYrdloGYsLwcCxO0lzW4ck8Ey7BL27i0dcqQzYztHjpx5HABPEHSAJgriGVpBbHPTvPNOYx8HoOE4AgK+E4QjDGrhtGlf4TEEhbf9IohtbqrJkxv7OAANR9CB79hLmiWVPMIgiIJ4hlYQ29xUs2adqcF5++0WFi9zAAAUHklEQVQz01TD2TU6s2Y1v20AJHHWFXxq+EiCH48wACRJTz4pfeELZ/5/4a/TlpYz//vv/y597nPNbxdgkHr6b4IOfMsON3YBLCEHvvXkk9I99xQXJnd0SH/3d4QcoAEIOhUi6AQPS5oRGLncmdVV77xzpiZn1iyWlAMNwvJyGIklzQiUSETq7JS+9KUz/0vIAXyBoANf8npJcyKRKPleqVQqkKu+ACCMCDrwnVJLmpsZdpp9AKffEPQAmILl5fAdPyxpdtorJkwrv0xd4g8ghBp45pbvcagnqhXmwyw5tRyAX3CoZ4VYdYVahHnlF0v8AfgBq64Al4R95Vc8Hs9/9mg0SsgBEDgEHaAEr1d++UHYgx6A4KMYGXDAYZalj+GQzP/szWbqQbaAHxB0AAd+WPnlJYJec7HKDXAPQQdwUO6v5zB08GEPes0W9u0MADex6goAfIJVboAzDvWsEEEHgN+FeTsDoBSWlwOAAVjlBjQeQQej4twjwH1sZwC4g2JkjIoVIYC7WOUGuIegg1GxIgSoTrX74rDKDXAPQQcVKQw7q1evZkUIUEa1o6Bh384AcBM1OqgY5x4Bv1Oudk2SOjs7i2psGAUFvBG4oPPtb39b06ZN09lnn61rr71We/fu9bpJoRGmFSEUYGM09qjN8J8TO9DMnTs3X1Aci8UIOYBXrADZsmWLNWbMGOuRRx6xXn31Veuee+6xzj33XOsXv/hFRc9Pp9OWJCudTrvcUvMkk0lLkpVMJh1vm6bU5zP9c6M6lfx3EY1GLUlWNBr1qplA4NXTfwcq6PzxH/+xdffddxfdd8UVV1grV66s6PkEndqEtdMPW7hDbeyfCzvQFP58lPsegMqFIuhkMhkrEolYTz75ZNH93/zmN63Zs2dX9BoEndr09vaW/AWdTCat3t7e5jaoieioUAmnURuCMtA49fTfgVl1dezYMeVyOV144YVF91944YV69913HZ+TyWSKtlAfHBx0tY2mCvOKkHg8nl9lRgF2dapdYh1UpWrX2BcH8IfAFSO3tLQU3bYsa8R9tjVr1qi9vT3/1dHR0Ywmok5+KgQOUwF2o41WrBuJRDxqWeOU2s14165dJffFSSaT7IsDNFPjB5jcUcvU1Ycffmil0+n815EjR5i6CgC/1AQx9VCslilMk6+hX35OgTAIRY2OZZ0pRl60aFHRfVdeeSXFyAbyuoOkExup1mtiap1TmGvXgGYLTdCxl5c/9thj1quvvmotW7bMOvfcc63Dhw9X9HyCTrB42UHSiTmrNYCyxBpAPUITdCzLsjZu3GhdfPHFVjQata655hprz549FT+XoBM8dJD+U20ANXVEB0Dz1NN/B64Y+Rvf+IYOHz6sTCajF198UbNnz/a6SXAJhcD+VM1RIKWKdfm3BNAsgQs6CIfCbfSdOkiOYfBOpQHU6Wwne9URYQdAswRmHx2ER2HI2bVrl1KpVNEeJLt3784v30VzDQ8v9m1p5L4wuVyu5BJr+/sA4DoXptJ8ixqdYCgsBB5e7Dp37lzqPDzCSjQAXgnFzsgIj8IpqcKRHHuHYk6A9gYjNACCqMWyLMvrRjTL4OCg2tvblU6n1dbW5nVzUIVYLJavCyk81gMAYL56+m+KkeF7rL4CANSKoGMQP50R1SgsTwYA1IOgYxDTDlFkeXLwmBi2AQQbxcgGKSzctW87hYWgoPjVG4lEQpFIxPHnJZVKKZfLlQwsdtiWipebF/4cAkBTNXwNmI+FZXk5W+6jHvUuI/f6QFYA5gnVWVf1CEvQsSzOiEJ96g0rhG0AjVRP/83ycgPZ0wT2KqUgTlvBe/X+HLElAIBGYXk58lilhEap5vDO4dgSAIBfEHQMwiolNFKtYYWwDcBPWHVlEFYpoVGqObyz3PMKH1/J8wGg0Qg6Bim3RwmdCypVT1ghbAPwG4qRgQL17CFjCq4BAL+pp/9mRAcowIZ3jAwCMAtBByhg2u7SABB2TF0BDtiLCAD8o57+m6ADlMCGdwDgD2wYCDQYG94BgBkIOsAwbHgHAOagGBkowIZ3AGAWgg5QgA3vAMAsFCMDAABfoxgZCIlEIlGyViiVSrFjMQAMQ9ABAsTeuXl42LFriyKRiEctAwB/okYHvsI5S+Wxc7P/8DML+BsjOvAVRixGF4/H80veY7EYIcdj/MwCPmeFSDqdtiRZ6XTa66agjGQyaUmyksmk422cEY1GLUlWNBr1uimhx88s4K56+m+CDnzJ7ijsztzLDqO3t7fk+yeTSau3t7e5DbL8dX1wBv8mgHvq6b+ZuoIvxePx/PEL0Wi07mmZelYr+W1qohk7N7O6q3qN/pkF0CAuBC/fYkQnOBr913GpqYRKpxj8MjVR7+fw2/uYhBEdwD1MXVWIoBMMboWKel/XDx1ZM6fR/BLugoBrBbiLoFMhgo7/uT2SUG9YCVsBsB/Cnd8x+gW4r57+myMg4CvN2JMkFovl6ygymUzFz7NrY+w6jLAs6a71eoUF++gA7qur/2547PIxRnRQ6whFWKcmGNEB4AdMXVWIoBNutYYVL6cmvFzaHtZwB8B/6um/OQICoeB0TILTcQpOcrmc4zSVfTuXy7nV7PzS9uHtK/w8bqjnegGAnxB0EAr1hJVy9RVud/ZenW3lZbgDgEaiGBkIgLAWQgOAVF//TdABAoLVTwDCqp7+myMggABIpVL5kJPNZht63AMAmIygA/hcM862AgBTUYwM+BirnwCgPgQd1I2dYd3D6icAqA9BB3Xzaq+XMPByaTsAmICgg7p5tdcLAACjYXk5Giase70wdQcA7mJ5OXwhHo/nQ040Gg1FyJF+N3U3fBWUHfwikYhHLQMABCbofOtb39J1112nsWPH6rzzzvO6OXAQ1r1e4vH4iCXfTN0BgD8EJuhks1ndfvvtWrRokddNgYOw7/VSGHZisRghBwB8InA1Ot/5zne0bNkyffDBB1U/lxodd5QavQjjqAbHNABA49XTfxu96iqTyRR1NoODgx62xlzs9XKG09RdWAIeAPiV0UFnzZo1WrVqldfNMB57vYwcvbJvS+G5BgDgR57W6CQSCbW0tJT9OnDgQM2vf//99yudTue/jhw50sDWA2eUOqYhbHVKAOBHno7oLFmyRHfccUfZx1xyySU1v34sFlMsFqv5+UAlmLoDAP/yNOhMmjRJkyZN8rIJQN2YuqsNGy0CaIbALC8fGBjQwYMHNTAwoFwup4MHD+rgwYM6efKk100DUAM2WgTQDIEpRu7p6dE///M/52//0R/9kSTp+eefV2dnp0etAlArzkgD0AyB20enHuyjA/hPWM9IA1C5evpvgg4Az7HRIoByONQTQGCF9Yw0AM1B0AHgmbCfkQbAfYEpRgZgllIbLUpiV2kADUPQAeAJNloE0AwUIwMAAF+jGBkAAMABQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMAAIxF0AEAAMYi6AAAAGMRdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjBSLoHD58WAsXLtS0adN0zjnn6NJLL1Vvb6+y2azXTQMAAD52ltcNqMRrr72moaEhbd68WZdddpl+8pOf6C//8i916tQpPfTQQ143DwAA+FSLZVmW142oxdq1a7Vp0ya99dZbFT9ncHBQ7e3tSqfTamtrc7F1AACgUerpvwMxouMknU5rwoQJZR+TyWSUyWSKniOduWAAACAY7H67prEZK4DefPNNq62tzXrkkUfKPq63t9eSxBdffPHFF198GfB16NChqjODp1NXiURCq1atKvuY/fv3a/r06fnbR48e1fXXX6/rr79ejz76aNnnDh/R+eCDD3TxxRdrYGBA7e3t9TU+xAYHB9XR0aEjR44wBVgnrmXjcC0bg+vYOFzLxkmn05o6dap+85vf6LzzzqvquZ5OXS1ZskR33HFH2cdccskl+f9/9OhRzZkzRzNmzNDDDz886uvHYjHFYrER97e3t/ND1wBtbW1cxwbhWjYO17IxuI6Nw7VsnNbW6heLexp0Jk2apEmTJlX02Lfffltz5szRtddeq8cff7ymDwsAAMIlEMXIR48eVWdnp6ZOnaqHHnpI7733Xv57H//4xz1sGQAA8LNABJ1nn31Wb775pt58801NmTKl6HvVlBjFYjH19vY6TmehclzHxuFaNg7XsjG4jo3DtWyceq5lYPfRAQAAGA2FLgAAwFgEHQAAYCyCDgAAMBZBBwAAGCu0Qae7u1tTp07V2WefrcmTJ+vOO+/U0aNHvW5W4Bw+fFgLFy7UtGnTdM455+jSSy9Vb2+vstms100LnG9961u67rrrNHbs2Kp3/gy7b3/725o2bZrOPvtsXXvttdq7d6/XTQqk/v5+zZ8/XxdddJFaWlr09NNPe92kQFqzZo0+85nPaPz48brgggt066236vXXX/e6WYGzadMmXX311fkNF2fMmKEf/vCHVb9OaIPOnDlz9G//9m96/fXX9cQTT+jQoUP6whe+4HWzAue1117T0NCQNm/erFdeeUUbNmzQP/3TP+mBBx7wummBk81mdfvtt2vRokVeNyVQ/vVf/1XLli3TX//1X+ull17SrFmzdMstt2hgYMDrpgXOqVOn9OlPf1r/+I//6HVTAm3Pnj1avHixXnjhBe3cuVMfffSRbrrpJp06dcrrpgXKlClT1NfXpwMHDujAgQOaO3euFixYoFdeeaWq12F5+f+3detW3XrrrcpkMhozZozXzQm0tWvXatOmTXrrrbe8bkogfec739GyZcv0wQcfeN2UQPiTP/kTXXPNNdq0aVP+viuvvFK33nqr1qxZ42HLgq2lpUVPPfWUbr31Vq+bEnjvvfeeLrjgAu3Zs0ezZ8/2ujmBNmHCBK1du1YLFy6s+DmhHdEp9P777+tf/uVfdN111xFyGiCdTmvChAleNwMhkM1m9eKLL+qmm24quv+mm27Sf/3Xf3nUKqBYOp2WJH4v1iGXy2nLli06deqUZsyYUdVzQx107rvvPp177rmaOHGiBgYG9Mwzz3jdpMA7dOiQ/uEf/kF33323101BCBw7dky5XE4XXnhh0f0XXnih3n33XY9aBfyOZVlasWKFZs6cqauuusrr5gTOyy+/rHHjxikWi+nuu+/WU089pU9+8pNVvYZRQSeRSKilpaXs14EDB/KP/6u/+iu99NJLevbZZxWJRPTnf/7nVR0pYbJqr6V05kyym2++Wbfffru+9rWvedRyf6nlOqJ6LS0tRbctyxpxH+CFJUuW6Mc//rG+973ved2UQLr88st18OBBvfDCC1q0aJHuuusuvfrqq1W9RiDOuqrUkiVLdMcdd5R9zCWXXJL///bp6X/wB3+gK6+8Uh0dHXrhhReqHhYzUbXX8ujRo5ozZ45mzJihhx9+2OXWBUe11xHVmTRpkiKRyIjRm1//+tcjRnmAZlu6dKm2bt2q/v7+Eec0ojLRaFSXXXaZJGn69Onav3+//v7v/16bN2+u+DWMCjp2cKmFPZKTyWQa2aTAquZavv3225ozZ46uvfZaPf7442ptNWqgsC71/ExidNFoVNdee6127typ2267LX//zp07tWDBAg9bhjCzLEtLly7VU089pd27d2vatGleN8kYlmVV3U8bFXQq9aMf/Ug/+tGPNHPmTJ1//vl666231NPTo0svvZTRnCodPXpUnZ2dmjp1qh566CG99957+e99/OMf97BlwTMwMKD3339fAwMDyuVyOnjwoCTpsssu07hx4zxunX+tWLFCd955p6ZPn54fURwYGKBOrAYnT57Um2++mb/985//XAcPHtSECRM0depUD1sWLIsXL9Z3v/tdPfPMMxo/fnx+xLG9vV3nnHOOx60LjgceeEC33HKLOjo6dOLECW3ZskW7d+/Wjh07qnshK4R+/OMfW3PmzLEmTJhgxWIx65JLLrHuvvtu65e//KXXTQucxx9/3JLk+IXq3HXXXY7X8fnnn/e6ab63ceNG6+KLL7ai0ah1zTXXWHv27PG6SYH0/PPPO/4M3nXXXV43LVBK/U58/PHHvW5aoHz1q1/N/3f9sY99zJo3b5717LPPVv067KMDAACMRTEFAAAwFkEHAAAYi6ADAACMRdABAADGIugAAABjEXQAAICxCDoAAMBYBB0AAGAsgg4AADAWQQcAABiLoAMg0Hbs2KGZM2fqvPPO08SJE/Wnf/qnOnTokNfNAuATBB0AgXbq1CmtWLFC+/fv13PPPafW1lbddtttGhoa8rppAHyAQz0BGOW9997TBRdcoJdffllXXXWV180B4DFGdAAE2qFDh/Rnf/Zn+sQnPqG2tjZNmzZNkjQwMOBxywD4wVleNwAA6jF//nx1dHTokUce0UUXXaShoSFdddVVymazXjcNgA8QdAAE1vHjx/XTn/5Umzdv1qxZsyRJ+/bt87hVAPyEoAMgsM4//3xNnDhRDz/8sCZPnqyBgQGtXLnS62YB8BFqdAAEVmtrq7Zs2aIXX3xRV111lZYvX661a9d63SwAPsKqKwAAYCxGdAAAgLEIOgAAwFgEHQAAYCyCDgAAMBZBBwAAGIugAwAAjEXQAQAAxiLoAAAAYxF0AACAsQg6AADAWAQdAABgLIIOAAAw1v8DplfmbiUeNWMAAAAASUVORK5CYII=",
"text/plain": [
"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\");\n",
"\n"
]
},
{
"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": 3,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAG2CAYAAAB20iz+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VGX2wPHvzCST3gkhIfTesQCCoAiKSijKqru2tf12XQu6q7trTQgZFazsurourr2tu3YgiChSRLqC9F4CJKT3MvX+/nhnSAIJJFOSSTif55knuZPJvTfJZO6Z9z3vOTpN0zSEEEIIIdohfWufgBBCCCGEr0igI4QQQoh2SwIdIYQQQrRbEugIIYQQot2SQEcIIYQQ7ZYEOkIIIYRotyTQEUIIIUS7JYGOEEIIIdotCXSEEEII0W5JoCOEEEKIdqvNBjpz5sxBp9Pxxz/+sbVPRQghhBB+qk0GOhs3buT1119n6NChrX0qQgghhPBjbS7Qqaio4Oabb+bf//43MTExrX06QgghhPBjAa19As113333kZKSwuWXX85TTz11xseazWbMZvPJbYfDQVFREXFxceh0Ol+fqhBCCCG8QNM0ysvLSUpKQq9v3hhNmwp0Pv74Y37++Wc2btzYpMfPmTOH2bNn+/ishBBCCNESjh49SnJycrO+p80EOkePHuXBBx9k6dKlBAcHN+l7HnvsMR566KGT26WlpXTt2pWjR48SGRnpq1MVQgghhBeVlZXRpUsXIiIimv29Ok3TNB+ck9d9+eWXXHvttRgMhpP32e12dDoder0es9lc72sNKSsrIyoqitLSUgl0hBBCiDbCk+t3mxnRmThxItu2bat33x133EH//v155JFHzhrkCCGEEOLc02YCnYiICAYPHlzvvrCwMOLi4k67XwghhBAC2uDyciGEEEKIpmozIzoNWbFiRWufghBCCCH8mIzoCCGEEKLdkkBHCCGEEO2WBDpCCCGEaLck0BFCCCFEuyWBjhBCCNGOpKenYzKZGvyayWQiPT29ZU+olUmgI4QQQrQjBoOBtLS004Idk8lEWlraOVdgt00vLxdCCCFEfampqQCkpaWd3HYFORkZGSe/fq5oM72uvEF6XQkhhDhXuIIbo9GIxWJp00GOJ9dvCXSEEEKIdiooKAiLxYLRaMRsNrf26bjNk+u35OgIIYQQ7ZDJZDoZ5FgslkYTlNs7CXSEEEKIdqZuTo7ZbCYjI6PBBOVzgSQjCyGEEO1IQ4nHDSUonysk0BFCCCHaEbvd3mDisWvbbre3xmm1GklGFkIIIYRfk2RkIYQQQogGSKAjhBBCiHZLAh0hhBBCtFsS6AghhBCi3ZJARwghhBDtlgQ6QgghhGi3JNARQgghRLslgY4QQggh2i0JdIQQQgjRbkmgI4QQQoh2SwIdIYQQQrRbEugIIYQQot2SQEcIIYQQ7ZYEOkIIIYRotyTQEUII0Sakp6djMpka/JrJZCI9Pb1lT0i0CRLoCCGEaBMMBgNpaWmnBTsmk4m0tDQMBkMrnZnwZwGtfQJCCCFEU6SmpgKQlpZ2ctsV5GRkZJz8uhB16TRN01r7JFpKWVkZUVFRlJaWEhkZ2dqnI4QQwg2u4MZoNGKxWCTIOQd4cv2WQEcIIUSbExQUhMViwWg0YjabW/t0hI95cv2WHB0hhBBtislkOhnkWCyWRhOUhQAJdIQQQrQhdXNyzGYzGRkZDSYoC+EiychCCCHahIYSjxtKUBaiLgl0hBCinUhPT8dgMDR4sTeZTNjt9jZda8ZutzeYeOzattvtrXFaws9JoCOEEO2Eq84M1B/ZqDsS0padKUiTkRzRGAl0hBCinZA6M0KcTpaXCyFEOyN1ZkR7I3V0mkgCHSHEuULqzIj2ROroCCGEOMkXdWakoaZoqyTQEUKIdsRXdWakoaZoqyQZWQgh2glf1pmRRGfRVkmgI4QQ7YSv68zUDXaeeuopSXQWbYIkIwshhGgWSXQWLU2SkYUQQrQIaagp2hoJdIQQQjSJNNQUbZHk6AghhDgraagp2ioJdIQQQpyVNNQUbZUkIwshhBDCr50TycivvfYaQ4cOJTIyksjISEaPHs3XX3/d2qclhBBCCD/WZgKd5ORk5s6dy6ZNm9i0aRMTJkxg+vTp7Nixo7VPTQghhBB+qk1PXcXGxvL8889z1113NenxMnUlhBBCtD2eXL/bZDKy3W7nk08+obKyktGjRzf6OLPZXK+YVVlZWUucnhBCCCH8RJuZugLYtm0b4eHhBAUF8Yc//IEvvviCgQMHNvr4OXPmEBUVdfLWpUuXFjxbIYQQQrS2NjV1ZbFYyMrKoqSkhM8++4w33niDlStXNhrsNDSi06VLF5m6EkIIIdoQT6au2lSgc6rLL7+cXr16MX/+/CY9XnJ0hBCtbbcZjttgYlhrn4nwhNmhPga1qXmRtuucy9Fx0TRNGsoJIfxeqR0WVcCnZbDdDAkGGN8dDLrWPjPhjjwb3JMD3QPhhQTQyd/Rr7WZQOfxxx/n6quvpkuXLpSXl/Pxxx+zYsUKlixZ0tqnJoQQp9E0+KkGPiyFpZVgcY6dBwLDg6HcAdGGVj1F4aZDFthhhq1m6BsEd8e09hmJM2kzgU5ubi633norOTk5REVFMXToUJYsWcIVV1zR2qcmhBAnVTlgQbkKcHZbau/vb4RfRcK0CIiVAKdNGxUKqfGQng8vFarAdVRIa5+VaEybztFpLsnREUL4ykGLCm6+KFejNQDBOhXY3BgJg4JkisOb0tPTMRgMDTYSNZlM2O120tPTfXoOf81Vf+8EA2R2hSgJYH3mnGgBIYQQ/sauwbJKuP04XJkF75WqIKdLADwWB6u7w9MdYXCwBDneZjAYSEtLw2Qy1bvf1WXdYPB91JEeDz0DIdcOfy/y+eGEm9rM1JUQQviLcjv8t0yN4Byzqft0wIQwuCkSxoaCXgIbn3KN5KSlpZ3cdgU5DXVZ94VQPcyKh9uy1XPhhkjoH+Tzw4pmkqkrIYRoohwbvFsCH5dCpfOVM0oP10fCTVHQJbB1z+9c5ApujEYjFoulxYKcuu7PgW8qYUQwfNhZRu984Zyto9NcEugIIdzxSw28XQLfVIBzAIfeRrgzGqaEQ4gkAbSqoKAgLBYLRqOxVUqOHLfCVVlQo8FrneDy8BY/hXZPcnSEEMLLNA1+qIRbjsN1xyDTGeSMDIbXEyGzixrJkSCndZlMppNBjsViOS1npyV0DoTbotTnb5W0+OHFWci/qBBC1GHTILMcph+FO3NgfbVKZrw2Ar7qAh8mw2VhkoPjD+rm5JjNZjIyMhpMUG4Jt0Sr58nGGthR0+KHF2cgychCCIEq6LegHP5VDEes6r5QHfw6Eu6IhsR2nn9jtkOJGYrNUGKBajvU2Go/uloegEq81ulUZefQAAgPVLeIAIgwQqcQCPbx1aWhxOOGEpRbSqcAuDocFlbAO6XwfHCLHVqchQQ6QohzWpUD/lcGb5bACWcCTowebo2GW6Igpp3URrHa4Vgl5FTDiSrIqVIf86qhyAyVtrPvozligiApFBJDoXMY9IuGflEQ4qWrjt1ubzDx2LVtt9u9c6BmuD1aBTqZ5ZAWDxEyZ+IXJBlZCHFOKnfA+yXwTgkUO0crOhrU6M2NURDWRi9SmgYFNbCvFA6Xw+EKOFIOx6vAcZZX+wCdClCig9RITYhBBSbBBgiqE/BpGmioOkIVVhUkVVjVrcQMFY0ETXqgRyQMiIZBsTCiA4S2o5EyTYNJWXDYCn9LgJSI1j6j9uOcbeophBDNVemAD0rhjWIocQY4XQPhd9EqD6etdaOussHuYthTCntLYW8JFFsafmxYACSFqVGWTiHQyfkxLlgFOGEB3lkaXW5VI0Y5lerj4QrYVQz5NXCgTN0WZUGgHi6Mh3GdYGRHFVy1ZTodTAqD10tUfzMJdPxDG39aCSFE01Q64KNSeKMEipyzGj0D4b5YmByuRjPaghIz7CyGHc7b/rLTR2oMOugeAT0i1Mdu4dAtAuJaqA1FRCC8OO/0Fg0F1bCrBJ5/xsSxCjudbkhnbS6szQWjHkYnwPU9oWcbHnCfFK4CnZWVKq+prQXO7ZEEOkKIds3sgI/KVJKxK8DpGggzY2BKhP8HODU22F4Mmwtgc6GajjpVQoiaDurrzIPpGVl/qqk1uFo0QG3eTIcQWPGCiYWvpDF7dga3XgyrTsAPOZBdBStz1G1MAtzUu20GPEOCIN4A+XbV3XyENPtsdRLoCCHaJasGn5XBK0WqFxGoAOeeGJgeAYF+HOBkV8K6PNiQp0ZvbKeM2HQLh0Ex6jYwVgU6/qapLRp6RMJv+6iRqc8PwaocWJOrbmMS4Obe6jFthV6nupl/W6kKTUqg0/okGVkI0a44NFhUoZosZjmXiScGwP0xcG2kfwY4Dg12l8C6XFifB0cr6389PhjO6wDnxcHwOIhqQ/2Umtui4Ug5/OeAGuXRUEvZp3dXwZCvl6x7y/xieKEQrgqDfyS29tm0D9ICookk0BGi/dI0WFkFLxXCLmcybpwB/hADN0b6X66EzQHbitTIxdpctcTbxaCDIbEqQffCeOgc2rb7J7nTouFIOXy4H1afUNvJYfDk+dC1DbRXWFcFt2ZDUgCs7N7aZ9M+yKorIcQ5bVM1vFgIm5wVacN08PsYuC3av5aJ2x2wpVDloazLU8uxXUIDYEQ8XJQAF3RQBfjag4ZaNDSlkF+3CHj8PNiYDy9vUzWA/rQGHhoKF3dqgRP3QF/niFu2TRKS/YEEOkKINmu/BZ4rgOVVajtIp4r83R3jP4X+NE0t+16eraZj6i79jjKqlUZjEmBYnFpu3Z6cmpPj2oamVy0eEQ//uBjmboGtRfD0ZrihJ/y2r/+24YjRq6raVZoKdnoYW/uMzm0S6Agh2pxCu0oy/k8p2AEDcF0k3B+rSvH7g4Ia+P44LD2mVhS5RAbCuES4pJNKJDb46cXaU95s0RAdBE+PgLf3qoTl/x1UrSoeGOKfvz+dTjX63GeBYx4GOunppy/TdzGZTNjtdtLT090/wDnAT14ShBDi7KodqpLx/GKodGYXTgyDv8ZBTz9411xjU/k23x1XU1SuBMggvZqSGp8E53dofyM3DfF2iwaDHv6vv6oNNG8bfHtcBUB39PPaKXtVUoAKdHKsZ3/smTS0TB/qB5LizCTQEUL4Pc25kuq5wtp+VIOC4NE4uCi0dc8N4FAZLD6qRnCq61y/B8XAFZ1hbGLbr/rbXGcaZfCk2ebEzmol1gtb4ZODaqn9hM5u785nIp3BbKWHy32aukxfNO4c+9cTQrQ122vgqQL4yZlo3DkAHoqDKeGtm6Nhtqucm6+Pqmq/Lp1C1IV3YmfVaqE5vD1N0V6nPSZ0hqwKNYX19+21TUP9Sagr0HGc+XFNUTfYeeqpp5q0TF/UOgcGUIUQbVGeDR7LhRnHVJATooM/xcKSrjAtovWCnBNV8NZu+O1yeGmbCnIMOhjbCZ4ZAW9cCrf0aX6QA7XTFCaTqd79rnfwBkPzMqy9vT9/8tu+MKojWB1g+rn+8nx/EO68ulZ5IdABFey4Vq4ZjUYJcppBRnSEEH7FosG7JfBqUe2w//QI+HNc6yUaOzT4uQAWHVHLnV2zER1D4OoucEUyxHqhiJ+3pyna87SHXgd/GQYPr4UjFfD6Lnh0eGufVS2jMxA3e6lSnbvL9AWgnUNKS0s1QCstLW3tUxFCNGBVpaZNOqxpvfep24wsTfu5qvXOp9qmaYsOa9pdKzTt6sW1t8fXa9qaE5pmc/jmuBkZGRqgGY1GDdAyMjL8an/+ZH+JpqU4/y5bC1r7bGrNzVfP4Tn5nu/L9fdz/d1O3T4XeHL9lsrIQohWd8wKzxSo/kCgKhr/OQ5mtNIUVakZFmapEZwy56qZsAC4vDOkdIXkFqjO60414Zbcnz95ZbtKBu8bBfNG+0cV6afz4Z1SVZn74Tj399PYCFx7GZlrKqmMLIRok8wOeKNEdRav0VQ9nFuj4IFYiGiF9JFjlapOy/fHweLMrUgIgWu7q+mpkBZ6xWxsmqKx5GJNg1mzTVRV27l7ZjplVVBWBRYrmK3w0ZtqfwGBan833GHihttTCTBAYABEhEJ0uPMWBuGh/luMryE394Hvs1Vhxh9zVb5Ua3OtKvf0KePtZfrnIgl0hBCtYlUlZBTAEecVYVQIzOoAfVqhYeXeErVUeU1ubf5Nnyj4VQ+4OEHVcGkpZ6omrKGSi4/lw9QbU8krgfwSWPAfE2sWpjFqcgaWL+rvb8MSE+sXq6+NvCqVDUtMfPJOGll5MPKqhkcCDHpIjoduCdC1o/oYH+O/wU9MkGr8+d8DkJnlH4FOuTP+8LQFia+W6Z9LJNARQrSobKtaLu6apoo3wKMdYGp4y045aBr8UqiWKG8prL1/ZDxc11PVwGnpKZC6QY7ZYufBP5uYdlMq1+5TCcUXTc5g1OQMXv9HGplff8+MB5bXC2TGTU0lKgwiwyAyFBb/V33t5t9l8NvfqRGcaWNSeS8Z3n89jZ5JcM3NqZRVQkkFlFZCWaXqyXUkV91cQoJgcHc4vw/06woBfrZg66ou8L8D6m+aUwmJYa17PoXOQCfOz35P5yIJdIQQLcKmwful8LdC1QPIgGq6eX8sRLTgiImmweZC+Gg/7CxW9+l1MD5RBTjdI1ruXOoqq4Tj+XZ+c1cGIX1TWfy+iXWZaazfpUZeRhWCw2EnOEhdOY/vX8G/Hg7CarXw18cyyJidypyn09EMBh5wvtM/vM5eb2TIVTfn8vmp9ElW0x43jK9/HnY7FJVDVp4KdLJy4WgeVJth4x51Cw2Cob3gwr7QJ9k/cmISQlTV6Z8KYOlxuK1v655PkQQ6fkMCHSGEz22tgdR82OnMgb0wGNLjoV8LTlNpGmzIh4/3w55SdV+gXo0EzOihLpQtyWyB/dmw5yjsyYKcIjD2SSceFWSMuDKV4EBY8WUa/brAW6+m8tGbJp5+So34uArHGY1Gnn1GBTantgtwTXs01C6gsWkPgwHio9XtAmewYLfD4ROweb+6lVfBup3q1isJrh0LXRN89qtqsiuTVaCz7Dj8tk/rBmAFzkAnVgKdVieBjhDCZ8rt8FIRfFiqcl8i9aov1fWRLZfvoWmwqQDe3wv7y9R9RmeAc31PiAtumfMAlU+z7RBsPwSHctQUkYsOSIyDnkkqeOiRCLEzUzGZ1LTVx2/XVsQFGkxW9lXdHIMBenVWtxnjVID2817YuBsOZMML/4ML+8GU0RDbSiNiACM6QoBONVTNrlIVk1tDtQPynYFOcmDrnIOoJYGOEMInllXCrDzIdb7gT4+Ax+IgzsNXHU1r3jv1pzerJGOAYINaHj6jh0pg9TVNg8O5sO2guuUW1/96bCT07wL9ukDfZAhrYFQpNTW13ugN0Giy8qnBji/aBej16lz7JsOVI2DRWjWdtWkP/LIfLjsPJo0AYytcXYIM0D8athfDtqLWC3SynAn2UXqIlhGdVid1dIQQXlVoA1MBZFao7a6BkBEPF3uh+eY3R2FVDiSFwY4VC0gs+uWsfZw+PwTv7YUp3eD6HhDl4wBH01ROy+b9sHmfyndx0euhdxIM6QkDu0GHqLMHba5AxjV6AzSppkpL1s3JyoMvV8P+42o7KQ5+NwXiWuFl9oN9Kv9qfCL8tZUqJX9TAfefgKFB8FmX1jmH9kbq6AghWp2mwZflqvBfiUMlG98ZrWriBHuQbKxpYHPAv3bBujyYlKz6TR3uMZXvNmzDPPdvPPXoH08+3nXRnzBhAqBGcC5JhA7BtV/3RTPLnEI1lfPzfigqq70/KBAGdVfBzYBuKpG3qU4NYC677DJWrFhx2uNOranS0u0CunaEmdfC1oPwv+WQXQgvfwb3X6tyfVrS4Bj10ZWH1RoOqHiU7jJt5R+8XKXZr0kLCCF8I9uiaXccr23dMPWIpm2r9t7+C6tVG4b1ubX3LcnStBmf5Wq9b0o7rTT+hAkTGiyR7+3S+WVVmrZii6Y9+x9Nm/ly7e3hf2raW4s1bfM+TTNb3dt3Y+d6tp+htdsFFJdr2lPvq9/Dk29q2onCFjnsSflVqh3ElK81zWZv2WO73Jut/g/eKGqd47dHnly/ZURHCOE2TYP/lcGcAtWAM0gHM2PVSE6gh8nGdkdtob4DZWC2Q3ydHJYru8Ch8o4ETL+Lfzz/B556KqhePsqpuSveKpnvcMCuLFizHXYcUdugznVQd7igHwzqBkYP3827UxG3oZ+xoQRlX4oOh5kz4NUv1SjXy1/Afdeo6ayWEBusVtNZHZBfA528MGXaXK7VhYNbofilOJ3k6Agh3HLcCk/kwY/Vavv8YJjTEXoaPdtvkRne3K2mIC5PVhet/Gq4fQU8MxKGxYHVDoEGyK6EF7fCknf/xu53HidAs9fLRzk1v8WTIKe4onZJdXGdvJuuHWFkfzi/L4S38BL1UzXWIgJ8N2XXmIpqFewcL4CwYPjzr1suZ+fuVXC0Ep4eAed1aJljupTYYcQh9flPPSBSkpG9wpPrtwQ6Qohm0TT4bxnMrTOK81CsKv5n8HAUJzML3toNwzvA5C4wKBaC9CphN30TVNjghYvUY+2aOt7v3vyBPbY4Nv7lYmrKS04LZjxJytU0Vedm1VbYcVhtg8qzGTkARg+CxFjPfub2rLJGBTvH8qF/V7hnWsvUtnl0PWwtgr8Mg8uSfH+8ur6vhLtzoEcgLO3WssduzyQZWQjRIk7Y1CjOqiq1fWEwPNMReng4igNQbIaV2fC7AarGzakmd4U5W2DpMZWQjAamp0x8/Nb/uOSVzRzIKebNl06frnInKddsgQ27VYBTd0l4ryQYMwiG9W6d5dNtTVgw3HYlPPsf2J0Fm/bCiH6+P66zeDSWVuh3ucE5wjmylUf3RC35VxVCnJWmwcIKyMiHUgcYdfBwLNwe7b3Cfz/kQLlVBTmHyuGrw6ruTfcI1aTxwnhV+Xb+ThgYA287g5p7X1tMcGwAscH181FWrFjB999/32jOTkOKy2HFFli7E2qcK2eCAmHUABg3BBJk9KbZEmJUvZ3MdfD5KhjQ1fdTfCcDHceZH+cL652BzigJdPyGBDpCiDMqsqvCf0ucTTiHBMFzCdDbC6M4AA5NBUsGPSSHwY8n4PVdoB3dgmYw8mPiQH4ugD8Ohrv6w+FymLm0iMqgEfxh/lIKe17BpBg1xQUqiDk1yHHdD7VB0LJly06ew9F8WL4ZXnvZhMNuZ9TkdOKj4ZKhaooqpJk/qz/lyviDiefDz/tUcvKSDXDdpb49ntH5XDC38IhOqb02EVlGdPxHC7bSE0K0NT9UwpQsFeQEoGri/DfZsyDnSLkq4Jd5RHWZduW9VNngWKWqlTOtG/Q+spAPbhlE3M5Pya+Gzw+rxOTg7+ay9eMX6JjcnYihV/C7AXBn//q5H+PGjWt0xdKECROcQZCJXUfglS/g+Y/hn39XTTQ7RBu4eyo8cQtcOqz5QQ7U9pwymUz17neNKhkM9TNU09PTT3us6z6TyVQvKLLbIWO2iVmz0mkrAgyqHxbA+l21o2W+Ync+pwJa+Aq3ugocqP+PBBlG8BvypxBCnKbGAc8Vqm7jAL0C4YUEGOxhX6j/HYAP96su08uOw6IsGNcJbuoD07vBF4fg++Pw2Hnwq9RUdED6n27mlrd6UxAxnNlPPcPs1CfIyMjgiRv7NzptdqbRkm+/W8bMh0zMmpXG4vWqM/jGb0ysX5zGQ49k8OJcz5dfN7fnlF6vAqOyIvi/m1IpK4IjOwy886n6/hkTMnjjCbBUQ+YaE4vXpzH5ogz+WQ56AxgCICwKImMhItb5MQ4SukJEjMc/jlf06wIdoyGvRLWLGDvEd8eqtqmPIS284ukHZ+7apa2wpF00TgIdIUQ9O8zw8Ak44OzXc2uUasTpSXVjUI0W1+WpKajLOqsRnE8PwtdHISEUJnaGG3vDv3ZCTpWqo+MKCD7bvovdR/PZ4AxyUlNTae56UYcDthyApRvB0DOVUZNh/eI0fvr2KWxW7/aDgvrBTt2eU08+kUphNuQfg4JsKDgOSTWpTB4FL/wtjZ1r4aqRqeQdrd1XTSVUl8OSDc4gZ1QGV41IxW4Duw2sZvWYwmz1+MXr09HrDFw1MpXIOOjcG7r2h64D4PkXW2fqTKdTwc3nP8DqbXDxYN+twKpxTlkFt+AVTtPUiA54p92J8B5ZXi6EAFSuzFsl8FIhWIF4AzzbEcZ5qTHi7hJ4eC38cyx0c3a4LqqB/xxQicgfTFBTDY+sh1KL6ix+SaJq93DDa6vIWvIGBav/2+wl4g5N9Zz6er0aTQCVYHzJULjuMt/3g3Itbw8MNPLpy2ZyDqmRmVPpdPD9VhNfrEgjMMCI1WZh5u8zCDDCvFdqawGlPZnBY4+mYreDw6amsmxWqCyBsiJ1m/+Oife/SCPlogyuHFEbvH2z0UTmujQe+0sGzzzn28KBDakyQ+pbYLXBozdCko9q3Ny/Gg6Ww+wLVEfzlrDTDNOPqnILm3p4/sZA1OfJ9Vv+FEII8mxwZzY86wxyrgiDRV29F+SAGqHpHgG5dS7yscFwdRcID4Q3dqv7/jIUekXCK9vhiQ1w7woz1YXZFK3/8uQS8aZwaLBlPzz7Ebz7jQpyQoPg6pEw+3bYsuz0pefeoGlQlAtblsNNU9UxDHojVquF+e+YsFRDoBGSesGQcXDZb+CGh+EPz8Pny1MxGlWQYzQaeXl+Ki/9I/XkORqNRmabUgkOg7BINU0VHQ8dkqDbQBgyFi6eBu99nkpGRgaZ69I4rJk4fyKs3KmCnMmjMkg0p7L4Tcg/evafpzkayjVyMZlMPDcnnZ6JavtAtnePXVd+jfrYwcOp1uZY6mxie0moBDn+RqauhDjHLa+ER/PU6qpgHTzZAW6I9P60QvcItQpmRzEMiYUQ56tPUpgq6rY+T01vdQiBmYNhajd47YPPWf3q8zx002RSK8qatERc02DnYVi0TlXlBZVQfNl5cOlw9fmp+TJN2e+ZOOxwbD8c2AJHdkF5Uf1ppmsnprJyh4kPvkpjyMXwzLOp6BvIH2mo7g/gVi2gulNnr7yups4eeSiDcQNSObgVDvynJ+kYAAAgAElEQVSibt0GwthrITah2T/2aVxJ2HWP7/q5XL/vXkmqCOOBbBg31PNjnqraBj+/nY5Ob6Dj5S236u0b56rEK8O9ulvhDV7tuuXnpKmnELVq7JqWkVe/Eed+c9O+d9asWac1iay2adqeYk17Ys7ftFmzZtX7ms2hPn52UNOuW6ppWwvq72/5cU37v5WaVlhTe587TS0PZmva3z6tba7559c0bdFaTav0cL8NsVo07eA2Tfv2A017/VFNe3lm7S1ltNrXH27P0AqyNc3hOPsxGmvG2dB9zWnQaTQaNUAzGo0n7ys4rmlL3tG0fzygzvfVhzRt98Ym7/KMztZUdO+x2oafrt+LNx0p07Q+t3jnb9xU+83qf2jAPk0rtXl118JJmnoKIZrlsAUezK2t+XF7FPy5g8ovaIpT37lnHoF394G9NI/yPnfSs2gjOVWQGKqW+rpG8mf0UKutPj8MYYHQ0znVrtepJox1NaepZU4RLFoD25w9hgINcMkwuPx8CDulnok7zTJdNAccPwB7NsL+LWCpqf1aSDj0HAo9hkBepJ1RVzf9GE1tONrcBp2NVYaOS4Irb4NRV8PKTyFrNyx9H6wWGDzmjLs8q8aSsF33d3PmzJRWQrUZQr08vXSsEnrflEpscNNXvXnqK2fvszGh0tvKH0kyshDnmAXlkJan+lTF6FXxv/Fu5OK4LhwPz3uX0vN/S8i2L3jD9BC3ZcxHN3ASXcPhkWFqCkzTVM6MQQ8HSuHFbRAeoJp2dghWTTwHxsAfBjavX1Zppaq4u36XOoZOp6oYXz0KYrw4hVCUC7s3wN5NUF6nJURYFPQaBr2GqpybhqajmqKhAoOu+4DTplqaMv3S2PTcqRd6zaGCnW2r1fa4GTB8vHs/R11n6jH26L+hqgYeuwkSvdzV/IN98NF+1Sak/EvvNXVtjEODy45Atg3+lgApEV7dvXA6J5p6zpkzh88//5zdu3cTEhLCmDFjePbZZ+nXr+mNUyTQEeeyageYCuCTMrU9Mhhe7ASdPBjXNZlMvLfHTtIlv2HV3QNOXkiWHlM1cVw1cuwONWrjyvvZUQxfZ8GBMqiwwsiOKi+nqSw2WLEZlv4EFucy+KE9Ycpo6OSlNg02q8q52b4Gsg/U3m8Mgd7Dof+FKrjR+WHiaWNBTaPBjgY/fgWbv1fbF02BEZM8P35jAcacj1SV5Humq5YQ3pTxkypjcPcAmN7ds6auTbG2Cn6bDRF6WNu9tkK38K5zoqnnypUrue+++xgxYgQ2m40nnniCSZMmsXPnTsLCvLg0RIh2aL8FHjgB+yygA+6PgftimzZ6cqZ2BhqgNwRSmXOAkMiYk48ZkwDZlfB9NlyaBJ2d/6LVNpWEPChG3UrMqpJsbFDTfg5NU60EFqxRfakAuiXAteM4uZrHU0W5sONHNYJT46yLotNDtwHQfyT0GAwBgd45lq80d3pOp4OLp0NgEGz4GtYtUqu6Bl7U/GM3JdE7MlQFOuWVbvxwZ3HAGcj3jGx86s6bPnc+DyeHS5Djt7ycL9Ri8vLyNEBbuXJlk79HkpHFueirMk0bul8lS44+qGlrKpv3/WdL3u1y9e+1iR/lauHdBtV7zPYiTfvLWk17f6/a3pCrabe+u1l75Pn5jR7n1CTmuo7katqL/6tNNE59S9M27NY0uxcSWh12TTu0Q9O+fLV+UvFbaZq2frGmlRd7foy2Ym2m+tn//bimmaua971NTfT+1wL1N/xxu7fOWjlRpWlXL9a0lK81LTX9zEnR3lBg07SBzv+tzdVe261owDmZjFxaqmrTx8Y2Pk5tNpvrDVWWlZX5/LyE8BdmBzxdAP9xPu1Hh8BLCdChmf/1Z2pnAPC70cns69SR62a/TdqdF518zKAYCA2AYue/YKUNyiO7886/ZhNWndvolMqpKqvVUvE129UIkjEQrrgALhuuPq+rtLQUq9Xa5J/NaoZ9P6vpqVLnUnR0qopw/5HQpS/o9VBjg5qCM+6q3ehxAfz8A+TnwbefqoTlpiovL+fRRx/lnnvuoaCg9hd2zz33UFVVRXl5OQUFBZSXQnUFlJVAgRd/rz/mgKUUSjNfxPTB3Hrn4jqHtLQ0qqqqePjhhz0+3rslUF0EA4KgczkUVJz+mMDAQKKiojw+lnBfm8nRqUvTNKZPn05xcTE//PBDo49LT09n9uzZp90vOTqivTtqhZknVDsHHXBvDMxs4lRVY07NuwBOTk8sz4aXtkKfPR8x70+3kDF7NqmpqTy7RQU6c0epfZyogn+/2LQkWYcD1u6EhWtV4irABX1h+sUQ3UCicWlpKfPnzz/5JuhMrFbIz4LcLJWLA2AwQHxn6NgVgs7xEv4l+bBvs5rSGnKx938fP+2FwlIY3BOSvJiMvCFPVUQ27F9FYpiesWPHnvaY1atX43A4uOSSSzw6lkODD0uhXIMJodC/kenXqKgo7r77bgl2PHROJCPXdd9995GZmcnq1atJTk5u9HENjeh06dJFAh3Rrq2ohIdzocyhVlW9mOC9CseuxE6DwcCsWbPqBSbPbFadyYO2fUl0wXbu+fOTzP1FVT6+qkv9/ZwtWfVoHvx3OWTlqe2kOLjuUujdufFzKygoYN68eYSEhBAa2vCVuaYSju+HE0dU+wSA4DCVVJzQFQI86MrenmiaGuUqyYXEnir52pt++AXyy2BEP+jipRYNmgYfH4Byq3q+dfVx4b49NfBZOYTo4P5YCGwgP6eqqorq6mr+9Kc/0aFD/X4X69evZ+7cufz000/k5uYSHR1Nz549GTNmDC+++KJvT74NOieSkV1mzpzJggULWLVq1RmDHFAvykFBTcxyFKKNs2vwShG84lz+PDwIXu4EiV5KnD01sfNUDw+FV3fAur7XYDjvGh5cA/1jYGyn0/eVmpp6ssaK0Wg8GeTUWNRy8VVb1YUr2AgpF6lmkIYmJnqGhoYSHl7/KlddAYd3wonDajl1oB4iElSTy/hkNT3lDxYuXIheryclJeW0r2VmZuJwOJg6dWqLnEvfIfBLMVQUQGiod39HDgMYgyA2BsK9FJDkV0N1AIQEQp+ODQce3vSLFfRhMCr0zKUMqqtPb2yWmZnJtGnTGD9+PM899xyJiYnk5OSwadMmPv74Y78NdM60MMFXFae9oc0EOpqmMXPmTL744gtWrFhBjx49WvuUhPAbJXY1irPKuUrolih4rAMYvdTGoSkraYIMcO8g1brhWCVEG+G8Rpo2nho0ZWSYmH5TKp+uUrVxAM7vAzPGQaQHo1HVlXBkJ+Qc5GS38+iOagVVbCffdc92l16vZ8GCBQD1gp3MzEwWLFjAtGnTWuxcYhLU6jJLDZQWQky8d/araapQIKh2HN5yyLn6qUu474OcE1Y4aFWFMEeGnPXhp3nuuefo0aMH33zzDQEBtZfh3/zmNzz33HPeO1Eva0qLD3/UZgKd++67j48++oivvvqKiIgITpw4Aaj5z5AQN55pQrQTO81wXw4cs6leVaZ4uMaLM7N1X8TufCiVpUcbr9AbbIA+UerWlP2lpqbyeKqJWbPSWLweRl6VSocouH68Z/VVapwBTvYhNYIDKrDpMQiivHTBbiqHHezVYKsEazXYq8BW5SxwqAd06qNOB+P6p2Crpl6wUzfIaWikx1f0BujQWY2C5R/1XqBjtYHNubo9xItVkQ86k+57tEDBvnXOQZoBQRDtRpHIwsJCOnToUC/IcdH7y/BiA860MMEXxRi9pc0EOq+99hoA48ePr3f/22+/ze23397yJySEH1hYDo/nQY0GXQLg1UT14utNdrud2bMzGH1HKjN/BIsdOoc3rWXCqeq+KD7xZCo/bgdbciqjJsP6xWn07gwv/DMVo5uvTDYrHNwGRUdVMjOokYkeg1WXb1/SNLCWqdVZNYVgLnJ+LEYVC2qiXqQwKkEFO5mLFmN32JiS0rJBjktckgp0SvO9t88KZ5AQFKhadXhDiRnyatQIS3cfBzpldtjsTI4f7eZ77NGjR/PGG2/wwAMPcPPNN3P++ecTGOjnxZmcztbiwx+1yWRkd0llZNFe2DR4oRDeLFHbl4SqpeNRPuizU2WFV3bAihy1PTQW/jxMtW5oLtcc/x9mpvKf72H/cXV/twTI3WwiNMi9OX6bBVYuKODFF+cRZIgjyBhOdLzqOxXjYbJrY3kzmgO++jQTc6mD0clTqcpWozYN0kNAKASEQECY+qgzqH1oDkADza6WRptLQLPCqzvvw6HZ0OsCmDnsVSJ6QcwACOvcctWYK0pgwxI1hXXJr7yzzwPZsHE3JMSojvLesC4XNuRDt3BVDdmXllTA6iroFgj/F33m6c+KigoKCwtPS0YuLCzkmmuuYfVq1XcjMDCQESNGMHXqVO6///7Tcsz8ka8rTp/qnEpGFuJcV2yHP56ANc53xn+IgT96uHS8MftLYc4WyKlSLRxu6Q3X93L/WGlp6azYAnM/AqsdjAGQMhouHQr6G5r/jlBzwK4NsC4Tck+AzQaxsTBwBMQleicHp27ezFUTUyjPgvJDsGxdJmuzFzAqfhplrnxTAwTHQnAHCIqF4Dj1MTC86cGJ5oCFX2Ti2GHDoA/A7rCx7ngmI60plO6GwAjocAHEDvZ9wBPivN7arKrmUKAXRgtLnLk0MV4aedE02OMM+PtHe2efjal0wAbn3/qSUPefX3Fxcfzwww9s2rSJZcuWsWnTJlasWMFjjz3G/Pnz2bhx42mrtPxJS1Sc9iYJdIRoQ/aZ4ffOfJwQHTybAFf74M2fpsGiLPj3LjV6FB8MjwxXjTfdlVcCH34Hh5wjQ/26wG8ugzg3y4scPwA/fAb5x9R2eDT0HAy9B4E3B2wnXZpCZY6aSspdByPjU9iQl8n6/AVclDiNy85LITQRQjtDaEfQezgDsfjrTDKX1ubkLFqUycKFCwhJgGHGFKzlkLMCyg5A1xQw+HBJvCEAgoLBXKNWrnkj0Cl2FtVrqBaSO3KroNSqEpB7+njaam0VWDRIDIC+Xvi9X3jhhVx44YUAWK1WHnnkEebNm8dzzz3nt0nJTVmY4G8k0BGijVhWCQ+fUF3HkwPgtcTGi5R5otIKf98Oq1W+P6MT4I9DIMLNC7hDg9Vb4as1KhE1KFD1pho90L13xGWFqgHl/i1q2xgMI66EpIFw4hXvLIO2m6H0AJTugcpjMFBLoSIe1uctYFP+YuyajavGT2P6r1O8uuy6ocTjKVNS0OlUoBUxBUZ3TSF3DVQehcNfQLdpairMVwKMKtCx2zzfl80ORc6k4VgvBaM7nKM5PSO8l/PTkAo7rHWO5lzqwWhOYwIDA5k1axbz5s1j+/bt3t25lzSUeNzYwgR/IoGOEH5O0+D1EnixULVAGBWi6uPE+uBF/UAZzNkM2VUQoIM7+8P0bu6/qBeVw0ffwV7nqEvfZLhponsXOasFfvoWfl6mLro6HQwaAxdNhpAIz1sJaA4V1BRtV1NTWp0c69AkmHZJCj+9vhib3UZAQADX3uj95GCHw9Hg6irXtsPhIG4YhHSCIwugOhcOfgbdrwGjj9I63E3ibCi3Kb9UBb7H92SyosLBtGme1QSy2GGvsxD2YC91rW/Myiowa9A5AAZ5+AYjJyeHxMTTu9Du2rULgKSkJM8O4CPNbRbrLyTQEcKPmR1qVdUC53D/jZGQGg+BPsjHWXoU/rkTLA7oGAKPDYd+HuQ8bNwDn6xQRQCNATDtYlX4T+/GuR/aDqs+hbIitZ3cB8bNUMufPWWtgpJdULxdJQO7BMVCdD+I6gvGKDXa4gpybDYbmZmZXl8JdaZigHWPFZoAPX4Fh78ESxEc+gS6XwtBvshRcUU6zRy5aqgmUG4RHN6RyaHtC+jf1fOaQHtKweqAWCMk+bBtR6m9NjfninDPR3OuvPJKkpOTmTp1Kv3798fhcLBlyxZefPFFwsPDefDBBz0/aR8400IBfxzJcZFARwg/VWCDe3JgixkMqADnZh+0y7HY4bWd8I1z1GVEPPx5KES4mYNQZVYBzk971Xb3TnDLFdDRjYtwWZEKcA45R/LDY2DctdBrmOcXm6pcKNyspqhwvhHVGyG6P8QMUgnFrmOcOqXk2gZaZdk3qKTnntc5g50SOPat2vb2lIqrDlFzA1TX76Xu72nFMhXkjLvM8+XymgbbnIHvoFjfFn9cUameIj0CoZcXVoE/+eSTfPXVV8ybN4+cnBzMZjOJiYlcfvnlPPbYYwwYMMDzg4iTJNARwg/tNcPdzqTjSD280glG++Ada241PP0z7C9Tb9hv6QM39HJv1AXgwHF4/1s1ZaXXwZUjYdKFTW/f4GK3w5blsOFrteJHr4fhE2DklZ4lxGoOKD8MBZuh6njt/SEJEDMEovqA4ZQLWUN5Mw1dxFuDMRK6z4B970J1DlRlq+Xn3qI5oMY5kmF0o5xA3d9TZuZi7HYbPQZP47oZnv++jlZCQY1KQh7gw9VWBTb4yVk3Z2KYdwKqG264gRtuuMHzHYkmkUBHCD+zqhIecCYddwuEfydCDx+srNlcAM9ugTIrRAbCX4fD+W6uaLU7YMkGWLpJvdPuEAW3ToIeDfS5Opu8LFj2HyhwBiKde8P46yH29JSGJnPY1fRUwWawOHuBYYDovhA3HELOUEywKXkzrckYDtED1NRb/ibvBjpmswp2dDr3O5inpKSwePFibDYbOn0Ao8alEOyF5/NmZ07WwGgI9uGVbFmlqvfYzwjdpelrmySBjhB+5MNSyMhXL6wjglWl4xgvJx1rGnx2CN7Zo47TJwoePw8S3Fy5U1wO734DB53LxkcNgF9dQrMvZlYLrF+sRnI0DYJDYey10H+k+++iHXYo2akCAKuzfos+GGIHQdwwVd/mbJqaN9OaOlwAxTuh4ghU5585cGsOszM3LNiDpp6ZmZnYbDb0+gAcDhtHd2XCcM9+b4U1cKRCjUIO92G5mWwrbDODDrjCg55ronVJoCOEH3Bo8FydSsczIsDU0XtNOV3Mdvj7ttoqx5OS4d6BYHQzmNp+CD74Dqpq1LLx30yAC/o2fz/H9sGyj9TScYC+F6hk41A366I47OrCX7ARrM6LdUCYCghiBvq29kxrCIqCyJ5Qtl/dvBXoVDqXgge7uaLLNe036appWKNSOLwjkx9XLCAu0rMgcaOzJUXPSIjy0d9S0yDT+dwZGgSd2kaHBtEACXSEaGVmB/w1DxY7X1T/GAv3xng/uTK/GkzOfByDDu4eACld3TuOzQ4L18ByZy2bLh3h9qsgvpnJ0lYLrFkAW1ep7bAouOzXqjeVOzQHlOyF47/UjuAEhDsrCQ8CfTt+xQvrooKc6jzv7bPUGXhGurF0u25uU+d+KezOgjGXpjCkp2e5TUU1sM+5Om6ED/uXbTPDEata4XiF/3dkEGfQjv/thfB/xXa1suqnGggE5iTAdB9Ud91RrJKOSywqH+eJ82BInHv7KiqHd5bAYWdBwfHDYeqY5hdryzkI330IJc5354PGwNjpYHRjCk3TIP9nOPYNBFVBSIAKcOIvUCuo2nOA4+IaxanOd3ZG90KgXOrMg4lyY3rIldt01dUpLPhR3dcrCToP9Sy3aUO+WvHeKxLifVQo0eJQPa1AtXpwp0O58B/nwL+/EP7pqBXuyoZDVojQw6s+Wln17TH4x3bVyqFnBKRe4H4+zq4seO8bqKyBkCC4eSIM7dW8fditsO5r2LxMXZDDomDiTdDNzRW1pXth95tw6CdVByckAhJGQdxQz9sxtCXBHVThwOB47wQ5VjNUOaeuIt0Iil25TVm5YLaq50uic2TI3WmruqM5I304mrOyCsocEKOHsT6sOi1ahgQ6QrSCXWYV5OTbISkA3kiEPl5u5+DQ4N298MlBtT22Ezw0xL0VKg4Nvt0Ei9epd9NdOsKdVzW/T1VRDnzzXu2Kqv4j4ZIZ7q3oMRfD3nfh2FK1rQ9Qq496jXFvBKKt0wdAeFfv7a/QmccVFune0nKXPc76TD0TPW/PsS5PPf96RvhuNCffprqTg+ojF9hCneKF70igI0QLW1cF95yACodasvpmEiR4+T+xxgbPb4W1uWr7xl5wcx/36uNUmeH9pbDjsNoeM0itqgpsxjlrGmxfDT98qUZ0gsNgwo3Qa2jzz8dhU+0P9n8ENucFKWkCxKbALx+AwXlRbqgFAYDeamDh0gU4HI4zrqg617mCUU+qTxeUQmGpet719nDZe161yi/TARcleLavxmgaZJar4oD9jDDQgwBP+A8JdIRoQV9XqMacVtTy8X8lQqSX5/8LamD2T6pvVYAO/jQELnPzIpNdCG9kqgtWgAFuGA8XDWzePqrLVV0cV3Xjrv3h8pvVlFVzFfwMO+erhpagCvwNuAdi+p/e68rVgmD37t3079+flJQUei+/COPWGBzFHTF3L2XJwi+wYmlXAU9jAR6oBOGmBHgOe+2IjieBzh7n36lbJzV15QlX0N43Cjr4KADZZYb9VlWJfLIkILcbEugI0UI+LIXZzkTKSWHwUgIEeXlY/GAZzNoEhWa17PbJ82FQjHv7+uWAqnJssUJsBNw1WU1ZNcexvbD0PbVMWW+Ai6fBsEtB18yf21wKu1+H7OVq2xgFfe+A5Msb31fdqrzH92Vz67G/oLOEMu9oKpMGTuHygusp2l/GwQt/at7J+LmGekxB/VVQZ1OUqxqnBgW7t+IKoKIajjkTzfsmu7cPl6MVtXVzLmrmc7CpzI7alY9jQyFOro7thvwphfAxTYPXimGesy/PzVGQ2kEt8famn/Lhmc1QbYcuYTD7QujkRu6LQ4NvNsDXG9R232S1dDy8GTkRDgds+gY2LFE/f0wCXHk7xDdzdEDTIPt72PVvsJYBOug2Ffrc0rRifymT1YX+l693klzUhz/m/oaEKyJZzPt8vfMr7uv+OOOixrKfdWfcT9VasOWDVgWOSufHKiAADNHqpo8CQxQEdoeAMyTKemPE5Yw/cwPtKRpqY3EmOc68rviuzQ9KXXYeVn+/TrEQ48FKQocDVjlHl4bEQpSXc9lcllVCiQOi9WqllWg/JNARwoc0DeYWwlvOQoD3xcCDPmhAuOQovLJDBSlDY9VITrgbK47MFlUA8JcDavvSYXDN2Ob1qqoqV6M4R/eo7QGj4NLrmt+jqioHtr+iGm8CRPSAwQ+ojuJno7Pr0QwO0KkLffVyPVVaBUX2fH5avBqbzcb0qdPRRdYQfTiZ6KxESrrmNLo/83awHmn4a45iNRVZl7EvhI6BwF6n/629MeJyGg2VvOJUN9hxtV9oapBjqa7Nz0nq2fxTATWac8hZfmBwD/f24bKtWI1QhhhglI9Gc45aYa2zp9e0CO+PtIrWJYGOED5i0+DJPPjMWbju8Q5wh5ebD2oavLcP/usMTCYkwYND3FspUlQOry9UeTkGPfz6subn4xw/AN+8raaqAgJh/A0q0GkOzQFZmbDnLbCbVUfx3jdBjxlnr4ejtxsYuOYSwi3RVMUWU9Ari/9ueo89NTtIDOhKR0Mi+yw7MAYYmTxlMvknDhF5Ip4O+7ufMdAJHgGBPUAfBvpQ9VEXCpoFHKVgLwV7CdgLwHoILHvVLaAjhE2CoP61+/LGiEs9Dj3oHWA3gMFe7ziuICcgIKDJ+z1xWD2vIuMg3M3na93RnA5u5GK5VNtgvbMA4kUdfdPTyqbBV+UqVhweDH19NGIkWo8EOkL4gFmDh07A0kqV2PhMR5gR6d1j2Bzw8nb4zvnu+8Zeqvu4O6NFR3Lh9UVQXgURofB/k6FHM5poahpsWQE/fqUClZgEuPpOiGtmI87qfNg2DwqdFZdjh6pRnLCks3+vdZWRGxb9lcq4Yiq6FJC4vR/lW638uHUdF6dcxMEj27kr5s/8WP0dFpuFzMxMNdoTVU54QQwGcyD2oFPHZpSQ4U07/4ULF0KEnvEhKdT8DLY8KP0AwifD8uLaaSl3R1xOnfYKO9yf0GO90NkD2Ve6kXztIH1SwtACbCd7TAUEBGCz2U7+vGficMDx/erzRDdHYsoqa0dzBnk4mrMmF2rsKvnY3Vyzs1lRCSdsEKpTy8lF+yOBjhBeVuWAe3Pgx2pV7fjvnbxfQr7aBnO2wKZ8tXR35iC4sot7+9qyXyUdW22QFAe/n6qSj5vKaoHv/wN7nTm9/S6E8b8GYzPeGbtycXb+C2yVoA+CfndAtylNyxHRzGD+PJhdfVZybNw2wsPD2Zi1hmGHL+WqK6+inBLm7n2ENzp/zaOXmZi7PPXkiEpCj1503DMWe6Ct6SfcgIULF7Jv3z727NmDfhpc/ZcUKr+D6vXw8kcvcci6p960lDsjLnWnvW7smUbY0T5UdzrCwcMHseYGMajjlURsDGJlwQcsWFQ7QuQaMXIdtzH5x6C6EgKN0Kmbe7+HLQfU3zOpQ/NbgtR1olJV9AYY74UaPA05boVVzhIF0yIgTKas2iUJdITwolI7/C4HNteod4ivJcIYLyc2lphh1k+qQmyQHh47D0a6kbugabDsZ1iwRm0P7KaSjpvTdby0EBa/oXI69HrVbXzoJc0bVbKUw/aXIdfZJiCqHwx7GMKcK3XS09MxGAykpqae9r0mkwm73c5jN6fjyDJQ1bn85NcqdeUUJ+Rg0ZlPTgvl6fZz09EHKZ6Ux/LD37BgwQKmjryD+K7Z6FDTF+7S6/Xs2bOHfv361QYVU1N4bbUKcnqG9WPyVfVzcpo74uL6+tLMZdw26X2KR6/ms03zWZC5gBum3krPIXdxYksQyYXjeHDKHAamxNb7vjMFO5oGWbvV5537gMGNHK/cYsguUH//4c2smF2XQ6ttPDsgGpJ80DncrsEX5eAAhgTBYKmZ025JoCOElxTa4PZs2G2BKD28kaTm/L0ptwqe2AjZVapnVfoF0N+NIX27Az5ZAWt2qO1LhsK145qXdHx0Dyx5G2qqVNuFq++Azr2bdx5FO+CX56AmH3QG6H0z9ACb+FkAACAASURBVLxeLUV3MRgMpKWlAZCamopmAZPpKVauWc7y75czO2M2+nDQJ9npffh8PsneyLDqPtxW/DsqOhbS9dCzXDp6GoEpFezSVjLivRncUzKbmwbdy6Gog1xQMp5DIzeg6T0Jc+oHE65gZ9GiRTgcDnoG9eOO0IeoWgthYzktJ6epIy6ur8doXckvOc7rr/yNn/O+ZdrUaUycMpZK9rDuly2MiruViYEzyatcgD2svN5+G+sxVZIH5UXqd5/cp/k/v0ODX5zTXr2SINKD4GR7EeTVqED+Yh8VB1xdVTtllSJTVu2aBDpCeEGeDW7Lhv0WiDfAO0neT2rMKldBTqFZ9aoyXQjJbrxAm62qKeeOw+qd94xxanVVc/yyCn74XOXjdOwKKXdBeDMCLs0OBz6B/R+ofYQmwfC/QlTf0x/rGslJS0tjwPYxjLJMpOSEme/Xfc+ECRN48slUdDoI+m01m574geFZ4xkQNZStv1pCVUwJHfb14NKfp3F4y2ayh+9i2/RviTuUTPSxRPrGDWTblCWUdyo4/cBuOHXkxOFwoNfrmXnTQ1QsAssuWFF6euJxU0Zc6hqTMoSydwoYEz+DncWrSZmSAg4d6DXGThsObEP7diBBhZ2oCqsd5Wpsv5oGB7epzxN7uNfy4VCOSmgPMHi20qrconJzQFVADvVBv7JsK3xfqT6fHA7h0rSzXZNARwgPZVvht9lwxAoJBvigM3RvxvRPU+wpgbRNUG6FruHw9AiIc+NiVF4F8xdCVp5q4XD7lTCkGUuIHXZY9RlsW622+42ACb9RK6yaylwEv7xQm3CcdBkMug8CzjDF9/i9qVy0ehq2tQHkBRQRVZHINZf+ii+//4ynM57iyVlP8vIPc5m7ZS5LB+0ib8IeCnuqsrzHz9tB9PFOROZ2IFuDmuhyjp+3i+whe9AC3OugfSYpKSknR3JABTvfHc7kIlKwnQB7gqPBxOOzjbic5NCR+fUiqo+t476B/8ShOfhm0TKunDKx3sNs4SUYS+OoYt9Zz7kwR3Uq1xugWzNX2oEqKrnVWXtncI+mT3+emlytabA8W3UPr96cyZGDDoZN827VaqsDPilTbR4GBsEwmbJq9yTQEcIDR6xw23E4boPkAHi3M3T18jvQzQVg+lmtPukXpQoBRroRSOWXwGsLVDuHsGD4/ZTmrawyV8HXbzvr4+hgzFQ4f2Lz8nGKtsGWuaohpyEIBt4LnS+vv4+GcnIcFhg+dhjX/jiJBH0SD8Sl89Bf7uVvE02kpqViesaExWLh+jE3EeOIoTB298nv1Wl6QkoiKUnO4WQijg6fBDkAL7300smRHIfDQb9+/chcvoDqcLiMFFLGTsXQSLXhpiQkZ3696OSIUOn5Kxm//gbMFdewc8ECel0RAQ49hppQggo6U9Z3y1n3pzngkHM0J7kPBLuRU7b9kKrBFBnavCrIp9YU2l0ChyugYEMmuesXMNidmkJn8U2laqYboYfpEd6vaSX8jwQ6QrjpiAVuPg65dugRCO8mQaKXg5z1ufD0ZlXrY3gcpJ4PIW781/4/e+cdHlWV/vHPzCSTSiqBJPTeqyBFRIogEAHbqriWta26igpYcJcJIcEVK7afva1t14IKShEBJYDSQXpvgZCeTDJJpp/fH2eGhDBJ5t4MTe/neebJnZl7z5x7M3Pv957zvu83K0+KHEslxEfBAxOgiYKpJpcL5r0KhdkyI2fU7coMOYWAo/Nhz3vywhrZCvo85dttu2ZMDkBQPDyw7GbWlP+EwWBgYqPbcL9o5ckPTaSlpeGwOzAajXQY1pqw7yNJ2t4Jgz2YH7Z+w6i8mxBGQUkLT3TrWbywvfTSS6cCkqdOnXoq9qZTp06s2LsA3HBNbkqtQqc+ChbEs3jxj0wcfw3jxo3FqjuKY1QR+cvDucR+C+JHByFhRoTBRXmrPdjjc+ttMy8LyorBECR9yJRSVAb7PSUO+nZUlh1VfcrO7oL8DinkrV9I/jqVNYXqYa8N1noKA16nZVn9adCEjoaGCo454NZsKXLaG+HjZEgI8K9p9Ul49neZHTK4KTzZC4JVxBLsPyELAdoc0DwB7h+vPFDUYJDZVOsXw/j7IEHBXbvLCjteq/KpSh4O3SZDUC1TBtVjcrzPrxw7khWZKxg5fCTLVizjrYf+Q/DCTjzzl5fRuXW4ceO2u1m9ejXdhi6n66YhhO6O428VXbFFW9hxwzIq48x19rOh1gwLFy48TeTA6RfythGdEMItCyupIGpXP2LcLRh79VWMSxl76nVXRBlxE8r4/Ie7aeMaTO+uPXFEFeFoVPf+gnSSP+AZ9GnVRXlsjtsNG/ZIIduiiSwQqJSUlBSEgO+/X4BOvwjhdnL1+MCLnFIXzCuVy4PCoINWGPBPgyZ0NDQUcswBt52QGRttg+GTZGgc4F/Szyfgxe0yk2VYEkzrqSwjysv2Q/DhEnC6oH0zuPdqCPNMewlRNWwv3PXXq+k+GDr2VXYxrMiFzRlQdki23/leaDXB93RB5SYIagpBiaeLnVmzZuFyuRgxYgTLli8D4P7X7+CzVcvomN2XW4fcQ7vRzUhNTSUzM5Pw8HDi4+axZ+tBLr28H71v8C/opKHWDG533bE3Zb+4GREyHr2K2jI6RzBR+3tRMHApKYkpGCrDCc1tSVBFI1yh5dhi87jq6lEAVOApk10jicyXkDuyG2yVsO/EQsrC3LTupiweZk8WFJeBMUiO5qil9eAUdAulyDEYghh/dWBFjvCkklcISAqC0VqW1Z8KTehoaCjgqEfknHTK6apPmgVe5CzJgtd2yOvUlc2kpYMaA9D1u+Hz5VIsdW8Dd46RAchenA55AQgO9t+4UYnIKdoBm2dLM05jNPR+CuJ9THdZt0LudHBbQBck7RJi/yHFzuzZs7Hb7RgMBpYvXw7IjK3Z/87gnV0fMK/dWsY3uZFrHh1B1NGmfPX5VyxZsgSACRMm0DvF/8jahloz1DXaM/bKFAp+k8sGFUInbstQrE2PY03MIsgSTdzmKzBUygrIQucm7GRrzF02YI/Lr9qoHo+tSgtk7YGdRxay4/ACWnRWFg9TWi5jcwD6dKgS0EoptcPX8xci3E70hiBcLv9qCilhXSXst0OQDm6IgmAtLudPhSZ0NDT85IQDbveInLYekdMkwL+gxcfgNU9tm5SW8EBXWflYKau2yzo5AJd2hkkjTx8Ryj0Kv34PlRZwWOW0VPNOyt3Fa+PEctj+CggnRLWHviYI8+Ho7SqGwrnSADPmbnAcg8LnoeBZ+C7hXex2O6HGUKx2KxkZGZhMJnQGcDnd/D31LtpFNqXpoqac+Cuk7LuXLaE7WVO5rNYpqPpoiBlmXbg8+kMXBnqF00OGigjCsttS3koGWEfv7oczwkxR71U4o4oxFiQSu+0yoncNIH/QIjD4DrKuvm9CQItGKWw/JEXOeIVTRW43rNsj/ybGQetEZftUvZ23vljIybULaDdkAo/dmsLiRf7XFPKHbAcs8aSSXxUBTbWr3p8O7V+uoeEH3jo52Z6RnE+bBT4mZ0lWlciZ2Ar+3kVdRsjPW+HbVXL5il6yEGB1sWQpge/egC6Xyho4xbkyXfzEAWnA2U5hTZ3qCAEHPoMDn8vnTS+TVY4NtVzcHcehcgPET4XgZvKh+ydseGw/IT8nkZ6eLkd2Zs3GlCqns0wmE2mzZuIshKK50lU8YjR80vo5/vPaqxgMBlwul+pRASXWDP7G9dg8SWBGFfVldK4gLG13YCxJIOnHSejcBvIuX4AzUgac2BvnUNZhK1H7+qBzG6Rrex37BlLs6HWLcAsnY0ZP4GqFU0V7sqDQLEcI+3dWn7n0zlcLObh6AUkDJ3DfTSnodcprCtVFpRv+Z5bB/J2NMCBMdVMaFzGa0NHQqIdClxQ5Rx1VKeSBFjk/ZkmDTmiYyFm6AX5YK5ev7AvjB5/Zzp710KQ5DL2+6rXmHWHTT7BtFaBTllHlxeWAHS9XBR23/Qt0vKPuaTFXkefiXy1A94XFGfy++RAz2r1E+/5XAzBj5gyEXpwKUJ7xlIm8x8G2E5r9F15YlEFqaipDhw5l4sSJrFy5UvWFUok1g79xPXav0OmiqCsAOBuZMXfdQEhhEmHZbXAH23AbbfJNT6q8K7QSoQO9IwRXsG9jUi+jRqbww/dS5Bj0QVx7vbLjU2KBHZ6aOX07yFIFajhRDicsbhIGTGDStSlEV5v68rumUB1443KK3BCrh+uj1I2Oalz8aEJHQ6MOzC64y1PxuKkBPm4mgxkDydLjVSJngkqRIwQsWgc/bpDPxw6AMf19t2MIhgqLDEI1hgA6aNERwiJlteMdq6WVg5J6Ko5y2DQLindIK4duD0KLMfVvF9IJHFlg3w+h3eVrLpeLwff3JdEWi2UxRI4Bdxk83NZE3MNJ5LuOowuCJk9DUJL0u0pNTWX69OnoPbnNakcFlFoz+BPX4yoGZw6ggxClAbuegGIR5MTaNAt7dAEGW1iV0PH8fyOOdsIZYcYVbqm3yc8+WohbONHrg3C5lcXDuNywdpeM+0purH7KqtIJPx6HhAHj6RIDnWLOXKeh01ZrKmCXTWrom6MhTEsl/9OiCR0NjVqwuOHubHmyjPeInBYBrpOz4gS8sl1ez8a3gvtUjuQsXl8lciYMhisvqX3diGjpaVScA4ltZLVjdNA4WVY5/u+zsHst9Bnh32fbimGDSWZWBYVDn39B4z7+bRuUCBGjwPwRRIyQgbppaWkAFL4I1i0eZ/IDUPQa3HTNPTSe7tnWU+zQ5XKRnp7OAw88wNy5c0+1rXRUwFfgsT+Cqb64nsr1cr3gNqBX6v/k/S4IHegE7tBK3KGVp97W20IIy2lFWE5LckbMq7e5Lz5byNotC+jeZgK335PC6nXK4mF+PyhHdIzB0L+Tuu+qW8BPx8HigBgjXKGgaKW/HLHDUq/FQyNodhZsJDQuHjSho6Hhg0o3/D0bfrdBjF56V7UNsK1D5kl4aZsUOSkt4X6VImfJevkAuGYIjKhHZHTsC4e2wffvwDUPyQBkt1ummEc3hhadwFJ/CRZApo9v+BdUZIMxFvqnQ5RC1+r4R+DYOCj9GqJvBb2nvomhCdj2yeWwPtD4XxAx7MztvcKooOBMvyplAbbqrRlqi+txV1YJnfDBfncF3HrCj7fDFp8jTTl13qEdTsumCrLEEH68PaWdttQ7mjP/u4WsyJQiZ9zYFKIbKxv5OlEA+6SrBgO6QJjKOjRbC2T14yAdjG0BxgD7TFlc8GWpdCXvFQKXahYPf3o0oaOhUQObgH+chA1WiNTDB8nQOcDFxX7Nhed+lyfj0c1ldpUakfPjBjllBTDxMt8ipzgPKspk1lVCczl1NXg8rHLCvJdh9O3QtkfV+k47Z9Rg8UXZUdgwA2yFENYU+j8NEcnK9yG4JcRNhpJ3QB8GjSYCBrBugkZXA5678UaBtTw6g7rSw+sTTLXF9VSulSNSQU3AqGDaKnbbYCKOdKGixX4qkg9jj8vBHWI7zcICwB6fS0H/ZQijvc72hIC8LDfd20xgYJ8U2nQ/c9/qEnLlVli3Wy53bAHNGvu/L9XJKa8y7ByaBAkBDg52C+ljVeqW5roTNIsHDTSho6FxGi4BU3JgdSWE6+D9JOgR4DvCTfkwZ4s8KQ9Phsnd1QVJLt8MCz2Bx+MHw8i+Z66zdwNs/EnWzCktkkX/ht8EUfFS7GxtBAvfg879ZYyO3QonD8PIW+r+bPN+KXIcZdLOof9sCI1Xvg9eYu8DZx4UvSFHdtxlMjU98WX/a/ycL2qL6xF2GLBLiojwYf7vh6EigpD8ZCztdmAsTiB2+yAqE49SmXwEe0weIsgFgM5uRO80ypGcGiM9NTm+H9oljEefCF0HSvPO6tQl5NyeuBy7A2IbQS+FI3ZeKp2w+LgU9x2joZsCCxJ/WVEOBx1g1Mm4nJAL/LujcW7QhI6GhgchYGY+/FQuT5RvJUHfAN9x7iqG2R7vqiGJMFVlMcBV22H+Grl89UAY5SMm58R++PkLGHYTJLWR01Hz/0/G4vS4HGKbwpBroU0P2LJcxu2ER8H1D0Okj+BQL+Z9sH4GOC0Q3Qn6pYOxkfJ9qEnjJyDqeplJBRD1l4a3ebapL66nPBxGtU0hpHtdrZxOkCUGZ2Qplta7cPYwE3moG5GHuhFSmERF8iGsicdwRBcRdaAXwaWxFAxcWqfIKS2qsnlo10vGaClh2yFpCBtkgMHd1FXodrtl+YQyT1zOiOTAj7Rst8IvFXJ5QqRWL0ejCu2roKHh4aUi+KJUXjNebAqDVLg418XBUpi5EWwu6NsYHu+l7qKxfk9VMcDR/WB0/zPXsVXIkZw+I+RojTf+puMlkH+iar1gI7TuKrOuDEHSvNNQR8xEdZET0wX6ZUBwjePky33cS0ZGBi6X61RcTXV0IRDSVT4CRUP9q+qjtrieUd1SsCwDgZvICcpGpewx+ZR22IorTEbTWtrupLzFPqL29iXySFdCC5KxxeXSaF9vCgYtrrMthx12rpH//4Tm0p1cCcfzYc8xuXxpF2jk52+i5nFfnQtZ5WDUQ9juhfy4t2HHvSbZDvimTC4PCYfeWr0cjWpoA3saGsA7xfBWsVyelQBjAuyFc9wCMzZAuRO6xsKMPhCs4tf3+0H4XNo9MbQnpAz0vZ61Ul7kmrSQz70X2rhEyM+So1duOQOC2yVFDtTtPF1d5MR2hf4+RA5UuY9nZGSc9ro3DdxQTUmlpaWdsV719X0JIiV469wsXLjwtNe9IzF6JVbbPhg/fvwZIkc4oWw+DI9IYdwV4zG2UtamMNqxx+cggpxySkroEMEOzN3XkX/ZQtxGG1F7+1KZdARr0+O1tyNkzaTKcgiLgM6XKhtFKas8PS6nZRP/t61+3PcWw9ZC+Xr47oUsW9Tw414dixs+M4NDQAcjjFaa2abxh0cb0dH40zO/DJ73nIgfj4dJKryI6qLICjM2gtkO7aNg1iUQ6vnlKRn92H8cPloiY3sGdIHrhtZ+4YqOlwUBvULH7ZYiJjJWLut0st5NpQWO7pJTGsEhtbdXeljG5HhFTr90mUruC1/u416RMystnamDq/bVK4qqb+fd99TUVNLT031/iJ801L9KDeU/S8sHfQREjlbZiK76X3EqONwZaaaszU7CsltT0uPXOpvI2gP5x6XI7TpYjt75i8MJq7fLv42jobfCuJzqxz0xD+L7p2DcuZA1KwJ73J1CVj42u6GxAW7UigJq+EATOhp/an6rgKc8WSB3x8DfAxwgWeGA1I2QVwnJ4ZDRDyKq1fTw90KflQ/v/CALtvVsC5NG1H5C94oar8gR7qqRmpAw+T7IAOUvX4Dk9vJuvzbKj3sCjy0Q07lukeOlutjxGnPOmpXOfZUmTt4PcQ/LTKu6RJHX/qGhnC3/Kl/YD0NFplyOHA/6QE1/ev7XOpeBmB0DqUw6giu8vNbVi07CwW1yuX1vKXzBv6m8lKvHs3YXmC0QapRxOWoGYIaNTmFzARz/dQF5GxbhdgX2uAsBP5TBEQeE6OQNilYUUMMX2tdC40/LPhs8mAMOYFwkPNGArCFfONwy8PhQGcQaIaM/RNdIUzeZTKSnp5821VPzQp9fAm/NB5sD2jeDO67yfeEpOCFtGGq+Vz0+xBgCtnIZw/P92zL4eNStte9DZT6s/xfYi6FRW/9ETvV9MxqN2O12jEYj97tMlH0N6MHY3vcxCAkJCajI8ZKSknIq9bs+/yq1uCtlxhgCQvtWVXpWRS3p/cLgorjvSgr7r6h108oy2PGbFAJJbU+Py/FnKm/nYVkzR6+HIT0gXEXWocsNi7Mg+pIUdIYg3K7AH/d1lbDRKjXgjVFa8LFG7WhCR+NPSY4T7jkJZW7oFwrPNQnskLdbyGKAWwsh1ABp/SCpjqme2i70peXwxnxYNi+NPZkZ3JsijRSrk5GRwbNPrGThu7Bvi0wRrw1jqJyiWvQBmPPhL1NrX9deKkdyrPkQ0UymkAcriF3KyMg4JXLsdjv/fj0DdND0OYgce+YxqC6KAilywHedm0Ai3FA2D9xmMMRDpNLruVt++XQOz3CfDlkN2YfgcUQVVxUQrIHTAdvXyFpIUfHQ6ZLTpyNTUlKYMGHCaWKn+lRej/4p7Dwi1+3fSU5bKUUIWHkSsiugaMNChCvwx/2AHRZ56iOOjoBOAa5zpfHHQtPAGn86ytxwbzacdELbYHgjKfD1Nv6zT57sDTr4Vx/oUM8Fw2QynZri8V7o7Q45XVVYChHhBn76OpUXep45xfXVhxu5e8LDGONh8zJwWKFDHwjzkfKtM0BJAbiccNMT8jVfcUIum/SuevfnDPRhLt76KI2QOlLOa1J9VOqhpiZmPZXBK0WpRF4Jsyf6jkeqLooyMjIUiZ20tDRsNpvPINeXXnqJvXv3+u1fpYbyZUiHcgNE3VRV3dkfwrPaEZbdlqCKRjgjzDiiirC03VXlZwWE5CcRZImhosV+GaRcg++//x6dTk+LqBQsJVLQdr9M1supmV1W21TewCEprPCkoXdsAW1UWjNsLoAdxZC/fiF56/z3DfOXHAf81yzr8fQOlVlWGhp1oY3oaPypcAiYfBL22GXw4vvJEBvgEvRLsuArj7vzoz3gkoT6t6l5oU9Pz+DjpXAsT7pDf/G+7ymu9FlPc9uN0+g5KJrrHpYu5JtXwI5foaz4zM+Jawq9r4Dx98tMHDgzS0q4YduL8Mb3GXy6P5UWVxoIU5BxU13kPNrTROHz8FCciaeuSefpL2vPxkpPT8dms52xn/5gMBiYM2cOq1evPu11r8jp1KnTaXVuao5qNITKTVVxOVHXQrCC6tChuS2I2XYZrrByrE2zEAYXoTmtaJI5gchDVXn24VkdiNrXB0Ol75QivV7P998vIPPXhegN0GOINGWtLbus5lTesBEprNouywskxSsPPvayzwxrciGvhsjxfmZDj3upCz4xy+rlrYNholb5WMMPtBEdjT8NQoApD9ZUQpgO3k2C5gE2+9uUD697Ct79tT2MbFb/NjVjcrzPB4yDQeNM3JMCCTG+A3zTZ6Vz1z1DsVlloPEVN0B4I9i+WrqTdxskCwMCFOXI9PLLrzv982sGBN+QaOKlj6TImf6PdGa/qGwayWuy+dhVJk78Tb4Wcxf8+ykTYRny/dr23Vd//BnZMZlMVFRUMGfOHMLDw7n22mtlarNH5EydevocnVLDz9qwH4YyOUhB+DAI7a1s++jdl2Bpu5PSLpsAOXVlLEkgLLsNEUc6Y7BGYO66AXO39VgLsnA28m1C1rNDCofawI7DC2jSAqIb151dVnMq752PF5LUIYWYSPXBxzkVsMyT7Z4Y6maQSt+w2nC44XNPhlWCAW6JhmBN5Gj4gU4I4YerzR+D0tJSoqOjMZvNREVFne/uaJxj3iiCuUVyGPOtJBge4HobR8pg2m9Q6ZICZ2qP+u82fV3oM7fBE09lsG5RKvc/nM6br5x+oQ8JCTk1+mOz2Xw1y85fYd1iaNUFel0BlhJY+TWk3AONaxFf3r4E6Yw4hZ3H7kzn+Q/Uxco4jkPWteAugYgxkPiK76J5aosL+qKgoIDrr7+ezMzMUxfws5VdBeDIhZL3QFRCSA9ZyVlJYUCd3UjjDVdibZJFWYftp72nt4YReaQzkUe6UtRnZZ31cgpOSGErBGRXLmTV2gV17n91ATRmbApvfLCQHRsW0KH3BB64M4UIFcHHZjt8dRAqXNCmkTSpDXTM2xelsNMmb1Luj4X4i+A23WKxUFhYyJQpU2jcWKVBWIAJ5G/uXNKQ67c2daXxp2CZRYocgNSEwIscsx3SN0mR0zMOHu7u35C6d/TDe9LZmwXfZMKlY0zcfl86TWNdp63vK5alOsJzo9xtMAy/GU4cgMx58MO7sipybSIH4KHrTATppcgJDjKqFjluG+RMliInpCc0fb52AZCWllbriI3JZFJ8wh0yZAgGg+GsZlcBOAvA/JEUOUEtoNF1yj25hNGOI9JM+PH2BFmiTgs8dodWUtp5C/aoIsJOtqm1jZJ8OU3pzbD669/qzi6rLnLGjkth/W5IaJtCux4T2L91Ab8sVz6lVOGA+UekyEkIhauaB1bkCCEDj3fawIBMI78YRM6FipKCnn8ULiqhk5mZyfjx40lOTkan0/Hdd9+d7y5pXATsssE0T62c26LhrwEuCOhwwezNkFMJJ79II2hxhs+qx74q/Va/0OeXwIeL5d1r/07w0ZunX+j9iWXR6eWFAaBNNxg8EU4chJ6Xw6Cra98HayFMnZSB020n2GDE4TxTRPlLwWyw7QB9LCS+BvoAm6LWxerVq3G5XGctuwrAVQglH0jj0aCmEHMb6FVOgZa33IvOFUzMtsEEl8adkWVlj8s9QwR5sZhhW6asbB2fDJ36waJFdWeXVbes2HYQjuZKQX7bJBk/o3RKye6C749BiR2igmFCKzAG+DqZWQFrK+Xy9VHQVkHhQ40z8aekxR+Ni0oXl5eX06tXL+68806uv/76890djYuAk074ezZUCBgcBk8FePRYCHhtJ+wshvAgGNHcwNOzUgkxKKv0W2mXGVYVNmjVFG4ecfqIUK2xLOLMWBbvduVm+HUBtOwsqyTXhtsBD0/I4KPtqdzdP523Vpp45oUMRTEyXkq/g9L/IdPIX1AWmNtQXnzxRTIzMxkzZsypGJ1AZ1e5iqH4Q3CXgiEBou9sWFFAR2wBhf1/Im7TcJqsmoClzS4qkw/jMlrRO4OJONYRS5tdZ5h2Wsvh919kOnl0Y+lKv3ixbxd1qNp/b+bV3qzTPayS4pUfI2+tnNxKCDPAxFanF8MMBJsqpckuyFpXPc+haP4j4zPe7w8qciAAQscb4qM7B6HvY8eOZezYsfWvqKEBWN3w4EnIdUF7I7yWGPjgxe+OwLITcmj0qd5wySgTMSHKKv0KAZ/9BLnFEB0BVRwW8wAAIABJREFU9/ioleNyuXjmn98xtNdElv8XWnSCZu3AlGo69X5NnA7pVD7u7rr34dEbMnj311Ru75rOK9+ZCApTFxDsyIJ8uTqxD0LEULl8LmICMjIymDNnDkOHDmXUqFGAb/uHhuAsgJIPPbVyEmSAtSEAnmiO6GJyR3xD5KFuRO3tQ/jx9niHcOyx+ZR12Hba+vZK2PqLDDaPiJKjdUt+rNtFvfrzwzmwZb9sq2dbaJOovM9uIb/3Ry3Ss218K4gNsAjZa5P2LACXh8NgLY08oPgqafFHRbXQef/995k7dy7798tfTIcOHXj00Ue55557Ata5hmKz2U4L1iwtLT2PvdE4lwgBqfmw3QYxengnCaICPKS+vRDe3yuX7+1SlUau9G5p2SbYdkg6md+TIsVOTVIGp7F1hQwqLsqB3KOwI1xOR5lSTadic6oT3bh+kZO7Dgp3ubi1QzrPvW8iotoIjLe/vkRUTYSA/JkyZiW0P8Q9VPXe2faz8vZx+vTp6PV6fvjhB0JCQhg1ahRXXHEFdrsdq9WKxWLhp59+wu12c9VVVylq35ErR6rc5bIgYPSNUKkDLA3u+iksTdaR02QdjfLaIxA4QyzYI4pwWxyn1rHbYMcqKC+FkAhod4msmG21WmnXrh12ux2LpapT3v3fuXMnVquV7pdcxfo98v/VLglaxINF4T4IAb/lylFMPXBFc4h0K2+nLk44ZIaVU0D3UBgsAtv+uaKiouJ8d6FWGlq76mJCldAxmUzMnTuXyZMnM2jQIAB+++03pkyZwpEjR5g9e3ZAO6mWZ555hlmzZp3vbmicBz4zw7dl8kT8ciK0CPCQeoEVntkq72yHJ8vYhOr4e7e0+xj8sFYu33CFnLaqSVEu7FwDo2+HVp7SKgd/h70b4adPYcg10LanfH3Xb1BRBn1GVDmS10ZlAWx/CW7tkEbra6CpDyd0f098loVQsQp0Rmjyb1mYsGYbZ9PPKi0tDbPZzNtvv43L5WLJkiVUVFQwZMgQ+vbtC8D8+fPJzMxk6NChFBYWnto2MzMTvV7PkCFDzmh39erVOEvd9MoZClYwNIbIcTImhcIzVq+TzMxM7EV6Lh825Izii6tXr8btdsu+GTwNO4Fq2eROh/yfV5RJK492XaHcKh/9+vXDarWett9eKioqOHjwILGNm7FqUyFCQHJjaBYDRUXK9gFgVzFs82w3sAlE2KFQ4bGoixIXfFsKVqBlEAzSQZE9cO2fa6KjowkODvAJqIHUVtIClE1VXyyoEjpvvvkm7777LpMmTTr12oQJE+jZsyeTJ0++YITOU089dVr9jNLSUlq0aHEee6RxLthQCU8XyOXH4+GyAA95O9zwzBZ5sWvTCCb7yLDy526psBT+86O8Qx7UFS6rxRvJaZOVjMOrBVG36wWNYmHbKtiwFIxhkNQajuwCcwF0HgCRdQRdCxf8/hw4yiCqPXS8U92xAHCVQsHTcjn2H2BsfeY65yImIDo6mvvuu4+77rqLF198kTlz5jB48GCmTZt2Kn5n+vTpTJs27dQ2zz33HDqdjszMzFPrernuuutYtWoV98VP529RUwi9FJq+DAYfFaf9oWCHm3dWz6GtGMyLP0zD6MmQrd63KVOm+NzWboXFH0BMOwiNgKvvraqP5GXKlCm17vd9D05H32oaDid0aw2TRsoRRKUsPwHrdkMb4PYOMK5VvZsoIs8pq5aHOqFPiKxaHnFRpcycSXBwMNHRAc6AaACBql11USFUEBMTI/bt23fG63v37hXR0dFqmlQMIL799ltF25jNZgEIs9l8lnqlcb7JcQgx6JAQ7fcL8chJIdzuwH/GGzuFGLtIiBuWCnHCcub76enpAhDp6ek+nwshhNMpxAtfCDH5VSGe/58Qdkftn1dWLMSn/xZiW6YQbtfp7508IsQXLwixZoF8bq0QIv9E/ftw6FshFo0V4sfrhLAcr3/9uih4QYj97YU4MloIt7XudY1GowCE0Wisc72ZM2eedryqk56eLmbOnFlvv7zH3fuZvtrzrjNixIjT1hk+XD4fGDZC7G8vxIl7hXBV1PuRPnG7hdjzgTzet3ao/7tRE1uFEF++KMSrk4V4+0kh8ur5f9Xc7ymPp4vH3pTftf/7Tgi7U91+LD8uxLhF8rv/4R51bdRFgVOIMUflb3fkESEK6vhNaKgnEL+t80FDrt+qhM5DDz0kpkyZcsbr06ZNE//4xz/UNKkYTeho1MTqFuL6Y/JEOe6oEOWu+rdRyooT8kQ/dpEQa3PPfL+2C1fN179bLS88T7wtRGE9X0eXU4gf/yPEx+lC5B498/0NPwrx3r+EKC+tvY3qJ7fybCF+vEZeeI8uatjJzVEgxIEeUuhYltW9rj/Co+a69R3H+vBHWNUUOwaD4TSRkzdTCLfKi67bJcTON+WxXjRWiINfKzsO1hoix9f/3xfV9/uJt+R37ZV5Qljt6vYjM1uIFM/3/vUdgb+BKHUJMdHz2x1yWIhjKvup8celIddvv6euqk8B6XQ63nvvPZYuXcrAgXJif+3atWRlZXH77berGFfyD4vFwoEDB049P3z4MFu3biUuLo6WLVuetc/VuDCpns0jBMzKh99tEKWHzh9k8ByBrfB5tAxe3SGXb24HA3z4P9UsAOilemDv7mOwfLN8/ZYREFdPkU+9QcbnfPUSLP4QRt0KTVtVxeBEN5aPumJyvAHBQsBVFSZcNojrCR9tyGDmTPUBwSXveCoD94TwEbWvpzQmIBBxPf4GW9b8LJfLhQEDn7RYTuMZsk6OGtxO2D4Xsn8GdNDtQWg5DkzX+xe/ZauA+W/KwPPQcLjmQUjwY+a95n6vXJDBpLtM3Hc1hKgIFVmfB8/9Lk00RzeHB7oG1l/K6ob7s2VBwDgD/Cc58DF1Gn9y/FVEw4YN8+sxfPhwxWrLX37++WeBzLs87XHHHXf4tb02ovPHovrd/X9L5N1gx/1C3D1D2V2/P1Q4hLh3pbyjfWqdEE6Vd7TmciH++a68w/7iZ2Xb2m1CLHhLiDcfk6M4h3cIkXtMiM+eEWLpJ/Vv7z1et3ZIF0uuEWLG4w07To48IQ5084zmrKz/c9WMzigZ/airbX8+y6A3nHZemXGX+u+Po1KIDalyFGdxihDHVyjbp4oyIf77rBzJeWd6/dNVNdue9mS6ePIdIQaMk89TZ6rbl835QkxYIr/3c7ao/97Xhs0txN9PyN9u74NC7Khn6lPjz8s5n7q6WNGEzh8P74k94dF00X6/EOOfCrzIEUKIudvkyf7W5UKUqDwZu91CvLVAipxnPhPCpnI6ZO0iIT6fIwXPZ89I8eMPtlIh7ujuucgGKRMOvih6S4qcY9fVPZVRMyag+vOa02a+ptH8jeup3oYSYeV2CXF5RzltpUeKnWGDRqg+PtYiIdY8LEXOkolC5K6rvQ+++lRWLMQns6XIefcpIfIVipypT6SLJ96W37MXvxTCNFPdb2JTvhATPSJn1kYhHAGeCra7hXggW4qcbgeEWKsyBkrjz4EmdPxEEzp/PMpdQrSbKk/keoV3/f6yKlue7MctEmJbQQPa2SYvPlP+T4gT+Q3rU2mhEIU5QhTn+h8vsfMNefEN1isTDr5wu4U4MkoKHfOXyrZVIkTUjOgoCbZ0FAhxWVJV4HHONCFmzfAdoOwPliwhfr5THuefbhKiaJey/S7OE+KjmVLkfGCS/2N/mTlzppjyeLp43BOT89JXQlTYfO93fVQXOWkb1Qcw14bDLcTkk1LkdNkvRKaPoH4NjepoQsdPNKHzx+OJHHmy1Cm86/eXnAoh/rK04ZkmOUVCTH1DXoBWbAlc//zFfEiIRSlVWT9Kp4JqUrFBipwDPYVwlSnf3p+RDTXTT0qw/CLElBayzUERI4T5yyrRWD1A2V+BULRDiptFY6XYsWSd/n59AuzxqTPFe/+SIuc/6UKYC5Xtz4ET4lR21ctfC2G1Kdvey9kWOU63EFOriZyfNZGj4Qea0PETTej8sfjaLE+W8Y8E5uJdE6dbiGm/yRP+I2vUD927XFWp5K9/K4TrLKS818f6f6lLba6NXJMUOjnT1feprtGaQGVd+cJpFiLnSdn/ybEzxdT26cK613f//BU5WUuFWDxeipw1D8vpKyWcOCCzql6dLKcjLQpPUbuPCjHNI6Rf+0Z9dtXGvBoiJ8DTVS531c1J5/1CLNNEjoafnJOsKw2NC4kDdpllVfR6BkWvnJ0Kn18fklVgwwwwvTcEqSxc9svv0iU6zAh/vRL0Z98W7jSKd8OrX2bw6f5UTI8HpkiYbZf86/WzUkNd1aP9yV5TQ/kKyEsFVy6gA9PUNOIf8+2w7nK5MBh8+4Z4PbpmmtLY8wEc+Va+3nQw9HwMghT4Ph3YCks/lkUhE1vD+PtllpW//H4APvpRmmx2biltRIwqzuzrcuHpLdJ2YUAT+Gcf6WMVKNwC/pkH35SBAXgpEUb6sDvR0Ag0mtDRuKhIS0tD6A1s/JuJE69JkTNrVpXI8V4gGyp29pvhU4/x4QNdIVFldeUCMyz0WDxMHAIxCk0gLSWy+nG/UWBUaZp48L/gFi4eHJtO+nMNFw5CgOOgXA5up65PUHf6d11lAdT8T50FUPAMWKS/JcGtockzENav9m3q8+iaOSOdjWlQsEm+3v4W+dApEAdbf4FV3wJC2nhcdTsEGf3fft1u+Hy5/J/0bg+3j4YgFZ5ua3JgzlZwCbisKTzRO/AiZ0Y+zPOInBebwtgAGKJqaPiDJnQ0Liq8F5+4IggTLp5MSyc11XeNFbV3/TYXvPB71Ul/ZDN1fRUC/rcCHE7o0FzaPChlzXzYtwmKcyFFhV9uyT7I3wi3dUrj8nd9r6NUOLjywG0B9GBUaQFwrrx2hBPM/4WiueAuk32OuRPiHgF9WN3b1lXLZ8bUdEbkmig4DvoQ6DkVki5X0C83rFkAW1bI5z2GwNAbQK9AXPyyFb5ZJZcHdoWbhyvb3svKbHh+mxQjVyTBYz3V2UPUhluAKR++KpXecy82hRSVNhoaGmrQhI7GRcXgqSbiCqHolVTuMaUzZ6ZvkdOQi+Wn+yGrHOJCfPtY+cuGPbDvOAQHwc0jlLeTc1iKHHTQX5nZ9ikOz5N/k4dDRJK6NmrizJZ/g5qCLkT59ufKa6fiNyiYDfZ98nlId0hIg9Be/rfhy6PrsTvTufywiXIrhCZA31SIVjCy5bDDT59IY1aAQePhkiv9/34IAd//Jl3vAYb3hmuGqPueLsmC13bIwkFXNoNHeoAhgFOrrmrTVZrI0ThfaEJH46Kh2AUz8iDuIRN9w+C9jFQ+fjawBpH7zfBsRhroDXzygomoGtMI3umx+ioul1vhuzVyecylkKDQ008I+PV7udzlUmiiwovWVgS5v8rl1tcq37429J5Kzm6Luu3PVvyNF/t+KJwL5T/J5/oYiH8Uom4+3VXdX6rHEgUHGRmRY8IFxPeCXk9yhhN5XVjMsPBdyDsmK16PvAU69/d/e5cL/rsC1u+Rz68eCKP6qRM53xyG9zztjG0BD3YLbPyYS8D0PPjOM131QlO4WhM5GucBTehoXBQIAWn5UOCCdsEw/2kTUS/VX0ZfCU43vLId0BvY/2kqSzvCAB+xGf7YJfzwG1gqITFO3nEr5ehuOHFA2joMGKd8e4CspdKlPKYLRLVV14YvDPHyr7sMhE35qE6g42+8OI5D0atQNh/pV2CA6Fsg7mEwKBAjNfHGEgXrjTicdj4/IGOJOtyiTDjlZcEP70C5WTqQj7sHmikYCbLa4YPFsOeYFCQ3j5BTVkoRAj47AJ973HRuaAN3dgqsrYNTwJO5sMAiRc7cRC0mR+P8oQkdjYuCb8pgkUV+YZ9vCs897Z+PkRK+PAiHyqDPHSb+0la9z9LRHPjV44n1l2HKg0OFG371BM32HAqNYhXuCFLgZC2Wyy1VCqXa0Ecj/xFOcBZCcHJg21eKMweK3wHz/wCHfC1iFMRPAWOHhrXt/b/f3i2dm1ua+OJYBv/ZmUrHI2Ay+P99O7AFfvoUnA6IS4Sr/y79yfzFXA7v/ABZeTKj6s4x0K2N8v1xC3h3N8w/Kp/f0RFubBtYkeMQ8HguLPT8XucmwhhN5GicRzSho3HBc8IBGfly+ZE4+O75wAeyZlngf55Moge6wrArTYQGnR6b4Y/IEQLmrZIxD/07QQcVgcyHtkNhNhjDZLaVGor3gDUfgiIgcYi6NmpDp4PgVjLzqvxHGdx7PrDuAPNHULaIUwInbDDET1UWh1Mbs1IzSMtI5dYOUuREd4I3PzLR7m3/Y4mEG9Ythg0/yuctO8OYOyGknkDo6mQXwtvfQ3EZRITC/eOhVaLy/XG4Ye42+OWkfH5/F5jQWnk7deEWMDkHlpfLi8sriTBaEzka5xlN6Ghc0LgFPJUH5QL6hkLu6xmkzQxsIKsQ8PpOOdzeP0Fmnnjb8cdlujpbD8CRHHnXPWGwom6c6stGT2xJz8vlFIcacj3xQU0uBYOKgOH6iLkD8lNlLEz4CPXZV0oRLihfDiUfgnVj1euh/WQmVfjAwHxOyR44/L2LWzukc0sHE+1uhPZ/BX2Q/7FEtgpY+gkc2Smf9xkBg8fL2Bx/2XNMTldZ7dAkBu4bDwkqpuEqnbJGzuYCGWw8pQeMUJlNWBd6HQwOg9UV8HoiDNPq5GhcAGhCR+OC5nMz/FYJoTp4tgl85A58IOuyE7C9CEIM8I+uVcP4ddV58YXDBQs8wb8j+kK0ijvZ4/tkoGpQMPS+Qvn2IMVS7m9yuakKseUPUTeBZRFUroW86dDsM//rx6SlpWEwGHwey9qCve0HoPQ7sMyXU1UABEHkOIj5G4T2UNb/2vrgdsHUmzMo+N3FX9unEZoAvR6DuBrt1yd6C7Nh4ftgzgdDMIy4CTpfqqyPv+6AL1eC2w3tkuGecRChYCTIi9kGMzfBPrP8jv+rD/RLUN6Ov9weA1dGQHLw2fsMDQ0laEJH44IlywHPFcrlJ+KhtTHwgaxldnjfk3lya3to6ikMqKbOy5rtUFgKUeEwso/irgCwebn823UQhKnMULEcgcpcWd+l8SXq2qgPnV4W3Dt2tRxZyXkIEmZBkB8X0PoK8XmDvR0noHypDC627azaXh8D0ZNkoHGQiimc2vpQdhimTsrgvXVyuirpCuj2IAQrFKx7N8KK/4HTDo3iYNzdyrLmXG6Yv0bWyQHo1wkmjYRgFRljJ8vBtBGyKyAqGNL6QecGBGZ7EaLuuB5N5GhcSGhCR+OCRAhp8VAp4NJQ+KvC9Gx/+ewAlDqgVSRMbC1fU1Pnxe6oqmsydgCEKKhu66UkD47tAXTQZ7iavfG046kbE9tZmRWBUoKbQ5NZkDtdpnJXboDG/4RG4+rOxKqzEN/t6dzvMnF0FDiOVNsoCMKHQtS1EDFcXf2e2vogXDCpo4mM9Aw+2SsDj//9qonk4cqCdF0OWD0ftmXK5807wpi/QZgCoVRpg4+WwO5j8vnYS2V5AjXBwntKYNYmMNshIRRm94cWDYiXOWCHCD0kBQU2eFlD42yjCR2NC5Il5bCyAoKBjCZnxx/qmAV+8FxQ7utS5WWlps7Lqu1QWgFxUTCgi7r+7PDWvOkKUfHq2gCY82oapQcM/PNa/6eG1NJoosxsyntK+l/lPQH5MyCkN4RdCmH9IaSHHAESbsAp42yevMWEfb8n2DttNna3nUfi07njVxNmz3HAAKF9oNHVEDkWDHEB6fIpTCYT1kKYOSuVDN1snMLO/cPTmfu1iVCFn2UugCUfyhRygH6jZVkAJZWK80rg3R8gt1gWmbz1SuijMmvst1x4bivY3NAuCmZdAnEqRW+OEx44CeVusAroaIQZjaFl8Ln3bdPQUIMmdDQuOIpdkO7JsrovFtqqGB3xh/f2yGDngU2gd7VUX6XTY5X2aqM5l6rzGnI6YPc6udz9MuXbV8dRYuDT/akk/grP36uuDpASQrpC86+h+D0w/wdchWBdLx/FdWx3ByaeRYqcYIw8FGsiuA2EXwZhl0HYADCcpQJzzgrY/xlcfsBEkEfkGIONvLHcpHi04sBWWP452K3SjPPKW6FNd2Vt7DkmR3IqbNIP7d4UaNFEWRtevj8Kb+2SmX/9EuCp3hCm8kyf54S7s6FPKNwQBSUuOZ08NRcejNNMOTUuDjSho3HB8WxBVWHA+1XUkPGHTfmwMR+CdHB354a1tXqbrITcJEbGU6jh8A6wlkNkLLRSUQSuOjc2M1HZAV74KJWotsrrAKlBFwxxD0Ds/eA4JKexKtfLhyvXx/rh8KYjAwd2jEFG7E47n12bwcznAt+36gghM9J2vQ22Qvj8QIYUOZ6A89mz/a/H5LTD6u9g+2r5PKkNXPU3ZXWPhIDlm6WlgxDQqql0H49WISBcbnh3Dyzw1MgZ0wIe7Now36pCF1S64Y5o6OCZLuwfJsXPJyUQrYd+KgKkNTTOJZrQ0bigWFcpHY51wDNNICSA5oJeXKKq9P34VtCsAXelDhes9HgWje6n/qJyYIv82/ESdcaMXoQARwXc0t5E2xuU1wFqKDodGNvJR/TNsj+iHDB4qggHyWmsjIwMXvQR7K1vFFhTz+qUn4Rdb1S5jX+Vk8Gn+9XVYyrMhh8/ln8B+l4JA1PAoGA0z2qXzuNbPRWKB3SBG4fJaSulVDjh2a2wwTMS+reO8JcAFAI84QCLG2I9+2VzyzidfzeB+07Ct2XQI+Ts/E41NAKFJnQ0LhhsAlLz5PLNUdDnLN0pLj8ORy0QGQw3KyjB74uNe2RsTkwkXNKx7nVrS2l22OGN9zNwOl3c9Fhag/rjtiPtD4AZM0zMebFhNhlqUsGro9OBrkYA7Lky9fTiqIBDX8Dhb6WbuS4Ilhgz+HCL8j4IN/yeKStXu5wyM270rdBSYVxWXjG8twhyiqQ4vn4oXKbSQDavEtI2wZEyKTim9YIhKrPRanJZuKy780EJPNFYtu8Ucjp5emMZuzMpCrqfxaB3DY2GogkdjQuG94rhkAMaG+CxBgTj1oXVBR/vl8s3t4NGDYj/cXumHQCG9a7/br62tOrHH83gh19TuXZYOgkqzDur47JWLT/zfMNtMvxNBVfUx7Ns6ulFuOD4T7DvE7B7goUa94Wu98Pa91yk91LWB4sZln/myYwDWneDkZMgPOrUjsGqVXDyJCQlweWX+/xSbDsEn/4kR3SiIuDusdBGpbP87mLI2AwldogNgZl9oaPK9HGLG74plQHHvUNlgc4wPfw9Vv42e4ZWWTkIIeNzOhnh5wpN6Ghc4Ig/EWazWQDCbDaf765o1OCkQ4juB4Rov1+I+aVn73PmHRJi7CIh/vazEHZnw9rafVSIya8K8fhbQlTa/NsmPT1dACI9Pf205+MGpItV3zasP0IIYS0RYtFYIW7t4PtzvM+VUFuf1bR1rsjfLMSqB+WxWDRWiF/uFiJ3rRBut7r29m0S4p0nhXh1shBvTBViW2aNtubNE6J5cyE8s3UC5PN5806t4nQK8U2m/M5MflWIuV8JUWJRv49Ls4QYv1h+nx9cJURehfq2dluFGHpYiJuyhJhwTIhO+4VY6enbIZsQT+YIceURIbZWVm1T4ZLrfqOdTjXOAQ25fmsjOhoXBC8VyjvJfqEw/ix541id8PUhuTypvboCbNVZ4zHu7N8JQv0cGao+PeKNnbl2eDrDu5toUc/Ulz8Yo+B/h2XsyYxpgZka8tXncxHvo4aSvbDvIyj0xE0FRUL7SdDqatCrKGJXWQa/fF0VQ9WkBYy6TRpznuKbb+CGG6S8qc6JE/L1r7+meNR1fLhE2oOAdLQfP1hdhp7LDe/vhe+OyOeDm8K0nuozqzZUwiM5MqvqgVg5inN3tnQeHxoBbYwwKRrK3HDHCZjVBBIMsMMGJx3yfQ2NCxmdEDV/nX9cSktLiY6Oxmw2ExUVVf8GGueE7Va47rhcntdcDpGfDb45LIOQE8PgnaFVdXPUUFoOqR/J8vzTJ0GyAidqgJCQkFPTSi/cZ0Ong7/PkUaeDeXOfmk4Sw289LmJhH6nv9eQOjrV+2yz2Rre0QBSdhT2f1xlfaELgpYpUuQYVf7UD2yBX76CSosMoO43CvqPqTEb5XJB69Zw/LjvRnQ6HInNmTnzMBa7gTAj3HIl9FIZG1bmkEHHmwvk87+2l6JdbT0bqxteK4JQvSzlYPS0k5YPyUGyVs7AMIgxQJkL5hbBT+UQrpPp6+kJMDBc3WdraCihIddvbURH47wiBPzbc9Ke2OjsiRybq2o05+Z2DRM5AOt2S5HTOlG5yKnpobVkfQa3X28KiMgBeGhCGnlroTwbajoyqB2FUer7da4oOwqHvoTsX5BXXj00GwkdboGwpurarCiDlV/J+jgA8Ulw5V+hSUsfK69aVbvIARCC4JNZJO5chXXwMO4aC41VVvk+XAazN8PJCulZNa1nw4OOQ/Vwveea4RU5bxVJj7lLQuGoAxKD4PZouCYKUhPgnhh5qKMM0EjLttK4CNC+phrnlcUW2GiFMB1MO0sByAA/Z8uAzYTQwLg2b/LYLAxSWPOmehCvzWbjnlvSWbQulR83ZjS8Ux6i28u/OasD017NPqenp5OamkpGRuD6rJSSvbApHVY/ANk/A0IamF7+BvScok7kCCGLNn72bylydHpZ4fimx2oROSADj/2gf/RJptygXuSszIapv0mR0zQMXhzYMJFjrzaO39ZYVZRzuxUWWqTz+LvJsKq1FDPLymVZBpAWEM2CNZGjcfGgjehonDds7irTzntj5QkUGp7SXBO3gG8Py+WJrRs+mpNTBNmFst5NTwVTEL7Sqq8dYSL7IHz4ZSptugcmrbrZaDjwORTvAPNBiG5ACv25TgWvCyGgcAsc+hoKPaMt6CDxMmh7Y5XAU0NJHvz8pXSPB2jcDEbe4ocZZ5J/6VKDRiapOtu63PDhPjntCtAnHp7sDVEq42KyHPD0qqzcAAAgAElEQVR0gRQ6TQwwIgJGV4uJa2uET5rJqSqnkKnl/cPgq1KocEMjg+ZzpXHxoQkdjfPGhyVwwimHxu+plhIb6JTmtbmQVQ4RQTCmecP77R3N6dISIhRMtflKqy7JgzGXmujUP3Bp1WGNIfFyOLkSjs6HnlPVt3WuUsHrwmmF7BVwdAFYPN5kOgMkD4e2f4HIBqTku5zSMX7Dj3I5KBguHQu9h/tZ/O/yy6F5cxl47CPcUeh06Jo3l+sppMQGc7bCtiL5/C9t4faOUnyoYZcNbj8BwyKgbTCsqYAXC2G7rWo0NVwnR7JAVg0XQnpdXRkhRY6GxsWIJnQ0zgvFLnjLU9vksXiZ6eGlLndrNdk+8zx3w1e3gnAVmTc12eKpw1NfgcCa+BqFMnvik1JTTcTUDKhpAK2vkUInewUkj4DGvdW1o9T3K5BUnIRji+D4j+CwyNcMYdB8FLS5Vn0MjpesvbDyayj2WFS07AzDboRoJTFXBgO88grihhsAHTqqxI7Q6dABvPyyspLJwM5imLMFCm0QZoApAYjHWV0Bl4bB803kqMwt0fBdmRQ7cQa4s9rNhkvIOJwXCqW57qsBKkCooXE+0ISOxnnhvWIoF9A1xHc6eaBSmg+YYXeJvDsd36rh/S40S5dpvR66tW5YWy6XNPMECA2wOWJMJznikf0zbHkaBr0IkbXFmVxAuGyQ86sUN0Xbql4PS4RWE6TICW7gsSorkh5V3mDjsEi4/Fro2E/5tIzTBYsTryPvrq+5bt4jxJZUBSbrmjeXIue66/xuTwj49gh8sFdOubaIgH/1gZYBMDfdZ4d8Z9U+xhjgukYybXxuIXQ2wqBwmYn1iVkKo2MO+DAZ2msp5BoXMZrQ0TjnFDjliRTgkbjaU2NNJtMpkaPWwuAHz1THZYkQF6Kyw9XY7WmvTSKENbA9l71qOTgAI0016f4IVOZC8S7YmAaDXoIQlVVzzybCLYOLs1fI7ClnuecNHTTuA63GQ0I/j1dWA3A5YPPPsHGpNOTU6aDH5TBgnHQdV8rJIvhkKRzPB3pdR+gNE7nBsIqQgrorI9dGuQPmbodfPSNMVyTBw93V18epSc8QGaOzxwadPd/dKI/YOWiH90tk1mOEHrqFQJEL3k6SmVkaGhczmtDROOe8XQyVAnqFwPA6LjANTWkus8MvHtPFQIzmAOz2OEN3CcDoiMMzmoMO9Gfhl2gwQl8T/DoFKnNgczr0nQkhKjN/AolwQ8kemRmWsxqsBVXvhTWFZqOg+ZUQ1iQAnyXg4FZYswBKPcHvye1g6A2QoCIDz+2Gn7fCwrVyRCciFG4aDr3bG4Bhqvp40Az/3iqzqoJ0cF8XGNcysIG/fUPhnWJYYoFmQVUxN808tXI+LJF+cxHA4HD50ND4I6AJHY1zSqET/lsqlx+Jr/1EXjMmR4mztJdfToLdDW0bQReFIxm+Mr+EgAPZsH5JBoajLkb3T1PWaA2Ex3zzbGaxGKOh3yz4bZoUFivvhNYTofV1YAzAdIgSHBZZsbhgM+StB1th1XuGMGg6UAqc+J5VAbENJeeInKY66amhFBEFl01UN00FkFsMny2rqnDctRVMGgnRKqfThIDvj8pClk4BTcLgn73V+1XVbNu7j0JIP6r7YmXdquQgGBtZJXY6GuUUVqlLxutoaPyR0ISOxjnlY7O8a+wRAkNqKZAXqJTm5Sfk31HNlV/UfGV+5Zth5fwM1i1KZeylys0saxLmNUh0g7W86nmgiWwBA56B7a9A6QE4+AUcWSADeltfA8E+PjcQKf7OCjDvh6IdUtyY91aJO4CgcGgyQGaINe4rR6ACRWkR/PY97Nvk+axg6DMS+o4EYwOmHMsr4WiOtPy4dggM7KpeqJbZ4eUd8JtnqmpgE5jSo2FGs4VOORUcWy0N3C3A28XbYuCwA14rhgIXXNtIrru0XIqfWE3kaPwB0YSOxjnD7JJCB6Qjcm0XiECkNB+3wD6zPOlfocIZ2pewSp0pRc6YG9OZObPhGUeGIBmEbC2H8tKzJ3QAotrB4Fcgbx3s/xTKDslaO0fmQ0J/OYoS1wPCk+X/RUmKvxBgL4WKbLAclfE25r1Qdgxwn96PiOZS1DTuC/F9wBDg2KSKMhmDs301uF2ADrpcCgNTILKWUZLDOVBRCUmNIa6eUa62yXDzSOjcAmIbMCK2q1haOeRb5VTVXZ1hYquGje69ViSnpSrc0CUEUiIhpZH8Dbg9TqN6naxuHKKD5eXwf8XQ1SgrIL+bDNGa0NH4A6J5XWmcM14ulCfWzkaY30K9P48//GcffHEQ+ifArH71r18b3gu7N05owLh0HnvCxA1XBKafnz8DhSdhwgPQqktg2qwP4YbcX2H/Z1KYVMcYA1Ft5ZTXe5kZvP5DKo/elM6U20y8/GkGc/+XyuSJ6fx9pAlHGVTkSIHjrPD9WaEJENNZBhU37tPwlPDasFXClhWw9WdweIK8m3eEIRMhoZY6OyUW+PQnyMqXU09lFXBVfxjcDYxnITgcZNr2Vwfh0wNSfCSHywKAHRoYN/VqoZwSnt4YgpEF/orcsv7N5Di5jtMzsuOtw5PtgE0eoTXI42eloXGhonldaVzwlLvhU89ozoN1ZFoFil+9LtHJDWuneuZXUJCRS8eYaKbQ26ouYppIoZN37NwJHZ0eEodIy4SinVD0OxRukzE89hI5zQQwDhMlHeDlL1J5/cvZOIWdWzukM9ZuImvxme2GJkBEMkR3lOnt0Z0hNO7s7ovdCttXyaJ/Vo/YatISBo2HFh3rHiFZs0OKDdNtYHfA74fgxw1ymyt6Bb6vuZXwwu+yRg7IkcbJ3Rpe28nqhs1W6Tw+0TPK1C8MviiFeaXQ2CDdx4M8x6LYJaeokoPlQ0Pjj44mdDTOCV+VgtkNrYNhVIBrxtTkuEVWQg7SyRGdhlAz82v9kgweuiZwhfJadISDv8vidf2vClizfqHTQ3wP+eiArGFTdkQaZTrKZJr3U+Umvpg8G4fLTrDByJMPmwgKk8HDweEQlgQRSbLOTSBjbOrDVgnbMuUIjlfgxDaFgVdDu551CxwhwGqH/cehZVOIDAPCYHhvyC2CrQegaSx0DmDdoZXZ8PpOKHfKAoD/6AYjkgMXiF7kgvxqM7pNguCGRnK6+FOzTCfvEwrzy2R9nAdjobVWG0fjT4JWIUHjrOMU8FGJXL4rRn0Je3/5LU/+7RUPEQ24Y60ej2K12hh8dTrrFqXy4VuBM7Ns0Un+PXkYHLaANasKQ4gciWkxGtpeDx1vh3kFGTj+v707D4+6vNc//p4sM9l3CPsiAWQRK6AcFBC0oraiaGtrrUtb9RR39LQ/l54ETOzBFuvSVilaD7ZHe/C0Llh3quJSpCyKIIKyGiAsAbKTdeb5/fFMSIgBMpmZzJL7dV1zhQzJ5GES8r3nWT4ftw16je4GXiwr4uTrbHfwQTMhd4ItRNhVIafuMPzrdfjTXFjxqn0/o4ftLn7lPZB3avvhYd8hW6AR7N8nuuzSVaZ3X1Rjk3179jfsLM+GHdAYgO4WhxvhN+vgV5/akHNyBvx+EpzbN3AhJ8YB/eNhe4OdrWnWOx4uSoW0GPjAW5uozA2f1kG5p/3HEolGmtGRoHur2va0yo61pzyC7V/eoPNvfuwHaXvyq/IwjJuej8fAA/9VQFJCYFogpPeA1Cxbrbd4EwwJwpJJZwXiiH+gVB2Ctctgw0ctgTAz186CDR1rK1W35/118NYqSE+xQee88TByECQ67YzNio0w7TSIj7MzPb2z4KTesH0PHKqAXD+W3tYfhN+sh/219hXl94fAD/Labyrrzyk3pwN+lAE/3A0TquCq9JYQ9Y0EGOy0bRxuy7YfNz1ZS1bSvSjoSNA9692b84O04FdZPdwIm7yzR+P92EvT9uRXnXeD6+QZ+XxrQuCaWTocMPQ0u8fk47fhpBMsu3SVcOlaXrobPnkbvvy45Wh6dh84fToM+caxAw7AsrXw0edw+dmQkw5rNsPz70NFNZwzFsYOg1Wb4NOtcOoQW/wvPg4mjrIFAZs6+S1ucNvN8C/tsP2ieiXCf4yBUccJTf42sj09EWZnwa8OQm4cnJ1sT1YBDIiz1cjrjb1PIUe6GwUdCarNDbCyDmKB73VBRd71h1pOs+T6Udm17avn5qWN+LjAX+BPnWr3m+zdAVvX2eWXUAtl13KPG7Z/Bus+gF1fttzfbxiMPQcGjDh+GPQYu7l43TY4ZTCcmmfv79sDtpbYENM/F4b2hTFDYMk/bdCJ9/42jHFAQjxU1NjP8cXWCpi/Doq9TUgv6A/XnwxJJ/hNG4hGtjdlQUkT3LPfFgY8PRESHfC3KrsvzhUGAVokFBR0JKgWe2dzzkmG3l3w0/aJt9ruN7ID+7itg06gpaTDaefAqjdh+csweGTg68v4KhRdy2urYcNyWP9PqPaeTHI4IO80G3B6HmdzcE2tPRIeH2eDSmwMlByAs8fYv3d77H2ZKbYx64frIK+PPU7+2xdstePJp9jNyR9vhp6ZMMSH9hCNHlvO4Lmt9gh5phNuOwUm+NDCIhCNbO/vCakx8GY1LCiDzBgYlwh3BfCkoEikUdCRoKk39pQH2GWrrrDukH37jQD/Ynd7l02CtZF67Ln2Il9RCkufhfOvCVwbhHBmPLB7C3y+wnYTd3sDZUIyjJoIoydB2nGWfCpq4K/LbGuG5AS77DRyIKQmwbD+8MYq6JEBfXJg+QYoq7Z7dPYchJ37bbD50QV2SevJV+0JrAMVdrnL1cGwubUSHl4H27w/62flwi2jIL0TFZg70si2yg2uGLs3pz135dilqr3eTuWjAtDMViSSKehI0CyrsUfKc2O7pkFgXRMUey82vva2OpHmTuWHg3QyypkA510Df18Amz+G1EzbkylaVZfBxpU24FS26nnVcwCMmWw3GMedIGhU18LTb7S0Y1i7Fd7+2B4bv+o8+M5keOR5eOo1+31rcsNPLrTBZ96zdvMxwLB+cNPFUHLQBqfTT7azPyfSdhYnLR5uHAlTend+n9WJGtmur4Pb98L0FFsc8Fhy4uxNRBR0pAM6eyLkBW/zzotTg3+kHGBrle04kOWC7ITAPnZyq6DjMcEpeDhgOJx7JSx9xm5MTsmAUwNUgTkc1B+2NYO+/NjuvWkOGvEuGDYORk6EXB86dm/fA3sPwc+/D1lpdqbmo8/hlY/s7M2Zo+DWy+wMTW293YcD9nvojG/ZYA72VFa6Dy04vqyAR9fD9lazODePggw/Zk+Od8rtP/8zn0Xl8OBBaMQuTd2SBSndYNZPxF8KOnJCnTkRsq8JlnkLuV3WRctWm72nrYYFYdNzkrcBqTFQVw9JAQ5SzU4+A6rL4aNX4P0X7NcbM+X4p4vCWWM97PjcNtfcscHbf8qrzxAbbvJOtWHHV01u+32IbzXzM+Yk2LUfXlthg05W6tf7V63YYPtUDerl+9esa7LtG17abkN1Wrwt/je5l3+n5U50ym1JJVT81L5/fjL8sqdCjkhHKejICXXmRMgr3tmVcQmQ10XF5L7ynnQ5KQjBKj7W7gGpqbMzBAOCFHQAxp1nw876D+GDF+CLVXD296DXwOB9zUCqqbSnpravh51fgrux5e+ye8PQcTDsNFtDyB9NHjuztvcgpPaz9yUnwPiTYf12+GAdTB5jw2J5DTQ02OWtD9fb2jkup/27jgaUTw7A7z6DvbX2/am94d9HtD+L4+ss6LFOuZ33H/k8WQbb6t30dsC9OXBlWniUIBCJFAo60iG+ngh5xRs6ZnRBgcBme7wzSL2DtB+of0/YVAzF3k2sweJwwNnfhaxe8NGrsH8n/PUhuzl34gxIDHILDV+53bBvh21j8dVG2NemUWhatt1zM2wc5PjQe2x3KcTFQlpyyx6p5uUuhwPGDoUlH9rvycDclkacPdIhr689Sj7pFPuxu0vhvU/hUJXdaDxmSMvjnEhFAzy1Cf6x2/v4CXaz8enHOVHl6yzo18oZGPjtIVhYBq6b8hkZD4/0sq0cRMQ3CjrSYR05EQKwowE+q7e1cy7wYd+Dv5qDTp8gBZ0BrYJOsDli7JJV3jfgwyV2VmfDcrvHZew5dokruQvqErXH3WgL+e3ZZmdsSra0dAxvljsQBo+GwafYWRxfZiAOVtrj3rsPQEay/dzLp8KQVr2hPN7j4lNOtYUBRw6yfw/25FST23bqbv74vL52H08fH8oOeAz8Yxc89QVUNdrHmzEQrhlm6+Icb9YGYOrUqZ2qi7OjAf5jH6zzbnz/Xhr8IgeStFQl0ikRF3Qef/xx5s+fz549exg1ahSPPPIIkydPDvWwuoUTnQhp9pp3NufMJNv2oSs0uuFAnf1zIGZ02ruI9fe+gl+0oIgvPzh2Sf5ASkqD6Vfb2Zz3/mo7nS//Oyx/BXoPspWUh4zxfxnoWJoaoXw/HCyBvV/Z2ZrSXUfvtQF7HLz/cFvUb9AoWxuoM9weWLoaklzws+/Zwn8v/RNe+hC+Oc5uKPZ4WgLM9PGw5gs7W5Pkgt7ZNqA0NkGvVsfSE5y+hZwdVfDYhpZO44NT7SzOiMyWj+nIrM0555zT4VlQY2xxv/tL4bCB9BhbF6crXyyIRKOICjrPPfccs2fP5vHHH+ess85i4cKFXHjhhXz++ecMGBDAVsPyNb70PXrb20Dwgi5cYqlstOX2YxyQHoA9Qe1dxPL6wOo3i/jo1QJOG3r8kvyB1jcPvv//4IuVtt/T3h22Eeie7fDPJZCUalsjZPe2t6ze9j5Xoj26fqyaPO5GOFwFNRVQXel9WwZl++DQPqg80LJc1FpCMvQaBH2H2oCT0zswdX8O19mWDDMm2mPgYGdz3lhpA9Apg+3GbI9pmdW5bDIsXWML/00eAzv22qPiF57h+9evbYLFW+GF7fbIeEIsXDUULhn49SPnHd271pFZ0HI35O+HN7z/d85IgAdzbWNOEfFPRAWdhx56iOuuu47rr78egEceeYQ333yTBQsWMG/evBCPLnr50vfoQFPLlPvULgw61d4Nrylxgdmo2d6/76EHbciZ8K1CJl3UtU0tAWJj7SmlkRNtGNm23raM2L3FhpXDX9h9Ml/jAKfLBh6PxxblczfZkNNeiGnLmQhZuXY5KnegDThp2cHZEFtda/fkJLfa7N0zA07Lg5eX25mbaad5j/d7v/7wAbb55ofrobQc0pPh6umQmtjxr2sMvLfH7sU56P35nZgLs0ZAj+M8zon2rnVkFvSDGrh7P+x321/Is7Ph+oyuKckg0h1ETNBpaGhgzZo13H333UfdP336dJYvXx6iUXUPvvQ9ete7T+YUF/Tswp+uam9F3eQAvgJu7yJ24+2FxA3JZ/WXcNGZwamn0xEpmXYPz5gp9gj3wb1wqMQubR0osTMydTXeSsMGGursrT0xsXa/T3Ka9206ZPS0m6Gzcu3yWVed8umdbY/vF++HU05qOVY/uLfdZ7NhB5w52lYt3lhs6xsNyIWMFLhoot0YHevjcun2KliwAT7zLlP1SoSfjoAJHdxwfqy9ayeaBa31wK8PwjPeNiknxdtZnFOCeKJPpDuKmKBz4MAB3G43ublH//bJzc1l79697X5OfX099fUtpWwrKyuDOsZo5Uvfo/e9U+9Tu6AScmuHm4NOgH+i217EHn0wn//8byivtkssp+UF9ut1RrzLHj1v7/h5UyM01EK9N+jExEJsHMTFed/G2xmbYASZzhaaPOsUO3Nz1mi7gRjsBuNeWTYAHa63Hcj/tswGoCu/2RI4fQk5VQ3w7BZ4pdguhbli4HtD4DuDwenD47Q3awMcdxZ0TxNs+Uk+270zkdekw8+zIUEbjkUCLuL+Wzna/EY2xnztvmbz5s0jPT39yK1///5dMcRuyxj4l7fGiD8tH+bOnXvkYtFWUVFRuxfH5h/kDqzE+KTtReyBeUVHGkW+8pHdBB3O4uLtjExmT1t1uEdfO0uTlm1nblxJwZutad7n1PZ72TyzEXuMVDJljF2VWvYp1Lc6zZWTZht1xsXapptnjrahyNdZtSYPvLwDrn8fXv7KhpyzcmHhFPhBnu8hpznQ1NfXU1hYSEFBAe+88067s6D/7z/zueDuQv6vzM32RtseZVEfyO+hkCMSNCZC1NfXm9jYWPPCCy8cdf9tt91mpkyZ0u7n1NXVmYqKiiO3nTt3GsBUVFR0xZC7nS/rjMnbbMzoLcbUeTr/OIWFhQYwhYWFHbrfGGM+KTXmwteMufGDzn/dE42j+f38OYXm3j8ac+tvjXn9X4H7euFuzpw57T73xtjnZs6cOe3e395zeKzHabZqkzE/W2DMO58YU1Nn73v+fWOees2Y+obOjd/jMWblPmP+/T37s3Lha8bMet+YNaWdezxff04/qzPm21/Z/yN5m42ZvceY8qbOfW2R7qaioqLT1++IWbpyOp2MGzeOpUuXcumllx65f+nSpVxySfvdD10uFy6XKmx1lZXePSBjE8DlxyxBZyoxN2/cbO4y7q8TbcCedRvE5+Xz1mpbfM6Xo8uRqjOtQNrb53TPLwoZe14+C5bArIvbn1EaP9x2GH/vU1i50Z54OlABPzi3pTCgL7ZX2no4Hx+w76fF23o45/frWAPP9nR071qjgT+UweOHoAnIjIH7esKFOjYu0jWCELyCZvHixSY+Pt489dRT5vPPPzezZ882ycnJZseOHR36fH8SoZzYnXvsK9XfHgzM4zW/MnY6nSecBfiizL5Cv+rtwHztE81eFBTMMY+/ZGd18v/bmENVgfm64a6zMzTN38O4OKe57bf2ebv1t8Zs23Psz2lqMmZXqTEfrjfmH2uMaezE7MeBWmMeXmfMt7wzODNeN+aPG42p7uSskK821BlzcXHLLM7NJcYcaOyary0STfy5fjuM6cgB0/Dx+OOP8+tf/5o9e/YwevRoHn74YaZMmdKhz62srCQ9PZ2KigrS0rqo02Q3cv5XsK0RnuwduKPlLpfryB6Z1hvL2zpYB1e/a/drLJne+VfpvqiphUeeh31ldqPs7d85+lh0tGqewWnet3SsWTZjbD2bn91dxOKnCoiJdeJxNzDhW4X88Pp8Jp0CYwb7fkqqI2qb4Pnt9lbv3Uc1uRf8aBj07oKyB/UeeKwMniyzszgZMTCnB3w7RX2qRDrDn+t3xAUdfyjoBE+1B8Zus5uBPxoEOQFYFO3oBRVscbeZb9q3f54KOT7UUPHHoSp4+K9QUQMn9YabZoIzYhaEO+94AfRgJaz50lYsXvKXIv71mq09dN538tm+oohnnzxxG4TOanTD6zvhua1Q5t3IPCIDrj/56KrGwfRJra2Ls817ouqCZCjoAT26wc+FSLD4c/3Wfz0JiE31NuT0igtsyDnnnHN4++23v1aDpO3x5FgHZLmgtA7213Vd0MlKhRsvhkefh217YMESuOo8yI7iHN3ecerb7sjn023w8Zf2eQBY+YYNOZdcVcj9hfmMGAixP8lneP+vF5r0l9vAshL4n82w33vyr1ci/GQ4nNWra2ZRajzwyEH4U4X9v5ATC3N7wPnaiyMSUgo6EhBbvK+ehweg/ULrkPPOO+8cVU22oKCAZcuWHTm+29qAFBt0tlXCyC569Q7QJwf+/SJY8HfbMfuBv8DMSXDmqOhbpmi98fjm2fn87C77/msr4PQL7PfIAQztBzWD3EwvKKTwvhMXmuwsY2DFfvjzl/CVt8dalsseE5/eD+K76Mj2shqYUwol3npOl6bCvTmQ0UW93kTk2BR0JCB2eKfpBwegMnHr0yxtZ3Jah5y2swHDM2DNAdhUDhe1UzwvmIb0hbuusF23t+2B596FtVvsKaGs1K4dS7AUFhYxZ04BP7yhENfQfOY+DSkj8pnwLVjxWgEZKfCzu/IZOwwyU4BL5x7zsfydyTEG/rXfFvzb6q0DmhIHlw+xHcYTuihglDbB/QdaGtn2jYPCHjClC9ufiMjxaY+OBMQNJbDssJ2q/2EnO1cfS0f36qwphfzV0CcJ/nh2YMfQUR4PvLcOXllu94u44m1rgn8bAa4AzHZ1tYOV8OUu+HInPPn7uTR5Yjmj1czNSX3gG3nw5t+KcMYGv6O7MbByPzzTKuAkxNpw892TILWLmmB6DPy10rZwqPTYgpU/zoDbsiBJhf9EAk6bkTtIQSd4zv0KihvhmT4wIQjtHzpy+qqqEb7/D/vnP0/tun067dlXZmd3dni7k7ji4fSTYdJou9QVjoyxm6u377HhZvMuG3RaS3LByQNg5CAYMQBSu6jVh8c7g7N4K2z29oZqDjiXDQ5Mx/qO2lwP+aWwxls3arQL7u8Bo7rBiTuRUNFmZAkpY2Cfd29CnyC8ou5IB2iwr+ZHZcKGMnh/r70AhkpuJsz+DvzzM9vKoLTcdtf+cL09nTV2GAzvZ1sZhGofT00t7D4AO/bZQPbVPqg6fPTHxMTAwFwY2hdGDoSBvbrm6H4zt4F/7rWnqLZX2ftcsXBxCAJOrQceL4M/eo+MJzlsp/Gr0yEuyvZiiUQTBZ0o0tkmiv6q9EC9d16wZ4D3RpyoA3RbU3vboLOsJLRBB2xImDzG9mPavMuGnnXb7B6e5pNJ6ckwrB8M6w8DekJOOsQH8H+lx2OPvh+stJWF9xy0tW32HILKmvbH3C/Hdgof1s8uTSWEYMmt0WO/h3/bBju940z0zuBcOgjSu7jg+QfezcY7vYH+m8mQnxOcYC8igaWgE0U6U6I/EPZ7f/lnxNgO0IFyojYMrd9vNqkXLNgIWyqhuNqexAq1GAcM729vFTWwchNsKrZLRBU1sOoLewO77yUzDXqk29CTkmiDRqLTvnU57ce4jQ0xHo9tUllbb2doaursreowlFXZpajjtcXISrMBa1Ave+vXw786QP6G7cON8OpO23TzoHeFMiUeZg6EGYO6bg9Os5JG+OUBeMsbtnrFQUEOnBcGP1ci0jEKOlGkMz2iAuGg9xSxg5YAABxnSURBVJRwdoBnczraS6i1dBec0cMeOf7bNrhzTGDH5K/0ZDhvnL01NHn3w+yEzbth3yGobYBDlfb2xc7AfM2YGHvyKzvNVnDune29ZbU/W+NPWOls2D5Ub8PNq8VQ4w3O2S64eBB8uz8kdXHAaTCwqBweOwS1BmKxS1S3Z0OKNhuLRBQFnSjTXhPFYIYcsBcCgOQAXwCO98r/eP+eK4bYoPP2bvjOYBgYpse7nXEtMz1g9zpV18L+cns7WAGH66Gu4eibw2FniWJjbIiJcUCiy87+JCd4b4n2iHd2OmQk24/rKH9mBn0N25srYMkOeH8PNHl/jvonw3dOgml9uq4OTmsfHYa5pS2Vjccn2NOEw9UfWCQi6dRVlOpoj6hAeL0abtsLZyTAs/2C+qU67P6PYfk+mJgL+WNDPZrIc6y9UR0NzccrCeD22O/Nkq/g87KWzxmRYY+IT+hpw1tXK26EX7VapsqOhbuyYWZq9BV+FIk0Ol7eQd0l6PjSIyoQnq+0vX2mJMFTfYL2ZXxSXA03fQAe4IEzYEx2qEcUefz9OWobtg/Vwxs74fXilv03cQ6Y3BsuGQjDMoL0DzmBwx74Qxk8VW6XrGKAK9PhjixIU2VjkbDgz/Vbq81RpvUr7/r6egoLCykoKKCoqChoX7M5KYfTi94BKXC+d0lo/qdQEdxJraiUn59/JOQ4nU6fQk7bkgDTbyriR+/CM5ttyEl3wg+GwNNT4eenhibkGAOvVMH5xbCgzIacMxPh7/1tp3GFHJHooD06UaQzp5QCwelNOA1hNjd4/cnw2SF7PPnBdXDf+NAsiUSqjtYvau/zCgoK+MEdhWRels/ShUUsXVDA0CqYcVM+Fw2ESbkQH8Ig8XEtzDsAa70BuF+c7U31zWQtU4lEGwWdKNKZU0qB4ArToJMYB3efBncstz2w/rYNvjck1KOKDL7WLwJbvfjme4r4w68KGH51IeXn5VNeBSOuymdIGrzxeAEZw2BaEJdRT+SrRnjwALzh3YeT5IAbMuH6DEjQ/LZIVFLQiSKdPaXkr+agUxtmQQdgcCrMGgm//cx2uE6NhwsHhHpU4c3XmcFdNfaE2zu7YflON0OvKuSkH+STlwbn9YOpfSD1/HyKegUvbJ9IhdsuT/25HBqxa/bfTYPbs6CnfguKRDX9Fxe/ZXmXIA42hXYcx3J+P9sA8tVi+N0GKG+wR9DbW6IIVXXpcNKRmcGqBvhgrw04G8tbPua0H81lWl+Y3g+GtNkvGMywfSz1HnimwoacCm/hxLMS4Z4cHRcX6S4UdMRvud6folK3rYUSbn1/HA64aaStsPvcVvifzTbs/HTE1/fshKq6dDg5VpCraoQJP87ngz1w5Tu2OjPY53BcDnyzrz0a7gyDTbweA3+vhocOQok3gA91ws+zYWqS9uGIdCcKOuK3nFhbOdYNHHDbMvnhxuGAa4dBphP+sBH+/hWU18Mdp0BCq/GGqrp0uKposF3D/7kXPjnQUtQP4KRUOKevXZrKCpPZEWNg2WEbcDY12PtyY23zzUtTIVYBR6TbCcNLkkSaWIed1Slpgl2N4Rl0ml08CNKc8Jt1dullWxXMGgHjerR8TCiqS4eTA3Xw0T5b1G/9ITs70mxwKkzuBZN6Q7/k0I2xPZ/UwfwDsKrOvp8SA7My4Zp0SNRGY5FuSwUDJSBuKLGvpAty4OoQFX7zxfqD8OtPWwrXTehpj6P3bXXx7srq0qHk9sCmClhdCmtKbUPU1k5KhTNzbbgJhyapbW1psDM4S70nqZwO25fqp5mQGQbLaCLiP3+u32H82lsiyQiXDTobG0I9ko45JRv+MBn+sgVe/souz6wphZmD7Ebl3zzQuRoykcAY2H3Yhr21B+2SVHWbjeQjMuCsXraFRu+k0IzzRLY1wO8PwSvVtmhlDHBZKtyWBb27uAmoiIQvBR0JiJHePRobImjiIzkebhgBF/SHJzZ6a+1sh4ceKGLDnwq48xeF/Ob+jtWQCWfGwP5a+KzMBpu1B1pmspqlxsPYHBjfw77NDJM9N+0pbrQBZ0mVbfEBcF4y3JFtNxyLiLSmoCMBcar3wripHsrdkBFBSwb9U6BwPKwshZvusSFn6FWFbDwrn//4CM7/UT75nuBWlw6kBjdsr4KNZfB5uX3bNtjEOWBEJpySZU9MDcsI/426OxthYZntrdY8AXVuMtyaBaPCOJiJSGgp6EhA9I6HPKfdL/HPw/Dt1FCPyDcOh92n8+1+bk6/q5B+38vnX/ttjZiN5RA3IZ9J/w6flrrZVA5D08MjGFQ02Aam2yptraCtlfZ9d5udd7EOW9fm1Gw4LduGHFeEhNHtDbbp5pIqe7IPYHISzM6CMQkhHZqIRABtRpaAeeCA7QB9aSr8OjfUo/HfoXpbEO+tXbC75ui/S4mDkzNhUAoMSoWBqfYUUqDDg8fYY/D7amF/HZTW2krEu2pgVzVUNrb/eWnxMDwDRmbCyAwYmgEJERJsmn1Zbwv9vVbdskQ1KRFuzoLxiSEdmoh0MX+u3wo6EjAfHYZrSmyl5A8HQXwYzHgEgjGw57DdtPvJQfj0INS0UwXaga0nk+mCDJet2ZPpsj234mPAGWPfxsfYANPoabk1eOxjVjbYWZqqRvu2tPbo2jWbn5mLIyaWvCtbls96JtqTUev+p4iMODfziuaSkxC5RfHW18HjZfCPVuHynCS4KQtO1QyOSLekU1cSFsYn2uKBB9z2InVhGB5F7gyHA/ok29u3B9rj2JsrYWsF7KiGHVV2uaiq0e6Fabsfxl8xQHaCDTRN6bEsXVDA5F5wb34+fZPtTE1RURHPP2qLGvaI0NmOj2vhsTJ4/7B93wGcnww3ZrVsdhcR8ZWCjgRMvAMuT7PLDYsrOhd0IqHXVGwMnJxhb82MsTMwB+qgrN57a7Bva5uOnrlp9Ng9M/ExENdqpicl3i45pTptUcO0eOiRYENOXHPBu3/Lp6i33Rg9KC3yKzd7DLx7GJ4qayn0FwvMSLXF/oZEwCmqSPiZFenOFHQkoL6fZjeOLq+1m0gH+3ihitReUw6HXa7K6IKZh2io3FzvgeerYFE57PDuM4oHLkuDGzJhYATVwYnUn1mRbsN0IxUVFQYwFRUVoR5KVLt+tzF5m425a2/nPr+wsNAAprCwsN33xXI6nQYwTqcz1EPpsINNxjx20Jh/22Z/RvI2GzN2qzG/LjVmT2OoR9d5+pkVCS5/rt8KOhJwH9faC9jQzcZsrOvcYzRfKJov5qG8YMyZM+eYX7+wsNDMmTOnawdkwuv56Ygv64z5xT5jRm1pCThTthuzqMyYaneoRxcYkfY9EYkkCjodpKDTdW7dYy9mP97d+ccI5IyFP2HlWK/OQ/WqvStmDwIR7po8xrxdbcyPdrWEm7zNxlxSbMxLlcY0eAI23LARibNsIpFAQaeDFHS6zo4GY0Z4L2zvVvv++YF+dexvWAmXpYmuCl3+fJ2yJmMWHjLm7O0t4WbYZmNuKjFm1WFjPFEYcIzRjI5IMCnodJCCTtf6r1J7kZu4zZhSH/ZfBCtU+Pu44XAh68plNF+fr8+9y1OjWy1Pjd9qzLxSY4obAjassBQuQVgkWinodJCCTtc67Dbmwq9alrDcHXglH+wZC3/DSndbmjjR81XtNua5cmMuKz56eWrGV8b8rcKY2ijZf3M84ba0KRKN/Ll+qzKyBNXmerhsF9QZ+H/Z9ujw8XRFTRKXy0VDQwNOp5P6+o5X92s+Lux0OiPySHdntX2+jIF19fC3Svh7FdR4f4PEY5tsXp0Bp0dwZWZfqY6OSPD5df0OeOwKY5rRCY3F5fZV/smbjXm/E/t1AqmzMzrddWmi7fN16b2F5ltfHT17c+4OY544ZMyBCD4eLiLhTUtXHaSgExoejzF3ek9hnbLFmDWHQzOOzoaVUC5NhPJoe/O/77r8QjN7jzE9Ztv3s24vNKO22O/pipro3VwsIuHDn+u3KiNL0DkcMC8Xyj22j9ENe+DpPnBKFzZobK9NQusKw63fb8vtdre7TNX8vtvtDtawQ1J11xi4dW4RjxUW0O+OQt67Jh+qIf3mfLJi4IuHCrgxE+6fE/3LdiIS+RR0pEs4HfD7XvDjElhTB1fvhoW9YUJS13x9f8LK8fZXBHuPTnthLFi9rbY1wCtV8Go1rCpzk3V7IQk35ZMZY/uWfTcNRj+Yz/0ZwQ13IiKBpM3I0qWqPTBrD/yr1qbse3Lg6vTus3G1s4K1EXpvE7xWBS9Xw4ZW+7JdDvhmsm2uOTnJBlURkVDx5/qtoCNdrs4Dd++3MwcAM1Lg/p6QFHP8z+vuOntarK2djbC0GpbW2Nm15l8AccCkJPh2CpybAqn6fohImPDn+q2lK+lyCTHwcC6cmgC/OgB/r4YvGuA3uXByF3T/jkRFRUVHQk5DQwNFRUUdntExBjY2wDs1Ntx83iYjjU+Ai1Lt8lRWbBAGLyISQgo6EhIOB/w4A0a54Pa98GUDXLoTrsuAWVmQotmEI9ruyWl+H469R6jewIrDNtgsq4F9rbbUxABnJML0ZPhmCvTWbwERiWL6FSchdUYivNwf5pTai/LCcnihCu7MhstSIaab7w3x5bTYrkZ7qu2Dw/DR4ZZCfgCJDpiYaIPNucmauRGR7kNBR/zmb2XYHnHwWC94uwbmHYTiRrhnP/y5HH6aCeenQFw3DTzHOy1W54Eva90UlcKHh2Fb49GfmxtrQ825yTAhEVyaJRORbkhBR/wWiFovDoedbZicDM+Uw2Nldl/J7H3Q9yBcm2GPN3e3DbKtA2JpE3xcBytr7e2LH+ZjgBUV9u9jgdMSYEqSPSk10qUZMRERBR3xWyBrvbgccF0mXJpmA88zFbC7Cf7rADx60B53/lYKnJ4Y3bM8xsBXjbC6DlbV2rfFjV//uCHxdvnvrCQ4MxFStSQlInIUHS+XgAlGrZc6Dyypgv8uP3ppJjPGzgCdnwwTQ1znxd+lO2NsPZv19d5bHXxWDxWeoz/OAQx12lNSExJtwMnRSxUR6QZUR6eDFHSCL1C1XtoyBlbU2qPo/6iGslYhIMkB30iAsQkwLtH+uStPbR1r9qrt/cbAfjfsaIQtDfBFPWxusCfOKj1ff1ynA8a4YHyiDTdjEzRjIyLdU7eoo/PLX/6SV199lbVr1+J0OikvLw/1kKQNf2q9nIjDYWduJiZBYQ+7nPNmjS18t98Ny2vtjTJ7fHqoE/KcMMT79qR4GBhva/gEWuuluzoP/PDufObfX8TT9xdw3t2F7Ls+nxnFdunp8DFeVsR6xzwmAUa7bMAZ6lJFYhERf0XMjM6cOXPIyMhg165dPPXUU50KOprRCZ5j1XoJdD+mtjzGzo6srrMbddfUwq6mY398Wgz0iIXcOHvaKyvGVmROirEzQ0kxEO8NF83/MQzQaOCwB2q9bw977NLSIXfLbeOjRex/pADindDYQNbthWTdcvS/PQboG2cD2HCnDTPDnTDYafcniYjI13Wrpaunn36a2bNnK+iEkY4u3XSVvU22+u/WBhuCtjXA1kaoamd5KNC2jHRBYwMxTiezdtTTKw76xcMg761vvGZpRER81S2Wrjqjvr7+qH0ilZWVIRxN9PKnM3gw9Iqzt3OSW+4zxgad/W7Y3wSlbtjXBBVuu5zUPEtTY8DdTvSP9872JLZ6mx5rC+813555oIiHGluW7nr9MXBLdyIi0jlRHXTmzZvHfffdF+phRL3jnSgKlwu9wwFpsfaW5wz84xcVFfHQfb61aRARkeALafm1uXPn4nA4jntbvXp1px//nnvuoaKi4sht586dARy9iHWsNg2FhYUUFBRQVFQU4hGKiHRfIZ3RueWWW7jiiiuO+zGDBg3q9OO7XC5cLrXDluAKt6U7ERFpEdKgk5OTQ05OTiiHIOK3SFi6C0f+FloUEemIiOkcVFxczNq1aykuLsbtdrN27VrWrl1LdXV1qIcmIp3Q3COt7dJe81JgbKyqI4qI/yJmM3JBQQF/+tOfjrx/2mmnAfDuu+8yderUEI1KRDorkD3SRESOJeLq6PhDdXREwk8weqSJSHTpVgUD/aGgIxKegtUjTUSigz/X74jZoyMi0am9HmkiIoGioCMiIdN6T059fb1qD4lIwEXMZmQRiS7HKrQIqKq0iASMgo6IhIQKLYpIV9BmZBEREQlr2owsIiIi0g4FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUUtARERGRqKWgIyIiIlFLQUdERESiloKOiIiIRC0FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUUtARERGRqKWgIyIiIlFLQUdERESiloKOiIiIRC0FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUUtARERGRqKWgIyIiIlFLQUdERESiloKOiIiIRC0FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUUtARERGRqKWgIyIiIlFLQUdERESiloKOiIiIRC0FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUUtARERGRqKWgIyIiIlFLQUdERESiloKOiIiIRC0FHREREYlaCjoiIiIStRR0REREJGop6IiIiEjUioigs2PHDq677joGDx5MYmIiQ4YMYc6cOTQ0NIR6aCIiIhLG4kI9gI7YtGkTHo+HhQsXkpeXx2effcYNN9xATU0NDz74YKiHJyIiImHKYYwxoR5EZ8yfP58FCxawbdu2Dn9OZWUl6enpVFRUkJaWFsTRiYiISKD4c/2OiBmd9lRUVJCVlXXcj6mvr6e+vv6ozwH7hImIiEhkaL5ud2puxkSgLVu2mLS0NPPkk08e9+PmzJljAN1000033XTTLQpuW7du9TkzhHTpau7cudx3333H/ZhVq1Yxfvz4I++XlJRw9tlnc/bZZ/PHP/7xuJ/bdkanvLycgQMHUlxcTHp6un+D78YqKyvp378/O3fu1BKgn/RcBo6ey8DQ8xg4ei4Dp6KiggEDBlBWVkZGRoZPnxvSpatbbrmFK6644rgfM2jQoCN/LikpYdq0aUycOJEnnnjihI/vcrlwuVxfuz89PV0/dAGQlpam5zFA9FwGjp7LwNDzGDh6LgMnJsb3w+IhDTo5OTnk5OR06GN3797NtGnTGDduHIsWLerUP1ZERES6l4jYjFxSUsLUqVMZMGAADz74IKWlpUf+rlevXiEcmYiIiISziAg6b731Flu2bGHLli3069fvqL/zZYuRy+Vizpw57S5nScfpeQwcPZeBo+cyMPQ8Bo6ey8Dx57mM2Do6IiIiIieijS4iIiIStRR0REREJGop6IiIiEjUUtARERGRqNVtg87FF1/MgAEDSEhIoHfv3lx99dWUlJSEelgRZ8eOHVx33XUMHjyYxMREhgwZwpw5c2hoaAj10CLOL3/5S84880ySkpJ8rvzZ3T3++OMMHjyYhIQExo0bxwcffBDqIUWk999/nxkzZtCnTx8cDgcvvfRSqIcUkebNm8fpp59OamoqPXv2ZObMmXzxxRehHlbEWbBgAWPGjDlScHHixIm8/vrrPj9Otw0606ZN4//+7//44osveP7559m6dSvf/e53Qz2siLNp0yY8Hg8LFy5kw4YNPPzww/zhD3/g3nvvDfXQIk5DQwOXX345N954Y6iHElGee+45Zs+ezS9+8Qs++eQTJk+ezIUXXkhxcXGohxZxampqOPXUU/n9738f6qFEtPfee4+bb76ZFStWsHTpUpqampg+fTo1NTWhHlpE6devHw888ACrV69m9erVnHPOOVxyySVs2LDBp8fR8XKvl19+mZkzZ1JfX098fHyohxPR5s+fz4IFC9i2bVuohxKRnn76aWbPnk15eXmohxIRJkyYwNixY1mwYMGR+0aMGMHMmTOZN29eCEcW2RwOBy+++CIzZ84M9VAiXmlpKT179uS9995jypQpoR5ORMvKymL+/Plcd911Hf6cbjuj09qhQ4d49tlnOfPMMxVyAqCiooKsrKxQD0O6gYaGBtasWcP06dOPun/69OksX748RKMSOVpFRQWAfi/6we12s3jxYmpqapg4caJPn9utg85dd91FcnIy2dnZFBcXs2TJklAPKeJt3bqV3/3ud8yaNSvUQ5Fu4MCBA7jdbnJzc4+6Pzc3l71794ZoVCItjDHceeedTJo0idGjR4d6OBFn/fr1pKSk4HK5mDVrFi+++CIjR4706TGiKujMnTsXh8Nx3Nvq1auPfPzPf/5zPvnkE9566y1iY2O55pprfGopEc18fS7B9iS74IILuPzyy7n++utDNPLw0pnnUXzncDiOet8Y87X7RELhlltuYd26dfzv//5vqIcSkYYPH87atWtZsWIFN954I9deey2ff/65T48REb2uOuqWW27hiiuuOO7HDBo06Mifm7unDxs2jBEjRtC/f39WrFjh87RYNPL1uSwpKWHatGlMnDiRJ554Isijixy+Po/im5ycHGJjY782e7N///6vzfKIdLVbb72Vl19+mffff/9rfRqlY5xOJ3l5eQCMHz+eVatW8eijj7Jw4cIOP0ZUBZ3m4NIZzTM59fX1gRxSxPLludy9ezfTpk1j3LhxLFq0iJiYqJoo9Is/P5NyYk6nk3HjxrF06VIuvfTSI/cvXbqUSy65JIQjk+7MGMOtt97Kiy++yLJlyxg8eHCohxQ1jDE+X6ejKuh01MqVK1m5ciWTJk0iMzOTbdu2UVBQwJAhQzSb46OSkhKmTp3KgAEDePDBByktLT3yd7169QrhyCJPcXExhw4dori4GLfbzdq1awHIy8sjJSUlxKMLX3feeSdXX30148ePPzKjWFxcrH1inVBdXc2WLVuOvL99+3bWrl1LVlYWAwYMCOHIIsvNN9/MX/7yF5YsWUJqauqRGcf09HQSExNDPLrIce+993LhhRfSv39/qqqqWLx4McuWLeONN97w7YFMN7Ru3Tozbdo0k5WVZVwulxk0aJCZNWuW2bVrV6iHFnEWLVpkgHZv4ptrr7223efx3XffDfXQwt5jjz1mBg4caJxOpxk7dqx57733Qj2kiPTuu++2+zN47bXXhnpoEeVYvxMXLVoU6qFFlJ/85CdH/l/36NHDnHvuueatt97y+XFUR0dERESiljZTiIiISNRS0BEREZGopaAjIiIiUUtBR0RERKKWgo6IiIhELQUdERERiVoKOiIiIhK1FHREREQkainoiIiISNRS0BEREZGopaAjIhHtjTfeYNKkSWRkZJCdnc1FF13E1q1bQz0sEQkTCjoiEtFqamq48847WbVqFW+//TYxMTFceumleDyeUA9NRMKAmnqKSFQpLS2lZ8+erF+/ntGjR4d6OCISYprREZGItnXrVq688kpOOukk0tLSGDx4MADFxcUhHpmIhIO4UA9ARMQfM2bMoH///jz55JP06dMHj8fD6NGjaWhoCPXQRCQMKOiISMQ6ePAgGzduZOHChUyePBmADz/8MMSjEpFwoqAjIhErMzOT7OxsnnjiCXr37k1xcTF33313qIclImFEe3REJGLFxMSwePFi1qxZw+jRo7njjjuYP39+qIclImFEp65EREQkamlGR0RERKKWgo6IiIhELQUdERERiVoKOiIiIhK1FHREREQkainoiIiISNRS0BEREZGopaAjIiIiUUtBR0RERKKWgo6IiIhELQUdERERiVoKOiIiIhK1/j9JIPmKA72TjwAAAABJRU5ErkJggg==",
"text/plain": [
"Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"p(x⋅∈S|m) ≈ 0.20987514863080037\n"
]
}
],
"source": [
"using HCubature, LinearAlgebra# Numerical integration package\n",
"# Maximum likelihood estimation of 2D Gaussian\n",
"μ = 1/N * sum(D,dims=2)[:,1]\n",
"D_min_μ = D - repeat(μ, 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}(undef,100,100); B = Matrix{Float64}(undef,100,100)\n",
"density = Matrix{Float64}(undef,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": "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\", read(f, String))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.1.0",
"language": "julia",
"name": "julia-1.1"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.1.0"
}
},
"nbformat": 4,
"nbformat_minor": 1
}