{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"Open\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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Algorithme EM" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**2) On va maintenant appliquer l'algorithme EM aux données précédentes.** \n", " - **Ecrire une fonction `EM(K,nbpas,data)` prenant en entrée les données 'data', un entier K (nombre de gaussiennes du mélange) et un nombre de pas nbpas.**\n", " - **Appliquer la fonction précédente aux données pour les paramètres du GMM sous-jacent pour un nombre de classes K donné.**\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On fixe le nombre de pas de EM et un intervalle pour les valeurs de K que l'on va tester. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "nbpas=int(5e1) # nombre de pas dans EM\n", "rangeK=range(1,5) # valeurs de K testees pour le nombre de classes (dans la selection de modele)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sélection de modèle\n", "\n", "Pour la sélection de modèle, il faut calculer la log-vraisemblance et la complexité du modèle.\n", "\n", "Une manière classique de choisir le nombre de composantes $K$ dans un intervalle d'entiers est de minimiser en $K$ \n", " $$-\\ell(x;\\widehat{\\theta}_K)+ \\frac{\\operatorname{dim}_K \\times \\log n}{2},$$\n", " avec \n", " + $\\widehat{\\theta}_K$ l'estimateur de $\\theta$ issu de l'algorithme EM dans le modèle de mélange avec $K$ classes,\n", " + $\\ell(x;\\hat{\\theta}_K)$ la log-vraisemblance de l'échantillon observé $(x_1,\\dots , x_n)$ \n", " + $\\operatorname{dim}_K$ le nombre de degrés de liberté du modèle de mélange utilisé avec $K$ classes (voir le cours).\n", "\n", "Le K optimal pour la sélection de modèle est l'argmin de cette fonction.\n", "\n", "**3) Ecrire une fonction `logLikelyhood(K,alpha,mu,Cov,data)` permettant de calculer la vraisemblance du GMM donné par les paramètres (K,alpha,mu,Cov) sur les données `data`.** \n", "\n", "**4) Implémenter la sélection de modèle et l'appliquer aux données précédentes pour déterminer le K optimal. Commentez le résultat.**\n", "\n", "**5) Afficher les données en les colorant en fonction de leur classe telle qu'estimée par EM.**\n", "\n", "**6) Afficher la fonction $-\\ell(x;\\widehat{\\theta}_K)+ \\frac{\\operatorname{dim}_K \\times \\log n}{2},$ en fonction de $K$**\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8VNXdx/HPjx0VQSAoAoIEtNaN0rRarXvBureVKtWKWivuj1qt+4JEqEu1uFUfrVrXqtVqqeJCtdZq3YKK+qhFQJBNWVREcSHwe/44N2USEjJJZnJm7nzfr9e8MnPnzMzvcvWbkzPn3mPujoiIpFeb2AWIiEh+KehFRFJOQS8iknIKehGRlFPQi4iknIJeRCTlFPRS1MxsgJm5mbXLou0RZvZsjj8/5++ZvO+hZvZErt9XSpOCXlqNmc0ys6/NrGed7a8lYT0gTmVx1ffLyt3vcvfhMeuS9FDQS2t7D/hZzQMz2xroHK8ckfRT0EtruwMYlfH4cOD2zAZm1tXMbjezRWY228zOM7M2yXNtzey3ZrbYzGYC+9Tz2pvNbIGZzTOzi82sbWNFmVknM7vTzJaY2Sdm9rKZbdjU9zSzb5jZZDP7yMz+Y2YHZTzX2cyuSPZpqZk9a2adgWeSJp+Y2Wdm9r26Q0JmtkNS09Lk5w4Zzz1tZpVm9pyZLTOzJ+r+1SSlTUEvre0FYH0z2yIJy4OBO+u0uQboCgwEdiH8Yjgyee5oYF/gW0AFMKLOa28DqoFBSZvhwC+zqOvw5DP7AT2AY4EvmvKeZrYuMBm4G+hF+Mvl92a2ZdLkt8C3gR2A7sAZwCpg5+T5bu6+nrs/X+d9uwOPAFcntV0JPGJmPTKaHUL4N+oFdABOz2KfpUQo6CWGml79MOAdYF7NExnhf7a7L3P3WcAVwGFJk4OACe4+x90/An6T8doNgb2AU9z9c3dfCPwOGJlFTSsIITrI3Ve6+xR3/7SJ77kvMMvdb3X3and/BXgAGJH8RfIL4GR3n5d8xr/d/assatsHeNfd70je90/Jv9t+GW1udfdp7v4FcB8wJIv3lRLR6EwFkTy4gzBcsSl1hm2AnoQe6eyMbbOBPsn9jYE5dZ6r0R9oDywws5ptbeq0X1tN/YB7zKwb4a+Mc5v4nv2B7czsk4xt7ZL37gl0AmZkUUtdG1N7P6H2vwnABxn3lwPrNeNzJKUU9NLq3H22mb0H7A0cVefpxYTedX/grWTbJqzu9S8gBDIZz9WYA3wF9HT36ibWtAK4CLgomf0zCfhP8jPb95wD/NPdh9V9IunRfwmUA1Prfnwj7zuf8O+RaRPgsUZeJwJo6EbiOQrY3d0/z9zo7isJQw/jzKyLmfUHfsXqcfz7gP8xs75mtgFwVsZrFwBPAFeY2fpm1sbMys1sl8aKMbPdzGzrZOjoU8Ivm5VNfM+Hgc3M7DAza5/cvmNmW7j7KuAW4Eoz2zj5Uvl7ZtYRWEQYqx/YQHmTkvc9xMzamdnBwDeTzxNplIJeonD3Ge5e1cDTJwGfAzOBZwlfbt6SPHcT8DihV/wK8Jc6rx1FGPp5C/gYuB/onUVJGyVtPwXeBv7J6l8uWb2nuy8jfFE7ktAL/wC4FOiYNDkdeAN4Gfgoea6Nuy8HxgHPJTN+tq/zvksI4/+nAUsIX+Lu6+6Ls9gvEUwLj4iIpJt69CIiKaegFxFJOQW9iEjKKehFRFKuIObR9+zZ0wcMGBC7DBGRojJlypTF7l7WWLuCCPoBAwZQVdXQTDsREamPmdU9Y7peGroREUk5Bb2ISMop6EVEUk5BLyKScgp6EZGUU9CLiKScgl5EJOWKO+hnzIBTToEVK2JXIiJSsIo76N95B666Cm67LXYlIiIFq7iDfu+9YbvtoLISvspmjWURkdKTVdCb2Swze8PMXjOzqjrPnW5mbmY9k8dmZleb2XQze93Mhuaj8OTDQ8i//z784Q95+xgRkWLWlB79bu4+xN0rajaYWT9gGPB+Rru9gMHJbTRwfS4KbdAPfgA77wzjxsEXX+T1o0REilFLh25+R1i/MnM9wgOA2z14AehmZtms2dk8Nb36BQvg+vz+ThERKUbZBr0DT5jZFDMbDWBm+wPz3H1qnbZ9gDkZj+cm2/Jn551h2DC45BL47LO8fpSISLHJNuh3dPehhGGZE8xsZ+Bc4IJ62lo929ZYgdzMRptZlZlVLVq0KOuCG1RZCYsWwbXXtvy9RERSJKugd/f5yc+FwIPALsCmwFQzmwX0BV4xs40IPfh+GS/vC8yv5z1vdPcKd68oK2v0uvmN22472GcfuOwyWLq05e8nIpISjQa9ma1rZl1q7gPDgZfdvZe7D3D3AYRwH+ruHwATgVHJ7JvtgaXuviB/u5Bh7Fj4+GOYMKFVPk5EpBhk06PfEHjWzKYCLwGPuPtja2k/CZgJTAduAo5vcZXZGjoUfvITuPJK+OijVvtYEZFCZu5rDJ+3uoqKCs/ZUoJvvgnbbANnnQXjx+fmPUVECpCZTcmc8t6Q4j4ztj5bbQUjR8LVV8PChbGrERGJLn1BD3DhheHkqUsvjV2JiEh06Qz6zTeHUaPg97+H+WtM+BERKSnpDHqACy6A6mqN04tIyUtv0G+6KRx1FNx0U7jomYhIiUpv0AOce274efHFcesQEYko3UHfrx8ccwzccktYjUpEpASlO+gBzj4bOnQIZ82KiJSg9Ad9795wwglw551h6UERkRKT/qAHOOMM6NwZxoyJXYmISKsrjaAvK4NTToF774XXX49djYhIqyqNoAc47TTo2jWcNSsiUkJKJ+g32CCE/UMPwZQpsasREWk1pRP0ACefDN27h7NmRURKRGkF/frrhy9mJ02C55+PXY2ISKsoraAHOPFE6NULzj8/diUiIq2i9IJ+3XXDSVRPPglPPx27GhGRvCu9oAc49ljo0yf06gtghS0RkXwqzaDv1Clc8OzZZ2Hy5NjViIjkVWkGPYRLGPfvD+edp169iKRa6QZ9hw5hmuXLL8PDD8euRkQkb0o36CEsNzhoUBirX7UqdjUiInlR2kHfrl240NnUqfCXv8SuRkQkL0o76AFGjoRvfjNcA2flytjViIjknIK+bdvQq3/rLbjnntjViIjknIIe4MADYdttQ+BXV8euRkQkp7IKejObZWZvmNlrZlaVbLvczN4xs9fN7EEz65bR/mwzm25m/zGzPfNVfM60aROWGpw+HW6/PXY1IiI51ZQe/W7uPsTdK5LHk4Gt3H0bYBpwNoCZfRMYCWwJ/BD4vZm1zWHN+bHffvCd74TA//rr2NWIiORMs4du3P0Jd68Z53gB6JvcPwC4x92/cvf3gOnAd1tWZiswg8pKmD0bbrkldjUiIjmTbdA78ISZTTGz0fU8/wvg0eR+H2BOxnNzk221mNloM6sys6pFixY1peb8GT4cdtwRLr4YvvwydjUiIjmRbdDv6O5Dgb2AE8xs55onzOxcoBq4q2ZTPa9f4xoD7n6ju1e4e0VZWVkTy84TsxDy8+bB//5v7GpERHIiq6B39/nJz4XAgyRDMWZ2OLAvcKj7fy8YMxfol/HyvsD8XBWcd7vuCrvvDuPHw+efx65GRKTFGg16M1vXzLrU3AeGA2+a2Q+BM4H93X15xksmAiPNrKOZbQoMBl7Kfel5VFkJCxfCddfFrkREpMWy6dFvCDxrZlMJgf2Iuz8GXAt0ASYn0y5vAHD3/wPuA94CHgNOcPfiOuV0hx1gr73gssvg009jVyMi0iLmBXCJ3oqKCq+qqopdRm1VVaunW2rZQREpQGY2JWPKe4N0ZmxDKirgRz+CK66Ajz+OXY2ISLMp6Nfmootg6VK48srYlYiINJuCfm222QYOOggmTIDFi2NXIyLSLAr6xowZA8uXhy9mRUSKkIK+MVtsAYceCtdeCx98ELsaEZEmU9Bn48ILw4XOfvOb2JWIiDSZgj4b5eVw5JFwww0wZ07j7UVECoiCPlvnnQfuMG5c7EpERJpEQZ+t/v1h9Gi4+WZ4773Y1YiIZE1B3xTnnAPt2oWzZUVEioSCvik23hiOPz4sNzhtWuxqRESyoqBvqjPPhM6dw1mzIiJFQEHfVL16wUknwZ/+BG++GbsaEZFGKeib49e/hi5dwlmzIiIFTkHfHN27w6mnwgMPwKuvxq5GRGStFPTNdeqpsMEGcMEFsSsREVkrBX1zde0ahnAefhhefDF2NSIiDVLQt8RJJ0FZmVagEpGCpqBvifXWg7POgsmT4ZlnYlcjIlIvBX1LHXcc9O4devUFsP6uiEhdCvqW6tw5XBrhmWfgySdjVyMisgYFfS4cfTT066devYgUJAV9LnTsGEL+hRdg0qTY1YiI1KKgz5UjjoCBA8O8evXqRaSAKOhzpX37sOTgK6/AQw/FrkZE5L8U9Ll06KGw+eahV79qVexqRESALIPezGaZ2Rtm9pqZVSXbupvZZDN7N/m5QbLdzOxqM5tuZq+b2dB87kBBads2XL74zTfhvvtiVyMiAjStR7+buw9x94rk8VnAk+4+GHgyeQywFzA4uY0Grs9VsUXhpz+FrbcOwzjV1bGrERFp0dDNAcBtyf3bgB9lbL/dgxeAbmbWuwWfU1zatAlLDU6bBnfdFbsaEZGsg96BJ8xsipmNTrZt6O4LAJKfvZLtfYA5Ga+dm2yrxcxGm1mVmVUtWrSoedUXqgMOgG9/OwzjrFgRuxoRKXHZBv2O7j6UMCxzgpntvJa2Vs+2NeYbuvuN7l7h7hVlZWVZllEkzEKv/r334NZbY1cjIiUuq6B39/nJz4XAg8B3gQ9rhmSSnwuT5nOBfhkv7wvMz1XBRWOvveB734PKSvjyy9jViEgJazTozWxdM+tScx8YDrwJTAQOT5odDvw1uT8RGJXMvtkeWFozxFNSzELIz50LN90UuxoRKWHZ9Og3BJ41s6nAS8Aj7v4YcAkwzMzeBYYljwEmATOB6cBNwPE5r7pY7L477LorjB8Py5fHrkZESlS7xhq4+0xg23q2LwH2qGe7AyfkpLpiV9Or32knuP56OO202BWJSAnSmbH59v3vw/DhcMklsGxZ7GpEpAQp6FtDZSUsXgzXXBO7EhEpQQr61vDd78J++8Hll8Mnn8SuRkRKjIK+tYwdG0L+d7+LXYmIlBgFfWsZMgRGjAhBv2RJ7GpEpIQo6FvTmDHw2WdhCEdEpJUo6FvTllvCIYeEL2U//DB2NSJSIhT0re3CC+Grr8J0SxGRVqCgb22DB8OoUeEEqnnzYlcjIiVAQR9DzVKD48fHrkRESoCCPoYBA+Coo8LFzmbPjl2NiKScgj6Wc88Nq1FVVsauRERSTkEfS9++cOyx8Mc/wvTpsasRkRRT0Md01lnQoUNYclBEJE8U9DFttBGcdFJYRPytt2JXIyIppaCP7de/hnXXDWfNiojkgYI+tp494dRT4c9/hqlTY1cjIimkoC8Ev/oVdOsWzpoVEckxBX0h6NYtLDP417/Cyy/HrkZEUkZBXyhOPhl69AhnzYqI5JCCvlB06QJnngmPPQbPPRe7GhFJEQV9ITnhBNhwQzj//NiViEiKKOgLyTrrwDnnwD/+AU89FbsaEUkJBX2hGT06XB7h/PPBPXY1IpICCvpC06kTnHce/Pvf8PjjsasRkRRQ0BeiI48MlzJWr15EciDroDeztmb2qpk9nDzew8xeMbPXzOxZMxuUbO9oZvea2XQze9HMBuSn9BTr0CFMs6yqgokTY1cjIkWuKT36k4G3Mx5fDxzq7kOAu4Hzku1HAR+7+yDgd8CluSi05Bx2WFh2sGY1KhGRZsoq6M2sL7AP8IeMzQ6sn9zvCsxP7h8A3Jbcvx/Yw8ys5aWWmHbtwoXOXn8d7r8/djUiUsSy7dFPAM4AMruWvwQmmdlc4DDgkmR7H2AOgLtXA0uBHnXf0MxGm1mVmVUtWrSomeWn3MEHw5ZbhmvgrFwZuxoRKVKNBr2Z7QssdPcpdZ46Fdjb3fsCtwJX1ryknrdZ4xtFd7/R3SvcvaKsrKyJZZeItm3DoiTvvAN33x27GhEpUtn06HcE9jezWcA9wO5m9giwrbu/mLS5F9ghuT8X6AdgZu0Iwzof5bLokvLjH8OQISHwV6yIXY2IFKFGg97dz3b3vu4+ABgJPEUYh+9qZpslzYax+ovaicDhyf0RwFPumiPYbDULiM+YAbffHrsaESlC7ZrzInevNrOjgQfMbBXwMfCL5OmbgTvMbDqhJz8yJ5WWsn32ge22g7Fj4ec/h44dY1ckIkXECqGzXVFR4VVVVbHLKGyTJ8Pw4XDddXD88bGrEZECYGZT3L2isXY6M7ZY/OAHsNNOMG4cfPFF7GpEpIgo6IuFWRirnz8fbrghdjUiUkQU9MVkl11Cz/6SS+Czz2JXIyJFQkFfbCorYeFCuPba2JWISJFQ0Beb7bcPs3Auvxw+/TR2NSJSBBT0xWjsWPjoI5gwIXYlIlIEFPTFaOjQcMbsFVeEwBcRWQsFfbG66CJYtiyEvYjIWijoi9XWW4erW151FejqnyKyFgr6YjZmTDh56lKt7SIiDVPQF7PNNw8rUV13HSxYELsaESlQCvpid8EFUF0N48fHrkRECpSCvtgNHAi/+AXceCO8/37sakSkACno0+C8ZF32iy+OW4eIFCQFfRr06wejR8Ott8LMmbGrEZECo6BPi3POgXbtwlmzIiIZFPRp0bs3nHAC3HFHWExcRCShoE+TM8+Ezp3DWbMiIgkFfZqUlcHJJ8M998Abb8SuRkQKhII+bU47DdZfHy68MHYlIlIgFPRp0717CPsHH4QpU2JXIyIFQEGfRqecEgL/ggtiVyIiBUBBn0brrw9nnAGTJsHzz8euRkQiU9Cn1YknQq9e6tWLiII+tdZdF846C/7+d/jnP2NXIyIRKejT7NhjYeON4fzzwT12NSISSdZBb2ZtzexVM3s4eWxmNs7MppnZ22b2Pxnbrzaz6Wb2upkNzVfx0ojOneHcc+Ff/4LJk2NXIyKRNKVHfzLwdsbjI4B+wDfcfQvgnmT7XsDg5DYauL7lZUqzHXUU9O+vXr1ICcsq6M2sL7AP8IeMzccBY919FYC7L0y2HwDc7sELQDcz653DmqUpOnYMIf/SS/DII7GrEZEIsu3RTwDOAFZlbCsHDjazKjN71MwGJ9v7AHMy2s1NttViZqOT11Yt0uLW+TVqFJSXh8Bftarx9iKSKo0GvZntCyx097qnWXYEvnT3CuAm4Jaal9TzNmuMGbj7je5e4e4VZWVlTSxbmqR9+7CQ+GuvhTNmRaSkZNOj3xHY38xmEcbhdzezOwk99QeSNg8C2yT35xLG7mv0BebnpFppvp/9DLbYIsyrX7kydjUi0ooaDXp3P9vd+7r7AGAk8JS7/xx4CNg9abYLMC25PxEYlcy+2R5Y6u4Lcl+6NEnbtqFX/9ZbcO+9sasRkVbUknn0lwAHmtkbwG+AXybbJwEzgemEIZ3jW1Sh5M6IEbDNNiHwq6tjVyMiraRdUxq7+9PA08n9Twgzceq2ceCEHNQmudamTVhq8Ec/CitRHXlk7IpEpBXozNhSs//+UFERAv/rr2NXIyKtQEFfasygshJmzYJbbmm0uYgUPwV9KdpzT9hhB7j4Yvjyy9jViEieKehLkVkI+Xnz4MYbY1cjInmmoC9Vu+0WbuPHw/LlsasRkTxS0Jeyykr48EO47rrYlYhIHinoS9mOO8IPfwiXXgrLlsWuRkTyREFf6saOhSVL4KqrYlciInmioC913/kOHHAA/Pa38PHHsasRkTxQ0Evo1S9dCldeGbsSEckDBb2E698cdBBMmACLF8euRkRyTEEvwZgxYZrl5ZfHrkREckxBL8EWW8Ahh8A118AHH8SuRkRySEEvq114YbjQ2SWXxK5ERHJIQS+rDRoERxwB118Pc+fGrkZEckRBL7Wdfz64w7hxsSsRkRxR0Ett/fvD0UfDzTfDe+/FrkZEckBBL2s655ywGlVlZexKRNLl66/hnXfgb3+DK66AY46BP/857x/bpKUEpUT06QPHHw9XXw1nnQWbbRa7IpHisWoVzJkD774L06bVvr33Xni+Rs+eMHBg3kuysMRrXBUVFV5VVRW7DMm0cCFsumlYX/auu2JXI1JY3GHRohDedQN9+vTaC/qsu27oLA0eHH7W3AYPhu7dW1SGmU1x94rG2qlHL/Xr1QtOOgkuuywM5Wy5ZeyKRFrfp5+uDvK6gb506ep27dtDeXkI8D33rB3ovXuHxX4iUo9eGrZkSejV77lnq4wjikTx1VcwY0b9YZ558qAZbLJJ7R55zf3+/aFd6/eb1aOXluvRA049NVz07LXXYMiQ2BWJNM/KlfD++6sDPDPQZ8+uPW7eq1cI7733rh3o5eXQuXO8fWgB9ehl7T75JHxZ9P3vw8SJsasRaZh7WDEts0deE+jTp4cZLzW6dKk9vFIT6IMHQ7du8fahidSjl9zo1g1OPx3OPRdefBG22y52RVLqPvlkzSGWmseZK6V16BDO9t5sM9h339q98w03jD5u3prUo5fGffZZGKsfOhQefzx2NVIKvvhi9bh53UBfuHB1OzMYMKD+3vkmm0DbttF2oTXkvEdvZm2BKmCeu++bsf0a4Eh3Xy953BG4Hfg2sAQ42N1nNa18KSjrrRfm059+OvzrX7DTTrErkjSorg7j43XDfNq0MA89sxO60UYhwPffv3agDxwIHTvG24ci0ZShm5OBt4H1azaYWQVQd0DrKOBjdx9kZiOBS4GDW1qoRHbccWG5wfPPh3/8o6T+7JUWcIcFC+oP85kzYcWK1W27dg3hvdNOtcN80CBYf/2GP0MalVXQm1lfYB9gHPCrZFtb4HLgEODHGc0PAMYk9+8HrjUz80IYI5LmW2edME5/0knw1FOwxx6xK5JC8tFH9Z8J+u678Pnnq9t16hSCe8st4cc/rh3oPXuqA5En2fboJwBnAF0ytp0ITHT3BVb74PQB5gC4e7WZLQV6ALXWqDOz0cBogE022aRZxUsrO/rocALVeefB7rvrf8pSs3x5mL1SX+98yZLV7dq2Dd/pDB4Mu+xSO8z79g3XUZJW1WjQm9m+wEJ3n2JmuybbNgZ+Cuxa30vq2bZGb97dbwRuhPBlbPYlSzQdO4aQP+YYePTRMM9Y0mXFinA9lvp653XXKOjTJ4T3gQfWDvNNNw0zXqRgNDrrxsx+AxwGVAOdCGP0XyW3mgs6bALMTMblHwfGuPvzZtYO+AAoW9vQjWbdFJEVK+Ab3wjTLquq1KsvRqtWwbx59Z8JOnNmOLmoxgYbwOabr3km6KBB4Ut6iSpns27c/Wzg7ORNdwVOz5x1k2z/zN0HJQ8nAocDzwMjgKc0Pp8i7dvDBReElageeiiMs0rhWbky9MBnzKh9qwn3L75Y3bZz5xDeQ4bAQQfVDvQePeLtg+RMk+bRNxL0NdMrOwF3AN8CPgJGuvvMtb2vevRFproattoqhP7UqRpzjWX58tADnzlzzUCfNav2jJb27cN88/p65xtvrGNYpLLt0euEKWmee+6Bn/0s/DxYs2fzwh0WL64d4JmhvmBB7fZdu4brsZSXh/nlNffLy8OXoCk/eagUKeglv1atgm23Db3GN9+McuW+VKiuDicH1e2R1wR65in9EL4AzQzwzEDv3l3fmZQYXetG8qtNm3BVy5/8BO6+G0aNil1R4fr88zV74zW32bND2Nfo0CHMWikvDycOZQb6ppsW7dUTJS716KX53KGiIlxk6p13wjhwKXIP119paIjlww9rt99gg/qHV8rLQ49d4+WSJfXoJf/MwgLi++wDf/xjOKEqrVasCNczb2iIJfPsT7MwJl5eHq6aWDfUN9gg3n5ISVKPXlrGHXbYIczLfvfd4r7A1LJlDQ+xvP9+7fnlHTvW3yMvLw+rDXXqFG8/pGSoRy+to6ZXP2wY3HQTnHhi7Ioa5h6WhmtoiGXRotrte/QIwb3ddnDIIbXDvHdvDbFI0VCPXlrOHXbbDf7znxCY66wTr5avvw5fcNY3xDJzZph7XqNNG+jXr+FZLF27xtsPkSyoRy+tp6ZXv/POcP31cNpp+f28pUvr75HPmBGmKmau/9m58+rwHjasdqAPGKBrskhJUI9ecmf4cHj11XBRrJZcB2XVqnAyUENDLJlXSgQoK2t4FstGG2luuaSWevTS+iorYfvt4Zpr4Oyz1972q6/Cafp1h1hmzAi/KL78cnXbtm3DsnDl5TBiRO1AHzhQi1KINEI9esmt/faD554LYb1qVcOzWObOrb1U3DrrrNkbrwn0/v1Ld46+yFqoRy9xjB0bFhHv3bv2FRIBNtwwhPcuu6wZ6r16aYhFJE8U9JJb3/pWWFt2+vQ1e+e6frlIFAp6yb18z7oRkSbRGR8iIimnoBcRSTkFvYhIyinoRURSTkEvIpJyCnoRkZRT0IuIpJyCXkQk5QriWjdmtgiY3cyX9wQW57CcmLQvhSkt+5KW/QDtS43+7l7WWKOCCPqWMLOqbC7qUwy0L4UpLfuSlv0A7UtTaehGRCTlFPQiIimXhqC/MXYBOaR9KUxp2Ze07AdoX5qk6MfoRURk7dLQoxcRkbVQ0IuIpFzRBL2Z3WJmC83szQaeNzO72symm9nrZja0tWvMRhb7sauZLTWz15LbBa1dY7bMrJ+Z/cPM3jaz/zOzk+tpU/DHJcv9KIrjYmadzOwlM5ua7MtF9bTpaGb3JsfkRTMb0PqVNi7LfTnCzBZlHJdfxqg1G2bW1sxeNbOH63kuv8fE3YviBuwMDAXebOD5vYFHAQO2B16MXXMz92NX4OHYdWa5L72Bocn9LsA04JvFdlyy3I+iOC7Jv/N6yf32wIvA9nXaHA/ckNwfCdwbu+4W7MsRwLWxa81yf34F3F3ff0f5PiZF06N392eAj9bS5ADgdg9eALqZWe/WqS57WexH0XD3Be7+SnJ/GfA20KdOs4I/LlnuR1FI/p0/Sx62T251Z1wcANyW3L8f2MOs8FZmz3JfioKZ9QX2Af7QQJO8HpOiCfos9AHmZDyeS5H+zwp8L/lz9VEz2zJ2MdlI/tT8FqHXlamojsta9gOK5LgkQwSvAQuBye7e4DFx92r/VJA2AAACFUlEQVRgKdCjdavMThb7AnBgMix4v5n1a+USszUBOANY1cDzeT0maQr6+n77FeNv/1cI16/YFrgGeChyPY0ys/WAB4BT3P3Tuk/X85KCPC6N7EfRHBd3X+nuQ4C+wHfNbKs6TYrmmGSxL38DBrj7NsDfWd0rLhhmti+w0N2nrK1ZPdtydkzSFPRzgczf5n2B+ZFqaTZ3/7Tmz1V3nwS0N7OekctqkJm1J4TjXe7+l3qaFMVxaWw/iu24ALj7J8DTwA/rPPXfY2Jm7YCuFPhwYkP74u5L3P2r5OFNwLdbubRs7Ajsb2azgHuA3c3szjpt8npM0hT0E4FRySyP7YGl7r4gdlFNZWYb1YzNmdl3CcdoSdyq6pfUeTPwtrtf2UCzgj8u2exHsRwXMyszs27J/c7AD4B36jSbCBye3B8BPOXJt4CFJJt9qfN9z/6E71cKiruf7e593X0A4YvWp9z953Wa5fWYtMvVG+Wbmf2JMPOhp5nNBS4kfDmDu98ATCLM8JgOLAeOjFPp2mWxHyOA48ysGvgCGFmI/xMmdgQOA95IxlEBzgE2gaI6LtnsR7Ecl97AbWbWlvDL6D53f9jMxgJV7j6R8EvtDjObTug1joxX7lplsy//Y2b7A9WEfTkiWrVN1JrHRJdAEBFJuTQN3YiISD0U9CIiKaegFxFJOQW9iEjKKehFRFJOQS8iknIKehGRlPt/E6E56QobnTgAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## En utilisant la librarie sklearn\n", "\n", "Refaire la même chose en utilisant la fonction `sklearn.mixture.GaussianMixture` de la bibliothèque `scikit-learn`. Voir [la page suivante](http://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html)\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "import sklearn.mixture" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }