kenjisato/2015-07-06.ipynb

## 2015-07-06.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Policy Iteration\n",
    "\n",
    "Stachurski (2008) Chapter 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.interpolate import interp1d\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## ラムゼーモデル\n",
    "\n",
    "再掲"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Ramsey:\n",
    "    \"\"\"One-sector Ramsey model\"\"\"\n",
    "\n",
    "    def __init__(self, A, α, ρ):\n",
    "        self.A, self.α, self.ρ = A, α, ρ\n",
    "\n",
    "    def f(self, x):\n",
    "        \"\"\"Production function\"\"\"\n",
    "        return self.A * x ** self.α\n",
    "\n",
    "    def U(self, x):\n",
    "        \"\"\"Utility from consumption\"\"\"\n",
    "        return np.log(x)\n",
    "    \n",
    "    def u(self, x, y):\n",
    "        \"\"\"Reduced form utility\"\"\"\n",
    "        return self.U(self.f(x) - y)\n",
    "\n",
    "    def is_feasible(self, x, y):\n",
    "        \"\"\"Checks feasibility\"\"\"\n",
    "        return self.f(x) >= y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 解析解\n",
    "\n",
    "再掲"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def analytic_solutions(A, α, ρ):\n",
    "    c0 = (1.0 / (1.0 - α * ρ) / (1.0 - ρ) * np.log(A) +\n",
    "          α * ρ * np.log(α * ρ) / (1.0 - α * ρ) / (1.0 - ρ) +\n",
    "          np.log(1.0 - α * ρ) / (1.0 - ρ))\n",
    "    c1 = α / (1.0 - α * ρ) \n",
    "    \n",
    "    def value_func(x):\n",
    "        return c0 + c1 * np.log(x)\n",
    "    \n",
    "    def policy_func(x):\n",
    "        return ρ * c1 * A * x**α / (1.0 + ρ * c1)  ## Fixed this line\n",
    "        \n",
    "    return value_func, policy_func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 線形近似\n",
    "\n",
    "前回より実用的に ... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class PLApprox:  # Refactored this class\n",
    "    def __init__(self, a, b, N, upsample=10):\n",
    "        self.a, self.b, self.N = a, b, N\n",
    "        self.grids = np.linspace(a, b, upsample*N)\n",
    "        self.centers = np.linspace(a, b, N)\n",
    "    \n",
    "    def proj(self, f):\n",
    "        return np.array([f(c) for c in self.centers])\n",
    "    \n",
    "    def inj(self, a):\n",
    "        return interp1d(self.centers, a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Policy Iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "任意の $x \\in X$, に対して $(x, f(x)) \\in \\mathbb{D}$ を満たす関数 $f$ を policy function という. \n",
    "\n",
    "Policy function $f$ と十分大きな $T$ について, その policy function の value\n",
    "\n",
    "$$\n",
    "    v_f(x) = \\sum_{t=0}^\\infty \\rho^t u(f^t(x), f^{t+1}(x)) \n",
    "    \\approx \\sum_{t=0}^T \\rho^t u(f^t(x), f^{t+1}(x))\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def value_of_policy(f, model, apx, T):\n",
    "    fc = apx.inj(f)\n",
    "    X = apx.centers\n",
    "    Y = fc(X)\n",
    "    vf = np.zeros_like(X)\n",
    "    for i in range(T):\n",
    "        vf += model.ρ**i * model.u(X, Y)\n",
    "        X, Y = Y, fc(Y)\n",
    "    return vf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Policy の更新則: $f\\to v_f \\to \\hat f$ \n",
    "\n",
    "$$\n",
    "  \\hat f (x) = \\arg\\max_y u(x, y) + \\rho v_f(y)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def greedy(vf, model, apx):\n",
    "    vc = apx.inj(vf)\n",
    "    Y = apx.grids\n",
    "    fhat = np.empty_like(vf)\n",
    "    for i, x in enumerate(apx.centers):\n",
    "        X = np.ones_like(Y) * x\n",
    "        with np.errstate(invalid='ignore'):\n",
    "            maximand = np.where(model.is_feasible(X, Y),\n",
    "                                model.u(X, Y) + model.ρ*vc(Y),\n",
    "                                -np.inf)\n",
    "        fhat[i] = Y[np.argmax(maximand)]\n",
    "    return fhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Policy Iteration Algorithm\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{array}{cccccc}\n",
    "f_0        &          & f_1        &          & f_2        & \\cdots & \\to & f^* \\\\\n",
    "\\downarrow & \\nearrow & \\downarrow & \\nearrow & \\downarrow &        &     &     \\\\\n",
    "v_0        &          & v_1        &          & v_2        & \\cdots & \\to & v^* \\\\\n",
    "\\end{array}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "いつものモデルで実験"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A, α, ρ =1.1, 0.4, 0.9\n",
    "ramsey = Ramsey(A, α, ρ)\n",
    "apx = PLApprox(a=1e-3, b=5.0, N=500, upsample=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1106df588>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEACAYAAACnJV25AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHDNJREFUeJzt3Xl0VeW9//H3l3kmQCABAglDmFEEFFCGg/NAUVSqOFFt\n1Vat9opDS68l99dae2urrdOtYqtiLSpwURQRETgIiiAyhTGEMWHIAGQAEsjw3D9IkZ8GzLxP9vm8\n1spynz0937MXftaznvPsvc05h4iI+EcdrwsQEZGqpWAXEfEZBbuIiM8o2EVEfEbBLiLiMwp2ERGf\nqVcdJzWzt4GeJR8jgCzn3HnV0ZaIiPz/qiXYnXM3/3vZzP4EZFVHOyIi8l1WnTcomZkBu4HRzrnt\n1daQiIicUt1j7COANIW6iEjNqfBQjJktAKJL2TTZOfdByfIE4F8VbUNERMqv2oZizKwekAoMdM7t\nO8M+elCNiEg5OefsbNurcyjmUmDzmUL935xz+nOOKVOmeF5DKPzpOuha6Fqc/a8sqjPYbwKmV+P5\nRUSkFNUy3RHAOXdndZ1bRETOTHeehohAIOB1CSFB1+Ebuhbf0LUon2qdx/69jZs5L9sXEaltzAzn\n4Y+nIiLiAQW7iIjPKNhFRHxGwS4i4jMKdhERn1Gwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iIzyjY\nRUR8RsEuIuIzCnYREZ9RsIuI+IyCXUTEZ6rtDUoiIlJ1jh3NY/v2PWXaV8EuIlIDDmUcZvv2PeSl\n7iNr7z4yDmaQXFxI96QtWM4hsuo41vSM54aP59A0/xjpkW2Z/oMb+NvTT9LyRDEbu/bkkZ/9tExt\n6Q1KIiJlUFRYxN60gxxK3k3Grj2kHTjAxoI84rdswWVlklt8ghV9enHT3Nk0zT9KdkQEU394O288\nmUBEfjGb4rpw76RJ/Gvyf5DdsAFJHTsz/epruf392Rxv0pxDbduzvt85DN9/kPqtI7HoaPI7dmRA\nVBTRXTsR1SmaevXrlekNSgp2EQkbxUXFZO7PZN/23ezdtYfNWYeI2rmTgsx0cvNz+SquE1cv/oQm\nx3I51rghr42/hef/9CQRxwvJbNOBCVMS+OTBezjcsD7J0dFMveEW7pr1DnmNTwbz6nMHMnxfBg0i\no6gb3YG8zjGcH9ORjvFxRERGUKdu5X/WVLCLiG8dPpjN6k3bKN6xm+yUFA6nHWBli4YMW/kl9XMP\nUuQKefeqMfzmleeIOH6c/KYtufm3T7P2xz/icKM6JEW14z/v+wVPvPgMRxs3I7NVOz4dOZrLN2zB\nIiKhXXsy42IZHNGGtnGd6BgfR+uo1l5/bQW7iIS+woJC9iansDtpB2sO7CNy+w5OpO0nL+cgi3vH\nc8O8OTQ5lkMdc/zXA48wc/Ik2uQXkdu0BRc/+yyzJt3H4UYNORDRhn+Mv5WbPvmEwhZtKIhsz5a+\nfRlOfVp16kyb2BgiunUmtkNbr79ypXgW7GZ2AfACUB8oBO5zzn1Vyn4KdhGfycrMYuf6raQmJfPF\n0Sx6JG6gOHM/lpPJzItH8+CbrxKRd5TGxY7rnvkb2267jaMNYG+LJkz43Z945vdPkNu4GTktWvPe\nFWO4Zv1m6kZG0yCqAwe7dWVIdHs69upK+9gOVTK0Udt4GexB4Cnn3Hwzuwp4zDk3upT9FOwiIS49\nNZ1d67dwIHk7S/Nz6bV+PRzcT/2sDF4fO5bHX32B1kePEHGigEFv/i8Hxo3lUEPIaNKIn07+LY+/\n/Dx5TSM43jKSRSMCXJaZQ7OOnWgdFwdxnbigT3eatmjq9desNcoS7NU13XE/0LJkOQLYW03tiEg5\nncg/wfb1SaRu2kLwYDpdEhOx9BQaHDrAW1dfw0PTXiEqN4c2ecc5/4132XjHrURZAfUaN+SLx6cw\nOHEpBY2ac6JlJF3zitl/7U/Ii+lMQbdufN0znkZ5BcTUrUMMsALgoYdPtX2XV186zFRXjz0WWAY4\nTt7dOsw5l1LKfuqxi1SBwoJCdmzYxtJNW2iUlMyJlJ0UH0jhw4HncOO892ifmUbksaNc/+wrfPbg\nz6mfn0N6kwbc/9hvuHv6m7gGTTnRqh3Lh17EhflFtI2JJSq+G5G9uxPXKQqzs3YQpQZVa4/dzBYA\n0aVs+jXwIPCgc262mY0H/gFcVtp5EhISTi0HAgECgUBFSxLxndT9mexavYGMTZvJ2bGNJS0bcf7K\nz2mXvpeII1n89oFHePqFvzAgJY3WDY1//P5pHp79N5qdKOJI89a0KOhL1jkjsXbtOdE9ntdje9Bi\n116atWhKWyAI8NjkU+2pRx16gsEgwWCwXMdUV489xznXomTZgCznXMtS9lOPXcLSkewjbP0qkaXJ\n22i8ZTP1diVTL2Mv748azvj5czhn9y7aHz3BFc++yO9feJoWhzLIbNqSaeNuZkjSdiKaRNAwJo6c\nnj0ZFN+NPuf2plnLZl5/LakBXo6xJ5vZKOfcEuBiIKma2hEJKUWFRexJ2k3yytWkb9jAl3acrpvW\nEbMnmVZHsnhu4t38fPZMRm1Iol3TugQfepQrV39K09w88lpHE1tYn5yrf8S+bt1peF5/VvSLp85P\nf3bq/Nd4+N2k9qiuHvtg4EWgIZDHyemOa0rZTz12qVW270glacVqshITyduRxBdREfTdtJYBSRto\nf+QIv3g0gXvmzKbX1vXsb96MaWNvoue+NKLrNqFhp64c79eHc3v3pv+5PWnQqIHXX0dqId2gJFIO\n+cfyWbz0Kw6sWYvbsA7bt5PlPbrSL3kTl6z9mo5HjnP/pF8yfNUX9N6cyMEWbVh80cV0dPWJbR1N\ndL/+xA0ZQNf4WK+/iviYgl3kNPnH8lny+SpS16zDEtdSJ2Uby3t2pXvqbsYvW0L7o0VMuvduOuzd\nQ++kJHJat2fjwKFEtYykX5du9BgyiNheXcLyphgJHQp2CSuFBYUsW7qK7avXYOvXYCnJrIttT5uc\nw9y2cAEdjhaScNst1M07yqANiWS36ciucwcTGdOFwf3602/4YP0AKSFPwS6+k5NzlJVLVpCxfDnH\nNq9jc/P61CvMZ+L8D4nLOcEzN4wjs0ljRn65nOzIjuzpP4hmnbsxrE9f+o+8QMEttZ6CXWql3Kxc\nvv54CSkrlpOcnU5uUR5jgh/TOTub+aOuZH3nGCa8/y5praPY2v986NKLYd16cO4lI2gX087r8kWq\nlYJdQlZxUTHJ67eyZNESklN30G/5IqLSU9jRvSfBISP57Qt/ZFdEBKv6DWZvj74Mq9uMzkOGct6l\nF6nXLWFNwS6ey805ysL5S8hdupSCLWvIIp9V5w7ghReeA2DugH7MHzqKizduo36v84i68CIGXTGS\nNq1aeFy5SGhSsEuNKSwoZNUny1i9JMjK4qOM+3AGXQ+mkdeqA/9x3/1Mee6/SWvXmeweAzh6wTDG\nXTSM7uf01AwTkXJSsEuVKyosYumyr9n/yQKOrf2Supl7ePXm2/j4l49xsFFdVsfE8M9xExiTepj2\nw0dxwdgrQuKtMyJ+4eUjBcQH8o/l89mseczYuZWRH84gbv8Oumcd45bX/smbbz3LgXadyOt6LmOa\nRnF4025iu3cmFhjndeEiYU7BLgBk7s/k87ffI3PZQmb368aj015m8N5M4prVZ++Dj5PbtgOpl15P\n2+uvY9/gvnDjeK9LFpEz0FBMGDp0OIfPpr/HwQUf0irpa16c+GP++5Wp1MtJY1tUDJ9edQMjouO4\n5Pof0D6ug9flishpNMYuFBcVM2dekJS5H9B6VZD41GSmPDSZ8fNm0/xwNrm9h1DnqrGMvfZyWmsm\nikjIU7CHoUNph/j79BmcWL6IISsXMyAtkyl33U3bvSl0zoMWwy9n2K3j6dilo9elikgFKNjDwO7N\nO3nrnXc5mrSGHwQ/on9mLn+6dgzZTSMY0LwD5910E/2HD/S6TBGpIgp2H0pas5lZb/2T9Mzd3DV3\nFnE5+bwaGMG27v0ZFRnH6Lsn6rZ6ER9TsPtAemo6M6a+xrqMndw/4w3icvN5d/B5rBh6KWM79eTS\nuybQpHkTr8sUkRqiYK+FCgsKeX/qdObsWMc9b7/KuenZBLt35p/X3srNbWK5/N7bFeQiYUzBXkuk\nbk/hj2+/w/B3pjJyRzIZTRryy/sfZUJefS5/4B4NrYjIKQr2EPbuBwvJePVF+q9ezHnpWdz8q8lc\ntnU3Q2/7CUOvCXhdnoiEKD1SIMQsmD6H1Kl/ZVDiF8y++0H6FR9j94330WvSQ8xVr1xEqoh67NWo\nuKiYBW/MZP/rz7OqTzxtjhylx9YkGt94N2Me+oneUi8i5aYeu0fWLFrBjDdfIb/uCX4+423Sew9g\nRN/hjLv3dho0qO91eSLicwr2KrJzRyov/c//8IO3X6JvZjadB59P1k0/JfblN+iiZ46LSA2qlsQx\ns3PNbLmZrTezOWbWvDra8VpxUTFz/vp3ZvWLp96gHqxs25TtY+6ibsohfrp0Bb984E69SEJEaly1\njLGb2VfAw865pWZ2J9DFOfebUvarlWPsaXsO8OEjj/DGkHOY/NrL7IvpyYW/e5peg/t6XZqI+JyX\nY+zxzrmlJcufAh8D3wn22uZ/Pwqy47nfcdeSRbTvGM0N5wzk4tVbNG4uIiGluoJ9o5ld65x7HxgP\ndKqmdmrEnL/+nTovTmFfz35s6xrPnl89xdWjzve6LBGRUlU42M1sARBdyqbJwF3Ac2b2BDAHOHGm\n8yQkJJxaDgQCBAKBipZUpYqLinnp//2Jz/P28fuXX2RRYCzjn3+RBzqX9pVFRKpHMBgkGAyW65hq\nn8duZj2AN51zQ0rZFnJj7EWFRbzzmz/Q/bWnaFjkePmBJ/jzpIdo3LSx16WJiHj3SAEza+ucyzCz\nOsDrwCLn3Oul7BdSwT7rqed5I38/NyxegOt1Ibc897RuIhKRkOLlj6cTzOz+kuVZpYV6KPls9gKO\nPjqRQRnp7LjjF4yfv4wmjRt6XZaISIWE9SMFdm/eyVuPP8CCwGhun7+Q69+aTkRkhGf1iIh8n7L0\n2MPy7pmiwiJeu+ch6p/fna57dvLkiCu4a/48hbqI+ELY9dhXLficpxbNZuS61cReejPXPXxPjbYv\nIlIZegjYaYqLinnth7dy7bx3GHblOO7413u0imjhdVkiIlUuLIJ9w+pN/Hr2NB77aj6bps3lkRuv\n8rokEZFq4/sx9ref+AOtRvUnJv0gvdftYqRCXUR8zrdj7EWFRfz+jon8eM50lj36R344ZVK1tCMi\nUpPCdlZM2p4DvD58EB+fP4C0D5cr1EUkrPiux75l1UaOjrmA3W2iGR1cRau2rar0/CIiXgq7WTFz\nPgoy76N/cUGPgUxcvEQvuRCRsOSb5Fs0fQ597rga17Qtd362VKEuImHLF0MxS2bOo/udY/jk+ru5\n842/VUFlIiKhybOnO5ZVVQT7J+8vJGF3IhNXbeDeaa9WUWUiIqHJ97NiNq/cQPeJV3Llpu0KdRGR\nErW2x56xL4PdA7uQ2Pt87ly8uIorExEJTb7tsRcXFXPH1Of4aOhIJn660OtyRERCSq3ssf99zHVE\n7klk2KKvaadH7YpIGPHlPPb3nnmFMYvnsHtGUKEuIlKKWhXsSTtSeG/fRrLv+RUTrx7pdTkiIiGp\nVo2xL59wHa3T05n47JNelyIiErJqTY995pN/5ZJNawm8PtPrUkREQlqt6LEfOpxDg5kvsfjHjxPb\nu4vX5YiIhLRaMSsm4Sf3sjw2mvlP/FcNVCUiErp8MY99T9JufjZ9Ko90PsfrUkREaoUKB7uZjTez\njWZWZGYDv7XtV2a2zcy2mNnllSlwyR03sqRbDy6beENlTiMiEjYq02NPBMYBn52+0sz6ADcBfYAr\ngZfMrELtLP8qkfeuvIpeL75ViTJFRMJLhYPdObfFOZdUyqZrgenOuQLn3C4gGbigIm1sfOR++iWu\n45wRgypapohI2KmO6Y4dgC9P+5wKdCzvSfbuSOW6r5aRPGNRlRUmIhIOzhrsZrYAiC5l02Tn3Afl\naOeMU18SEhJOLQcCAQKBAACzfj2Zbh07cM01gXI0IyLiL8FgkGAwWK5jKj3d0cwWA5Occ6tLPv8S\nwDn3h5LPHwNTnHMrSjm21OmOxUXFjP7L01xCC34z6WeVqk9ExE9q8iFgpzcyB/iXmT3DySGYeGBl\neU4259mpvPzUFLrtz62i8kREwkdlpjuOM7MUYCgw18zmATjnNgHvApuAecB95X0274lpz/LF4NHU\nb1C/ouWJiIStkLvzNDvnCPPGXEr3J/7M4Msu8qgyEZHQVCvvPJ32zEvMGDFSoS4iUkEh93THDjOn\ncnWXfl6XISJSa4VUsBcWFDJyx3b2/GWa16WIiNRaIRXs/5g2kyOXXMrDlw7zuhQRkVorpII9beGH\npEd18LoMEZFaLaR+PB21dC5DWnb2ugwRkVotZKY7Zu7PpGFcW45s3Uv7OPXaRURKU5N3nlbaa2/P\n4sCECfxZoS4iUikhE+zuyyAN6xR5XYaISK0XMmPsvdZ+TnyD1l6XISJS64VMj71P2n5SrxrjdRki\nIrVeSPTYV21IJuG+h7nwukq9HlVERAiRYN846z2GfPU5DRo18LoUEZFaLySCvejzT2mVned1GSIi\nvhASwd5q7zZOdD/H6zJERHwhJIJ96s23cnzYCK/LEBHxBc+DvbiomMfefYeRo4Z7XYqIiC94/kiB\njV+uo+3oAbQ9VuxZHSIitUWteIPSluAydrZs6nUZIiK+4fkNSkuOHOTArT/mHa8LERHxCc+Dvdva\nVXQuOu51GSIivuH5UExUyjYiGrbwugwREd/wvMcemZVBQddeXpchIuIbFe6xm9l4M9toZkVmNvC0\n9a3NbLGZ5ZrZ8993nqd+9gsKzhtQ0TJERORbKjMUkwiMAz771vp84D+BR8pykien/o3BfdRjFxGp\nKhUeinHObYGTcyq/tf4Y8LmZxZflPAN3puL6dK9oGSIi8i3V+eNpme58OtTI9FRHEZEqdNYeu5kt\nAKJL2TTZOfdBVRRw4dBR3J6QAEAgECAQCFTFaUVEfCEYDBIMBst1TKUfKWBmi4FJzrnV31o/ERjs\nnPv5WY510wf04OY1WytVg4hIuKjJRwqU1shZG/63xnn5VVSCiIhA5aY7jjOzFGAoMNfM5p22bRfw\nZ+BHZrbHzM447eVos1YVLUFEREpRmVkxs4HZZ9gWV9bzfHTJFdxS0SJEROQ7PH+kwOBdqV6XICLi\nK54He4t6Db0uQUTEVzwP9not23hdgoiIr3ge7CmxnbwuQUTEVzwP9viGenuSiEhV8jzYIyMjvS5B\nRMRXPA/2ltHtvC5BRMRXPA/2Fp3ae12CiIivVPpZMZVq3MwdO3KMxk0be1aDiEhtUpZnxXge7F62\nLyJS29TkQ8BERCREKNhFRHxGwS4i4jMKdhERn1Gwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iIzyjY\nRUR8RsEuIuIzCnYREZ9RsIuI+Eylgt3MxpvZRjMrMrNBp62/zMxWmdn6kv+OrnypIiJSFvUqeXwi\nMA54GTj9+bsZwBjn3AEz6wvMB2Iq2ZaIiJRBpYLdObcFTj4f+Fvr1572cRPQ2MzqO+cKKtOeiIh8\nv5oYY78B+FqhLiJSM763x25mC4DoUjZNds598D3H9gX+AFx2pn0SEhJOLQcCAQKBwPeVJCISNoLB\nIMFgsFzHVMmr8cxsMTDJObf6tHUxwELgR8655Wc4Tq/GExEph5p+Nd6phswsApgLPH6mUBcRkepR\n2emO48wsBRgKzDWzeSWbHgC6AVPMbE3JX2QlaxURkTKokqGYCjeuoRgRkXKp6aEYEREJAQp2ERGf\nUbCLiPiMgl1ExGcU7CIiPqNgFxHxGQW7iIjPKNhFRHxGwS4i4jMKdhERn1Gwi4j4jIJdRMRnFOwi\nIj6jYBcR8RkFu4iIzyjYRUR8RsEuIuIzCnYREZ9RsIuI+IyCXUTEZxTsIiI+o2AXEfEZBbuIiM9U\nONjNbLyZbTSzIjMbeNr6C8xsTcnfejO7qWpKFRGRsjDnXMUONOsFFAMvA5Occ6tL1jcGjjvnis0s\nGtgARDnniko5h6to+yIi4cjMcM7Z2fapV9GTO+e2/LuRb63PO+1jYyC7tFAXEZHqUS1j7CXDMRuB\njcDD1dGGiIiU7qw9djNbAESXsmmyc+6DMx3nnFsJ9C0ZrvnYzILOuezS9k1ISDi1HAgECAQCZShb\nRCQ8BINBgsFguY6p8Bj7qROYLea0MfZSti8EHnPOfV3KNo2xi4iUQ1nG2KtqKOZUI2YWZ2b1SpZj\ngXhgWxW1IyIi36My0x3HmVkKMBSYa2bzSjaNANaa2RpgBnCPcy6n8qWKiEhZVHooplKNayhGRKRc\nanIoRkREQoSCXUTEZxTsIiI+o2AXEfEZBbuIiM8o2EVEfEbBLiLiMwp2ERGfUbCLiPiMgl1ExGcU\n7CIiPqNgFxHxGQW7iIjPKNhFRHxGwS4i4jMKdhERn1Gwi4j4jIJdRMRnFOwiIj6jYBcR8RkFu4iI\nzyjYRUR8psLBbmbjzWyjmRWZ2cBStnc2syNmNqlyJYqISHlUpseeCIwDPjvD9meAuZU4f1gJBoNe\nlxASdB2+oWvxDV2L8qlwsDvntjjnkkrbZmbXATuATRU9f7jRP9yTdB2+oWvxDV2L8qnyMXYzawY8\nBiRU9blFROT71TvbRjNbAESXsmmyc+6DMxyWADzrnDtmZlbJ+kREpJzMOVe5E5gtBiY551aXfP4M\n6FSyOQIoBp5wzr1UyrGVa1xEJAw5587aaT5rj70cTjXinBt5aqXZFCC3tFAvS3EiIlJ+lZnuOM7M\nUoChwFwzm1d1ZYmISEVVeihGRERCiyd3nprZlWa2xcy2mdnjXtQQCszsH2aWZmaJXtfiNTPrZGaL\nS25622BmD3pdk1fMrJGZrTCztWa2ycye8romr5lZXTNbY2ZnmrQRFsxsl5mtL7kWK8+4X0332M2s\nLrAVuBTYC3wFTHDOba7RQkKAmY0AjgDTnHP9va7HS2YWDUQ759aWTJn9GrguHP9dAJhZk5KZZfWA\nZcAjzrllXtflFTN7GBgENHfOjfW6Hq+Y2U5gkHPu0Nn286LHfgGQ7Jzb5ZwrAN4GrvWgDs8555YC\nh72uIxQ45w4459aWLB8BNgMdvK3KO865YyWLDYC6wFn/R/YzM4sBrgZe5bSJGmHse6+BF8HeEUg5\n7XNqyToRAMwsDjgPWOFtJd4xszpmthZIAxY758L5Lu5ngUc5OXU63DngUzNbZWZ3n2knL4Jdv9bK\nGZUMw8wEHirpuYcl51yxc24AEAOMNLOAxyV5wszGAOnOuTWotw5wkXPuPOAq4P6S4dzv8CLY9/LN\nDUyULKd6UIeEGDOrD8wC/umce8/rekKBcy6bkw/TG+x1LR65EBhbMrY8HbjYzKZ5XJNnnHP7S/6b\nAczm5ND2d3gR7KuAeDOLM7MGwE3AHA/qkBBS8viJvwObnHN/8boeL5lZpJlFlCw3Bi4D1nhblTec\nc5Odc52cc12Am4FFzrk7vK7LC2bWxMyalyw3BS7n5FN2v6PGg905Vwg8AMzn5NMf3wnjmQ/TgS+A\nHmaWYmZ3el2Thy4CbgNGl0zlWmNmV3pdlEfaA4tKxthXAB845xZ6XFOoCOeh3Chg6Wn/Lj50zn1S\n2o66QUlExGf0ajwREZ9RsIuI+IyCXUTEZxTsIiI+o2AXEfEZBbuIiM8o2EVEfEbBLiLiM/8HQSAn\nwMLcfG4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11012bb38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = apx.proj(lambda x: 0.5 * ramsey.f(x))\n",
    "\n",
    "fhat = df\n",
    "for _ in range(10):\n",
    "    vf = value_of_policy(fhat, ramsey, apx, T=100)\n",
    "    fhat = greedy(vf, ramsey, apx)\n",
    "\n",
    "for _ in range(3):\n",
    "    plt.plot(apx.centers, vf)\n",
    "    vf = value_of_policy(fhat, ramsey, apx, T=100)\n",
    "    fhat = greedy(vf, ramsey, apx)\n",
    "    \n",
    "v, h = analytic_solutions(A, α, ρ)\n",
    "plt.plot(apx.grids, v(apx.grids), ':')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x110131160>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEACAYAAACuzv3DAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXJyEJYQkYhixBQcRtC4gDvMgwDgRFBdxW\nqbWgqNW6MfanbW1rxT1BnMUKyhKJMoIiqICAjICEPQRkb8j4/P5IqmkaMsg4ubnv5+ORB/fc873n\nfLgPfOfr95zv95i7IyIi4Ssq6AJERKRkFOQiImFOQS4iEuYU5CIiYU5BLiIS5hTkIiJhrtAgN7NE\nM1tiZsvM7P589ieY2UQzm2dmC83spjKpVERE8mUF3UduZtHAUqArsB6YBfRz99RcbZKAOHd/0MwS\ncto3cPeMsixcRESyFdYjbw+kufsqd08HRgA987T5EaiV87oWsFUhLiJSfqoUsr8xsDbX9jrgrDxt\nXgemmNkGoCZwdemVJyIihSmsR16U+fsPAfPcvRFwBvCimdUscWUiIlIkhfXI1wNNc203JbtXnts5\nwJMA7r7czFYCrYHZuRuZmRZ1ERE5Au5uBe0vrEc+G2hlZs3NLBboA4zN02YJ2RdDMbMGZIf4isMU\nox93HnvsscBrqCg/+i70Xei7KPinKArskbt7hpkNBJKBaGCou6ea2W05+18F/gy8aWbzyf7F8Ed3\n31aks4uISIkVNrSCu38KfJrnvVdzvd4C9Cj90kREpCg0szMAoVAo6BIqDH0Xv9B38Qt9F8VT4ISg\nUj2RmZfXuUREKgszw0t4sVNERCo4BbmISJhTkIuIhDkFuYhImFOQi4iEOQW5iEgFdSjzUJHaFToh\nSEREyt7s2TDyowxSd3zHhrgp/Bg3lc1xM4v0WQW5iEhAsjyLOevn88gbU0lZPQU7djoJjZtxYlxn\n2sfdzglxI7iLuoUeRxOCRETKibuz6KdFTF05lamrpjJ1xTQObKtPvT2defjazlxxZoh61ev912eK\nMiFIQS4iUso2bIBHHoH53zsHqi1jd8IUdh89lT0JKURl1KT6T52pvrkzuxeEeOL+RgwYAHaYqFaQ\ni4iUszdHreOu5ybT7PxJbKo2hSiLpl3CBbRN6Ey7hM40qt4Ms+zgPvZYOOqogo+nIBcRKSXukJX1\ny5//+dm2bwdfrk1h8spJjJo7iV0ZWwg160Kf9l3o0qILx9U5Djtcd7sIFOQiIqVgyhS4/Xb44QeI\nij2ANZsJLSaR1WISnrCY6PXnUmVNV0LNuvDeP0/n6Lqld2e3glxEpAQyM+FPT2Ty0qh59Bg0mXWx\nk5i5bian1D+Fri260uW4Lpzd5GziqsSVWQ0KchGRYjhwAObNc1bvXs6MjZN4b+ZkdtadQot6DUg8\noStdWnQh1DzEUVULGdguRQpyEZEi2Lx3M89/MplnRk/iUNPJEJVOnW1d6Xp8F/762y40rd04sNoU\n5CIi+TiUeYgZa2eQnJZM8vLPSN24HF95Prd07sbAi7twYsKJJbpAWZpKJcjNLBEYQvbDl99w96fy\n7L8XuDZnswrQBkhw9x152inIRaTMrV4Nn34KBw9mb2ffZeL8lLWMJRnJLM34jJWZX1AvqjUnRF/I\nxundqXeoA++9HUOjRsHWnp8SB7mZRQNLga7AemAW0M/dUw/T/lLgLnfvms8+BbmIlKl33oF77oFL\nL4XYWjtYGzOFtTHJrIn5jCxL59j0Czk2szvNMrpSjaMxgzZt4OabIaqCLiFYlCAvbK2V9kCau6/K\nOeAIoCeQb5AD1wD/KmadIiIlcuAA3DEok+QFs7jqxWTm7/mMBZsWcG6zc/nNcd25sOUg2iS0qTDD\nJaWtsCBvDKzNtb0OOCu/hmZWDbgQ+H3plCYicnhjxsB749ewNjaZxYc+48AxkznhtKZUr92dx9s+\nznnNzqNqlapBl1kuCgvy4oyF9ACm5x0bFxEpLQcyDjA57QsefvNTFqdPJK7pVk6t1o1b6l3KvZc/\nR6NaxwRdYiAKC/L1QNNc203J7pXnpy+FDKskJSX9/DoUChEKhQotUEQi26odq/h02adMSJvAtJVf\nwOZTabz/Ij7/w7t0bHUmUVZBB7ePUEpKCikpKcX6TGEXO6uQfbGzC7AB+JZ8Lnaa2VHACqCJu+8/\nzLF0sVNEDuvQIXjoIZg46SB7j57O7gYT2N3wUzLjtlD9x4uovuFi9i3oxiN/qMs991Tci5OlrcQX\nO909w8wGAslk33441N1Tzey2nP2v5jTtBSQfLsRFRAry7ZK19Bv8KXsbT2DPVVM5rmYbLm54MZ2O\neZuT6vyK6KgozKBOHWjYMOhqKx5NCBKRcpWRAfsOpPPV2q+YuHwCYxZ/ypptP3J69UTuuewiLmp1\nIQnVEoIus8IojdsPRURKxZrtP/Kbv0wgZd0EMo+dTNT2VkSvvJimB94geXBbunWJDrrEsKUeuYiU\nCXdn7sa5jFs6jtGp41m4Po16u7vzyNWX0Pv0C2lQo0HQJYYFrbUiIuXixx+z1+o+kLGfOdsnM3Pr\neL7eNp7YqHjOrNaDOf/qwfXnn8dfnoghWh3vYtHQioiUuQ8m/Ej/v40n5pRx7KyTQo09Z3L01h4c\nv3Uy1Q+0JisOXrkfevYMutLKSz1yESmW/wyZjF06jmHTx7Fu7wo6N72QWzr2ILFlInXj6wZdYqWi\nHrmIHJGDB2H4cNi1K3v1wENZ+1mWOZklWeNZkjWeGOJJ2NaD6AV/J/W182jdMibokiOaeuQi8l9+\n+gmuuAKia26h6mnjWR4zmjXRU6ifdSYtM3vQMqsHCbSmRg24806oWTPoiis3XewUkWKZ+PUqrv3T\nGGq2G832+O/oelxXerXuxSUnXKIhk4AoyEXksLKy4JlnnKmLF7Cq6sesrzmanayjY/0e3HvJ5XQ9\nrivxMfFBlxnxFOQi8j8yszKZuvwrBr4wmlVVR1OzltO+5uW0q9mLy9uew+mn6tJZRaKLnSICwP70\n/UxaMYnRS0Yzduk49m9uTLP0Xky752PaNzut0j5wIVKoRy5SyezeDQMGwLzUXWypO57dTT5iX8PP\nid12JtXWXE76wp4MuKY5Tz4ZOSsIhjMNrYhEmCWrt5N45zi8zUi21kyhbb1OdG/amy5NenB0fAJR\nUdl3mdSrF3SlUlQKcpFKbO9eSE+HLfu28MnyMXywYBRfr5/OCbEX8HCvK+nR+lJqV60ddJlSQhoj\nF6mE0tOh/92beG/OaGgzkoyG3xKzpjt1f7yRl677gN/drBu7I4165CJhYsPuDYyY/xFPfjyKXdXm\n0vOki+h3+pUktkykemz1oMuTMqKhFZEwdPAgzJmTPXSy+cA6vtgykunbRrJ6/yLi1/TglKgrGftM\nd2rGR8YT4iOdglwkzGzeDF16bmL7MSPZ2/wD9lRbSP3tPWm49SoSdnUhsVscgwbpbpNIoiAXCRPb\n9m9j2IyPeWzkCDIbzKL3qZfQ95S+dD++O3FV4oIuTwJUKkFuZonAELIfvvyGuz+VT5sQ8AwQA2xx\n91A+bRTkIkBmJrzxBqzeuIvFmWNZyAjWRn1JldXduOKEvrx678VUi6kWdJlSQZQ4yM0sGlgKdAXW\nA7OAfu6emqtNbeAr4EJ3X2dmCe6+JZ9jKcgl4u3ct4+L7/qEJVVGsLf+JI61TpwR3ZdT4y7jzJNq\ncsklQVcoFU1p3H7YHkhz91U5BxwB9ARSc7W5Bhjl7usA8gtxkUh2KPMQyWnJvPf9CEZ9/wm1a7Tn\nyT596XPaG9SJrxN0eVIJFBbkjYG1ubbXAWfladMKiDGzqUBN4Fl3f6f0ShQJL7t2wYMPOYt2zWD1\nUe+yofaH1DxwEtGp/bjkqGf4cHh9YvQcBilFhQV5UcZCYoBfAV2AasBMM/va3ZflbZiUlPTz61Ao\nRCgUKnKhIuHg25WpXPHYe+xs9h61m1QjVOd6zqs9h/qxx9Lgcjj7bND6VFKQlJQUUlJSivWZwsbI\nOwBJ7p6Ys/0gkJX7gqeZ3Q/Eu3tSzvYbwER3H5nnWBojl0rpx90/MmLhCN6a9y6L12zk5Kx+DLvr\nOs5oeLpWFZQSK40x8tlAKzNrDmwA+gD98rQZA7yQc2E0juyhl38eScEiFd3evXDddTBnwW72Hvsx\ne49/l0MJs6i6qhc+/2/c3j3EkH9Gq9ct5arAIHf3DDMbCCSTffvhUHdPNbPbcva/6u5LzGwi8D2Q\nBbzu7ovLunCR8rZjVwbn3fgZ+1q+y452E2jfoBNXtLyVbs3GEB8TT/XqUEfXLiUAmhAkkg932LEj\nuwc+d/1C/v3DcD5IfY86tGBwr+vpc8pVJFRLCLpMiQBa/VDkCGRmQt+btzJ21b/w04bj1TdSd+0N\nXN8khdf/2lrT46XCUY9cJEdGVgaf/DCRO4cOZ0PVSVxx6sXc8uub6NKiC9FR0UGXJxFKPXKRw1i1\nChYtyn6S/Kp9C5m8dTjTtr9Lld3HcfSGm1j9/Bs0qquHMkh4UJBLxHnnHRj0wDaO6f4+6+sN52DM\nRppuu4FfbZtGqzqtefo1qK7lvSWMaGhFIoa788R7Kfx54uvEnDSBS1pfxM1n3KyhE6nQtIytRLT9\n++H99+GHDRuZlT6c73woe3ZW5e5O/XnwkuuoG1836BJFCqUgl4h18FAmnX6TzIrab7AnYSonWW86\nxN7KvX3PokULzdaR8KEgl4izZucahn43jKenDCP6QEP+elV/rju9LzXj9EBiCU+6a0UqvalT4cWX\n01kZO4619V5nR41vqb3mGlptHseXH55OjRpBVyhS9tQjl7A1btp6+v7tdWLPfp1G8ceRWO+3nFvn\nSmrExdO5M1oqVioF9cil0nF3pqycwt+mvsTnaVO5+JJr+OuVyZxS/5SgSxMJjIJcwsKOAzt4Zupb\n/CPlZQ7tjyVm3u95/vrhDOivsW8RBblUaN/9+B0vzXqJkYtH4T9cTL/mb/CnQedSu7ZRTc8nFgEU\n5FIBuMPGjdmLVbnD/vQDjF/5b95Z8hI/7d9In5a3ccaXSzn1uPo896iesCOSly52SqAyMuD662Hi\nRIhL2MC+k15mb5vXiN3yK2ouGUD8+ouoEhXN+efDa69BFXU9JMLoYqdUePfdB8sPfMuFrz3LZys+\n5YZTr+GO9l/QOqF10KWJhA0FuZS7rCz48qt0np44is8znqV+p430aTyQV3q8SO2qWnFQpLgU5FKu\ntu7byqVJrzHbXiQhqiUv9vsjN3a4TItWiZSAglzKxcLNC3n262d5f95IqmzsxZRHx9Ox1RlBlyVS\nKRQa5GaWCAwh++HLb7j7U3n2h4AxwIqct0a5+xOlXKeEmT17YMgQZ9G+FL6J/jubo+Zy8v7fE//B\nUqYn1+fEVkFXKFJ5FBjkZhYNvAB0BdYDs8xsrLun5mk6zd0vK6MaJczsO5BBu5tHsfG4vxNTfQ+d\nqtzLr2M+Ir5OVTqPhRNPDLpCkcqlsB55eyDN3VcBmNkIoCeQN8h1Z6+w99Behs4dxiPjnyGmeSOG\n3fQoPU/sQZTpacUiZamwIG8MrM21vQ44K08bB84xs/lk99rvdffFpVeiVFSHDsEDD8C8ZZtZ1eAF\n1jV8mepbOnLMsvf49qOzOeqooCsUiQyFBXlRZvB8BzR1931mdhEwGjghv4ZJSUk/vw6FQoRCoaJV\nKRWOO/QbkMbs2H+w7awP6HR0Hx5t/BXH1jiBs8+G+PigKxQJTykpKaSkpBTrMwXO7DSzDkCSuyfm\nbD8IZOW94JnnMyuBX7v7tjzva2ZnJbFo8yJuHvZnvtuVzB/Ov50/nHcH9avXD7oskUqpNGZ2zgZa\nmVlzYAPQB+iX5yQNgM3u7mbWnuxfDtvyHkjC2/79cMXAOUxJf5L0hjOos/Quvn7iZdqeWivo0kQi\nXoFB7u4ZZjYQSCb79sOh7p5qZrfl7H8VuBK43cwygH1A3zKuWcrZ9DXT6ffyE2xtuIhHQ/fxmzPe\npUHdakRrDo9IhaBFs+R/bNgA+/Y50zdM4qWFT7Bi6zqiZzzAwvduoF7duKDLE4koWjRLiu2tt5zb\nn/mErI5/IqvKHuoueog2B/sy7I0q1KsbdHUikh/1yAXIfoTa0GnJ3P7vwRx3wkH+3H0wl7e5XPeA\niwRMPXIp0N69MHWqM2PjZD7cMph1P+3kxlZJvDaotwJcJIyoRx6hMjOhbe9prGg+GK+2kVO3JdH1\nmKtJGhytJ/CIVCBF6ZEryCPQV2u+4sa3BrN+z2pe7juY606/hipR+p8zkYpIQysCwPz5MHYs/LBr\nLl/GPchWW0rUl4+y4M3radkiJujyRKSE1COv5LZvhzbnLKf2FY+yIXYq3ao+wrlV+3NJYiyt9TQ1\nkQpPQysRbtOezXR65P9YfdT7PHzBXdx99t3UiK0RdFkiUgwaWolQTz+/mxe+e5q1jZ7n6O3Xs+Th\nJTSvVy/oskSkjOges0rkUOYhBr33An9c14oWv07jnfNmk/bCEIW4SCWnHnkl4O6MWTqGeybex4aF\nx/PcRRMZ0FvPwxSJFAryMHbvvTDsk/ns6HA3VN+Mff4i913enQG9g65MRMqTgjxMjZy4kVc2PEr8\nTWN5pmMSt5zRnypRVahaNejKRKS8KcjDyK5dsHz1AV5fMIRXF/yDyzrdxJs3LaV21dpBlyYiAVKQ\nh4mVK532N41iZ/v7qLrzDB5o/TVP/q5l0GWJSAWgIA8DS35ayrnP3kGVCzYw8YahXNDigqBLEpEK\nREFeQWVmwvjkvbz+wxNM2fE6ddY9TNpTA4mP05R6EflvCvIKyN256rGPGJ9+N00yO9Jz1/c88ngj\n4vVwHhHJh4K8gvlh6w/cMvJOZu5ay4gb3ubKtqGgSxKRCq7QmZ1mlmhmS8xsmZndX0C7dmaWYWZX\nlG6Jld8HH8DtAw/S7r7HOf25c1g6oTt/bjZPIS4iRVLgollmFg0sBboC64FZQD93T82n3efAPuBN\ndx+Vz7G0aFY+li+HX/X6ivir+1O/Skv61nqRtq2a0q0besCDiJTKolntgTR3X5VzwBFATyA1T7s7\ngJFAuyMrNTLt2L+TbkMexK8cwwu9n6V3m96Y0ltEiqmwoZXGwNpc2+ty3vuZmTUmO9xfznlL3e5C\nrFsHp/cZTb2kU9i5O4Mldy7kypOuVIiLyBEprEdelFAeAjzg7m7ZSXTYNEpKSvr5dSgUIhQKFeHw\nlcumPZv59V8HkHHi9zzX4V1u6XI+sbFBVyUiFUVKSgopKSnF+kxhY+QdgCR3T8zZfhDIcvencrVZ\nwS/hnUD2OHl/dx+b51gRP0Y+avEobv1oINELb2TV8CRqaGEUESlEaYyRzwZamVlzYAPQB+iXu4G7\nH5frhG8C4/KGeCTLzISO3bfydd078IZzqDX1Iz4ffjY1lOEiUkoKDHJ3zzCzgUAyEA0MdfdUM7st\nZ/+r5VBjWBv00jjmtP8dA8+7micumEuNuGpE6XEeIlKK9MzOMvLdop08+MUgJi/7klcuepNbu3UK\nuiQRCUNFGVpR37AMDH79K9oNPYO538bz0mnzFeIiUqY0Rb8UZWRlMHjSk/w17WX+2e017rrosqBL\nEpEIoCAvBT/9BG+NWc0rm69l17aq9Nz8HXdd1CjoskQkQmhopYQOHYKTr/6AR1a3I2FLT35X/TPe\nfVkhLiLlRz3yEtifvp/uz97J3rOmMf32CbRt1DbokkQkAinIj8D27TDo8TQ+qX4le1e14fNBc2jb\nqGbQZYlIhNLQyhHo/8+P+aDGOSTW78+nt75Px/YKcREJjnrkxZCemc5d4x/g432jGNV3PL3atQ+6\nJBERBXlRrFgB/fpvYkGbK8ncV4vrYufQq93RQZclIgJoZmeR9LhtDtMaXs61J9/MvW0fo/mxUURH\nB12ViESC0lg0K+INmfwvPqlzJ8MvfYUb2vUOuhwRkf+hID+MGV9n0v0vj7C3xQgG1pvEDe1OD7ok\nEZF8aWglH7sO7qLZPddwTLM9TBvwIfVr1Au6JBGJUFo0q5iysmDijHWc9I+OZO1ozNy7P1eIi0iF\np6GVXPrd9T2j4i6l2aaBfPnH+6gaq2doikjFpyDPMfzLSXwYfw3Drnqem9r2CbocEZEii/gg//57\n+Ntnw/n39vvpFz2Km9p2DLokEZFiieggX70azr33aeysF7i79jT+dOeJQZckIlJsEXvXirtz7qOD\nWRo1kvn3fk6TWk2CLklE5H+Uyl0rZpZoZkvMbJmZ3Z/P/p5mNt/M5prZHDO7oCRFl4dlaVmcfO8g\nvt3+CeOu+EIhLiJhrcAeuZlFA0uBrsB6YBbQz91Tc7Wp7u57c16fCnzs7i3zOVaF6JGnZ2ZwzO9u\nIa7hCj68fDzn/OqooEsSETms0pii3x5Ic/dVOQccAfQEfg7y/4R4jhrAliOqthykZ6aT+Po17I/a\nw+pHk6keWy3okkRESqywIG8MrM21vQ44K28jM+sF/AU4BuheatWVojffyuCer65lX/p+BrceTfXY\nuKBLEhEpFYUFeZHGQtx9NDDazDoC7wCt82uXlJT08+tQKEQoFCpSkSV1KCODAZOv46Qz9/Bm4kec\n3FohLiIVU0pKCikpKcX6TGFj5B2AJHdPzNl+EMhy96cK+MxyoL27b83zfiBj5AcOZtLlxRtYkLaV\nzc+NpmqVquVeg4jIkSqNu1ZmA63MrLmZxQJ9gLF5TnK8mVnO618B5A3xoIwek0W1frcwa9FPjL/+\nY4W4iFRKBQ6tuHuGmQ0EkoFoYKi7p5rZbTn7XwV6AzeYWTqwB+hbxjUXibtz5/j7aNVhGXPv/pxq\nMfFBlyQiUiYq5YSggwfh+lefYvSKd1id9AXH1K5bLucVESltEfmEoKwsOP2moaxo+grvJk5XiItI\npVfpgvzp8WNIO/YRFv5hGifWaxx0OSIiZa5SBfmfh8/h8bRbuSFmAifWOyHockREykWleULQZ1+v\n59FFvehhr/LUne2CLkdEpNxUioudew/tpcX/deIUu4opf3qgTM4hIhKEolzsDPsg/+jjLO6eeRWb\n19Qi7elhNG6sx7OJSOVR6R++7A63vvkUVvNHvnrwFYW4iESksL7Y+faXk9l54vN8f+csmhyl9VNE\nJDKFbZDf/dg6nt1/HTfWeI8mR+k2QxGJXGEZ5PsOHuL5TVdxZ9dBDLmywj+QSESkTIXdGHlWFlw7\n7GGqU49/9v5j0OWIiAQu7HrkXW6dzBcJ7/N24nyiLOx+D4mIlLqwCvJNu7Yxre5NfHjNm/Q+IyHo\nckREKoSwuY98zx7ngpevZuWCxvz09pBSrExEpOKqVKsfnt3/36Q1WswrF74TdCkiIhVKWAT5mi1b\nWdTkLib99iMuaKWn/IiI5Fbhh1YyMuDMP93I1vW12TD02TKoTESk4qoUQyu//VsyPxycxic3Lgy6\nFBGRCqlC37+3a98B3trye166+BW6dqoRdDkiIhVSkYLczBLNbImZLTOz+/PZf62ZzTez783sKzM7\nraSFHTgAHe4aQu2Dp3LL+YklPZyISKVV6NCKmUUDLwBdgfXALDMb6+6puZqtADq5+04zSwReAzqU\npLDx0zbwQ8I/mHrt1yU5jIhIpVeUHnl7IM3dV7l7OjAC6Jm7gbvPdPedOZvfAE1KWtgTMx+ibdSt\ndDy5ZUkPJSJSqRXlYmdjYG2u7XXAWQW0vwWYUJKinnpnFt/v/YwJFy0tyWFERCJCUYK8yPcMmlln\n4DfAufntT0pK+vl1KBQiFArle5y/zn6Q3578OIkX1CzqqUVEKoWUlBRSUlKK9ZlC7yM3sw5Akrsn\n5mw/CGS5+1N52p0GfAQkuntaPscp0n3kz4+fyt1T+rPzyVSqx8cU/W8iIlIJldaj3mYDrcysuZnF\nAn2AsXlO1IzsEL8uvxAvqqws54HkRxlwcpJCXESkiAodWnH3DDMbCCQD0cBQd081s9ty9r8KDAbq\nAC+bGUC6u7cvbjEvTEwmvcp2nr6pX3E/KiISsSrMFH13qHtPiEuP+S3v/PGacqlJRKSiK62hlXIx\nLPlb9sSs5LVBVwVdiohIWKkwQf7U9L/TpdrdxMdpbFxEpDgqxKJZ42csJy1jKiOvezPoUkREwk6F\n6JE/MuZFOla/ldNO1MJYIiLFFXiP/GDGQRZGvcNHnbWmiojIkQi8R37t/42mytbTuPSc44MuRUQk\nLAUa5BkZMHrt6/zjmv5EBf4rRUQkPAUan5/MWAEN5tP/vMuDLENEJKwFGuQvTn+LkzKvJa5KXJBl\niIiEtcCC3N35aue/ue50TccXESmJwIJ8Rtoi9mfs43c9ir0ki4iI5BJYkP/u+X9zZuxV1KpV4BIC\nIiJSiEDuIz90CBb7KD69fmgQpxcRqVQC6ZGP+3INUTU307WNhlVEREoqkCB/d2YyLa07Uaabx0VE\nSiqQJJ2xaSKJrRKDOLWISKVT7kG+d386m2tM5vZuF5b3qUVEKqVyD/L3U76j6sFjOaFx/fI+tYhI\npVSkIDezRDNbYmbLzOz+fPafaGYzzeyAmf2hoGONnj2TVlXPOdJ6RUQkj0KD3MyigReAROAkoJ+Z\ntcnTbCtwB/CPwo733eavuaDV2UdQqoiI5KcoPfL2QJq7r3L3dGAE0DN3A3f/yd1nA+kFHSg9HTbF\nzuSaTh2OuGAREflvRQnyxsDaXNvrct4rtk+nbyCq6h7aHdfqSD4uIiL5KEqQe2mdbNQ3X9OEDphp\nWr6ISGkpyhT99UDTXNtNye6VF1vy2GeoWT2apKQkQqEQoVDoSA4jIlJppaSkkJKSUqzPmHvBHW4z\nqwIsBboAG4BvgX7unppP2yRgt7s/nc8+r3fHZfQ/6waevLZ3sYoUEYlUZoa7FziMUWiP3N0zzGwg\nkAxEA0PdPdXMbsvZ/6qZNQRmAbWALDMbBJzk7ntyH2tH7CLObXXyEf51REQkP4X2yEvtRGbOw1XZ\n/dBualQLZNFFEZGwU5QeebnO7Kyy+3iFuIhIKSvXIK+Z2aI8TyciEhHKNcjrxRxbnqcTEYkI5Rrk\nTWs2L8/TiYhEhHIN8hPqNy/P04mIRIRyDfJTmjYrz9OJiESEcg3y4xsdXZ6nExGJCOUa5E0Sapfn\n6UREIkJaoA08AAAD6klEQVQ5B3mt8jydiEhEKNeZnVlZjhY+FBEpugo3s1MhLiJS+sr94csiIlK6\nFOQiImFOQS4iEuYU5CIiYU5BLiIS5hTkIiJhTkEuIhLmFOQiImGu0CA3s0QzW2Jmy8zs/sO0eS5n\n/3wzO7P0yxQRkcMpMMjNLBp4AUgETgL6mVmbPG0uBlq6eyvgt8DLZVRrpZGSkhJ0CRWGvotf6Lv4\nhb6L4imsR94eSHP3Ve6eDowAeuZpcxnwFoC7fwPUNrMGpV5pJaJ/pL/Qd/ELfRe/0HdRPIUFeWNg\nba7tdTnvFdamSclLExGRoigsyIu6NGLe5bDKZ0lFEREpeBlbM+sAJLl7Ys72g0CWuz+Vq80rQIq7\nj8jZXgKc7+6b8hxL4S4icgQKW8a2SiGfnw20MrPmwAagD9AvT5uxwEBgRE7w78gb4kUpREREjkyB\nQe7uGWY2EEgGooGh7p5qZrfl7H/V3SeY2cVmlgbsBW4u86pFRORn5faEIBERKRtlPrOzKBOKIoWZ\nDTOzTWa2IOhagmZmTc1sqpktMrOFZnZn0DUFwcyqmtk3ZjbPzBab2V+CriloZhZtZnPNbFzQtQTJ\nzFaZ2fc538W3BbYtyx55zoSipUBXYD0wC+jn7qlldtIKzMw6AnuAt9391KDrCZKZNQQauvs8M6sB\nzAF6ReK/DTOr5u77zKwKMB24192nB11XUMzsHuDXQE13vyzoeoJiZiuBX7v7tsLalnWPvCgTiiKG\nu38JbA+6jorA3Te6+7yc13uAVKBRsFUFw9335byMJftaVKH/4VZWZtYEuBh4g/+9rTkSFek7KOsg\nL8qEIolwOXdFnQl8E2wlwTCzKDObB2wCprr74qBrCtAzwH1AVtCFVAAOTDKz2WbWv6CGZR3kupIq\nBcoZVhkJDMrpmUccd89y9zPInhHdycxCAZcUCDO7FNjs7nNRbxzgXHc/E7gIGJAzNJuvsg7y9UDT\nXNtNye6Vi2BmMcAo4F13Hx10PUFz953AJ0DboGsJyDnAZTljw/8CLjCztwOuKTDu/mPOnz8BH5M9\nVJ2vsg7ynycUmVks2ROKxpbxOSUMmJkBQ4HF7j4k6HqCYmYJZlY753U80A2YG2xVwXD3h9y9qbu3\nAPoCU9z9hqDrCoKZVTOzmjmvqwPdgcPe7VamQe7uGWTP+kwGFgMfROJdCf9hZv8CZgAnmNlaM4vk\nyVPnAtcBnXNur5prZolBFxWAY4ApOWPk3wDj3H1ywDVVFJE8NNsA+DLXv4vx7v7Z4RprQpCISJjT\no95ERMKcglxEJMwpyEVEwpyCXEQkzCnIRUTCnIJcRCTMKchFRMKcglxEJMz9P3UddRBfXju+AAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110131128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(apx.centers, fhat)\n",
    "plt.plot(apx.centers, h(apx.centers))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## なぜ収束するのか？\n",
    "\n",
    "$$\n",
    "f_{n+1}(x) = \\arg\\max_y u(x, y) + \\rho v_n(y)\n",
    "$$\n",
    "\n",
    "に対応する作用素\n",
    "\n",
    "$$\n",
    "  T_{n+1} w (x) = u(x, f_{n+1}(x)) + \\rho w(f_{n+1}(x))\n",
    "$$\n",
    "\n",
    "を定義. \n",
    "\n",
    "3つの性質に注意:\n",
    "\n",
    "$$\n",
    " T_{n+1} v_n = v_{n+1}\n",
    "$$\n",
    "\n",
    "$$\n",
    " T_{n+1} v_n = Tv_n\n",
    "$$\n",
    "\n",
    "$$\n",
    "  v \\le w \\Longrightarrow T_{n+1}v \\le T_{n+1}w\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 価値の単調性\n",
    "\n",
    "\\begin{align}\n",
    "v_n(x) &= u(x, f_n(x)) + \\rho v_n(f_n(x)) \\\\\n",
    "       &\\le u(x, f_{n+1}(x)) + \\rho v_n(f_{n+1}(x)) \\\\\n",
    "       &=: T_{n+1} v_n(x) \\\\\n",
    "       &\\le T_{n+1}^2 v_n(x) \\\\\n",
    "       &\\le \\cdots \\\\\n",
    "       &\\to v_{n+1}(x)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bellman作用素 $T$ との関係\n",
    "\n",
    "\\begin{align}\n",
    "    T^n v_0 \\le v_n \\le v^*\n",
    "\\end{align}\n",
    "\n",
    "#### 証明\n",
    "\n",
    "$n=0$\n",
    "のときは自明に成り立つ. \n",
    "$n$ まで成り立つとして, \n",
    "\n",
    "\\begin{align}\n",
    "    T^{n+1} v_0 = T(T^n v_0) \\le Tv_n = T_{n+1} v_n = v_{n+1}\n",
    "\\end{align}"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}