Yet another VAE note

vae.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume the following generative process"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "z_i \\sim p(z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "x_i \\sim p_{\\theta}(x|z_i)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Giving the following joint distribution\n",
    "\n",
    "$$\n",
    "p_{\\theta}(x,z) = p_{\\theta}(x|z)p(z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We wish to optimize $p_{\\theta}(x)$ the probability of the data under our model, with respect to the model parameters $\\theta$\n",
    "\n",
    "$$\n",
    "p_{\\theta}(x) = \\int_z p_{\\theta}(x|z)p(z)\n",
    "$$\n",
    "\n",
    "Which is equivialent to the expectation of $p_{\\theta}(x|z)$ under $p(z)$\n",
    "\n",
    "$$\n",
    "p_{\\theta}(x) = E_{z \\sim p(z)} p_{\\theta}(x|z)\n",
    "$$\n",
    "\n",
    "Now we could stop here and optimize this w.r.t. $\\theta$. That is entirely theoretically possible. Assume we have $N$ training samples $x$. For each of the N training samples $x$, we'd approximate $p_{\\theta}(x)$ by sampling $M$ $z$'s, and take the average of $p_{\\theta}(x|z)$. Given the approximate $p_{\\theta}(x)$ we'd take the gradient of it w.r.t. $\\theta$ and take a gradient step in the direction that maximizes $p_{\\theta}(x)$ and repeat this $T$ times (or until $\\theta$ converges). In math:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\Delta\\theta_t = \\nabla_{\\theta} \\left[\\frac{1}{N} \\sum_{i=0}^{N-1} \\frac{1}{M} \\sum_{j=0}^{M-1} p_{\\theta}(x_i|z_j \\sim p(z)) \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\theta_{t} = \\theta_{t-1} + \\alpha\\Delta\\theta_t\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this is theoretically possible, it is computationally infeasible for anything but very simple $p_{\\theta}(x|z)$ and $p(z)$, since we'd have to draw a very large amount, $M$, samples of $p(z)$ in order for the average of $p_{\\theta}(x|z)$ to be a good approximation of the expectation $E_{z \\sim p(z)}p_{\\theta}(x|z)$. It is important to understand why we'd need a large number of samples $M$ if we were to draw them from the prior. A computational example at this point is instructive. We'll assume the following generative process:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(z) = \\mathcal{U}(-5, 5)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(x|z) = \\frac{1}{2\\sqrt{2\\pi(0.5)^2}}e^{-\\frac{(x-z)^2}{2(0.5)^2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "i.e $z$ is drawn from a uniform distribution between -5 and 5, and $p(x|z)$ is a gaussian with mean $z$ and variance 0.5. Let's assume we only have 1 training sample $x=0$. We'll see how efficiently we can approximate $p(x=0|z)$ by sampling $z$ from the prior and averaging $p(x=0|z)$. Note we're not trying to optimize $\\theta$, we're trying to approximate $E_{z \\sim p(z)} p(x=0|z)$, which is the inner loop of our optimization objective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean with 100 samples: 0.087408\n",
      "Mean with 100000 samples: 0.099855\n",
      "Standard deviation of 1000 means using 100 samples each: 0.020751\n"
     ]
    }
   ],
   "source": [
    "z = np.random.rand(1000,100)*10-5 # z: 1000x100 uniform randomly distributed numbers between -5 and 5.\n",
    "px = norm.pdf(z, scale=0.5) # px: p(x=0|z) for all of the 1000x100 z's above.\n",
    "\n",
    "print \"Mean with %d samples: %f\"%(z.shape[1], px.mean(axis=1)[0])\n",
    "print \"Mean with %d samples: %f\"%(np.product(z.shape), px.mean())\n",
    "print \"Standard deviation of %d means using %d samples each: %f\"%(z.shape[0], z.shape[1], px.mean(axis=1).std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot what's happening"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHuCAYAAABH3t/lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QVOWd7/HPl0FFdGBAUyiQDCO4pXcuEk1FSS2babBW\nUUmk8kcKTLDWixusrN5UjDduxCtDYvyRyuancgm7lLcImaCuWyurUfS6jlupmOiNYkQHMC6OgHO1\nVn6IbJlhmO/9o2eG+dE90z1z+jznnH6/qqaq+/SZ098+M92ffp5znueYuwsAAMRrXOgCAACoRgQw\nAAABEMAAAARAAAMAEAABDABAAAQwAAABjBjAZrbRzN41sz8Ms85PzOwNM9tuZp+MtkQAALKnlBbw\nA5IuL/agmV0haba7nytplaT1EdUGAEBmjRjA7v5rSQeHWeVqSZt61v2dpMlmNi2a8gAAyKYojgHP\nkLS33/39PcsAAEAR4+N8MjNj3ksAQNVxdxu8LIoW8H5JH+93f2bPsmJFpOpnzZo1wWvI+g/7mH2c\nlR/2M/u40E8xpQaw9fwUslXStZJkZvMlHXL3d0vcLgAAVWnELmgza5GUk3SGmb0taY2kkyW5u29w\n91+Z2ZVm9kdJRyVdV8mCAQDIghED2N2vKWGdG6MpJ3lyuVzoEjKPfVx57ON4sJ8rL0v72Ibrn478\nycw8zucDACA0M5MXOAkr1rOgAQCVNWvWLLW3t4cuoyrV19frrbfeKnl9WsAAkCE9ra3QZVSlYvu+\nWAuYizEAABAAAQwAQAAEMAAAARDAAAAEQAADAILr7OxUY2Oj3n13+IkUH3vsMS1btiymqiqLAAYA\nBLdhwwY1NTVp2rThr2a7ZMkSvf7669qxY0dMlVUOAQwA0L590pe+JP3FX0jf+Y7U1RXv869fv14r\nVqwoad1ly5bpZz/7WYUrqjwCGACqwLZt0he+IC1fLr300sDHDh6UPvUp6cEHpV//WrrnHumv/3ro\nNrq6pD17pMOHR1dDQ0OD7rnnHjU2NuqMM87QypUr1dnZqbffflt79uzRJZdcIknq6OhQbW2tJk2a\npEmTJum0005TTU1N33ZyuZwef/zx0RWRIEzEAQAZUmgyiH/+53zr9j//M39/4sR80F54Yf5+S4u0\napX04YcnfqemRvroI2l8z3yJu3dLixblw7qrS7rjDmn16vJqa2hoUG1trZ588klNnDhRS5Ys0aJF\nizR//nzdeuutevXVVwv+3pe//GVJ0ubNmyVJBw8e1JlnnqnDhw/r9NNPL6+ICmIiDgDAAN/5zonw\nlfK3f/Sj8rZx9dXSO+/kf7ezU7r7bunf/q38Wm666SZNnz5ddXV1Wr16tVpaWnTo0CHV1tYWXP/e\ne+/Vrl27tHHjxr5ltbW1cncdOnSo/AIShAAGgIw7fnzosv7HeK+4It8q7u3lnThRWrHiROvXPd8C\n7t+46+qSXn65/FpmzpzZd7u+vl4dHR2aOnWqPvjggyHrPvHEE/rpT3+qRx99VKecckrf8iNHjsjM\nVFdXV34BCUIAA0DG3XRTPlR7nXpqvsu515Qp0u9/Ly1bJi1YIH3rW9Lf//2Jx82ks84auM2TTpLO\nOaf8Wvbu3dt3u729XdOnT9fcuXO1Z88edXd39z22a9cuXXfddXr44Yc1ffr0Adtoa2vTrFmzEtX9\nPBpcDQkAMm7lSmncOGndOumUU/LHbz/72YHrzJwp9RxiLeihh6TFi/Ot5GPH8l3SS5aUX8v999+v\nq666SqeeeqruuusuLVu2TDNmzNC5556rF154QfPnz9eRI0e0dOlSffe739VnPvOZIdt47rnndMUV\nV5T/5AnDSVgAkCGVvBrSe+9J27dLZ56ZP4HLhpxWNLyGhgbdcMMN2rRpkzo6OrR06VKtW7dOEyZM\n0Lp167Rjxw6tW7dOzz33nBYtWqTTTjtNkuTuMrO+buoLLrhAv/jFLzR37tyoX+KYlHsSFgEMABmS\n5MsRNjQ0aOPGjVq0aNGQxzo7O3XRRRfpmWeeGXYyjscee0ybN2/Wli1bKlnqqJQbwHRBAwCCO/nk\nk0ua3WrJkiVaMpq+7wTiJCwAQCys3D7rjKMLGgAyJMld0FnHRBwAAKQAAQwAQAAEMAAAARDAAAAE\nQAADABAAAQwACK6zs1ONjY169913Y3vOtWvX6tvf/vaI691yyy1av3595M9PAAMAgtuwYYOampqG\nnQVrNH74wx/q7LPPVl1dna6//nodO3as7G3ccsstuuuuu9TV/xJSESCAAQDBrV+/XitWrIh0m9u2\nbdP3vvc9Pfvss2pvb9ebb76pNWvWlL2ds846S+eff762bt0aaX0EMADghObmim26oaFB99xzjxob\nG3XGGWdo5cqV6uzs1Ntvv609e/bokksukSQdO3ZMF154oe677z5JUnd3txYsWKA777yzrOfbtGmT\nVq5cqfPOO0+TJ0/WHXfcoQceeKDgug899JBqa2s1adIkTZo0SRMmTBgwZ3VTU5Mef/zxUb7ywghg\nAKgmIwXs2rVj38YwWlpa9PTTT+vNN9/Url27dOedd2rHjh0655xzNG5cPpJOOukkbd68WWvWrNHO\nnTt19913q7u7W6tXr5Yk/fKXv9SUKVM0depUTZkyZcDtqVOnat++fZKk1157TfPmzet77nnz5um9\n997TwYMHh9T1xS9+UUeOHNEHH3yg/fv3a/bs2brmmmv6Hj///PP1yiuvjPp1F0IAA0A1KRawzc0n\nri9olv8pFrSlhHQRN910k6ZPn666ujqtXr1aLS0tOnTokGprawes19jYqNtvv11Lly7VD37wA23e\nvLlvLunly5fr4MGDOnDggA4ePDjg9oEDBzRz5kxJ0ocffqjJkyf3bXPSpElydx05cqRofe6u5cuX\na+HChbr++uv7ltfW1urQoUOjft2FEMAAUA0GB+zgcG1ulnrnMXbP/1SgO7o3HCWpvr5eHR0dmjp1\nat+1fvu79tpr1d7eriuvvFLnnHNO2c91+umnD9ju4cOHZWZDwr6/2267TUePHtWPf/zjAcuPHDmi\nurq6smsYDgEMANVgcMAWC9diJyn1BvhwIV6CvXv39t1ub2/X9OnTNXfuXO3Zs0fd3d0D1v3qV7+q\nz33uc9q2bZt+85vf9C1vaWkZcLy296d3WW8XdGNj44Bu4+3bt2vatGmaMmVKwdq2bNmiBx98UI88\n8ohqamoGPNbW1jagOzsS7h7bT/7pAACVMuLn7Jo1UTzJqH5t1qxZfsEFF/i+ffv8/fff9wULFvjt\nt9/u7u7z5s3z559/vm/dTZs2+Zw5c/zo0aPe0tLis2fP9qNHj5b1fE8++aSfffbZ/vrrr/uBAwc8\nl8v5bbfd1vd4c3Ozr1271t3dX3rpJf/Yxz7mr7zySsFtXXbZZf7www8P+3zF9n3P8iGZSAsYAKpJ\nFN3KoxjK0+uaa67RZZddpjlz5ujcc8/tO7HqK1/5ijZt2iQp30q++eab9fOf/1wTJ07U8uXL9elP\nf1pf//rXy3quyy+/XN/85je1cOFCNTQ0aPbs2Wou8vq3bt2qQ4cOacGCBX0t6auuukqS1NHRoba2\nNi1dunTUr7sQrgcMABmS5OsBNzQ0aOPGjQOG9/Tq7OzURRddpGeeeSbyyTiKWbt2rcxMd9xxx7Dr\n3XLLLZozZ45uuOGGYdcr93rA48usFwCAyJ188snasWNH6DIK+v73v1+R7RLAAIBY9A4jSoqFCxcG\nfX66oAEgQ5LcBZ115XZBcxIWAAABEMAAAARAAAMAEAAnYQFAhtTX1yfuZKdqUV9fX9b6nIQFAEAF\ncRIWAAAJQgADABAAAQwAQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAADABAAAQwA\nQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAAD\nABAAAQwAQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDA\nAAAEQAADABAAAQwAQAAlBbCZLTaznWa228xuLfD4JDPbambbzexVM/uryCsFACBDzN2HX8FsnKTd\nki6V9I6kFyUtc/ed/db5lqRJ7v4tMztT0i5J09y9a9C2fKTnAwAgS8xM7m6Dl5fSAr5Y0hvu3u7u\nxyRtkXT1oHVcUm3P7VpJ7w8OXwAAcEIpATxD0t5+9/f1LOvvPkn/xczekfSKpK9FUx4AANkU1UlY\nl0t62d2nS7pQ0v1mdnpE2wYwFmYDfwAkwvgS1tkv6RP97s/sWdbfdZLuliR3f9PM9kg6T9L/Hbyx\n5ubmvtu5XE65XK6sggGUoVDgmkmciwFUTGtrq1pbW0dcr5STsGqUP6nqUkkdkl6QtNzd2/qtc7+k\n99x9rZlNUz5457n7gUHb4iQsIC4jtXZ5LwKxKHYS1ogtYHc/bmY3SnpK+S7rje7eZmar8g/7Bkl3\nSvrfZvaHnl/75uDwBRCjUrqaaQkDQY3YAo70yWgBA/Eo9Vgv70eg4sYyDAlAmpRzohUnZQHBEMBA\nNXCntQskDAEMZAktWiA1CGAg60Zq+RLaQBAEMFBN6IYGEoMABrKCliyQKgQwkGWltngJbyB2BDBQ\nbeiGBhKBAAYAIAACGMgCupCB1CGAgawqt6uZEAdiRQAD1YjjwEBwBDAAAAEQwEDa0XUMpBIBDGTR\naLuYCXMgNgQwUK04DgwERQADABAAAQwAQAAEMJBmHLMFUosABrJmrMd2x/GxAMSBdxpQzQqFNSdn\nAbEggAEACIAABgAgAAIYAIAACGAgrTgDGkg1AhjIkqhOoCLcgYojgIFqx1nPQBAEMAAAARDAAAAE\nQAADABAAAQykESdJAalHAANZEfXJVIQ8UFEEMADOhAYCIIABAAiAAAYAIAACGACAAAhgAAACIICB\ntCl0dnJNTfx1ABgTAhjIgq6uymyXoUhAxRDAAPIYigTEigAGACAAAhgAgAAIYAAAAiCAAQAIgAAG\n0oSzkoHMIICBtKv02cuEPlARBDCAExiKBMSGAAYAIAACGACAAAhgAAACIIABAAiAAAbSYtas0BUA\niJB5jGc9mpnH+XxAphQaDlSJ91NczwNUCTOTuw95Y9ECBjAyxgIDkSOAAQxEaxeIBQEMAEAABDAA\nAAEQwAAABEAAAwAQAAEMpEFzc+gKAESMccBAGsQ9NpexwEBkGAcMYGwYCwxEigAGMBStXaDiCGAA\nAAIggAEACIAABgAgAAIYAIAACGAAAAIggIGkY/gPkEkEMJBGDBMCUo8ABlBYoZCnNQ5EhgAGACAA\nAhgAgAAIYAAAAiCAAQAIgAAGACAAAhhIslwudAUAKsQ8xvGEZuZxPh+QeoWG/cT5Hgr9/EAGmJnc\nfcibiRYwgPKM42MDiALvJADFFWrt0gIGIkEAAwAQQEkBbGaLzWynme02s1uLrJMzs5fNbIeZPRtt\nmQAAZMuIJ2GZ2ThJuyVdKukdSS9KWubuO/utM1nSbyRd5u77zexMd/+PAtviJCygHEk4CSoJNQAp\nNpaTsC6W9Ia7t7v7MUlbJF09aJ1rJD3i7vslqVD4AgCAE0oJ4BmS9va7v69nWX9/JmmqmT1rZi+a\n2YqoCgQAIIvGR7idiyQtknSapOfN7Hl3/+PgFZubm/tu53I55ZhoACis33sFQHq0traqtbV1xPVK\nOQY8X1Kzuy/uuf+3ktzd7+23zq2SJrj72p77/yDpCXd/ZNC2OAYMlCopx16TUgeQUmM5BvyipDlm\nVm9mJ0taJmnroHUelbTAzGrMbKKkSyS1jbVoAAk1YULoCoDUG7EL2t2Pm9mNkp5SPrA3unubma3K\nP+wb3H2nmW2T9AdJxyVtcPfXK1o5gHi4D20F/+lPYWoBMoS5oIGkSlLXb5JqAVKGuaABAEgQAhgA\ngAAIYAAAAiCAAQAIgAAGACAAAhhIokJnHTc1xV8HgIphGBKQREkb9lOonvp66a23Yi8FSJtiw5AI\nYCCJkhbAUjJrAlKAccAAACQIAQwAQAAEMAAAARDAAAAEQAADABAAAQwAQAAEMJA0hYb7FFoGINUI\nYCANurtDVwAgYgQwgNLQMgciRQADKA2tcCBSBDAAAAEQwAAABEAAAwAQAAEMAEAABDAAAAEQwECS\n1NWFrgBATMxjvKC2mXmczwekTtIvep/0+oAEMjO5+5A3Dy1gAGNDqx0YFVrAQJKkoYWZhhqBBKEF\nDABAghDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAADScHF7YGqQgADScb4WiCzCGAAY0frHSgb\nAQygPLTKgUgQwAAABEAAAwAQAAEMAEAABDAAAAEQwAAABEAAAwAQAAEMJEGhcbSTJ8dfB4DYmMc4\nps/MPM7nA1KjUAAn+b2StnqBgMxM7j7kTUMLGED51qwZuozZsICy0AIGkiCNLco01gwEQAsYAIAE\nIYABAAiAAAYAIAACGACAAAhgAAACIICB0HK50BUACIBhSEBoaR3Ok9a6gZgxDAlA5Y0fH7oCIDUI\nYACjU6i1e/x4/HUAKUUAAwAQAAEMAEAABDAAAAEQwAAABEAAAwAQAAEMAEAABDAQEhexB6oWAQwk\nDbNJAVWBAAYQLVr1QEkIYACjR2sdGDUCGACAAAhgAAACIIABAAiAAAYAIAACGACAAAhgAAACIICB\nUAqNl508Of46AARhHuM4PjPzOJ8PSLRCAZzG90dWXgdQIWYmdx/yRqEFDGBsmpqGLmM2LGBEtICB\nULLUcszSawEiRgsYAIAEIYABAAiAAAYAIAACGACAAEoKYDNbbGY7zWy3md06zHqfNrNjZvaF6EoE\nACB7RgxgMxsn6T5Jl0tqlLTczM4rst49krZFXSSQOePofAKqXSmfAhdLesPd2939mKQtkq4usN5N\nkv5R0nsR1gdkU6EhOgzbAapKKQE8Q9Lefvf39SzrY2bTJS119/8liRH4QLU55ZShy8aPj78OIEWi\neof8SFL/Y8NFQ7i5ubnvdi6XUy6Xi6gEAMF89NHQyTiOHw9TCxBYa2urWltbR1xvxJmwzGy+pGZ3\nX9xz/28lubvf22+df++9KelMSUclfcXdtw7aFjNhAVI2Z47K4msCIlBsJqxSArhG0i5Jl0rqkPSC\npOXu3lZk/Qck/Yu7/1OBxwhgQMpmWGXxNQERKBbAI3ZBu/txM7tR0lPKHzPe6O5tZrYq/7BvGPwr\nkVQMAECGcTEGIIQsthaz+JqACHAxBgAAEoQABgAgAAIYiBtdtQBEAAOoJMb5A0VxEhYQtyy3gLP8\n2oBR4iQsAAAShAAGACAAAhgAgAAIYAAAAiCAAQAIgAAGACAAAhiIU6FhOjU18dcBIDgCGAitqyt0\nBQACIIABRGfNmqHLCrX6ATATFhCrapgpqhpeI1AGZsICACBBCGAAAAIggAEACIAABgAgAAIYAIAA\nCGAgLoXODq6vj78OAInAMCQgLtUyPKdaXidQIoYhAYhHU9PQZUzGAQxBCxiISzW1DKvptQIjoAUM\nAECCEMAAAARAAAMAEAABDABAAAQwAAABEMBAHAqdFVzo2rkAqgbDkIA4VNuwnEKvt6ZG6uqKvxYg\nMIYhAYhPoS8Xx4/HXweQYAQwAAABEMAAAARAAAMAEAABDABAAAQwAAABEMBApXEpPgAFEMBACFke\nAwygJAQwgPjQGwD0IYABVAatfGBYBDAAAAEQwAAABEAAAwAQAAEMAEAABDAAAAEQwEAlFRp209QU\nfx0AEsc8xqECZuZxPh8QXKEArqb3QLFxv9W0D1D1zEzuPuTNQAsYQOUQtEBRBDAAAAEQwAAABEAA\nAwAQAAEMAEAABDAAAAEQwEClFBqCU18ffx0AEolxwEClVPsY4F6MBUaVYxwwgDAIWqAgAhgAgAAI\nYAAAAiCAAQAIgAAGACAAAhiohFmzQlcAIOEYhgRUAkOQBmJ/oIoxDAlAODU1Q5cVGx8MVAkCGEDl\ndXWFrgBIHAIYAIAACGAAAAIggAEACIAABgAgAAIYiBpn9wIoAQEMxIExr4XV1YWuAAiGiTiAqDHp\nRHHsG1QhJuIAACBBCGAAAAIggAEACIAABgAgAAIYiBJDkACUiAAGKo2zfAEUQAADiE+hLyP0GqBK\nEcAAAARQUgCb2WIz22lmu83s1gKPX2Nmr/T8/NrM5kZfKgAA2THiTFhmNk7SbkmXSnpH0ouSlrn7\nzn7rzJfU5u6HzWyxpGZ3n19gW8yEhWxjpqeRsY9QZcYyE9bFkt5w93Z3PyZpi6Sr+6/g7r9198M9\nd38racZYCwYAIMtKCeAZkvb2u79Pwwfs9ZKeGEtRQCpxMhGAMoyPcmNmtlDSdZIWRLldILXoWi3N\nuHFSd3foKoBYlRLA+yV9ot/9mT3LBjCzCyRtkLTY3Q8W21hzc3Pf7Vwup1wuV2KpADLBfWhvAV9U\nkCGtra1qbW0dcb1STsKqkbRL+ZOwOiS9IGm5u7f1W+cTkp6RtMLdfzvMtjgJC9nFyUWlY1+hihQ7\nCWvEFrC7HzezGyU9pfwx443u3mZmq/IP+wZJ/1PSVEnrzMwkHXP3i6N9CQAAZMeILeBIn4wWMLKM\nVl3p2FeoImMZhgQAACJGAANRKNSiq6+Pvw4AqUEXNBAFulTLU2zMNPsMGUQXNIDkIGgBAhgAgBAI\nYAAAAiCAAQAIgAAGxoqLMAAYBQIYqAROMhodvsygihDAAMLgSwqqHAEMAEAABDAAAAEQwAAABEAA\nA2NR6KQhTiQCUAICGIhad3foCtKNLzCoEgQwgHA4ExpVjAAGACAAAhgAgAAIYGC0OFYJYAwIYCBK\nHNOMxqxZoSsAKs48xg8MM/M4nw+oqEItYP6/R4d9iQwzM7n7kH9yWsAAAARAAAMAEAABDABAAAQw\nMBqFjlmeckr8dQBILQIYiMpHH4WuIL0KfXlhmBcyjgAGEB5fXlCFCGAAAAIggAEACIAABsrFNYAB\nRIAABqLANYDHrr5+6DK+2CDDmIoSKBfTJlYO+xYZxFSUAAAkCAEMlIMuUQARIYCBsaKLtLL40oOM\nIoABJAdfZlBFCGAAAAIYH/cT0puENOuW1P9f2CWN4386UuxjVAtawECJumUiBwBEhQAGRskljRPH\nLKP2nJoG3Dflv/wAWRP7RBziAwspNbgFTABXDvsa2cJEHAAAJAYBDJSA478AokYAA6NAl2hltWvg\nhRk4Dows4mIMQCm4SED82OfICC7GAABAghDAwEiYPSY5+FsgQwhgYDToCq089jEyjgAGACAAAhgY\nDl2eACqEAAbKRddoWHwpQkYQwACSiy87yDACGCimuTl0BQAyjIk4gGKYCCIZinU587dASjARB4B0\nImiRUQQwAAABEMBAIYW6PZuahi5DOJwNjZTjGDBQCMd/k4e/CVKKY8BAqWhZAYgBAQyUgpZWMvFl\nCSlGAANIB74EIWMIYKA/WlTpM2FC6AqAUSGAgZHQ8kqOQn+LP/0p/jqACBDAQC9avwBiRAADw6H1\nmw58eUIKEcAA0oUvRcgIAhiQmOQhC2gFI2UIYADpw5cjZAABDBRqOa1ZE38dGDtawUgR5oJGdeNa\ns+nGoQOkAHNBA6XiAzzdaAUjJQhgVC8+qNOv2Jcl/rZIAQIY1am5ufByWr/pw98MKcUxYFQnjh1m\nC8fykWAcAwZ60T2ZPXRFI4UIYFQXWkrZRQgjZQhgVA/CN/sIYaQIAYzqQPhWD0IYKUEAI9vMCN9q\nNFwIE8RICAIY2TTSBy3hm33D/Y0JYiQAAYxsKeWDlfCtHiP9rQliBFRSAJvZYjPbaWa7zezWIuv8\nxMzeMLPtZvbJaMsMp7W1NXQJmTemfdz7AVrqB2mVhm9V/x+X8jcv9/+oiKrezzHJ0j4eMYDNbJyk\n+yRdLqlR0nIzO2/QOldImu3u50paJWl9BWoNIkt/7KQquI8HfyAW+ylHlYavxP9x2X/7Uf7vVf1+\njkGW9vH4Eta5WNIb7t4uSWa2RdLVknb2W+dqSZskyd1/Z2aTzWyau787ZGtp7O5ZuzZ0BdlXyX1c\nxcGLfnr/D6L6DCq2HT4vKi8j+7iULugZkvb2u7+vZ9lw6+wvsA4QL3fCF0Pxf4GEKKUFHKkUtn+V\nje9ayVaRfZzG3pYKWpuRVkPSsZcrLyv7uJQA3i/pE/3uz+xZNnidj4+wTsHJqAEAqEaldEG/KGmO\nmdWb2cmSlknaOmidrZKulSQzmy/pUMHjvwAAQFIJLWB3P25mN0p6SvnA3ujubWa2Kv+wb3D3X5nZ\nlWb2R0lHJV1X2bIBAEi3WK8HDAAA8pgJqwxm9g0z6zazqaFryRoz+56ZtfVM5PKImU0KXVNWlDKR\nDkbPzGaa2b+a2Wtm9qqZ/ffQNWWVmY0zs5fMbPBh0FQigEtkZjMl/aWk9tC1ZNRTkhrd/ZOS3pD0\nrcD1ZEIpE+lgzLok3ezujZI+I+lv2McV8zVJr4cuIioEcOl+KOl/hC4iq9z9/7h7d8/d3yp/Jj3G\nrm8iHXc/Jql3Ih1ExN3/n7tv77n9oaQ2MQ9C5HoaQVdK+ofQtUSFAC6BmX1e0l53fzV0LVXiv0l6\nInQRGVHKRDqIiJnNkvRJSb8LW0km9TaCMnPiUuwTcSSVmT0taVr/Rcr/oW+XdJvy3c/9H0OZhtnH\nq939X3rWWS3pmLu3BCgRGDUzO13SP0r6Wk9LGBExs6skvevu280sp4x8BhPAPdz9LwstN7P/KmmW\npFfMzJTvGv29mV3s7u/FWGLqFdvHvczsr5TvYloUS0HVoZSJdDBGZjZe+fD9ubs/GrqeDPpzSZ83\nsyslnSqp1sw2ufu1gesaE4YhlcnM9ki6yN0Phq4lS8xssaS/k/RZd38/dD1ZYWY1knZJulRSh6QX\nJC1397aghWWMmW2S9B/ufnPoWrLOzJokfcPdPx+6lrHiGHD5XBnp/kiYn0o6XdLTPcMM1oUuKAvc\n/bik3ol0XpO0hfCNlpn9uaQvSVpkZi/3/P8uDl0Xko8WMAAAAdACBgAgAAIYAIAACGAAAAIggAEA\nCIAABgATZ0NiAAAAqUlEQVQgAAIYAIAACGAAAAIggAEACIAABjLMzFb1m53p383smdA1AchjJiyg\nCvRcLOAZSfe6+69C1wOAFjBQLX4i6V8JXyA5uBwhkHE9l3n8uLt/NXQtAE4ggIEMM7NPSfqGpAWh\nawEwEF3QQLb9jaQpkp7tORFrQ+iCAORxEhYAAAHQAgYAIAACGACAAAhgAAACIIABAAiAAAYAIAAC\nGACAAAhgAAAC+P8i9F9W8DJb5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104721e90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(z, np.ones(z.shape)/10., color='b', label='p(z)')\n",
    "plt.scatter(z, px, marker='+', color='r', label='p(x=0|z)')\n",
    "plt.xlabel('z')\n",
    "plt.ylim([0,1])\n",
    "plt.xlim([-5,5])\n",
    "_=plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most drawn z's will result in $p(x|z)$ very close to 0 and thus not contribute meaningfully to the expectation of $p(x|z)$. Just eyeballing from the plot we can see that picking z at random from -5 to 5 will hit a nonzero region of p(x=0|z) about 40% of the time. So if we only drew 10 z's there'd be a good chance we would not hit a nonzero region at all, and would estimate the expectation to be 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can we do better? If we could draw $z$ from $p(z|x=0)$ instead all z's drawn would contribute meaningfully to the expectation. However, if we did that, then the we'd no longer be computing $p(x=0)$, because we'd only be drawing z's that have high probability of causing $x=0$. So we'll have to re-weight $p(x=0|z)$ with how likely it was to draw $z$ from $p(z|x=0)$ compared to the prior $p(z)$. This is also known as importance sampling.\n",
    "\n",
    "In general we don't know $p(z|x)$, so we'll draw from something else $q(z|x)$. In this case, let's pick $q(z|x) = \\mathcal {N}(0,1)$, i.e. a normal gauss. Let's see how that looks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHuCAYAAABH3t/lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+QVfV9//HXmx+CPxYWtGP4kSwgfmu+BGnMGM2EhpVM\n/J1KmqYDttjw1apjYzqmfmMV4+4ak0imadKoDCVh7KChJNZMtbHG5Etd2zRa861iQQGtwRWQr04D\nCLEju8D7+8fdXXb3nrv317nnc865z8fMzuw99+w5n733nPM678/5Ze4uAACQrDGhGwAAQDMigAEA\nCIAABgAgAAIYAIAACGAAAAIggAEACKBsAJvZOjN708z+Y5Rxvm1mr5jZZjP7rXibCABA/lRSAd8v\n6aJSb5rZJZLOcPczJV0naU1MbQMAILfKBrC7/0zS/lFGuULS+v5x/03SZDM7PZ7mAQCQT3EcA54h\nadeQ13v6hwEAgBLGJTkzM+O+lwCApuPuNnJYHBXwHknvHfJ6Zv+wUo3I1E9HR0fwNuT9h8+Yzzgv\nP3zOfMZRP6VUGsDW/xPlUUlXSZKZnS/pgLu/WeF0AQBoSmW7oM1sg6R2Saea2euSOiSdIMndfa27\n/6OZXWpm/ynpHUkrGtlgAADyoGwAu/uVFYzzuXiakz7t7e2hm5B7fMaNx2ecDD7nxsvTZ2yj9U/H\nPjMzT3J+AACEZmbyiJOwEj0LGgDQWLNmzVJPT0/oZjSltrY2vfbaaxWPTwUMADnSX22FbkZTKvXZ\nl6qAeRgDAAABEMAAAARAAAMAEAABDABAAAQwACC43t5ezZs3T2++OfqNFH/0ox9p6dKlCbWqsQhg\nAEBwa9eu1aJFi3T66aM/zfbyyy/XSy+9pK1btybUssYhgAEA2r1b+oM/kH77t6Uvf1k6ciTZ+a9Z\ns0bLly+vaNylS5fqr//6rxvcosYjgAGgCTzxhPS7vystWyY999zw9/bvlz70Ien735d+9jPp7rul\nP/7j4mkcOSLt3Cm9/XZtbZg9e7buvvtuzZs3T6eeeqquvvpq9fb26vXXX9fOnTt13nnnSZL27t2r\nlpYWTZo0SZMmTdLJJ5+ssWPHDk6nvb1djz32WG2NSBFuxAEAORJ1M4i///tCdfvf/114fdJJhaD9\n4AcLrzdskK67Tvr1r4//zdix0rvvSuP675f48svS4sWFsD5yRLrjDmnlyuraNnv2bLW0tOjHP/6x\nTjrpJF1++eVavHixzj//fN1yyy3asmVL5N/94R/+oSTpwQcflCTt379fp512mt5++22dcsop1TWi\ngbgRBwBgmC9/+Xj4SoXfv/Wt6qZxxRXSG28U/ra3V/ra16R//ufq23LjjTdq+vTpam1t1cqVK7Vh\nwwYdOHBALS0tkeOvWrVKO3bs0Lp16waHtbS0yN114MCB6huQIgQwAOTc0aPFw4Ye473kkkJVPNDL\ne9JJ0vLlx6tf90IFPLS4O3JEev756tsyc+bMwd/b2tq0d+9eTZ06VQcPHiwa9/HHH9c999yjRx55\nRBMmTBgcfujQIZmZWltbq29AihDAAJBzN95YCNUBJ55Y6HIeMGWK9O//Li1dKi1cKN16q/Sd7xx/\n30x6z3uGT3P8eGnOnOrbsmvXrsHfe3p6NH36dM2fP187d+7UsWPHBt/bsWOHVqxYoYceekjTp08f\nNo1t27Zp1qxZqep+rgVPQwKAnLv6amnMGGn1amnChMLx2499bPg4M2dK/YdYI/3gB9LFFxeq5L6+\nQpf05ZdX35b77rtPl112mU488UR99atf1dKlSzVjxgydeeaZevbZZ3X++efr0KFDWrJkib7yla/o\nIx/5SNE0nnrqKV1yySXVzzxlOAkLAHKkkU9DeustafNm6bTTCidwWdFpRaObPXu2rr/+eq1fv157\n9+7VkiVLtHr1ak2cOFGrV6/W1q1btXr1aj311FNavHixTj75ZEmSu8vMBrupzz77bH3ve9/T/Pnz\n4/4X61LtSVgEMADkSJofRzh79mytW7dOixcvLnqvt7dX55xzjjZt2jTqzTh+9KMf6cEHH9TGjRsb\n2dSaVBvAdEEDAII74YQTKrq71eWXX67La+n7TiFOwgIAJMKq7bPOObqgASBH0twFnXfciAMAgAwg\ngAEACIAABgAgAAIYAIAACGAAAAIggAEAwfX29mrevHl68803E5tnV1eX7rzzzrLj3XzzzVqzZk3s\n8yeAAQDBrV27VosWLRr1Lli1+OY3v6lp06aptbVV11xzjfr6+qqexs0336yvfvWrOjL0EVIxIIAB\nAMGtWbNGy5cvj3WaTzzxhL7+9a/rySefVE9Pj1599VV1dHRUPZ33vOc9ev/7369HH3001vYRwACA\n4zo7Gzbp2bNn6+6779a8efN06qmn6uqrr1Zvb69ef/117dy5U+edd54kqa+vTx/84Ad17733SpKO\nHTumhQsX6q677qpqfuvXr9fVV1+ts846S5MnT9Ydd9yh+++/P3LcH/zgB2ppadGkSZM0adIkTZw4\ncdg9qxctWqTHHnusxv88GgEMAM2kXMB2ddU/jVFs2LBBP/3pT/Xqq69qx44duuuuu7R161bNmTNH\nY8YUImn8+PF68MEH1dHRoe3bt+trX/uajh07ppUrV0qS/vZv/1ZTpkzR1KlTNWXKlGG/T506Vbt3\n75Ykvfjii1qwYMHgvBcsWKC33npL+/fvL2rX7//+7+vQoUM6ePCg9uzZozPOOENXXnnl4Pvvf//7\n9cILL9T8f0chgAGgmZQK2M7O488XNCv8lAraSkK6hBtvvFHTp09Xa2urVq5cqQ0bNujAgQNqaWkZ\nNt68efN0++23a8mSJfrLv/xLPfjgg4P3kl62bJn279+vffv2af/+/cN+37dvn2bOnClJ+vWvf63J\nkycPTnPSpElydx06dKhk+9xdy5Yt0wUXXKBrrrlmcHhLS4sOHDhQ8/8dhQAGgGYwMmBHhmtnpzRw\nH2P3wk8DuqMHwlGS2tratHfvXk2dOnXwWb9DXXXVVerp6dGll16qOXPmVD2vU045Zdh03377bZlZ\nUdgPddttt+mdd97RX/3VXw0bfujQIbW2tlbdhtEQwADQDEYGbKlwLXWS0kCAjxbiFdi1a9fg7z09\nPZo+fbrmz5+vnTt36tixY8PGveGGG/TJT35STzzxhH7+858PDt+wYcOw47UDPwPDBrqg582bN6zb\nePPmzTr99NM1ZcqUyLZt3LhR3//+9/Xwww9r7Nixw97btm3bsO7sWLh7Yj+F2QEAGqXsdrajI46Z\n1PRns2bN8rPPPtt3797tv/rVr3zhwoV+++23u7v7ggUL/Omnnx4cd/369T537lx/5513fMOGDX7G\nGWf4O++8U9X8fvzjH/u0adP8pZde8n379nl7e7vfdtttg+93dnZ6V1eXu7s/99xz/hu/8Rv+wgsv\nRE7rwgsv9IceemjU+ZX67PuHF2UiFTAANJM4upVruJRnwJVXXqkLL7xQc+fO1Zlnnjl4YtW1116r\n9evXSypUyV/4whf0wAMP6KSTTtKyZct07rnn6qabbqpqXhdddJG++MUv6oILLtDs2bN1xhlnqLPE\n///oo4/qwIEDWrhw4WAlfdlll0mS9u7dq23btmnJkiU1/99ReB4wAORImp8HPHv2bK1bt27Y5T0D\nent7dc4552jTpk2x34yjlK6uLpmZ7rjjjlHHu/nmmzV37lxdf/31o45X7fOAx1XZXgAAYnfCCSdo\n69atoZsR6S/+4i8aMl0CGACQiIHLiNLiggsuCDp/uqABIEfS3AWdd9V2QXMSFgAAARDAAAAEQAAD\nABAAAQwAQAAEMAAAARDAAIBU6O3t1bx58/Tmm2+WHXfg0YXlbNmyRR/96EfrbVpDEMAAgFRYu3at\nFi1aVNGdsCq9pnj+/PmaMmWKHnvssXqbFzsCGACQCmvWrNHy5ctjn+6VV16pNWvWxD7dehHAAIBB\nnd2dDZ3+888/rw996EOaPHmyli5dqmXLlumOO+7Qrl27tHPnTp133nmSCg9AGPrIwZNPPrnoEYED\nbrjhBv3e7/3e4OtbbrlFn/jEJwZft7e3a9OmTerr62vo/1YtAhgAmki5gO16qqvuaZTS19enT33q\nU/qjP/oj7du3T5/5zGf08MMPSyocq50zZ87gsd1p06bp0KFDOnjwoA4ePKhPfepTWrZsWeR0v/GN\nb2jr1q1av369/uVf/kX333//4JOVJGn69OkaP368duzYUVO7G4UABoAmUipgO7s7ZV2F46rWZbIu\nKxm0lYR0lGeeeUZHjhzR5z//eY0dO1af/vSnde6550qSDhw4oJaWlsi/W7VqlXbs2KF169ZFvn/i\niSfqgQce0E033aSrrrpK9957r6ZNmzZsnJaWFh04cKCmdjcKAQwATWBkwI4M1872TnlH4T7G3uHy\nDldne6fi9MYbb2jGjBnDhrW1tUmSpkyZokOHDhX9zeOPP6577rlHjzzyiCZMmFBy2ueee67mzJkj\nd9dnPvOZovcPHTqk1tbWOv+DeBHAANAERgZsqXDtWNQR/ff9AT5aiJczbdo07dmzZ9iw119/XZJ0\n9tln65e//KWOHTs2+N6OHTu0YsUKPfTQQ5o+ffqo077vvvvU29ur6dOna9WqVcPee+ONN9TX16ff\n/M3frKq9Defuif0UZgcAaJRy29mOJzvqn0dnbdvy3t5eb2tr829/+9ve19fnDz/8sI8fP96/9KUv\nubv7ggUL/Omnn3Z394MHD/pZZ53l3/3udyOnNWbMmMHfd+zY4VOmTPEtW7b4K6+84lOnTvXNmzcP\nvr9hwwa/7LLLampzNUp99v3DizKRChgAmkgc3cqlquRyxo8frx/+8Ie6//77deqpp+qhhx7Spz/9\n6cH3r7322sGTp5577jm9/PLLuummmzRp0qTBM6JHOnr0qJYvX65bb71VH/jABzR37lx95Stf0VVX\nXTV41vP3vvc9XX/99TW1uZF4HjAA5EjWnge8YsUKvfe979Wdd96p3t5enXPOOdq0aVPZm3GMGTNm\nWHd1KVu2bNH111+vf/3Xf42rySVV+zzgcQ1vEQAAFTjhhBO0devWisat5k5YSYRvLeiCBgAEU2mQ\njtTRUVs3eJrQBQ0AOZK1Lug8qbYLmgoYAIAACGAAAAIggAEACICzoAEgR9ra2mo+sQn1GbitZqU4\nCQsAgAbiJCwAAFKEAAYAIAACGACAAAhgAAACIIABAAiAAAYAIAACGACAAAhgAAACIIABAAiAAAYA\nIAACGACAAAhgAAACIIABAAiAAAYAIAACGACAAAhgAAACqCiAzexiM9tuZi+b2S0R708ys0fNbLOZ\nbTGzz8beUgAAcsTcffQRzMZIelnSxyW9IekXkpa6+/Yh49wqaZK732pmp0naIel0dz8yYlpebn4A\nAOSJmcndbeTwSirgD0t6xd173L1P0kZJV4wYxyW19P/eIulXI8MXAAAcV0kAz5C0a8jr3f3DhrpX\n0v80szckvSDpT+NpHgAA+RTXSVgXSXre3adL+qCk+8zslJimDaAeZsN/AKTCuArG2SPpfUNez+wf\nNtQKSV+TJHd/1cx2SjpL0v8dObHOzs7B39vb29Xe3l5VgwFUISpwzSTOxQAapru7W93d3WXHq+Qk\nrLEqnFT1cUl7JT0raZm7bxsyzn2S3nL3LjM7XYXgXeDu+0ZMi5OwgKSUq3ZZF4FElDoJq2wF7O5H\nzexzkn6iQpf1OnffZmbXFd72tZLukvQ3ZvYf/X/2xZHhCyBBlXQ1UwkDQZWtgGOdGRUwkIxKj/Wy\nPgINV89lSACypJoTrTgpCwiGAAaagTvVLpAyBDCQJ1S0QGYQwEDelat8CW0gCAIYaCZ0QwOpQQAD\neUElC2QKAQzkWaUVL+ENJI4ABpoN3dBAKhDAAAAEQAADeUAXMpA5BDCQV9V2NRPiQKIIYKAZcRwY\nCI4ABgAgAAIYyDq6joFMIoCBPKq1i5kwBxJDAAPNiuPAQFAEMAAAARDAAAAEQAADWcYxWyCzCGAg\nb+o9tjuGzQKQBNY0oJlFhTUnZwGJIIABAAiAAAYAIAACGACAAAhgIKs4AxrINAIYyJO4TqAi3IGG\nI4CBZsdZz0AQBDAAAAEQwAAABEAAAwAQAAEMZBEnSQGZRwADeRH3yVSEPNBQBDAAzoQGAiCAAQAI\ngAAGACAAAhgAgAAIYAAAAiCAgayJOjt57Njk2wGgLgQwkAdHjjRmulyKBDQMAQyggEuRgEQRwAAA\nBEAAAwAQAAEMAEAABDAAAAEQwECWcFYykBsEMJB1jT57mdAHGoIABnAclyIBiSGAAQAIgAAGACAA\nAhgAgAAIYAAAAiCAgayYNSt0CwDEyDzBsx7NzJOcH5ArUZcDNWJ9Smo+QJMwM7l70YpFBQygPK4F\nBmJHAAMYjmoXSAQBDABAAAQwAAABEMAAAARAAAMAEAABDGRBZ2foFgCIGdcBA1mQ9LW5XAsMxIbr\ngAHUh2uBgVgRwACKUe0CDUcAAwAQAAEMAEAABDAAAAEQwAAABEAAAwAQAAEMpB2X/wC5RAADWcRl\nQkDmEcAAokWFPNU4EBsCGACAAAhgAAACIIABAAiAAAYAIAACGACAAAhgIM3a20O3AECDmCd4PaGZ\neZLzAzIv6rKfJNeh0PMHcsDM5O5FKxMVMIDqjGGzAcSBNQlAaVHVLhUwEAsCGACAACoKYDO72My2\nm9nLZnZLiXHazex5M9tqZk/G20wAAPKl7ElYZjZG0suSPi7pDUm/kLTU3bcPGWeypJ9LutDd95jZ\nae7+XxHT4iQsoBppOAkqDW0AMqyek7A+LOkVd+9x9z5JGyVdMWKcKyU97O57JCkqfAEAwHGVBPAM\nSbuGvN7dP2yo/yFpqpk9aWa/MLPlcTUQAIA8GhfjdM6RtFjSyZKeNrOn3f0/R47Y2dk5+Ht7e7va\nudEAEG3IugIgO7q7u9Xd3V12vEqOAZ8vqdPdL+5//eeS3N1XDRnnFkkT3b2r//V3JT3u7g+PmBbH\ngIFKpeXYa1raAWRUPceAfyFprpm1mdkJkpZKenTEOI9IWmhmY83sJEnnSdpWb6MBpNTEiaFbAGRe\n2S5odz9qZp+T9BMVAnudu28zs+sKb/tad99uZk9I+g9JRyWtdfeXGtpyAMlwL66CDx8O0xYgR7gX\nNJBWaer6TVNbgIzhXtAAAKQIAQwAQAAEMAAAARDAAAAEQAADABAAAQykUdRZx4sWJd8OAA3DZUhA\nGqXtsp+o9rS1Sa+9lnhTgKwpdRkSAQykUdoCWEpnm4AM4DpgAABShAAGACAAAhgAgAAIYAAAAiCA\nAQAIgAAGACAAAhhIm6jLfaKGAcg0AhjIgmPHQrcAQMwIYACVoTIHYkUAA6gMVTgQKwIYAIAACGAA\nAAIggAEACIAABgAgAAIYAIAACGAgTVpbQ7cAQELME3ygtpl5kvMDMiftD71Pe/uAFDIzuXvRykMF\nDKA+VO1ATaiAgTTJQoWZhTYCKUIFDABAihDAAAAEQAADABAAAQwAQAAEMAAAARDAAAAEQAADacHD\n7YGmQgADacb1tUBuEcAA6kf1DlSNAAZQHapyIBYEMAAAARDAAAAEQAADABAAAQwAQAAEMAAAARDA\nAAAEQAADaRB1He3kycm3A0BizBO8ps/MPMn5AZkRFcBpXley1l4gIDOTuxetNFTAAKrX0VE8jLth\nAVWhAgbSIIsVZRbbDARABQwAQIoQwAAABEAAAwAQAAEMAEAABDAAAAEQwEBo7e2hWwAgAC5DAkLL\n6uU8WW03kDAuQwLQeOPGhW4BkBkEMIDaRFW7R48m3w4gowhgAAACIIABAAiAAAYAIAACGACAAAhg\nAAACIIABAAiAAAZC4iH2QNMigIG04W5SQFMggAHEi6oeqAgBDKB2VOtAzQhgAAACIIABAAiAAAYA\nIAACGACAAAhgAAACIIABAAiAAAZCibpedvLk5NsBIAjzBK/jMzNPcn5AqkUFcBbXj7z8H0CDmJnc\nvWhFoQIGUJ9Fi4qHcTcsoCwqYCCUPFWOefpfgJhRAQMAkCIEMAAAARDAAAAEQAADABBARQFsZheb\n2XYze9nMbhllvHPNrM/Mfje+JgIAkD9lA9jMxki6V9JFkuZJWmZmZ5UY725JT8TdSCB3xtD5BDS7\nSrYCH5b0irv3uHufpI2SrogY70ZJfyfprRjbB+RT1CU6XLYDNJVxFYwzQ9KuIa93qxDKg8xsuqQl\n7n6BmQ17D0DyrKuyG2F4R0yhP2GCdPjw8GHjxklHjsQzfSCHKgngSnxL0tBjwyXX/s7OzsHf29vb\n1d7eHlMTAFQavCPHrzuI3323+GYcR4/WN00go7q7u9Xd3V12vLJ3wjKz8yV1uvvF/a//XJK7+6oh\n4/xy4FdJp0l6R9K17v7oiGlxJyxAiv3OUa13t+rtw2/X0aAYQpi7YQGRSt0Jq5IAHitph6SPS9or\n6VlJy9x9W4nx75f0D+7+w4j3CGBAijWsqq16y6k5iAlgIFLNt6J096OSPifpJ5JelLTR3beZ2XVm\ndm3Un9TdWgAViTt8GzVNAMV4GAMQQgzVYrmgLFfJ1vv3xROkAgai1NwFHXMjCGBAqjusRgvPaoMz\ntmkRwEAkAhhIkzrCqlRg1nMSVSwhTAADkQhgIC1SFr6xTZ8ABiLxPGAg4xodvqNNq+YTs7jOHyiJ\nChhIWg2VYhLhG8v8qIKBIlTAQEYlHb6NnjaAAgIYyKBQAck1wkB8CGAgxUIGXuzHgwEMQwADKRWi\n67nSeRHCQP0IYCCFJt41MXJ4iK5njgcDjUEAAyl0+OjhomFpC0KqYKA+BDCQpKjLdMaOHT5KCoON\nrmggfgQwENqRI2VHSUP1m4Y2AHlCAAMpksWKclibOzoiRsje/wQkgTthAUka5U5RaTjruRJR7RzW\nRu6GBQzDnbCADEpb+ErS5AmTi4ZlsXIHQiOAgRTIUoAd+PMDkcOz9D8AaUAAAymVxup3QJrbBmQF\nAQwElqfKMU//C9BoBDAQUGd79PAsVJhZaCOQZgQwkJSIs4O7FhWPlvVgs4grkQAUI4CBQKxDUsZ7\nbLO+swCERAADKZKLQLOIKpibcQBFCGAggDxUvwNysdMABEAAAymRqyCLqoIBDEMAAwnLU/U7IFc7\nD0BCCGAgBXIZYFTBwKgIYCBBeax+B+RyJwJoIAIYCCzXwUUVDJREAANJMMt19Tsg1zsTQMwIYCCg\npgisgSp43LjQLQFShQAGEtDZrtxXvwNK7lQcPZpsQ4CUI4CBBBTd89mbpPodwLFgoAgBDCB23uFS\nE+1fALUggIEGsy5rmu7nURnPCwaGIoCBpLnkXaEb0XjeJapgYBQEMNBAVHzF+EyAAgIYSFKTVL8D\nqIKB0ghgoEFGrfScVAKaHQEMJKXJqt8BUVUw3dAAAQw0BAEjqnygDAIYSEKTVr+jYScFzY4ABtBw\nnIwFFCOAgZhR2VWOzwrNjAAGGo3uZ0lUwcBIBDAQIyq66vGZoVkRwEAjRVW/i0Y+Gql50BMAHEcA\nAzGpuJLr7m5oO1KPbmhAEgEMNFTTV3wVXgtMNzSaEQEMIFFNv1MC9COAgRhQwdWPzxDNhgAGGsQ7\nONhZCp8NQAADdaNyi8/EuyaGbgKQGAIYaAAqvPKiPqPDRw8HaAkQBgEM1GHU6tci3mtra1xjAGQK\nAQzEbNTq97XXEmtHqkXtnAy8RZc+mgQBDKCxRrkWmK56NDMCGKgRlVrj8NmiGRDAQIyo6KrHZ4Zm\nRQADNaBCA1AvAhiICZVcvNjJQd4RwEAjzJoVugWZws4LmhEBDFQpqjIzjRjW01P8hxU+GQjHUQUj\nzwhgIAbHOo6FbkK6jR1bPGzEtcBUwWg2BDBQBSqyGh05UvOfdnZ3xtcOIEXME+wWMzNPcn5A3KIC\nOLJyi7rTU7Mv+xV+JhV/xkBGmJncvWjBpgIGKkT1CyBOBDBQByqzZLDzgzwigAGkCjs1aBYEMFCB\nqiqwUZ70g9pRBSNvCGCgRlVVas1+AlYpra2Rg6mC0QwIYKAMKq+YRO2EvP128u0AUoIABmpAhRYG\nO0PIEwIYQCqxk4O8I4CBUVBxpQ/fCfKCAAaqRGWWHD5r5BkBDJRQU6XFJUgAKkQAA1WoqSLjEqS6\nTBg7oWgY3dDIAwIYiDDxromhm5BPUTsjZXoN3r393QY1BgiLAAYiHD56uGgYxyPThSoYWVdRAJvZ\nxWa23cxeNrNbIt6/0sxe6P/5mZnNj7+pAJoVOz/Io7IBbGZjJN0r6SJJ8yQtM7OzRoz2S0kfc/cF\nku6S9J24GwokhcoqO/iukGWVVMAflvSKu/e4e5+kjZKuGDqCuz/j7gP3lHtG0ox4mwmERQUWHt8B\n8qaSAJ4hadeQ17s1esBeI+nxehoFhFJXRcUlSACqEOtJWGZ2gaQVkoqOEwNZVVflxSVIlRlT+6aI\nbmhk1bgKxtkj6X1DXs/sHzaMmZ0taa2ki919f6mJdXZ2Dv7e3t6u9vb2CpsKIBfci3sLKtxR8Q4n\ncJF63d3d6u7uLjueeZkF38zGStoh6eOS9kp6VtIyd982ZJz3Sdokabm7PzPKtLzc/IBQojbsi9oW\nqfuz3RVOICIYWN6j1fFZlQpgjhEjrcxM7l604Jbt93H3o5I+J+knkl6UtNHdt5nZdWZ2bf9oX5I0\nVdJqM3vezJ6Nse1AMBWHLxJD0CIvylbAsc6MChgpFUtVRQVcuTo/K6pgZEnNFTDQrNiYpxffDfKA\nAEbTi+WknqiKrq2t/ukCyC0CGIgQS4X12mv1T6OZxHAdNWdII0sIYDS19r9pD92E5hTDsXG6oZF1\nBDCa2lM9TxUNY8OebVTByAoCGEBmsbOELCOA0bSolPKL7xZZQAADQ9RUUfEQhqCogpFVBDCaUsMr\nJG7AUZsYd2bG3VnJre6BcAhgoB+VVMJi3EmJ+u6O+tHYpg80AgGMpsPxQQBpQAADovrNK3a2kGYE\nMJoKG+T8YicKWUMAo+nVteGOOmmIs6JThZ0upBUBDMTt2LHQLci2OnZgqIKRJQQwmgaVUAoldLkW\n3z3SiAD5TgjoAAAJ4ElEQVRGU6Niyh++U2QFAYymQAUElgGkDQGMplV3pcTJVqlFFYwsIICRe4lW\nPtyCMh6zZjVkslTBSBMCGE2JCilFonZaenrqnyzfMVKOAEauUfFgJJYJpAUBjKZDZdQ8+K6RZgQw\ncotKB6WwbCANCGA0ldgqoqgzoCdMiGfaiBVVMNKKAEYuBalw3n03+XnmRdTOS4Mv86IKRmgEMJoG\nlVCKNXjnhe8eaUQAI3eobFAplhWERACjKVABgWUAaUMAI1cSqWh4BnCuUAUjFAIYuZdI5cMzgOvX\n1lY8LOYdG6pgpAkBjNygksm4114LNmuWHYRAACPXqHgwEssE0oIARi5EVTBjbWwDZkSllFdUwUga\nAYzMK7XhPHLHkWQawCMIG6sBOz1UwUgDAhi5xAY2owLvzFAFI0kEMDKNDSZqVWonrfXu1oRbgmZl\nnuAep5m5RGWCGHWYNDSDXVJX45axY7Ki2Y1hmY5V4p9xwssQmpHJ3YuqBSpgZNfIDWeDjQwG5JSp\nsGwBDUYAIz8SrlyofhvjKS0a9tpU2PlpmC465hAGAYxsiqp+D08O0hTE6wJ1h89DqmAkgABG9kSF\nr0u6+0CI1iAPqIIRAAGM7Eug65njv83Ail9SBaOBCGBkS8InXpXC8d/G6tHwBzM0/DiwJHUdK66C\nCWE0UOKXISU5P+RLqWt+E7npRtTdmFiWGyvQZx61nHFjF9TDjMuQkENsGJEEbviCRiCAkQlBN4A8\ngCE9EvguSu3UEcKIGwGM1Bt357jI4UGrX7qfGy/gZ0zPCpJAACP1jvrRomFsIBECVTDiRAAj1YJv\n8Oh+blp0RaPRCGCkVtCznkdD93NYCe4UBV/WkGsEMDKFDWKTSenODlUw4kAAI5VSsYHr7AzdAqQA\nXdFoFG7EgdRJTdczN99Ih1Jdzgl/F6lZLpE53IgDmcBGDkVSstPDMoi4EcBIPTZ8SDO6olErAhip\nkaoNWVS356JFxcMQToBLxDgejDhxDBipkLquZ47/pk+KvpPULa9INY4BI7VStzHj5hsog0oYcSCA\nEVTqwrcUqt90CrizRAijXgQwgslM+CIdUrgTxLKKehDACCK14Uv3c/ZMnBh09lHLLFUwKkEAI3Gp\nDd9SUlh5Na2o7+Lw4eTbUQFCGOUQwEhUqsOX6hc1Gu14cPvftCfbGGQGlyEhMakOXylVl7lgFCm5\nNWWU0are1CznSByXISGo1IcvsiMFQVvKaMszXdIYiQBGw2UifKl+sy8lhxAIYVSKAEbDdHZ3ZiN8\nkT0p3zkihFGJcaEbgHzK1LGwqMqpoyP5dqB+ZqkJ54HlPGpdGBiWunUBieIkLMQu8+ErpWYjjjIy\ncuggU+sEYlfqJCwCGLEp17WWyg1NRjbgKCFDO1CZXD8QC86CRkNlcuOSkpN2UIdSQZvC77bcOsCx\n4ebDMWDUJZPBK0mdndHDU1g5oQz3VAZuFO/wUdcZjg03F7qgUZNK9tZTvRGh6zlfMtQVPSDz6xAq\nxjFgxKLSbrJUbzgyuLFGBTL4veZifUJZBDDqkpsNRQY30qhCRr/f3KxfiEQAo2rVnhSS+o1DRjfO\nqFJGv+fcrW8YRACjIrWciZmJDUFGN8qoUYa/79yug02MAEZJtV7+kImVfrSzY1kW8y3j332u18sm\nQwBDUjzXGmZiBS93WQrLYXPIwXJQ7zqbifU15wjgJhT3hf2ZWJEruR6UZbC55GSZaMr1OSfqCmAz\nu1jSt1S4c9Y6d18VMc63JV0i6R1Jn3X3zRHjZC6Au7u71d7eHroZRZK6a04SK2ldn3G1N2DI2PIX\nl7Qux4lJaDlJ6nNOYv1Pa0BncVkuFcBl74RlZmMk3Svp45LekPQLM3vE3bcPGecSSWe4+5lmdp6k\nNZLOj631AcX5ZWflVnNJr3iRn3Ej7mzUpOErZXOjFatq75ZVY9Wc1Oc8dB1t1HYlrRV3npblSm5F\n+WFJr7h7jySZ2UZJV0jaPmScKyStlyR3/zczm2xmp7v7myMnZp3ZCKFB3VJXZ1c800rLvx6xHvjQ\nfzHEd9QV02ccpYmDF0MMLAdx7dyVmk4jl+UIQ5duG/kUzbRscxTjtr87xm1yYJUE8AxJu4a83q1C\nKI82zp7+YUUBnKYFomJZbPOAcmGbZwQvosQdxCkyct0uCuTBNxrelMbOMydfXfIPY+hMfI71eyp0\nA+KVxmW3IfsEOdzA1qMr4cqsWaXqU05VY2KUk21yJQG8R9L7hrye2T9s5DjvLTNO5EFoAACaUSXP\nA/6FpLlm1mZmJ0haKunREeM8KukqSTKz8yUdiDr+CwAACspWwO5+1Mw+J+knOn4Z0jYzu67wtq91\n9380s0vN7D9VuAxpRWObDQBAtiV6Iw4AAFBQSRc0+pnZn5nZMTObGroteWNmXzezbWa22cweNrNJ\noduUF2Z2sZltN7OXzeyW0O3JGzObaWb/ZGYvmtkWM/t86DbllZmNMbPnzGzkYdBMIoArZGYzJX1C\nUk/otuTUTyTNc/ffkvSKpFsDtycXhtxI5yJJ8yQtM7OzwrYqd45I+oK7z5P0EUl/wmfcMH8q6aXQ\njYgLAVy5b0r636EbkVfu/n/c/Vj/y2dUOJMe9Ru8kY6790kauJEOYuLu/2/g1rvu/mtJ21S4DwJi\n1F8EXSrpu6HbEhcCuAJm9juSdrn7ltBtaRL/S9LjoRuRE1E30iEcGsTMZkn6LUn/FrYluTRQBOXm\nxKXkb8SRUmb2U0mnDx2kwhd9u6TbVOh+HvoeqjTKZ7zS3f+hf5yVkvrcfUOAJgI1M7NTJP2dpD/t\nr4QREzO7TNKb7r7ZzNqVk20wAdzP3T8RNdzMPiBplqQXzMxU6Br9dzP7sLu/lWATM6/UZzzAzD6r\nQhfT4kQa1BwquZEO6mRm41QI3wfc/ZHQ7cmhj0r6HTO7VNKJklrMbL27XxW4XXXhMqQqmdlOSee4\n+/7QbcmT/kdefkPSx9z9V6HbkxdmNlbSDhWeZrZX0rOSlrn7tqANyxkzWy/pv9z9C6HbkndmtkjS\nn7n774RuS704Blw9V066P1LmHkmnSPpp/2UGq0M3KA/c/aikgRvpvChpI+EbLzP7qKQ/kLTYzJ7v\nX34vDt0upB8VMAAAAVABAwAQAAEMAEAABDAAAAEQwAAABEAAAwAQAAEMAEAABDAAAAEQwAAABEAA\nAzlmZtcNuTvTL81sU+g2ASjgTlhAE+h/WMAmSavc/R9DtwcAFTDQLL4t6Z8IXyA9eBwhkHP9j3l8\nr7vfELotAI4jgIEcM7MPSfozSQtDtwXAcHRBA/n2J5KmSHqy/0SstaEbBKCAk7AAAAiAChgAgAAI\nYAAAAiCAAQAIgAAGACAAAhgAgAAIYAAAAiCAAQAI4P8DE/uBT4vSy6kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10c021510>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "qz = norm.pdf(z)\n",
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(z, np.ones(z.shape)/10., color='b', label='p(z)')\n",
    "plt.scatter(z, px, marker='+', color='r', label='p(x=0|z)')\n",
    "plt.scatter(z, qz, marker='+', color='g', label='q(z|x)')\n",
    "plt.xlabel('z')\n",
    "plt.ylim([0,1])\n",
    "plt.xlim([-5,5])\n",
    "_=plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean with 100 samples: 0.101230\n",
      "Mean with 100000 samples: 0.100134\n",
      "Standard deviation of 1000 means using 100 samples each: 0.007207\n"
     ]
    }
   ],
   "source": [
    "z2=np.random.randn(1000,100)\n",
    "px2=norm.pdf(z2, scale=0.5)*(1./10)/norm.pdf(z2)\n",
    "\n",
    "print \"Mean with %d samples: %f\"%(z2.shape[1], px2.mean(axis=1)[0])\n",
    "print \"Mean with %d samples: %f\"%(np.product(z2.shape), px2.mean())\n",
    "print \"Standard deviation of %d means using %d samples each: %f\"%(z2.shape[0], z2.shape[1], px2.mean(axis=1).std())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see the standard deviation of the means is ~ 3 times smaller, with the same number of samples. In general, the closer $q(z|x)$ is to $p(z|x)$ the fewer samples we'll need."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Ok, back to our general problem. We'll use importance sampling, by sampling from $q(z|x)$ instead of $p(z)$ to reduce the number of samples we need.\n",
    "\n",
    "$$\n",
    "p(x) = \\int_z p(x|z)q(z|x)\\frac{p(z)}{q(z|x)}\n",
    "$$\n",
    "\n",
    "$$\n",
    "p(x) = E_{z \\sim q(z|x)} p(x|z)\\frac{p(z)}{q(z|x)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Jensen's inequality\n",
    "\n",
    "$$\n",
    "\\log p(x) = \\log E_{z \\sim q(z|x)} p(x|z) \\frac{p(z)}{q(z|x)}\n",
    "$$\n",
    "\n",
    "\n",
    "$$\n",
    "\\log p(x) \\geq  E_{z \\sim q(z|x)} \\log \\left( p(x|z) \\frac{p(z)}{q(z|x)} \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using $\\log(ab) = \\log(a)+\\log(b)$\n",
    "\n",
    "$$\n",
    "\\log p(x) \\geq  E_{z \\sim q(z|x)} \\left[ \\log \\left(\\frac{p(z)}{q(z|x)} \\right) + \\log p(x|z) \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split the expectation\n",
    "\n",
    "$$\n",
    "\\log p(x) \\geq  E_{z \\sim q(z|x)} \\left[\\log \\frac{p(z)}{q(z|x)} \\right] + E_{z \\sim q(z|x)} \\log p(x|z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$log(a) = -log(1/a)$\n",
    "\n",
    "$$\n",
    "\\log p(x) \\geq  -E_{z \\sim q(z|x)} \\left[\\log \\frac{q(z|x)}{p(z)} \\right] + E_{z \\sim q(z|x)} \\log p(x|z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\log p(x) \\geq  -D_{kl}(q(z|x) || p(z)) +  E_{z \\sim q(z|x)} \\log p(x|z)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\log p(x) \\geq  E_{z \\sim q(z|x)} \\log p(x|z) - D_{kl}(q(z|x) || p(z)) \n",
    "$$\n",
    "\n",
    "So, if we maximize the right hand size, we're maximizing a lower bound of the log probability of the data, which is exactly what we want to do. Provided $q(z|x)$ and $p(z)$ are simple distributions we can even calculate the KL term analytically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, what's missing for this to be an equality and not a lower bound? We'll add a term C to the right hand side and see if we can't derive what C should be in order for the equality to hold.\n",
    "\n",
    "$$\n",
    "\\log p(x) =  E_{z \\sim q(z|x)} \\log \\left[ p(x|z) \\frac{p(z)}{q(z|x)} \\right] + C\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = \\log p(x) - E_{z \\sim q(z|x)} \\log \\left[ p(x|z) \\frac{p(z)}{q(z|x)} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = - \\int_z \\left[ q(z|x)\\log \\left[ p(x|z) \\frac{p(z)}{q(z|x)} \\right]\\right] + \\log p(x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = - \\int_z \\left[ q(z|x)\\log \\left[ p(x|z) \\frac{p(z)}{q(z|x)} \\right] - \\log p(x) \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = - \\int_z q(z|x)\\log \\left[ \\frac{p(x|z)p(z)}{p(x)q(z|x)} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "p(z|x) = \\frac{p(x|z)p(z)}{p(x)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = - \\int_z q(z|x)\\log \\frac{p(z|x)}{q(z|x)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = \\int_z q(z|x)\\log \\frac{q(z|x)}{p(z|x)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C = D_{kl}(q(z|x) || p(z|x))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's insert $C$ back in our optimization objective"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\log p(x) = E_{z \\sim q(z|x)} \\log p(x|z) - D_{kl}(q(z|x) || p(z)) + D_{kl}(q(z|x) || p(z|x))\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\log p(x) - D_{kl}(q(z|x) || p(z|x)) =  E_{z \\sim q(z|x)} \\log p(x|z) - D_{kl}(q(z|x) || p(z))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So if we're maximizing the right hand side we're maximizing the log probability of the data or minimizing the KL divergence between $q(z|x)$ and $p(z|x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}