{ "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