expectation maximization

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fda7baad710>]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VfWd//HX52ZfyEIWtiSENci+hEXFKlItbmCtVdC6\n0tFOta2j047+Zmr768z8ZmwdbavW1tYFtK2i1Uq1FTesOxB2kC0gSxBISEIgCSHb9/fHvTBpCibA\nTc5d3s/HI4/cc+73ej9fvb7vyfd8z/macw4REYkOPq8LEBGR7qPQFxGJIgp9EZEootAXEYkiCn0R\nkSii0BcRiSIKfRGRKKLQFxGJIgp9EZEoEut1Ae1lZ2e7wsJCr8sQEQkry5cv3++cy+moXciFfmFh\nISUlJV6XISISVsxsR2faaXhHRCSKKPRFRKKIQl9EJIoo9EVEoohCX0Qkiij0RUSiiEJfRCSKhNw8\n/dOxp+YwizdWUFV3hMLsFM4flktyfER1UUTktERMIn68rZI5v/6Ytkv+ZibH8f1Lh3PF+DzvChMR\nCSERE/pj8zP4py8O5ZLRfeiXkcTqXQe4//VN3LlgNVvKa/nel4owM6/LFBHxVMSM6SfGxfDt6UMY\nlJNKYlwMkwdm8ft/mMKcSQU8+s5Wfv3eNq9LFBHxXMQc6R9PbIyP/7x8JAcbmvivv2xkZL90zhqU\n7XVZIiKeiZgj/RPx+YyfXDmaAVkp3LVgNTWHm7wuSUTEMxEf+gDJ8bH8dPZY9h1s4ME3NntdjoiI\nZ6Ii9AFG52VwzeQC5n+0nY17D3pdjoiIJ6Im9AHuuqCI1IRY7l+ko30RiU5RFfqZKfHMnTqQNzfs\nY93uGq/LERHpdp0KfTObYWabzKzUzO4+zvMJZvZc4PklZlYY2B9nZvPMbK2ZbTCze4Jb/sm7aWoh\naYmx/PRNHe2LSPTpMPTNLAZ4BLgIGA7MMbPh7ZrNBaqdc4OBB4H7Avu/CiQ450YBE4Bbj34heCUt\nMS5wtF9OafkhL0sREel2nTnSnwSUOue2OecagWeBWe3azALmBR6/AEw3/+WvDkgxs1ggCWgEPD+L\n+rUpBcTH+njyg+1elyIi0q06E/r9gF1ttssC+47bxjnXDNQAWfi/AOqAPcBO4H7nXNVp1nzaslIT\nuHxsX/6woowD9Y1elyMi0m26+kTuJKAF6AsMAO4ys4HtG5nZLWZWYmYlFRUVXVyS301nD6ChqZXn\nlu3quLGISIToTOjvBvLbbOcF9h23TWAoJx2oBK4BXnPONTnnyoEPgOL2b+Cce8w5V+ycK87JyTn5\nXpyCM/qkMaF/Js+V7MK1vTWniEgE60zoLwOGmNkAM4sHZgML27VZCNwQeHwl8LbzJ+lO4HwAM0sB\npgAbg1F4MFxdnM+2ijqW76j2uhQRkW7RYegHxuhvBxYBG4AFzrn1ZvYjM5sZaPY4kGVmpcCdwNFp\nnY8AqWa2Hv+Xx5POuTXB7sSpumR0H5LjY1hQoiEeEYkOnbrLpnPuz8Cf2+27t83jBvzTM9u/rvZ4\n+0NFSkIsl47uwytr9nDvZSNITYjom46KiETXFbnHc1VxPvWNLby2bq/XpYiIdLmoD/0J/TPJy0zi\nT6s/87oUEZEuF/Whb2ZcNqYv75fup6pOc/ZFJLJFfegDXDa6Ly2tjj+v3eN1KSIiXUqhD5zRpweD\nclI0xCMiEU+hj3+IZ+aYfizdXsXemgavyxER6TIK/YDLxvTBOXhVQzwiEsEU+gEDc1Ip6tWD19dr\n6qaIRC6FfhtfGtGLZdurqKw94nUpIiJdQqHfxoUjetPq4K2N5V6XIiLSJRT6bYzom0a/jCQN8YhI\nxFLot2FmXDiiF+9u2U/dkWavyxERCTqFfjsXDu9NY3Mr727unsVcRES6k0K/nYmFmWQmx7FIQzwi\nEoEU+u3ExviYfkYv3t5YTnNLq9fliIgElUL/OKYPy+VgQzMrdh7wuhQRkaBS6B/H1CHZxMUYb23c\n53UpIiJBpdA/jh6JcUws7MlizdcXkQij0D+B84flsnlfLWXV9V6XIiISNAr9E5g2LBdAR/siElEU\n+icwMDuF/lnJvK3QF5EIotA/ATNjWlEuH26t5HBji9fliIgEhUL/c5w/LJcjza18tG2/16WIiASF\nQv9zTB7Yk+T4GA3xiEjEUOh/joTYGM4enM3ijRU457wuR0TktCn0O3D+sFx2HzjM5n21XpciInLa\nFPodmFYUmLq5SUM8IhL+FPod6J2eyBl90jRfX0QigkK/E6YV5VCyo5qDDU1elyIicloU+p0wbVgu\nLa2O97do6qaIhDeFfieMy88gLTFWQzwiEvYU+p0QG+PjnKE5vLO5gtZWTd0UkfCl0O+kaUW5VBw6\nwid7DnpdiojIKVPod9K5Q3MAeEdTN0UkjCn0OymnRwKj89JZvKnC61JERE6ZQv8knFeUy8qd1VTX\nNXpdiojIKVHon4RpRTm0Onh3i472RSQ8dSr0zWyGmW0ys1Izu/s4zyeY2XOB55eYWWGb50ab2Udm\ntt7M1ppZYvDK716j8zLomRLPOxriEZEw1WHom1kM8AhwETAcmGNmw9s1mwtUO+cGAw8C9wVeGws8\nA3zDOTcCOA8I28taY3zGuUNz+OvmClo0dVNEwlBnjvQnAaXOuW3OuUbgWWBWuzazgHmBxy8A083M\ngAuBNc651QDOuUrnXFgvQ3VeUQ5VdY2sKTvgdSkiIietM6HfD9jVZrsssO+4bZxzzUANkAUMBZyZ\nLTKzFWb2vdMv2VtfGJKDz9AsHhEJS119IjcWmApcG/j9ZTOb3r6Rmd1iZiVmVlJREdphmpkSz7iC\nTM3XF5Gw1JnQ3w3kt9nOC+w7bpvAOH46UIn/r4J3nXP7nXP1wJ+B8e3fwDn3mHOu2DlXnJOTc/K9\n6GbTinJYU1ZDxaEjXpciInJSOhP6y4AhZjbAzOKB2cDCdm0WAjcEHl8JvO386wsuAkaZWXLgy+Bc\n4JPglO6d8wILq/x1c2j/VSIi0l6HoR8Yo78df4BvABY459ab2Y/MbGag2eNAlpmVAncCdwdeWw08\ngP+LYxWwwjn3avC70b1G9E0jt0eCVtMSkbAT25lGzrk/4x+aabvv3jaPG4CvnuC1z+CfthkxzIzz\ninL4y7q9NLe0Ehuja9xEJDworU7RtKJcDjU0s2Knpm6KSPhQ6J+is4dkE+szDfGISFhR6J+itMQ4\nigsztZqWiIQVhf5pmFaUy8a9h9hb0+B1KSIinaLQPw3ThvmnbupCLREJFwr90zAkN5V+GUka1xeR\nsKHQPw1Hp26+v2U/jc2tXpcjItIhhf5pmlaUS11jCyXbq7wuRUSkQwr903TW4CziY3wa4hGRsKDQ\nP03J8bFMHthTt1oWkbCg0A+CaUW5lJbXsquq3utSREQ+l0I/CI5O3Xxrwz6PKxER+XwK/SAYkJ3C\nkNxUFq1X6ItIaFPoB8mMkb1Z8mklVXWNXpciInJCCv0g+dKI3rQ6ePMTHe2LSOhS6AfJiL5p5GUm\n8dr6vV6XIiJyQgr9IDEzvjSiN+9v2c+hhiavyxEROS6FfhDNGNmbxpZWzdkXkZCl0A+i8QWZZKcm\nsGidhnhEJDQp9IMoxmdcOKIXizeV09DU4nU5IiJ/R6EfZBeN7E19Y4tW1BKRkKTQD7IzB2aRnZrA\ny6s+87oUEZG/o9APstgYH5eO7sPbm8o5qFk8IhJiFPpdYObYvjQ2t+qEroiEHIV+FxiXn0F+zyQW\nrtYQj4iEFoV+FzAzZo7pywel+6k4dMTrckREjlHod5GZY/rR6uDVNTraF5HQodDvIkW9e3BGnzRe\nWFHmdSkiIsco9LvQ7In5rNt9kHW7a7wuRUQEUOh3qcvH9iM+1seCkl1elyIiAij0u1R6chwzRvTm\njyt367YMIhISFPpdbPbEfA42NPOa5uyLSAhQ6HexKQOzyO+ZxHPLNMQjIt5T6Hcxn8+YPbGAj7ZV\nsmXfIa/LEZEop9DvBrMn5hMf6+PJD7d7XYqIRDmFfjfISk3gy2P78eKKMqrrGr0uR0SimEK/m9w0\ntZCGplZ+v2yn16WISBRT6HeTYb3TOHtwFvM/3EFTS6vX5YhIlOpU6JvZDDPbZGalZnb3cZ5PMLPn\nAs8vMbPCds8XmFmtmf1zcMoOT1+fOpC9Bxt4UbdmEBGPdBj6ZhYDPAJcBAwH5pjZ8HbN5gLVzrnB\nwIPAfe2efwD4y+mXG97OK8phdF46D71dqqN9EfFEZ470JwGlzrltzrlG4FlgVrs2s4B5gccvANPN\nzADM7HLgU2B9cEoOX2bGHV8cQln14Yg+2q890szWilrWltWw9NMqVu6sZtPeQ+w72EBrq/O6PJGo\nFtuJNv2AtlcWlQGTT9TGOddsZjVAlpk1AP8CXACccGjHzG4BbgEoKCjodPHhaFpR7rGj/SvG5xEX\nE96nVXZU1lGyvZqVu6pZW1bDjqp6DtSfeJnI+BgfeZlJDO+bxpi8DCYN6Mmofun4fNaNVYtEr86E\n/un4IfCgc642cOB/XM65x4DHAIqLiyP6UPDo0f7NT5XwuyU7ueGsQq9LOinOOVbtOsBf1u3lzQ37\n2FZRB0BqQiyj+qVz8ag+5Gcm0yc9keT4GJLiY2hucdQ3tlBV30hZdT079tezcucBXlmzB4Ds1ATO\nH5bDzDH9OGtQlr4ARLpQZ0J/N5DfZjsvsO94bcrMLBZIByrx/0VwpZn9GMgAWs2swTn38GlXHsam\nFeVy5sAsHnxzM7PG9iUjOd7rkjpUcegIL64o4/nlZZSW1xIXY0wZmMV1U/pz9uBsBuWkEnOSYV1+\nqIEPSvfz1oZy/rJ2LwtKysjvmcTVxflcO7k/mSmh/+9FJNyYc59/YB0I8c3AdPzhvgy4xjm3vk2b\n24BRzrlvmNls4Arn3FXt/jk/BGqdc/d/3vsVFxe7kpKSU+lLWNmw5yCX/Pw9rj+zkB/OHOF1OSdU\nWn6IX7/7KS+t3E1jSysT+mfy1Ql5XDy6D2mJcUF7n4amFl7/ZB/PLt3Jh1srSY6P4drJBXz9nIH0\nSksM2vuIRCozW+6cK+6oXYdH+oEx+tuBRUAM8IRzbr2Z/Qgocc4tBB4HnjazUqAKmH165Ue+M/qk\nMWdSAU9/vIMrJ+Qxsl+61yX9jU8+O8gDb2zmzQ37SIj1cdXEPG48awCDc1O75P0S42KYOaYvM8f0\nZcu+Q/zina08/v6nPP3xDm45ZyC3njuIlISuHo0UiXwdHul3t2g50gc4UN/IhQ++S8+UeBbePpX4\nWO9P6pZV1/PA65t5adVu0hLjuOnsQq6b0p+s1IRur2VHZR33v76ZP63+jNweCdxz8TAuH9uPzzs/\nJBKtOnukr9D32Fsb9jF3Xgm3TxvMP3+pyLM6qusaeWRxKfM/2gEGN51dyDfPHUx6cvCGcE7V8h3V\n/Psrn7Bq1wHOHZrD/7tiFP0ykrwuSySkKPTDyF0LVvPSyjLm3TyJc4bkdOt7NzS18MQHn/LoO1up\nO9LMV8bn8U8XDKVviIVqa6tj/kfb+fGiTRjwfy45g2smFeioXyRAoR9G6hubueIXH7L3YAN/un0q\n+T2Tu/w9W1odLyzfxYNvbGHvwQamD8vlezOGUdS7R5e/9+nYVVXPPS+u5f3S/Vw0sjf//ZXRpCd5\n/9eIiNcU+mFmR2Udlz30Pr3SEnnu1jPp2UXTFZ1zvLWhnPte28iW8lrG5mdwz0XDmDwwq0veryu0\ntjp+8/42fvzaJnqlJfLQNeMYX5DpdVkinups6Ht/5lAA6J+Vwq+uK2ZnVT3XP7GEA/XBv+/+x9sq\n+cqjH/L1+SU0tzoevXY8L33zrLAKfPCvRnbLFwbx/DfOxAyu+uVHPPPxDq/LEgkLOtIPMYs3lnPr\n08vJ65nEUzdOoiDr9Id61pbVcP/rm/jr5gp6pSXwnelD+Wpx+N8CAqDmcBN3PLuSxZsq+NqUAn5w\n2YiI6JfIydLwThhb+mkVtzxdQkur499njWTW2L4nfcLSOce7W/bz2Ltb+aC0kozkOL553iCuP7OQ\nxLiYLqrcGy2tjh8v2siv/rqNyQN68ujXJnTZ8JhIqFLoh7mdlfX804JVLN9RzeQBPbnji0OZMrBn\nh+G/q6qel1ft5sUVu9m2v45eaQncfPYA5kwuCOoVtKHopZVl/Msf1pKXkcS8myd1ywlxkVCh0I8A\nLa2O3y7Zwc/fKmV/7REKs5I5d2gOo/Iy6JWWQHyMj/qmFsqqD7Np70E+3Fp57AZokwb05KrifGaO\n6RsSF311l5LtVcydV+JfiP7GiSF3pbNIV1HoR5D6xmb+vHYvL6/aTcn2ag43tfxdm5T4GCYO6MnU\nwdlcOLx3UM4FhKst+w5xwxNLOdjQzC+/NoGpQ7K9Lkmkyyn0I1RTSytl1YcpP9hAc6sjMc5Hv4xk\neqUl6EKlNvbWNHDjk0vZWlHLQ3PGMWNkH69LEulSCn2JejWHm7jpyaWsLqvhgavGMGtsP69LEuky\nmqcvUS89KY75cydT3D+TO55bxfMluzp+kUiEU+hLREtNiOWpmyYxdXA2331hjS7ikqin0JeIlxQf\nw6+vL+aLZ+Tyb39cx/yPtntdkohnFPoSFRLjYvjFtRO4YHgv7n15Pc8u3el1SSKeUOhL1IiP9fHw\nNeM4ryiHe15ay4sryrwuSaTbKfQlqiTExvDLr03grEFZ/PPzq3llzWdelyTSrRT6EnUS4/xj/MX9\ne/KdZ1exaP1er0sS6TYKfYlKyfGxPHHTREbnpXP771aweFO51yWJdAuFvkSto9M5h/bqwT8+s5yS\n7VVelyTS5RT6EtXSk+KYd/Mk+qQncfNTy9iw56DXJYl0KYW+RL3s1ASenjuJ5PhYrn9iKTsr670u\nSaTLKPRFgLzMZJ6eO4mmlla+9vgSyg81eF2SSJdQ6IsEDOnVgydvnMj+2iNc//hSag43eV2SSNAp\n9EXaGFeQya+um8DWilrmPrWMw41/v3aBSDhT6Iu0c86QHH569TiW76zmm79dTlNLq9cliQSNQl/k\nOC4Z3Yf/uHwkizdV8N3nV9PaGlrrToicqlivCxAJVddO7k91XSP3v76ZzJR47r10uFYnk7Cn0Bf5\nHLdNG0xlXSNPfrCd7NQEbps22OuSRE6LQl/kc5gZ379kOAfqm/jJok1kJsdzzeQCr8sSOWUKfZEO\n+HzGj68czYH6Rv7tj2vJTI7jolFaaF3Ck07kinRCXIyPX1w7gXEFmXzn2VV8WLrf65JETolCX6ST\nkuJjePyGYgqzk/mH+SWsKTvgdUkiJ02hL3ISMpLjmX/zZDKS47nxyWVsq6j1uiSRk6LQFzlJvdMT\neebrkzHguseXsrdG9+mR8KHQFzkFA7JTmHfzJGoON3Hd40s4UN/odUkindKp0DezGWa2ycxKzezu\n4zyfYGbPBZ5fYmaFgf0XmNlyM1sb+H1+cMsX8c7Ifuk8dv0EdlTWc9NTy6hvbPa6JJEOdRj6ZhYD\nPAJcBAwH5pjZ8HbN5gLVzrnBwIPAfYH9+4HLnHOjgBuAp4NVuEgoOGtQNj+fM5bVuw7wj8+soLFZ\n9+mR0NaZI/1JQKlzbptzrhF4FpjVrs0sYF7g8QvAdDMz59xK59xngf3rgSQzSwhG4SKhYsbIPvzn\nl0fx180VfPcF3adHQltnLs7qB+xqs10GTD5RG+dcs5nVAFn4j/SP+gqwwjl35NTLFQlNcyYVUFXX\neOyq3R9cpvv0SGjqlityzWwE/iGfC0/w/C3ALQAFBbrEXcLTN88bRFVdI4+//ylZKfF8a/oQr0sS\n+TudGd7ZDeS32c4L7DtuGzOLBdKBysB2HvAScL1zbuvx3sA595hzrtg5V5yTk3NyPRAJEWbGv158\nBleM68f/vLGZZz7e4XVJIn+nM6G/DBhiZgPMLB6YDSxs12Yh/hO1AFcCbzvnnJllAK8CdzvnPghW\n0SKhyucz7rtyNOcPy+X7L6/jD8vLvC5J5G90GPrOuWbgdmARsAFY4Jxbb2Y/MrOZgWaPA1lmVgrc\nCRyd1nk7MBi418xWBX5yg94LkRDiv0/PeM4elM13X1jNy6va/2Es4h1zLrRmGhQXF7uSkhKvyxA5\nbYcbW7jxyaWU7Kjm4TnjdGdO6VJmttw5V9xRO12RK9JFkuJjeOLGiYzNz+Bbv1/JG5/s87okEYW+\nSFdKSYjlyZsmMqJfOt/87XIWbyr3uiSJcgp9kS6WlhjH/JsmUdS7B7c+vZz3tlR4XZJEMYW+SDdI\nT47j6ZsnMzA7hbnzSli8UUf84g2Fvkg3yUyJ5/f/MIWhvVK55ekSXlu31+uSJAop9EW6UWZKPL/9\n+hRG9Uvntt+t0HRO6XYKfZFulp4Ux/y5k5nQP5M7nlvFgpJdHb9IJEgU+iIeSE2IZd5Nk5g6OJvv\nvbCGeR9u97okiRIKfRGPJMXH8Ovri7lgeC9+sHA99y/aRKhdLCmRR6Ev4qHEuBgevXY8syfm8/Di\nUv7lD2tobtFCLNJ1uuXWyiJyYrExPv7rilHk9kjg52+XUlnbyMPXjCcpPsbr0iQC6UhfJASYGXde\nWMS/Xz6StzeVc+1vPqayVusNSfAp9EVCyHVT+vPoteNZ/9lBZj3yAZv2HvK6JIkwCn2REDNjZB8W\n3Homjc2tXPGLD3h7o27UJsGj0BcJQWPyM3j59rMZkOO/bcNv3tummT0SFAp9kRDVJz2JBbeeyYwR\nvfmPVzdw54LV1Dc2e12WhDmFvkgIS46P5ZFrxnPnBUP546rdXP7IB2ytqPW6LAljCn2REOfzGd+e\nPoT5N09if20jMx96n1fX7PG6LAlTCn2RMHHOkBxe+dZUhvbuwW2/W8G9L6+joanF67IkzCj0RcJI\n34wknrvlTOZOHcD8j3Zw6UPvs253jddlSRhR6IuEmfhYH9+/dDhPz53EoYYmvvyLD3hkcSktrZrd\nIx1T6IuEqXOG5LDoji9w4fDe/GTRJq7+1UeUluskr3w+hb5IGMtIjufha8bxwFVj2FJey8U/e4+f\nvbmFxmbdtE2OT6EvEubMjCvG5/HmnefypZG9efDNzVzy8/dYtr3K69IkBCn0RSJETo8EHpozjidv\nnEh9Ywtf/eVHfPv3K/nswGGvS5MQotAXiTDThuXyxp1f4FvnD2bR+r2c/z/v8OAbmzncqOmdotAX\niUjJ8bHcdWERb911LtPP6MXP3trCuT9ZzPyPtnOkWeEfzRT6IhEsLzOZR64Zz/PfOJPCrBTufXk9\n037yDr9bspMmrdAVlSzU7txXXFzsSkpKvC5DJOI453i/dD//8/pmVu06QF5mEl+fOoCvFueTkqBF\n9MKdmS13zhV32E6hLxJdnHO8s6mCh97ewoqdB0hPiuPayQXccFYhvdISvS5PTpFCX0Q6tHxHNb95\nbxuL1u8lxmdcOKI3cyYWcNagLHw+87o8OQmdDX39TScSxSb0z2RC/wnsqKxj3oc7eHFlGa+u2UN+\nzySuLs7nivF59M1I8rpMCSId6YvIMQ1NLbz+yT6eXbqTD7dWAv4vhktG9eHiUX3ona7hn1Cl4R0R\nOS07K+v505rPeGXNHjbsOQhAcf9Mzj8jl2lFuQzr3QMzDQGFCoW+iATN1opaXl2zh7+s23vsC6B3\nWiLnDs3hnKHZTCrsSa5OAntKoS8iXWLfwQb+uqmCdzaX897m/Rw64l+3t39WMhMLezKpsCfj+2cw\nIDuVGJ0M7jYKfRHpck0traz/7CDLPq1i6fYqSrZXUV3fBEBSXAzD+6Yxsm8aI/qlc0bvNAbmpOia\ngC4S1NA3sxnAz4AY4DfOuf9u93wCMB+YAFQCVzvntgeeuweYC7QA33bOLfq891Loi4Sv1lbH1opa\nVpfVsP6zGtbvPsj6z2qoa3Pfn95piQzKTWFgdioDc1Lon5VM34wk+mYkkZYY52H14S1oUzbNLAZ4\nBLgAKAOWmdlC59wnbZrNBaqdc4PNbDZwH3C1mQ0HZgMjgL7Am2Y21Dmnm3+IRCCfzxjSqwdDevXg\nygl5gP+LYHtlHZv3HWJrRR1bK2rZWlHHH1ft5lBD89+8vkdCbOALIJHe6Ulkp8aTlRJPz9QEslPi\nyUpNoGdKPJnJccTG6C4yp6Izf2dNAkqdc9sAzOxZYBbQNvRnAT8MPH4BeNj8p/VnAc86544An5pZ\naeCf91FwyheRUOfzGQNzUhmYk/o3+51zVNQeoaz6MJ8dOPrTwO7A4zVlNVTXN3KiVSCT42PokRhL\nakIsPRLj6JEY6/9JiCM1MZakuBgSYn0kxsWQGOcjIS7G//jYPv/+uBgfsT4j9thvI8ZnxPp8xMYY\nsT7/dpzPFxEXrHUm9PsBu9pslwGTT9TGOddsZjVAVmD/x+1e2++UqxWRiGFm5PZIJLdHIuMLMo/b\npqXVUXO4icraI+yvbaSqrpHKuiNU1TVS29DMoYZmDh1p8v9uaOazA4c51NBM7ZFmGppaTviFceo1\nQ5zPR4zP8Jm/DxbYb/a/+/zfDUe3wWjTPrDPF3itz98AnxnTinL410uGB7fodkLijIqZ3QLcAlBQ\nUOBxNSISKmJ8Rs+UeHqmxDOk18m91jlHU4vjSHMLDU2tNDS1HHvcdl9TSyvNrY6WVn/7llb/dnOL\nC/z+3+ePPm5udTjnaHXgHLQGzo22Onds2wVqOLbt8LfH//jY6wOvw0Hv9K6/+rkzob8byG+znRfY\nd7w2ZWYWC6TjP6HbmdfinHsMeAz8J3I7W7yIyImYGfGxRnysjx66hOCYzpwJWQYMMbMBZhaP/8Ts\nwnZtFgI3BB5fCbzt/NOCFgKzzSzBzAYAQ4ClwSldREROVodH+oEx+tuBRfinbD7hnFtvZj8CSpxz\nC4HHgacDJ2qr8H8xEGi3AP9J32bgNs3cERHxji7OEhGJAJ2dp6+JriIiUUShLyISRRT6IiJRRKEv\nIhJFFPoiIlEk5GbvmFkFsOM0/hHZwP4gleOlSOkHqC+hKFL6AerLUf2dczkdNQq50D9dZlbSmWlL\noS5S+gHqSyiKlH6A+nKyNLwjIhJFFPoiIlEkEkP/Ma8LCJJI6QeoL6EoUvoB6stJibgxfRERObFI\nPNIXEZGwN4zqAAADdUlEQVQTiJjQN7MZZrbJzErN7G6v6zkZZvaEmZWb2bo2+3qa2RtmtiXw+/hL\nC4UQM8s3s8Vm9omZrTez7wT2h2NfEs1sqZmtDvTl/wb2DzCzJYHP2XOB242HBTOLMbOVZvZKYDss\n+2Jm281srZmtMrOSwL5w/IxlmNkLZrbRzDaY2Znd0Y+ICP02i7dfBAwH5gQWZQ8XTwEz2u27G3jL\nOTcEeCuwHeqagbucc8OBKcBtgf8O4diXI8D5zrkxwFhghplNAe4DHnTODQaqgbke1niyvgNsaLMd\nzn2Z5pwb22Z6Yzh+xn4GvOacGwaMwf/fpuv74V/OK7x/gDOBRW227wHu8bquk+xDIbCuzfYmoE/g\ncR9gk9c1nkKfXgYuCPe+AMnACvxrQ+8HYgP7/+ZzF8o/+Fetews4H3gFsDDuy3Ygu92+sPqM4V9d\n8FMC51W7sx8RcaTP8RdvD/cF2Hs55/YEHu8FTnKFUG+ZWSEwDlhCmPYlMByyCigH3gC2Agecc82B\nJuH0Ofsp8D2gNbCdRfj2xQGvm9nywPraEH6fsQFABfBkYMjtN2aWQjf0I1JCP6I5/9d+2EyzMrNU\n4A/AHc65g22fC6e+OOdanHNj8R8lTwKGeVzSKTGzS4Fy59xyr2sJkqnOufH4h3NvM7MvtH0yTD5j\nscB44FHn3DigjnZDOV3Vj0gJ/U4twB5m9plZH4DA73KP6+kUM4vDH/i/dc69GNgdln05yjl3AFiM\nfwgkw8yOLjMaLp+zs4GZZrYdeBb/EM/PCM++4JzbHfhdDryE/ws53D5jZUCZc25JYPsF/F8CXd6P\nSAn9zizeHm7aLjZ/A/7x8ZBmZoZ/veQNzrkH2jwVjn3JMbOMwOMk/OcmNuAP/ysDzcKiL865e5xz\nec65Qvz/b7ztnLuWMOyLmaWYWY+jj4ELgXWE2WfMObcX2GVmRYFd0/GvJd71/fD6hEYQT4xcDGzG\nP+76r17Xc5K1/x7YAzThPwKYi3/M9S1gC/Am0NPrOjvRj6n4/xxdA6wK/Fwcpn0ZDawM9GUdcG9g\n/0BgKVAKPA8keF3rSfbrPOCVcO1LoObVgZ/1R/9fD9PP2FigJPAZ+yOQ2R390BW5IiJRJFKGd0RE\npBMU+iIiUUShLyISRRT6IiJRRKEvIhJFFPoiIlFEoS8iEkUU+iIiUeT/A7wZzsmGIzqTAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fda7bb15b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.mlab as mlab\n",
    "%matplotlib inline\n",
    "\n",
    "ls = np.linspace(0, 60, 1000)\n",
    "\n",
    "multimodal_norm = mlab.normpdf(ls, 0, 5) + mlab.normpdf(ls, 20, 10)\n",
    "\n",
    "plt.plot(ls, multimodal_norm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import mixture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 2)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# concatenate ls and multimodal to form an array of samples\n",
    "# the shape is [n_samples, n_features]\n",
    "# we reshape them to create an additional axis and concatenate along it\n",
    "samples = np.concatenate([ls.reshape((-1, 1)), multimodal_norm.reshape((-1,1))], axis=-1)\n",
    "samples.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "gmix = mixture.GaussianMixture(n_components = 2, covariance_type = \"full\")\n",
    "fitted = gmix.fit(samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitted.means_.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 19.69808581,   0.03686826])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitted.means_[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}