{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "cmap = plt.get_cmap('tab10')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_pdf(x, mu, var):\n",
    "    D = len(x)\n",
    "    e = 1e-8\n",
    "    ln_pdf = - 0.5 * ( np.dot(x - mu, x - mu) / (var + e) + D * np.log(np.abs(var)) + D * np.log(2*np.pi) )\n",
    "    return np.exp(ln_pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "D = 2\n",
    "a = 10\n",
    "b = 100\n",
    "w = 0.5\n",
    "#theta = np.random.multivariate_normal(mean=np.zeros(D), cov=np.identity(D))\n",
    "theta = np.array([1,1])\n",
    "\n",
    "X = np.zeros((N,D))\n",
    "X[:int(w*N)] = np.random.multivariate_normal(mean=theta, cov=np.identity(D), size=int(w*N))\n",
    "X[int(w*N):] = np.random.multivariate_normal(mean=np.zeros(D), cov=np.identity(D)*a, size=N-int(w*N))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGAhJREFUeJzt3X+MnVWdx/HPl3EqXdSwq6VkaWlLqYzd+gN3QEllwxaqbGFpovsHGMXQxAajwCbWH0g0G5JNzGpUzG6yKVizBpQs4G43Mv4oAV00ttsplBXa4acWRqWOEpUsgwztd/+4d+rtzJ25d+5z7nPOeZ73KyHtnd6599xh5jPnfs/3nMfcXQCA8p0QewAAUFcEMABEQgADQCQEMABEQgADQCQEMABEEiSAzexkM7vTzMbM7KCZnRficQGgyl4R6HFukvQdd/87M1sk6U8CPS4AVJYV3YhhZq+R9JCkM5xdHQDQtRAz4DMkTUj6qpm9WdI+Sde5+/+13snMtkraKkknnXTSXw4NDQV4agBIz759+37t7ks63S/EDHhY0m5J6919j5ndJOn37v7puT5neHjYR0dHCz0vAKTKzPa5+3Cn+4VYhBuXNO7ue5q375T01gCPCwCVVjiA3f1ZSc+Y2VnND10o6UDRxwWAqgvVBXGNpNuaHRBPSboq0OMCQGUFCWB33y+pY70DAPBH7IQDgEgIYACIhAAGgEgIYACIhAAGgEgIYACIhAAGgEgIYAAYG5Hu3tb4s0QEMIB6GxuR7toi7b258WeJIUwAA6i3J++VpiYbf5+abNwuCQEMoN5Wb5AGFzf+Pri4cbskoQ7jAYA8DW2S3rOjMfNdvaFxuyQEMNI1NhLlhwI1NLQpyvcYJQikKeLCCFAWAhhpirgwApSFAEaaIi6MAGWhBow0RVwYAcpCACNdkRZGgLJQggCASAhgAIiEAAaASAhgAIiEAM5ZpCP0AIRBAOeKnWJA9gjgXLFTDMgeAZwrdooB2WMjRq7YKQZkjwDOGTvFgKxRggByQddL5RDAQA7oeqkkAhjIAV0vlUQAAzmg66WSWIQDckDXSyURwEAu6HqpHEoQABAJAQwAkRDAABAJAQwAkRDAABAJAQwAkRDAABAJAQwAkRDAABAJAQwAkQQLYDMbMLMHzexboR4TAKos5Az4OkkHAz4eAFRakAA2s2WSLpF0S4jHA4A6CDUD/pKkj0s6OtcdzGyrmY2a2ejExESgpwWAfBUOYDO7VNKv3H3ffPdz9+3uPuzuw0uWLCn6tACQvRAz4PWSLjOzn0m6XdIGM7s1wOMCQKUVDmB3v97dl7n7SkmXS7rX3d9XeGQAUHH0AQNAJEEvSeTu35f0/ZCPCQBVxQwY/TM2It29rfEngFkIYPTH2Ih01xZp782NPwlhYBYCGP3x5L3S1GTj71OTjdsAjkMAoz9Wb5AGFzf+Pri4cRvAcYIuwgHHDG2S3rOjMfNdvaFxG/03NsLXPCMEMPpnaBMhUKbpuvvUpLT/1sYvQL7+SaMEARSVSrcHdffsEMBAESl1e1B3zw4lCKCIdrPOWG/7qbtnhwAGili9oVFvnZpMY9ZJ3T0rBDBQBLNOFJB/ANN2g9iYdaJHeQcwbTfAvHYdOKz7H5/Q+WuWaOPapbGHgxny7oKg7QaY064Dh3XtNx7U1358SNd+40HtOnA49pAwQ94BTNsNMKf7H5/Q5NQRSdLk1BHd/zjXYkxN3gE8vQByzgcpPwAznL9miRYPDkiSFg8O6Pw1XIsxNebupT/p8PCwj46Olv68aINFzEqjBhyHme1z9+FO98t7EQ7FsIiZjV6DdOPapQRvwvIuQaAYFjGzwGJadRHAdcYiZhZYTKsuArjOKraIuevAYX1m58OVmyGymFZdLMKhEqbfpk9OHdHiwQF9+YqzK1X7ZDEtLyzCoVbavU2vUlDNtZhGMOeNEgSiCVkyqOPbdBbn8scMGFG0lgzuGB0vXDLYuHapvnzF2bWaDVZ91t+TzPramQEjin6s7G9cu1Q3bl5XmxCq46x/XildnaRLBDCiIDyKm571X3neisotOvYkw752ShCIoo4lg150WmTreadbZm/Vu5La1Um6QBsakKi+tda1bkEfXFyJHvBjEvnF0m0bGiUIIFF92wGX4Vv1rg1tki75fDa/UAhgIFF9q5OzBT0ZlCCQlV42HuS8WaFvY0/krXpVdVuCIICRjV5qolXfoow0UQOum7ER6e5tWfQ+9qqXmigniSFlBHAVZNiA3oteaqL0GyNl9AFXQbtV7QrW9XrpHabfGCkjgKsgwwb0XvWy8YDL8iBVBHAVTB+szqo2kBUCuCqGNhG8QGZYhAOASAhgAIiEAAaASAhgAIiERThUUs7nP6A+Cs+AzWy5md1nZgfN7BEzuy7EwIBecbFK5CJECeJlSR919zdIerukD5vZ2gCPC/SE8x8qpsLnnBQOYHf/pbs/0Pz785IOSjqt6OMCveL8hwqp+DknQWvAZrZS0tmS9rT5t62StkrS6aefHvJpgeMkc/4DZ+4WV/FzToKdB2xmr5L0A0n/6O7fnO++nAeM3Cx4Ua9f112rW6hnev26bs8DDjIDNrNBSXdJuq1T+AK5aT3U/Y7R8e4Ode/HzK01jPbfmk0YFVLxc05CdEGYpK9IOujuXyg+JCAtPS3q9eO6a1W+mOZ8MrvQ5kKE6IJYL+n9kjaY2f7mf9X7SqH65lht72lRb3rmds4Hw81UuZhm5XBNOEDqWGsMubGj0GPVrQacqVJrwEDf9Tt4OtRsQx3q3lpPvm33IV19wZn62LvO6v4BOh07SkBnpR5nQVS4kbsWyugFLentfWs9+YhL//qDJ8Pt1Kt4z2wVVT+A+abMXxmLT/2o2bZx/polGrA/3j5y1PX1PYf0mZ0PFw/iui7SZaz6Acw3Zf7KWnwqYbV949qluvqCMzVwQiOFFw2coB898Zsw51awSJed6teAa3TByio5fqGqWr2gH3vXWXrL8pN1/+MTeua5F3Tfo422tukWt55rzRXvma2ienRBsDCRldaFqsWDA91tfAj43PN1KIQ+5rKfr5UjOePptguiHgGMrHxm58P62o8PHbt95XkrdOPmdX1/3k5h2K+w7EdQxvwlhu4DuPo1YGQn1mlmnXa89euYy41rl+rGzeuCBiRHcuaBAEZypk8zu/K8FaXO3DoFf07HXL76xMF5bxe168DhMJ0bNVf9RThkKdTGh4U+53zHWCZzzGUXnn9xat7bRfR0OBHaIoCBFp2Cf+a/p7rQdf6aJbpjdPxYDTjkbL1deSOl154TArgTOihqoZcgLToTLBre831+P2fr/Qz3uqELYj6ZHgaNhem1Y6BIt0bRLoXYXQ6pzvxTQRdECOyiq4VeOwZmLsq9+sTBrheminYpxO5y6EfnRh0RwPNha2ct9Nrd0NqtseUdq7Tjhz/tektx0Y6KnDoyMDdKEJ1QA85Hgf9XRd9S91KO6GcNGHGxEw71ErleH7smi7RwIDvqJfLly3PqEUY6CGCkaaHlhAROvYuxeaQdShP5oASB9LSUE44MnKjbV/yDThl+d+cw6UO9PrcwoxSSBtrQkK+WcsLAkRc19di93R1WHvhA9ekwC3JYeklit6dhYQhgpKel/e8FX6QfHn1jlDDJMcxoT8sLNWCkp3llh6f3fkuffezPdc/Rs6OESY5bblkMzAs1YCQtdg029vMjT/QBA8hLhTY9sQgHIB/TnS97b278OTYSe0SlIICBHnFViIBqevAVAQz0IMcWtaTV9OArAhjoQY4taklrdr7onA/W6txt2tCAHuTYopa8oU21Cd5pBDBmya71qsTV89avzVz9tp2+fl1/fSO9riz+n1cEbWg4TnZnCZR4DGU3X5tO9+n665vY68LC0IaGnmRX2yxx9bybr02n+3T99U3sdaE/CGAcJ7uzBEpcPe/ma9PpPl1/fRN7XegPShCYJbt6YGK1UmrAYCsyUGcV2tabIy5JBHQhx5lfxzG3LuDtv7VWfbW5oQaM2spxN1tXY67ptt4cEcCorW5X/2ee+RDzDIiuxlzTbb05ogSB2upmN1trj+wdo+Pa8o5V2vHDnx67XXbPbFc78Ka39VIDTh4BjNrq5uoRM2ec9xx4dtYMNEgAz1g0m6vO2/UVL2q4rTdHBDBqZ2a4zRegM2ecF609VU8/99OwZ0DMWDTbf+4XdO1//9mcs+xOY0Y+CGDUysySQqcSQrsZ51uWnzxrBlqom2LGotmLj96jyan3SAo8y0ZyCGDUSrtFrE7hNnPGOfP2QkN9ltUbGu1izXMfTjzrIi0+PMBJazUQJIDN7GJJN0kakHSLu382xOMiIBrzJYU7RrJ1xttLqB9nxqLZW4Y26cun5defjIUrvBPOzAYkPSZpo6RxSXslXeHuB+b6HHbClazEk7VyUHTzxczTw1o7IzhNDFK5O+HOlfSEuz/VfOLbJW2WNGcAo2TtGvPLCuAEZ95zLWJ1G8wzZ7zPvzg1q06c4w47lC/ERozTJD3Tcnu8+bHjmNlWMxs1s9GJCY67K1WsxvyMrnS7kF1x7U4P27h2qW7cvO5Y+Oa2ww5xhAhga/OxWXUNd9/u7sPuPrxkCYsKpYp1va2MtsQu5Ezc6c6IK89b0bbcMPOxvr7nUP8GjqyFCOBxSctbbi+T9IsAj4uQhjZJl3y+3DJARltiF3ombuuMt91jLRr444/Wj574zYJmwVzuvj5C1ID3SlpjZqsk/VzS5ZLeG+BxkbuMtsR2vcOsy8daf+Zrdd+jjVn0S0eOdt0ZUbilDVkpHMDu/rKZfUTSd9VoQ9vh7o8UHhmqIaMtsSF3mL33bSu0+6nnFtzu9vU9h/qz1blMCS68pipIH7C7j0hKd4UFKFkvM+pdBw7rR0/85tjtRQMn5LcJg7OIF4SdcEBgrS1oN25e1/Xn3f/4hF46cvTY7fVnvja/2W/MlscMcR4wEFCRFrSZC4HvfduKfg2zfzJaeE0BM2AgoCLbkkMuBEaT0cJrCghgIKCiZ03066jJUnfmZbTwGhsBDASU4iyW1rZ0EcBAYKkdmF74tDb0DYtwQMUtdJcfysMMGOgg95PNUiyLoIEABuZRlfppamURNFCCAOaxkFPSgIUigIF5UD9FP1GCAOZB/bQPOKznGAIY6ID6aUAc1nMcShAAypPRVVLKQAADKA+H9RyHEgSA8nBYz3EIYADl4rCeYyhBAEAkBDAAREIAA0AkBDAAREIAA0AkBDAAREIAA0AkBDAAREIAA0AkBDAARMJWZFRO7tdwQ32kMwMeG5Hu3tb4E+jR9DXcvvbjQ7r2Gw9q14HDsYcEzCmNAJ4+pHnvzY0/CWH0iGu4ISdpBDCHNCMQruGGnKRRA169oXF5kqlJDmlGIVzDDTkxdy/9SYeHh310dPT4D9btQn11e71IH9+TwZjZPncf7ni/ZAK4TlovTDi4uPYXJkQC+J4MqtsATqMGXDfUvJEaviejIIBj4MKESA3fk1GksQhXN1yYEKnhezIKasAAEBg1YABIHAEMxMDWe4gABsrH1ns0EcBA2Wj5QhMBDJSNli800YYGlI2WLzQVCmAz+5ykv5X0kqQnJV3l7r8NMTCg0oY2zR+8nMtQC0VLELskrXP3N0l6TNL1xYcE1ByLdLVRKIDd/Xvu/nLz5m5Jy4oPCag5FulqI+Qi3BZJ357rH81sq5mNmtnoxERGVymgXxNlY5GuNjpuRTazeySd2uafbnD3nc373CBpWNK7vYu9zdlsReaIPsRCDThr3W5F7rgI5+4XdXiiD0i6VNKF3YRvVtq9FeSHAWXotEiHSihUgjCziyV9QtJl7v5CmCElhLeCAPqoaB/wP0t6paRdZiZJu9396sKjSgX9mgD6qFAAu/uZoQaSLN4KAugTtiIDQCQEMABEQgADQCQEcArY7AHUEgEcG/v+gdoigGNj3z9QWwRwbGz2AGqLA9ljY7NH/3GuAhJFAKeAzR7903qg0v5bOVAJSaEEgWqjxo6EEcCoNmrsSBglCFQbNXYkjABG9VFjR6IoQQBAJAQwAERCAANAJAQwAERCAANAJAQwAERCAANAJAQwAERCACMdXBkENUMAIw1cGQQ1RAAjDZxahhoigJEGTi1DDXEYD9LAqWWoIQIY6eDUsjxwiadgKEEA6B6LpUERwAC6x2JpUAQwgO6xWBoUNWAA3WOxNCgCGMDCsFgaDCUIAIiEAAaASAhgAOmpycFMBDCAtNSo15gABpCWGvUaE8AA0lKjXmPa0ACkpUa9xgQwgPTUpNeYEgQAREIAA0AkBDCA+dWkJzcGAhjA3GrUkxtDkAA2s21m5mb2uhCPByARNerJjaFwAJvZckkbJT1dfDgAklKjntwYQrShfVHSxyXtDPBYAFJSo57cGAoFsJldJunn7v6QmXW671ZJWyXp9NNPL/K0AMpUk57cGDoGsJndI+nUNv90g6RPSXpnN0/k7tslbZek4eFhX8AYAaCSOgawu1/U7uNm9kZJqyRNz36XSXrAzM5192eDjhIAKqjnEoS7/0TSKdO3zexnkobd/dcBxgUAlUcfMABEEuwwHndfGeqxAKAOmAEDQCQEMABEQgADQCQEMABEQgADQCQEMABEQgADQCQEMABEQgCjPS5DA/QdAYzZuAwNUAoCGLNxGRqgFAQwZuMyNEApgh3GgwrhMjRAKQhgtMdlaIC+owQBAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQCQEMAJEQwAAQSeEANrNrzOxRM3vEzP4pxKAAoA5eUeSTzeyvJW2W9CZ3/4OZnRJmWEDT2Ij05L3S6g3S0KbYowGCKhTAkj4k6bPu/gdJcvdfFR8S0DQ2It21RZqalB74N2nVBdLwVQQxKsPcvfdPNtsvaaekiyW9KGmbu++d475bJW1t3lwn6eGenzg9r5P069iDCCiJ17PyZFv+2sXHv6ty19Gf/c6fem7Sf7fAh0viNQXE60nbWe7+6k536jgDNrN7JJ3a5p9uaH7+n0p6u6RzJP27mZ3hbVLd3bdL2t58zFF3H+703Lng9aSvaq+J15M2Mxvt5n4dA9jdL5rnST4k6ZvNwP0fMzuqxm+yiW4HCgB1VbQL4j8lbZAkM3u9pEWq1tsIAOibootwOyTtMLOHJb0k6QPtyg9tbC/4vKnh9aSvaq+J15O2rl5PoUU4AEDv2AkHAJEQwAAQSdQAruI2ZjPbZmZuZq+LPZYizOxzZjZmZv9rZv9hZifHHlMvzOzi5vfYE2b2ydjjKcLMlpvZfWZ2sPkzc13sMYVgZgNm9qCZfSv2WEIws5PN7M7mz89BMztvrvtGC+AZ25j/QtLnY40lFDNbLmmjpKdjjyWAXZLWufubJD0m6frI41kwMxuQ9C+S/kbSWklXmNnauKMq5GVJH3X3N6jRe//hzF/PtOskHYw9iIBukvQddx+S9GbN89pizoCruI35i5I+Lin7lU13/567v9y8uVvSspjj6dG5kp5w96fc/SVJt6vxSz9L7v5Ld3+g+ffn1fjBPi3uqIoxs2WSLpF0S+yxhGBmr5H0V5K+Iknu/pK7/3au+8cM4NdLOt/M9pjZD8zsnIhjKczMLpP0c3d/KPZY+mCLpG/HHkQPTpP0TMvtcWUeWNPMbKWksyXtiTuSwr6kxqTlaOyBBHKGGhvRvtosq9xiZifNdeeifcDzCrWNORUdXs+nJL2z3BEVM9/rcfedzfvcoMZb39vKHFsg1uZjyX5/dcvMXiXpLkl/7+6/jz2eXpnZpZJ+5e77zOyC2OMJ5BWS3irpGnffY2Y3SfqkpE/Pdee+qdo25rlej5m9UdIqSQ+ZmdR4u/6AmZ3r7s+WOMQFme//jySZ2QckXSrpwpR/Mc5jXNLyltvLJP0i0liCMLNBNcL3Nnf/ZuzxFLRe0mVmtknSiZJeY2a3uvv7Io+riHFJ4+4+/c7kTjUCuK2YJYjKbGN295+4+ynuvtLdV6rxP+GtKYdvJ2Z2saRPSLrM3V+IPZ4e7ZW0xsxWmdkiSZdL+q/IY+qZNX67f0XSQXf/QuzxFOXu17v7subPzOWS7s08fNX8mX/GzM5qfuhCSQfmun9fZ8Ad9LqNGeX4Z0mvlLSrOavf7e5Xxx3Swrj7y2b2EUnflTQgaYe7PxJ5WEWsl/R+ST9pHgUrSZ9y95GIY8Js10i6rflL/ylJV811R7YiA0Ak7IQDgEgIYACIhAAGgEgIYACIhAAGgEgIYACIhAAGgEj+H/DrbUn+rxikAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.scatter(X[:int(w*N),0], X[:int(w*N),1], s=10, color=cmap(0))\n",
    "plt.scatter(X[int(w*N):,0], X[int(w*N):,1], s=10, color=cmap(1))\n",
    "plt.xlim(-6,6)\n",
    "plt.ylim(-6,6)\n",
    "plt.tight_layout()\n",
    "#plt.savefig('data.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1 1]\n",
      "[ 1.17613383  1.07267182] 0.0327877614257\n"
     ]
    }
   ],
   "source": [
    "max_iter = 100\n",
    "e = 1e-8\n",
    "mu_prior = np.array([0.0,0.0])\n",
    "var_prior = b\n",
    "mu_new = np.copy(mu_prior)\n",
    "var_new = np.copy(var_prior)\n",
    "\n",
    "mu = np.zeros((N,D))\n",
    "var = np.ones(N) * 999.0\n",
    "scale = np.sqrt((2 * np.pi * var)**D)\n",
    "evidence = np.zeros(max_iter)\n",
    "\n",
    "for i in range(max_iter):\n",
    "    for n in range(N):\n",
    "        var_cav = var[n] * var_new / (var[n] - var_new)\n",
    "        mu_cav = mu_new + (mu_new - mu[n]) * var_cav / var[n]\n",
    "        \n",
    "        z_norm = ( (1 - w) * gaussian_pdf(X[n], mu_cav, (1+var_cav)) \n",
    "                  + w * gaussian_pdf(X[n], np.zeros(D), a) )\n",
    "        \n",
    "        rho = 1 - w * gaussian_pdf(X[n], np.zeros(D), a) / (z_norm+e)\n",
    "        var_new = ( var_cav - rho * var_cav**2 / (var_cav + 1) \n",
    "                   + rho * (1 - rho) * var_cav**2 * np.dot(X[n] - mu_cav, X[n] - mu_cav) / (D * (var_cav + 1)**2) )\n",
    "        mu_new = mu_cav + rho * var_cav * (X[n] - mu_cav) / (var_cav + 1)\n",
    "        \n",
    "        var[n] = var_new * var_cav / (var_cav - var_new)\n",
    "        mu[n] = mu_cav + (var[n] + var_cav) * (mu_new - mu_cav) / var_cav\n",
    "        scale[n] = z_norm / ( np.sqrt((2 * np.pi * var[n])**D) * gaussian_pdf(mu[n], mu_cav, (var[n]+var_cav)) )\n",
    "        \n",
    "print(theta)\n",
    "print(mu_new, var_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGV5JREFUeJzt3XuQXGWZx/HfM7dkJpeZXAjZIgSSgAIbEoQRErO77opaRLJYayxES7S0ytRaK4W7ChqQqv1jMRpWgQKrtiJiFYjlugEvJUHFUrcKYbImCEQIKolc5ZIMmSRkksyln/1jpsNk0j3TM+d0P6e7v5+/MjNn3vN2debX73nO+77H3F0AgMpriO4AANQrAhgAghDAABCEAAaAIAQwAAQhgAEgSCoBbGYdZrbZzJ42s51mtjKNdgGgljWl1M6tkn7q7h80sxZJbSm1CwA1y5IuxDCzmZIel7TYWdUBACVLYwS8WNIeSd82s+WStku62t0PjTzIzNZJWidJ06ZNu+Css85K4dQAkD3bt2/f6+4njXdcGiPgTkldkla5+1Yzu1XSAXe/odjvdHZ2+rZt2xKdFwCyysy2u3vneMelcRPuRUkvuvvW4a83Szo/hXYBoKYlDmB3f0XSC2b21uFvXSzpqaTtAkCtS2sWxFWS7hmeAbFb0idSahcAalYqAezuj0kat94BAHgTK+EAIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDBSl8u5enr75O7RXQEyrSm6A6gtuZzr2nufUNfubq1YPEcb1y5TQ4NFdwvIJEbASNWBI/3q2t2t9tZmde3u1oEj/dFdAjKLAEaq2lubtWLxHO0/3K8Vi+eovbU5uktAZlGCQKrMTBvXLtOBI/1qb22WGeUHoJjURsBm1mhmvzOzn6TVJqpTQ4Opo62F8AXGkWYJ4mpJO1NsDwBqWioBbGYLJF0q6Y402gOAepDWCPgWSddKyhU7wMzWmdk2M9u2Z8+elE4LANUrcQCb2RpJr7n79rGOc/dN7t7p7p0nnXRS0tMCQNVLYwS8StJlZvaspO9JepeZfSeFdgGgpiUOYHdf7+4L3P10SVdI+qW7fzRxzwCgxrEQAwCCpLoQw91/LenXabYJALWKETAABCGAASAIAQwAQQhgAAhCAANAEAIYqGI8/qm6sR8wUKV4/FP1YwQMVCke/1T9CGCgSvH4p+pHCQKoUjz+qfoRwEAVyz/+CdWJEgQABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIEgmApg9TQHUo/C9INjTFEC9Ch8Bs6cpgHoVHsDsaQqgXoWXINjTFEC9Cg9giT1NAdSn8BIEANQrAhgAghDAABCEAAaAIAQwMoeVkagXmZgFAeSxMhL1hBEwMoWVkagnBHCNqfbLd1ZGFlbt7ysKowRRQ2rh8p2VkSeqhfcVhTECriG1cvmeXxlZjeFbjpFqrbyvOBEBXEO4fI+VH6muue0hXbP5CeVy6YQw72vtogRRQ7h8j1VopJrGHie8r7WLEXCNqebL92pXzpEq72ttYgQMpISRKiaKAAZSxNaqmAhKEAAQhAAGUsJiCUxU4gA2s1PN7FdmttPMnjSzq9PoGDBS1sOtXFPQUNvSqAEPSPqcuz9qZjMkbTezB939qRTaBqpiJVi5pqChtiUeAbv7y+7+6PC/D0raKemUpO0CeZErwUodeVd6sUTWrwhQmlRnQZjZ6ZLeJmlrgZ+tk7ROkhYuXJjmaVHj8uGWHwFXaiVYfuT9yK5uXXBah26+/Dw1NhYesxSagpbLeVmmpFXDFQFKk1oAm9l0SfdK+qy7Hxj9c3ffJGmTJHV2dvKxjZJFza89cKRfj+zq1qGj/bp/x8uSTLd86LyiYTdyClo5Q5JyR+1IZRaEmTVrKHzvcff70mgTGKncK8EKXdK3tzbrgtM6dPDogKZPadb2514vufxRzrIJe0PUjsQjYBv6i/iWpJ3u/vXkXQIqq9ho1cx08+XnSTJtf+51rVwyt+SwK2fZZCJXBOUqgyAdaZQgVkm6UtIOM3ts+HvXufuWFNoGym6sS/rGxgbd8qHzJhxi44Vk0mAsZcUdteLsSxzA7v6QJN5VVK3xRquTXV5c7PcqFYzj1YoZHcdjLwjUvUrf5MsH48ypTXr4mb3af7hPs6ZNSf08Y32wMDrOBgIYUOU20cnlXO6uC0+fpZ8++aok6cYtO7Vx7fLUA3CsDxZmUmQDe0EAFZIfdf7j7b9R/6Br5tQmLZrbpkd2va4X9vWWZVFFsdkjzKTIBkbAwLBy10RHjjoffb5HFy6ao+3P7VNTo+kj3+zSyiVzK1YKYO/ibGAEDKgym+mMHHWuXDJHN1++XN/91EUaGMypo62lIsusR8535ikb8RgBA3pz1duMqY16ZFfxmmixUXIpo+f8qLPncJ9MQ+WBhbPbtHLJ3Ioss87lXNdsflwP7+rWqjPmlKXujIkhgAFJM6Y0qanR9KfX3tCCWW2aMeXEP41iMwfe3DNiry44bbZuvnx50T0jJOnLW54+ro1KlQJ6evv0wO9fUd/AoLbseEXXrT5bs6enP/sCpaMEAUg6eHRAA4M5nTFvhgYGczp4dOCEY4otLx4aPe/VG0cHdf+Ov+hfv/9Y0RJGoTYqVQrIN58/D5WHeAQwoKH67Molc3Xo6EDRJcfFZg60tzbr/IWz9MaRfk2f0qTtz/UUreVGzj7oaGvR6qXzddL0KVq9dD7TzjKAEgSg0mYFFDvGXWpuNLU0Nap/0LVi8ezjgrWvb1B/eO2gzjp5ug715/TVD5yrg0cHKj77wMx00weXM/MhQwhgYFgpizFGH5PLuV7Y16utf96nRXPb1NPbry9devaxcOvrG9Tb/uNBHeobVGODNG96i1adOU8b1y6Tu7T/cF9Fw5CnNmcLAYy6N9n5vyNvvjU1Nmj/4QG944y5x+0JvO3513Wob1CSNJiTprYM1X57Dvfpxvt36uFd3XrHkjm6/tKzNWuSdWD2dKheBDDqWpI9EfI31DraWtTT26fvfmqFFs5uO/Y0jGvvfUK/+dNramwYCt/GBulIX79WnTlPnvOhGQn9g/rhYy/p4We6terMExdijBeu7OlQ3Qhg1IzJjAST7IkwcrOblUvmakFHq/YfHjp/vt3Z06dKkr52+Xl6+2mzdKg/p/bWZr1+6KiODgxq0CW5NGNqwwnnHy9c8+WPR3btPW4hByWG6kEAoyZMdiSYZOP0kTflZkxp0hfu23Gsna9+4Nxj7a46c55WLpkrM1NHvnmXmhsaJM/JJR08ktOqM4+ffTHWh8PI59U1NTaop7dvQhvGIxsIYGReKSPbyY5kk+6JkL+p1dPbd8L5R7Y7+vV8+YGnZQ2mFmvU6qXzdcOas0+YCzzWh8Ob5Y9m9fT2H1f+QPVgHjAyrdQ9GpLMry11IcRYj4Ifef5PP/tLHb3lazLTsQ+BVzds0J7bbpc0FJ5b//y6Fs1p06y2Zt2w5mzNmjblhPPnPxx+ctXf6KYPLisYzvl9JQjf6sQIGJlW6si23Lt7jVfiyJ9//+E+Hb2lS/vuvlsy6eT16/Xqhg3ad9fdmvWxK+Xux41s87Mmio3yi00bYzez2kAAI9MmUqMt5xzXUj4IGhpMs6ZNkV+3XjJp3113a99dd0uSZn3sSp28fv2xoBwZnoODrn/7n8ePPfiz1Pp1qc+FI6SziwBGpkWO9EaG10Q+CMxMJ69ffyx8JR0XvtJQeM6c2qx9vX369x8/qft3vKzpU5r1yK69qc1kYIpa9hHAY2D0kA0Rq7cKhVepHwTurlc3bDjue69u2HBcCOfbf/iZveo5PLSHxBtH+/XOt5yU2kwGHjuUfdyEK6ISG3Qjuya7a1k+fPM137N2PqVZH7tS++66W69u2HDsBt7IWQyS1NbSpEvP/Svd8qHlqX3Y89ih7GMEXEQ1jB4YoZfPZOcHm5kaZ8w8ruZ78vr1kqTGGTOPvU8j21+9dL6+dOmb09DSel+5UZd9Vo4HAY6ns7PTt23bVvHzToS765rNb16Cjp4GFI36XvklCUJ3P+53Rn9drH3e19pgZtvdvXO84xgBF5H10UM1jNCzbryAzZcccjkfc9ey0e0UarcSS6NRfagBjyHLDy2kvpdMqTX+8Y4b/fOBgVyidttbm3XRotnad6hPFy2azfta4xgBV6msj9CzrtSR5nj7MQxthjN0M61rd7de2n84Ubv5iqBr6JFB7jw6qJYxAq5iWR6hZ12pVxDFjsuPYD/yzS41NZp6eod+fuqs1kTt5pcpz57Woq7dr5f9MfWIxU041K1Sb7IVOq6nt09rbntI7a3NBfcCnmy7Wb/5i9JwEw4YR6kLPAodN3ov4JGb4YzX7sjgHX0cpaX6QgADkzBWUI41As7lXNdsflwP7+rWqjPmaOPa5WM+AYO53rWNGjAwASO3pCxUgx9v1kRPb58e+P0r2nPwiLbseEU9vX1Ff7fUGRWoXgQwUKJSpq4Vmt0wUj6r35wfXPx3C82oQG0hgIESjReu0vizKzraWrR66XydNH2KVi+df1wNePTvljqjAtWLWRBAiUqdoVDKk4zHqhFTA65+zIIAUjbeDIWxZjeMNNYsidE/i9iKE5VDACNUtY3wigUim+hgMqgBI0wt7blcSn04bWM9JBTVgQBGmIjQksoTXJXeHGmyH14Roc0HRXGUIBBmspueJ5G0VFCsZFLpFWyT2bYyokxCaWZsBDDCRCy7nWhwjQxcd40ZJpW8YTaZD6+IvYbZ33hsBDBCVfou/0SCa/To7br3nZWZMJnMh1fEFUfEOasJAYy6MpHgGj16M6msYTLROcAT/fCKuOJgc6GxEcCoO6UG1+jRW0dbS9nCZPRo+yv/dK6++IMdiWqnhQI8Yl4xc5mLSyWAzewSSbdKapR0h7t/JY12gUiFRm9mKjoPOEkwl7IPxERCjJtf1SHxNDQza5T0DUmrJZ0j6cNmdk7SdoEsKOWpI2nMZ057H4ioKX6YmDRGwBdKesbdd0uSmX1P0vslPZVC28CkVWqVXRp3+guNtpOUO7j5VR3SCOBTJL0w4usXJV00+iAzWydpnSQtXLgwhdMCxZX7EnxkuKcVdmnuA8HNr+qQRgAXemdPuAZz902SNklDu6GlcF6gqHLOPy0U7lkMO25+ZV8aS5FflHTqiK8XSPpLCu0Ck1bOpcGFwp0nVGMy0hgB/1bSmWa2SNJLkq6Q9JEU2gUmrZyX4KWUHKptlzfESBzA7j5gZp+R9DMNTUO7092fTNwzIKFyXYKXsi8wU8BQilTmAbv7Fklb0mgLqISkI9Sxwj1fopg5tUkPP7NX+w/3ada0KUm7jBrEdpSoO+Xeh7i9tVkXLZqtZ7t71XO4Xzdu2VnVex2jfAhg1J1yL1IwM13/vrPV3tqsRXPb1LX7dRZCoCACGHWnEpunz5rWolVnzNX+wwOTPgcbmdc+noqMulSJWQpJzsGNvOpW6lORGQGjLlVi3m6Sc7CXQ30ggIEMqvQz5hCD/YCBDMrPNe7p7RPrOGoXAQxk2JcfeJo6cA2jBAFkFHXg2kcAAylLa/oYdeDaRwkCSFGa08fY07f2MQIGUpR22YBtLmsbAQykiLIBJoISBJAiygaYCAIYSBmPAkKpKEGg7rHpDaIwAkZdY9MbRGIEjLrGYgdEIoBR15i1gEiUIFDXmLWASAQw6h6zFhCFEgQABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMCoK7mcq6e3T+4e3RUgWQCb2U1m9rSZPWFmPzCzjrQ6BqQtl3Nde+8TWnPbQ7pm8xPK5QhhxEo6An5Q0lJ3Xybpj5LWJ+8SUB4HjvSra3e32lub1bW7WweO9Ed3CXUuUQC7+8/dfWD4yy5JC5J3CSiP9tZmrVg8R/sP92vF4jlqb22O7hLqXFOKbX1S0n8X+6GZrZO0TpIWLlyY4mmB0piZNq5dpgNH+tXe2iwzi+4S6ty4AWxmv5A0v8CPrnf3Hw0fc72kAUn3FGvH3TdJ2iRJnZ2dFN8QoqHB1NHWEt0NQFIJAezu7x7r52b2cUlrJF3s3FoGgJIlKkGY2SWSviDpne7em06XAKA+JJ0FcbukGZIeNLPHzOy/UugTANSFRCNgdz8jrY4AQL1hJRwABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMELxhArUszS3owQmJP+Eiq7d3VqxeI42rl2mhga2iET9YASMMDyhAvWOAEYYnlCBekcJAmF4QgXqHQGMUDyhAvWMEgQABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABCGAAdSMXM7V09snd4/uSkmaojsAAGnI5VzX3vuEunZ3a8XiOdq4dpkaGiy6W2NiBAygJhw40q+u3d1qb21W1+5uHTjSH92lcRHAAGpCe2uzViyeo/2H+7Vi8Ry1tzZHd2lclCAA1AQz08a1y3TgSL/aW5tllu3yg5TSCNjMPm9mbmZz02gPACajocHU0dZSFeErpRDAZnaqpPdIej55dwCgfqQxAr5Z0rWSqmPeBwBkRKIANrPLJL3k7o+XcOw6M9tmZtv27NmT5LQAUBPGvQlnZr+QNL/Aj66XdJ2k95ZyInffJGmTJHV2djJaBlD3xg1gd393oe+b2bmSFkl6fLjgvUDSo2Z2obu/kmovAaAGTXoamrvvkDQv/7WZPSup0933ptAvAKh5LMQAgCCpBbC7n87oF0im2jaTQTKshAMyoho3k0EylCCAjKjGzWSQDAEMZEQ1biaDZChBABlRjZvJIBkCGMiQ/GYyqA+UIAAgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCAEMAAEIYABIAgBDABBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAQhgAEgCAEMAEEIYAAIQgADQBACGACCEMAAEIQABoAgBDAABCGAASAIAQwAQQhgAAhCAANAEAIYAIIQwAAQhAAGgCCJA9jMrjKzP5jZk2a2MY1OAUA9aEryy2b2D5LeL2mZux81s3npdAsAal/SEfCnJX3F3Y9Kkru/lrxLAFAfEo2AJb1F0t+a2Y2Sjkj6vLv/ttCBZrZO0rrhL4+a2e8TnjtL5kraG92JFNXa65Fq7zXxerLtraUcNG4Am9kvJM0v8KPrh39/lqQVkt4u6ftmttjdffTB7r5J0qbhNre5e2cpHawGvJ7sq7XXxOvJNjPbVspx4wawu797jJN8WtJ9w4H7f2aW09An2Z5SOwoA9SppDfiHkt4lSWb2Fkktqq3LCAAom6Q14Dsl3Tlcz+2T9PFC5YcCNiU8b9bwerKv1l4TryfbSno9VlpeAgDSxko4AAhCAANAkNAArsVlzGb2eTNzM5sb3ZckzOwmM3vazJ4wsx+YWUd0nybDzC4Z/j/2jJl9Mbo/SZjZqWb2KzPbOfw3c3V0n9JgZo1m9jsz+0l0X9JgZh1mtnn472enma0sdmxYAI9axvzXkv4zqi9pMbNTJb1H0vPRfUnBg5KWuvsySX+UtD64PxNmZo2SviFptaRzJH3YzM6J7VUiA5I+5+5na2ju/b9U+evJu1rSzuhOpOhWST9197MkLdcYry1yBFyLy5hvlnStpKq/s+nuP3f3geEvuyQtiOzPJF0o6Rl33+3ufZK+p6EP/ark7i+7+6PD/z6ooT/sU2J7lYyZLZB0qaQ7ovuSBjObKenvJH1Lkty9z917ih0fGcD5Zcxbzex/zeztgX1JzMwuk/SSuz8e3Zcy+KSkB6I7MQmnSHphxNcvqsoDK8/MTpf0NklbY3uS2C0aGrTkojuSksUaWoj27eGyyh1mNq3YwUnnAY8prWXMWTHO67lO0nsr26Nkxno97v6j4WOu19Cl7z2V7FtKrMD3Mvv/q1RmNl3SvZI+6+4HovszWWa2RtJr7r7dzP4+uj8paZJ0vqSr3H2rmd0q6YuSbih2cNnU2jLmYq/HzM6VtEjS42YmDV2uP2pmF7r7KxXs4oSM9f5Ikpl9XNIaSRdn+YNxDC9KOnXE1wsk/SWoL6kws2YNhe897n5fdH8SWiXpMjN7n6Spkmaa2Xfc/aPB/UriRUkvunv+ymSzhgK4oMgSRM0sY3b3He4+z91Pd/fTNfQmnJ/l8B2PmV0i6QuSLnP33uj+TNJvJZ1pZovMrEXSFZJ+HNynSbOhT/dvSdrp7l+P7k9S7r7e3RcM/81cIemXVR6+Gv6bf8HM8ruhXSzpqWLHl3UEPI7JLmNGZdwuaYqkB4dH9V3u/s+xXZoYdx8ws89I+pmkRkl3uvuTwd1KYpWkKyXtMLPHhr93nbtvCewTTnSVpHuGP/R3S/pEsQNZigwAQVgJBwBBCGAACEIAA0AQAhgAghDAABCEAAaAIAQwAAT5f5BfI+7gqRvmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "plt.scatter(mu[:,0], mu[:,1], s=5.0, color=cmap(0), alpha=0.8)\n",
    "plt.scatter(mu_new[0], mu_new[1], color=cmap(3), marker='x')\n",
    "#plt.scatter(theta[0], theta[1], color=cmap(3))\n",
    "plt.xlim(-6,6)\n",
    "plt.ylim(-6,6)\n",
    "plt.tight_layout()\n",
    "#plt.savefig('ep_result.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}