{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\"Open\n", "\n", "# TP Algorithmes Stochastiques, SEM, N-SEM, EM" ] }, { "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", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from TP_SEM_EM import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\renewcommand{\\si}{\\sigma}\n", "\\renewcommand{\\al}{\\alpha}\n", "\\renewcommand{\\tta}{\\theta}\n", "\\renewcommand{\\Tta}{\\Theta}\n", "\\renewcommand{\\Si}{\\Sigma}\n", "\\renewcommand{\\ld}{\\ldots}\n", "\\renewcommand{\\cd}{\\cdots}\n", "\\renewcommand{\\cN}{\\mathcal{N}}\n", "\\renewcommand{\\R}{\\mathbb{R}}\n", "\\renewcommand{\\p}{\\mathbb{P}}\n", "\\renewcommand{\\f}{\\frac}\n", "\\renewcommand{\\ff}{\\frac{1}}\n", "\\renewcommand{\\ds}{\\displaystyle}\n", "\\renewcommand{\\bE}{\\mathbf{E}}\n", "\\renewcommand{\\bF}{\\mathbf{F}}\n", "\\renewcommand{\\ii}{\\mathrm{i}}\n", "\\renewcommand{\\me}{\\mathrm{e}}\n", "\\renewcommand{\\hsi}{\\hat{\\sigma}}\n", "\\renewcommand{\\hmu}{\\hat{\\mu}}\n", "\\renewcommand{\\ste}{\\, ;\\, }\n", "\\renewcommand{\\op}{\\operatorname} \n", "\\renewcommand{\\argmax}{\\op{argmax}}\n", "\\renewcommand{\\lfl}{\\lfloor}\n", "\\newcommand{\\ri}{\\right}\n", "$\n", "\n", "\n", "## Rappels de cours\n", "### Mélange de gaussiennes\n", "On rappelle qu'une mixture gaussienne a une densité de la forme \n", "$$x\\mapsto \\sum_{k=1}^K\\alpha_kf_{\\mu_k,\\sigma_k}(x)=\\sum_{k=1}^K\\alpha_k\\frac{1}{\\sqrt{2\\pi \\sigma_k^2}}\\mathrm{e}^{-\\frac{(x-\\mu_k)^2}{2\\sigma_k^2}}.$$\n", "\n", "On commence par définir une fonction qui permet d'évaluer la densité théorique d'un mélange de gaussiennes sur une grille x. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Calcul de la densité théorique d'un GMM sur une grille x \n", "def densite_theorique(K,mu,sigma,alpha,x):\n", " y=np.zeros(len(x))\n", " for j in range(K):\n", " y+=alpha[j]*sps.norm.pdf(x,loc=mu[j],scale=sigma[j])\n", " return y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation des données\n", "\n", "Pour simuler une mixture gaussienne, on commence par fixer $K\\geq 1$, et par choisir un paramètre $\\theta:=(\\alpha_k,\\mu_k,\\sigma_k)_{1\\le k\\le K}$, avec $\\mu_1, \\ldots, \\mu_K\\in \\R$, $\\sigma_1, \\ldots, \\sigma_K>0$ et $\\alpha_1, \\ldots, \\al_K>0$ tels que $\\al_1+\\cd+\\al_K=1$. On \n", "tire au hasard des données $(z_1,x_1), \\ld,(z_n,x_n\n", ")$, copies indépendantes du vecteur aléatoire $(Z,X)$ de loi $\\p_\\theta$ donnée par :\n", "- $Z\\in \\{1, \\ld, K\\}$ de loi $(\\al_1, \\ld, \\al_K)$,\n", "- pour tout $k\\in \\{1, \\ld, K\\}$, conditionellement à $Z=k$, $X\\sim \\cN(\\mu_k, \\si_k^2)$. \n", "\n", "Notons que le passage de $X$ à $Z$ se fait avec la formule de Bayes : pour tout $x\\in \\R$ et $k\\in \\{1, \\ld, K\\}$, \n", "$$\\label{eq:Oleron1}\\p_\\theta(Z=k | X=x)\\;=\\;\\f{\\al_kf_{\\mu_k,\\si_k}(x)}{\\sum_{l=1}^K\\al_lf_{\\mu_l,\\si_l}(x)}.$$ \n", "avec $$f_{\\mu,\\si}(x):=\\ff{\\sqrt{2\\pi}\\si}\\me^{-\\f{(x-\\mu)^2}{2\\si^2}}.$$\n", "\n", "Les algorithmes SEM, $N$-SEM, EM ont pour but d'estimer les paramètres $\\mu_1, \\ld, \\mu_K\\in \\R$, $\\si_1, \\ld, \\si_K>0$ et $\\al_1, \\ld, \\al_K>0$ ainsi que les valeurs des $z_i$ par la seule observation des $x_i$. \n", "\n", "1. Ecrire une fonction `echantillon(K,alpha,mu,sigma,n)` permettant de simuler $n$ réalisations indépendantes d'une telle loi. $K$ est le nombre de gaussiennes, *alpha* le vecteur des paramètres du mélange, *mu* les moyennes et *sigma* les écart-types (ce sont donc des listes de $K$ élements). La fonction doit renvoyer à la fois les Z_i et les X_i, par exemple sous forme de deux listes. Essayer de coder cette fonction SANS BOUCLE FOR.\n", "\n", "2. Tester cette fonction avec K=3, alpha=(.4,.3,.3), mu=(-4,4,0), sigma=(1,1,1) et n = 1000 et tracer sur le même graphe l'histogramme des réalisations ainsi obtenues et le mélange correspondant." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1oklEQVR4nO29eZhU1bX3/1nQgCAyRFtRQBsQxWZQoWmgBUFwAMSg8V7n4ZqBcNU4xTcxufnd++a+r7l5ExOveeOrj3Ocx0SRoKAMCo1ANzSDUExBGloRGhBBZe71+2NXSdN00ae7T9UZen2e5zynqs4+e686Xf09+6y991qiqhiGYRjxpVnQBhiGYRiZxYTeMAwj5pjQG4ZhxBwTesMwjJhjQm8YhhFzcoI2oDZOOOEEzcvLC9oMwzCMyLBw4cKtqppb27FQCn1eXh6lpaVBm2EYhhEZRKQ83TFz3RiGYcQcE3rDMIyYY0JvGIYRc0zoDcMwYo4JvWEYRswxoTcMw4g5JvSGYRgxx4TeMAwj5pjQG4ZhxBwTeqNW8vJAxNtm0SoMI9yEMgSCETzl5eA1+ZhIZm0xDKNxeOrRi8hoEVklImtF5L5ajvcSkY9EZK+I3FvjWAcReV1EVopIQkSG+GW8YRiGUTd19uhFpDnwMHARUAGUiMgkVV1Rrdh24A7g8lqqeAh4V1X/SURaAm0abbVhGIbhGS89+kJgraquU9V9wMvA+OoFVHWLqpYA+6t/LiLtgPOBJ5Pl9qnqDj8MNxrJwYOwcCGsWOHdR2MYRiTxIvSdgY3V3lckP/NCd6ASeFpEykTkCRE5traCIjJBREpFpLSystJj9UaDWLUKzj4bCgqgd28YMwa2bw/aKsMwMoQXoa9tqM1rFzAH6A88oqrnAl8DR/j4AVT1MVUtUNWC3NxaY+cbflBZCZdcAlu2wNNPw+9/DzNnwve+BwcOBG2dYRgZwMusmwqga7X3XYDPPNZfAVSo6vzk+9dJI/RGlrj1Vti8GWbPdj16gJNOgptugt/9Dn75y2DtMwzDd7z06EuAniLSLTmYeg0wyUvlqvo5sFFEzkx+NApYcZRTjEzy4Yfw+utOzFMiD3DjjXD55fCb38CmTYGZZxhGZqhT6FX1AHA7MBVIAK+q6nIRmSgiEwFEpJOIVAD3AL8SkYrkQCzAT4AXRGQpcA7wmwx8D8ML//Vf0KkT3Hvvkcd+/3vYswcefDD7dhmGkVFEQzjjoqCgQC1nrM8kEpCfD//rf8GvflV7mauvhqlTYeNGpN1x9VowFcKfkWE0KURkoaoW1HbMQiA0Ff7v/4VWreDHP05f5q674Msv4ZVXsmaWYRiZx4S+KbBvnxPvK6+Eo81oGjwYzjgDXnghe7YZhpFxTOibAu+/7+bJX3vt0cuJwA03wKxZdGVDdmwzDCPjmNA3BV56CTp2hIsvrrvs9dcDcA0vZ9gowzCyhQl93Nm/HyZNgiuugJYt6y7fvTucey6X8XbmbTMMIyuY0Medjz6CnTth3Djv51x2GUXMhW3bMmeXYRhZw4Q+7rz7LuTkwMiR3s8ZN47mVME772TOLsMwsoYJfdx5910oKoL27b2fM2AAm+gEkydnzi7DMLKGCX2c+fxzKCuD0aPrd16zZrzHRTBjhq2EMowYYEIfZ6ZNc/v6Cj0wkwtcpMsVFprIMKKOCX2c+eADN63y7LPrfeosRrgXM2f6a5NhGFnHhD7OzJ4NQ4dCs/r/mdfTDU47zYTeMGKACX1c+fxzWLMGhg1reB0jRringqoq38wyDCP7mNDHlTlz3L4xQn/BBW4u/ccf+2NThMjLcxEhvGx5eUFbaxhHx4Q+rsyeDW3aQP/+Da8jdZP46CN/bIoQ5eVuwpGXrbw8aGsN4+iY0MeV2bNdNEovYQ/S0a2bi3bZBIXeMOKEJ6EXkdEiskpE1orIETlfRaSXiHwkIntF5Ij0RSLSXETKRMRW4GSDb76BpUvdQqnGIOJuFvPm+WOXYRiBUKfQi0hz4GFgDJAPXCsi+TWKbQfuAB5IU82duDSERjYoK4ODB6GwsPF1DRkCq1a5MMeGYUQSLz36QmCtqq5T1X3Ay8D46gVUdYuqlgD7a54sIl2AS4EnfLDX8EJJidsPHNj4ugYPdvsFCxpfl2EYgeBF6DsDG6u9r0h+5pX/Bn4GHHWOnohMEJFSESmtrKysR/XGEZSUQJcuLhF4YykocPPwzU9vGJHFi9BLLZ95CoAiIuOALaq6sK6yqvqYqhaoakHu0dLdGXVTUuJPbx7guOOgTx/z0xtGhPEi9BVA12rvuwCfeaz/POC7IrIe5/IZKSLP18tCo37s2OEWSvkl9ODqWrjQApwZRkTxIvQlQE8R6SYiLYFrgEleKlfVX6hqF1XNS543Q1VvaLC1Rt2Ulrq9n0Lfv79bOLVxY91lDcMIHTl1FVDVAyJyOzAVaA48parLRWRi8vijItIJKAXaAVUicheQr6o7M2e6USupgdiCAv/qTC26WrQITj3Vv3oNw8gKdQo9gKpOAabU+OzRaq8/x7l0jlbHLGBWvS006kdJCfTsCR06+Fdnv35uQHbRIrj8cv/qNQwjK9jK2LhRVta4sAe10aYN5Oc7oTcMI3KY0MeJL7+E9esbFH++Tvr3N6E3jIhiQh8nli51+0wJ/aZNbjMMI1KY0MeJJUvcPlNCD841ZBhGpDChjxNLlsDxx8Mpp9R6uD4x1k87rcbJ55zj9gvrXPtmGEbI8DTrxogIS5a43rzUtpj5UIz1BnHccW42T+qpwTCMyGA9+rhw8KDLBJUJt02Kvn1h2bLM1W8YRkYwoY8La9bA7t2ZF/o1a1y8e8MwIoMJfVzI5EBsin79nO9nxYrMtWEYhu+Y0MeFJUsgJwfOOitzbfTt6/bmvjGMSGFCHxeWLHEi36pV5tro3h1atz40X98wjEhgQh8Xli51rpVM0ry5i01vPXrDiBQm9HFg506oqIDevTPfVt++1qM3jIhhQh9x8vKgsP1KAC7/5Vn1WwTVEPr2hcpK2LzZh8oMw8gGJvQRp7wcFvzFCf2bK89ClbTb+vU+NJhyD5n7xjAigwl9HEgkoEUL6NEj822lZt6Y+8YwIoMnoReR0SKySkTWish9tRzvJSIficheEbm32uddRWSmiCREZLmI3Omn8UaSRMKFJ8jJQkSL3Fw46STr0RtGhKhTGUSkOfAwcBEuUXiJiExS1eqrZrYDdwCX1zj9APBTVV0kIscBC0XkvRrnGo0lkTjU084G/fpZj94wIoSXHn0hsFZV16nqPuBlYHz1Aqq6RVVLgP01Pt+kqouSr3cBCaCzL5YbALRgH/zjH5ldKFWTvn3d6tgDB7LXZhhocEQ4wwgWL0LfGdhY7X0FDRBrEckDzgXmpzk+QURKRaS0srKyvtU3WXqyxgU0y6bQ9+kDe/bAJ59kr80g+ewzuPFGaN/ebRMmwPbtQVsVXr7+Gu67D7p1gy5d4M477XoFjBehry3mbb26NiLSFngDuEtVd9ZWRlUfU9UCVS3Izc2tT/VNmrNIJF9kUejz892+KcS8Wb0aBg6E11+Hq6+G730PnnkGCgrg00+Dti58bNsGQ4bA737nnvzOOw8efhgGD4YNG4K2rsniRegrgK7V3ncBPvPagIi0wIn8C6r61/qZZ9TFWSTcJPkzz8xio8mbStyFftcuGDcO9u+HBQvg8cfh6afhgw9g61Z3bO/eoK0MDwcPwhVXuJvjlCkwaRK88oq7Xps3w+WXuwirRtbxIvQlQE8R6SYiLYFrgEleKhcRAZ4EEqr6x4abaaTjLBJuJVSbNtlrtF0790ged6H/2c/c+Mfrrx8+2D1kCLz4IixeDL/+dWDmhY4HHoDZs90NcfToQ5+fdx688IJLQ/kf/xGcfU0ZVa1zA8YCq4F/AP+W/GwiMDH5uhOu578T2JF83Q4YinPzLAUWJ7exdbU3YMAANbyxiHNUx4zJfsOXXKLav7+quuVYcaMPS1WbNVO98870hW6+WbVFCz2NT7JkVYjZtEn12GNVL79ctaqq9jK33KKak6OaSGTXtiYCUKrpNDzdgSA3E3qPHDyoX9Na9Z57st/23Xertm6tevBgLIV+CqNVO3ZU3bYtfaGNG1VbtdKnuTlrdoWWW291Ir56dfoyW7a4m8H112fPribE0YTeVsZGmfJy2rA7uwOxKfLznb/Vl7gKIWPJEsbwrnPdfOc76ct16QK3385NPOtcPE2VTZucu+aHP3QL99KRmwu33govveQylRlZw4Q+yiQCmHGTIs4zb/7wB3bRFiZOrLvs3XdzkObw//5f5u0KK4884tZU/PSndZe95x5o2RIefDDzdhnfYkIfZYIU+rjOvNmyBV56iSf5AXToUHf5zp15nX+CJ5+Er77KuHmhY88eePRRNwPp9NPrLt+pE1x1FTz/fNO8XgFhQh9lEgk2c+LR3QuZomNHOPnk+An9iy/CgQM8xgTPp/yZ2+HLL+HVVzNoWEh5+20Xtvr2272f8+Mfu6mrr7ySObuMwzChjzKJBAkC6M2n6N07fkL/zDMwcCAJ8j2fMpciFzn0hRcyZ1dYee45OOUUGDXK+zlDhrjfzpNPZs4u4zBM6KOKavBCn5+fFPqYxIBZvNjl3r355nqeKHDDDTBzZtNaLbt1K7zzDlx3nUsz6RURuP56+OgjWy2bJUzoo8qWLfDFF6ykV3A25OfD119zKtH5Z83LS5+B64Fzn2cfLTj+9mvqn43r+uvdzfellzJhdjh55RU3CHvjjfU/96qr3P611/y1yagVE/qokhyIDbxHD+QTHfdNeXmaDFxVyr3d/0bLsRexTY+v/6zRnj2hsND5+JsKL7/sXDANSUrfowcMGGB++ixhQh9VTOj95eOPYd06F4+lnpx2mnsi+NmCK6GsjFNlQ9qnhrw83y0PhspKmDsXrryy4XVcfTWUlDSdKKgBYkIfVRIJaNuWCroEZ8Pxx8OJJ8ZD6P/2N6fEl11W71PXr3dPBb9b6dI0bPjz22nz9paX+2x3UEyeDFVVMH583WXTkbpJTPIUOstoBCb0USWRgF69qD2KdBbp3TseQv/mm242SKdODa/jzDPhjDPgrbd8MytIjjae8db332QDXZEB5zb8SaV7d/cb/vvffbbcqIkJfVRJJIJZKFWT/Hwn9FHOvlRe7iIrNsBtcwTjx8OsWW5efcRJO57x9TeMb/0ep94+HlVp3JPKpZe6MMa2eCqjmNBHkZ073TS+kAh9e3ZGe1rh5Mlu3xg3RIrx4138+qlTG19XWJk508U5aoCb6wguvRT27YP33298XUZaTOijyMqVbh8SoQeivXBq2jSX9u6MMxpf16BBLt3ge+81vq6wMm0aHHMMnH9+4+saOtTlNzD3TUYxoY8iQca4qUnUhX7/ftdDvegif+rLyYGRI50YRtmddTTeew+GD3di31hatHDXfurU+F6vEGBCH0USCfcP0qNH0JZAbi5bOR6WLw/akoaxYIGLu+KX0IOra8OGeIbi3bjR/f78vF6jRrl6m3Ko5wxjQh9FEgm3QCcnJ2hLQITl9D70lBE1pk2DZs1cL9wvLr7Y7ePovkl9p9R39INUnJzp0/2r0zgMT0IvIqNFZJWIrBWR+2o53ktEPhKRvSJyb33ONRpAWGbcJFlBvuvRR/HR+733oKDA3wigPXo4n/+0af7VGRamTXNTUPv08a/Onj2hc2eYMcO/Oo3DqFPoRaQ58DAwBsgHrhWRmqH9tgN3AA804FyjPuzd6x5xwyb0O3bA558HbUr9+PJL57rx0w2R4qKLnO9//37/6w6Kqio3O+bii91ker8Qcb36GTNcG4bveOnRFwJrVXWdqu4DXgYOm4emqltUtQSo+auu81yjnqxZ4/4ZQiT0y+mdfBExP/3MmXDwYOaEftcuKC31v+6gSCRg2za44AL/6x450kXD/Phj/+s2PAl9Z2BjtfcVyc+84PlcEZkgIqUiUlpZWemx+iZImGbcJFlBRGfezJgBrVu7FbF+k5p6+OGH/tcdFHPnuv155/lfd2qMxPz0GcGL0Nf2jObVGev5XFV9TFULVLUgNzfXY/VNkJUr3aPumWcGbcm3bOYk5+OOmtAXF7t57y1b+l/3iSe65f1xEvriYjjhBG8pA+tL166u3pkz/a/b8CT0FUDXau+7AJ95rL8x5xq1kUi4cIlt2gRtSTXEzaePkuvmq69ckpFM9E5TDB8Oc+Y491AcmDsXior89c9XZ9gw14b56X3Hi9CXAD1FpJuItASuAbyGm2vMuUZthGzGzbfkR2zmzYIFToAzLfQ7d7obStSprHTjQ2muVypUc2O27z89FLZt46zmq+IV0jkE1Cn0qnoAuB2YCiSAV1V1uYhMFJGJACLSSUQqgHuAX4lIhYi0S3dupr5M7KmqglWrwin0vXvDF1/A5s1BW+KN4mKnLpnwz6dI+ek/+CBzbWSLjz5y+6KiWg+nQjU3Zntq1VAAEo8XxyukcwjwNI9eVaeo6hmq2kNV709+9qiqPpp8/bmqdlHVdqraIfl6Z7pzjQZSXu6CSYVR6KMWCqG42N2cOnTIXBudO7s59XEQ+rlz3WrsAQMy10bPnpCb69xdhq/YytgoEcIZN9/SO0JTLKuqXA81Te/UV84/H2bP/tbvXB8XR6hcF8XF0L+/m6WUKURckDMTet8xoY8SYRb6Tp1c7zgKPfrly53vPJP++RTDh8P27d/eAOvj4giN62LfPpfyLxs3xqFD3YLATZsy31YTwoQ+SiQSbtqen8v1/UKSM2+iIPTFxW6fDaEfNuzwNqNIWZlbkZ2N6zXU+ekjfb1CiAl9lAjrjJsUvXtHY+ZNcTGcdJJLZZdpunVzbaUWG0WRlO2ZHLhOce65zj1k7htfMaGPCqrhF/r8fLdEPuwrm4uLMzsfvDqpmT2pWStRZO5cN2BwyimZb6tFC7eIzXr0vmJCHxW2bHHTF8Mu9BDqAdlObIJPPsmOGyJFURGsXev+hlFD9dBCqWwxeDAsXswx7M5emzHHhD4qpAZie/UK1o6jkZp5E2I/fREZjNeSttGkSEaxV19eDp99lt3rNXgwHDjAuZRlr82YY0IfFcI84ybFKae4/J8hFvrzKHYp8Pr3z16jAwY4l0QU/fQpm7PZox80yO2Yn702Y44JfVRIJKBtW+jSJWhL0iNyaEA2pJxHMQwcmJlAZulI3ViiKvRt2/qbaKQuOnWCU09lMPOy12bMMaGPComEc9tkYwCxMYR5iuXu3fRnUXZ7pymKilxs+n37st92Y5g71/Wws522cvBg69H7iAl9VAj7jJsU+flu1k0YZ96UlNCCA9n1N6coKoI9e2Dx4uy33UCOJRnhM4gb46BB5FEevaxlIcWEPgrs3AmffhoNoQ/zgGxqyl5QPXqIlPtmEPNd6IYgboxJPz3zrVfvByb0UWDlSrePgtCHObhZcTEJesHxx2e/7VNOcYFuIiT0Rcx1rsKU6GaT/v3ZT44JvU+Y0IeQvLzDg1vdPMjNuDnzirOOCHx12mnB2noEXbrAcceFb0A2GcismAB6pykitnCqiLmZj/CZjtatWcLZMM8GZP3AhD6ElJcfHtzqLz9PQIsWrNrf44jAV+vXB21tDcIa82bVKti+nbkE4LZJUVQEFRWwcWPdZYOmqoohZCnCZxrmM8gFU4tLhq4AMaGPAitWuByx2Z750FDCmFYw6Z8PtEcfJT99IkEHvgxU6Ocx2KV8DFunIYJ4EnoRGS0iq0RkrYjcV8txEZE/JY8vFZH+1Y7dLSLLReRjEXlJRI7x8ws0CVasiIZ/PkXv3m65/9atQVtyiGRi69WcEZwN/fq5gF1RcN9kM8JnGuZjA7J+UafQi0hz4GFgDJAPXCsi+TWKjQF6JrcJwCPJczsDdwAFqtoHaI7LG2t4ZfduF5slv+YlDzFhHJBNBTIjwHUILVq4xVpR6NHPncsWcl2GrIBYQ0/o2NGE3ge89OgLgbWquk5V9wEvA+NrlBkPPKuOeUAHETk5eSwHaC0iOUAb4DOfbG8arF7tBhKjJPR9+7r90qXB2pEildg6QDfEtwwZ4uK77w55wK65c914RqAL9AQKC21A1ge8CH1noProUUXyszrLqOqnwAPABmAT8KWqTqutERGZICKlIlJaGcbFNkGR6hVHyXXTubNLjhIWoZ8bQCCzdBQVwYEDsHBh0JakJ3ljDHTgOsWgQe5/4KuvgrYk0ngR+tpu6TUzS9RaRkQ64nr73YBTgGNF5IbaGlHVx1S1QFULcnNzPZjVREgkoFkzOCNA33IdHJEHtZkwc3s/5j2+NBw5UIuLXWybgoKADKjG4MFuH2Y/fdK20Ah9VVW4b4wRwIvQVwBdq73vwpHul3RlLgQ+UdVKVd0P/BXC8OuJECtWwOmnQ6tWQVuSltryoF5wRz8Gt1mGHqwKPgdqcbGLIHlMCOYBnHii83uH2U9fXAwtWlBKCG6MAwe6vfnpG4UXoS8BeopINxFpiRtMnVSjzCTgpuTsm8E4F80mnMtmsIi0EREBRgEJH+2PP1GbcZPi7LPhm29g3bpg7di71wUTC4PbJkVq4VRYUy7OnQsDBrCXENwYc3NdyscFC4K2JNLUKfSqegC4HZiKE+lXVXW5iEwUkYnJYlOAdcBa4HHg1uS584HXgUXAsmR7j/n9JWLL/v1uEDFKA7Ep+vVz+yVLgrVj4UIXMTIMA7EphgyBzZtDuNoNd61KSsJ1vQoLrUffSDzNo1fVKap6hqr2UNX7k589qqqPJl+rqt6WPN5XVUurnfsfqtpLVfuo6o2qujczXyWGrF3rBu6iKPT5+W5sIegB2SADmaUjzBmnysrcU1CYrldhoVtR/JlN2GsotjI2zKRm3ERR6Nu0gZ49wyH0p58OJ50UrB3V6dMHjj02nH76IDJK1UUqqFpJSbB2RBgT+jCTEvozzwzWjoZy9tnBCn0qsXWY/PPgQlkUFoazR19cDN26wckn1102W5x7rrtm5r5pMCb0YSaRcHMSjz02aEsaRr9+bjB2165g2l+71s0JD5vQg/PTL1kCX38dtCWHUK22gjhEtG7tfks2INtgTOjDzIoV0XTbpEgNyC5bFkz7YfTPpygqclEZS0vrLpstypMZncJ4vQoLneumqipoSyKJCX1YOXjQJRyJ4tTKFCmhD8p9U1zsYqmH8RqGceFUGP3zKQYNcpnWVq0K2pJIYkIfVj75xM1+iHKP/tRToX37YIW+qMjN/gkbxx/vVjuHaUB27lxo2/ZQrKIwUVjo9ua+aRAh/A8wAOefh2gLvYgbkC0ry37bW7e6azhsWPbb9krYFk4VF7snjebNg7bkSHr1cpnLbEC2QZjQh5VU4o5evYK1o7EMGOAGHQ8cyG67qZ7y0KHZbbc+FBW5G9I//hG0JW7AfOnScLptwD2VDRxoPfoGYkIfVj7+2OVfDSJfp58MGOBC8iayHPli9uzwBDJLx5Ahbh8GP/2CBW6gM6xCD859s2RJ+EM8hxAT+rCybFk4faX1JSW02Z5dMmeO6wGGIZBZOvLznTsiDH764mLnakstTgojgwa5J8PFi4O2JHKY0IeQHPa7HnAchL5nTydm2Qwz+803rr0w++fB+cIHDQpHj7642P3ewvwEaQOyDcaEPoScwWoX0CwOQt+smVvZmE2hLylx1y/M/vkUQ4a4p7egFpWB6yXPnRv+63XKKc6daQOy9caEPoT0JbnAKA5CD859s3gxzcnSgOycOW4fZn9ziqIi5xsPMo7LsmUug1PYhR5cr9569PXGhD6E9GWZe6yP+oybFAMGwJ495JOlZOFz5rjAYR07Zqe9xpDyiQfpvkndGKMg9IMGuVlKW7cGbUmkMKEPIX342C2mCXFWqXoxYIDbkQX3zcGDzg0Rdv98io4d3crdIAdk58xxi9u6dq27bNCk/PQWybJemNCHkL7EZMZNiuSAbFaEftkyt1Q+Cr3TFEOGwLx5wSycUnVCH5XrVVDgxn3MfVMvPAm9iIwWkVUislZE7qvluIjIn5LHl4pI/2rHOojI6yKyUkQSIjLEzy8QO3btojufxEvomzWD/v0ZSBZ6YbNnu30YI1amo6gItm+H1auz3/b69S6hR1SuV9u2blqqDcjWizqFXkSaAw8DY4B84FoRqbkufwzQM7lNAB6pduwh4F1V7QWcjeWMPTqpFbFxEnqAQYM4h8WwZ09m25k1y4V2Pu20zLbjJ6mFU0G4b6Lkn0+RGpANS+iICOClR18IrFXVdaq6D3gZGF+jzHjg2WRKwXlABxE5WUTaAecDTwKo6j5V3eGf+TFkWcxm3KQoKqIV+2DRosy1UVUFM2fCyJGZayMT9Orlgpx9+GH2254zxwWe6907+203lEGDYNu24BPPRwgvQt8Z2FjtfUXyMy9lugOVwNMiUiYiT4hIrVk0RGSCiJSKSGllZaXnLxA7li3jK451vdI4kY1e65Il8MUX0RP6Zs1g+HD3NJJt5sxxrqMwBjJLhy2cqjdehF5q+azmM1O6MjlAf+ARVT0X+Bo4wscPoKqPqWqBqhbk5uZ6MCumLFvGx/QJZ2jdxnDiiazh9MwK/YwZbn/BBZlrI1OMGOH85evXZ6/N7dtdcpsouW3ATZ1t3dqEvh54UZMKoPq8qy5AzXTs6cpUABWqmho5eR0n/EZtqB4S+hgylyIn9Jnyrc6Y4dwgp5ySmfozyfDhbv/BB9lrMzVwHVKhP+00F37niK1FDrN3D2Duf8//9rO4PQD7jRehLwF6ikg3EWkJXANMqlFmEnBTcvbNYOBLVd2kqp8DG0Ukld16FGRr1UwEqaiAbdso49ygLckIcymCzZtdUhW/2b/f+bhHjfK/7mzQpw985zvZFfoZM1zPOKSBzNavd32C2rZhPx1EUatF6L79qLosiEZ66hR6VT0A3A5Mxc2YeVVVl4vIRBGZmCw2BVgHrAUeB26tVsVPgBdEZClwDvAb/8yPGcmofIs5J1AzMsVckiEJMuG+KSlxy/ij5p9PEYSffuZM15uP4sK8wkKXgS2o7GURw5MjWFWnqOoZqtpDVe9Pfvaoqj6afK2qelvyeF9VLa127uKk772fql6uql9k5qvEgLIyEGEJZwdtSUZYQT60a3coabdP5OXBr86bQRXC8VcOr/1xP7mFetbliBHuaScb3dMtW9wMryiOZ8ChpxDz03siZiN+EaesDHr25GvaBm1JRqiiuetB+txrLS+H/z10Ks36n8s2PT7t475qdsc6682IEW6fDfdN6m8Q1SegU0+FE080ofeICX2YKCtzIX3jzMiRsHKlW43pEx34wrmDLr3UtzoDIeWnz4b7ZuZMlycgGYcocog4942tkPWECX1Y+OIL1zU955ygLcksqR5kaiqkD1zMNLdYauxY3+oMhJSffubMzK/6nDHDtZWTk9l2MsmgQa7T8OWXQVsSekzow0IqPVrce/Rnn+0iNvoo9GOZ4laWDhzoW52BcdFFzr+0Zk3m2vj0UxdXJ6r++RSFhe6GmO00lRHEhD4slJW5fdyFvlkzJzDTp/vTa62qYgzvwOjR0VrdmY5LLnH7d9/NXBvTp7t9VP3zKVI3dvPT14kJfVgoK3MLfU48MWhLMs/IkbBhg0sg0VhKSzmRyui7bVJ07+5yEWRS6N95Bzp1gn79MtdGNujY0V2refOCtiT0mNCHhbKy+PvnU4we7fZ//3vj65o0iYM0O9QTjgOjR8OsWRzDbv/rPnAApk6FMWPiEWajKLna+oioLEZ1YvCXjgG7d7tBpbi7bVL06OGyKr39duPqUYXXXmMmFzgffVy45BLYvZthzPa/7vnz3cD/mDH+1x0Ew4bB1q30YmXQloQaE/owUFbmUuDFYTDRK5dd5uaLN2bGxLJlsHo1r/HP/tkVBoYPh1atuISp/tc9ZYoby7joIv/rDoJkysihzAnYkHBjQh8GUoNJTU3oU26EhvLaa9CsGX/jCv/sCgPHHgvnn+8Gmf3mnXdcNqkOHfyvOwhOPx1OPDEzTz8xwoQ+DJSUQOfO0Yy62FCGDHHulrfeatj5SbcNw4dTSQwHsMeOJZ8ErF3rX52ffeaeHuPitgG3cGrYMBP6OjChDwMLFhxKptBUaN4crrgCJk2Cb76p//nz58OqVXDttf7bFgauSD6lvPGGf3W++abbX3aZf3WGgWHD6MZ6F/3VqBUT+qDZvt312pqA0NeML37BE9fBV19xzbGTjgg+Vmd88SefhDZt4Oqrs2F69jntNEoo8CT0eXnpg7hV36bf9gYJeiF9ensqH+oAcNVJxdOfY376dJjQB01Jids3AaGvGV985sHh0LkzL4974YjgY0cN4PjVV/Dyy3DVVS4aZkx5gyvd72PDhqOWKy9PH8Tt221LJaOazeKsX/1T3WWjEACuOmefzS7aHkqkYhyBCX3QLFjguk9RDS7VGJo1g+uvdwOE9XnsfuUVJ/Y/+EHmbAsBb3Cle/HXvza+srfecvGArryy8XWFjZwcl+vAevRpMaEPmpISl/6uffugLQmGH//YCdCjj3orX1UFf/gD9O3rZo/EmLX0dN/z9dcbX9mrr7r1C2fHM9fBHIa66bY7dgRtSijxJPQiMlpEVonIWhE5Irl3MoXgn5LHl4pI/xrHm4tImYhM9svwWKDaNAdiq9O9O4wbB489Bnv21F1+8mRIJODnP3dPQnHnuutcopbGBDnbsAHefx9uuCG212w2w9z/k89JbeJCnUIvIs2Bh4ExQD5wrYjk1yg2BuiZ3CYAj9Q4ficuDaFRnQ0bXA7VpjR/vjbuuQcqK+GRmj+bGlRVwa9/7UYJ4zoIW5ObbnIurmeeaXgdzz3nRPDmm30zK2zMZxC0bJndnLsRwkuPvhBYq6rrVHUf8DIwvkaZ8cCzyZSC84AOInIygIh0AS4FnvDR7niQyp06ZEiwdgTNiBFupeb99x99pewzz8CiRfCb30Q7jnp9OOUUF/vmL39xq6fri6q7biNGQLduflsXGvbQ2v0fpSJzGofhReg7Axurva9Ifua1zH8DPwOqjtaIiEwQkVIRKa2srPRgVgyYMwfato1+FEE/+K//clNNf/GL2o9//jncd58LYhXXufPpuOUWF0O+IREtZ81y03f/5V/8tip8jBrlFoRt2xa0JaHDi9DX5tSrGSqu1jIiMg7YoqoL62pEVR9LJhEvyM3N9WBWDJgzx/VCmkrv9GgMGAB33+3cN6++evix3budq+arr5wvP6Z+5rR897tu5fQf/1j/cx98EE44wU1FjTsXXpictzszaEtChxehrwC6VnvfBaiZ8DNdmfOA74rIepzLZ6SIPN9ga+PEjh1ulkBqsYfhXDfDhsENN3AXD8LXX8OKFS6a4+zZ8MQT0Lt30FZmn5Yt4c47XVauhXX2mQ6xerUbvP7Xf4XWrTNnX1gYONDlwTX3zZGo6lE3IAdYB3QDWgJLgN41ylwKvIPr2Q8GFtRSzwhgcl3tqSoDBgzQ2DNliluXMn36EYcgAHvCwhdfqI4de/janeOOU33xxbSnxPV6Hfa9duxQbddO9aqrjl6uOjffrNqqleqmTZkwL1R8ew3GjVM9/fRAbQkKoFTTaGqdPXpVPQDcDkzFzZx5VVWXi8hEEZmYLDYleTNYCzwO3OrLXSjOzJnj4r0MGhS0JeGiQweYPJmRTIf//E/4859dz7Sp+eVr0r493H67c2t5yZH68cfw7LPunE6dMm9fWBg1yo1J1LGauKkhmuls8w2goKBAS+Oe8Hf4cOd7riXfpYg/6VSjTH2uQVyv1xHfa+dO6NnTbR9++G2GqCPKqcLFF7vf1rp18UrKkoZvr8GyZW5yw1NPuUHsJoSILFTVgtqO2crYINizx/0Txnxlp+Ez7dq52UnFxfDQQ+nLPf64WyD12982CZE/jD59XN5l89Mfhgl9EHz0kRP7Cy4I2hIjatxyC4wf71YG17Y4aP58uOsuNwNl4sQjj8cdEee+ef99t8DOAEzog2H6dOefHz48aEuMqCHi3BI9esCll8LzzyNUOb/Fq6+6hWcnnwwvvND0pqGmGDPGrThftChoS0KDCX0QTJ/upoI11UBmRuP4znfcVMt+/eDGG9nEyW4F7dVXO//9Bx8490VTZfRod5P7+9+DtiQ0mNBnm507XcTKUaOCtsSIMief7AZkX3qJKYx1g6/PPutcN126BG1dsOTmwuDBbg2BAZjQZ5W8PLis/Qdw8CAX3D8q+pl9jGDJyYFrruH7PO1i4dx4o62yTnHppW4a6uefB21JKDChzyLl5fD2ndPhmGOYuXtI9DP7GEZYGTfO7d95J1g7QoIJfbZ57z0X9uCYY4K2JNTUzC8bi9ym9cTrNYjr928U/fo5F5b56QEX3sDIEnl84mK3/PCHQZsSeuypxq5BoxCBsWPhpZdg715o1SpoiwLFevRZ5FKSvYvUY6VhGL5Q29PPmMcuh127GHfMe4d9npcXtLXZx4Q+i4xjMpxxhpsCZxiGb6xff+RY1zt7R0GHDky+6bXDPi8vD9ra7GNCny2++ooLmGm9ecPIFi1bwuWXw1tvOfdNE8aEPltMn04r9pnQG0Y2+ed/dukp33svaEsCxYQ+W0yaxJe0s0QjhpFNLrzQhb5+7bWgLQkUE/pssH8/vPkmb3MZtGgRtDWG0XQw9w1gQp8d3n8ftm/nFa4O2hLDaHpcfbVz3zThOfWehF5ERovIKhFZKyL31XJcRORPyeNLRaR/8vOuIjJTRBIislxE7vT7CwRNXl7dC1qeGfsKO2jPqlMvDtpcw2h6XHihiw309NNBWxIYdQq9iDQHHgbGAPnAtSKSX6PYGKBncpsAPJL8/ADwU1U9C5dL9rZazo005eW1hzH4dtuzl39p/yYdbr6c1eVNe9GGYQRCTg7cdJMLh9BEY9946dEXAmtVdZ2q7gNeBsbXKDMeeDaZo3Ye0EFETlbVTaq6CEBVd+Fyznb20f7wM3Wqe2y82tw2hhEYt9wCBw/Cc88FbUkgeBH6zsDGau8rOFKs6ywjInnAucD8elsZZZ59Fk44wT0+GoYRDGeeCUOGJN03MUwwXAdehL62NDU1r9RRy4hIW+AN4C5V3VlrIyITRKRUREorKys9mBUBtmxxo/033WSzbQwjaH74Q0gkGE4tKRhjjhehrwC6VnvfBfjMaxkRaYET+RdU9a/pGlHVx1S1QFULcnNzvdgefp57Dg4cgB/8IGhLDMO49lo4/nju5CiJ1WOKF6EvAXqKSDcRaQlcA0yqUWYScFNy9s1g4EtV3SQiAjwJJFT1j75aHnZU4Ykn3ONifqzGnw0jmrRuDT/+MeN5Cz75JGhrskqdQq+qB4Dbgam4wdRXVXW5iEwUkVSa+SnAOmAt8Dhwa/Lz84AbgZEisji5jfX7S4SS6dNh5UqYMCFoSwzDSHHrrVTRDP7856AtySqiGr6BiYKCAi0tLQ3aDE+IuM77EYwd67LQl5c3+VjYhhEmXpTruO64ya5Xf/zxQZvjGyKyUFULajtmK2MzQSLh5uzedpuJvGGEjGdO/iVVu77if5/wYJ2LHeMSu96EPhP8/vcuVeDEiXWXNQwjq0z7rA/N/vmf+NVxf0K3bjvqgse4xK43ofebNWvc3PmJEyEus4cMI278+7/Drl3whz8EbUlWMKH3m//8Txcx7+c/D9oSwzDS0acPXHcdPPhgk0jOa0LvJ4sXw4svOt98p05BW2MYxtH4P/8HmjWDe+8N2pKMY0LvF6rwk5/Ad74Dv/xl0NYYhlEXXbrAL34Bb7zhQonHGBN6v3jxRZgzB37zG+jYMWhrDMPwwr33wumnw49+BDtrjc4SC0zo/WDTJrjjDigshO9/P2hrDMPwyjHHwF/+Ahs2wN13B21NxjChbzTqgiV9842bbdO8edAGGYZRH4qK3OSJp56CV14J2pqMYELfSO7jtzBlips7f+aZQZtjGEZD+J//E847z8WtX7QoaGt8x4S+Mbz9Nvfzby4q3m23BW2NYRgNpWVLNyh7wgkwfrxz5cQIE/qGMmsWXHUVi+jvolRKbSH5DcOIDCedBJMmuYVUI0bAxo11nhIVTOgbwrRpMG4cdO/OaN6FNm2CtsgwDD845xz3/71tG5x/PvksD9oiXzChrw+q8MgjLjLl6afD+++zjROCtsowDD8pLHTz6nfvZp4MYby8FfngZyb0XqmshO99D269FS6+GD78EE4+OWirDMPIBAMHQkkJx53bk7e4HP2XW9AvdkQ2+JkJfV188w088IDrwf/97y4I0uTJ0K5d0JYZhpFJunaFuXPdSvfnnoMePeCPf3SaEDFM6NOxerX7A3ftCv/jf8DQobBkCdxzj4uPYRhG/GnVCu6/H0pLoaAAfvpTFzrh3nth2bI0WYfChyfFEpHRIrJKRNaKyH21HBcR+VPy+FIR6e/13NCwZQu89Rb87GfQr5+bE//b38Lw4fDBB643f9ZZQVtpGEYQnHMOTJ3qXLYXXggPPeR0okcPuOMOruElWLcutMJfZypBEWkOrAYuAipwycKvVdUV1cqMBX4CjAUGAQ+p6iAv59ZGg1MJHjwIe/bA3r2173fsgK1bD20bN7r48WvWuPcALVrA4MHOH3/lla5Hf9TrE9q/rWEYmWLzZjcV8623YMYM2L3bfX7ccc7Ne/rpTjtycw9t7du7GXqtW7utTRsXgqFFC7fl5DhvQQOnah8tlWCOh/MLgbWqui5Z2cvAeKC6WI8HnlV315gnIh1E5GQgz8O5/tG2rRN0LzRv7kIJ9+wJV1zhevCDB8OAAe7iG4ZhpOOkk1wgtB/9CA4c4JwWH7P40fmwfDmsXetClk+efOgGUJ96P//cd3O9CH1noPrKgQpcr72uMp09nguAiEwAJiTffiUiqzzYVhsnAFvrLHXwIHz6qdtmzWpgUw6PN2BvdmUfs6t+mF31o8nYJX5kDt28+QREGmrXaekOeBH62mSsprMiXRkv57oPVR8DHvNgz1ERkdJ0jy9BYnbVD7Orfphd9aOp2eVF6CuA6o7qLsBnHsu09HCuYRiGkUG8zLopAXqKSDcRaQlcA0yqUWYScFNy9s1g4EtV3eTxXMMwDCOD1NmjV9UDInI7MBVoDjylqstFnEdKVR8FpuBm3KwFvgFuOdq5Gfkmh2i0+ydDmF31w+yqH2ZX/WhSdtU5vdIwDMOINrbE0zAMI+aY0BuGYcScWAq9iJwjIvNEZLGIlIpIYdA2pRCRnyRDQiwXkd8FbU91ROReEVERCUXsZRH5vYisTIbV+JuIdAjQllCG8hCRriIyU0QSyd/UnUHblEJEmotImYhMDtqW6iQXdL6e/G0lRGRI0DYBiMjdyb/hxyLykoj4tnIzlkIP/A74taqeA/x78n3giMgFuJXB/VS1N/BAwCZ9i4h0xYWqCFMOtfeAPqraDxdK4xdBGJEM5fEwMAbIB64VkfwgbKmFA8BPVfUsYDBwW4hsuxNIBG1ELTwEvKuqvYCzCYGNItIZuAMoUNU+uMkr1/hVf1yFXoFUHOH2hGfu/r8Cv1XVvQCquiVge6rzIPAz0ixoCwJVnaaqB5Jv5+HWYQTBt2FAVHUfkArlETiquklVFyVf78KJVudgrQIR6QJcCjwRtC3VEZF2wPnAkwCquk9VdwRq1CFygNYikgO0wUfdiqvQ3wX8XkQ24nrNgfQEa+EMYJiIzBeRD0RkYNAGAYjId4FPVXVJ0LYche8D7wTUdroQH6FCRPKAc4H5AZsC8N+4jkNVwHbUpDtQCTyddCs9ISLHBm2Uqn6K06oNwCbcWqRpftXvZWVsKBGR94FOtRz6N2AUcLeqviEiV+Hu3heGwK4coCPuEXsg8KqIdNcszHGtw65fAhdn2obaOJpdqvpWssy/4VwUL2TTtmp4DuURFCLSFngDuEtVdwZsyzhgi6ouFJERQdpSCzlAf+AnqjpfRB4C7gP+vyCNEpGOuKfEbsAO4DURuUFVn/ej/sgKvaqmFW4ReRbnHwR4jSw+PtZh178Cf00K+wIRqcIFV6oMyi4R6Yv7cS0RF52tC7BIRApV1f8weh7tqmbfzcA4YFQ2bohp8BIGJDBEpAVO5F9Q1b8GbQ9wHvDdZPjyY4B2IvK8qt4QsF3g/pYVqpp66nkdJ/RBcyHwiapWAojIX4EiwBehj6vr5jNgePL1SGBNgLZU502cPYjIGbhYQIFG9lPVZap6oqrmqWoe7h+hfzZEvi5EZDTwc+C7qhpk/rbQhvIQd3d+Ekio6h+DtgdAVX+hql2Sv6drgBkhEXmSv+uNInJm8qNRZCpsev3YAAwWkTbJv+kofBwkjmyPvg5+BDyUHNTYw6Hwx0HzFPCUiHwM7ANuDrCXGgX+DLQC3ks+bcxTVT+CwdaLgEJ5eOU84EZgmYgsTn72S1WdEpxJoecnwAvJm/Y6kiFbgiTpRnodWIRzU5bhYzgEC4FgGIYRc+LqujEMwzCSmNAbhmHEHBN6wzCMmGNCbxiGEXNM6A3DMGKOCb1hGEbMMaE3DMOIOf8/WjlI3LmMubMAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Simulation des Données \n", "K=3\n", "alpha_real=[.4,.3,.3]\n", "mu_real=[-4,4,0]\n", "sigma_real=[1,1,1]\n", "n=1000\n", "\n", "Z,X=echantillon(K,np.array(mu_real),np.array(sigma_real),alpha_real,n)\n", "\n", "plt.hist(X,bins=round(3*n**(.3)),density=1,histtype=\"step\",color = 'b') # remplacer density par normed si matplotlib <2.2\n", "\n", "x=np.linspace(-8,8,num=1000)\n", "y=densite_theorique(K,mu_real,sigma_real,alpha_real,x)\n", "plt.plot(x,y,\"r\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimation des paramètres dans le cas où l'on connaît les $z_i$\n", " L'estimateur du Maximum de Vraisemblance (EMV), étant données les observations $(z_i,x_i)_{1\\le i\\le n}$, fournit les estimations suivantes pour le paramètre $\\tta=(\\al_k,\\mu_k,\\si_k)_{1\\le k\\le K}$: pour tout $k=1,\\ld, K$, \n", " $$\\widehat{\\al}_k=\\f{|\\{i\\ste z_i=k\\}|}{n},$$ \n", " $$\\widehat{\\mu}_k=\\ff{n\\widehat{\\al}_k}\\sum_{i\\ste z_i=k} x_i,$$ \n", " $$\\widehat{\\si}^2_k=\\ff{n\\widehat{\\al}_k}\\sum_{i\\ste z_i=k} (x_i-\\widehat{\\mu}_k)^2.$$\n", " \n", " 3. Ecrire une fonction *estimation(K,Z,X)* qui calcule l'EMV dans le cas où l'on connait les $(z_i,x_i)_{1\\le i\\le n}$. Cette fonction doit retourner les trois listes de taille K (alpha,mu,sigma). Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Vrai alpha et alpha estimé :\n", "[0.4, 0.3, 0.3] [0.386, 0.282, 0.332]\n", "\n", "Vrai mu et mu estimé :\n", "[-4, 4, 0] [-3.981 3.997 0.042]\n", "\n", "Vrai sigma et sigma estimé :\n", "[1, 1, 1] [0.996 1.056 0.98 ]\n" ] } ], "source": [ "(alpha_estime,mu_estime,sigma_estime)=estimation(K,Z,X)\n", "\n", "print(\"\\nVrai alpha et alpha estimé :\")\n", "print(alpha_real,alpha_estime)\n", "\n", "print(\"\\nVrai mu et mu estimé :\")\n", "print(mu_real,np.around(mu_estime,decimals=3))\n", "\n", "print(\"\\nVrai sigma et sigma estimé :\")\n", "print(sigma_real,np.around(sigma_estime,decimals=3))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Algorithmes\n", "\n", "Dans la suite de ce TP, on souhaite appliquer les algorithmes SEM, N-SEM et EM aux données $x_1, \\ld, x_n$ pour estimer $\\tta$. Estimer alors les $z_i$ : pour tout $i$, on choisit $$z_i:=\\argmax_k \\p_\\theta(Z=k | X=x_i).$$\n", "Comparer alors les valeurs de $\\tta$ ainsi que des $z_i$ à leurs vraies valeurs.\n", "\n", "### Algorithme SEM\n", "- initialisation arbitraire des classes $z_i\\in \\{1,\\ldots,K\\}$\n", "- répéter un grand nombre de fois:\n", " - a) calculer $\\widehat{\\theta}$\n", " - b) pour tout $i=1,\\ldots,n$, tirer la valeur de $z_i$ selon la loi $$\\mathbb{P}_{\\widehat\\theta}(Z_i=\\bullet | X_i=x_i)$$ \n", "\n", " 4. Ecrire une fonction `SEM(nbpas,K,x)` qui implémente l'algorithme précédent avec L étapes. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "SEM :\n", "\n", "Vrai alpha et alpha estimé :\n", "[0.4, 0.3, 0.3] [1. 1. 1.]\n", "\n", "Vrai mu et mu estimé :\n", "[-4, 4, 0] [ 3.885 -3.938 0.059]\n", "\n", "Vrai sigma et sigma estimé :\n", "[1, 1, 1] [1.153 1.014 0.833]\n" ] } ], "source": [ "L=int(100)\n", "print(\"\\nSEM :\")\n", "alpha_est,mu_est,sigma_est =0,0,0\n", "(alpha_est,mu_est,sigma_est)=SEM(L,K,X)\n", "\n", "print(\"\\nVrai alpha et alpha estimé :\")\n", "print(alpha_real,alpha_est)\n", "\n", "print(\"\\nVrai mu et mu estimé :\")\n", "print(mu_real,np.around(mu_est,decimals=3))\n", "\n", "print(\"\\nVrai sigma et sigma estimé :\")\n", "print(sigma_real,np.around(sigma_est,decimals=3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithme SEM à $N$ tirages\n", "- dupliquer $N$ fois le jeu d'observations $(x_1, \\ld, x_n)$, qui devient donc $(x_{i}^{(j)})_{1\\le i\\le n,\\, 1\\le j\\le N}$\n", "- appliquer SEM à ce jeu de données étendu \n", "\n", "5. Ecrire une fonction `NSEM(N,L,K,x)` qui implémente l'algorithme précédent. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "N-SEM :\n", "\n", "Vrai alpha et alpha estimé :\n", "[0.4, 0.3, 0.3] [1. 1. 1.]\n", "\n", "Vrai mu et mu estimé :\n", "[-4, 4, 0] [ 3.94 -3.947 0.066]\n", "\n", "Vrai sigma et sigma estimé :\n", "[1, 1, 1] [1.104 1.014 0.869]\n" ] } ], "source": [ "L=100\n", "N=10\n", "print(\"\\nN-SEM :\")\n", "(alpha_est,mu_est,sigma_est)=NSEM(N,L,K,X) \n", "\n", "print(\"\\nVrai alpha et alpha estimé :\")\n", "print(alpha_real,alpha_est)\n", "\n", "print(\"\\nVrai mu et mu estimé :\")\n", "print(mu_real,np.around(mu_est,decimals=3))\n", "\n", "print(\"\\nVrai sigma et sigma estimé :\")\n", "print(sigma_real,np.around(sigma_est,decimals=3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Algorithme EM\n", "- initialisation arbitraire du paramètre $\\theta_0$, \n", "- étant donné le paramètre $\\theta_t$, répéter pour $t=0,1,2,\\dots$ et jusqu'à convergence\n", " - Calcul de la matrice \n", " $$\\left[\\p_{\\tta_{t}}(Z=k | X=x_i)\\right]_{1\\le i\\le n, \\, 1\\le k\\le K}=\\left[\\f{\\al_k^tf_{\\mu_k^t,\\si^t_k}(x_i)}{\\sum_{l=1}^K\\al_l^tf_{\\mu_l^t,\\si_l^t}(x_i)}\\right]_{1\\le i\\le n, \\, 1\\le k\\le K}$$ \n", " - Calcul de $\\tta_{t+1}$ : pour tout $k=1,\\ld,K$, \n", " $$\\al_k^{t+1}=\\ff n\\sum_{i=1}^n \\p_{\\tta_{t}}(Z=k | X=x_i),$$\n", " $$\\mu_k^{t+1}=\\ff{n\\al_k^{t+1}}\\sum_{i=1}^n \\p_{\\tta_{t}}(Z=k | X=x_i)x_i,$$\n", " $$(\\si_k^{t+1})^2=\\ff{n\\al_k^{t+1}}\\sum_{i=1}^n \\p_{\\tta_{t}}(Z=k | X=x_i)(x_i-\\mu_k^{t+1})^2.$$\n", " \n", "6. Ecrire une fonction `EM(L,K,x)` qui implémente l'algorithme précédent. Appliquer cette fonction aux réalisations de la question précédente et comparer les valeurs estimées avec les valeurs réelles du mélange." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "EM :\n" ] } ], "source": [ "print(\"\\nEM :\")\n", "L=100 \n", "(alpha_est,mu_est,sigma_est)=EM(L,K,X) " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "Vrai alpha et alpha estimé :\n", "[0.4, 0.3, 0.3] [0.33462064 0.31353918 0.35184018]\n", "\n", "Vrai mu et mu estimé :\n", "[-4, 4, 0] [-4.09 2.332 0.687]\n", "\n", "Vrai sigma et sigma estimé :\n", "[1, 1, 1] [0.934 2.286 2.523]\n" ] } ], "source": [ "print(\"\\nVrai alpha et alpha estimé :\")\n", "print(alpha_real,alpha_est)\n", "\n", "print(\"\\nVrai mu et mu estimé :\")\n", "print(mu_real,np.around(mu_est,decimals=3))\n", "\n", "print(\"\\nVrai sigma et sigma estimé :\")\n", "print(sigma_real,np.around(sigma_est,decimals=3)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## EM en dimension supérieure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On se place cette fois en dimension deux. Chaque gaussienne du mélange est déterminée par sa moyenne (un vecteur de $\\R^2$ et sa matrice de covariance (une matrice $2\\times 2$).\n", "\n", "7. Ecrire une nouvelle fonction `échantillon2d` permettant d'échantillonner $n$ fois le mélange de gaussiennes de dimension 2 ayant les caractéristiques suivantes : \n", "- K=3\n", "- alpha=(.5,.3,.2)\n", "- les 3 gaussiennes ont les moyennes $\\mu_1=(-5,-5)$, $\\mu_2 = (5,5)$, $\\mu_3=(0,0)$. On pourra rassembler ces moyennes dans une liste $mu=[[-5,-5],[5,5],[0,0]]$\n", "- les 3 gaussiennes ont la même covariance $$\\begin{pmatrix} 1 & 0.5 \\\\ 0.5 & 1\\end{pmatrix}.$$ \n", "\n", "8. Appliquer la fonction précédentes pour réaliser une liste $X$ constituée de $n$ observations du mélange (donc $X$ est une liste de vecteurs de dimension 2).\n", "\n", "9. Ecrire une nouvelle fonction `EM2d` permettant d'estimer les paramètres du mélange à partir des observations $X$ et comparer le résultat obtenu avec les paramètres réels du mélange." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "attention les classes estimées peuvent être dans un ordre différent des classes réelles !\n", "\n", "Vrai alpha : [0.5, 0.3, 0.2]\n", "\n", "Estimation de alpha : [0.23196642 0.2688371 0.49919648]\n", "\n", "Vrai mu : [[-5, -5], [5, 5], [0, 0]]\n", "\n", "Estimation de mu: [[-0.02430919 0.18438086]\n", " [ 4.9245234 4.91653035]\n", " [-5.09633677 -4.9550138 ]]\n", "\n", "Vraie covariance : [array([[1. , 0.5],\n", " [0.5, 1. ]]), array([[1. , 0.5],\n", " [0.5, 1. ]]), array([[1. , 0.5],\n", " [0.5, 1. ]])]\n", "\n", "Estimation de la covariance : [[[1.2292515 0.64568701]\n", " [0.64568701 1.28573749]]\n", "\n", " [[0.90058606 0.33126681]\n", " [0.33126681 0.82164365]]\n", "\n", " [[0.95682149 0.50140257]\n", " [0.50140257 1.13658659]]]\n" ] } ], "source": [ "d=2\n", "K=3\n", "alpha=[.5,.3,.2]\n", "mu=[[-5,-5],[5,5],[0,0]]\n", "covariance=[np.array([[1,.5],[.5,1]])]*K\n", "n=int(5e2)\n", "L=100\n", "\n", "# on échantillonne les données\n", "(z0,x0)=echantillon2d(K,alpha,mu,covariance,n) \n", "# on applique l'algorithme EM et on vérifie que les paramètres estimés sont proches des vrais paramètres\n", "(alpha_em,mu_em,covariance_em)=EM2d(3,L,x0)\n", "\n", "print(\"attention les classes estimées peuvent être dans un ordre différent des classes réelles !\")\n", "print(\"\\nVrai alpha :\",alpha)\n", "print(\"\\nEstimation de alpha :\", alpha_em) \n", "\n", "print(\"\\nVrai mu :\",mu) \n", "print(\"\\nEstimation de mu:\",mu_em)\n", "\n", "print(\"\\nVraie covariance :\",covariance)\n", "print(\"\\nEstimation de la covariance :\", covariance_em)\n" ] } ], "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.10.4" } }, "nbformat": 4, "nbformat_minor": 2 }