add variation inference training given to core team

trainings/variational_inference_core_team/.gitignore

@@ -0,0 +1,4 @@
+.ipynb_checkpoints/
+.idea/
+__pycache__/
+data/

trainings/variational_inference_core_team/expectation_maximization.ipynb

@@ -14,12 +14,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAFlCAYAAADxilWiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3RUdf7/8eedNNILKbQgLSF0pEuXSBOpFlC/KK7r6lrBtf/c9csXXdsiWLCCiKIioiIgJXRFQAjSQ5eSACkE0kmd+/sjyIoECJCZO0lej3M4x9yZzPuFejJ5zf3czzVM00RERERERETEWWxWBxAREREREZHqRUVUREREREREnEpFVERERERERJxKRVREREREREScSkVUREREREREnEpFVERERERERJzK3arBoaGhZoMGDawaLyIiIiIiIg60adOmE6ZphpX1mGVFtEGDBsTHx1s1XkRERERERBzIMIzDF3pMS3NFRERERETEqVRERURERERExKlUREVERERERMSpVERFRERERETEqVRERURERERExKlUREVERERERMSpVERFRERERETEqVRERURERERExKlUREVERERERMSpVERFRERERETEqVRERURERERExKlUREVExGnsdrvVEURERMQFqIiKiEiFW7VqFddffz0tWrTgmmuuoWbNmnh5eTFixAgANm/erFIqIiJSjblbHUBERKqOU6dOERwczOHDh7n//vtp2bIlfn5+Z/94eXlht9t57rnn+OSTT/D09CQgIAA3Nzero4uIiIgT6YyoiIhUCNM0GTJkCLt37+buu+9m1KhRtGzZkgYNGhAaGkqNGjUwDAObzcaiRYuIiIhg0qRJtGjRgq+++srq+CIiIuJEKqIiInJVEhISGDVqFEVFRSxfvpyYmJhyf+/48eOZMmUKzz//PC+//LIDU4qIiIgrUREVEZErdvz4cQYMGEDv3r3x8PDA09Pzsr7fMAxiY2NZvXo1M2bMYMKECQ5KKiIiIq5E14iKiMgVycvLY8iQIdx///088MADV/VaderUYdWqVfTp0wcvLy+eeuqpCkopIiIirkhFVERELpvdbmf06NE0b96c5557rkJes1atWqxatYoTJ05QUlKCzWbDMIwKeW0RERFxLVqaKyIil+2FF14gLS2NDz/8sELLYnh4OM2bN2f8+PFMnz69wl5XREREXIthmqYlgzt06GDGx8dbMltERK7Oli1bqFevHqGhoQ55/ZMnT+Lt7Y2Hhwfu7lq8IyIiUhkZhrHJNM0OZT2mM6IiIlJu69at4+mnn6Zt27YOK6EAISEheHp60rlzZ7Zs2eKwOSIiImINfcwsIiLl1qpVq8veGfdKubm5MW7cOO644w7i4+Px8fFxylwRERFxPJ0RFRGRcnnnnXc4duwY7du3d9rMO++8kzZt2mgXXRERkSpGRVRERC4pISGB8ePHExwc7NS5hmHw3nvvMX/+fBYuXOjU2SIiIuI4KqIiInJRpmkyduxYnn/+ecLCwpw+PygoiE8//ZS//vWvpKamOn2+iIiIVDwVURERuaj58+eTmJjIgw8+aFmGXr16MWHCBAoLCy3LICIiIhVHRVRERC6ooKCAxx9/nMmTJ+Ph4WFplnvvvZcaNWqwePFiS3OIiIjI1VMRFRGRC5o2bRrNmzenf//+VkcBID09nQ0bNlgdQ0RERK6SYZqmJYM7dOhgxsfHWzJbRETKp7CwkMzMTEuuDb2YxMREIiMjrY4hIiIiF2EYxibTNDuU9ZjOiIqISJnGjx/P/v37Xa6E5ufn07lzZ7Zs2WJ1FBEREblCKqIiIlKmDh06UK9ePatjnKdGjRr885//ZNy4cVi1qkdERESujoqoiIicJy4ujtjYWAICAqyOUqb77ruPtLQ0vv/+e6ujiIiIyBVQERURkXMcPHiQO+64g6KiIqujXJC7uztvvPEGTz75JAUFBVbHERERkcukIioiIueYPHky9957L/7+/lZHuah+/foRHR3NO++8Y3UUERERuUwqoiIiclZGRgafffYZjzzyiNVRymXixIl4enpaHUNEREQuk7vVAURExHV8+OGHDBo0yCU3KSpLTEwMMTEx7N27l+joaKvjiIiISDnpjKiIiACl9wx96623+Mc//mF1lMuSnZ3N3XffTX5+vtVRREREpJx0RlRERABIT09n9OjRtG3b1uool8Xf35+1a9daHUNEREQug86IiogIpmkSEBDAyy+/bHWUKxYbG8vu3butjiEiIiLloCIqIiLs3LmTAQMGWB3jihmGQWxsLC+++KLVUURERKQcDNM0LRncoUMHMz4+3pLZIiJyvvz8fGrUqGF1jCuWlZVF48aN+fnnn7VxkYiIiAswDGOTaZodynpMZ0RFRKq53bt388ILL1TqEgoQEBDAI488wksvvWR1FBEREbkEFVERkWpu0qRJGIZhdYwK8eijj/LDDz+wf/9+q6OIiIjIRaiIiohUYzk5OcyePZsHHnjA6igVIigoiIceeogZM2ZYHUVEREQuQrdvERGpxr755ht69OhBrVq1rI5SYZ577jk8PDysjiEiIiIXUa4zooZhDDAMY49hGPsNw3jmIs/raBhGiWEYt1RcRBERcZRPPvmEu+++2+oYFcrLy4tff/2VF154weooIiIicgGXLKKGYbgBU4CBQHPgdsMwml/gea8CSyo6pIiIVLyDBw+yfft2brrpJqujVLhGjRrRr18/q2OIiIjIBZTnjGgnYL9pmr+ZplkIzAKGlvG8R4BvgNQKzCciIg7i6+vLjBkz8PLysjpKhQsJCaF169YsWrTI6igiIiJShvIU0bpA4h++Tjpz7CzDMOoCw4H3Ky6aiIg4it1uJz09nUGDBl31a2XnF7E/NYeMvEKsujd1WfLz87nzzjtJSUmxOoqIiIj8SXk2KyprT/8//6YxGXjaNM2Si90CwDCMvwF/A6hfv355M4qISAVLSkpi7NixLF68+Ipv3ZJbUMwnaw/x/uoDZOcXA+DhZhDi60mYnxd1grx5pE8UreoFVmT0cgsLC2PkyJG88847TJgwwZIMIiIiUjbjUp9eG4ZxHfC/pmn2P/P1swCmab78h+cc5L+FNRTIA/5mmubcC71uhw4dzPj4+KtLLyIiTldQXMIXvxzhnRX7Sc8tpEdUKH1iwsk8XcTJ3EJO5haSnlvInuRssk4XMa5vNA/0aoybzfn3Kt23bx9du3bl0KFD+Pr6On2+iIhIdWYYxibTNDuU9Vh5zohuBKIMw2gIHAVGAXf88QmmaTb8w7BPgAUXK6EiImKd3Nxcbr31VubNm4e7e/nv4mWaJnM2JTFp2V6OZeTT/ppgXr2lNa3qln3GM/N0Ea8u2s3rS/awck8qk25rS2SIT0X9NcolKiqKHj16MH36dB5++GGnzhYREZELu+Q1oqZpFgMPU7ob7i5gtmmaOw3DeMAwjKpxB3QRkWrk22+/xWazXVYJBZiycj9PztlGsI8n79x+Le/e2e6CJRQg0NuDl4a35IXBzdl1LIsBb/7IN5uSnH4d6ZNPPskbb7xBcXGxU+eKiIjIhZXrtxDTNBcCC/90rMyNiUzTHHP1sURExFFmzJjBAw9c3ueIs+MT+U/cXga0rMULg5tjK+d1pYZhcGOr2rSNDGL8/AT+8fVWNieeYsLQlld8berluu666xg9ejRZWVmEhIQ4ZaaIiIhcXHl2zRURkSri8OHDbNmyhcGDB5f7e1buTuXZb7fTqWEIzw9qVu4S+kd1grx598523N4pkpnrj/DlhsRLf1MFGj9+PIBL7eorIiJSnamIiohUI5999hkjR44s971DtyRm8ODnv9IkzI9XRrTCw+3K3zbcbAaP9Imic8MQXpi3g62JGVf8Wldi0KBB7N+/36kzRUREpGwqoiIi1YRpmsyYMYMxY8aU6/mHTuRyz/QNBPt6MGlkG3y9Lu+a0rK42QwmDG1JTV8v/v75Jk7mFl71a5bXmjVriIqKcto8ERERuTAVURGRasI0Td588006dChzF/VzpGUXcNfHGzCBN0deS02/8p1BLY9AHw9eHtGKtOwCHpu1mRK7c5bLurm5MXbsWHbu3OmUeSIiInJhKqIiItXEjh076N279yU3CTJNkyfnbCU1K5+Jt7ahfs2Kv+VK8zoBPNGvKT/tO8Gby/ZW+OtfSM2aNZk4caLT5omIiEjZVERFRKqJSZMmcfDgwUs+b2lCCqv2pPFA78a0vMjtWa7W0LZ1uKl1bd5asZ8Vu1McNuePHnzwQebOncvx48edMk9ERETKpiIqIlJNTJ8+nRYtWlz0OacLSxg/P4HGYb7c2r6eQ/MYhsGT/ZvSNMKfsbO2kJKV79B5UHpG9Pbbb+fdd991+CwRERG5MBVREZFq4N///jfz5s275PPeW7WfoxmnebJ/U9yvYofc8qrh4caLw1uSV1jCpKXOWaL72GOP8cEHH3D69GmnzBMREZHzqYiKiFRxpmny4Ycf0qBBg4s+79CJXN7/8Tf6t4jg2vrBzgkH1A/x4eb29Zgdn8jelGyHz4uOjqZTp07MnDnT4bNERESkbCqiIiJV3KZNm/D09KRVq1YXfI5pmoyfvxOPM/f6dLZ7uzXEx9OdVxbtdsq8t956ixEjRjhlloiIiJxPRVREpIqbM2cOt95660V3y122K5WVe9L4a49GhPlX3K1ayivQx4O7u17Dit2prD1wwuHzGjVqxJEjR9i8ebPDZ4mIiMj5VERFRKow0zSZM2cOt9xyywWfk19Uwvj5O2kU6sttHRy7QdHF3NYhkogAL15euAu7E+4teujQIRITEx0+R0RERM6nIioiUoVt2bIF0zRp27btBZ/z3qoDJJ06zRNO2qDoQmp4uPFAr8ZsP5rF/G3HHD5v+PDhDB48mLy8PIfPEhERkXOpiIqIVGERERF8+OGHF1yWm5qVz3urD9CveQTtr3HeBkUXMqBLine truncated
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFlCAYAAADmqMVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3RVVd6H8WenQxIgjdBJ6IQOCV0QkCbSBAQRsSFiLzOWcXREGSs6ioogIIqgIl2kS7eAEgQpUkSKdEKv6fv9A2ZehgEJkGTf3Hw/a7HIveec3CdrGbm/e5qx1iIiIiIiIiLimo/rABERERERERHQgCoiIiIiIiIeQgOqiIiIiIiIeAQNqCIiIiIiIuIRNKCKiIiIiIiIR9CAKiIiIiIiIh7Bz3XAxURGRtqYmBjXGSIiIiIiIpLNVq5cedBaG3WxZR45oMbExJCYmOg6Q0RERERERLKZMWbHpZbpEF8RERERERHxCBpQRURERERExCNoQBURERERERGPoAFVREREREREPIIGVBEREREREfEIGlBFRERERETEI2hAFREREREREY+gAVVEREREREQ8ggZUERERERER8QhZGlCNMe2MMZuMMVuMMc/8yXoJxpgMY0z3K91WRERERERE8rfLDqjGGF9gKNAeiANuNcbEXWK914G5V7qtiIiIiIiISFb2oNYHtlhrt1prU4HxQOeLrPcwMBk4cBXbioiHs9a6ThARERERL5eVAbUksPO8x7vOPfcfxpiSQFdg+JVue9736G+MSTTGJCYlJWUhS0Ryw9atW7nrgUcJKRxO3RY3smnfCdavX8+ZM2dcp4mIiIiIl8nKgGou8tyFu1LeAZ621mZcxbZnn7R2hLU23lobHxUVlYUsEckpmZmZTJw6ndpNWlGlZl2mJO4k8rY3SKp5O23fWUqrO/7Ck6NmM+3b1YwdO851roiIiIh4Cb8srLMLKH3e41LAngvWiQfGG2MAIoEbjTHpWdxWRDzIgePJPPjmOL4e+TohdW7kuuef4ub65WlbLZrktEy+/S2JpeUGMmf7Yaau/YnDU19izcYtDH55oOt0EREREcnjsjKgrgAqGmNigd1AL6D3+StYa2P//bUx5hNghrV2mjHG73LbiojnWLPhN3o8+z7EtWXA2xPoUqckVYsX+s/y0CC4uW4pbq5bipMp6Xy/5SCvFQrl/ZHPkJKSzJDBr3LugyoRERERkSt22QHVWptujHmIs1fn9QVGW2vXG2MGnFt+4Xmnl902e9JFJDudSknnySm/ciwzkNG96lCvbNifrh8S6EfbasUoF3UD95nBjBr3DKmpqQwb8paGVBERERG5KsYTr8wZHx9vExMTXWeI5BuHjh6jSbd7Sa7Tkzd61ad5pSs7D3zz/hMMGLWEHeOe5d7bevD2Gy/nUKmIiIiI5HXGmJXW2viLLcvKIb4i4sXOJKdQp1k7DhPKO51rX/FwClApOpTh/Zpzv3mNRanH2br/ODFRIfj4ZOU6bCIiIiIiZ+ndo0g+lpGRQULbbhw8Yxn4xjt0rHXRu0BlSaXoUIbdcx22cAma3Xo/740YnY2lIiIiIpIfaA+qSD52452PsWXrNp741xhub1zumr9fpehQ3u9dh37HD/IjZbHW6nxUEREREcky7UEVyafmLlvN/Elj6Pn0WzzSplq2fd9K0aE83CGBxYlr6HH3Q9n2fUVERETE+2lAFcmHrLW8NX0lxVvfw997NMn2vZxd65akaoVYZnw1hRmz5mbr9xYRERER76UBVSQfGjt/JZvSI3nmsQcoXMA/27+/n48Pf+tUh0Kt7uPOe+8jOTk5219DRERERLyPBlSRfCY5LYNnBr5C4YPr6Fr36i+KdDk1ShWm582dSA4tyZPPvZhjryMiIiIi3kMDqkg+M3LJFvyb3s3Lj96JXw7fBubB6ytQ+sb7GTF8GJs2bcrR1xIRERGRvE8Dqkg+8tvO/fz11tY0KhVE/XIROf56YcEBPNKpIcENevD0y//K8dcTERERkbxNt5kRyUduffApAktW4clOdXPtNTvXLsm0Dr3ZfSqNE8lphAZl/zmvIiIiIuIdtAdVJJ/4aulKVs2fRr/HnqVUWMFce11fH8PTN8ax/8ABasQ3IjMzM9deW0RERETyFg2oIvlAZqZlwMOPU6xJdx7qEJ/rr1+tRGG6NK6GbXIvB0+m5vrri4iIiEjeoAFVJB/4YPJ8krZt4Pln/kpwoJsj+29vVBZCIhnw7CtYa500iIiIiIhn04Aqkg+MnjyHkk260qlejLOGMuEFua5qcWZ+MYpFS7911iEiIiIinksDqoiXW7/nGIdjWvLI40/g62OctvRpGEtIvS489Y+XnXaIiIiIiGfSgCri5R5+9p+kblpCp1olXadQu3QR6t7QhTUrf2LDho2uc0RERETEw2hAFfFih06msL1QTW5q3ZJCBdzf3sUYQ5+mlShYqz1Pv/iq6xwRERER8TAaUEW82KtjviblzCnubpfgOuU/WlUtSmyzrsz5eioHDhxwnSMiIiIiHkQDqoiXSsvI5MM3XiDG9yixkcGuc/7D39eH3tfXoFCLe9mw+4jrHBERERHxIBpQRbzUu+PncOboQR67u5frlP/RpXZJIuvcwGfLt3P69GnXOSIiIiLiITSginipIe++R+mmXWhaOdp1yv8oVMCfjjVLMH7YYBZ9/6PrHBERERHxEBpQRbzQgp83s/uXb7mvXz98jNtby1xKz4TSRNz4OBttCdcpIiIiIuIhNKCKeKHnXnuHQtWacUvTqq5TLql0eEGuqxTJ6y8+x+y581zniIiIiIgH0IAq4mUOHE9my9FMbr69H8GBfq5z/tStCWVILxjBK+8Md50iIiIiIh5AA6qIl/l4ySZC6nbkga7Xu065rDplilCpURt+XPINx48fd50jIiIiIo5pQBXxIpmZlsHPPU7Z05soFVbQdc5lGWPo2KAKfqWqM+rTz13niIiIiIhjGlBFvMiP2w4T3PYx7ujR0XVKlrWvXpyQ6i0ZNnK06xQRERERcUwDqogXGfzhp7BvAy3j8s6VcUuGFaBBsxvY8ftmtm7d6jpHRERERBzK0oBqjGlnjNlkjNlijHnmIss7G2PWGGNWG2MSjTFNz1u23Riz9t/LsjNeRP5fSnoGcz7/kBrFggny93Wdc0VurF2aiK7PczA90HWKiIiIiDh02QHVGOMLDAXaA3HArcaYuAtWWwDUstbWBu4GRl2wvIW1tra1Nj4bmkXkIsbN/ZGUI/u5+5a8c3jvv7WqGk1w6SoM/XIO1lrXOSIiIiLiSFb2oNYHtlhrt1prU4HxQOfzV7DWnrT//64yGNA7TJFcNvyjMUTUaE6DClGuU65Y4QL+NCkfydcTP2Pvvv2uc0RERETEkawMqCWBnec93nXuuf9ijOlqjNkIzOTsXtR/s8A8Y8xKY0z/S72IMab/ucODE5OSkrJWLyIAHDuTyi9LZtKmc3f8fPLmqeU31ixOoRv/wubjebNfRERERK5dVt4Jmos89z97SK21U621VYAuwKDzFjWx1tbl7CHCDxpjml3sRay1I6y18dba+KiovLcHSMSlDybMw2ZmcleXG1ynXLXG5SMJ8c3glnbNOXPmjOscEREREXEgKwPqLqD0eY9LAXsutbK1dilQ3hgTee7xnnN/HwCmcvaQYRHJRku2HKJ8+7upVqKw65SrFuDnQ5uaZThlCvDlpCmuc0RERETEgawMqCuAisaYWGNMANALmH7+CsaYCsYYc+7rukAAcMgYE2yMCT33fDDQBliXnT+ASH63+/ApNiQXok/v3pz7Ncyz2tcoRoFqLXl3+EeuU0RERETEgcsOqNbadOAhYC6wAZhgrV1vjBlgjBlwbrVuwDpjzGrOXvG357mLJkUD3xljfgF+AmZaa+fkxA8ikl99MOkbDkwcSNvqxVynXLMaJQtTIaEl61YnsnfvXtc5IiIiIpLL/LKykrV2FjDrgueGn/f168DrF9luK1DrGhtF5E+sTo7g+kfepkx4Qdcp18wYQ4e6Mawv34BhH33CS8/9zXWSiIiIiOQiXS5TJA9bu/0A3332Du1qlnKdkm3aVS9GwarNGfv5l65TRERERCSXZWkPqoh4pjdGfk7a/t9pHRftOiXblAkvSL1GTThaqqjrFBERERHJZdqDKpJHWWuZNW0S1Zt3ICIk0HVOtmpdrSR/pBRgzpJlrlNEREREJBdpQBXJo77fsJMjW37mjt63uE7Jds0rRZGWtIOPvpx++ZVFRERExGvoEF+RPGrYZ1MIKhVH2zrlXKdku5JhBajR4DpsQX/XKSIiIiKSi7QHVSSPWjB7BlUbtCQ0yDuHuOsrRbHwy5G8/+Eo1ykiIiIikks0oIrkQVsOnCSlcGlu7tLJdUqOaV45Ct9C0YwaM851ioiIiIjkEg2oInnQnLW7KdygOx0bV3edkmMqFg0htk5Tfl21gqNHj7rOEREREZFcoAFVJA8a8uqLBG9dRHShINcpOcYYQ6uaZQgoXYOJU6e5zhERERGRXKABVSSPOXA8mbR6vbilV2/XKTmueaUogio2ZvTYL12niIiIiEgu0IAqksd8OmcZpzYto3XNMq5TclzNUkUoXqMxK5d/y8mTJ13niIiIiEgO021mRPKYcZ99RsDRE5SLCnadkuN8fQwtapXneNdnSLeua0REREQkp2kPqkgecjIlnQ3LF9KqbQeMMa5zckXzylFQogZTFq90nSIiIiIiOUwDqkge8uWCRNJPHqFnh5auU3JNQkwYgRmn+Ofzf8Na7UYVERER8WY6xFckDxkzfiJFqjSkTtkI1ym5JtDPl2a1KrIq7EUyLfjmjx3HIiIiIvmS9qCK5BGp6ZmsXPoNjVq0xdcnf01pzSpFsXvTKvo/+qTrFBERERHJQRpQRfKIZb8nQcFwbu3awXVKrmtSPpLAwkUZ/+lHpKenu84RERERkRyLine truncated
       "text/plain": [
        "<Figure size 1152x432 with 1 Axes>"
       ]
@@ -37,8 +37,8 @@
     "\n",
     "np.random.seed(654)\n",
     "# Draw samples from two Gaussian w.p. z_i ~ Bernoulli(phi)\n",
-    "generative_m = np.array([stats.norm(2, 1), stats.norm(5, 1.8)])\n",
-    "z = stats.bernoulli(0.75).rvs(100)\n",
+    "generative_m = np.array([stats.norm(-1, 1), stats.norm(5, 1.8)])\n",
+    "z = stats.bernoulli(0.9).rvs(100)\n",
     "x_i = np.array([g.rvs() for g in generative_m[z]])\n",
     "\n",
     "# plot generated data and the latent distributions\n",
@@ -72,8 +72,8 @@
     "## Continuous form\n",
     "\n",
     "$$ \\begin{eqnarray}\n",
-    "\\ell(\\theta) & = & \\int \\log P(x | \\theta, z) P(z) dz \\\\\n",
-    "\\log P(x; \\theta)& = & \\int \\log P(x | \\theta, z) P(z) dz\n",
+    "\\ell(\\theta) & = & \\int_z \\log P(x | \\theta, z) P(z) dz \\\\\n",
+    "\\log P(x; \\theta)& = & \\int_z \\log P(x | \\theta, z) P(z) dz\n",
     "\\end{eqnarray}$$\n",
     "\n",
     "\n",
@@ -86,6 +86,9 @@
    "source": [
     "**Jensens inequality:**\n",
     "*Only for concave functions*\n",
+    "\n",
+    "let $f$ be a concave function.\n",
+    "\n",
     "$$ E[f(X)] \\le f(E[X]) $$\n",
     "\n",
     "**Jensens equality:**\n",
@@ -128,7 +131,7 @@
    "source": [
     "# Now what?\n",
     "\n",
-    "Jensen's inequality holds to equaltiy if and only if $X = E[X]$. This will be the case if we choose a distribution that leads to:\n",
+    "Jensen's inequality holds to **equaltiy** if and only if $X \\iff E[X]$. This will be the case if we choose a distribution that leads to:\n",
     "\n",
     "$$ \\begin{eqnarray}\n",
     "\\frac{P(x|z) P(z)}{Q(z)} &=& \\frac{P(x, z)}{Q(z)} = c\n",
@@ -138,17 +141,17 @@
     "\n",
     "$$ P(x, z; \\theta) \\propto Q(z) $$\n",
     "\n",
-    "A next requirement is that $\\int Q(z) dz = 1$ (requirement of probability distribution).\n",
+    "A next requirement is that $\\int_z Q(z) dz = 1$ (requirement of probability distribution). And because $Q(z)$ is a probability distribution, it must integrate to one; $\\int_z Q(z) dz = 1$. So ideally we want to take joint distribution $P(x, z; \\theta)$ and transform it so that it is proportional to itself, but also sums to 1 over all values of $z_i$. It turns out we can find that by normalizing the joint distribution.\n",
     "\n",
     "$$ \\begin{eqnarray}\n",
-    "\\int \\frac{P(x, z; \\theta)}{\\int P(x, z; \\theta) dz} dz &=& 1 \\\\\n",
+    "\\int_z \\frac{P(x, z; \\theta)}{\\int_z P(x, z; \\theta) dz} dz &=& 1 \\\\\n",
     "\\end{eqnarray}$$\n",
     "\n",
     "This will lead to:\n",
     "\n",
     "$$ \\begin{eqnarray}\n",
-    "Q(z) &=& \\frac{P(x, z; \\theta)}{\\int P(x, z; \\theta) dz}  \\\\\n",
-    "&=& \\frac{P(x, z; \\theta)}{\\int P(x; \\theta)} \\\\\n",
+    "Q(z) &=& \\frac{P(x, z; \\theta)}{\\int_z P(x, z; \\theta) dz}  \\\\\n",
+    "&=& \\frac{P(x, z; \\theta)}{\\int_z P(x; \\theta)} \\\\\n",
     "&=& P(z|x; \\theta)\n",
     "\\end{eqnarray}$$"
    ]
@@ -160,40 +163,40 @@
     "## Gaussian mixture model\n",
     "\n",
     "### Expectation\n",
-    "Given a Gaussian mixture model, observed data points $x$ and current parameters $\\theta_t$ find $Q(z; \\theta_t) = P(z|x; \\theta_t)$.\n"
+    "Given a Gaussian mixture model, observed data points $x$ and current parameters $\\theta_t$ find $Q(z; \\theta_t) = P(z|x; \\theta_t)$."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array([ 4.28244273,  4.66972656,  3.35275075,  2.36158676,  3.74751768,\n",
-       "        9.75240814,  5.15756521,  4.80213993,  5.90145058,  1.91271217,\n",
-       "        4.71273852,  2.36241133,  1.61084671,  6.71697406,  3.11442326,\n",
-       "        5.15422491,  5.09580485,  2.11274369,  6.76574344,  1.48354695,\n",
+       "array([ 4.28244273,  4.66972656,  3.35275075,  5.65085616,  3.74751768,\n",
+       "        9.75240814,  5.15756521,  4.80213993,  5.90145058,  4.84288191,\n",
+       "        4.71273852, -0.63758867, -1.38915329,  6.71697406,  7.00596187,\n",
+       "        5.15422491,  5.09580485,  5.20293865,  6.76574344, -1.51645305,\n",
        "        2.65304672,  6.81892055,  2.77103106,  5.06418788,  5.72517362,\n",
-       "        3.54313443,  5.67858554,  1.03222311,  1.63983108,  5.51107055,\n",
-       "        0.65124788,  4.97843073,  4.32989217,  2.28812021,  4.51797044,\n",
-       "        4.52825271,  6.93340533,  4.94273792,  2.54969132,  0.78634279,\n",
-       "        4.81172184,  5.19576022,  6.11967179,  3.5235398 ,  0.75307613,\n",
-       "        4.23428727,  7.15230773,  2.83345813,  8.60608829,  4.27073679,\n",
-       "        7.75388195,  6.79259227,  0.62404892,  4.26260727,  4.09584051,\n",
-       "        3.26993383,  4.86431429,  8.38977374,  6.48961063,  4.09783105,\n",
+       "        3.54313443,  5.67858554, -1.96777689,  4.35169594,  5.51107055,\n",
+       "        2.57224619,  4.97843073,  4.32989217,  5.51861638,  4.51797044,\n",
+       "        4.52825271,  6.93340533,  4.94273792,  2.54969132,  2.81541702,\n",
+       "        4.81172184,  5.19576022,  6.11967179,  3.5235398 , -2.24692387,\n",
+       "        4.23428727,  7.15230773,  6.50022463,  8.60608829,  4.27073679,\n",
+       "        7.75388195,  6.79259227, -2.37595108,  4.26260727,  4.09584051,\n",
+       "        0.26993383,  4.86431429,  8.38977374,  6.48961063,  4.09783105,\n",
        "        4.40498267,  3.27433241,  7.44164281,  7.62286463,  3.21797707,\n",
        "        4.42538064,  7.06135769,  3.94733606,  1.91385988,  5.23226903,\n",
-       "        3.61125312,  1.86866663,  5.71568321,  2.50687193,  2.26842535,\n",
-       "        3.75523264,  0.71738199,  8.02274291,  2.46995812,  5.35222504,\n",
-       "        2.72716235,  4.15793385,  5.33299762,  1.00091742,  7.04351757,\n",
-       "        6.14452194,  2.62636128,  4.55251749, -0.38101307,  4.27776483,\n",
-       "        1.57442618,  3.75797097,  0.52259715,  1.59952851,  6.13345271,\n",
-       "        1.42972786,  2.82725842,  4.57300419,  2.59096887,  3.45567689])"
+       "        3.61125312, -1.13133337,  5.71568321,  5.91236947,  5.48316563,\n",
+       "        3.75523264,  2.69128758,  8.02274291, -0.53004188,  5.35222504,\n",
+       "       -0.27283765,  4.15793385,  5.33299762, -1.99908258,  7.04351757,\n",
+       "        6.14452194,  2.62636128,  4.55251749, -3.38101307,  4.27776483,\n",
+       "        4.23396713,  3.75797097,  2.34067487,  4.27915132,  6.13345271,\n",
+       "        3.97351015, -0.17274158,  4.57300419,  6.06374396,  3.45567689])"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -250,7 +253,7 @@
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAHjCAYAAAA5Y6JnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbMElEQVR4nO3df7Dld13f8dfbXULk14aa1YHdyKZOQKJAwCXSagsa1A1aUzrqJOGHUmiMJRZrZ8r2JzvFtiJK0TEQI00DrSXjIC0rWUgdEBlwgGwAA0sauoQtWUPN8sMUQSdG3v3jnqWHm5u9Zzfnfs7eex+PmTt7v9/v55zzPpczWZ77/Z5zq7sDAAAAo3zDogcAAABgcxGiAAAADCVEAQAAGEqIAgAAMJQQBQAAYCghCgAAwFBCFAAAgKGEKAAAAEMJUQA2har691X1czOs+1BVfceImSaPN+tc+6pq3xrOMfR5A7C5CVEA1r2qenhV/UJVfaqqvlRVn6iqn546vj3JC5P8xgx398tJ/s2c5np0VXVV/VlVfaWq7pqOzpOcay6q6qlV9f7JPB+qqm+dHJrb8waA1QhRANa1qnp0kvclOTfJRUkeleQfJHllVb14suynkhzo7j+f4S73J/m+qnrMHMa7IMmx7n5Edz8syc8k+Q9VtfMU5nrQJo97IMmrknxTkjuS/MvJ4Xk+bwA4ISEKwHr3q0n+JMnzu/tIL3l/ktckeelkzcVJ/mD6RlX125Mzlce/uqqu6u6/SHJLkh+cw2wXJLl5avuDkz/PWGmuB5pppTuuql+qqv82tf3qqnpXVT2kql5XVa9b4Wa/kuQ3u3v/JH5vSPL0JJnz8waAE9q66AEA4FRV1eOSPC/J07u7lx3+VJJdk++flOT26YPd/RNT9/PzSV6Q5L9Odt2W5CkrPN7bk3zvA4zzvu7+kWX7nprkQ5PbnpXk32Yp9j690lyrzLTcq5J8qqouSPKMJHuSfG93/2WSf7jC7I9KckmSx0/t/oYkfzG1veLzBoB5E6IArGfPTnJnd394hWM7khydfH9Wki+tdAdV9bIsvU/zou7+wmT3l5Lc7xLVFUJzNRck+XuTx/hikt9P8nemonnFuR5gpuWzfL6qXpvkTUm2ZSlC7znBLBcleUiSW6vq+L6HJnnb1JoVnzcAzJsQBWA9257/H5vLPTfJuyfffzHJI5cvmFz2+uIk39/dn5869Mgkf/pgBquqhyZ5YpJzu/uBZrzfXCeYaSUfSfKKJM/r7jtXWbsryf7ufu7UYx1I8s6pNQ/6eQPALLxHFID17NNJHldVX/f3WVX9QJLdWXpPZJLcmq+/JDVV9TNJfjpLZx0/t+x+n5jkj5Y/WFW9Y9l7OKe/3rFs+Xcm+fIJIvR+c60y0/JZnpTk9UnemOTvn2jtxEOTfGXq9udm6We0f2rNis8bAOZNiAKwnt04+fMXquphVfXQqnp+kjcn+fGps4QHkjzz+I2q6ookVyV5dncfm77DyZnM70rye8sfrLsvnnwC7kpfFy9b/tQkh1aZ/2tznWim5apqR5LfTXJllt4P+qSqetbU8eur6vplN7s5yTOr6rFVdU6W3nv6L45f+nui5w0A8yZEAVi3uvvPsvTexyclOZLkz5P8fJJndveBqaVvSvKcqvrGyfYvJfm2LH3Yz/Ezmi+YHPvRJO/p7rse5HgXJPn4Kmum5zrRTF8z+dChA0leM/n0268keXWWPgjpuHOSvH/ZTd+dpXj9ZJZ+3c1/7u7fnDo+r+cNAKuq+3/IIACsT1X1E1n6dS7fsfxDfqrq3yW5u7tfu8p9fDDJi7t7tYici5OYa1+SdPe+VdadkaXLa588+QTdWecY+rwB2NyEKAAbyuR9loe6+72LnmWeZg1RAFgPhCgArAPH3wPa3e9Z7CQA8OAJUQAAAIbyYUUAAAAMtXVRD3z22Wf3rl27FvXwAAAArKFbbrnlc929faVjCwvRXbt25eDBg4t6eAAAANZQVf3vBzrm0lwAAACGEqIAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoAAAAQwlRAAAAhhKiAAAADCVEAQAAGEqIAgAAMJQQBQAAYCghCgAAwFBCFAAAgKFWDdGquq6q7q6qjz/A8aqqX6uqw1V1a1U9bf5jAgAAsFHMckb0+iR7TnD84iTnTb6uSPL6Bz8WAAAAG9WqIdrd703yhRMsuSTJm3rJB5KcVVWPmdeAAAAAbCzzeI/ojiR3Tm0fnewDAACA+5lHiNYK+3rFhVVXVNXBqjp47NixOTz0Gtu3bdETnN5W+/kcP34yP8d5/Myn7+N0+9/wRPNM/7xO5Wd3Mo/9YO/3VG4//bxO5f4eaP5Tvc/lx1faXqvXz8k+31lvczKPtxavMdhkdu29cdEjnJIHmvv4/lme16xrFvkzOtXHPpnbrcXzW6+vq2R9z8548wjRo0nOmdremeSulRZ297Xdvbu7d2/fvn0ODw3A3J3sPx6JWTaZRQcWwEYwjxDdn+SFk0/PfUaSe7r7s3O4XwCAuRGPAKePrastqKo3J3lWkrOr6miSVyR5SJJ09zVJDiR5TpLDSb6S5EVrNSwAAADr36oh2t2XrXK8k7x0bhMBsHns25bsu2fRUwAAg83j0lwAAACYmRAFANYV7/UEWP+EKAAAAEMJUQAAAIYSogCsH35nKQBsCEIUAACAoYQoAAAAQwlRANY3l+sCwLojRAEAABhKiAIAADCUEAUAAGAoIQoAAMBQQhQAAIChhCgAG5NP0z0t7dp746JHAOA0IEQBAAAYSogCsLE5MwoApx0hCsDGJ0YB4LQiRAHYHMQoAJw2hCgAAABDCVEAAACGEqIAAAAMJUQBAAAYSogCsDn58CIAWBghCgAAwFBCFAAAgKGEKAAAAEMJUQAAAIYSogAAAAwlRAEAABhKiAIAADCUEAUAAGAoIQoAAMBQQhQAAIChhCgATNu3bdETAMCGJ0QBAAAYSogCACdt194bFz0CAOuYEAUAAGAoIQoAy3mfKACsKSEKAADAUEIUAACAoYQoAAAAQwlRAAAAhhKiAAAADCVEAQAAGEqIAgAAMJQQBQAAYCghCgAAwFBCFAAAgKGEKAAAAEMJUQAAAIYSogAAAAwlRAEAABhKiAIAADCUEAUAAGAoIQoAAMBQQhQAAIChhCgAcD+79t646BEA2MCEKAAAAEMJUQA4Wfu2LXoCAFjXhCgAAABDCVEAAACGEqIAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAEiS7Np746JHAGCTEKIAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoADxY+7YtegIAWFdmCtGq2lNVt1fV4arau8LxbVX1u1X1R1V1qKpeNP9RAQAA2AhWDdGq2pLk6iQXJzk/yWVVdf6yZS9N8onufkqSZyX5lao6Y86zAgAAsAHMckb0wiSHu/uO7r43yQ1JLlm2ppM8sqoqySOSfCHJfXOdFAAAgA1hlhDdkeTOqe2jk33Tfj3JE5PcleRjSV7W3V9dfkdVdUVVHayqg8eOHTvFkQEAAFjPZgnRWmFfL9v+oSQfTfLYJBck+fWqetT9btR9bXfv7u7d27dvP+lhAYD52rX3xkWPAMAmNEuIHk1yztT2ziyd+Zz2oiRv7SWHk3w6ybfPZ0QAAAA2kllC9OYk51XVuZMPILo0yf5laz6T5KIkqapvSfKEJHfMc1AAAAA2hq2rLeju+6rqqiQ3JdmS5LruPlRVV06OX5PklUmur6qPZelS3pd39+fWcG4AAADWqVVDNEm6+0CSA8v2XTP1/V1JfnC+owEAALARzXJpLgAAAMyNEAUAAGAoIQoAAMBQQhQAAIChhCgAAABDCVEAAACGEqIAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoAAAAQwlRAAAAhhKiAAAADCVEAQAAGEqIAgAAMJQQBYBNYtfeGxc9AgAkEaIAAAAMJkQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoAKy1fdsWPQEAnFaEKAAAAEMJUQAAAIYSogAAAAwlRAEAABhKiALABrZr742LHgEA7keIAgAAMJQQBQAAYCghCgAAwFBCFAAAgKGEKAAAAEMJUQAAAIYSogAAAAwlRAEAABhKiAIAADCUEAUAAGAoIQoAAMBQQhQAAIChhCgAjLRv26InAICFE6IAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoAAAAQwlRAAAAhhKiAAAADCVEAQAAGEqIAgAAMJQQBQAAYCghCgAAwFBCFAAAgKGEKAAAAEMJUQAAAIYSogAAAAwlRAEAABhKiAIAADCUEAUAAGAoIQoAAMBQQhQAAIChhCgAAABDCVEAAACGEqIAAAAMJUQBAAAYSogCAAAwlBAFAABgKCEKAADAUEIUAACAoYQoAAAAQ80UolW1p6pur6rDVbX3AdY8q6o+WlWHquoP5jsmAAAAG8XW1RZU1ZYkVyf5gSRHk9xcVfu7+xNTa85K8roke7r7M1X1zWs1MAAAAOvbLGdEL0xyuLvv6O57k9yQ5JJlay5P8tbu/kySdPfd8x0TAACAjWKWEN2R5M6p7aOTfdMen+TRVfWeqrqlql640h1V1RVVdbCqDh47duzUJgYAAGBdmyVEa4V9vWx7a5LvSvLDSX4oyb+qqsff70bd13b37u7evX379pMeFgAAgPVv1feIZukM6DlT2zuT3LXCms9195eTfLmq3pvkKUk+OZcpAQAA2DBmOSN6c5LzqurcqjojyaVJ9i9b87Ykf6uqtlbVw5J8d5Lb5jsqAAAAG8GqZ0S7+76quirJTUm2JLmuuw9V1ZWT49d0921V9c4ktyb5apI3dPfH13JwAAAA1qdZLs1Ndx9IcmDLine truncated
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAHjCAYAAAA5Y6JnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAbtUlEQVR4nO3df5Dtd13f8dfbXCCi8qPm6mASvWknKlF+eo202oIFNUFLakcdiIIiNmKJxdpOibXVTNVWxR/oCKQRU6RaMo5ijXIRKf5g1AFzQfkRaOg1pOQaai6iFEGHRt/9Y8+1y2Zz9yQ5+z7Z3cdjZufu9/v97Dnvw5dJ8rzf7zlb3R0AAACY8nHrHgAAAICDRYgCAAAwSogCAAAwSogCAAAwSogCAAAwSogCAAAwSogCAAAwSogCAAAwSogCcCBU1X+sqm9bYt3vVdXnTMy0eL5l57q6qq7exTlGXzcAB5sQBWDPq6pPqKrvrao/rKoPVdU7q+qbNx0/nORZSf7TEg/3Q0n+/YrmenhVdVX9eVV9pKpu3xyd93Culaiqx1XV7yzm+b2q+vTFoZW9bgDYiRAFYE+rqocn+e0kFyR5cpKHJPmnSb6nqp6zWPYNSY51918s8ZA3JPniqnrECsZ7bJJT3f2J3f3gJN+S5Eer6rx7Mdd9tnjeY0l+IMknJ7klyb9dHF7l6waAMxKiAOx1P5bkj5N8XXff2ht+J8mPJHneYs2lSX5r8w9V1c8trlSe/uqqurK7/zLJm5N86Qpme2ySGzdtv2nx5wO3m+vuZtrugavqB6vqFzdtv7CqXl9VD6iql1TVS7b5sR9O8pPdfcMifq9P8vlJsuLXDQBndGjdAwDAvVVVn5Hka5N8fnf3lsN/mOTI4vtHJbl588Hu/ppNj/PtSZ6Z5L8udr0ryWO2eb5fSfJFdzPOb3f3V2zZ97gkv7f42Ycl+b5sxN57tptrh5m2+oEkf1hVj03yhCSXJPmi7v6/Sf7ZNrM/JMllST5z0+6PS/KXm7a3fd0AsGpCFIC97ClJbuvut2xz7NwkJxffPyzJh7Z7gKp6fjbep/nk7v7AYveHktzlFtVtQnMnj03yTxbP8adJfiPJP9oUzdvOdTczbZ3lT6rqRUlekeSh2YjQD55hlicneUCSt1XV6X0PSvJLm9Zs+7oBYNWEKAB72eH8/9jc6iuT/Pri+z9N8klbFyxue31Okn/Y3X+y6dAnJfmz+zJYVT0oySOTXNDddzfjXeY6w0zb+f0k353ka7v7th3WHklyQ3d/5abnOpbkVzetuc+vGwCW4T2iAOxl70nyGVX1Mf8+q6ovSXI0G++JTJK35WNvSU1VfUuSb87GVcf3b3ncRyZ569Ynq6rXbHkP5+av12xZ/rlJPnyGCL3LXDvMtHWWRyV5aZKfTvKNZ1q78KAkH9n08xdk43+jGzat2fZ1A8CqCVEA9rJXL/783qp6cFU9qKq+Lskrk3z1pquEx5I88fQPVdUVSa5M8pTuPrX5ARdXMj8vyeu2Pll3X7r4BNztvi7dsvxxSW7aYf6/metMM21VVecm+eUkz83G+0EfVVVP2nT85VX18i0/dmOSJ1bVp1XV+dl47+l3nr7190yvGwBWTYgCsGd1959n472Pj0pya5K/SPLtSZ7Y3cc2LX1FkqdW1ccvtn8wyd/Jxof9nL6i+czFsacl+c3uvv0+jvfYJO/YYc3muc40099YfOjQsSQ/svj0248keWE2PgjptPOT/M6WH/31bMTru7Px627+S3f/5Kbjq3rdALCjuuuHDALA3lRVX5ONX+fyOVs/5Keq/kOSO7r7RTs8xpuSPKe7d4rIlbgHc12dJN199Q7rHpiN22sfvfgE3WXnGH3dABxsQhSAfWXxPsubuvsN655llZYNUQDYC4QoAOwBp98D2t2/ud5JAOC+E6IAAACM8mFFAAAAjDq0ric+55xz+siRI+t6egAAAHbRm9/85vd39+Htjq0tRI8cOZLjx4+v6+kBAADYRVX1v+7umFtzAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGLVjiFbVdVV1R1W9426OV1X9eFWdqKq3VdXjVz8mAAAA+8UyV0RfnuSSMxy/NMmFi68rkrz0vo8FAADAfrVjiHb3G5J84AxLLkvyit7wxiQPq6pHrGpAAAAA9pdVvEf03CS3bdo+udgHAAAAd7GKEK1t9vW2C6uuqKrjVXX81KlTK3hq9pyrH/qxf04/736zjtc1+Zz79byxHOefu3Hkqldv+/3drVnm+JnWb3fs9L67e5ydnv/0mmXnXObxtq5b9mfuy+Mv8xrujXv7uNu9/vv6OPdlzaqea9m1q5hpN87nfXF/m4fVWUWInkxy/qbt85Lcvt3C7r62u49299HDhw+v4KkBYBPxuiftVswAcP+1ihC9IcmzFp+e+4QkH+zu963gcQHgnhOj91ubr1AJT4CD7dBOC6rqlUmelOScqjqZ5LuTPCBJuvuaJMeSPDXJiSQfSfLs3RoWAO6Vqx+aXP3BdU8BACzsGKLd/YwdjneS561sIgBgzzty1atz6/d/+brHAOB+ahW35gIAAMDShCgAcJ95zycA94QQBQAAYJQQBYDTfOIuAIwQogAAAIwSogAAAIwSogAAAIwSogAAAIwSogDAUvyKFgBWRYgCAAAwSogCAGfkSigAqyZEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAQAAGCVEAeCAOXLVq9c9AgAHnBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABglBAFAABg1FIhWlWXVNXNVXWiqq7a5vhDq+qXq+qtVXVTVT179aMCAACwH+wYolV1VpIXJ7k0yUVJnlFVF21Z9rwk7+zuxyR5UpIfrqoHrnhWAAAA9oFlrohenOREd9/S3R9Ncn2Sy7as6SSfVFWV5BOTfCDJnSudFAAAgH1hmRA9N8ltm7ZPLvZt9hNJHpnk9iRvT/L87v7rrQ9UVVdU1fGqOn7q1Kl7OTIAAAB72TIhWtvs6y3bX5bkD5J8WpLHJvmJqnrIXX6o+9ruPtrdRw8fPnyPhwUAAGDvWyZETyY5f9P2edm48rnZs5O8qjecSPKeJJ+9mhEBAADYT5YJ0RuTXFhVFyw+gOjpSW7Ysua9SZ6cJFX1qUk+K8ktqxwUAACA/eHQTgu6+86qujLJa5OcleS67r6pqp67OH5Nku9J8vKqens2buV9QXe/fxfnBgAAYI/aMUSTpLuPJTm2Zd81m76/PcmXrnY0ANiHrn5ocvUH1z0FAKzVMrfmAgAAwMoIUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQAAAEYJUQDYw45c9ep1jwAA95gQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBYA96shVr173CABwrwhRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARglRAAAARi0VolV1SVXdXFUnquqqu1nzpKr6g6q6qap+a7VjAgAAsF8c2mlBVZ2V5MVJviTJySQ3VtUN3f3OTWseluQlSS7p7vdW1afs1sAAAADsbctcEb04yYnuvqW7P5rk+iSXbVlzeZJXdfd7k6S771jtmAAAAOwXy4TouUlu27R9crFvs89M8vCq+s2qenNVPWu7B6qqK6rqeFUdP3Xq1L2bGAAAgD1tmRCtbfb1lu1DST4vyZcn+bIk/66qPvMuP9R9bXcf7e6jhw8fvsfDAgAAsPft+B7RbFwBPX/T9nlJbt9mzfu7+8NJPlxVb0jymCTvXsmUAAAA7BvLXBG9McmFVXVBVT0wydOT3LBlzS8l+ftVdaiqHpzkC5K8a7WjAgAAsB/seEW0u++sqiuTvDbJWUmu6+6bquq5i+PXdPe7qupXk7wtyV8neVl3v2M3BwcAAGBvWubW3HT3sSTHtuy7Zsv2C5O8cHWjAQAAsB8tc2suAAAArIwQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQQBQAAYJQLine truncated
       "text/plain": [
        "<Figure size 1152x576 with 1 Axes>"
       ]
@@ -263,22 +266,31 @@
    ],
    "source": [
     "torch.manual_seed(1)\n",
-    "g in better poste-rior a\n",
+    "\n",
+    "# Number of clusters\n",
     "k = torch.tensor([2])\n",
+    "\n",
+    "# Current values for theta\n",
     "mu = dist.Uniform(x_i.min(), x_i.max()).sample(k) \n",
     "std = torch.ones(k)\n",
+    "phi = torch.ones(k)\n",
     "\n",
-    "def expectation(x, mu, std, k):\n",
+    "def expectation(x, mu, std, phi, k):\n",
+    "    # P(z|x)\n",
     "    qz = torch.zeros((k, x.shape[0]))\n",
+    "    phi = F.softmax(phi, dim=-1)\n",
     "    \n",
     "    for z in range(k):\n",
-    "        qz[z] = dist.Normal(mu[z], std[z]).log_prob(x).exp() \n",
+    "        # Define the joint: P(z, x)\n",
+    "        qz[z] = (dist.Normal(mu[z], std[z]).log_prob(x) + torch.log(phi[z])).exp() \n",
     "    \n",
     "    # normalize zo that marginalizing z would lead to p = 1\n",
+    "    # make conditional by dividing: 1/ P(x)\n",
     "    qz = qz / qz.sum(0)\n",
+    "    \n",
     "    return qz\n",
     "\n",
-    "qz = expectation(x_i, mu, std, k)    \n",
+    "qz = expectation(x_i, mu, std, phi, k)    \n",
     "\n",
     "plt.figure(figsize=(16, 8)) \n",
     "plt.title(r\"$Q(z) = P(z|x; \\theta)$\")\n",
@@ -300,55 +312,62 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
-      "  after removing the cwd from sys.path.\n",
+      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:5: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
+      "  \"\"\"\n",
       "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
-      "  \n"
+      "  \n",
+      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
+      "  import sys\n"
      ]
     }
    ],
    "source": [
+    "torch.manual_seed(0)\n",
     "k = torch.tensor([2])\n",
     "mu = dist.Uniform(x_i.min(), x_i.max()).sample(k) \n",
-    "std = torch.ones(k)\n",
+    "\n",
     "mu = torch.tensor(mu, dtype=torch.float, requires_grad=True)\n",
-    "std = torch.tensor([1., 1.])\n",
-    "std = torch.tensor(std, dtype=torch.float, requires_grad=True)\n"
+    "std = torch.tensor(torch.ones(k), requires_grad=True)\n",
+    "phi = torch.tensor(torch.ones(k), requires_grad=True)\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "tensor(-540.7363, grad_fn=<AddBackward0>)"
+       "tensor(-941.1580, grad_fn=<AddBackward0>)"
       ]
      },
-     "execution_count": 31,
+     "execution_count": 14,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
+    "def log_likelihood(mu, std, phi, x):\n",
+    "    return torch.log(phi) + dist.Normal(mu, std).log_prob(x)\n",
+    "\n",
     "def get_elbo(mu, std, phi):\n",
     "    elbo = 0\n",
     "    for k in range(2):\n",
     "        z_i = torch.tensor(1.)  # prior on z\n",
     "        mu_i = mu[k]\n",
     "        std_i = std[k]\n",
+    "        phi_i = F.softmax(phi, dim=-1)[k]\n",
     "        \n",
-    "        for i in range(len(x_i)):\n",
-    "            elbo += qz[k, i] * ((dist.Normal(mu_i, std_i).log_prob(x_i[i]) + torch.log(z_i) - torch.log(qz[k, i])))\n",
+    "        elbo += (qz[k] * (log_likelihood(mu_i, std_i, phi_i, x_i) + torch.log(z_i) - torch.log(qz[k]))).sum()\n",
+    "\n",
     "    return elbo\n",
     "\n",
     "elbo = get_elbo(mu, std, phi)\n",
@@ -357,70 +376,75 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "optim = torch.optim.Adam([mu, std, phi])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "neg. ELBO: 147.22\tmu: tensor([2.8749, 5.8349])\tstd: tensor([1.4770, 1.5313])\n",
-      "neg. ELBO: 146.09\tmu: tensor([2.7886, 5.6798])\tstd: tensor([1.4575, 1.5878])\n",
-      "neg. ELBO: 145.73\tmu: tensor([2.7476, 5.5918])\tstd: tensor([1.4493, 1.6264])\n",
-      "neg. ELBO: 145.62\tmu: tensor([2.7310, 5.5414])\tstd: tensor([1.4473, 1.6521])\n",
-      "neg. ELBO: 145.59\tmu: tensor([2.7267, 5.5119])\tstd: tensor([1.4488, 1.6695])\n",
-      "neg. ELBO: 145.57\tmu: tensor([2.7280, 5.4938])\tstd: tensor([1.4517, 1.6818])\n",
-      "neg. ELBO: 145.57\tmu: tensor([2.7318, 5.4821])\tstd: tensor([1.4552, 1.6908])\n",
-      "neg. ELBO: 145.56\tmu: tensor([2.7362, 5.4739])\tstd: tensor([1.4586, 1.6978])\n",
-      "neg. ELBO: 145.56\tmu: tensor([2.7406, 5.4678])\tstd: tensor([1.4618, 1.7033])\n",
-      "neg. ELBO: 145.56\tmu: tensor([2.7447, 5.4631])\tstd: tensor([1.4646, 1.7079])\n"
+      "neg. ELBO: 313.53\tmu: tensor([3.0438, 6.6150])\tstd: tensor([1.0943, 1.0870])\tphi: tensor([0.5431, 0.4569])\n",
+      "neg. ELBO: 230.57\tmu: tensor([1.7261, 5.3456])\tstd: tensor([1.9472, 1.4770])\tphi: tensor([0.2626, 0.7374])\n",
+      "neg. ELBO: 224.19\tmu: tensor([0.4331, 5.1358])\tstd: tensor([2.2290, 1.5297])\tphi: tensor([0.1794, 0.8206])\n",
+      "neg. ELBO: 222.64\tmu: tensor([-0.5704,  5.0750])\tstd: tensor([2.0617, 1.5510])\tphi: tensor([0.1491, 0.8509])\n",
+      "neg. ELBO: 219.15\tmu: tensor([-1.3539,  5.0196])\tstd: tensor([1.0274, 1.5799])\tphi: tensor([0.1279, 0.8721])\n",
+      "neg. ELBO: 219.14\tmu: tensor([-1.3579,  5.0186])\tstd: tensor([1.0029, 1.5811])\tphi: tensor([0.1277, 0.8723])\n",
+      "neg. ELBO: 219.14\tmu: tensor([-1.3579,  5.0186])\tstd: tensor([1.0029, 1.5811])\tphi: tensor([0.1277, 0.8723])\n",
+      "neg. ELBO: 219.14\tmu: tensor([-1.3579,  5.0186])\tstd: tensor([1.0029, 1.5811])\tphi: tensor([0.1277, 0.8723])\n"
      ]
     }
    ],
    "source": [
     "lr = 0.01\n",
-    "for i in range(10):𝑝(𝑧𝑖)\n",
+    "for i in range(150):\n",
     "    \n",
     "    # expectation step\n",
     "    with torch.no_grad():\n",
-    "        qz = expectation(x_i, mu, std, k)\n",
+    "        qz = expectation(x_i, mu, std, phi, k)\n",
     "    \n",
     "    # maximization step\n",
-    "    for j in range(150):\n",
+    "    for j in range(100):\n",
     "        negative_elbo = -get_elbo(mu, std, phi)\n",
     "        negative_elbo.backward()\n",
     "        \n",
-    "        with torch.no_grad():\n",
-    "            mu -= mu.grad * lr \n",
-    "            std -= std.grad * lr\n",
-    "            std.grad.zero_()\n",
-    "            mu.grad.zero_()\n",
-    "            phi.grad.zero_()\n",
+    "        optim.step()\n",
+    "        optim.zero_grad()\n",
     "    \n",
+    "    if i % 20 == 0:\n",
     "        print(\"neg. ELBO: {:.2f}\"\n",
     "              \"\\tmu: {}\"\n",
     "              \"\\tstd: {}\"\n",
-    "         \"\".format(negative_elbo.item(), mu.data, std.data))\n"
+    "              \"\\tphi: {}\"\n",
+    "             \"\".format(negative_elbo.item(), mu.data, std.data, F.softmax(phi, dim=-1).data))\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7f7f3402bf28>]"
+       "[<matplotlib.lines.Line2D at 0x7f3a8b76cb38>]"
       ]
      },
-     "execution_count": 27,
+     "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFlCAYAAADmqMVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeZzO5f7H8dc1M8aIscQgxp49+9gisowtJSRbVDpFUTmd+knldFJJnTpp0SLFsZRUlN2MQgkxdlmHLMNgxpYZY9br98etjmowY+77/t4z834+Hh7MfX+/1+c95WHmM9f1vS5jrUVERERERETEaX5OBxAREREREREBNagiIiIiIiLiI9SgioiIiIiIiE9QgyoiIiIiIiI+QQ2qiIiIiIiI+AQ1qCIiIiIiIuITApwOkJlSpUrZypUrOx1DRERERERE3GzDhg3x1tqQzN7zyQa1cuXKREVFOR1DRERERERE3MwYc/By72mJr4iIiIiIiPgENagiIiIiIiLiE9SgioiIiIiIiE9QgyoiIiIiIiI+QQ2qiIiIiIiI+AQ1qCIiIiIiIuIT1KCKiIiIiIiIT1CDKiIiIiIiIj4hSw2qMaaLMWa3MSbaGPP0Fa5raoxJN8bcld17RUREREREJH+7aoNqjPEHJgJdgTpAf2NMnctc9yqwNLv3ioiIiIiIiGRlBrUZEG2t3W+tTQFmAT0yue5R4CvgxDXcKyIiIiIiIvlcVhrU8sDhSz6Oufja74wx5YGewAfZvfeSMR4yxkQZY6Li4uKyEEtERERERETykqw0qCaT1+yfPp4AjLLWpl/Dva4XrZ1krQ2z1oaFhIRkIZaIiPiqCxcu8Oabb5KRkcELL7xAvXr16Ny5M0OGDGHMmDHs3LnT6YgiIiLigwKycE0MUOGSj0OBo3+6JgyYZYwBKAV0M8akZfFeERHJI6y1HD58mNDQUBITE7lw4QKPP/44PXv25MiRIxw5coTo6GjatGnDlClT6N69u9ORRURExIcYazOd0PzfBcYEAHuADsARYD0wwFr782WunwossNZ+md17fxMWFmajoqKy+amIiIiTtm/fzvDhwwkODmbBggVXvHb9+vWkp6cTFhaGn58ffn469UxERCS/MMZssNaGZfbeVb8jsNamASNw7c67E5htrf3ZGDPMGDPsWu7N7icgIiK+beLEibRv356+ffvyzTffXPX6pk2b0qJFC9577z1eeeUVLyQUERGR3OCqM6hO0AyqiEjusXz5cgYMGMDatWupVKlStu5NSUnhwoULnD17lvPnz1OzZk0PpRQRERFfkaMZVBERkcs5cuQIAwcOZPr06dluTgECAwMpWrQoUVFRtGvXDu3iLiIikr+pQRURkWtWtGhR3nvvPTp27JijcXr27MnAgQN5+OGH8cWVPSIiIuIdalBFROSa/Pvf/+bUqVPceeedbhnvxRdfZMeOHcyaNcst44mIiEjuowZVRESuSZkyZShevLjbxgsKCmLatGmMHDmS2NhYt40rIiIiuYcaVBERyZatW7cye/ZsBg8eTLFixdw6dlhYGEOHDmX06NFuHVdERERyBzWoIiKSZcnJydx1112kpKR4rMZzzz3HG2+84bHxRURExHepQRURkSybNGkSNWrU4J577vFYjcDAQEqUKEHPnj05evSox+qIiIiI7wlwOoCIiOQOiYmJjBs3jsWLF3u8lp+fH8888wxly5b1eC0RERHxHZpBFRGRLHn77bdp06YNDRs29Eq9pk2bEhkZycKFC71ST0RERJynBlVERK7q3Llz/Oc//2Hs2LFerWutZdSoUWRkZHi1roiIiDhDDaqIiFxVkSJFWLlyJTVr1vRq3c6dOxMUFMS8efO8WldEREScoQZVRESu6MSJE7zzzjvUqVPH67WNMTzzzDOMGzcOa63X64uIiIh3qUEVEZErSk5OJjg42LH6d955JwkJCSxbtsyxDCIiIuIdalBFROSyjh8/TmBgIPfff79jGfz8/Bg9ejQvv/yyYxlERETEO3TMjIiIXNaYMWMICQlxvDns378/pUuXxlqLMcbRLCIiIuI5mkEVEZFM7d27lzlz5vCPf/zD6SgEBARw6623Mn/+fKejiIiIiAepQRURkUxNmDCBhx9+mOuvv97pKIDryJk5c+aQkpLidBQRERHxEC3xFRGRv0hISOCzzz5j69atTkf5XVBQEFOnTiUpKcnpKCIiIuIhmkEVEZG/mDVrFrfccguhoaFOR/mD9PR06tWrx+7du52OIiIiIh6gBlVERP4iPj6e4cOHOx3jL/z9/Rk8eDCvvvqq01FERETEA4wvHnweFhZmo6KinI4hIpIv+fpOuXFxcVSvXp2DBw9SrFgxp+OIiIhINhljNlhrwzJ7TzOoIiLyByNHjuSrr75yOsZlhYSE0L59e2bPnu10FBEREXEzNagiIvIHY8eOpWPHjk7HuKIhQ4bwySefOB1DRERE3EwNqoiI/C4yMpJffvnF55fOdunShZiYGA4fPux0FBEREXEjNagiIgK4nj196qmniIuLczrKVQUEBLBnzx4qVKjgdBQRERFxIzWoIiICwLp16zh37hwdOnRwOkqWFChQgCeffJK0tDSno4iIiIibqEEVEREAPvjgA4YOHYqfX+740hAQEEDt2rVJTk52OoqIiIi4Se74LkRERDzqzJkzzJ07l/vuu8/pKNnywAMP5IolySIiIpI1WWpQjTFdjDG7jTHRxpinM3m/hzFmqzFmszEmyhjT+pL3Dhhjtv32njvDi4iIe0yfPp0uXbpQunRpp6NkS0JCAg0bNlSTKiIikkdctUE1xvgDE4GuQB2gvzGmzp8u+xZoYK1tCAwBJv/p/XbW2oaXO4xVREScFRISwt///nenY2RbkSJF6NGjBzNmzHA6ioiIiLhBVmZQmwHR1tr91toUYBbQ49ILrLUJ1lp78cPCgEVERHKFCxcu0LdvX5o3b+50lGty//33M2XKFP73ZUhERERyq6w0qOWBSw+ai7n42h8YY3oaY3YBC3HNov7GAhHGmA3GmIdyElZERNxv4sSJjB071ukY16xNmzYkJiayYcMGp6OIiIhIDgVk4RqTyWt/+TG1tXYuMNcY0wZ4Eeh48a1W1tqjxpjSQKQxZpe19vu/FHE1rw8BVKxYMav5RUQkh5544olcvROun58f9913H1OmTCEsTE+SiIiI5GZZmUGNAS49CT0UOHq5iy82n9WMMaUufnz04u8ngLm4lgxndt8ka22YtTYsJCQki/FFRCQndu7cyccff0xQUJDTUXLk3nvv5ccff9QyXxERkVwuKw3qeqC6MaaKMSYQ6AfMu/QCY8yNxhhz8c+NgUDgpDGmsDEm+OLrhYFOwHZ3fgIiInLtpk2bxt69e52OkWMVK1Zk06ZNXPxSJCIiIrnUVRtUa20aMAJYCuwEZltrfzbGDDPGDLt4WW9guzFmM64df/te3DSpDLDKGLMFWAcstNYu8cQnIiIi2WOt5fPPP6dv375OR3GLn3/+mVdeecXpGCIiIpIDWXkGFWvtImDRn1774JI/vwq8msl9+4EGOcwoIiIesGHDBvz9/WnUqJHTUdyiTJkyeeZzERERya+yssRXRETyoN9mT/PKstiQkBBatWrF0aOX3SZBREREfJwaVBGRfMhay+zZs/PM8t7fTJkyhWeffdbpGCIiInKN1KCKiORDCQkJDBgwgJtuusnpKG7Vq1cv5s2bR2pqqtNRRERE5BqoQRURyYeCgoJ45ZVX8szy3t+EhoZSvXp1li9f7nQUERERuQZqUEVE8hlrLQ0bNuTIkSNOR/GI3r1789VXXzkdQ0RERK6BGlQRkXzGGMNPP/1E+fLlnY7iEb1792bu3Lmkp6c7HUVERESySQ2qiEg+M2PGjDy9023VqlWpXLkyW7dudTqKiIiIZJMaVBGRfCQ9PZ0nn3wSP7+8/c//jz/+qDNRRUREcqG8/R2KiIj8wQ8//ED58uW58cYbnY7iUdZaRo0aRUZGhtNRREREJBvUoIqI5CPffPMNPXv2dDqGxwUGBlK5cmWSk5OdjiIiIiLZoAZVRCSfsNYyb948br/9dqejeMWwYcNISEhwOoaIiIhkgxpUEZF8YufOnaSmplK/fn2no3jF7t27CQsLw1rrdBQRERHJIjWoIiL5xK+//soTTzyBMcbpKF5Rs2ZNChUqRFRUlNNRREREJIsCnA4gIiLe0aJFC1q0aOF0DK8xxtC7d2+++uormjZt6nQcERERyQLNoIqI5APx8fE0adIk3y137dGjB/Pnz3c6hoiIiGSRGlQRkXygZMmSfPXVV/lmee9vwsLCiI+P55dffnE6ioiIiGSBGlQRkXxg/vz5FCtWzOkYXufn50fXrl1ZsmSJ01FEREQkC4wvLvcKCwuz2tRCRMQ9kpOTKVOmDHv27KF06dJOx/G6M2fOEBwcjL+/v9NRREREBDDGbLDWhmX2nmZQRUTyuJUrV1KnTp182ZwCFC9enKlTp5KYmOh0FBEREbkKNagiInnc/Pnzuf32252O4ajjx49z9uxZp2OIiIjIVeiYGRGRPMxay/z581mwYIHTURz1zDPPkJGR4XQMERERuQrNoIqI5GHbt2/HGEPdunWdjuKoI0eOULt27Xx3zI6IiEhuowZVRCQP279/P4MGDcp3x8v8Wbly5UhLS2Pr1q1ORxEREZErUIMqIpKH9ejRg7Fjxzodw3HGGG677TYWLlzodBQRERG5AjWoIiJ51LFjxxg8eLDTMXzGbbfdlu+fxRUREfF1alBFRPKo4OBgHnzwQadj+Iy2bdvy888/Ex8f73QUERERuQw1qCIiedSePXto3ry50zF8RlBQEO3atWPJkiVORxEREZHL0DEzIiJ5UGpqKrfeeivR0dGEhIQ4HcdnvPrqq5QoUcLpGCIiInIZWZpBNcZ0McbsNsZEG2OezuT9HsaYrcaYzcaLine truncated
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFlCAYAAADmqMVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3iV9f3/8ecnJwOSMLPYK4FAUGREJOy9tALFPVqrdW9UtEOroq36q9b6tUqpVq2iFgdDg4QZQkREtuywNyRhhATIOp/fHwkWFSXc505ykrwe18Vlcu77fp9XiCHnfe7PMNZaRERERERERCpbQGUHEBEREREREQE1qCIiIiIiIuIn1KCKiIiIiIiIX1CDKiIiIiIiIn5BDaqIiIiIiIj4BTWoIiIiIiIi4hcCKzvA2URGRtpWrVpVdgwRERERERFx2fLly7OstVFnO+aXDWqrVq1YtmxZZccQERERERERlxljdv7UMQ3xFREREREREb+gBlVERERERET8ghpUERERERER8QtqUEVERERERMQvqEEVERERERERv6AGVURERERERPyCGlQRERERERHxC2pQRURERERExC+UqUE1xgw3xmwyxmwxxjz2M+ddbIwpNsZccb7XioiIiIiISM12zgbVGOMB/gGMABKAa40xCT9x3vNAyvleKyIiIiIiIlKWO6jdgS3W2m3W2gLgQ2DUWc67F/gEOOTgWhEREREREanhytKgNgV2n/H5ntLHvmOMaQqMASae77UiUjVYays7goiIiIhUc2VpUM1ZHvvhK9WXgUettcUOri050ZjbjDHLjDHLMjMzyxBLRMrCWsvGAzlMXbmHzOP55339+vXrefDBB4mMjOS2224DYNu2bXi9XrejioiIiEgNF1iGc/YAzc/4vBmw7wfnJAIfGmMAIoGRxpiiMl4LgLV2EjAJIDExUbdqRHxwqrCYr7ZmM3/jIeZvPMTeoycBCPIYhnZsxPXdW5AUG0Hpz+yPnDhxgilTpvCvf/2L7du385vf/IalS5fSsGFDAG6++WbefvttiouLOXToEElJSRX2tYmIiIhI9VWWBvUboK0xpjWwF7gGuO7ME6y1rU9/bIx5G/jcWjvNGBN4rmtFxD1er+XZmRuY/PVOThV6qR3koXfbSO4ZGEf7RnX4fM1+Pl6+h+Q1+2kTGca13Vtw1cXNqVc76Ht1vvjiCz7++GPGjx/PpZdeSmDg9/+pSE1NBSA9PZ2xY8cyceJExowZU1FfpoiIiIhUU+dsUK21RcaYeyhZndcD/Ntau84Yc0fp8R/OOz3nte5EF5EzFXstj36yho+X72FMl6aM7tKUS1o3pFaQ57tzurRowCPD4pn57X4mf72LZ2du4L/LdvPR7Uk0CAtm5cqVLF68mLvvvpuxY8ee8zl79+7NF198wYgRIzDGMHr06PL8EkVERESkmjP+uPBJYmKiXbZsWWXHEKkyir2WRz5azacr93L/oLY8MLjtTw7fPdOijExueWcZCY3rMvm3l3D40H7Wrl3LiBEjzuv5V6xYwYgRI5g0aRKjRmmhbhERERH5acaY5dbaxLMdK8siSSLix4qKvYybsopPV+7loSHteHBIuzI1pwB92kbxf9d2YcWGLVw0+JdERjc67+YUoGvXrsycOZPbbruN5OTk875eRERERATUoIpUaYXFXu7/7yqmr9rH+OHx3Duo7XnXuLhxMHbms2QHRjD+07V4vc5GVXTr1o3k5GRat2597pNFRERERM5CDapIFVVU7OW+D1aSvGY/vx/Znrv6x513jdzcXEaOHMm1V4zmz3/6I5+t3seTn61zvOdpYmIiCQkJ3HPPPWiYvoiIiIicr7Ks4isifug/X+3ki7UH+OOlHfhtnzaOatxxxx106NCB559/HoAjJwqYlLaNhmHBPDC4neNst9xyC+3bt3d8vYiIiIjUTGpQRaqggzmneGnOZvq1i+KW3s6G1Kanp7Nw4UI2btz43ZzV341oz+G8Al6em8ElrSNIio1wVLtLly6sWLGC7du3l2k1YBERERER0BBfkSppwufrKSj28vSojmVeEOmHvF4vr7/+OmFhYd89ZozhmdEX0LR+bZ76bB3FDuejAgQGBnLnnXeSlZXluIaIiIiI1CxqUEWqmEUZmXy+Zj9394+jZUTYuS84i3Xr1pGUlMRll132o2O1gjz8fmQHNh44zoff7HKcs1OnTlx//fX87ne/c1xDRERERGoWNagiVcipwmIen7aW1pFh3N7P2bxTgL/+9a+sWbPmJ4+PvLAR3Vs35MXZmzl2stDx8zz55JPMnDmTJUuWOK4hIiIiIjWHGlSRKmRS2jZ2ZJ/g6VEdqRXkcVQjPz+ft956i27duv3kOcYY/vSLBI6cKOCVeRlO41KvXj1eeOEF7r77boqLix3XEREREZGaQQ2qSBWxMzuPVxds4bJOjenTNspRja1bt3LBBRdQUFBwznM7NqnHNRc3553FO9hyKNfR8wFcd911hIeH8/777zuuISIiIiI1gxpUkSrAWssT09cR7Ang8csSHNd5+OGHufnmmwkODi7T+Q8Njad2kIdnktc7fk5jDFOmTOHaa691XENEREREagY1qCJVQMq6AyzcnMm4Ie2IqVvLUY158+axevVqHnzwwTJfExkewv2D25K6KZMFGw85el6AmJgYNmzYwNNPP+24hoiIiIhUf2pQRfyctZYXZ28mPqYOv0pq6ahGcXExDzzwAC+++CK1ap1fg/urpFa0iQwr2dqmyOvo+QEaN25Mz549HV8vIiIiItWfGlQRP5e6OZOMQ7nc3q8NgR5nP7LJycnUrl2b0aNHn/e1wYEB/PGyDmzLymPy1zsdPT9AZGQkPXr0YOnSpY5riIiIiEj1pgZVxM+9sWgbMXVDuKxTE8c11q5dy7hx4zDGOLp+QHw03Vs15M307RR7reMce/bs4Re/+AWnTp1yXENEREREqi81qCJ+bN2+Y3y5JZuberYmOND5j+vvf/97rrnmGsfXG2O4uXcr9hw5yZz1Bx3Xad++PYmJiUyePNlxDRERERGpvtSgivixNxdtJzTYw3XdWziuMW7cOD7//HOfswxJaETT+rV568vtPtV58MEH+dvf/oa1zu/EioiIiEj1pAZVxE8dOHaKGav3cVVic+qFBjmu8+ijj9KrVy+f83gCDL/u2ZKvtx9m3b5jjusMGjSIgIAA5syZ43MmEREREale1KCK+Km3F+/Aay239G7tuMbMmTPJycmhQYMGrmS6OrEFtYM8vLN4h+MaxhgefPBBXnrpJVcyiYiIiEj1oQZVxA/l5Rfx/tc7GX5BI5o3DHVUo6ioiDvvvJOcnBzXctULDeKXXZsybdU+snPzHde57rrriI6OpqCgwLVsIiIiIlL1qUEV8UNTlu0m51QRv+3TxnGN6dOn07x5c7p16+ZiMvhNr1YUFHn5YOkuxzVCQkL4z3/+42IqEREREakO1KCK+Jlir+XfX26nW8sGdG3hfGjuK6+8wn333edishJx0XXo0zaSd5fspLDY67hOfn4+CQkJ5ObmuphORERERKoyNagifiZl3QF2Hz7JrX2czz1dtWoV27ZtY8yYMS4m+5+be7XmYE4+M7/d77hGSEgIK1asIDw83MVkIiIiIlKVqUEV8TP/WrSNlhGhDElo5LjGK6+8wl133UVQkPPVf39Ov3ZRtI4M420fFksCCAsLY9y4ceTnO5/PKiIiIiLVhxpUET+ybt8xVu46yk09W+EJMI7rdOzYkVtvvdXFZN8XEGD4dVJLVu46yqrdRx3X8Xg8rFy5kuTkZBfTiYiIiEhVpQZVxI98umIvQR7D6M5NHdc4fPgw48aNIzIy0sVkP3ZFYnPqhATy1pfbfapz44038u6777qUSkRERESqMjWoIn6iqNjL9FX7GNg+mgZhwY7rXHfddXzzzTcuJju78JBArkxsTvKa/RzOc75dzNixY5k/fz6HDx92MZ2IiIiIVEVqUEX8xKItWWTl5vPLrs18qjNz5kwuvvhil1L9vCsTm1HktSSv2ee4Rr169Rg+fDhTpkxxMZmIiIiIVEVlalCNMcONMZuMMVuMMY+d5fgoY8waY8wqY8wyY0zvM47tMMZ8e/qYm+FFqpOpK/ZSPzSIAfHRjmtMnDiRr7/+GmOcz189Hx0a1yU+pg7TVjlvUKFkmO97773nUioRERERqarO2aAaYzzAP4ARQAJwrTEm4QenzQMustZ2Bm4G3vjB8QHW2s7W2kQXMotUO8dPFZKy7gC/6NSE4EBnAxu8Xi/PPvssdevWdTndzxvVpQnLdx5h9+ETjmsMGzaMyy67DK/X+b6qIiIiIlL1leWVcHdgi7V2m7W2APgQGHXmCdbaXGutLf00DLCISJl98e0B8ou8/LKr88WR0tLSiIiIoGPHji4mO7fLL2oCwPRVex3XCAoK4rHHHuPYsWNuxRIRERGRKqgsDWpTYPcZn+8pfex7jDFjjDEbgWRK7qKeZoHZxpjlxpjbfAkrUl19unIPrSPD6Ny8vuMakydP5vrrr3cxVdk0axBK99YNmbpyL/97n+r8rV69mjFjxriYTERERESqmrI0qGebzPajV6HW2qnW2vbAaGDCGYd6WWu7UjJE+G5jTN+zPokxt5XOX12WmZlZhlgi1cOeIydYsu0wv+zS1PHc0fz8fD799FOuueYal9OVzejOTdmamce6fTmOa3Tq1Il58+a5mEpEREREqpqyNKh7gOZnfN4M+MkVUay1aUCsMSay9PN9pf89BEylZMjw2a6bZK1NtNYmRkVFlTG+SNU3bWXJ0NjRXZwP7505cyYXXnghzZs3P/fJ5WDkhY0I8pjvvhYnjDGsWrWKV155xcVkIiIiIlKVlKVB/QZoa4xpbYwJBq4BZpx5gjEmzpTe+jHGdAWCgWxjTJgxpk7p42HAUGCtm1+ASFVmreXTlXu5pHVDmjcLine truncated
       "text/plain": [
        "<Figure size 1152x432 with 1 Axes>"
       ]
@@ -464,22 +488,22 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7f7f344ad438>]"
+       "[<matplotlib.lines.Line2D at 0x7f9520175f60>]"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6IAAAFlCAYAAADxilWiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3iV9f3/8eednZCQxYYsAjJERghhKihi2QkqAk6sWq04UNviV2y12tqfrVqLLSqiIi5wYMABUpbIJmxQIBCSMEMgIWSvc//+OEBRGUnOfZJzktfjuriUc+77fb9QriTv81mGaZqIiIiIiIiI1BaPug4gIiIiIiIiDYsaUREREREREalVakRFRERERESkVqkRFRERERERkVqlRlRERERERERqlRpRERERERERqVVedfXgJk2amNHR0XX1eBEREREREXGiTZs2nTBNs+mF3quzRjQ6OpqUlJS6eryIiIiIiIg4kWEYGRd7T1NzRUREREREpFapERUREREREZFapUZUREREREREapUaUREREREREalVakRFRERERESkVqkRFRERERERkVqlRlRERERERERqlRpRERERERERqVVqREVERERERKRWqREVERERERGRWqVGVERERERERGqVGlERkQZuz7F8cgvL6jqGiIiINCBedR1ARETqxr7jBfz16x9Yvicbb0+DQR2acVNca67t2AxfL8+6jiciIiL1mBpREZEG5lRRGa8uSeX9dRkEeHvy+1914FRRGclbj/DfH7II9vdmRNeW3NY7kitbBdd1XBEREamH1IiKiDQQFZU2PliXwT+XpJJfUs74hEgeH3IFTQJ9AZgytCOr95/ki82H+GLzYT7bdIh5v+1Hl9ZqRkVERMRahmmadfLg+Ph4MyUlpU6eLSLSEP0xeSfvr8ugf7twnh7RmU4tG1/02uz8Ukb/exWeHgZfPTyAkACfWkwqIiIi9YFhGJtM04y/0HvarEhEpAFYuTeb99dlcHf/aD64p/clm1CApkG+vH57T46fLuXROVux2ermQ0sRERGpn9SIiojUc3nF5Uz5fDuxTRsxZWhHDMOo0n3dI0J4ZnRnvtubzb+Wpjo5pYiIiDQkakRFROq557/6gazTJbx8S3f8vKu3G+6tCZHc3LMN/1qayrLdWU5KKCIiIg2NGlERkXpsyQ9ZfLbpEA8Oakf3iJBq328YBn9J6kLnlo2ZPGcrmSeLnJBSREREGho1oiIi9VRuYRlPzttBxxZBPDK4fY3r+Hl78sbtPTEMg/s/2ERxWaWFKUVERKQhUiMqIlJP/WnBLvKKy3jllu74eDn25T4yPIBXx3Xnx6OneWf1AYsSioiISEOlRlREpB76evtRvtx2hEeua0/nVpfeIbeqru3YjEEdmvLOqgOUlGtUVERERGpOjaiISD1TUl7JMwt20rVNML8dFGtp7d8OjOVkYRmfphy0tK6IiIg0LGpERUTqmQVbj3CioIwnh3bEy9PaL/MJMWHERYbw5so0KiptltYWERGRhkONqIhIPWKaJu+sPkCH5kH0jQ23vL5hGPx2UDsO5Rbz1fajltcXERGRhkGNqIhIPbIuLYfdx/K5u380hmE45RmDOzbjiuaBvL5iP6ZpOuUZIiIiUr+pERURqUfeXX2A0ABvknq0dtozPDwMHhgYy56sfJbvOe6054iIiEj9pUZURKSeOJhTxH9/zGJCQiR+3p5Ofdaobq1oHeLP6yv2O/U5IiIiUj+pERURqSfeW5OOh2FwR98opz/L29OD+66OYWN6LhvTc5z+PBEREalf1IiKiNQDhaUVzE05yLAuLWgZ7F8rzxzXK5KwRj68oVFRERERqSY1oiIi9cDnmw+RX1LB3f1jau2Z/j6eTOwXzdLdx9l97HStPVdERETcnxpRERE3Z7OZzFqdTrc2wcRFhtTqs+/sG0WAj6dGRUVERKRa1IiKiLi571KzSTtRyN39Y5x2ZMvFhAT4cHPPNnyz8xj5JeW1+mwRERFxX2pERUTc3Lur02kW5Mvwq1rWyfMTu7eirMLG4l1ZdfJ8ERERcT9qREVE3Ni+4wWs3JvN7X2i8PGqmy/pcZGhtA7xZ8G2I3XyfBEREXE/akRFRNzYvM2H8PQwmJAQWWcZDMNgVLdWrNp3gpMFpXWWQ0RERNyHGlERETdlmibf7DhK37bhNA3yrdMso7u1otJm8s3OY3WaQ0RERNyDGlERETf149F80k8W1dna0PN1ahlEu2aBfLlV03NFRETk8tSIioi4qW92HMXDgF9d2byuo2AYBqO7tWJDeg5HThXXdRwRERFxcWpERUTc0NlpuX3ahhMeWLfTcs8a3a0VAF9t16ioiIiIXJoaURERN7QnK5+0E4UuMS33rOgmjejaJli754qIiMhlqREVEXFD32y3T8sd2qVFXUf5idHdWrHz8GnSsgvqOoqIiIi4MDWiIiJuxjRNvt5xlN4x4TRxkWm5Z43s2grDQKOiIiIicklqREVE3Ezq8QL2Zxcy/CrXGg0FaBHsR0J0GAu2HcE0zbqOIyIiIi5KjaiIiJv5evtRDAN+5WLTcs8a3b0VadmF7Dpyuq6jiIiIiItSIyoi4ma+2XGUhOgwmgX51XWUCxrepSVeHgZfanquiIiIXESVGlHDMIYahrHHMIx9hmE8eYnrehmGUWkYxs3WRRQRkbNSs/JJPV7gUrvl/lxoIx+ubt+EL7cdwWbT9FwRERH5pcs2ooZheAL/AYYBnYEJhmF0vsh1LwLfWh1SRETsvtlxDMOAYS46Lfes0d1bcSSvhC0Hc+s6ioiIiLigqoyIJgD7TNNMM02zDJgDJF7guoeBz4HjFuYTEZHzfLPjKL2iwmjW2DWn5Z51XcfmeHoYLNutbwkiIiLyS1VpRFsDB8/7/aEzr51jGEZrYAzwhnXRRETkfPuOF7AnK98ld8v9uWB/b3pGhrJ8d3ZdRxEREREXVJVG1LjAaz9f9PMqMMU0zcpLFjKM3xiGkWIYRkp2tn44ERGpjoU7jgIwtIvrrA/dsGED7733HgDjx4/nxx9/JCMjg8mTJ2Pb+RU7M4+TdbqkjlOKiIiIq6lKI3oIiDjv922An2+FGA/MMQwjHbgZmG4YRtLPC5mmOcM0zXjTNOObNm1aw8giIg3Tol3H6BkVSovgup+Wu2PHDpKSkrjpppvOvfbss88SFRWFj48PkZGRnD6wk6OzHuXtL5bUYVIRERFxRVVpRDcC7Q3DiDEMwwcYDyw4/wLTNGNM04w2TTMa+Ax40DTNZMvTiog0UMfzS9h15DSDOzWr0xz79u3jtttuY8iQIQwcOJDU1FTuuusuADp27EhAQAAtW7bk8ccf59uvkokaMpHnH7mT559/noqKijrNLiIiIq7jso2oaZoVwEPYd8P9EfjENM1dhmE8YBjGA84OKCIi8P3eEwBc077uZpPk5eUxcOBAOnXqRGpqKo899hh+fhcfnTUMg5tuvoWoe6axYsV3zJo1q/bCioiIiEvzqspFpml+A3zzs9cuuDGRaZoTHY8lIiLnW5maTZNAHzq3bFwnzz99+jTBwcHs3LmT0NDQKt93bcemzE0J5S9vfkxCTBhLliyhuLiYUaNGOTGtiIiIuLqqTM0VEZE6ZLOZfJ96gqvbN8XD40L7xznf2LFj2bx5c7WaUID+7Zrg5WGwct9JPD09ady4MWFhYU5KKSIiIu5CjaiIiIvbeSSPnMIyBl5R+9Ny8/PzKSoqIjk5mbi4uGrfH+TnTXx0KMvPnCeakJBA3759+fvf/05+fr7VcUVERMRNqBEVEXFxK/faj7sa0L5JrT63rKyMMWPG8Prrr+Pv71/jOoM6NGP3sXyO5dmPcfHw8CA1NZV77rkH0/z5aWAiIiLSEKgRFRFxcSv3nqBL68Y0CfSttWfabDYmTpxIUFAQkydPdqjWtR3sO/2u2HP83GuvvfYaaWlpvPrqqw7VFhEREfekRlRExIWdLilnc2Zure+W++qrr5Kens5HH32Ep6enQ7WuaB5Iy2A/VuzJPvean58fn3/+OS+++CIrV650NK6IiIi4GTWiIiIubM2+k1TYTK6pxfWhR44c4YUXXmDWrFkOTck9yzAMBnVoxqp9JyirsJ17PSoqitmzZzNhwgSOHDni8HNERETEfagRFRFxYStTs2nk40lcZPV2q3XEH/7wB+6//36uuOIKy2oO6tCUgtIKNmXk/uT1G264gYcfflijoiIiIg1Mlc4RFRGR2meaJiv3ZtOvXRN8vGrvc8OnnnqKqKgoS2v2b9cEb0+DFXuP0zc2/CfvPfnkkwDk5eURHBxs6XNFRETENWlEVETERR04Ucih3OJam5ZbXl7OM888Q2xsLI0aNbK0dqCvF72iw1ixO/uC7//www+MGjXK0meKiIiI61IjKiLior47c2zLwFraqKisrIxmzZrh4+PjlPrXdmjGnqx8jpwq/sV7nTt3ZsmSJU55roiIiLgeNaIiIi5q5d5sosMDiAwPcPqzjhw5wv79+5k0aRKGYTjlGYM62Bvq83fPPZ9hGAwbNoycnBynPF9ERERchxpREREXVFpRybq0nFqblvv444/zySefOPUZ7ZoF0jrEn+/2Hr/g+97e3kRHR/OXv/zFqTlERESk7qkRFRFxQSnpuRSXV9bK+aFLly5l/fr1PPXUU059jmEYDGjXhHVpOdhs5gWvefbZZ5k9ezZpaWlOzSIiIiJ1S42oiIgLWrk3G29P4xc7zFrNZrPx2GOP8corrxAQ4PwpwH1jw8krLueHo6cv+H7z5s159NFHnd4Ui4iISN1SIyoi4oK+25tNfFQYjXyde8rWwoUL8fT0JCkpyanPOetsY712/8mLXvP444/z/fffs379+lrJJCIiIrVPjaiIiIvJOl3C7mP5tbI+9MUXX+QPf/iD0zYo+rnmjf1o26QRa9Mu3og2atSI5557junTp9dKJhEREalLine truncated
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6gAAAFlCAYAAADmqMVrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZfrG8e+bSSGETkgA6aH3EpDQIUAEC/JDRSy7rqvYGyqyRVdFdtW1sq6i667uKsiiSNHQQgcRaQKGThJqkN4Jqe/vjwkuKkg4c5KZJPfnuriSmTnnmTvJRTLPnLcYay0iIiIiIiIi/hbk7wAiIiIiIiIioAZVREREREREAoQaVBEREREREQkIalBFREREREQkIKhBFRERERERkYCgBlVEREREREQCQrC/A1xIZGSkrVevnr9jiIiIiIiIiMtWr159yFpb7UKPBWSDWq9ePVatWuXvGCIiIiIiIuIyY8zOiz2mIb4iIiIiIiISENSgioiIiIiISEBQgyoiIiIiIiIBQQ2qiIiIiIiIBAQ1qCIiIiIiIhIQ1KCKiIiIiIhIQFCDKiIiIiIiIgFBDaqIiIiIiIgEhAI1qMaYq4wxW4wx240xo37huI7GmFxjzA2Xe66IiIiIiIiUbpdsUI0xHuDvwACgOTDMGNP8Ise9BMy+3HNFRERERERECnIFtROw3Vqbaq3NAiYCgy5w3EPAZOCAg3NFRERERESklCtIg3oFsPu823vy7/uBMeYKYDAw7nLPFZHiwVrr7wgiIiIiUsIVpEE1F7jvp69U3wCestbmOjjXe6Axw40xq4wxqw4ePFiAWCJSINbC/g2w7r9w6sClj/+JjRs38thjjxEZGcnw4cMBSE1NJS8vz+2kIiIiIlLKBRfgmD1A7fNu1wLSf3JMLDDRGAMQCQw0xuQU8FwArLXvAe8BxMbG6lKNiC+yMyBtCWybDVtnw/H8gQxBIdDsGujwG6jfA8yF3kOCM2fOMGnSJP7xj3+QlpbGb37zG1asWEGVKlUAuPPOO/nwww/Jzc3lwIEDxMXFFdVXJiIiIiIlWEEa1JVAI2NMfWAvcDNwy/kHWGvrn/vcGPMh8KW1dqoxJvhS54qIi/LyIOlpWPlPyMmAkLLQoDf0eAKiW0Ly57B2PGyYAlUbQoc7oN3tEF7pR2VmzpzJZ599xsiRI7n66qsJDv7xr4qFCxcCsHTpUoYMGcK4ceMYPHhwEX2RIiIiIlJSXbJBtdbmGGMexLs6rwf4l7V2gzHm3vzHfzrv9JLnuhNdRH4kLxemPwxrP4bWQ6H1TVC3G4SU+d8xtWIh/mnYOA1W/Qvm/BHWfAR3zoKyVfj2229ZtmwZDzzwAEOGDLnkU3br1o2ZM2cyYMAAjDFcf/31hfgFioiIiEhJZwJx4ZPY2Fi7atUqf8cQKT7ycmHq/bB+IvQcBb1GXXT47o+kzIcJN0P1VvCraew+cJTk5GQGDBhwWU+/Zs0aBgwYwHvvvcegQVqoW0REREQuzhiz2lobe6HHCrJIkogEstwcmHKPtznt/Ufo/buCNacAMX3ghn+xd/Mqftu/BTWjIi+7OQVo3749M2bMYPjw4SQmJl72+SIiIiIioAZVpHjLzYbP74LvPoX4P0HPJy+7xJHoLvSfWo4mofvxfPGAdx6rAx06dCAxMZH69etf+mARERERkQtQgypSXOXmwGd3ehc86jcauo+47BKnTp1i4MCBXD3kVkb+6S+QPBlmjvRuTeNAbGwszZs358EHH0TD9EVERETkchVkFV8RCUQr/wGbpkP/MdDlQUcl7r33Xpo1a8ZLL73kveP0Ifj6LYiI9M5jdei3v/0tTZs2dXy+iIiIiJROalBFiqMT+2D+GGjYF+IecFRi6dKlLFq0iM2bN2POzVnt/wKcOQIL/wJ1u0L97o5qt2vXjjVr1pCWllag1YBFREREREBDfEWKp9m/h9wsGPjXgi+I9BN5eXm88847RERE/O9OY+Ca16BibZg1yrs6sEPBwcHcd999HDp0yHENERERESld1KCKFDcp82HD59D9cajSwFGJDRs2EBcXxzXXXPPzB0PCof9o2J8Ma/7tOGbr1q259dZb+d3vfue4hoiIiIiULmpQRYqT7LOQ+DhUiYGujzgu88orr7B+/fqLH9D8eu8Q3/kvQMYxx8/z7LPPMmPGDJYvX+64hoiIiIiUHmpQRYqTZWPhSCpc/QqElHFUIjMzkw8++IAOHTpc/CBj4Kq/eOejLnrZYVioWLEiL7/8Mg888AC5uc6HC4uIiIhI6aAGVaS4OJIKi1+BFv8HMX0clUhJSaFly5ZkZWVd+uAabaD9r2DFu3Bwq6PnA7jlllsoV64cEyZMcFxDREREREoHNagixYG1MONJ8IRCwp8dl3niiSe48847CQ0NLdgJfZ6GkLLeRZkcMsYwadIkhg0b5riGiIiIiJQOalBFioNNX8D2udDnD1ChhqMS8+bNY926dTz22GMFP6lcNej5FGxPgq1zHD0vQHR0NJs2beL55593XENERERESj41qCKBzlpYMAaimkPHux2VyM3N5dFHH+XVV1+lTJnLnLvaaThUbQizfwc5BRgafBE1atSgS5cujs8XERERkZJPDapIoNs+Fw5u9q7a6wl2VCIxMZHw8HCuv/76yz85OH9Y8eHtsOpfjp4fIDIyks6dO7NixQrHNURERESkZFODKhLolv0NytfwLo7kUHJyMiNGjMAY46xAo/5Qpwss/zvkOV+Nd8+ePVx77bWcPXvWcQ0RERERKbnUoIoEsn3rIW0RXHmP90qmQ7///e+5+eabnecwBjrfB8d2wZYZjss0bdqU2NhYxo8f7zyLiIiIiJRYalBFAtnXf4eQCOhwh+MSI0aM4Msvv/Q9S9OroWIdWD7OpzKPPfYYr7/+OtZa3zOJiIiISImiBlUkUJ1Ih+TPoP3tEF7ZcZmnnnqKrl27+p4nyAOd7oadS71Xdh2Kj48nKCiIpKQk3zOJiIiISImiBlUkUH3zLtg879Bah2bMmMGJEyeoXNl5g/sj7W/37ou64l3HJYwxPPbYY7z22mvuZBIRERGREkMNqkggyjwFqz+AZtdC5XqOSuTk5HDfffdx4sQJ93KFV4Y2N8P6T+H0IcdlbrnlFqKiosjKcr5tjYiIiIiUPGpQRQLRtx/D2eMQ95DjEtOmTaN27dp06NDBxWDAlfdCbqa3gXYoLCyM//znPy6GEhEREZGSQA2qSKDJy4Xlb0PtK6F2R8dlxo4dy8MPP+xisHzVmkBMH1j5T8jNdlwmMzOT5s2bc+rUKRfDiYiIiEhxpgZVJNBs+gKO7YS4Bx2XWLt2LampqQwePNjFYOe58j44uQ82TnNcIiwsjDVr1lCuXDkXg4mIiIhIcaYGVSTQfP0WVK7v3dbFobFjx3L//fcTEhLiYrDzNOwLVWLgG9+2nImIiGDEiBFkZma6FExEREREijM1qCKBZN962LPSO88zyOO4TIsWLbj77rtdDPYTQUFw5T3erHtWOy7j8Xj49ttvSUxMdDGciIiIiBRXalBFAsm6iRAUAq1vclziyJEjjBgxgsjISBeDXUDbWyCsAnzzjk9lbr/9dj766COXQomIiIhIcaYGVSRQ5ObAd59C4wQoW8VxmVtuuYWVK1e6GOwiwspDu9tgwxQ4fdhxmSFDhjB//nyOHDniYjgRERERKY7UoIoEitQFcPoAtBnmU5kZM2bQsaPz1X8vS9tbIS8HNnzuuETFihW56qqrmDRpkovBRERERKQ4KlCDaoy5yhizxRiz3Rgz6gKPDzLGrDfGrDXGrDLGdDvvsR3GmO/OPeZmeJESZd1ECK8Mjfo7LjFu3Di++eYbjDEuBvsF1VtCVHPvlV8f3H777Xz88ccuhRIRERGR4uqSDaoxxgP8HRgANAeGGWOa/+SweUAba21b4E7g/Z883tta29ZaG+tCZpGS5+wJ2PwltBwCwaGOSuTl5TFmzBgqVKjgcrhLaHUj7P4Gju5wXCIhIYFrrrmGvLw893KJiIiISLFTkCuonYDt1tpUa20WMBEYdP4B1tpT1lqbfzMCsIhIwW2cBjlnfRreu3jxYqpWrUqLFi1cDFYArW7wfvThKmpISAijRo3i+PHjLoUSERERkeKoIA3qFcDu827vyb/vR4wxg40xm4FEvFdRz7HAHGPMamPMcF/CipRY6//r3Vf0ig6OS4wfP55bb73VxVAFVKkO1O0K6yeBdf7e1Lp16xg8eLCLwURERESkuClIg3qhyWw/exVqrZ1irW0KXA+MPu+hrtba9niHCD9gjOlxwScxZnj+/NVVBw8eLEAskRLi2C7YscR79dTh3NHMzEw+//xzbr75ZpfDFVCrG+HQVti3znGJ1q1bM2/ePBdDiYiIiEhxU5AGdQ9Q+7zbtYD0ix1srV0MxBhjIvNvp+d/PABMwTtk+ELnvWetjbXWxlarVq2A8UVKgPX/9X70Ye/TGTNm0KpVK2rXrn3pgwtD80He/Vt9GOZrjGHt2rWMHTvWxWAiIiIiUpwUpEFdCTQyxtQ3xoQCNwPTzz/AGNPQ5C8baoxpD4QCh40xEcaY8vn3RwD9gWQ3vwCRYs1aWPdfqNsNKtd1XOaKK67g6aefdjHYZSpbxbt/63efQV6u4zKVK1dm9OjRZGdnuxhORERERIqLSzao1toc4EFgNrAJmGSt3WCMudcYc2/+YUOAZGPMWrwr/g7NXzQpGlhqjFkHrAASrbWzCuMLESmW9q6Bw9ugzVDHJTIyMmjRogXx8fEuBnOg1Y1w6ntIW+y4RIMGDWjSpAmzZunXhIiIiEhpVKB9UK21M6y1ja21MdbaMfn3jbPWjsv//CVrbYv8rWTirLVL8+9Ptda2yf/X4ty5IpJv3ScQXMY7RNahxYsXc/fdd7sYyqHGV0FYBZ/3RL3tttv46KOPXAolIiIiIsVJgRpUESkEOVmQPBmLine truncated
       "text/plain": [
        "<Figure size 1152x432 with 1 Axes>"
       ]
@@ -497,6 +521,7 @@
     "        self.mu = None\n",
     "        self.std = np.ones(k)\n",
     "        self.w_ij = None\n",
+    "    [1] https://www.ritchievink.com/blog/2019/05/24/algorithm-breakdown-expectation-maximization/\n",
     "        self.phi = np.ones(k) / k\n",
     "\n",
     "    def expectation_step(self, x):\n",
@@ -535,11 +560,12 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
-   "source": []
+   "source": [
+    "# References\n",
+    "[1] https://www.ritchievink.com/blog/2019/05/24/algorithm-breakdown-expectation-maximization/"
+   ]
   }
  ],
  "metadata": {

trainings/variational_inference_core_team/models.py

@@ -0,0 +1,63 @@
+import torch
+from torch import nn
+import torch.nn.functional as F
+
+
+class AutoEncoder(nn.Module):
+    def __init__(self, input_size=784, z_size=20):
+        super().__init__()
+        hidden_size = int((input_size - z_size) / 2 + z_size)
+        self.encoder = nn.Sequential(
+            nn.Linear(input_size, hidden_size),
+            nn.ReLU(),
+            nn.Linear(hidden_size, z_size)
+        )
+        self.decoder = nn.Sequential(
+            nn.Linear(z_size, hidden_size),
+            nn.ReLU(),
+            nn.Linear(hidden_size, input_size),
+        )
+
+    def forward(self, x):
+        x = x.view(-1, 784)
+        z = self.encoder(x)
+        x = self.decoder(z)
+
+        if self.training:
+            return x
+        else:
+            return F.sigmoid(x)
+
+
+class VAE(nn.Module):
+    def __init__(self, input_size=784, z_size=20):
+        super().__init__()
+        hidden_size = int((input_size - z_size) / 2 + z_size)
+        self.z_size = z_size
+        self.encoder = nn.Sequential(
+            nn.Linear(input_size, hidden_size),
+            nn.ReLU(),
+            nn.Linear(hidden_size, hidden_size),
+            nn.ReLU(),
+            nn.Linear(hidden_size, z_size * 2)
+        )
+        self.decoder = nn.Sequential(
+            nn.Linear(z_size, hidden_size),
+            nn.ReLU(),
+            nn.Linear(hidden_size, input_size),
+            nn.Sigmoid()
+        )
+
+    def reparameterize(self, mu, log_var):
+        std = torch.exp(0.5 * log_var)
+        eps = torch.randn_like(std)  # unit gaussian
+        z = mu + eps * std
+        return z
+
+    def forward(self, x):
+        x = x.view(-1, 784)
+        variational_params = self.encoder(x)
+        mu = variational_params[..., :self.z_size]
+        log_var = variational_params[..., self.z_size:]
+        z = self.reparameterize(mu, log_var)
+        return self.decoder(z), z, mu, log_var

trainings/variational_inference_core_team/requirements.txt

@@ -0,0 +1,5 @@
+torch==1.3.0
+torchvision>=0.4
+numpy==1.17.2
+matplotlib==3.1.1
+scipy==1.3.1

trainings/variational_inference_core_team/vae.ipynb

@@ -6,26 +6,416 @@
    "source": [
     "# Variational Autoencoders\n",
     "\n",
-    "A Bayesian approach\n",
+    "A *Neural Network* **and** a *Bayesian* approach\n",
     "\n",
-    "\n",
-    "## Suddenly a wild fashion-set appears...\n",
     "<img src=\"img/fashion-mnist-sprite.png\"/>"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Dataset and utilities\n",
+    "\n",
+    "Below we'll import some libraries and use some utility classes from `pytorch` and `torchvision`.\n",
+    "\n",
+    "<img src=\"img/auto-encoder.png\" width=\"900px\"> </img>\n",
+    "<small>[1]</small>"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import torchvision\n",
+    "from torchvision import transforms\n",
+    "import torch\n",
+    "from torch import nn\n",
+    "import torch.nn.functional as F\n",
+    "from torch.utils.data import DataLoader\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Load data\n",
+    "data_dir = './data'\n",
+    "t = transforms.ToTensor()\n",
+    "\n",
+    "fashion = False\n",
+    "if fashion:\n",
+    "    dset = torchvision.datasets.FashionMNIST\n",
+    "else:\n",
+    "    dset = torchvision.datasets.MNIST\n",
+    "\n",
+    "ds_train = dset(\n",
+    "    root=data_dir,\n",
+    "    train=True,\n",
+    "    transform=t,\n",
+    "    download=True\n",
+    ")\n",
+    "\n",
+    "ds_test = dset(\n",
+    "    root=data_dir,\n",
+    "    transform=t,\n",
+    "    train=False\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[torch.Tensor, int]"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# Show types in dataset\n",
+    "[type(a) for a in next(iter(ds_train))]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x7fd66b698710>"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAOS0lEQVR4nO3df4xU9bnH8c8jgqgQg7JQYsnd3kZNjcnd4kiuQQiXegnyDxDsTUlsaCTdxh9JMcRcszex/kgMMZdWjKbJ9oLQm15rFRBMzC1KSAyJVkdFBfF31rIFYYlKhSgt8Nw/9nCz4sx3lpkzc4Z93q9kMzPnOWfP47gfzsx8z5mvubsAjHznFN0AgNYg7EAQhB0IgrADQRB2IIhzW7mziRMnemdnZyt3CYTS19enQ4cOWaVaQ2E3s3mSVksaJem/3H1lav3Ozk6Vy+VGdgkgoVQqVa3V/TLezEZJelTSDZKulLTEzK6s9/cBaK5G3rNPl/SBu3/k7n+T9HtJC/JpC0DeGgn7pZL2Dnncny37GjPrNrOymZUHBgYa2B2ARjQS9kofAnzj3Ft373X3kruXOjo6GtgdgEY0EvZ+SVOHPP62pH2NtQOgWRoJ+yuSLjOz75jZGEk/krQln7YA5K3uoTd3P25mt0v6owaH3ta6++7cOgOQq4bG2d39WUnP5tQLgCbidFkgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCaGgWV7S/kydPJuvHjh1r6v7Xr19ftXb06NHktm+//Xay/tBDDyXrPT09VWuPPPJIctvzzz8/WV+1alWyfssttyTrRWgo7GbWJ+kLSSckHXf3Uh5NAchfHkf2f3H3Qzn8HgBNxHt2IIhGw+6StprZq2bWXWkFM+s2s7KZlQcGBhrcHYB6NRr2Ge4+TdINkm4zs1mnr+Duve5ecvdSR0dHg7sDUK+Gwu7u+7Lbg5I2SZqeR1MA8ld32M3sQjMbf+q+pLmSduXVGIB8NfJp/GRJm8zs1O/5H3f/31y6GmEOHz6crJ84cSJZf+ONN5L1rVu3Vq19/vnnyW17e3uT9SJ1dnYm6ytWrEjW16xZU7V20UUXJbedOXNmsj5nzpxkvR3VHXZ3/0jSP+XYC4AmYugNCIKwA0EQdiAIwg4EQdiBILjENQf9/f3JeldXV7L+2Wef5dnOWeOcc9LHmtTQmVT7MtRly5ZVrU2aNCm57bhx45L1s/FsUI7sQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAE4+w5uOSSS5L1yZMnJ+vtPM4+d+7cZL3Wf/vGjRur1s4777zktrNnz07WcWY4sgNBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIyz56DWddXr1q1L1p966qlk/dprr03WFy9enKynXHfddcn65s2bk/UxY8Yk65988knV2urVq5PbIl8c2YEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCHP3lu2sVCp5uVxu2f7OFseOHUvWa41l9/T0VK09+OCDyW23b9+erM+aNStZR3splUoql8tWqVbzyG5ma83soJntGrLsYjN7zszez24n5NkwgPwN52X8OknzTlt2l6Rt7n6ZpG3ZYwBtrGbY3f0FSZ+etniBpPXZ/fWSFubcF4Cc1fsB3WR33y9J2W3VibPMrNvMymZWHhgYqHN3ABrV9E/j3b3X3UvuXjobJ8MDRop6w37AzKZIUnZ7ML+WADRDvWHfImlpdn+ppPR1kAAKV/N6djN7XNJsSRPNrF/SLyStlPQHM1sm6c+SftjMJke6Wt+fXsuECfWPfD788MPJ+syZM5N1s4pDumhDNcPu7kuqlH6Qcy8AmojTZYEgCDsQBGEHgiDsQBCEHQiCr5IeAZYvX1619vLLLye33bRpU7K+e/fuZP2qq65K1tE+OLIDQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBCMs48Aqa+a7u3tTW67bdu2ZH3BggXJ+sKF6a8fnDFjRtXaokWLktty+Wy+OLIDQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBBM2Rxcrevd5807fU7Przt8+HDd+167dm2yvnjx4mR93Lhxde97pGpoymYAIwNhB4Ig7EAQhB0IgrADQRB2IAjCDgTB9ezBTZ8+PVmv9b3xd9xxR7L+5JNPVq3dfPPNyW0//PDDZP3OO+9M1sePH5+sR1PzyG5ma83soJntGrLsHjP7i5ntzH7mN7dNAI0azsv4dZIqnUb1K3fvyn6ezbctAHmrGXZ3f0HSpy3oBUATNfIB3e1m9mb2Mn9CtZXMrNvMymZWHhgYaGB3ABpRb9h/Lem7krok7Ze0qtqK7t7r7iV3L3V0dNS5OwCNqivs7n7A3U+4+0lJv5GU/kgXQOHqCruZTRnycJGkXdXWBdAeal7PbmaPS5otaaKkA5J+kT3ukuSS+iT9zN3319oZ17OPPF999VWy/tJLL1WtXX/99clta/1t3njjjcn6E088kayPRKnr2WueVOPuSyosXtNwVwBaitNlgSAIOxAEYQeCIOxAEIQdCIJLXNGQsWPHJuuzZ8+uWhs1alRy2+PHjyfrTz/9dLL+7rvvVq1dccUVyW1HIo7sQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAE4+xI2rdvX7K+cePGZP3FF1+sWqs1jl7LNddck6xffvnlDf3+kYYjOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EwTj7CFdryq1HH300WX/ssceS9f7+/jPuabhqXe/e2dmZrJtV/EblsDiyA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQjLOfBY4cOZKsP/PMM1Vr9913X3Lb9957r66e8jBnzpxkfeXKlcn61VdfnWc7I17NI7uZTTWz7Wa2x8x2m9nPs+UXm9lzZvZ+djuh+e0CqNdwXsYfl7TC3b8n6Z8l3WZmV0q6S9I2d79M0rbsMYA2VTPs7r7f3V/L7n8haY+kSyUtkLQ+W229pIXNahJA487oAzoz65T0fUl/kjTZ3fdLg/8gSJpUZZtuMyubWbnWedoAmmfYYTezcZI2SFru7n8d7nbu3uvuJXcvdXR01NMjgBwMK+xmNlqDQf+du5/6OtEDZjYlq0+RdLA5LQLIQ82hNxu8TnCNpD3u/sshpS2Slkpamd1ubkqHI8DRo0eT9b179ybrN910U7L++uuvn3FPeZk7d26yfu+991at1foqaC5RzddwxtlnSPqxpLfMbGe2rEeDIf+DmS2T9GdJP2xOiwDyUDPs7r5DUrV/Yn+QbzsAmoXTZYEgCDsQBGEHgiDsQBCEHQiCS1yH6csvv6xaW758eXLbHTt2JOvvvPNOXT3lYf78+cn63Xffnax3dXUl66NHjz7jntAcHNmBIAg7EARhB4Ig7EAQhB0IgrADQRB2IIgw4+x9fX3J+gMPPJCsP//881VrH3/8cT0t5eaCCy6oWrv//vuT2956663J+pgxY+rqCe2HIzsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBBFmnH3Dhg3J+po1a5q272nTpiXrS5YsSdbPPTf9v6m7u7tqbezYscltEQdHdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0Iwtw9vYLZVEm/lfQtSScl9br7ajO7R9JPJQ1kq/a4+7Op31UqlbxcLjfcNIDKSqWSyuVyxVmXh3NSzXFJK9z9NTMbL+lVM3suq/3K3f8zr0YBNM9w5mffL2l/dv8LM9sj6dJmNwYgX2f0nt3MOiV9X9KfskW3m9mbZrbWzCZU2abbzMpmVh4YGKi0CoAWGHbYzWycpA2Slrv7XyX9WtJ3JXVp8Mi/qtJ27t7r7iV3L3V0dOTQMoB6DCvsZjZag0H/nbtvlCR3P+DuJ9z9pKTfSJrevDYBNKpm2M3MJK2RtMfdfzlk+ZQhqy2StCv/9gDkZTifxs+Q9GNJb5nZzmxZj6QlZtYlySX1SfpZUzoEkIvhfBq/Q1KlcbvkmDqA9sIZdEAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSBqfpV0rjszG5D08ZBFEyUdalkDZ6Zde2vXviR6q1eevf2Du1f8/reWhv0bOzcru3upsAYS2rW3du1Lord6tao3XsYDQRB2IIiiw95b8P5T2rW3du1Lord6taS3Qt+zA2idoo/sAFqEsANBFBJ2M5tnZu+a2QdmdlcRPVRjZn1m9paZ7TSzQueXzubQO2hmu4Ysu9jMnjOz97PbinPsFdTbPWb2l+y522lm8wvqbaqZbTezPWa228x+ni0v9LlL9NWS563l79nNbJSk9yT9q6R+Sa9IWuLub7e0kSrMrE9Syd0LPwHDzGZJOiLpt+5+VbbsQUmfuvvK7B/KCe7+723S2z2SjhQ9jXc2W9GUodOMS1oo6Scq8LlL9PVvasHzVsSRfbqkD9z9I3f/m6TfS1pQQB9tz91fkPTpaYsXSFqf3V+vwT+WlqvSW1tw9/3u/lp2/wtJp6YZL/S5S/TVEkWE/VJJe4c87ld7zffukraa2atm1l10MxVMdvf90uAfj6RJBfdzuprTeLfSadOMt81zV8/0540qIuyVppJqp/G/Ge4+TdINkm7LXq5ieIY1jXerVJhmvC3UO/15o4oIe7+kqUMef1vSvgL6qMjd92W3ByVtUvtNRX3g1Ay62e3Bgvv5f+00jXelacbVBs9dkdOfFxH2VyRdZmbfMbMxkn4kaUsBfXyDmV2YfXAiM7tQ0ly131TUWyQtze4vlbS5wF6+pl2m8a42zbgKfu4Kn/7c3Vv+I2m+Bj+R/1DSfxTRQ5W+/lHSG9nP7qJ7k/S4Bl/W/V2Dr4iWSbpE0jZJ72e3F7dLine truncated
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "plt.imshow(next(iter(ds_train))[0].reshape(28, 28), cmap=\"Greys\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{'img': (torch.float32, torch.Size([64, 1, 28, 28])),\n",
+       " 'label': (torch.int64, torch.Size([64]))}"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "batch_size = 64\n",
+    "\n",
+    "# Create dataloaders utilities\n",
+    "dl_train = DataLoader(ds_train, shuffle=True, batch_size=batch_size)\n",
+    "dl_test = DataLoader(ds_test, batch_size=batch_size)\n",
+    "\n",
+    "# img and label\n",
+    "{k: (v.dtype, v.shape) for (k, v) in zip(['img', 'label'], next(iter(dl_train)))}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def loss_batch(model, criterion, x, opt=None):\n",
+    "    reconstruction_x = model(x)\n",
+    "    loss = criterion(reconstruction_x, x.view(-1, 784))\n",
+    "\n",
+    "    if opt is not None:\n",
+    "        loss.backward()\n",
+    "        opt.step()\n",
+    "        opt.zero_grad()\n",
+    "\n",
+    "    return loss.item(), len(x)\n",
+    "\n",
+    "def fit(epochs, model, criterion, optim, dl_train, dl_test, device='cuda', log_between_epochs=True):\n",
+    "    \"\"\"\n",
+    "    :param epoch: int\n",
+    "    :param model: torch.nn.Module\n",
+    "    :param criterion: torch.optim.Optimizer\n",
+    "    :param dl_train: torch.utils.data.DataLoader\n",
+    "    :param dl_test: torch.utils.data.DataLoader\n",
+    "    :param device: str\n",
+    "        'cuda' or 'cpu'\n",
+    "    \"\"\"\n",
+    "    for epoch in range(epochs):\n",
+    "        \n",
+    "        # train\n",
+    "        model.train()\n",
+    "        for i, (x, y) in enumerate(dl_train):\n",
+    "            loss, n = loss_batch(model, criterion, x.to(device), optim)\n",
+    "            \n",
+    "            if i % 150 == 0 and log_between_epochs:\n",
+    "                print(\"Loss: {:6.5f},\\t iteration: {:4d}\".format(loss / n, i))\n",
+    "        \n",
+    "        # evaluate test set\n",
+    "        model.eval()\n",
+    "        with torch.no_grad():\n",
+    "            losses, nums = zip(\n",
+    "                *[loss_batch(model, criterion, x.to(device)) for x, y in dl_test]\n",
+    "            )\n",
+    "        val_loss = np.sum(np.multiply(losses, nums)) / np.sum(nums)\n",
+    "        print(\"\\nEpoch: {:4d},\\t Test loss: {:6.5f}\".format(epoch, val_loss))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "class AutoEncoder(nn.Module):\n",
+    "    def __init__(self, input_size=784, z_size=20):\n",
+    "        super().__init__()\n",
+    "        # Implement autoencoder\n",
+    "    \n",
+    "    def forward(self, x):\n",
+    "        x = x.view(-1, 784)\n",
+    "        z = self.encoder(x)\n",
+    "        x = self.decoder(z)\n",
+    "        \n",
+    "        if self.training:\n",
+    "            return x\n",
+    "        else:\n",
+    "            return F.sigmoid(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "class AutoEncoder(nn.Module):\n",
+    "    def __init__(self, input_size=784, z_size=20):\n",
+    "        super().__init__()\n",
+    "        hidden_size = int((input_size - z_size) / 2 + z_size)\n",
+    "        self.encoder = nn.Sequential(\n",
+    "            nn.Linear(input_size, hidden_size),\n",
+    "            nn.ReLU(),\n",
+    "            nn.Linear(hidden_size, z_size)\n",
+    "        )\n",
+    "        self.decoder = nn.Sequential(\n",
+    "            nn.Linear(z_size, hidden_size),\n",
+    "            nn.ReLU(),\n",
+    "            nn.Linear(hidden_size, input_size),\n",
+    "        )\n",
+    "    \n",
+    "    def forward(self, x):\n",
+    "        x = x.view(-1, 784)\n",
+    "        z = self.encoder(x)\n",
+    "        x = self.decoder(z)\n",
+    "        \n",
+    "        if self.training:\n",
+    "            return x\n",
+    "        else:\n",
+    "            return F.sigmoid(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "On cuda!\n"
+     ]
+    }
+   ],
+   "source": [
+    "model = AutoEncoder()\n",
+    "optim = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "criterion = nn.BCEWithLogitsLoss()\n",
+    "epochs = 1\n",
+    "device = 'cpu'\n",
+    "\n",
+    "if torch.cuda.is_available():\n",
+    "    print('On cuda!')\n",
+    "    model.cuda()\n",
+    "    device = 'cuda'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Loss: 0.01081,\t iteration:    0\n",
+      "Loss: 0.00195,\t iteration:  150\n",
+      "Loss: 0.00163,\t iteration:  300\n",
+      "Loss: 0.00167,\t iteration:  450\n",
+      "Loss: 0.00165,\t iteration:  600\n",
+      "Loss: 0.00162,\t iteration:  750\n",
+      "Loss: 0.00165,\t iteration:  900\n"
+     ]
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/opt/miniconda3/lib/python3.7/site-packages/torch/nn/functional.py:1351: UserWarning: nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\n",
+      "  warnings.warn(\"nn.functional.sigmoid is deprecated. Use torch.sigmoid instead.\")\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Epoch:    0,\t Test loss: 0.67246\n"
+     ]
+    }
+   ],
+   "source": [
+    "fit(epochs, model, criterion, optim, dl_train, dl_test, device)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAOyklEQVR4nO3db4hV953H8c83o0biH6JxzA5qdrpNHjQRouZGlmQp2ZQtMQloH7TUQLEQ1jyI0IIPVmygBoKEZWvpg6VhukptaCyF1ijE7FakIWkeNLkGGzWS6Iqr0wzOqBhj/jmj3z6Yk2Vi5vzOeM+991zn+37BcGfO9/7mfL3jZ86d+7vn/MzdBWDyu6HqBgC0B2EHgiDsQBCEHQiCsANBTGnnzubNm+e9vb3t3CUQyokTJ3TmzBkbr1Yq7Gb2kKSfSeqS9F/u/mzq/r29varX62V2CSChVqvl1hp+Gm9mXZL+U9IKSXdKWm1mdzb6/QC0Vpm/2ZdLOubux939kqTfSFrZnLYANFuZsC+QdGrM1/3Zti8ws7VmVjez+tDQUIndASijTNjHexHgS++9dfc+d6+5e627u7vE7gCUUSbs/ZIWjfl6oaT3y7UDoFXKhP1NSXeY2VfMbJqk70ra3Zy2ADRbw1Nv7j5iZusk/Y9Gp962ufvhpnUGoKlKzbO7+x5Je5rUC4AW4u2yQBCEHQiCsANBEHYgCMIOBEHYgSAIOxAEYQeCIOxAEIQdCIKwA0EQdiAIwg4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBEHYgCMIOBFFqFVdMjLsn62bWpk4QWamwm9kJSR9KuixpxN1rzWgKQPM148j+z+5+pgnfB0AL8Tc7EETZsLukP5jZfjNbO94dzGytmdXNrD40NFRydwAaVTbs97v7MkkrJD1pZl+/+g7u3ufuNXevdXd3l9wdgEaVCru7v5/dDkraKWl5M5oC0HwNh93MZpjZrM8/l/RNSYea1RiA5irzavytknZmc8RTJL3g7v/dlK46UGqu/PLly8mxRfWPP/44WZ89e3ay3tXVlawDUomwu/txSXc3sRcALcTUGxAEYQeCIOxAEIQdCIKwA0FMmlNci04j/eyzz5L1Dz74IFl/5513cmv79+9Pjj1//nyyfvz48WS96N/24IMP5tYeeeSR5Nibb745Wb/hhnLHg9S0Y9GU47Rp05L1ot5uvPHG3NrUqVOTYyfjaccc2YEgCDsQBGEHgiDsQBCEHQiCsANBEHYgiEkzz37lypVS44sumfXyyy/n1t57773k2P7+/mT97NmzyXqRU6dO5db6+vqSYwcGBpL1ovcfzJw5M1lPzVcXzZOPjIwk60W93X13/kmZL7zwQnJsb29vsl72/QdVuP46BtAQwg4EQdiBIAg7EARhB4Ig7EAQhB0IYtLMsxfNexadG3377bcn64899lhu7aabbkqOHRwcTNbffffdZH3hwoXJ+p49e3JrFy9eTI49c6bcmpyXLl1K1lM/lxkzZiTHfvrpp8l60TUKDh48mFt7/fXXk2Nvu+22ZJ15dgAdi7ADQRB2IAjCDgRB2IEgCDsQBGEHgpg08+xF1/kuqk+fPj1ZX7x48TX39Lmic6NrtVqyPjw8nKwvW7Yst3bu3Lnk2I0bNybrRfP0qWuzS1JPT09ureg6AJs3b07Wd+3alaynrnFQ9N6FybgMduGR3cy2mdmgmR0as22ume01s6PZ7ZzWtgmgrIk8jf+lpIeu2rZB0j53v0PSvuxrAB2sMOzu/qqkq58LrpS0Pft8u6RVTe4LQJM1+gLdre4+IEnZ7fy8O5rZWjOrm1m96DpvAFqn5a/Gu3ufu9fcvdbd3d3q3QHI0WjYT5tZjyRlt+nTugBUrtGw75a0Jvt8jaT0HAiAyhXOs5vZDkkPSJpnZv2SfizpWUm/NbPHJZ2U9O1WNtkJpkxp/C0JReurF9WL5nxT9aL118v8u6Ry65jfddddyfqFCxdK7Tv13ol777231Pe+HhX+pN19dU7pG03uBUAL8XZZIAjCDgRB2IEgCDsQBGEHgpg0p7h2sqJpnLLLTaemz4qm1lo9xZSaVjx8+HBy7BtvvJGsF/X+6KOP5taKLv89GXFkB4Ig7EAQhB0IgrADQRB2IAjCDgRB2IEgmGfvAEXzxUWnwHby8sEjIyO5tZUrVybHFl3G+pZbbknWU5ei7uTHrFXi/YuBoAg7EARhB4Ig7EAQhB0IgrADQRB2IAjm2TtA0Tx7meWDi+boyyr6/qlz1s+ePZscW/S4rF6dd+HjUQsWLEjWo+HIDgRB2IEgCDsQBGEHgiDsQBCEHQiCsANBMM/eAYrmk4vOvU6dM150Tfqy142/dOlSsr5jx47c2uXLl5Nji85Xf/rpp5P1ybjschmFR3Yz22Zmg2Z2aMy2TWb2VzM7kH083No2AZQ1kafxv5T00Djbf+ruS7KPPc1tC0CzFYbd3V+VdK4NvQBooTIv0K0zs7ezp/lz8u5kZmvNrG5m9aGhoRK7A1BGo2H/uaSvSloiaUDST/Lu6O597l5z91p3d3eDuwNQVkNhd/fT7n7Z3a9I+oWk5c1tC0CzNRR2M+sZ8+W3JB3Kuy+AzlA4z25mOyQ9IGmemfVL+rGkB8xsiSSXdELSEy3sMbwy8/Bl5/CL5sIPHUr/nt+6dWvD+96wYUOyPnv27GQdX1QYdncf7woB+T9BAB2Jt8sCQRB2IAjCDgRB2IEgCDsQBKe4TgJllh8uuhT0Rx99lKxv2rQpWT9//nxubf78+cmx69atS9aLLrHNKa5fxJEdCIKwA0EQdiAIwg4EQdiBIAg7EARhB4Jgnn0SSM0nF82jF11q+sUXX0zWX3vttWR92rRpubWNGzc2PFZiHv1acWQHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSCYZ5/kiuaiz549m6xv27at1P5Xrx7v4sSj1qxZkxzLPHpzcWQHgiDsQBCEHQiCsANBEHYgCMIOBEHYgSCYZ5/kipZc3rlzZ7J+9OjRZH3p0qXJ+qpVq3JrM2fOTI5FcxUe2c1skZn90cyOmNlhM/tBtn2ume01s6PZ7ZzWtwugURN5Gj8iab27f03SP0p60szulLRB0j53v0PSvuxrAB2qMOzuPuDub2WffyjpiKQFklZK2p7dbbuk/OdrACp3TS/QmVmvpKWS/izpVncfkEZ/IUgad+EuM1trZnUzqw8NDZXrFkDDJhx2M5sp6XeSfujuFyY6zt373L3m7rXu7u5GegTQBBMKu5lN1WjQf+3uv882nzaznqzeI2mwNS0CaIbCqTcbPc9wq6Qj7r5lTGm3pDWSns1ud7WkQxRKXS765MmTybHPPfdcsj48PJys33fffcn6ihUrcmucwtpeE5lnv1/S9yQdNLMD2baNGg35b83scUknJX27NS0CaIbCsLv7nyTl/Qr+RnPbAdAqvF0WCIKwA0EQdiAIwg4EQdiBIDjF9TpQtOxy6nLQTzzxRHLssWPHkvXFixcn60899VSyPmUK/8U6BUd2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiCSdDrQNHloJ955pnc2iuvvJIc29XVlaxv2bIlWZ81a1ayjs7BkR0IgrADQRB2IAjCDgRB2IEgCDsQBGEHgmCevQOMjIwk6y+99FKy/vzzzzf8vYuWXL7nnnuSda79fv3gyA4EQdiBIAg7EARhB4Ig7EAQhB0IgrADQUxkffZFkn4l6e8kXZHU5+4/M7NNkv5V0lB2143uvqdVjV7Piq77XnTt9vXr1yfrFy5cyK0Vna++efPmZJ3rvk8eE/lJjkha7+5vmdksSfvNbG9W+6m7/0fr2gPQLBNZn31A0kD2+YdmdkTSglY3BqC5rulvdjPrlbRU0p+zTevM7G0z22Zmc3LGrDWzupnVh4aGxrsLgDaYcNjNbKak30n6obtfkPRzSV+VtESjR/6fjDfO3fvcvebute7u7ia0DKAREwq7mU3VaNB/7e6/lyR3P+3ul939iqRfSFreujYBlFUYdhs9rWmrpCPuvmXM9p4xd/uWpEPNbw9As0zk1fj7JX1P0kEzO5Bt2yhptZktkeSSTkhKrw0c2CeffJKs7927N1kfGBhI1lOnsc6dOzc5tmhqrugy1kXj0Tkm8mr8nySNd9Iyc+rAdYR30AFBEHYgCMIOBEHYgSAIOxAEYQeCsKLTL5upVqt5vV5v2/6uF0U/g6L68PBww2OnT5+erOP6UqvVVK/Xx72+N0d2IAjCDgRB2IEgCDsQBGEHgiDsQBCEHQiirfPsZjYk6f/GbJon6UzbGrg2ndpbp/Yl0Vujmtnb37v7uNd/a2vYv7Rzs7q71yprIKFTe+vUviR6a1S7euNpPBAEYQeCqDrsfRXvP6VTe+vUviR6a1Rbeqv0b3YA7VP1kR1AmxB2IIhKwm5mD5nZu2Z2zMw2VNFDHjM7YWYHzeyAmVV68n22ht6gmR0as22ume01s6PZ7bhr7FXU2yYz+2v22B0ws4cr6m2Rmf3RzI6Y2WEz+0G2vdLHLtFXWx63tv/NbmZdkt6T9C+S+iW9KWm1u7/T1kZymNkJSTV3r/wNGGb2dUkXJf3K3Rdn2/5d0jl3fzb7RTnH3f+tQ3rbJOli1ct4Z6sV9YxdZlzSKknfV4WPXaKv76gNj1sVR/blko65+3F3vyTpN5JWVtBHx3P3VyWdu2rzSknbs8+3a/Q/S9vl9NYR3H3A3d/KPv9Q0ufLjFf62CX6aosqwr5A0qkxX/ers9Z7d0l/MLP9Zra26mbGcau7D0ij/3kkza+4n6sVLuPdTlctM94xj10jy5+XVUXYx7s+VifN/93v7sskrZD0ZPZ0FRMzoWW822WcZcY7QqPLn5dVRdj7JS0a8/VCSe9X0Me43P397HZQ0k513lLUpz9fQTe7Hay4n//Line truncated
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.image.AxesImage at 0x7fd67021bc50>"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAANMElEQVR4nO3dXahd9ZnH8d9vYqPBFs0xRw1p9MQieHRwknKIQaU4lAm+XMRcODRKyaBMeqHSYi98mYtGQQzDtDUXQyGdxKTasRTamAgyNoSKKWjwKGc0meAcjWea1JjsEDBWhGryzMVZmTnGs9fZ7rX2S/J8P3DYe69nvTxs8svae//X3n9HhACc/f6q1w0A6A7CDiRB2IEkCDuQBGEHkjinmwebN29eDA0NdfOQQCoTExM6evSop6tVCrvtmyWtlzRL0r9FxLqy9YeGhjQ6OlrlkABKjIyMNK21/TLe9ixJ/yrpFklXS1pl++p29wegs6q8Z18q6Z2I2B8Rf5H0K0kr6mkLQN2qhH2BpANTHh8sln2O7TW2R22PNhqNCocDUEWVsE/3IcAXrr2NiA0RMRIRI4ODgxUOB6CKKmE/KGnhlMdfl/R+tXYAdEqVsL8m6Urbi2zPlvQdSdvraQtA3doeeouIz2zfJ+lFTQ69bYqIvbV1BqBWlcbZI+IFSS/U1AuADuJyWSAJwg4kQdiBJAg7kARhB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASlaZstj0h6SNJJyR9FhEjdTQFoH6Vwl7424g4WsN+AHQQL+OBJKqGPST9zvbrttdMt4LtNbZHbY82Go2KhwPQrqphvyEivinpFkn32v7W6StExIaIGImIkcHBwYqHA9CuSmGPiPeL2yOStkpaWkdTAOrXdthtn2/7a6fuS1ouaU9djQGoV5VP4y+RtNX2qf38e0T8Ry1dAahd22GPiP2S/qbGXgB0EENvQBKEHUiCsANJEHYgCcIOJFHHF2FSePXVV5vW1q9fX7rtggULSutz5swpra9evbq0PjAw0FYNuXBmB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEkGGdvUdlY9/j4eEeP/fjjj5fWL7jggqa1ZcuW1d3OGWNoaKhp7eGHHy7d9rLLLqu5m97jzA4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDO3qLnnnuuaW1sbKx022uuuaa0vnfv3tL67t27S+vbtm1rWnvxxRdLt120aFFp/b333iutV3HOOeX//ObPn19aP3DgQNvHLhuDl6QHH3yw7X33K87sQBKEHUiCsANJEHYgCcIOJEHYgSQIO5AE4+wtGh4ebqvWimuvvba0vmrVqtL6unXrmtYmJiZKt51pnH3//v2l9Spmz55dWp9pnH2m3huNRtPaVVddVbrt2WjGM7vtTbaP2N4zZdmA7R22x4vbuZ1tE0BVrbyM3yzp5tOWPSRpZ0RcKWln8RhAH5sx7BHxsqRjpy1eIWlLcX+LpNtr7gtAzdr9gO6SiDgkScXtxc1WtL3G9qjt0bL3UAA6q+OfxkfEhogYiYiRwcHBTh8OQBPthv2w7fmSVNweqa8lAJ3Qbti3Szr128qrJTX/jiWAvjDjOLvtZyXdJGme7YOSfiRpnaRf275H0h8l3dHJJlHuvPPOa1qrOp5c9RqCKmb6Hv/Ro0dL69ddd13T2vLly9vq6Uw2Y9gjotkVHd+uuRcAHcTlskAShB1IgrADSRB2IAnCDiTBV1zRMx9//HFpfeXKlaX1kydPltaffPLJprU5c+aUbns24swOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0kwzo6e2bx5c2n9gw8+KK1fdNFFpfXLL7/8y7Z0VuPMDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM6Ojnr33Xeb1h544IFK+37llVdK65deemml/Z9tOLMDSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKMs6Ojnn/++aa1Tz/9tHTbO+4onwn8iiuuaKunrGY8s9veZPuI7T1Tlq21/SfbY8XfrZ1tE0BVrbyM3yzp5mmW/zQiFhd/L9TbFoC6zRj2iHhZ0rEu9AKgg6p8QHef7TeLl/lzm61ke43tUdujjUajwuEAVNFu2H8m6RuSFks6JOnHzVaMiA0RMRIRI4ODg20eDkBVbYU9Ig5HxImIOCnp55KW1tsWgLq1FXbb86c8XClpT7N1AfSHGcfZbT8r6SZJ82wflPQjSTfZXiwpJE1I+l4He0Qfm2msfOvWrU1r5557bum2TzzxRGl91qxZpXV83oxhj4hV0yze2IFeAHQQl8sCSRB2IAnCDiRB2IEkCDuQBF9xRSUbN5YPzOzatatp7c477yzdlq+w1oszO5AEYQeSIOxAEoQdSIKwA0kQdiAJwg4kwTg7So2NjZXW77///tL6hRde2LT22GOPtdUT2sOZHUiCsANJEHYgCcIOJEHYgSQIO5AEYQeSYJw9uU8++aS0vmrVdD8u/P9OnDhRWr/rrrua1vi+endxZgeSIOxAEoQdSIKwA0kQdiAJwg4kQdiBJBhnP8udPHmytH7bbbeV1t9+++3S+vDwcGn90UcfLa2je2Y8s9teaPv3tvfZ3mv7+8XyAds7bI8Xt3M73y6AdrXyMv4zST+MiGFJyyTda/tqSQ9J2hkRV0raWTwG0KdmDHtEHIqIN4r7H0naJ2mBpBWSthSrbZF0e6eaBFDdl/qAzvaQpCWSdku6JCIOSZP/IUi6uMk2a2yP2h5tNBrVugXQtpbDbvurkn4j6QcRcbzV7SJiQ0SMRMTI4OBgOz0CqEFLYbf9FU0G/ZcR8dti8WHb84v6fElHOtMigDrMOPRm25I2StoXET+ZUtouabWkdcXtto50iEqOHTtWWn/ppZcq7f/pp58urQ8MDFTaP+rTyjj7DZK+K+kt26d+RPwRTYb817bvkfRHSXd0pkUAdZgx7BHxB0luUv52ve0A6BQulwWSIOxAEoQdSIKwA0kQdiAJvuJ6Fvjwww+b1pYtW1Zp388880xpfcmSJZX2j+7hzA4kQdiBJAg7kARhB5Ig7EAShB1IgrADSTDOfhZ46qmnmtb2799fad833nhjaX3y5w5wJuDMDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJMM5+BhgfHy+tr127tjuN4IzGmR1IgrADSRB2IAnCDiRB2IEkCDuQBGEHkmhlfvaFkn4h6VJJJyVtiIj1ttdK+kdJjWLVRyLihU41mtmuXbtK68ePH29738PDw6X1OXPmtL1v9JdWLqr5TNIPI+IN21+T9LrtHUXtpxHxL51rD0BdWpmf/ZCkQ8X9j2zvk7Sg040BqNeXes9ue0jSEkm7i0X32X7T9ibbc5tss8b2qO3RRqMx3SoAuqDlsNv+qqTfSPpBRByX9DNJ35C0WJNn/h9Pt11EbIiIkYgYGRwcrKFlAO1oKey2v6LJoP8yIn4rSRFxOCJORMRJST+XtLRzbQKoasawe/LnQzdK2hcRP5myfP6U1VZK2lN/ewDq0sqn8TdI+q6kt2yPFcsekbTK9mJJIWlC0vc60iEquf7660vrO3bsKK0z9Hb2aOXT+D9Imu7HwRlTB84gXEEHJEHYgSQIO5AEYQeSIOxAEoQdSIKfkj4D3H333ZXqgMSZHUiDsANJEHYgCcIOJEHYgSQIO5AEYQeScER072B2Q9L/TFk0T9LRrjXw5fRrb/3al0Rv7aqzt8sjYtrff+tq2L9wcHs0IkZ61kCJfu2tX/uS6K1d3eqNl/FAEoQdSKLXYd/Q4+OX6dfe+rUvid7a1ZXeevqeHUD39PrMDqBLCDuQRE/Cbvtm22/bfsf2Q73ooRnbE7bfsj1me7THvWyyfcT2ninLBmzvsD1e3E47x16Peltr+0/Fczdm+9Ye9bbQ9u9t77O91/b3i+U9fe5K+urK89b19+y2Z0n6b0l/J+mgpNckrYqI/+pqI03YnpA0EhE9vwDD9rck/VnSLyLir4tl/yzpWESsK/6jnBsRD/ZJb2sl/bnX03gXsxXNnzrNuKTbJf2DevjclfT19+rC89aLM/tSSe9ExP6I+IukX0la0YM++l5EvCzp2GmLV0jaUtzfosl/LF3XpLe+EBGHIuKN4v5Hkk5NM97T566kr67oRdgXSDow5fFB9dd87yHpd7Zft72m181M45KIOCRN/uORdHGP+zndjNN4d9Np04z3zXPXzvTnVfUi7NNNJdVP4383RMQ3Jd0i6d7i5Spa09I03t0yzTTjfaHd6c+r6kXYD0paOOXx1yW934M+phUR7xe3RyRtVf9NRX341Ay6xe2RHvfzf/ppGu/pphlXHzx3vZz+vBdhf03SlbYX2Z4t6TuStvegjy+wfX7xwYlsny9pufpvKurtklYX91dL2tbDXj6nX6bxbjbNuHr83PV8+vOI6PqfpFs1+Yn8u5L+qRc9NOnrCkn/Wfzt7XVvkp7V5Mu6TzX5iugeSRdJ2ilpvLgd6KPenpb0lqQ3NRms+T3q7UZNvjV8U9JY8Xdrr5+7kr668rxxuSyQBFfQAUkQdiAJwg4kQdiBJAg7kARhB5Ig7EAS/wseauFUg51ZyQAAAABJRU5ErkJggg==\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "with torch.no_grad():\n",
+    "    model.eval()\n",
+    "    for x, _ in ds_test:\n",
+    "        x_reconstruction = model(x.to(device)).cpu()\n",
+    "        break\n",
+    "\n",
+    "plt.imshow(x_reconstruction.view(28, 28), cmap='Greys')\n",
+    "plt.show()\n",
+    "plt.imshow(x.view(28, 28), cmap='Greys')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Variational Autoencoder\n",
+    "\n",
+    "## *A neural networks perspective*\n",
+    "\n",
+    "Let $Q(z|x)$ be the encoder, and $P(x|z)$ the decoder. Both are neural networks.\n",
+    "\n",
+    "We want to be able to generate new data, so we modify the autoencoder so that it has a random process at it's core.\n",
+    "\n",
+    "\n",
+    "$$ \\begin{eqnarray}\n",
+    "Q_{\\theta}(\\cdot) &=& \\{ \\mu_{\\theta}, \\log \\sigma^2_{\\theta} \\}\\\\\n",
+    "z &\\sim& \\mathcal{N}(\\mu_{\\theta}, \\sigma^2_{\\theta})\n",
+    "\\end{eqnarray}$$\n",
+    "\n",
+    "Now it's of course very easy for $E(\\cdot)$ to cheat and become an autoencoder by just learning $E_{\\theta}(\\cdot) = \\{ \\mu_{\\theta}, 0 \\}$. Without variance the random process isn't quite random.\n",
+    "\n",
+    "Therefore, we need a regularization term, the Kullback-Leibler-Divergence. \n",
+    "\n",
+    "$$ D_{KL}(Q(z|x) || P(z)) $$\n",
+    "\n",
+    "Where we don't want the $z$ values we draw from the stochastic decoder to deviate too much from a unit Gaussian.\n",
+    "\n",
+    "$$ \\begin{eqnarray}\n",
+    "z &\\sim& P(z) \\\\ \n",
+    "  &\\sim& \\mathcal{N}(0, \\text{diag}(\\sigma^2)) \n",
+    "\\end{eqnarray}$$"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "## Generative model\n",
+    "## *A Bayesian perspective*\n",
     "\n",
     "<img src=\"./img/plate-vae.png\"/>\n",
     "\n",
     "For each image $x_i$:\n",
     "\n",
     "$$\\begin{eqnarray}\n",
-    "z_i &\\sim& P(z) \\\\\n",
-    "x_i &\\sim& P(x|z_i)\n",
+    "z_i &\\sim& \\mathcal{N}(\\mu_i, \\sigma_i^2) \\\\\n",
+    "x_i &\\sim& P(x_i|z_i)\n",
     "\\end{eqnarray}$$\n"
    ]
   },
@@ -58,7 +448,246 @@
    "metadata": {},
    "source": [
     "## KL-divergence\n",
-    "There is a reason we choose to measure the 'closeseness' with the KL-divergence and not a real distance metric, such as the Wasserstein distance."
+    "There is a reason we choose to measure the 'closeseness' with the KL-divergence and not a real distance metric, such as the Wasserstein distance.\n",
+    "\n",
+    "$$\\begin{eqnarray}\n",
+    "D_{\\text{KL}}(Q(z) \\: || \\: P(z|x)) &=& \\int_{z} Q(z) \\log \\frac{Q(z)}{P(z|x)}\\text{dz} \\\\\n",
+    "&=& \\int_{z} q(z) \\log \\frac{q(z)p(x)}{p(z, x)}\\text{d}z  &\\qquad \\tiny{\\text{rewrite the posterior as}\\frac{p(z, x)}{p(x)} } \\\\\n",
+    "&=& \\int_{z} Q(z) \\log \\frac{Q(z)}{P(z, x)}\\text{d}z + \\int_{z} Q(z) \\log P(x) \\text{d} z \\\\\n",
+    "&=&  \\int_{z} Q(z) \\log \\frac{Q(z)}{P(z, x)}\\text{d}z +  \\log P(x) &\\qquad \\tiny{P(x)\\text{ is not parameterized by }z \\text{and} \\int_{z} Q(z) \\text{d} z = 1} \\\\\n",
+    "&=& \\log P(x) -\\int_{z} Q(z) \\log \\frac{P(z, x)}{Q(z)}\\text{d}z &\\qquad \\tiny{\\log A = - \\log \\frac{q}{A}} \\\\\n",
+    "&=& \\log P(x) - \\underbrace{ E_{z \\sim Q}[\\log \\frac{P(z, x)}{Q(z)} }_{\\text{ELBO}} ]\n",
+    "\\end{eqnarray}\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Reconstruction loss and KL-divergence\n",
+    "In the derivation of Variational Autoencoders, the ELBO is often written in two terms; a reconstruction error and a KL-term on the variational distribution $Q_z$ and the prior $P(z)$\n",
+    "\n",
+    "First we rewrite the joint probability $P(z, x)$ into conditional probability $P(x|z)P(z)$.\n",
+    "\n",
+    "$$ \\text{ELBO}  = E_{z \\sim Q}[\\log \\frac{P(x|z)P(z)}{Q(z)} ] $$\n",
+    "\n",
+    "Then we expand the equation and thereby isolating the reconstruction error $\\log P(x|z)$, i.e. the log likelihood.\n",
+    "\n",
+    "$$ \\text{ELBO}  = E_{z \\sim Q}[\\log P(x|z)] +  E_{z \\sim Q}[\\log \\frac{P(z)}{Q(z)}]$$\n",
+    "\n",
+    "If we rewrite the $E_{z \\sim Q}[\\log \\frac{P(z)}{Q(z)}]$ in the integral form $\\int_{z} Q(z)\\log\\frac{P(z)}{Q(z)}dz$, we can observe that this is the KL-divergence between the prior $P(z)$ and the variational distribution $Q(z)$. Resulting in an ELBO defined by the reconstruction error and $-D_{KL}(Q(z)||P(z)).$\n",
+    "\n",
+    "$$ \\text{ELBO}  = E_{z \\sim Q}[\\log P(x|z)] - D_{KL}(Q(z)||P(z))$$\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Implementation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 102,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class VAE(nn.Module):\n",
+    "    def __init__(self, input_size=784, z_size=20):\n",
+    "        super().__init__()\n",
+    "        hidden_size = int((input_size - z_size) / 2 + z_size)\n",
+    "        self.z_size = z_size\n",
+    "        self.encoder = nn.Sequential(\n",
+    "            nn.Linear(input_size, hidden_size),\n",
+    "            nn.ReLU(),\n",
+    "            nn.Linear(hidden_size, hidden_size),\n",
+    "            nn.ReLU(),\n",
+    "            nn.Linear(hidden_size, z_size * 2)\n",
+    "        )\n",
+    "        self.decoder = nn.Sequential(\n",
+    "            nn.Linear(z_size, hidden_size),\n",
+    "            nn.ReLU(),\n",
+    "            nn.Linear(hidden_size, input_size),\n",
+    "            nn.Sigmoid()\n",
+    "        )\n",
+    "        \n",
+    "    def reparameterize(self, mu, log_var):\n",
+    "        std = torch.exp(0.5 * log_var)\n",
+    "        eps = torch.randn_like(std)  # unit gaussian\n",
+    "        z = mu + eps * std \n",
+    "        return z\n",
+    "        \n",
+    "    def forward(self, x):\n",
+    "        x = x.view(-1, 784)\n",
+    "        variational_params = self.encoder(x)\n",
+    "        mu = variational_params[..., :self.z_size]\n",
+    "        log_var = variational_params[..., self.z_size:]\n",
+    "        z = self.reparameterize(mu, log_var)\n",
+    "        return self.decoder(z), z, mu, log_var\n",
+    "\n",
+    "\n",
+    "def loss_batch(model, criterion, x, opt=None):\n",
+    "    reconstruction_x, z, mu, log_var = model(x)\n",
+    "    loss = criterion(reconstruction_x, x.view(-1, 784))\n",
+    "    \n",
+    "    # see Appendix B from VAE paper:\n",
+    "    # Kingma and Welling. Auto-Encoding Variational Bayes. ICLR, 2014\n",
+    "    # https://arxiv.org/abs/1312.6114\n",
+    "    # 0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)\n",
+    "    kl = -0.5 * torch.sum(1 + log_var - mu.pow(2) - log_var.exp())\n",
+    "    \n",
+    "    loss += kl\n",
+    "\n",
+    "    if opt is not None:\n",
+    "        loss.backward()\n",
+    "        opt.step()\n",
+    "        opt.zero_grad()\n",
+    "\n",
+    "    return loss.item(), len(x)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 116,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "On cuda!\n"
+     ]
+    }
+   ],
+   "source": [
+    "model = VAE()\n",
+    "optim = torch.optim.Adam(model.parameters(), lr=0.01)\n",
+    "criterion = nn.BCELoss(reduction='sum')\n",
+    "epochs = 25\n",
+    "device = 'cpu'\n",
+    "batch_size = 256\n",
+    "\n",
+    "dl_train = DataLoader(ds_train, shuffle=True, batch_size=batch_size)\n",
+    "dl_test = DataLoader(ds_test, batch_size=batch_size)\n",
+    "\n",
+    "if torch.cuda.is_available():\n",
+    "    print('On cuda!')\n",
+    "    model.cuda()\n",
+    "    device = 'cuda'\n",
+    "    \n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 117,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "Epoch:    0,\t Test loss: 45651.31713\n",
+      "\n",
+      "Epoch:    1,\t Test loss: 44757.27022\n",
+      "\n",
+      "Epoch:    2,\t Test loss: 43418.09588\n",
+      "\n",
+      "Epoch:    3,\t Test loss: 43290.16544\n",
+      "\n",
+      "Epoch:    4,\t Test loss: 43172.56996\n",
+      "\n",
+      "Epoch:    5,\t Test loss: 43022.13056\n",
+      "\n",
+      "Epoch:    6,\t Test loss: 43347.33919\n",
+      "\n",
+      "Epoch:    7,\t Test loss: 42884.59440\n",
+      "\n",
+      "Epoch:    8,\t Test loss: 42916.98715\n",
+      "\n",
+      "Epoch:    9,\t Test loss: 42681.94180\n",
+      "\n",
+      "Epoch:   10,\t Test loss: 42751.86753\n",
+      "\n",
+      "Epoch:   11,\t Test loss: 42673.31136\n",
+      "\n",
+      "Epoch:   12,\t Test loss: 42635.39283\n",
+      "\n",
+      "Epoch:   13,\t Test loss: 42672.09193\n",
+      "\n",
+      "Epoch:   14,\t Test loss: 42560.42911\n",
+      "\n",
+      "Epoch:   15,\t Test loss: 42571.33200\n",
+      "\n",
+      "Epoch:   16,\t Test loss: 42495.00136\n",
+      "\n",
+      "Epoch:   17,\t Test loss: 42640.09606\n",
+      "\n",
+      "Epoch:   18,\t Test loss: 42661.70678\n",
+      "\n",
+      "Epoch:   19,\t Test loss: 42861.99558\n",
+      "\n",
+      "Epoch:   20,\t Test loss: 42472.55930\n",
+      "\n",
+      "Epoch:   21,\t Test loss: 42630.67527\n",
+      "\n",
+      "Epoch:   22,\t Test loss: 42511.97979\n",
+      "\n",
+      "Epoch:   23,\t Test loss: 42701.39958\n",
+      "\n",
+      "Epoch:   24,\t Test loss: 42459.78887\n"
+     ]
+    }
+   ],
+   "source": [
+    "fit(epochs, model, criterion, optim, dl_train, dl_test, device, log_between_epochs=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Generate new samples"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 171,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAQoElEQVR4nO3dW2zd1ZXH8d/Kxc7FCUmICQ6FuCQRDIwYLhZCZFQxVFMFHoA+dFQkqoyEJn0AqZX6MIh5KI9oNG3Vh1GldEBNRx2qSi2CBzQDQpVQQaowUUoSAkMmmFwwsUMCCZB71jz4z4wJPnuZ8z+3Zn0/kmX7rLNztv/xz3/7rLP/29xdAC5+c7o9AQCdQdiBJAg7kARhB5Ig7EAS8zr5YCtXrvTh4eFOPiSQytjYmA4fPmwz1WqF3cw2SvqppLmS/s3dHy/df3h4WKOjo3UeEkDByMhIw1rTv8ab2VxJ/yrpLknXSbrfzK5r9t8D0F51/ma/VdIed9/r7qcl/VrSva2ZFoBWqxP2KyTtn/b5geq2zzGzzWY2amajk5OTNR4OQB11wj7TkwBfeO2tu29x9xF3HxkcHKzxcADqqBP2A5KunPb5VyS9V286ANqlTthflbTezL5qZn2Svi3p2dZMC0CrNd16c/ezZvawpP/SVOvtSXff1bKZAWipWn12d39O0nMtmguANuLlskAShB1IgrADSRB2IAnCDiRB2IEkOrqeHc2pcwXgaGzdqwubzbh0elb1aCxaizM7kARhB5Ig7EAShB1IgrADSRB2IAlabx1Qt/1Vpx6NPX/+fLEeidpnc+Y0Pp+UarP5tyO09j6PMzuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEGfvQXq9rLPnDlTrJ86dapY//jjjxvWjh071vRYKZ77okWLivXly5c3rC1ZsqQ4tr+/v1ifO3dusV7qs2fs8XNmB5Ig7EAShB1IgrADSRB2IAnCDiRB2IEk6LPPUqmXXrePfvjw4WJ9z549xfrLL7/csLZz587i2LGxsWI96sNHffbrr7++Ye22224rjr3hhhuK9aGhoWK91Mfv6+srjp03rxyNqMcf9fFL9Xb16GuF3czGJB2XdE7SWXcfacWkALReK87sf+Pu5VMTgK7jb3Ygibphd0nPm9lrZrZ5pjuY2WYzGzWz0cnJyZoPB6BZdcO+wd1vlnSXpIfM7GsX3sHdt7j7iLuPDA4O1nw4AM2qFXZ3f696PyHpaUm3tmJSAFqv6bCb2WIzW/LZx5K+Ianc5wHQNXWejV8l6emqJzhP0n+4+3+2ZFZdUOfa7OfOnSuO/fDDD4v1HTt2FOvPP/98sf7KK680rO3bt6849vTp08V69BqBqB999OjRhrVPP/20OPbQoUPFetSHX7t2bcPaihUrimMHBgaK9fnz5xfr0XHpxlbWTYfd3fdK+qsWzgVAG9F6A5Ig7EAShB1IgrADSRB2IAmWuFbqXA46utRztIQ1ar29+eabxXqphXX11VcXx65evbpYj5Zynjhxoun62bNni2NLbTtJ+uCDD4r10hLYSy65pDg20s5LRUffi80+Nmd2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUgiTZ896l3WGR/1iz/66KNiPdpWecGCBcX6tdde27A2PDxcHLtmzZpifeHChcV6tHz3yJEjDWuffPJJcWy0zDS6jHVpGWqd7Z5bUW/X2BLO7EAShB1IgrADSRB2IAnCDiRB2IEkCDuQRJo+e111LiVd93LM0U46ixcvblhbv359cew111xTrEdzi9aU79+/v2EtulR06euS4ss5l+p1LvXciXo7cGYHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSTS9Nmjvmada3XPmVP+mRn1g6N128uWLSvWS/3oVatWFcdGa+VL18uX4rX80ZbQJVEvPFrP3t/f37DW19dXHButd4/+z7vRR4+EZ3Yze9LMJsxs57TbVpjZC2b2dvV+eXunCaCu2fwa/wtJGy+47RFJL7r7ekkvVp8D6GFh2N39JUkXXlvoXklbq4+3SrqvxfMC0GLNPkG3yt3HJal6f1mjO5rZZjMbNbPRycnJJh8OQF1tfzbe3be4+4i7j0QLOgC0T7NhP2RmQ5JUvZ9o3ZQAtEOzYX9W0qbq402SnmnNdAC0S9hnN7OnJN0haaWZHZD0Q0mPS/qNmT0oaZ+kb7Vzkr2gznr2qGcb9Yujnm5pf/aDBw8Wx548ebJYj/aej9akT0w0/qUvem1D9GdfdNxKvfTo/6TdffZu9OHDsLv7/Q1KX2/xXAC0ES+XBZIg7EAShB1IgrADSRB2IIk0S1zrKi31jFpvddswx48fL9ZL2yZHl3pevry8YDG6DHa0HXVpW+aVK1cWx0atuTrLVC/G1lqEMzuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJJGmzx71bOvWS6KebXQ55miZ6ZEjF14i8P9FvehSH1wqL5+V4j5+qU8f9bqjx46OS53/s4sRZ3YgCcIOJEHYgSQIO5AEYQeSIOxAEoQdSCJNnz0SbU1cqkdrl6Mtm6P6woULi/XSls5Rjz/qRdddz146NkePHi2OPXDgQLE+NjZWrJe2q47W8Uei75de3NKZMzuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJJGmz97Otc3RuuwFCxYU66tXry7WozXj/f39DWvR1x31i6P18FEfvnRsTp8+XRz71ltvFevRcV23bl3DWnTM636/1Bnfrh58eGY3syfNbMLMdk677TEzO2hm26u3u9syOwAtM5tf438haeMMt//E3W+s3p5r7bQAtFoYdnd/SVLj6x4B+LNQ5wm6h83s9erX/IYvNDazzWY2amajk5OTNR4OQB3Nhv1nktZKulHSuKQfNbqju29x9xF3HxkcHGzy4QDU1VTY3f2Qu59z9/OSfi7p1tZOC0CrNRV2Mxua9uk3Je1sdF8AvSHss5vZU5LukLTSzA5I+qGkO8zsRkkuaUzSd9s4x5aou592qT5vXvkwlvrgUry2emhoqFgvPX7U743mXlorL8VzP3nyZMNa1Effu3dvsb5ixYpivbRevm4vuxf3X4+EYXf3+2e4+Yk2zAVAG/FyWSAJwg4kQdiBJAg7kARhB5JIs8Q1El36t057K7oUdNSai8YPDAw0rC1evLg4NlomGrWYoi2fDx8+3LC2f//+4thS204qb1Ut1WuX1r0UdC+25jizA0kQdiAJwg4kQdiBJAg7kARhB5Ig7EASafrsdfuepUsi111GGl2K+uzZs8V6qVceLQNdunRpsR7NPeqFly4XHX1d0XbQx44dK9ZLy2+jr6vdfXS2bAbQNoQdSIKwA0kQdiAJwg4kQdiBJAg7kESaPns7t2yO/u1z584V69HWxVEfftGiRQ1r0ZbLUb3uawBKl3N+5513imOjtfLRJbYvu+yyhrXo6/pzXK8e4cwOJEHYgSQIO5AEYQeSIOxAEoQdSIKwA0mk6bNHzp8/33T9zJkzxbHRuutSL1qSTpw4UazPnz+/YS3qVUevEYge+9133y3Wt23b1rA2NjZWHLtkyZJi/aabbirWBwcHG9Yuxj56JDyzm9mVZvZ7M9ttZrvM7HvV7SvM7AUze7t6X96oG0BXzebX+LOSfuDufyHpNkkPmdl1kh6R9KK7r5f0YvU5gB4Vht3dx919W/XxcUm7JV0h6V5JW6u7bZV0X7smCaC+L/UEnZkNS7pJ0h8lrXL3cWnqB4KkGV+IbGabzWzUzEYnJyfrzRZA02YddjMbkPRbSd939/IzTtO4+xZ3H3H3kdITJgDaa1ZhN7P5mgr6r9z9d9XNh8xsqKoPSZpozxQBtELYerOpHsQTkna7+4+nlZ6VtEnS49X7Z9oyw1mKWkhRPWq9lZZyRu2tiYnyz8GofTU+Pl6sl+Z+6tSp4tho+W20LXI09zfeeKNhLTrmGzZsKNY3btxYrJe2ur4YW2uR2fTZN0j6jqQdZra9uu1RTYX8N2b2oKR9kr7VnikCaIUw7O7+B0mNfgx+vbXTAdAuvFwWSIKwA0kQdiAJwg4kQdiBJFjiWon68KU+e7QM9Pjx48V61EfftWtXsb53796Gtffff784NrqMdbQlc6TU67799tuLYx944IFife3atcX6nDmcy6bjaABJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEvTZK1FPtlQv9ZIlafny8oV3L7/88mJ93759xXppzXp0GetovfuyZcuK9XXr1hXrt9xyS8PanXfeWRwb9eEHBgaK9Whb5mw4swNJEHYgCcIOJEHYgSQIO5AEYQeSIOxAEhdNn73udcCjPntfX1/D2tKlS4tjo3XXl156abFe6lVL0j333NOwFq2Vj655H81tzZo1Tdej1xcsWLCgWJ83r/ztm/Ha8CWc2YEkCDuQBGEHkiDsQBKEHUiCsANJEHYgidnsz36lpF9KulzSeUlb3P2nZvaYpH+QNFnd9VF3f65dE60r6rnWqUfrpvv7+4v1aM34VVddVazffPPNxXod7Txudfvg9NG/nNm8qOaspB+4+zYzWyLpNTN7oar9xN3/pX3TA9Aqs9mffVzSePXxcTPbLemKdk8MQGt9qb/ZzWxY0k2S/ljd9LCZvW5mT5rZjNdeMrPNZjZqZqOTk5Mz3QVAB8w67GY2IOm3kr7v7sck/UzSWkk3aurM/6OZxrn7FncfcfeRwcHBFkwZQDNmFXYzm6+poP/K3X8nSe5+yN3Puft5ST+XdGv7pgmgrjDsNvWU5xOSdrv7j6fdPjTtbt+UtLP10wPQKrN5Nn6DpO9I2mFm26vbHpV0v5ndKMkljUn6bltm2CPa2UICOmE2z8b/QdJM380921MH8EW8gg5IgrADSRB2IAnCDiRB2IEkCDuQBGEHkiDsQBKEHUiCsANJEHYgCcIOJEHYgSQIO5CEuXvnHsxsUtK7025aKelwxybw5fTq3HpLine truncated
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "with torch.no_grad():\n",
+    "    z = torch.randn(20)\n",
+    "    x = model.decoder(z)\n",
+    "    plt.imshow(x.view(28, 28), Cmap='Greys')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# References\n",
+    "[1] https://lilianweng.github.io/lil-log/2018/08/12/from-autoencoder-to-beta-vae.html <br>\n",
+    "[2] https://jaan.io/what-is-variational-autoencoder-vae-tutorial/ <br>\n",
+    "[3] https://www.ritchievink.com/blog/2019/06/10/bayesian-inference-how-we-are-able-to-chase-the-posterior/ <br>\n",
+    "[4] https://www.ritchievink.com/blog/2019/09/16/variational-inference-from-scratch/ <br>"
    ]
   },
   {