{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"# Algorithme EM - Application à des données non simulées\n",
"\n",
"Dans ce TP, nous allons étudier des données issues de l'article suivant : \n",
"\n",
"- Bachrach LK, Hastie T, Wang M-C, Narasimhan B, Marcus R. *Bone Mineral Acquisition in Healthy Asian, Hispanic, Black and Caucasian Youth. A Longitudinal Study.* J Clin Endocrinol Metab (1999) 84, 4702-12. \n",
"\n",
"Ces données son disponibles [ici](http://www.math-info.univ-paris5.fr/~jdelon/enseignement/densitesOs.txt). Récupérez les et placez les dans le même répertoire que votre notebook. \n",
"Elles représentent des mesures relatives de densité minérale osseuse spinale sur des adolescents nord-américains. Chaque valeur est la différence de mesures prises sur deux visites consécutives, divisées par la moyenne.\n",
"\n",
"L'idée est de savoir si la population peut être décrite par un mélange de deux gaussiennes. \n",
"Le but est donc d'appliquer à ces données l'algorithme EM, vérifier que le $K$ optimal de la méthode de sélection de modèle vue en cours est bien $K=2$ et proposer un clustering de ces observations. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as sps\n",
"import scipy.linalg as spl"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On commence par créer le tableau de donnees \"data\" à partir du fichier."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"data=np.loadtxt(\"../datasets/densitesOs.txt\")\n",
"n,d = data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Données"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**1) De quelle dimension sont les données ? Affichez les données sous la forme d'un nuage de points.**"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXu4XVV16H/j7JMEX1U4+ARiULBXqF6BiD3ftRoFIrVW6MUH1RoUJU0FW2xtNB8+KFgesbdNK972RBpK+vBJVdrSD9FylN4cJeGlIqJAkQRQMYiPqgeSjPvHXMs99zprrT3XY++19j7j933723uv51hrzTXGnGPMOaaoKoZhGIYRM9G0AIZhGEa7MMNgGIZh9GCGwTAMw+jBDINhGIbRgxkGwzAMowczDIZhGEYPZhgMIwURWS4iPxGRTgPnvltEjq/pWLMi8pY6jmUsHswwGAuoUzH1Oc9SEfmGiOxKLH+piNwoIj8SkbtEZK23bpWI7IuUdvw5rW7ZVPUeVX2squ4NuI4VIqIiMlm3HMNkWM895bwqIocN+7xGNiNdkI2R54+B7wGPjReIyBLgU8B6YDOwErhWRL6sqrdEm92nqgcPW1jDWCxYi8EohIi8QkRuFpGHRGSbiDzXW/dOEblXRH4sIreLyHE5xzkU+B3gwsSqA4BfAv5eHduB24AjSsqrIvL7Ucvj+yLyARGZiNZNiMi7ReTbIvI9EdkqIo+P1vW0AiKXzPki8v+i6/usiBwYneaL0fdDUQtmWkQOE5EviMgPo/N+LEfGN0Qy7BaRcxLr/k5E3u/9X5VsYSW2PyFqhf1QRC4BxFv3TBH5j+g83xeRfxSRJ0Tr/h5YDvxLdA3ro+WfEJHvRMf7oogcmXPuN0b3+cci8l8i8npv3ekicpuI/EBErhaRp0fL43t3S3Te12Yd3xgiqmof+/R8gLuB41OWH42r4b8A6ACnRdsuA34Z2Ak8Ldp2BfDMnHP8K/BbwCpgV2LdPwFnRueYjs55SLRuFfAw8F3gv4C/AB6Tcx4FrsUZnOXAN4G3ROtOB+4AnoFrtfwzziDF8iswGf2fBe4EngU8Kvp/Udq20bKPAOfgKl/7AS/MkO8I4CfAi6L7+OfAnvj+A38HvN/bfsH98tYdCPwIeBWwBHh7dKz4eg8DTojO80ScQduU99yje/S4aJ9NwM0Z535MdO5fjv4/FTgy+n1ydJ+fjfNSvBvYlnhGhzVd7u3T/ViLwSjCGcCMqn5ZVfeq6uXAPPCrwF6c8jhCRJao6t2qemfaQUTkt3BK9FMZ5/kI8N7o2NcB56jqzmjdN4Dn4RTPS4FjcMo0j4tV9UFVvQen3H47Wv564M9V9S5V/QmwATg1J1Zwmap+U1V/Bnw8kiOLR4Cn4wzlz1X1PzO2exXwr6r6RVWdB94D7OtzPVm8HPi6qn5SVR/BXet34pWqeoeqXqOq86r6AO6+vTjvgKq6RVV/HMl2LvA/41ZVCvuAXxGRR6nq/ap6a7T8d4ELVfU2Vd0DXAA8L241GO3DDINRhKcDfxS5kR4SkYeAQ3DK7w7gbJzy+J6IfFREnpY8gIg8BtgIvC3tBCLyP4CPAWuApcCRwHoR+Q0AVf2Oqn5dVfep6n/hYhGv6iP3Tu/3t4FYrqdF//11k8CTM47zHe/3T/FiIymsx7lxrheRW0Xk9IztnubLp6r/DezOOW4eyWOp/19EnhQ9l3tF5EfAP+BaGamISEdELhKRO6Pt745WLdgnkvu1wDrgfhH5t+hZgis3f+mVmQdx9+agktdpDBgzDEYRdgJ/qqpP8D6PVtWPAKjqP6nqC3GKQIGLU45xOM71cp2IfAfnvnlq5MdeAfwKcLuqXh0p/9uBfwN+PUMmxfOjZ3CI93s5cF/0+75IVn/dHpybqggLUhRHBuwMVX0arsb8fzN63tzvyycijwamvPX/DTza+/+UHDmSxxJ6r/3CSNbnquov4WI8/r1LXsfrgJOA44HH454bZNzv6JmdgGvNfQP4cLRqJ/C7iXLzKFXdlnMtRoOYYTCyWCIi+3mfSdyLvk5EXiCOx4jIb4jI40Tkl8V1M10G/Bz4Gc69lORrOGX1vOjzFpwifh5OgdwEHB4dS0TkmcArgFvgF8HX5dG6Q4CLgM/0uZY/FpH9o+3/ANciAeeyeruIHCoij8W5OD4WuTuK8ADOjfKMeIGIvFpE4p5TP8Ap3bT78UngFSLyQhFZCpxH73t5M/ByETlARJ6Ca5Vl8W/AkSLyv6Pn9fv0GpLH4eIZD4nIQbheYT7f9a8h2n4e14J5NO7+pCIiTxaRV0YtwvnoPPH1/g2wIQ5ci8jjReTVOec1mqbpIId92vfBuQw08Xl/tO5EYDvwEK6G+gmcAnkucD3wY5yr4F+JAtF9zrWKhcHn1+AMyI+BXbiWx0S07g+Be3GunJ3AB4HH5RxfcQryLpyC+z9AJ1o3gYtl7MQp938A9o/WrWBh8Pkt3nHfCPyn9/+86BgP4WIuGyM5f4ILWq/NkfE04J5IvnPwgsC4wPXHcIHdr+ACyqnBZ+/5fBP4IXAJ8AW6wecjgRsimW4G/sg/Fq51cE90De/Auco+Ez2Hb+Pce6mBYlwr4QvReR+K7tcR3vo3AF+NrmMnsMVbty4qSw8Br2m6/NtHkejBGMZYIiIKHK4uBmIYRgDmSjIMwzB6MMNgGIZh9GCuJMMwDKMHazEYhmEYPYxkEr0DDzxQV6xY0bQYhmEYI8UNN9zwfVV9Yr/tRtIwrFixgh07djQthmEYxkghIt/uv5W5kgzDMIwEZhgMwzCMHswwGIZhGD2YYTAMwzB6MMNgGIZh9GCGwTAMw+jBDENbmJuDCy9034ZhGA0ykuMYxo65OTjuOHj4YVi6FD7/eZiebloqwzAWKdZiaAOzs84o7N3rvmdnm5bIMIxFjBmGNrBqlWspdDrue9WqpiUyDGMRY66kNjA97dxHs7POKJgbyTCMBjHD0Bamp80gGIbRCsyVZBiGYfRghsEYHNYF1zBGEnMlGYPBuuAaxshiLQZjMFgXXMMYWcwwGIPBuuAaxshiriRjMFgXXMMYWcwwGIPDuuAaxkhSiytJRE4UkdtF5A4ReVfK+j8Uka+LyFdE5PMi8nRv3Wki8q3oc1od8hiGYRjlqWwYRKQDfAj4deAI4LdF5IjEZjcBK1X1ucAngY3RvgcA7wNeABwLvE9E9q8qk2EYhlGeOloMxwJ3qOpdqvow8FHgJH8DVb1WVX8a/f0ScHD0+2XANar6oKr+ALgGOLEGmQzDMIyS1GEYDgJ2ev93RcuyeDPw70X3FZG1IrJDRHY88MADFcQ1DMMw8qjDMEjKMk3dUOR3gJXAB4ruq6qbVXWlqq584hOfWEpQwzAMoz91GIZdwCHe/4OB+5IbicjxwDnAK1V1vsi+xgCwdBWGYWRQR3fV7cDhInIocC9wKvA6fwMROQqYAU5U1e95q64GLvACzquBDTXIZORh6SoMw8ihcotBVfcAZ+GU/G3Ax1X1VhE5T0ReGW32AeCxwCdE5GYRuTLa90HgfJxx2Q6cFy0zBomlqzAMI4daBrip6lXAVYll7/V+H5+z7xZgSx1yjBVzc4MbNRynq4hbDJauwjAMDxv53EYG7eqxdBWGYeRghqGNpLl66lbelq7CMIwMLLtqG7HMpIZhNIi1GNqIuXoMw2gQMwxtxVw9hmE0hLmSDMMwjB7MMBiGYRg9mGEwFmLpMgxjUWMxBqOXKmMoBjkozzCMoWGGYZwpo6jLjqGw/EuGMTaYYRhXyirqsukyhjEozzCMoWCGYVwpq6jLjqGw/EuGMTaYYRglslxDacuLKurkMYrW9m1QnmGMDWYYRoUs11DW8iKKuq74gA3KM4yxwLqrtpVkl9GsORTy5laYnoYNG/ora5ufwTAMD2sxVGFQ3TPTavBZrqE6fPsWHzAMw8MMQ1kG2T0zrQa/YUO6a6gO377FBwzD8DDDUJZBds/MqsFn+fDr8O1bfMAwjAgzDGUZpPvFavCGYTSIGYayDFp5V6nB1xn7sDQXhrHoMMNQhTa6X8rEPvLGR1iaC8NYdFh31XGjaNfTWPm/5z3u28+oat1YDWNRYoZh3Cg6X3Se8re5pw1jUWKupHGjaOwjL4he9FgWjzCMsUBUtWkZCrNy5UrdsWNH02KMD3UodItHGEbrEZEbVHVlv+2sxWDUE0S3tNuGMTbUEmMQkRNF5HYRuUNE3pWy/kUicqOI7BGRVyXW7RWRm6PPlXXIYzSAxSMMY2yo3GIQkQ7wIeAEYBewXUSuVNWve5vdA7wReEfKIX6mqs+rKocxBPJcTnWN67A4hWE0Th2upGOBO1T1LgAR+ShwEvALw6Cqd0fr9tVwPqMJkjGETZtg9+6FeZuqKHOLUxhGK6jDlXQQsNP7vytaFsp+IrJDRL4kIidnbSQia6PtdjzwwANlZV28JNN4F8WPIczPw5lnpo99qIKNmzCMVlBHi0FSlhXp6rRcVe8TkWcA/yEiX1XVOxccUHUzsBlcr6Ryoi5S6qiJ+91aJyac8t63r95As6X/NoxWUIdh2AUc4v0/GLgvdGdVvS/6vktEZoGjgAWGwSiI76sv0mMoy8fvxxCmpuDss+tX4JY80DBaQR2GYTtwuIgcCtwLnAq8LmRHEdkf+KmqzovIgcD/AjbWINPiJi0eEFIT79ey8GMIz3nOYBR4G/NPGcYio7JhUNU9InIWcDXQAbao6q0ich6wQ1WvFJHnA58C9gd+U0T+RFWPBJ4NzERB6QngokRvptGhTb1pki2E3bvDauJFWhamwA1jbKllgJuqXgVclVj2Xu/3dpyLKbnfNuA5dcjQKEV8+MMwILGvfn4eRJzrJ0SRm4/fMAwsiV49hPamyctkWpa03kbT08591Om4APHZZ4edK/bxn3++dRU1jEWMGYY6CB31W3d3zDxDs3u3Mwp+zyHDMIwALFdSHYT2pqnbVbN1K/z856C6MCZQ5lyhLrE2xVOyGAUZDaOlmGGoixAffp1pI7ZuhUsvdUYBYHKyWspsCAs+j8Lo5FGQ0TBajBmGYVNX2oi4pQAuwPymN6XnL4KuGyn+n1WbDmlltDGLavJ62iijYYwQZhhGjVjp+UZhv/1gzZqF26bVnCG7Nh3Symhbz6W0a2ybjIYxYphhGDV8pTc56VoKa9akK/GsYHdWbTrEL9+20clp17hhQ7tkNIwRwwzDqFFEMWfVnNOWFfHLt2lwW9Y1tklGwxgxzDC0mby8RSFKL8uIpC0bVb9821owhjEG2JzPbSDNAAy7Z4315DGMscfmfB4VshTysGvwbal52/gDw2gcMwx1UVahZRmAJnrWNO2Xt1aLYbQCMwx1UGXEcF7wdLHNoTyqcQ7DGDPMMNRBlRHDeQZgsc2hbOMPDKMVmGGog9ARw/PzLqnd/Hyv8RiUC2fUauBtiXMYxiLHDEMdhCi0qSlnFMB9T02lHyvP9ZPVeylr+yo18KZcUE3HOQzDMMNQG2kKzVeuu3fDxIQzChMT7n+SPNdPvG5+3qX3vuQSN71mv6k4y9TAR80FZRhGrZhhGBRp8y4vW1Y+QZ3vitq3D848E97ylu728/Nw7rnuUzZOERuye+4ZLReUYRi1YoYhjyrulKSSD5l3Oc/1s2pVd0Y2cMe98UaXL0nVLf/c5+C668rV8H1DNjnpzgUWBDaMRYgZhiyqulPSlHxW7d03QHk9lC65xLUU9u51xmDHDqfEV650v/3Z2kJk9c/rGzKAM86A5cstCGwYixAzDFlU7dET6t9PM0AbNqRvu3atiyuce65rHezb5+Q7+mj46lerzda2aVOvIcvK2OrvP+jg9CiNwTCMMcIMQxZFevRUSXZX1ABNTzvDcN11vUp8zZpqs7WFuLr86x10cNoC4IbRGGYYsqhS46/qcior26DPG5M313RdNDUGw1ophmGGIZdB1PjTzlGmS2nV/v7J80J4Wo/LLuvOINfpDCY4PexR0PE82lu2uGdprRRjEWOGoSqhCiyvJtrUoC7/vBdeGGbgZmdhzx73WwROP30wsvuGa2pq4bzVdZI2j7Z10zUWMWYYqhJS42+zvzw2WFNTYQYuaQjT5pquKkuyB9eg713c6vPn0bZuusYiphbDICInAn8JdIBLVfWixPoXAZuA5wKnquonvXWnAe+O/r5fVS+vQ6ah0q/G38acRXNzsHEj/Mu/OIW4bJnrmbR7dzNzPs/NueM98ggsWdK9R8O4dyHzaFvswVhEVDYMItIBPgScAOwCtovIlar6dW+ze4A3Au9I7HsA8D5gJaDADdG+P6gqVzDDeOHbljU0VsIPP9xdNj/vjEJWV1mfQbi+tm7tyvPww+7/sOal6Gfs2tziM4wBUEeL4VjgDlW9C0BEPgqcBPzCMKjq3dG6fYl9XwZco6oPRuuvAU4EPlKDXP3xX3gRNx7gzW924wXqpG1ZQ2dnXc3cZ2KieYOVxrBiDXnGro0tPsMYIHUYhoOAnd7/XcALKux7UNqGIrIWWAuwfPny4lKmkRzte/317gODMQ5V51aoYlj8/Vetcu6auIbe6cCHPtSssluzxvV2SotdDCvWkEXbWnyGMWDqMAySskzr3ldVNwObAVauXBl6/HziF/5nP+tdfsUV9RsGKK/cq7oy0vafnXXuGug/yrkMRa91ehquvTZ7nyZr7W1r8RnGgKnDMOwCDvH+HwzcV2DfVYl9Z2uQKYz4hd+4ET796e7yU06p/1y+cu50XDfPUIVcVSmm7b9hw2BTWZQxZHmtqqbnlrB5IoxFxEQNx9gOHC4ih4rIUuBU4MrAfa8GVovI/iKyP7A6WjY8pqfhU5+CmRlYvdp9l2ktzM25sQBzc+nrk8p5ZsYpz6ztfWKl2OmUc2VU3b8oaYaoKrERP//8Yi2m2Ei95z3h99swFjmVWwyqukdEzsIp9A6wRVVvFZHzgB2qeqWIPB/4FLA/8Jsi8ieqeqSqPigi5+OMC8B5cSB66KxdW959FFJDjpVzPIiqSDqJNFdGkVrwsF0hg/LJl6m1W+DYMApTyzgGVb0KuCqx7L3e7+04N1HavluALXXI0RghyidWzlu3uiDrnj3FlKavFIu4anwDEtIVtQ6a9sknA+0WODaMQtjI5zooonyWL4e/+qv+A8nyCK0FN9n/vimffNo1W+DYMAphhqEOhp0WI9QQVXWjjOJo36xAe7wORudaDKMhzDDUxTDTYoS6aqr25BnF0b5p1zyq12IYDWGGYVjU7esOcdVU8fW3JWhbZjxE8ppDM8cahgGYYRgeg0w+l3fMsr7+NgRt6xoP0YZrMYwRwgzDMKk7IBuqOP3U2qFB76Z7FkF9rZY2XIthjBBmGEaZEMUZG4/5edi3zyXLW7YsrPbd9GjfOmv6TV+LYYwQdYx8NpoiZERzbDz2RYlt9+2rbzRyWfqNEo+ZnnZzRBx3nPs2xW4YQ8FaDEVoW/fNPBdJcmY2v8XQpJ+96OC8s8922153HTznOe247yG0rawYRgHMMITS1i6PaS6SpKzxzGxpMYZhK7AicYO29IwqSlvLimEEYoYhlFFSUklZs2Zmq1uBpRmZ5LIicYNR7U20dWs3J1bby4phpGCGIZQ2K6myytc3IPPzcO657lN2IqCkkYF0wxPaQ2gUexPNzblcWBpNGdLptKusGEYAZhhCaauSyqr1FxkZHccfPvc558sv03LwjczPf+5qzcuXp7eyivQQSiYPbNP9T5NndtYlSAQ3Xezpp7dDVsMogBmGUNqmlGKyXFxZsQf/GmIDcu65zij4PZaKXuOqVa52vHevqy1fdplLFlhXK6ttfvsseZKtNX+KUsMYEcwwhNA2peSTlRsozdefdg3T084wXHddNQU+Pe1qxzMzzjDs2eNiG5s2ualSTzml2j0rG+MZlEHPM8h1tyzbWikxxhYzDCG0OfDsz/MA8NWvdrt4+lOI5l1DEWWWp6TWrIHLL+8amIcegve9z52zanfTMjGeQRr0PHnqHEzX5kqJMbaYYYhJKrwmJnupUjOMFbKIcwnt2+cU8syMW7dpU/419FNmc3PO+GzZ4o6bpqR8A/PQQ/Bnf9YdWDc/Xz3td9Ga+CAN+rBiTm2ulBhjixkGSO/3H9e6YwU4aCWweTOcdZZTAKEpK2J85TEx4VoK8fShcZfJm26C005z269ZU1xBH3dctwsmdAPMyePE/1/84q5RgOK9c/JcX6GEGvSyBnkYaTba3BvOGFvMMMDCWtkVV6RP9jIoJTA3B2ee2e3NEteuY9n6Kayk8ti0yRmCeArRTqd3OtGiAdH4/sRGAdzvLVvSjYzfMyfm7W8Pv39zcy7uEfeWKltTDqnVt91V09becMZYY4YBFirWU06pHowtwuzswtr11FS4wspSHnFs4frr4TOfKT/gyr8/4Axm/J12rFWrXMsl3k4EnvCEsHOlJf2r8gz61eqbcNWUmWPCDIIxRMwwQLpifc5zuv/BJX0bVI1t1SrnPpqfd4rwkktcj548hZXW9TTLrXPuud3a/uRkcSXr35+pqV43W9axfMOwZEn6dlnjAOKkfxMTcPzxTn6o/gzSzjdsV03bWyiGgRmGLknFGv+fm3PK4pFHnIKrUqPMqimmGaa5uWyFVUS5zM721tzf9KZy8vv3xzeaacfyW0BZg7yyRkrfc083RtLpuNYbVFemeTGLYU6gZMFkYxRQ1ZH7HHPMMTo01q3zw7jufxm2bVN91KNUOx33vW1b2D4XXLBw2wsucMcB933BBfnHWLpUVcR9h5w379yh+/a71uQ1rFvX3WfJEvc9MeGWrVsXfr1ZFLlnVcm7/jLlwDBqAtihATrWWgzDokxNMdmKSabSDnV/iPR+96OquyOkFp504UD3/sStjTgmAtXdPcN0GdU1ZsRYHLRwAKMZhn6sWeN69FRNcVBVMWWl0u5XmOIeQvFo5BCDFGLEqs41nVSQ0B2LMTnp5I3HS6xZ0w2kF3l5qo6DKEKRcS9NBpNbqIQWNS2NOdViGETkROAvgQ5wqapelFi/DNgKHAPsBl6rqneLyArgNuD2aNMvqeq6OmSqjelpuPba6i9TVcWUVNZZqbSTFDFIoS2SugpzUkEmDUXyXpUZe1FlHESVc7WxVdBSJbSoaWnMqbJhEJEO8CHgBGAXsF1ErlTVr3ubvRn4gaoeJiKnAhcDr43W3amqz6sqx0CpS6FUOU7ZFkeoQerXIvFrmoMqzGkdAMoyzBcu7VyDHPdSlpYqoUVNSwcw1tFiOBa4Q1XvAhCRjwInAb5hOAk4N/r9SeASkVCH95DIamK3pemdpeBD5AsxSHktkjSjUbUwD/q+Zr1wgzhvS1/uBZSVsy3vwDjS0phTHYbhIGCn938X8IKsbVR1j4j8EJiK1h0qIjcBPwLerarX1SBTMbKa2EXnJw55uFVesrRgdF2ugaTSmJrqjhtIMxpp3WtDr8uXe3LSdaE96qiwmEkoWV2A67hfWenLhxG/qHLsMnKa+2nwtHAAYx2GIa3mr4Hb3A8sV9XdInIM8GkROVJVf7TgJCJrgbUAy5cvryhygqwmdmjTO/Tlqfslq9s1EOdSOuqo3kFsaS0EvzAnx3p88INdJR/LmTaIbe/ebqI/VTegrWieqDySMtaRZiMvdjGs+EVV41Bkf3M/LUrqMAy7gEO8/wcD92Vss0tEJoHHAw9G/WrnAVT1BhG5E3gWsCN5ElXdDGwGWLlyZdLwVCOriR3a9A59ebK284O+RWrNdbkwksoH+rcQfLZu7XYrffhheOtb3e9Ox3WRjXM0JSeziZPyxaOyqyjskOurI83GsBVl04p5VNxkRq3UYRi2A4eLyKHAvcCpwOsS21wJnAbMAa8C/kNVVUSeiDMQe0XkGcDhwF01yFSMrCZ2aNPbf3k6HTd6d25u4fZpL1ma0pqc7M6jkBZLgO7v0MBy3jbJyeshu4UwN9ebmmJuDm68sfd4+/a5YyXHI8RKLb6vcRrvPXt6FbbvxipTo89qoSTTbJRRsMNWlE0r5pb6wI0BEzIKrt8HeDnwTeBO4Jxo2XnAK6Pf+wGfAO4ArgeeES0/BbgVuAW4EfjNkPMNdeRzKNu2uRG6S5a4UcbLlqWPak2OKPZH5PofETcydmbGHXfpUrfd0qXu2KEjZ/uNtI1HRsfnjeVOG/mcPNbMjPuemOjKPDnZHbW8dGn4/ZiZ6X6XGRkc3/+0e1P3aOP161UPO8x9D4Mqo9ANw4NhjnxW1auAqxLL3uv9/jnw6pT9rgCuqEOGgRLa82frVudnB9cCSJuvIElcI4xbDCJd8zA/79Jxx/MoQ3YtPIt+rojZ2YW5lOLl/fL8xOnJ45r4ypVw883dVN9nn+3iE7G8WffNd6fddFO+uy3tGaTNF5HWQqmj1rt5M2zc6H5v3AjPfCasXVv+eCG0JThpvZMWDTbyuR91Bv+yjuVnLvXnUYgzlMbKTsQFd32/fT/XQj9XRHL9UUdlX29y22R68qOPhhtu6Bq4m2/uyp+Votu/L/Pzbr+JCbc86W7LegaxwfLvU/Ja61KuV1yx8P+gDUMbsN5JiwozDP0oEvxbs8b5zOPeOcn0GVnHSiqtOP2Dn+I67tYZHzO05tavtpxcn3e9acdKpif353wOnddidrbbYoo544xujOXCCxfK5N+DZIwnjs+A27doUD+PU06Bz3629/9ioOkguDFcQvxNbfsMNcZQ1D+d5w+uM8PqoKjqj0/KGyL/tm0uNhE70SYmerOf+jItXap68sndmEssY9p5/fjHxISLP6xb141nlL2nMzOqq1e778VCSKyqjnLa1nhKW+UqCIExhsaVfJnP0IPPdRaKrGO1qeA1IcvMjAtU+wo8GfiOg8siXSOSlUI7K6jvGx9Le+0Ifd55ZbeO4H7acdrwXoxRqnQzDKPEsGpjdVO3wVy3Lr01EJNU9nHPrbzWWdxiSPsMel6GUaAOpVfXXBd5c3Q0qZCHOZfHgAk1DBZjgOZ7W6T5b+Plyak02xL0qzO9RDyeYe/ebsA9bbCbH0vwYy5p5/XjIf6c1zFVBrk1XV7qpI7YQV1jLZLHgXbENZoeS9IAZhgG1duiSFK+tDxFsUx5irJJ6lAoad1MNUqNkdWzqEgVqljhAAAXjklEQVS3U78r7NVXuwB3pwNvfzs84QnFFXtsxOJeY8M21GUMUr996lB6dXUHTh4HejszNKWQ67q+EapQmGEoo+D6PeC0JHFxL5mQeYd9mWJFCU5ZTk0tPF9VyhTYOhRKWjfTZcvc2Iebbkrfp1+301h5Q7c1EfJihz7TrLESg6ZMBSZknzqVet35raA9o66rXt+odfcN8Te17VNrjKFMr6Mi8xn7vvDQuYuT51i/vhuYTY7o9UcND2p+5rx94/OWiTf45/YDzsleSMlAdPLc/rJly7r3PWSO6zi2kRXXiCkS3xgERef5vuCCYnNlF31+bY17tZWWxCmwGEMgRWtMIS2MtCRxWTmIQmSanXWuJN+dBAtzLJXJTFrFJeS7asrUhtLufTILapx99fLLe9Ohv+Ql3fNde61bHl9LzCOPLBzzkPwd2grIGisxrFrf1FR38F9W2dm82bW2br/d/V+yxMkK+eWt6PMbtdrvICjayg5tYbfF3RRiPdr2abRXUlYNO60ffVpNtGrNOj5OWnfMMjWRNvVKiWv8fnfUtGtbt6533bHHdu9tssXg515K5pnya9QhrYCmasl+D6slS9LHT8zMLLxnIu4a+8lc9Pm1pPbbGGXfmX7lZwjdYrEWw4DIquVmxQ7SJrEvWhPIatUsXVo9lXQdPua6em1s3equJ+aII+DOO/un/9i+3d3/z3/etR78GIPfIkrmmYLwXk6Q72euo6aXdYz4GuJUI7t3L9w3maoD3HUedVT/lB1Fn9/UlCtvqouml04PZVvZ/eIULRpdboYhpsiLnXzA/dJIxL+zUmeXKVRpOZbKXl/VwFpdAcwkL3oRXHrpwuPGqUdi5R4r+tnZ9LmWfReQn2dqzRr3SQari1LWtZIsD6E5qtIUcTJVR0xouUjO8Z23/dlnu7IuAi97Wb+rHD8G1X21Td1iQ5oVbfvU7kqqIw1Ev/2TgdZ+wc46zz/oJmpdLpakOyQv5YQ/ErrfdeUFyZtypSXP2y9QHHKPZ2ZUn/3sXldcXgDeH23e79rTAtpNBOGbpo4OHyHHH9D9xFxJBajahAupMee5NJLnK+qWmJ3tupTm59NTa+etT56v7PzNVQORu3c7F0XsGktzmcTkueqytvX/x9TRfE8bh9JvoqHkeSG/thjSqlu71iU1/LVf66ZS37Mn/Zrm5lxK9z173P+0cuFv63e/7nS6kzFlleFB0lSAdhhB97q6/VZkcRsGf0rNOgb55D3Q5KhdjVJRJ8/n95cXgXe8Ay6+OP/cU1NdY7Nv38KxDnnrk4V906ZiI62LKNaQwVbLlg2uKR06uLDo7HFJt17e/csqc7Fbq6rCm53tHeHd6aTfx7inW3K7vBnwYmNzxhnuOx6tXvezynO5DkI5hxqaFsUABk5Is6Jtn1pcScmmvN80HFT2zJB+/xdcsLBXTtpMYf7+F1zQm0U06YbIW590g6xeXcwtEuqKKbJdaFO6iBsob1vfPTBIt9K2bc61I9LtMRWSebbo+IK4B9PkZHYZTtsur8ddSE+8sjJnydXpuI9/v1av7pblOnpE1VWGRgQsiV4fsl7gIn5un7p8g9u2LUz8NjGxUJGlTbFZJsZQ9Fih155cVrWLY9o5ihyzn/GsQ8ZYzqwKx8kn9z7XuJtt6LFCjWloWUxWVPKUbpFjVlGe/nNK67LsP8M6lHPRZ95Ul+WaMMPQj6wCvHp1b2Fcvbr8scqyfv1Cw5BXy/dbDmVq4mUVSxZp96OMkutXky9y30MMftGacd71J+VeunShwusXuA1RWnWUPb+WHpe3eN6LtFHnedRRAViyJN0wxC3piQn3XtahnMegFVAEMwwhpL3wRVsM/WpaZVm/3h0zrWbU9sKcpRxCjFfaYLS450zRmmzS3RYrFpFsJbtuXa8yrGIs/PuQNmivaGqUtHOV7RHly+63ZCYmXEvGV84hqUXi465b5/YVcc+vTNmcmXHuLf95LVlSb0++pNwj3AooQqhhWNzB57SAcTwY6IorXN/wvMFB/lzFVQaZpXHxxXDyyelBsUGNG6jS28PfNy+gu2FD9jGSPbc0CqKqdjOudjpwzz3ufPHz69eTZulSeNvbeo+XDNAnt/enUE1LiR4SAE2m0di7txvAhfQMsj7xc47HWcRy+s+oaN/3tPvy6U9313c6bu7u7du7y+LUIiFpMuJ3QaQ3CJ4mexZx76o4SB+Pr4DB9EbKGms0roHlEEKsR9s+rZnBza+tFWne1lFDqbuWU6UVkuU6KhrQTbo04s+yZe44RcYtJFtx/v+JiYWpIvJaOXnpSIokp4trwlmz1PW7t/H4l2RguUhZSMp+2GG99zqOexRNRtgvRUtdrdxB1u6baomXiQmVBHMlpVDmxuYVljIFqU6fcJ0FuIpvOG/fMsE9X4nH+X5Cj5XmL08G1dMGGPZ7zqGD4nyDmFbWivZ4S1O48WfJkuLPPil7Mp7lGxvfrdbvupL3PelK6tdba9CB7aKuv34uvjIKOm2/IrG3Gt55MwxJyt7YfoWlaCGpo+dLHcdIUuT+FFWUdRnPoj73ZCsulnvduoU9lNLiC/1kzHrJs3rO+OuzkuGlnWfp0nTDkNW7Kku+rHX9jFXcgvCDv1m939at68YH/JZGlWcaU6bTRb9z+PuHyFJWj2S1qkPjkzW982YYkpS9sXXXztvaYoiPW/YlK6KIqsjS71gh9ybZwWD9+vDWQh793CnJrpiTk2HHXbeuNxCb1SmhyD0oQjKbbd47VLRlUKSjQvK6Qt2URdyE/Z53skNB3JrNkjlLhng+69Dut9ZiaKjFULdi6ydL22IMoYxCv++Ql9t/IbMG9RU1giEtBr+3T16NP+24vkIMVV5la5f+NaYZhrqMUpZyDjHURVw/accr26PLj79MTrr7EzKWyB/c6OecCo1P1vAuDdUwACcCtwN3AO9KWb8M+Fi0/svACm/dhmj57cDLQs5Xe4whtDAPS8k1pfRDKPLiD6plU5WkXOvXd4PC/RSHr/w7HdfVM01x5cUY0s4VInNomahanpNKbGam+3/Jkv7utqLlN7l9VYUfKlPZ8um34PzutHkuodigxLGXPEMywPd/aIYB6AB3As8AlgK3AEcktnkr8DfR71OBj0W/j4i2XwYcGh2n0++ctfdKKhLUHLSSa6sy9QktuG1uXfgK3Pf7r1+f73NOuoPiXlN1uMnytqmqbNP+Z5WzZAshNn7DfDZVFH7Rc+UZ8jz5fOMQ9xTLkjk0PpI0IDXf72Eahmngau//BmBDYpurgeno9yTwfUCS2/rb5X0aSbs9iIBvGsM6T5JBvPij0LpIupWSL3eaQk2OzM0bMFe290oZX3ro8fq5YZKGodMZfgVlmIaoSPpxX75k9+k84xLqDkzeez+GUQOhhqGOAW4HATu9/7uAF2Rto6p7ROSHwFS0/EuJfQ9KO4mIrAXWAixfvrwGsT1CBowNaxKNJibrGFQ64SID8ZrIXDk3B9df35t5Nn4lsyb+mZ6GSy6Bt761O1gt7TmF3NOswVRbt3bnoX74YTfYssq9Sbu3eYMQ16yBD3+4e32q9TyPKpNh1X18f5/Q9ONp8oVmxC2SgbcNhFiPvA/wauBS7/8bgA8mtrkVONj7fyfOMHwI+B1v+d8Cp/Q7Z2NzPo9rjKENLp+6WgyhsqW5A2JfcWgLJ697a0g357TrnZnp7dnUzx8dci/yck1lrQ+JhQwi7lGWvOPnyZl0DZYZG1KUfj2b/PhOQ66kOloMu4BDvP8HA/dlbLNLRCaBxwMPBu7bHoY1icawJ+so0kppQ+uiDtniWrSrkLgUDvvtFz7FZZH5N9LuadrkSQBnndWtqYu4eajjFBF+aox+pM2zkZzmM76GCy9c2KLYsKGbliLtXvgpMDod14pau7b/vNXJ1CIhzzukJZDV4uxXJvz5q+PrGPS7t2qVm5Nl71533i1butPKTk872ZtOyxFiPfI+uJjBXbjgcRx8PjKxzZn0Bp8/Hv0+kt7g8100EXxuI2Vq5VVq8oMKKA8Tv7dIP9n8GmbcfXAQtdise5qWrDFvnEPRGndWv/mQQYghgdg0Wfv1tCkTNwm97qzt+o2r8Dse1D3/Sh5FymqNMKwWg7qYwVm4wHEH2KKqt4rIeZEQV+JcRH8vInfgWgqnRvveKiIfB74O7AHOVNW9qSdaTJSplVetyYe2UqrGQAaVpGxuDi67rNsCyJq5LGZQiQiT58g6bto0pvEMdn4tPG0WtZA4Q/I5Qfb+af7vZEsg7fjxFJ/gvvNiIcn7HXo9odtlPc+88hofO076lzeVbN2sWQOXXz7cWGIRQqxH2z5j32IoUysfZk2+rt42ddbQ8/y2VWXOq/WXnekvq1ti1vnK3Dv/WKH7J1sCeT2SknGIfrGMvBZE3jmK9hhKuw9p8aCq9zTr2orINeTxStjI5xGmbIFt+/iHukfkJpeHujGKBG+z9isyb0eWYiwaZMw6TpEAcMj4iWRaj7wuk/0UZdmgsL9vv2lK++HLkMxqW0TRl3WHpR2noQGsZhhGnTpjDA0WxAVy1N0fP7k+7TrLGqS8/UJn+ivj+w4lqzZd9XknpyCt0pe+ynXW1Qr2jxO3KLPKT175SsqTlkol1Ng1VIEzw2A42taSqOLSKTtTXpMthizlVoeRnJzsnt/PElv1eWe5ucpQRZ66ym58HL9rclr5KdrFONliSPufLOsNd94ww2A4RiH2EHLc2KUQkomyTtmqxhiquFLyyOp/n4y1xBPvFKXOZ9kGP3wcZ8ibHjTtWYW4yvyWQnzv455OafEWazGYYWicYRXEQZ4n+cLVORH8KA1Y9I+T5X+Pa/t+a2YAA6UKydomQlw9RQP2/r5xS6vTKTdP+YAxw2B0GUZBHGTLZBBGp20utn6kuXeynmsy307duZxCZB2le5tFmYwAcYeCycmwaWiHTKhhqGPks9F2hjGSepA5ngYx5iBvpGzTo07T2LrVjS0A9711K/z1X6fLuGaNG9MRb79kSblcTmWpI+9VG55D0TI9O9sdzazqRq0vX96+shSAGQajHgY9YKyscctSMGkv/SCV5SBJXuP0NFx7bTeFRpxuwWeQSQvrGATZhucQUqb9e5+87rT7PiqENCva9jFXkhFE0e6tbU734bspsuZSLpLaY9DunqqB9VHoMBESrG4ZmCvJWPT0qxUnWyFNpDyPmZvLr+FnJVfzr3HvXpiZcakW+tWy29jCi2vfU1PDeQ5VWyZp5SuZpn1EMcMwKNrgIx0Gbb7Ooop+GPmT0pibg5e8pBsT2LIl3bWTpmzja4znb1ANdw0NO4tvHmnZYEOy3FahqjutyYrEgDHDMAja4iMNoYpib/t1llH0TSjLWEHFPPJIuJKKr3HrVhdw3rNnNJVUUknv3u1q34OkqmJvqiIxBMwwDIJBBvbqZBBN6bZdZ5tqxVnECiqvF1EeRWcTayNN1L7rUOyjUL5KYIZhEBQt5E25Y6wp3Q5CehGFHmdUlVRTte9RvmcDRFygerRYuXKl7tixo2kx8glV9k26Y+o4d5tjDEZ91Pmcmywzi7y8isgNqrqy33bWYhgUoTWRJt0x1pQePOOgiOqsQMQTAY1qRWiRYIahaZp2x5hiHxzjooiqVl78+zAx4Y6zb9/CYw3aiI5CTKwlmGFomjHu2bDoGRdFVLXy4t8HVWccRHqPNQwjOqhK2Di0ChOYYWgDVmsfT5puDdZFv8pLP8WYvA9pYxSGYUQHUQkbl1ZhAjMMhjEoxqk1mFV5CVGMIfdhWEa07krYuLQKE5hhGBZj2Nxc9IQ803FvDYYqxn73YVBGdNDv3bi0ChOYYRgGY9rcXNTYM3XUqRjrNqLDeEbj1Cr0MMNQF3k1kzFtbi5a5ubg3HPdSOW03jWLiTYrxmG9d2PYKjTDUAf9aiZj2txclMTPOjYKExPj/UxH2V1m711pzDDUQb+aSZtrVUYx4mcdG4Xjj3eth3F8pqPuLrP3rjRmGOogpGbS1lqVUYzksx5XowDj4QK1964UlQyDiBwAfAxYAdwNvEZVf5Cy3WnAu6O/71fVy6Pls8BTgZ9F61ar6veqyNQIVjNZPCymZ22umEVLpSR6IrIReFBVLxKRdwH7q+o7E9scAOwAVgIK3AAco6o/iAzDO1S1UEa8kUiiZxjjgHWzHiuGlUTvJGBV9PtyYBZ4Z2KblwHXqOqDkWDXACcCH6l4bsMwBo25YhYlExX3f7Kq3g8QfT8pZZuDgJ3e/13RspjLRORmEXmPiEjWiURkrYjsEJEdDzzwQEWxDcMwjCz6thhE5HPAU1JWnRN4jjRlH/uvXq+q94rI44ArgDcAW9MOoqqbgc3gXEmB5zYMwzAK0tcwqOrxWetE5Lsi8lRVvV9EngqkBY530XU3ARyMczmhqvdG3z8WkX8CjiXDMBiGYRjDoaor6UrgtOj3acBnUra5GlgtIvuLyP7AauBqEZkUkQMBRGQJ8ArgaxXlMQzDMCpS1TBcBJwgIt8CToj+IyIrReRSgCjofD6wPfqcFy1bhjMQXwFuBu4FPlxRHsMwDKMiNuezYRjGIiG0u2rVFoNhGIYxZphhMAzDMHoww9AEc3Nw4YXu2zAMo2VYEr1hM+oZKw3DGHusxTBs0jJWGoZhtAgzDMMmzljZ6VjGSsMwWom5kobNYkrbbBjGSGKGoQksY6VhGC3GXEmGYRhGD2YYDMMwjB7MMBiGYRg9mGEwDMMwejDDYBiGYfRghsEwDMPoYSTTbovIA8C3h3CqA4HvD+E8RTG5wmmjTGByFaGNMsFoyvV0VX1ivwOMpGEYFiKyIyR3+bAxucJpo0xgchWhjTLBeMtlriTDMAyjBzMMhmEYRg9mGPLZ3LQAGZhc4bRRJjC5itBGmWCM5bIYg2EYhtGDtRgMwzCMHswwGIZhGD2YYYgQkS0i8j0R+Zq37AARuUZEvhV9798SuT4gIt8Qka+IyKdE5AlNy+Ste4eIqIgcOEyZ8uQSkbeJyO0icquIbGyDXCLyPBH5kojcLCI7ROTYIct0iIhcKyK3RfflD6LljZb5HLmaLvOpcnnrh17u82SqXOZV1T4uzvIi4Gjga96yjcC7ot/vAi5uiVyrgcno98XDlitNpmj5IcDVuMGHB7bkXr0E+BywLPr/pJbI9Vng16PfLwdmhyzTU4Gjo9+PA74JHNF0mc+Rq+kynypX9L+Rcp9zryqXeWsxRKjqF4EHE4tPAi6Pfl8OnDxUoUiXS1U/q6p7or9fAg5uWqaIvwDWA430aMiQ6/eAi1R1Ptrmey2RS4Ffin4/HrhvyDLdr6o3Rr9/DNwGHETDZT5LrhaU+az7BQ2V+xyZKpd5Mwz5PFlV7wf3EIAnNSxPGqcD/960ECLySuBeVb2laVkSPAv4NRH5soh8QUSe37RAEWcDHxCRncCfARuaEkREVgBHAV+mRWU+IZdPo2Xel6st5T5xryqXeZvac4QRkXOAPcA/NizHo4FzcM39tjEJ7A/8KvB84OMi8gyN2tgN8nvA21X1ChF5DfC3wPHDFkJEHgtcAZytqj8SkWGLkEpSLm95o2XelyuSo/Fyn/IMK5d5azHk810ReSpA9D10N0QWInIa8Arg9S1Qcs8EDgVuEZG7cc38G0XkKY1K5dgF/LM6rgf24ZKMNc1pwD9Hvz8BDDX4DCAiS3AK5R9VNZal8TKfIVfjZT5FrsbLfca9qlzmzTDkcyXuBSb6/kyDsvwCETkReCfwSlX9adPyqOpXVfVJqrpCVVfgCubRqvqdhkUD+DTwUgAReRawlHZkxLwPeHH0+6XAt4Z5cnFNg78FblPVP/dWNVrms+RqusynydV0uc95htXL/LAi6G3/AB8B7gcewT3gNwNTwOdxL+3ngQNaItcdwE7g5ujzN03LlFh/N830Skq7V0uBfwC+BtwIvLQlcr0QuAG4BecXPmbIMr0QFyz9ileOXt50mc+Rq+kynypXYpuhlvuce1W5zFtKDMMwDKMHcyUZhmEYPZhhMAzDMHoww2AYhmH0YIbBMAzD6MEMg2EYhtGDGQbDMAyjBzMMhmEYRg//H00Yi4y28xMAAAAAAElFTkSuQmCC\n",
"text/plain": [
"