Reparameterization Trick.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory behind Reparameterization Trick - somemore additions to the original version by Goker Erdogan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will understand the reparameterization trick used by Kingma and Welling (2014) to train their variational autoencoder.  \n",
    "\n",
    "Assume we have a normal distribution $q$ that is parameterized by $\\theta$, specifically $q_{\\theta}(x) = N(\\theta,1)$. We want to solve the below problem\n",
    "$$\n",
    "\\text{min}_{\\theta} \\quad E_q[x^2]\n",
    "$$\n",
    "This is of course a rather silly problem and the optimal $\\theta$ is obvious. We want to understand how the reparameterization trick helps in calculating the gradient of this objective $E_q[x^2]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the expected value of $E_q[x^2]$ ?\n",
    "<br> It is nothing but variance + (mean)^2. Therefore, in this case, it is $\\theta^2$ + 1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to calculate $\\nabla_{\\theta} E_q[x^2]$ is as follows\n",
    "$$\n",
    "\\nabla_{\\theta} E_q[x^2] = \\nabla_{\\theta} \\int q_{\\theta}(x) x^2 dx = \\int x^2 \\nabla_{\\theta} q_{\\theta}(x) \\frac{q_{\\theta}(x)}{q_{\\theta}(x)} dx = \\int q_{\\theta}(x) \\nabla_{\\theta} \\log q_{\\theta}(x) x^2 dx = E_q[x^2 \\nabla_{\\theta} \\log q_{\\theta}(x)]\n",
    "$$\n",
    "\n",
    "For our example where $q_{\\theta}(x) = N(\\theta,1)$, this method gives\n",
    "$$\n",
    "\\nabla_{\\theta} E_q[x^2] = E_q[x^2 (x-\\theta)]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " What is the variance of $E_q[x^2 (x-\\theta)]$ in this case? <br>\n",
    "It can be found using some complicated, yet straightforward math. $Var(\\mu^2 \\epsilon + 2\\mu \\epsilon^2 + \\epsilon^3) \\le \\mu^4Var(\\epsilon) \\le \\mu^4$ <br>\n",
    "Therefore, variance of the sample mean of $E_q[x^2 (x-\\theta)]$ will be at least $\\frac{\\mu^4}{\\sqrt{n}}.  $ We can indeed verify that this variance increases with $\\mu$. Note that the pathwise gradient does not have this problem.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reparameterization trick is a way to rewrite the expectation so that the distribution with respect to which we take the expectation is independent of parameter $\\theta$. To achieve this, we need to make the stochastic element in $q$ independent of $\\theta$. Hence, we write $x$ as\n",
    "$$\n",
    "x = \\theta + \\epsilon, \\quad \\epsilon \\sim N(0,1)\n",
    "$$\n",
    "Then, we can write \n",
    "$$\n",
    "E_q[x^2] = E_p[(\\theta+\\epsilon)^2]\n",
    "$$ \n",
    "where $p$ is the distribution of $\\epsilon$, i.e., $N(0,1)$. Now we can write the derivative of $E_q[x^2]$ as follows\n",
    "$$\n",
    "\\nabla_{\\theta} E_q[x^2] = \\nabla_{\\theta} E_p[(\\theta+\\epsilon)^2] = E_p[2(\\theta+\\epsilon)]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above case, it is easy to note that the sample mean will have a variance of $\\frac{4}{\\sqrt{n}}$,which looks independent of $\\mu$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us compare the variances of the two methods; we are hoping to see that the first method has high variance while reparameterization trick decreases the variance substantially."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With mean = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.95569048662\n",
      "3.9150255039\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "N = 1000\n",
    "theta = 2.0\n",
    "eps = np.random.randn(N)\n",
    "x = theta + eps\n",
    "\n",
    "grad1 = lambda x: np.sum(np.square(x)*(x-theta)) / x.size\n",
    "grad2 = lambda eps: np.sum(2*(theta + eps)) / x.size\n",
    "\n",
    "print(grad1(x))\n",
    "print(grad2(eps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot the variance for different sample sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.81095352  4.03655932  3.96279395  4.0045522   4.00203414]\n",
      "[ 4.04018422  3.9878957   4.00037079  3.99667789  4.00066778]\n",
      "\n",
      "[  1.32731049e+01   9.12063666e-01   8.03604324e-02   1.01701641e-02\n",
      "   7.13464452e-04]\n",
      "[  3.96371303e-01   3.66080170e-02   5.87655498e-03   3.65440156e-04\n",
      "   3.89156655e-05]\n"
     ]
    }
   ],
   "source": [
    "Ns = [10, 100, 1000, 10000, 100000]\n",
    "reps = 100\n",
    "\n",
    "means1 = np.zeros(len(Ns))\n",
    "vars1 = np.zeros(len(Ns))\n",
    "means2 = np.zeros(len(Ns))\n",
    "vars2 = np.zeros(len(Ns))\n",
    "\n",
    "est1 = np.zeros(reps)\n",
    "est2 = np.zeros(reps)\n",
    "for i, N in enumerate(Ns):\n",
    "    for r in range(reps):\n",
    "        x = np.random.randn(N) + theta\n",
    "        est1[r] = grad1(x)\n",
    "        eps = np.random.randn(N)\n",
    "        est2[r] = grad2(eps)\n",
    "    means1[i] = np.mean(est1)\n",
    "    means2[i] = np.mean(est2)\n",
    "    vars1[i] = np.var(est1)\n",
    "    vars2[i] = np.var(est2)\n",
    "    \n",
    "print(means1)\n",
    "print(means2)\n",
    "print()\n",
    "print(vars1)\n",
    "print(vars2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1dd308df828>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHxxJREFUeJzt3XlwnPWdJvDn263WLUu2DlvW0RLYYBsjbCykNiR4ggcw\nCQ4QwGCQgN3ZeLO1R2Z3kyzZrSKhikolNakJs9mtsM4xY2zjEA4nhAHCFY4ZfMm2AGOBZce2DsuS\nLFuyZN3q7/7RbVmSdbS63+5f99vPp6pLfbzHo7fkR6/e9/X7E1UFERHFPofpAEREZA0WOhGRTbDQ\niYhsgoVORGQTLHQiIptgoRMR2QQLnYjIJljoREQ2wUInIrKJhEiuLCcnR0tKSiK5SiKimLd///4z\nqpo703QRLfSSkhLU1NREcpVERDFPRE4GMh0PuRAR2QQLnYjIJljoREQ2EdFj6EREszU0NISmpib0\n9/ebjhJ2ycnJKCwshMvlCmp+FjoRRbWmpiZkZGSgpKQEImI6TtioKjo6OtDU1ITS0tKglsFDLkQU\n1fr7+5GdnW3rMgcAEUF2dnZIf4mw0Iko6tm9zC8K9fuMiUL/l/oz+MV7x0zHICKKajFR6B/Wt+On\nb36B1vP2PylCRPb23nvv4aOPPgrLsmOi0B+qLMaIV7Fjb4PpKEREQRseHmahu7PTsOaqXOzY24Ch\nEa/pOEQUR06cOIGlS5fim9/8Jq655hrcdttt6OvrAwDU1tbC4/GgrKwM99xzD86dO3fZ/I899hi+\n9a1vobKyEhs2bMAzzzyDn/3sZ1ixYgU+/PBDS7PGzGWL1R43/t2zNXj7cCvuuDbfdBwiMuDJP36G\nw6fOW7rMZQvn4Afrr5l2mvr6euzYsQO//OUvsWHDBrz00kuoqqrCI488gp///OdYs2YNnnjiCTz5\n5JN4+umnL5u/qakJH330EZxOJ374wx8iPT0d3/nOdyz9PoAY2UMHgK8syUNBVgq27QnoHjVERJYp\nLS3FihUrAACrVq3CiRMn0NXVhc7OTqxZswYA8Oijj+KDDz6YdP77778fTqcz7DljZg/d6RA8VFmM\nv/vTFzja1oNFeemmIxFRhM20Jx0uSUlJo8+dTufoIZdApaWlWR1pUjGzhw4AD9xQBJdTsJ176URk\nWGZmJubOnTt6HHzr1q2je+vTycjIQHd3d1gyxVSh56Qn4Y7l+XhxfxN6B4dNxyGiOLdlyxZ897vf\nRVlZGWpra/HEE0/MOM/69euxc+fOsJwUFVW1dIHTKS8v11AHuNh34izuf2YXfvyNa/FgRbFFyYgo\nWtXV1WHp0qWmY0TMZN+viOxX1fKZ5o2pPXQAKHfPxZIFGdi6+yQi+cuIiCjaxVyhiwiqPG58duo8\nahs7TcchIooaMVfoAHD3ygKkJyVg626eHCUiuigmCz09KQH3rCzAq5+04OyFQdNxiIiiQkwWOgBU\nedwYHPbihZpG01GIiKJCzBb61QsyUFE6D9v3NMDr5clRIqKYLXTAd3+XhrO9+KC+3XQUIoojP/rR\nj0xHmNSMhS4ivxGRNhE5NOa9vxORz0XkExHZKSJZ4Y05uduvWYCc9CRs48lRIooQVcVTTz1lOsak\nAtlD/ycA6ya89xaA5apaBuAIgO9bnCsgiQkObKwowjuft6HxbK+JCEQUB06cOIGrr74ajzzyCBYt\nWoS+vj6sWLECDz/8sOlo48x4cy5V/UBESia89+aYl7sB3GdtrMBtrCjG//3zUezY24DvrVtiKgYR\nRcLrjwOnP7V2mQuuBe748YyT1dfXY8uWLXj22WeRnp6O2tpaa3NYwIpj6P8WwOtTfSgim0SkRkRq\n2tutP9a9MCsFa5fOx/P7GjEwPGL58omIAMDtdsPj8ZiOMa2Qbp8rIv8LwDCA7VNNo6qbAWwGfPdy\nCWV9U6n2uPHW4Va8ceg07lpREI5VEFE0CGBPOlwidQvcUAS9hy4ijwG4E8DDavimKl9alIOS7FSe\nHCWiiHC5XBgaGjId4zJBFbqIrAPwPQBfV1XjZyMdDt/9XfadOIe6FmuHpyIimmjTpk0oKyuLupOi\ngVy2uAPALgBXi0iTiPwNgP8DIAPAWyJSKyLPhDnnjO5bVYikBAf30onIciUlJTh0aPTKbfzkJz9B\nXV0dtm+f8mizEYFc5bJxkrd/HYYsIclKTcT66xZi58FmPH7HEmQku0xHIiKKqJj+n6ITVXvc6B0c\nwc6DzaajEBFFnK0K/bqiLJQVZmIbB78gspV4+fcc6vdpq0IHfHdhPNLag73Hz5qOQkQWSE5ORkdH\nh+1LXVXR0dGB5OTkoJcR0nXo0Wh92UI89ephbN19EpVXZJuOQ0QhKiwsRFNTE8LxHxOjTXJyMgoL\nC4Oe33aFnpLoxP3lRdjy0Qm0dfcjLyP433ZEZJ7L5UJpaanpGDHBdodcAODhymIMexXP7+XgF0QU\nP2xZ6FfkpuPLi3Pw3N4GDI94TcchIooIWxY64Ds52tLVj3c/bzMdhYgoImxb6GuX5CE/Mxlb+T9H\niShO2LbQE5wObKwoxof1Z3D8zAXTcYiIws62hQ4AD95QhASHYDv30okoDti60PPmJOP25Qvwwv4m\n9A9x8AsisjdbFzrgu79LV98Q/vjxKdNRiIjCyvaFXlk6D4vz0nlbXSKyPdsXuohv8IuPm7rwcWOn\n6ThERGFj+0IHgHuuL0BqopN76URka3FR6HOSXbh7ZQFe+fgUOnsHTcchIgqLuCh0AKiqdGNg2IsX\n9zeZjkJEFBZxU+jLFs5BuXsutu9pgNdr7/sqE1F8iptCB4Dq1W4cP3MB/3rsjOkoRESWm7HQReQ3\nItImIofGvDdPRN4SkXr/17nhjWmNdcsXIDstEVt38eQoEdlPIHvo/wRg3YT3HgfwjqouBvCO/3XU\nS0pwYsMNRXi7rhWnOvtMxyEistSMha6qHwCYOEDnXQC2+J9vAXC3xbnC5qGKYiiAHXsbTEchIrJU\nsMfQ56tqi//5aQDzLcoTdkXzUnHL1Xn47b5GDA5z8Asiso+QT4qqbyjuKS8bEZFNIlIjIjXRMshr\n1Wo32rsH8Obh06ajEBFZJthCbxWRfADwf51yWCBV3ayq5apanpubG+TqrLVmcS6K5qXw5CgR2Uqw\nhf4KgEf9zx8F8Adr4kSGwyF4uNKNPcfP4khrt+k4RESWCOSyxR0AdgG4WkSaRORvAPwYwK0iUg/g\nr/2vY8qG8iIkJjh4fxciso2EmSZQ1Y1TfLTW4iwRNS8tEXdem4+XDzTjf6xbgrSkGTcFEVFUi6v/\nKTpR1Wo3egaG8fvaZtNRiIhCFteFvrIoC8vy52DrrpPwXaxDRBS74rrQRQTVq934/HQ39p88ZzoO\nEVFI4rrQAeCuFQuRkZSArTw5SkQxLu4LPTUxAfeuKsRrn7bgTM+A6ThEREGL+0IHgCqPG0Mjit/V\nNJqOQkQUNBY6gEV56bjxymxs392AEQ5+QUQxioXuV+Vxo7mzD+99MeVdDIiIohoL3e/WZfORl5HE\nk6NEFLNY6H4upwMbK4rx/pF2NHT0mo5DRDRrLPQxNlYUwyGC7Xu5l05EsYeFPsaCzGTctmw+frev\nEf1DI6bjEBHNCgt9giqPG+d6h/Dapy0zT0xEFEVY6BPceGU2rshN48lRIoo5LPQJRARVlW4cbOjE\noeYu03GIiALGQp/EvasKkezi4BdEFFtY6JPITHHh7hUF+EPtKXT1DZmOQ0QUEBb6FKo8bvQNjeDl\nA02moxARBYSFPoXlBZlYUZSFrbs5+AURxQYW+jSqPW78pf0Cdh3rMB2FiGhGLPRpfK0sH1mpLl7C\nSEQxIaRCF5H/KiKficghEdkhIslWBYsGyS4nHigvwpuHW9F6vt90HCKiaQVd6CJSAOC/AChX1eUA\nnAAetCpYtHioshheVezY22A6ChHRtEI95JIAIEVEEgCkAjgVeqTo4s5Ow82Lc7FjbwOGRrym4xAR\nTSnoQlfVZgA/BdAAoAVAl6q+aVWwaFLtcaP1/ADePtxqOgoR0ZRCOeQyF8BdAEoBLASQJiJVk0y3\nSURqRKSmvb09+KQGfWVJHgqyUnhylIiiWiiHXP4awHFVbVfVIQAvA7hx4kSqullVy1W1PDc3N4TV\nmeN0CB6qLMZHxzpwtK3HdBwiokmFUugNADwikioiAmAtgDprYkWfB24ogssp2L6He+lEFJ1COYa+\nB8CLAA4A+NS/rM0W5Yo6OelJuGN5Pl7c34TewWHTcYiILhPSVS6q+gNVXaKqy1W1WlUHrAoWjapX\nu9HdP4xXam13MQ8R2QD/p+gslLvnYsmCDDy7i/d3IaLow0KfBRFBlceNwy3ncbCx03QcIqJxWOiz\ndPfKAqQnJXDwCyKKOiz0WUpPSsA3ri/Aq5+04OyFQdNxiIhGsdCDUOVxY3DYixdqGk1HISIaxUIP\nwlXzM1BROg/b9pyE18uTo0QUHVjoQar2uNF4tg/v18fm7QyIyH5Y6EG6/ZoFyElPwrZdPDlKRNGB\nhR6kxAQHNlYU4d0v2tB4ttd0HCIiFnooNlYUQwAOfkFEUYGFHoKFWSlYu3Q+nt/XiIHhEdNxiCjO\nsdBDVO1xo+PCIN44dNp0FCKKcyz0EH1pUQ5KslOxlSdHicgwFnqIHA7f/V1qTp5DXct503GIKI6x\n0C1w36pCJCU4eH8XIjKKhW6BrNRErL9uIXYebEZ3/5DpOEQUp1joFqn2uNE7OIKdB5tNRyGiOMVC\nt8h1RVkoK8zEVg5+QUSGsNAtVOVxo76tB3uOnzUdhYjiEAvdQuvLFiIzxcWTo0RkBAvdQimJTty3\nqhBvHDqNtu5+03GIKM6EVOgikiUiL4rI5yJSJyKrrQoWqx6uLMawV/H8Xg5+QUSRFeoe+j8AeENV\nlwC4DkBd6JFi2xW56fjy4hw8t7cBwyNe03GIKI4EXegikgngZgC/BgBVHVTVTquCxbIqjxstXf14\n5/M201GIKI6EsodeCqAdwD+KyEER+ZWIpFmUK6atXZKH/MxknhwloogKpdATAFwP4BequhLABQCP\nT5xIRDaJSI2I1LS3x8dwbQlOBzZWFOPD+jM4fuaC6ThEFCdCKfQmAE2qusf/+kX4Cn4cVd2squWq\nWp6bmxvC6mLLgzcUIcEh2M69dCKKkKALXVVPA2gUkav9b60FcNiSVDaQNycZty9fgBf2N6FvkINf\nEFH4hXqVy38GsF1EPgGwAsCPQo9kH9UeN7r6hvDHT06ZjkJEcSCkQlfVWv/hlDJVvVtVz1kVzA4q\nS+dhcV46D7sQUUTwf4qGkYigerUbHzd14eNGXtFJROHFQg+ze1YWIDXRyUsYiSjsWOhhlpHswt0r\nC/DKx6fQ2TtoOg4R2RgLPQKqKt0YGPbixf1NpqMQkY2x0CNg2cI5KHfPxbbdJ+H1cvALIgoPFnqE\nVK9240RHL/712BnTUYjIpljoEbJu+QJkpyVi6y6eHCWi8GChR0hSghMbbijC23WtONXZZzoOEdkQ\nCz2CHqoohgLYsbfBdBQisiEWegQVzUvFLVfnYcfeRgwOc/ALIrIWCz3Cqla7caZnAG8ePm06ChHZ\nDAs9wtYszkXRvBSeHCUiy7HQI8zhEDxc6cae42dxpLXbdBwishEWugEbyouQmODg/V2IyFIsdAPm\npSXizmvz8fKBZvQMDJuOQ0Q2wUI3pGq1Gz0Dw/j9wWbTUYjIJljohqwsysI1C+dg2+6TUOX9XYgo\ndCx0Q0QEVR43Pj/djf0nOdATEYWOhW7QXSsWIiMpAVt5cpSILMBCNyg1MQH3rirEa5+24EzPgOk4\nRBTjWOiGVXncGBpRPL+v0XQUIopxIRe6iDhF5KCIvGpFoHizKC8dN16Zjef2NGCEg18QUQis2EP/\nNoA6C5YTt6o8bjR39uG9L9pMRyGiGBZSoYtIIYCvAfiVNXHi063L5iMvI4knR4koJKHuoT8N4HsA\nprwXrIhsEpEaEalpb28PcXX25HI6sLGiGO8facfJjgum4xBRjAq60EXkTgBtqrp/uulUdbOqlqtq\neW5ubrCrs72NFcVwiOC5PRz8goiCE8oe+k0Avi4iJwD8FsAtIrLNklRxaEFmMm5bNh+/q2lE/9CI\n6ThEFIOCLnRV/b6qFqpqCYAHAbyrqlWWJYtD1R43zvUO4bVPW0xHIaIYxOvQo8jqK7NxRW4aT44S\nUVAsKXRVfU9V77RiWfFMRFBV6cbBhk4cau4yHYeIYgz30KPMvasKkezi4BdENHss9CiTmeLC3SsK\n8PvaZnT1DZmOQ0QxhIUehao8bvQPefHygSbTUYgohrDQo9DygkysKMrCVg5+QUSzwEKPUtUeN/7S\nfgG7jnWYjkJEMYKFHqW+VpaPrFQXL2EkooCx0KNUssuJB8qL8ObhVpzu6jcdh4hiAAs9ij1UWQyv\nKnbs5f1diGhmLPQo5s5Ow5qrcvHbfQ0YGpnyhpZERABY6FGvqtKN1vMDePtwq+koRBTlWOhR7itL\n8lCQlcKTo0Q0IxZ6lHM6BA9VFuOjYx042tZtOg4RRTEWegx44IYiuJyCbbt5cpSIpsZCjwE56Un4\n6rX5eOlAE3oHh03HIaIoxUKPEVUeN7r7h/FK7SnTUYgoSrHQY0S5ey6WLMjAs7t4fxcimhwLPUaI\nCKo8bhxuOY+DjZ2m4xBRFGKhx5C7VxYgPSkB23bxEkYiuhwLPYakJyXgG9cX4NVPW3D2wqDpOEQU\nZVjoMabK48bgsBcv1DSajkJEUSboQheRIhH5s4gcFpHPROTbVgajyV01PwMVpfOwbc9JeL08OUpE\nl4Syhz4M4L+r6jIAHgD/UUSWWROLplPtcaPxbB/er283HYWIokjQha6qLap6wP+8G0AdgAKrgtHU\nbr9mAXLSk3hylIjGseQYuoiUAFgJYI8Vy6PpJSY4sLGiCO9+0YbGs72m4xBRlAi50EUkHcBLAP5W\nVc9P8vkmEakRkZr2dh4isMrGimIIwMEviGhUSIUuIi74yny7qr482TSqullVy1W1PDc3N5TV0RgL\ns1Kwdul8PL+vEQPDI6bjEFEUCOUqFwHwawB1qvr31kWiQFV73Oi4MIg3Dp02HYWIokAoe+g3AagG\ncIuI1PofX7UoFwXgS4tyUJKdiq08OUpEABKCnVFV/wWAWJiFZsnh8N3f5al/rkNdy3kszZ9jOhIR\nGcT/KRrj7ltViKQEB7ZxiDqiuMdCj3FZqYlYf91C7DzYjO7+IdNxiMggFroNVHvc6B0cwc6Dzaaj\nEJFBLHQbuK4oC2WFmdjKwS+I4hoL3SaqPG7Ut/Vgz/GzpqMQkSEsdJtYX7YQmSkubOXJUaK4xUK3\niZREJ+5fVYg/HTqNtu5+03GIyAAWuo087HFj2Kt4fi8HvyCKRyx0GynNScOXF+fgub0NGB7xmo5D\nRBHGQreZKo8bLV39eOfzNtNRiCjCWOg2s3ZJHvIzk/H/3j+G/SfPoquX/9mIKF4EfS8Xik4JTgf+\n/c1X4Id/PIx7f7ELAJCbkYTFeelYnJeORXnpWJSXgcXz05GdlgjfTTOJyA5Y6Db02E2lWLt0Purb\nulHf2oP6Nt/jpQPN6BkYHp0uK9XlL/kMLPIX/uL56VgwJ5lFTxSDWOg2VTQvFUXzUnHLkvmj76kq\nTp/vR31rD476S/5oWzdeP9SCzjGHZtKTEnDlmD163959BgrnpsDhYNETRSsWehwREeRnpiA/MwU3\nX3Vp9ChVRceFQV/Rt/fgaGs36tt68MGRdry4v2l0uqQEB67M9e3Fjz18485OhcvJ0zFEprHQCSKC\nnPQk5KQnYfWV2eM+6+odwtH2bt8evf/wTc2Jc/hD7anRaVxOQUl2GhbPH3/4pjQnDckuZ6S/HaK4\nxUKnaWWmurDKPQ+r3PPGvX9hYBjH2i8duqlv7UFdSzfeOHQaXv/9wRwCFM9LHT0Ju8i/d39lbjrS\nkvijR2Q1/quioKQlJaCsMAtlhVnj3u8fGsHxMxfGHaM/2taD94+0YWjk0p0gC7JSRvfkF/lPxi7K\nzUBmqivS3wqRbcRGoe9+Bjj6FpAyd+ZHchbgjI1vy46SXU4szZ9z2XB4QyNenOzoxVF/yfvKvge7\n/9KBgeFL/6s1LyPpUtHPzxjdq+cllkQzi43m8w4BvWeBjmNA3zmgvwvANPf9TpoDpGQF9gtg7CMh\nKWLfUrxxOR3+k6jpABaMvj/iVTSf60O9f0/+YtFPvMRybqrr0jX0Y/bqeYkl0SUSyQERysvLtaam\nJvQFeUd8pd53Dujr9H8N8KEjUy/XlTr7XwIpcwFXCsBSsdTFSyzHnow91taDI23dU15ieekyS15i\nSfYiIvtVtXzG6UIpdBFZB+AfADgB/EpVfzzd9JYVerBUgYHuGUp/7C+Is5eejwxOvVxn0hRlP8Nf\nCUkZ/EUwSxcvsRzdm2+9dPimrXtgdLpkl+8Sy0vH6X1X3/ASS4pFYS90EXECOALgVgBNAPYB2Kiq\nh6eax3ihB0sVGOoN8K+ACX8xDPVOvVxxzvKvAf8viORMwMHLASea7BLLo209aO7sG53G5RSU5qRd\ndvjGnZ2KRKcDTofwEA5FnUALPZRj6BUAjqrqX/wr/C2AuwBMWegxSwRITPM9MgtnN+9QP9Af4GGh\nntNAe53vl8LA+ekC+Up9xr1/h/8hvq8Q//OxryebxnFpunGvZ5oHkyxjqnkwzXonm2e6LL6vmUmC\nVUVZWFU8d9xfPoFcYjmWQ4AEhwMOh/+r+O6R4xBBgkPgnPiQSd7zv5/glNH5HI7xX6ebz+n0z+8Q\nOB0OOB0Y/1UAp9MxOs3EZU+ZdZLMl83nnCEbf+lFrVAKvQDA2JEUmgBUhhbHhlzJgGsBkLFg5mnH\nGhkac54ggMe545f+QpjuhHE88Zd+mjhQJoKyCb8EdI7Aq4IRFXh1/FZTiP8rxn3FiO+hkMs/G12G\nQMcsTSG4+IfwxHnGvj+6TvW9M3YdYzNNRnX8Z9PNN/b12Om8AEYgGJpi2svWOaHUrah4mWwpAS74\n4mTTZZ5x0bP4Jmb7/fbe9lMsW33HLOeanbBf5SIimwBsAoDi4uJwr84+nC4gLcf3mA2vFxjo8p0r\nUK+vHdR/WeC41/6vk76HCa8nfj5mXqj/9RTTXDaPzjzNuNcBLndijgCWK1A41Qvn2G3k21Djt+ll\nhyU1sM8u+3w2n419qlBVeFX90X2vFRe/FS9U4f+l5P9s3HQ6+jkuTotLzxWA179NvP6Fei/mm2yZ\n8H2dLO9kuxI66bvTz3PZZrhsIp1i3qkXOuV6ZpNjiokDmTczLXOaqawRSqE3Ayga87rQ/944qroZ\nwGbAdww9hPVRIByOS4ddyBb8B5s4eAHNKJSfkX0AFotIqYgkAngQwCvWxCIiotkKeg9dVYdF5D8B\n+BN8ly3+RlU/sywZERHNSkjH0FX1NQCvWZSFiIhCwMNyREQ2wUInIrIJFjoRkU2w0ImIbIKFTkRk\nExG9fa6ItAM4GeTsOQDOWBjHKsw1O8w1O8w1O9GaCwgtm1tVc2eaKKKFHgoRqQnkbmORxlyzw1yz\nw1yzE625gMhk4yEXIiKbYKETEdlELBX6ZtMBpsBcs8Ncs8NcsxOtuYAIZIuZY+hERDS9WNpDJyKi\naURdoYvIOhH5QkSOisjjk3wuIvK//Z9/IiLXR0muvxKRLhGp9T+eiECm34hIm4gcmuJzU9tqplwR\n31b+9RaJyJ9F5LCIfCYi355kmohvswBzmfj5ShaRvSLysT/Xk5NMY2J7BZLLyM+Yf91OETkoIq9O\n8ll4t9foqCZR8IDvNrzHAFwBIBHAxwCWTZjmqwBeh++e/x4Ae6Ik118BeDXC2+tmANcDODTF5xHf\nVgHmivi28q83H8D1/ucZ8A1yHg0/X4HkMvHzJQDS/c9dAPYA8ETB9gokl5GfMf+6/xuA5yZbf7i3\nV7TtoY8OPK2qgwAuDjw91l0AnlWf3QCyRCQ/CnJFnKp+AODsNJOY2FaB5DJCVVtU9YD/eTeAOvjG\nxh0r4tsswFwR598GPf6XLv9j4kk3E9srkFxGiEghgK8B+NUUk4R1e0VboU828PTEH+xApjGRCwBu\n9P8Z9bqIXBPmTIEwsa0CZXRbiUgJgJXw7d2NZXSbTZMLMLDN/IcPagG0AXhLVaNiewWQCzDzM/Y0\ngO/BP2zrJMK6vaKt0GPZAQDFqloG4OcAfm84TzQzuq1EJB3ASwD+VlXPR3Ld05khl5FtpqojqroC\nvjGDK0RkeSTWO5MAckV8e4nInQDaVHV/uNc1lWgr9EAGng5ocOpI51LV8xf/DFTfSE4uEckJc66Z\nmNhWMzK5rUTEBV9pblfVlyeZxMg2mymX6Z8vVe0E8GcA6yZ8ZPRnbKpchrbXTQC+LiIn4Dsse4uI\nbJswTVi3V7QVeiADT78C4BH/2WIPgC5VbTGdS0QWiIj4n1fAt207wpxrJia21YxMbSv/On8NoE5V\n/36KySK+zQLJZWKbiUiuiGT5n6cAuBXA5xMmM7G9ZsxlYnup6vdVtVBVS+DriHdVtWrCZGHdXiGN\nKWo1nWLgaRH5lv/zZ+Abw/SrAI4C6AXwb6Ik130A/oOIDAPoA/Cg+k9rh4uI7IDvbH6OiDQB+AF8\nJ4iMbasAc0V8W/ndBKAawKf+468A8D8BFI/JZmKbBZLLxDbLB7BFRJzwFeLvVPVV0/8eA8xl6mfs\nMpHcXvyfokRENhFth1yIiChILHQiIptgoRMR2QQLnYjIJljoREQ2wUInIrIJFjoRkU2w0ImIbOL/\nA6hUTIMC4tooAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1dd3281cac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(vars1)\n",
    "plt.plot(vars2)\n",
    "plt.legend(['no rt', 'rt'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With Mean = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.87862307368\n",
      "7.9292520907\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "N = 1000\n",
    "theta = 4.0\n",
    "eps = np.random.randn(N)\n",
    "x = theta + eps\n",
    "\n",
    "grad1 = lambda x: np.sum(np.square(x)*(x-theta)) / x.size\n",
    "grad2 = lambda eps: np.sum(2*(theta + eps)) / x.size\n",
    "\n",
    "print(grad1(x))\n",
    "print(grad2(eps))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 8.26682027  7.82919365  7.99299786  7.97066258  8.00084004]\n",
      "[ 8.04052111  7.99936387  7.9954101   7.99919523  8.00066903]\n",
      "\n",
      "[  6.95956617e+01   6.44329020e+00   5.44355317e-01   3.08962462e-02\n",
      "   7.23940750e-03]\n",
      "[  4.18539528e-01   4.05795454e-02   4.69016086e-03   4.58622290e-04\n",
      "   4.22305479e-05]\n"
     ]
    }
   ],
   "source": [
    "Ns = [10, 100, 1000, 10000, 100000]\n",
    "reps = 100\n",
    "\n",
    "means1 = np.zeros(len(Ns))\n",
    "vars1 = np.zeros(len(Ns))\n",
    "means2 = np.zeros(len(Ns))\n",
    "vars2 = np.zeros(len(Ns))\n",
    "\n",
    "est1 = np.zeros(reps)\n",
    "est2 = np.zeros(reps)\n",
    "for i, N in enumerate(Ns):\n",
    "    for r in range(reps):\n",
    "        x = np.random.randn(N) + theta\n",
    "        est1[r] = grad1(x)\n",
    "        eps = np.random.randn(N)\n",
    "        est2[r] = grad2(eps)\n",
    "    means1[i] = np.mean(est1)\n",
    "    means2[i] = np.mean(est2)\n",
    "    vars1[i] = np.var(est1)\n",
    "    vars2[i] = np.var(est2)\n",
    "    \n",
    "print(means1)\n",
    "print(means2)\n",
    "print()\n",
    "print(vars1)\n",
    "print(vars2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1dd338c4cc0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtwVPeZ5vHv2y0JCSSQQAJkdGk8EGNwsMAKIJybDfY4\n4yQml3HisY2ytTVstmanZnbnsp6tLc+4ajblVO3OJJuq3ZQnmbV8iZ2MY8ceO/HEw9gTZ7hZBoGN\nIQE7QhIRIC7iYhC6vftHH4iQJdSSuvt0t55Plaq7T5/ufjglnv7p9K/PMXdHRESyXyTsACIikhwq\ndBGRHKFCFxHJESp0EZEcoUIXEckRKnQRkRyhQhcRyREqdBGRHKFCFxHJEXnpfLHy8nKPxWLpfEkR\nkaz35ptvHnf3irHWS2uhx2Ixmpub0/mSIiJZz8wOJbKedrmIiOQIFbqISI5QoYuI5Igx96Gb2XXA\n94csuhZ4EHgsWB4DWoG73f1U8iOKyFTW19dHR0cHPT09YUdJucLCQqqqqsjPz5/Q48csdHf/BVAH\nYGZR4DDwHPAAsNndHzazB4Lb/3VCKURERtHR0UFJSQmxWAwzCztOyrg7J06coKOjg4ULF07oOca7\ny2Ud8K67HwLuApqC5U3AhgklEBG5ip6eHubMmZPTZQ5gZsyZM2dSf4mMt9C/DDwVXJ/n7p3B9SPA\nvAmnEBG5ilwv80sm++9MuNDNrAD4LPAPw+/z+HnsRjyXnZltMrNmM2vu6uqaUMh/O3ic//PawQk9\nVkRkqhjPCP1TwE53PxrcPmpmlQDB5bGRHuTuj7h7vbvXV1SM+UWnEf3rL7v4Xz/9JZ2nL0zo8SIi\nmeK1115jy5YtKXnu8RT6PfxmdwvAC0BjcL0ReD5ZoYa7f00tg+48ua0tVS8hIpJy/f394Re6mc0A\nbgOeHbL4YeA2MzsArA9up0T17OmsWzKXp3a0cbF/IFUvIyLyAa2trVx//fX8/u//PsuWLeP222/n\nwoX43oKWlhbWrFnD8uXL+dznPsepUx+cuf2Vr3yFr371q6xevZq7776bb3/72/zt3/4tdXV1vP76\n60nNmtCxXNz9fWDOsGUniM96SYvGtTH+ed8OXtrTyedXVqXrZUUkgzz0j3t559dnkvqcS6+ZyV9+\nZtlV1zlw4ABPPfUUf/d3f8fdd9/ND3/4Q+677z42btzIt771LT7xiU/w4IMP8tBDD/GNb3zjA4/v\n6Ohgy5YtRKNR/uqv/ori4mL+9E//NKn/Dsiib4re/FvlXFsxg6YtrWFHEZEpZuHChdTV1QFw0003\n0drayunTp+nu7uYTn/gEAI2NjfzsZz8b8fG/+7u/SzQaTXnOtB5tcTIiEaOxIcZfvrCXlvZu6qpL\nw44kImk21kg6VaZNm3b5ejQavbzLJVEzZsxIdqQRZc0IHeALN1VRPC1Po3QRCd2sWbMoKyu7vB/8\n8ccfvzxav5qSkhLOnj2bkkxZVejF0/L4wsoFvLSnk66zF8OOIyJTXFNTE3/2Z3/G8uXLaWlp4cEH\nHxzzMZ/5zGd47rnnUvKhqMW/E5Qe9fX1PtkTXBw8do71f/Ov/MltH+IP1y1OUjIRyVT79u3j+uuv\nDztG2oz07zWzN929fqzHZtUIHWDR3GI+tricJ7Yfom9gMOw4IiIZI+sKHaCxIcbRMxf56d6jY68s\nIjJFZGWh37JkLtWzi/ThqIjIEFlZ6NGIcf+aWna0nkz6lwxERLJVVhY6wN311RTmR3hsa2vYUURE\nMkLWFnrp9AI21C3gRy2H6T7fG3YcEZHQZW2hQ/z4Lj19g/yguT3sKCIyhXzta18LO8KIsrrQr6+c\nyaqFs3ls6yEGBtM3n15Epi5356//+q/DjjGirC50iE9h7Dh1gVf3j3h+DRGRSWttbeW6665j48aN\nLFq0iAsXLlBXV8e9994bdrQrZM3BuUZz+7J5zJ9ZSNPWVtYv1WlNRXLaTx6AI28l9znnfxg+Nfbp\nHA4cOEBTUxOPPfYYxcXFtLS0JDdHEmT9CD0/GuHe1TW8fuA4B4+dCzuOiOSo2tpa1qxZE3aMq8r6\nETrAPatr+Na/HOTxra08dNcNYccRkVRJYCSdKuk6BO5kZP0IHaC8eBqfXl7JM292cLanL+w4IpLj\n8vPz6evLvK7JiUIH2Lg2xvu9Azy783DYUUQkx23atInly5dn3IeiOVPoddWl3FhdStPWVgY1hVFE\nkigWi/H2229fvv31r3+dffv28eSTT4aY6oMSKnQzKzWzZ8xsv5ntM7MGM5ttZq+Y2YHgsizVYcfS\n2FDLe13v8/ODx8OOIiKSdomO0L8JvOzuS4AbgX3AA8Bmd18MbA5uh+rO5ZWUFxfo+C4iMiWNWehm\nNgv4OPBdAHfvdfdu4C6gKVitCdiQqpCJmpYX5Z5VNWzef4y2E+fDjiMiSZLOM6uFabL/zkRG6AuB\nLuD/mdkuM/uOmc0A5rl7Z7DOEWDEb/WY2SYzazaz5q6urkmFTcTvra4hYsYT2w+l/LVEJPUKCws5\nceJEzpe6u3PixAkKCwsn/ByJzEPPA1YCf+ju283smwzbveLubmYjbm13fwR4BOLnFJ1w0gRVziri\njmXz+f4b7fzn9R+iqCCa6pcUkRSqqqqio6ODdAwIw1ZYWEhVVdWEH59IoXcAHe6+Pbj9DPFCP2pm\nle7eaWaVQMYcTGVjQy0vvdXJj1oOc8+qmrDjiMgk5Ofns3DhwrBjZIUxd7m4+xGg3cyuCxatA94B\nXgAag2WNwPMpSTgBqxbOZsn8Epq2tOb8n2kiIpckOsvlD4EnzWwPUAd8DXgYuM3MDgDrg9sZwcz4\nytoY+4+cZcevToYdR0QkLRIqdHdvcfd6d1/u7hvc/ZS7n3D3de6+2N3Xu3tGNedddQuYVZTPY1v1\n4aiITA05803R4YoKonzpI9W8vPcInacvhB1HRCTlcrbQAe5bXcugO09uaws7iohIyuV0odfMmc66\nJXN5akcbF/sHwo4jIpJSOV3oED+R9In3e3lpT+fYK4uIZLGcL/Sbf6ucaytm0KQPR0Ukx+V8oUci\nRmNDjN3t3bS0d4cdR0QkZXK+0AE+v3IBMwqiNG1pDTuKiEjKTIlCLynM54s3VfHSnk66zl4MO46I\nSEpMiUIHuL8hRu/AIE/v0BRGEclNU6bQF80t5mOLy3lyext9A4NhxxERSbopU+gAjQ0xjpzp4ad7\nj4YdRUQk6aZUod+yZC5VZUU0bW0NO4qISNJNqUKPRoyNDbXs+NVJ9nWeCTuOiEhSTalCB7i7vprC\n/IimMIpIzplyhV46vYANdQv4Ucthus/3hh1HRCRpplyhQ/z4Lj19g/yguT3sKCIiSTMlC/36ypms\nWjibx7cdYmBQp6gTkdwwJQsd4lMY209e4NX9GXNuaxGRSZmyhX77snnMn1moKYwikjMSKnQzazWz\nt8ysxcyag2WzzewVMzsQXJalNmpy5Ucj3Lu6htcPHOfgsXNhxxERmbTxjNBvcfc6d68Pbj8AbHb3\nxcDm4HZWuWd1DQXRCI9vbQ07iojIpE1ml8tdQFNwvQnYMPk46VVePI1PL6/kmTc7ONvTF3YcEZFJ\nSbTQHfhnM3vTzDYFy+a5+6Xzuh0B5iU9XRpsXBvj/d4Bnt15OOwoIiKTkmihf9Td64BPAX9gZh8f\neqe7O/HS/wAz22RmzWbW3NXVNbm0KVBXXcqN1aU0bW1lUFMYRSSLJVTo7n44uDwGPAesAo6aWSVA\ncDni/D93f8Td6929vqKiIjmpk6yxoZb3ut7n3949HnYUEZEJG7PQzWyGmZVcug7cDrwNvAA0Bqs1\nAs+nKmSq3bm8kvLiAh3fRUSyWiIj9HnAz81sN7ADeMndXwYeBm4zswPA+uB2VpqWF+WeVTVs3n+M\n9pPnw44jIjIhYxa6u7/n7jcGP8vc/X8Ey0+4+zp3X+zu6939ZOrjps7vra4hYsbj2w6FHUVEZEKm\n7DdFh6ucVcQdy+bz/TfaudA7EHYcEZFxU6EPsbGhltMX+ni+RVMYRST7qNCHWLVwNkvml/Dollbi\nMzFFRLKHCn0IM+Mra2PsP3KWN1pPhR1HRGRcVOjD3FW3gFlF+ZrCKCJZR4U+TFFBlC99pJqX9x6h\n8/SFsOOIiCRMhT6C+1bXMujO97a3hR1FRCRhKvQR1MyZzrolc3lqRxsX+zWFUUSygwp9FI1rYxw/\n18uP3+oce2URkQygQh/Fzb9VzrUVM3h0i745KiLZQYU+ikjEaGyIsbu9m5b27rDjiIiMSYV+FZ9f\nuYAZBVEe0xRGEckCKvSrKCnM54s3VfHink6On7sYdhwRkatSoY/h/oYYvQODPL1DUxhFJLOp0Mew\naG4xH1tczhPb2ugbGAw7jojIqFToCWhsiHHkTA8/3Xs07CgiIqNSoSfgliVzqSoromlra9hRRERG\npUJPQDRibGyoZcevTrKv80zYcURERqRCT9Dd9dUU5kd4bGtr2FFEREakQk9Q6fQCNtQt4Lldh+k+\n3xt2HBGRD0i40M0sama7zOzF4PZsM3vFzA4El2Wpi5kZNjbE6Okb5AfN7WFHERH5gPGM0P8I2Dfk\n9gPAZndfDGwObue0pdfMZFVsNo9vO8TAoE5RJyKZJaFCN7Mq4E7gO0MW3wU0BdebgA3JjZaZGtfG\naD95gVf3Hws7iojIFRIdoX8D+HNg6Ddr5rn7pWPLHgHmjfRAM9tkZs1m1tzV1TXxpBni9mXzmD+z\nUFMYRSTjjFnoZvZp4Ji7vznaOu7uwIj7INz9EXevd/f6ioqKiSfNEPnRCPeuruH1A8d5t+tc2HFE\nRC5LZIR+M/BZM2sFngZuNbMngKNmVgkQXE6ZfRD3rK6hIBrRURhFJKOMWeju/hfuXuXuMeDLwL+4\n+33AC0BjsFoj8HzKUmaY8uJp3Lm8kmfe7OBsT1/YcUREgMnNQ38YuM3MDgDrg9tTRuPaGO/3DvDs\nzsNhRxERAcZZ6O7+mrt/Orh+wt3Xuftid1/v7idTEzEz1VWXcmN1KU1bW4l/hCAiEi59U3QSGhtq\nea/rfX5+8HjYUUREVOiTcefySubMKKBJJ5IWkQygQp+EaXlR7llVw+b9R2k/eT7sOCIyxanQJ+ne\nNTVEzHh8m0bpIhIuFfokVc4q4o5l8/n+G+1c6B0IO46ITGEq9CTY2FDL6Qt9PN+iKYwiEh4VehKs\nWjibJfNLaNp6SFMYRSQ0KvQkMDMa18bY13mGN1pPhR1HRKYoFXqSbKhbwKyifJp0fBcRCYkKPUmK\nCqJ86SPVvLz3CEdO94QdR0SmIBV6Et23upZBd57crimMIpJ+KvQkqpkznXVL5vLUjjYu9msKo4ik\nlwo9yTY2xDh+rpcfv9U59soiIkmkQk+yjy4q59qKGTyq47uISJqp0JMsEjEaG2Lsbu+mpb077Dgi\nMoWo0FPg8ysXMKMgqlPUiUhaqdBToKQwny/eVMWLezo5fu5i2HFEZIpQoafI/Q0xegcGeXpHW9hR\nRGSKUKGnyKK5xXxscTlPbGujb2Aw7DgiMgWMWehmVmhmO8xst5ntNbOHguWzzewVMzsQXJalPm52\naWyIceRMD6+8czTsKCIyBSQyQr8I3OruNwJ1wB1mtgZ4ANjs7ouBzcFtGeKWJXOpKiviUX04KiJp\nMGahe9y54GZ+8OPAXUBTsLwJ2JCShFksGjE2NtSy41cn2dd5Juw4IpLjEtqHbmZRM2sBjgGvuPt2\nYJ67X/o65BFgXooyZrW766spzI/w2NbWsKOISI5LqNDdfcDd64AqYJWZ3TDsfic+av8AM9tkZs1m\n1tzV1TXpwNmmdHoBG+oW8Nyuw3Sf7w07jojksHHNcnH3buBV4A7gqJlVAgSXx0Z5zCPuXu/u9RUV\nFZPNm5U2NsTo6RvkH5o7wo4iIjkskVkuFWZWGlwvAm4D9gMvAI3Bao3A86kKme2WXjOTVbHZPLat\nlYFBnaJORFIjkRF6JfCqme0B3iC+D/1F4GHgNjM7AKwPbssoGtfGaD95gdd+MeIfMiIik5Y31gru\nvgdYMcLyE8C6VITKRbcvm8f8mYU8uqWVddfr82MRST59UzRN8qMR7l1dw+sHjvNu17mxHyAiMk4q\n9DS6Z3UNBdEIj2/VsdJFJPlU6GlUXjyNO5dX8sybHZy72B92HBHJMSr0NGtcG+PcxX6e3akpjCKS\nXCr0NKurLuXGqlk0bWkl/n0sEZHkUKGHoHFtjHe73ufnB4+HHUVEcogKPQR3Lq9kzowCmnQiaRFJ\nIhV6CKblRblnVQ2b9x+l/eT5sOOISI5QoYfk3jU1RMx4YptG6SKSHCr0kFTOKuK3l83j6TfaudA7\nEHYcEckBKvQQNTbEOH2hj+dbDocdRURygAo9RKsWzmbJ/BKath7SFEYRmTQVeojMjMa1MfZ1nuGN\n1lNhxxGRLKdCD9mGugXMLMyjaWtr2FFEJMup0ENWVBDlSx+p5uW3j3DkdE/YcUQki6nQM8D9a2IM\nuvPkdk1hFJGJU6FngJo501m3ZC5P7WjjYr+mMIrIxKjQM8TGhhjHz/Xy47c6w44iIllKhZ4hPrqo\nnGsrZuj4LiIyYWMWuplVm9mrZvaOme01sz8Kls82s1fM7EBwWZb6uLkrEjE2rqmlpb2b3e3dYccR\nkSyUyAi9H/gTd18KrAH+wMyWAg8Am919MbA5uC2T8IWbqphRENUURhGZkDEL3d073X1ncP0ssA9Y\nANwFNAWrNQEbUhVyqigpzOeLN1Xx4u5Ojp+7GHYcEcky49qHbmYxYAWwHZjn7pc+wTsCzEtqsinq\n/oYYvQODPL2jLewoIpJlEi50MysGfgj8sbufGXqfxw9EMuLBSMxsk5k1m1lzV1fXpMJOBYvmFvOx\nxeU8sa2N/oHBsOOISBZJqNDNLJ94mT/p7s8Gi4+aWWVwfyVwbKTHuvsj7l7v7vUVFRXJyJzzNjbE\nOHKmh5++czTsKCKSRRKZ5WLAd4F97v43Q+56AWgMrjcCzyc/3tR065K5VJUV0bSlNewoIpJFEhmh\n3wzcD9xqZi3Bz+8ADwO3mdkBYH1wW5IgGjE2NtSy/Vcn2dd5ZuwHiIiQ2CyXn7u7uftyd68Lfn7s\n7ifcfZ27L3b39e5+Mh2Bp4q766spzI/w2NbWsKOISJbQN0UzVOn0AjbULeC5XYc5fb4v7DgikgVU\n6BlsY0OMnr5BftDcHnYUEckCKvQMtvSamayKzebxbYcYGNQp6kTk6lToGa5xbYy2k+d57RcjzgoV\nEblMhZ7hbl82j/kzC3lUUxhFZAwq9AyXH41w7+oaXj9wnHe7zoUdR0QymAo9C3x5VQ0F0QiPb9Wx\n0kVkdCr0LFBRMo07l1fyzJsdnLvYH3YcEclQKvQs0bg2xrmL/Ty7syPsKCKSoVToWaKuupQbq2bR\ntKWV+MEtRUSupELPIo1rY7zb9T7/dvBE2FFEJAOp0LPIncsrmTOjQFMYRWREKvQsMi0vyj2rati8\n/yjtJ8+HHUdEMowKPcvcu6aGiBlPbNMURhG5kgo9y1TOKuK3l83j6TfaudA7EHYcEckgKvQs1NgQ\n4/SFPl7YfTjsKCKSQVToWWjVwtksmV/Co1sOaQqjiFymQs9CZkbj2hj7Os/QfOhU2HFEJEOo0LPU\nhroFzCzM0xRGEblMhZ6ligqifOkj1bz89hGOnO4JO46IZIAxC93M/t7MjpnZ20OWzTazV8zsQHBZ\nltqYMpL718QYdOd72zWFUUQSG6E/CtwxbNkDwGZ3XwxsDm5LmtXMmc6t183lezvauNivKYwiU92Y\nhe7uPwNODlt8F9AUXG8CNiQ5lySocW2M4+d6+clbR8KOIiIhm+g+9Hnu3hlcPwLMS1IeGaePLirn\n2ooZ+nBURCb/oajHJ0KPOhnazDaZWbOZNXd1dU325WSYSMTYuKaWlvZudrd3hx1HREI00UI/amaV\nAMHlqKekd/dH3L3e3esrKiom+HJyNV+4qYoZBVGatraGHUVEQjTRQn8BaAyuNwLPJyeOTERJYT5f\nuKmKF3d3cvzcxbDjiEhIEpm2+BSwFbjOzDrM7N8DDwO3mdkBYH1wW0K0sSFG78Ag33+jPewoIhKS\nvLFWcPd7RrlrXZKzyCQsmlvMxxaX88S2Q/yHj19LXlTfGROZavS/PodsbIjRebqHn75zNOwoIhIC\nFXoOuXXJXKrKimjSFEaRKUmFnkOiEeP+NbVs/9VJ9nWeCTuOiKSZCj3HfOkj1RTmRzRKF5mCxvxQ\nVLJL6fQCPrdiAU/taOdnv+xiRU0ZK2pKWVFTxg0LZjItLxp2RBFJERV6Dvrvdy5l8dwSdradYldb\nNy+9FT9KQ0E0wtJrZrKippSVQdEvKC3CzEJOLCLJYOk8hVl9fb03Nzen7fUk7tiZHna2dbOr/RS7\nDnWz53A3PX2DAFSUTGNlMIJfUV3K8qpSigo0ihfJJGb2prvXj7WeRuhTwNyZhdxxw3zuuGE+AH0D\ng/ziyNnLI/idbaf4p73xqY7RiHF9ZcnlEfyK6jJq50zXKF4kC2iELgCcOHeRlvbuyyW/u72b93vj\nx1ifPaOAFdWll3fVLK8upXiaxgIi6aIRuozLnOJprLt+Huuujx8JeWDQ+eXRs+xq62ZX2yl2tp1i\n8/74MdgiBh+aV3L5A9eVNWVcWz6DSESjeJEwaYQuCTt9vi++H76tm13t8aI/29MPwMzCPOpqyi7v\nj6+rLmVWUX7IiUVyg0boknSzpufzyevm8snr5gIwOOi8d/xc/APXYFfNNzcf4NIYYdHcYlZUl7Ky\nNj6SXzy3hKhG8SIpoxG6JNXZnj72dJwOdtPEi/7U+T4AiqflcWP1LFZUl7GytpS66jJmzygIObFI\n5tMIXUJRUpjPzYvKuXlROQDuzqET5y9/2Lqr/RT/91/fZWAwPpCIzZn+mxk1NWUsmV+iI0WKTJBG\n6JJ253v7eavj9OUR/M627ssn5ijKj/LhqllDSr6UuSWFIScWCZdG6JKxphfksfraOay+dg4QH8V3\nnLpw+YPWnW3dfPfn79E3EB9sVJUVXf7i08raMpZWzqQgT6N4keFU6BI6M6N69nSqZ0/nszdeA0BP\n3wB7f306mDbZTXPrSf5x968BKMiL8OEFs4K58fH98ZWzisL8J4hkBO1ykazRefoCLW2/+fLTnsOn\n6e2PH8Jg/szCK45Rc8OCWRTm6xAGkhu0y0VyTuWsIio/XMSnPlwJQG//IPs6z/xmRk37KX7y9hEA\n8qPG0sqZV3z5qapMByKT3DapEbqZ3QF8E4gC33H3q54sWiN0SbWusxfjc+KD/fG7209zoS9+CIPy\n4mmXP2hdWVPG8qpZTC/QmEYyX6Ij9AkXuplFgV8CtwEdwBvAPe7+zmiPmXChD8TnMRPJA42wZBz6\nBwb5RXAIg51tp2hp6+a94+8D8QORXTevhJW1pSytnEVhfoS8aIT8iJEXjZAXNfIjwWXUyItEiEaM\n/GH3XblehLyIEY2Y/hqQpEnHLpdVwEF3fy94waeBu4BRC33CXn4A3vhOcMPixR7Jg0g0/mPRMZZF\ngjeDqy0LHnfpfotOYllkyHMmc9mQ5SMuU4EMlxeNsOyaWSy7Zhb3rakF4NT7vbQMmVHzo12/5olt\nbcl/7YgNK/6R3yzyLi0b8maRF4nE30QuP+Y36/zmvivXyxv2ZhO99PpXee5Ly+NvVFfPqDeozDeZ\nQl8AtA+53QGsnlycUXzoU1A8H3wABvthcMjlVZddWj582SD094JfuHLZYP+w5xu4+jIfTMk/d+Ls\nyjeBK/4DDrk+0vIr/q+OsW5Yy0ddl1GWj7x+mRm3ALcEy31O/GBkgwAOl/5mdQcPFjjBcvf49WCl\n+HUfcj1+zYc9z+Vlbni/431DX2Pkx1zx2peeY9hrM8pf2M74y7cv+Lkw2goW34p26cpIK6TBmK+S\nohiTfdrzt/9PljZ8KilZRpPyHYhmtgnYBFBTUzOxJ1m8Pv6TaQYHg7IfWvzD3mASXTbYH3+DmNCy\nEd7Yhr7ZXPGf3kdYPtKyTFme+tc0II+RizGpUjCjzIFBj78ZXLp0H2TQ4280g8Ebw9D7By8tv3wZ\nX+fKZR98jDsMXn5juvq/xUfanmP88ye0dfyqNyfyFCM/SxKyz5xRmligSZhMoR8GqofcrgqWXcHd\nHwEegfg+9Em8XuaJRIAIRHVUQQmHEZ+RIAIwma/bvQEsNrOFZlYAfBl4ITmxRERkvCY8Qnf3fjP7\nT8A/ER8k/L27701aMhERGZdJ7UN39x8DP05SFhERmQQd4UhEJEeo0EVEcoQKXUQkR6jQRURyhApd\nRCRHpPV46GbWBRya4MPLgeNJjJMsyjU+yjU+yjU+mZoLJpet1t0rxloprYU+GWbWnMjRxtJNucZH\nucZHucYnU3NBerJpl4uISI5QoYuI5IhsKvRHwg4wCuUaH+UaH+Uan0zNBWnIljX70EVE5OqyaYQu\nIiJXkXGFbmZ3mNkvzOygmT0wwv1mZv87uH+Pma3MkFyfNLPTZtYS/DyYhkx/b2bHzOztUe4Pa1uN\nlSvt2yp43Woze9XM3jGzvWb2RyOsk/ZtlmCuMH6/Cs1sh5ntDnI9NMI6YWyvRHKF8jsWvHbUzHaZ\n2Ysj3Jfa7RU/M0lm/BA/DO+7wLVAAbAbWDpsnd8BfkL82P5rgO0ZkuuTwItp3l4fB1YCb49yf9q3\nVYK50r6tgtetBFYG10uIn+Q8E36/EskVxu+XAcXB9XxgO7AmA7ZXIrlC+R0LXvu/AN8b6fVTvb0y\nbYR++cTT7t4LXDrx9FB3AY953Dag1MwqMyBX2rn7z4CTV1kljG2VSK5QuHunu+8Mrp8F9hE/N+5Q\nad9mCeZKu2AbnAtu5gc/wz90C2N7JZIrFGZWBdwJfGeUVVK6vTKt0Ec68fTwX+xE1gkjF8Da4M+o\nn5jZshRnSkQY2ypRoW4rM4sBK4iP7oYKdZtdJReEsM2C3QctwDHgFXfPiO2VQC4I53fsG8CfA6Od\nQT6l2yvTCj2b7QRq3H058C3gRyHnyWShbiszKwZ+CPyxu59J52tfzRi5Qtlm7j7g7nXEzxm8ysxu\nSMfrjiVs/mYdAAABhklEQVSBXGnfXmb2aeCYu7+Z6tcaTaYVeiInnk7o5NTpzuXuZy79GejxMznl\nm1l5inONJYxtNaYwt5WZ5RMvzSfd/dkRVgllm42VK+zfL3fvBl4F7hh2V6i/Y6PlCml73Qx81sxa\nie+WvdXMnhi2Tkq3V6YVeiInnn4B2Bh8WrwGOO3unWHnMrP5ZmbB9VXEt+2JFOcaSxjbakxhbavg\nNb8L7HP3vxlltbRvs0RyhbHNzKzCzEqD60XAbcD+YauFsb3GzBXG9nL3v3D3KnePEe+If3H3+4at\nltLtNalziiabj3LiaTP7anD/t4mfw/R3gIPAeeDfZUiuLwL/0cz6gQvAlz34WDtVzOwp4p/ml5tZ\nB/CXxD8gCm1bJZgr7dsqcDNwP/BWsP8V4L8BNUOyhbHNEskVxjarBJrMLEq8EH/g7i+G/f8xwVxh\n/Y59QDq3l74pKiKSIzJtl4uIiEyQCl1EJEeo0EVEcoQKXUQkR6jQRURyhApdRCRHqNBFRHKECl1E\nJEf8fx8SP1GweOjaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1dd314a0dd8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(vars1)\n",
    "plt.plot(vars2)\n",
    "plt.legend(['no rt', 'rt'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Variance of the estimates using reparameterization trick is one order of magnitude smaller than the estimates from the first method!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are some additional resources I found interesting:\n",
    "\n",
    "<ol>\n",
    "  <li>Section 3.1 from [Yarin Gal's Thesis](http://mlg.eng.cam.ac.uk/yarin/thesis/thesis.pdf) contains a different and more rigorous derivation</li>\n",
    "  <li>[Miller's blog](http://andymiller.github.io/2016/12/19/elbo-gradient-estimators.html) reconciles various estimators like score fucntion estimators, pathwise estimators, etc. It was here I came across the relation between Control variates in linear regression and variance reduction in pathwise estimators. Mind = Blown </li>\n",
    "  <li> I am looking for a more rigorous formulation of the above relation between reparameterization trick and control variates in linear regression. If you have come across any, do let me know. you can mail me at kasa at the rate u dot nus.edu <li>\n",
    "  \n",
    "</ol>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.1 |Anaconda 4.4.0 (64-bit)| (default, May 11 2017, 13:25:24) [MSC v.1900 64 bit (AMD64)]\n",
      "3.6.1\n",
      "The Python version is 3.6.1\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "print(sys.version)\n",
    "import platform\n",
    "print(platform.python_version())\n",
    "print(\"The Python version is %s.%s.%s\" % sys.version_info[:3])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}