{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noisy Coin Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## License"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We turn on the ``division`` feature first of all, since this is just a plain good idea when working in Python 2.x, as it is a more sensible default for scientific computing and is the default in Python 3.x anyway."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get access to NumPy and matplotlib, IPython's ``%pylab`` magic command is quite useful. With the ``inline`` argument, all plots will be made a part of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.\n",
      "  warnings.warn(self.msg_depr % (key, alt_key))\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "try: plt.style.use('ggplot')\n",
    "except: pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to import from *Qinfer*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/metrics.py:51: UserWarning: Could not import scikit-learn. Some features may not work.\n",
      "  warnings.warn(\"Could not import scikit-learn. Some features may not work.\")\n",
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/IPython/parallel.py:13: ShimWarning: The `IPython.parallel` package has been deprecated. You should import from ipyparallel instead.\n",
      "  \"You should import from ipyparallel instead.\", ShimWarning)\n",
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/qinfer/parallel.py:53: UserWarning: Could not import IPython parallel. Parallelization support will be disabled.\n",
      "  \"Could not import IPython parallel. \"\n"
     ]
    }
   ],
   "source": [
    "# We need distributions to model priors.\n",
    "from qinfer import distributions\n",
    "# The noisy coin model has already been implmented, so let's import it here.\n",
    "from qinfer.test_models import NoisyCoinModel\n",
    "# Next, we need to import the sequential Monte Carlo updater class.\n",
    "from qinfer.smc import SMCUpdater\n",
    "# We'll be demonstrating approximate likelihood evaluation (ALE) as well.\n",
    "from qinfer import ale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "from scipy.special import betaln, gammaln"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we can compare to the analytic solution derived by Christopher Ferrie [cite]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def exactBME(k, K, a, b, gamma=1):\n",
    "    idx_k = np.arange(k+1)    \n",
    "    idx_K = np.arange(K-k+1)[np.newaxis].transpose()\n",
    "    \n",
    "    numerator   = (\n",
    "        gammaln(k+1) - gammaln(idx_k+1) - gammaln(k-idx_k+1) + gammaln(K-k+1) -\n",
    "        gammaln(idx_K+1) - gammaln(K-k-idx_K+1) + (idx_k+idx_K)*np.log(a-b) +\n",
    "        (k-idx_k)*np.log(b) + (K-k-idx_K)*np.log(1-a) +\n",
    "        betaln(idx_k+gamma+1,idx_K+gamma)\n",
    "    )\n",
    "    denominator = (\n",
    "        gammaln(k+1) - gammaln(idx_k+1) - gammaln(k-idx_k+1) + gammaln(K-k+1) -\n",
    "        gammaln(idx_K+1) - gammaln(K-k-idx_K+1) + (idx_k+idx_K)*np.log(a-b) +\n",
    "        (k-idx_k)*np.log(b) + (K-k-idx_K)*np.log(1-a) +\n",
    "        betaln(idx_k+gamma,idx_K+gamma)\n",
    "    )\n",
    "    bme = np.sum(np.exp(numerator))/np.sum(np.exp(denominator))\n",
    "        \n",
    "    var = np.sum(np.exp(\n",
    "            numerator - betaln(idx_k+gamma+1,idx_K+gamma) +\n",
    "            betaln(idx_k + gamma + 2, idx_K + gamma)\n",
    "        )) / np.sum(np.exp(denominator)) - bme ** 2\n",
    "    \n",
    "    return bme, var"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sequential Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's helpful to define a few constants, as we'll need to refer to them over and over below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "N_PARTICLES = 5000\n",
    "N_EXP = 250\n",
    "N_TRIALS = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a model to play with, using the prior $p \\sim \\mathrm{Uni}(0, 1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prior = distributions.UniformDistribution([0, 1])\n",
    "model = NoisyCoinModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to allocate an array to hold performance data. A record array is a rather convienent structure for doing so. First, let's define the fields in this array,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "performance_dtype = [\n",
    "    ('outcome', 'i1'),\n",
    "    ('est_mean', 'f8'), ('est_cov_mat', 'f8'),\n",
    "    ('true_err', 'f8'), ('resample_count', 'i8'),\n",
    "    ('elapsed_time', 'f8'),\n",
    "    ('like_count', 'i8'), ('sim_count', 'i8'),\n",
    "    ('bme', 'f8'),\n",
    "    ('var', 'f8'),\n",
    "    ('bme_err', 'f8')\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and then the array itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "performance = np.empty((N_TRIALS, N_EXP), dtype=performance_dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "true_params = np.empty((N_TRIALS, model.n_modelparams))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we run the experiments!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ALPHA = 0.1\n",
    "BETA = 0.8\n",
    "expparams = np.array([(ALPHA, BETA)], dtype=model.expparams_dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/ipykernel/__main__.py:8: RuntimeWarning: invalid value encountered in log\n",
      "/home/cgranade/anaconda/envs/qinfer-binder/lib/python3.5/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in log\n"
     ]
    }
   ],
   "source": [
    "for idx_trial in range(N_TRIALS):\n",
    "    \n",
    "    # First, make new updaters using the constructors\n",
    "    # defined above.\n",
    "    updater = SMCUpdater(model, N_PARTICLES, prior)\n",
    "    \n",
    "    # Sample true set of modelparams.\n",
    "    truemp = prior.sample()\n",
    "    true_params[idx_trial, :] = truemp\n",
    "\n",
    "    # Now loop over experiments, updating each of the\n",
    "    # updaters with the same data, so that we can compare\n",
    "    # their estimation performance.\n",
    "    for idx_exp in range(N_EXP):\n",
    "        \n",
    "        # Make a short hand for indexing the current simulation\n",
    "        # and experiment.\n",
    "        idxs = np.s_[idx_trial, idx_exp]\n",
    "        \n",
    "        # Start by simulating and recording the data.\n",
    "        outcome = model.simulate_experiment(truemp, expparams)\n",
    "        performance['outcome'][idxs] = outcome\n",
    "    \n",
    "        # Reset the like_count and sim_count\n",
    "        # properties so that we can count how many were used\n",
    "        # by this update. Note that this is a hack;\n",
    "        # an appropriate method should be added to\n",
    "        # Simulatable.\n",
    "        model._sim_count = 0\n",
    "        model._call_count = 0\n",
    "            \n",
    "        # Time the actual update.\n",
    "        tic = toc = None\n",
    "        tic = time.time()\n",
    "        updater.update(outcome, expparams)\n",
    "        performance[idxs]['elapsed_time'] = time.time() - tic\n",
    "        \n",
    "        # Record the performance of this updater.\n",
    "        est_mean = updater.est_mean()\n",
    "        performance[idxs]['est_mean'] = est_mean\n",
    "        performance[idxs]['true_err'] = np.abs(est_mean - truemp) ** 2\n",
    "        performance[idxs]['est_cov_mat'] = updater.est_covariance_mtx()\n",
    "        performance[idxs]['resample_count'] = updater.resample_count\n",
    "        performance[idxs]['like_count'] = model.call_count\n",
    "        performance[idxs]['sim_count'] = model.sim_count\n",
    "        \n",
    "        # Finally, record the ideal stats.\n",
    "        performance[idxs]['bme'], performance[idxs]['var'] = exactBME(\n",
    "            idx_exp + 1 - np.sum(performance[idxs]['outcome']), idx_exp + 1,\n",
    "            ALPHA, BETA\n",
    "        )\n",
    "        performance[idxs]['bme_err'] = np.abs(performance[idxs]['bme'] - truemp) ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff45cb02588>]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEGCAYAAABlxeIAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8FPW9//HXdxJIDISQhABJuBM0JXJREmxBooAFsac/\n6RFzjkd/R7S1FfV4t4q0Sltbb6BibbEqVaq2GqvSo61WK1osWgURLwHEICAXAwnhEm4hyXzPHwsZ\nwqWEZDczu/t+Ph4+yi47u5/9dMI7M9+Z79dYay0iIhL3HL8LEBGRYFAgiIgIoEAQEZH9FAgiIgIo\nEEREZD8FgoiIAAoEERHZT4EgIiIAJEbyzTdv3swLL7zA7t27uf766yP5USIi0koRPULo2rUrl19+\n+XFvV1ZWFoFqopN64VEvPOqFR73wtLYXzQqE2bNnc9lll3HjjTc2eX7p0qVce+21XHPNNcybN69V\nhRxM/wd71AuPeuFRLzzqhadNAmH06NFMmzatyXOu6zJnzhymTZvGzJkzWbhwIRs2bABgwYIFzJ07\nl61bt7aqOBERaTvNCoT8/Hw6dOjQ5Lny8nKys7PJysoiMTGRkSNHsmjRIgCKi4u5+OKLadeuHY8+\n+ihr1qwJ6xGEiIiEX4sHlaurq8nMzGx8nJGRQXl5eZPXdOzYkcsuu+yY71VWVtbkUKekpKSlZcUc\n9cKjXnjUC4964SkpKaG0tLTxcUFBAQUFBc3ePqJXGTXXkYreuHGjT9UES2pqKjU1NX6XEQjqhUe9\n8KgXnpycnFYFZIuvMsrIyKCqqqrxcXV1NRkZGS0uRERE/NXsQLDWcvBaOnl5eVRUVFBZWUl9fT0L\nFy6ksLAwIkWKiEjkNeuU0axZs1i2bBk1NTVMmTKFkpISRo8ezaWXXsodd9yBtZYxY8bQo0ePSNcr\nIiIRYoK0hOaBweWSkhKNIeyn86Me9cKjXnjUC09OTg6lpaXHPZh8QKAC4WAKhBDt7B71wqNeeNQL\nT05OTqu21+R2IiICKBBERGQ/BYKIiAABC4SysrImd9mJiMjxKS0tbfEkdxpUDjgNmHnUC4964VEv\nPBpUFhGRsFAgiIgIoEAQEZH9FAgiIgIoEEREZL9ABYIuOxURaR1ddhrDdEmdR73wqBce9cKjy05F\nRCQsFAgiIgIoEEREZD8FgoiIAAEOBLt7l98liIjElUAFwsGXndrnn/C3GBGRKNSay04Tw1xLqxy8\nDqgtX+5zNSIi0aekpKTF2wbqCKGJHVv9rkBEJK4ENxD27MHW1/ldhYhI3AhuIKR2gh3b/a5CRCRu\nBDcQOqXDdp02EhFpK8ENhLR0jSOIiLShwAaCSUvH6ghBRKTNBDYQ6KQjBBGRthSoQGiyHkJnjSGI\niByvmLwxzXRKx13+kc8ViYhEl9i8MS2tM+zY5ncVIiJxI7iBoMtORUTaVHADIS0UCAFd4VNEJOYE\nNhBMUjIkJGBf+gO2ssLvckREYl5gAwHAXHQFdmUZ9t03/S5FRCTmBeoqo0M5p52BW7sHVn/udyki\nIjEv0EcIACYjC1td6XcZIiIxL/CBQEYWKBBERCIuCgKhC1RX6WojEZEIC1QgNJm6Yj+TnAKJibCr\nxqeqRESiR0xOXdFERhZUV0HHTm1flIhIFInNqSsOlt5F4wgiIhEWFYFgMnWlkYhIpEVFIDSeMhIR\nkYiJjkDQKSMRkYiLikAwmV2xVZv8LkNEJKZFRSDQsw9s/BJbX+d3JSIiMSsqAsEkp0BWNqxb7Xcp\nIiIxKyoCAcD0Pwm7arnfZYiIxKyoCQT65cOqz/yuQkQkZkVNIJi8fOwXK/wuQ0QkZkVNIJCVDfX1\n2Ir1flciIhKToiYQjDGYM8/B/uU5v0sREYlJgQqEI812ejAz9t+wnyzGbtrYhlWJiESP1sx2amxA\nFxrYuPHI/+i7f3wCEhJwvvP/27Ygn6SmplJTo6m/Qb04mHrhUS88OTk5rdo+UEcIzWHyB2FXaXBZ\nRCTcoi4Q6HsSrCnHNjT4XYmISEyJukAwHTqGltXcsMbvUkREYkrUBQKA6XcS9vNluH/7E7a21u9y\nRERiQqCW0Gy2/vnYZx6BffswPfpC/mC/KxIRiXpRGQjm5GEwYizU1WHXr8YoEEREWi06TxmlZ+Jc\nOAX6nQRfagZUEZFwiMpAOMD07Iddr0AQEQmHqA4EcnvBpg1aOEdEJAyiOhBM+yTI7Aaa8E5EpNWi\nOhAATI8+WI0jiIi0WtQHAj37gcYRRERaLeoDwfTsg9VayyIirRb1gUCPvrB+NQGdtFVEJGpEfyCk\npYNxYOsWvysREYlqUR8Ixhjo2VfjCCIirRT1gQBgevbVOIKISCvFRCDQoy8oEEREWiVQgXCsNZWP\nxvTJw37xmQaWRSTutWZN5UDNdlpQUEBBQcHxb9gtF5KSYE059B0Q/sJERKJESUlJi7cN1BFCSxlj\nMKd8A7vkHb9LERGJWjERCABm2Ajsknd02khEpIViJhDo1R8aGmDjl35XIiISlWImEIwxmMGF2I8X\n+V2KiEhUiplAADCDirAfL/a7DBGRqBRTgUD+oNC8Rrtq/K5ERCTqxFQgmHbt4aRB2E8+8LsUEZGo\nE1OBAOCMGod9+Vnsvlq/SxERiSoxFwhmyHBMr37Y0jnYOq21LCLSXDEXCADmgh9gt2zGvfX72G3V\nfpcjIhIVYjMQUjuRcM106DsA+9knfpcjIhIVYjIQDjB9TwzNbyQiIscU24HQZwB2zed+lyEiEhVi\nOhDonQfrvsA2NPhdiYhI4MV0IJiUDpDeBb7S/EYiIscS04EA+08brdZpIxGRY4n5QGDA12CFrjQS\nETmWmA8EM6gIW7YEW1/vdykiIoEW+4GQngldusGq5X6XIiISaDEfCBCazsIufd/vMkREAi1OAqFI\nC+eIiBxDXAQCPfrC3t3Yygq/KxERCay4CATjOJiBQ7HLl/pdiohIYMVFIAAw8BRsWSgQ7L5aLaIj\nInKIuAkE87UhsOJjbEMD9tk5uA/+RDesiYgcJH4CoXMG9OiN++Mp2BUfYf7je7jPPoq11u/SREQC\nITGSb75o0SKWLFnCnj17GDNmDIMHD47kxx2Tc8PPYe0qSM+ETp2xC9+AT5dgBw6FT5fA4EKMMb7W\nKCLil4gGQlFREUVFRezatYsnn3zS90AwjgN9B3iPz5yA+4/XMDu2Yef+ElM0Ci65GpPYzscqRUT8\n0axAmD17NkuWLCEtLY0ZM2Y0Pr906VKeeOIJrLWMHj2aiRMnHnH7559/nrPPPjs8FYeRKRqFfX4u\ndu0qnP/5Ee5r8+CtVzBn/T+/SxMRaXPNGkMYPXo006ZNa/Kc67rMmTOHadOmMXPmTBYuXMiGDRsA\nWLBgAXPnzqW6upqnn36aU089lT59+oS9+NYyKR0wQ4dD5ww4eRjO+Zdg//oCtm6f36WJiLS5Zh0h\n5OfnU1lZ2eS58vJysrOzycrKAmDkyJEsWrSI3NxciouLKS4u5pVXXuHTTz9lz549VFRUcNZZZ4X/\nG7SSOf+7mLp9obGDXv2hdx7uIzMwQ4fjjAxevSIikdLiMYTq6moyMzMbH2dkZFBe3nT94gkTJjBh\nwoRjvldZWRllZWWNj0tKSkhNTW1pacfnkM9xL/8h+959k9qXnuGEXn1JHDi0beo4ivbt27ddLwJO\nvfCoFx71oqnS0tLGPxcUFFBQUNDsbSM6qNxcRyq6pqbGn2IS28Oo8dA+mZ1PPITpdxJ0zsD5Vokv\n5aSmpvrXi4BRLzzqhUe98KSmplJS0vJ/q1p8H0JGRgZVVVWNj6urq8nIyGhxIUFjikZBp85gXezr\nf8Ju2ex3SSIiEdXsQLDWNrmJKy8vj4qKCiorK6mvr2fhwoUUFhZGpEg/GMch4ZrbcS66AnPGBOxz\nj2MrNvhdlohIxBjbjFt1Z82axbJly6ipqSEtLY2SkhJGjx7Nhx9+2HjZ6ZgxY4562WlLbNy4MWzv\n1Vp29y7cJ2bB52U4k6/BDBneZp+tw2GPeuFRLzzqhScnJ6dV2zcrENrKgcHlkpKSQAXCAe7fX4WV\nn+JcdmObfaZ2do964VEvPOqFJycnh9LS0uMeTD4gEIPKB7T0S7QVc8ppuM/Pxdbtw7Rr73c5IiKH\n8WVQOR6ZTunQqx+Ufeh3KSIiYadAOE5m+CjcV58/7G5mu22LTxWJiISHAuE4mdPHYTKycB+7r/Gq\nK7vkHdybLsF95Y+4b7yMu/ANn6sUETl+CoTjZBwHc8m1sLUqdH9CxXrcp2bjXH4L9qP3sZ8uxj7/\nBLa+zu9SRUSOS6ACoaysrMlt10Fl2rXD+f5N2Ff+iHv3zZiJF2GGjSDhlntIuGY6dM+Fjxf7XaaI\nxKHS0tImUwEdj0BddnqwIF52eii7eSN0TMOkdGjyvPvOG9jFC0m4+rZWf4YuqfOoFx71wqNeeFp7\nH0KgjhCijemac1gYAJhhI2HN59iv1vlQlYhIyygQIsAkJWO+ORH7p9/7XYqISLMpECLEjPkWtnw5\ndmXLzuWJiLS1QN2pHEtMUjLO5P/BffgunCunYfrnN/l729CAfXo2dsNaTHoX6Hsi5rRiTOfMo7yj\niEhkBeoIIVquMmouc/IwnIuvxn34buyunY3PW9fFPv4AdmsVzvmXwNDhsGkD7u1X4b6/wMeKRSTa\n6SqjgHN//xvYsxvnu9eFHv/lOezS93Bu/DmmfVLj6+zqlbi/+gXOz36NOSEF0BUUB1MvPOqFR73w\n6CqjKGDOuxj7eRn2s0+wyz/Czn8Z5/JbmoQBgOl7ImbgUOwrz/lUqYjEMwVCGzBJyTiTJuM+/TDu\nnPtwvns9JqPLkV977oXYBa9ha/e2cZUiEu8UCG1l2EjIzMJ881zM14Yc9WUmMwvyvoZd9HYbFici\noquM2owxBufq2zHGHPO1zhln4/7p99ivnxn5wkRE9tMRQhtqThgAUHAKdOqMe/N3qft0SWSLEhHZ\nT4EQQMZJIOHq23AunMKe387CNjT4XZKIxIFABUKs3YfQaqd8HSc9E/vXFw5bkEdE5Eh0H0IMS9lW\nRc2v74JNG3C+dwNmUKHfJflG15t71AuPeuFp7X0IGlQOuISefUm4dQZ21QrcX/0ceudh+ufj/Nt/\n+F2aiMSYQJ0ykqMz/fNxbp2Bc+Y52PkvYzes9bskEYkxCoQoYrp0wwwpwoybiH35Wb/LEZEYo0CI\nQubMc0KnkP7yHO7f/hf3kXuxFesBsB8twn10Bnb3rlZ/jq2v02C2SBzRoHLAHW3AzFZX4j4+CxIT\nMXkDsa//KbSWc9UmzEmDsF+tw7l2OqZTerM+x366BKyFpGRwG6BPHu4D0yHphND7NPceigjS4KFH\nvfCoF57WDiorEAKuuTu7rd0Ln5dBdk/IyMK+/Cz23fmYIcOhvg5zwfdh9y7Ysxu6dIOabZDaGWMM\ndsk7uH94FLJ7wL7a0H+bNmKGjcRu/BJTPB6neHwbfNt/TT/4HvXCo154dJWRAKEJ9Dh5mPf42/+J\nm9UNKjdhV3wEv/8N9pMPoG4f1O4BkwA5PTE9+2KXvItzze2YPgMAsNbCpo2Q1R3z1TrcmT/C3bwR\nc+6FmHbtffqGIhJpgTpCKCsro6ysjJKSEh0h7BeO337sti24d92M+VYJ5vRvwp5dkJwCH76L3bYV\nkz8Yk9vr6Nvv2Ir7u1+B6+JMmYpp165V9bSUfhP0qBce9cKTk5NDaWkpBQUFFBQUHPf2gQqEgykQ\nQsK1s1trWzUOYOvrsY/NxG7agHPuhZihp7W6puOlH3yPeuFRLzw6ZSTN0tpBYZOYCD/4IebjxbhP\n/RrHdaFrNnbpP2HjOuy+WpwJkw5bO1pEoocCQZrNGANDinDSbsW97zZo3x7z9dEwuBBTW4v7q5+H\nVoI78fgPVUXEfwoEOW6mzwCcW++F9C6hwez9bFY33IfvwkyajDn1G5jkFP+KFJHjphvTpEVM9x5N\nwgDADDwF57qfYhf8FffGyaEb5LZU+lShiBwvBYKElenZl4Rb7sGZMRcysnB/+VNsfX2T19ja2sb/\ntXv3+FGmiByBThlJRJjkE+Df/xu77gvs7x/GrViPOfFkSErGvvgU9O4PVRWhO6Gn3otJa94d1SIS\nOTpCkIgxxuBcOAX75ReYEWOxqz/Hvr8A5xe/wTn3QpxbZ2JO/ybugz/B7tJlgyJ+030IARdL11hb\na8FajOM0ec7+8XHssqU41/3kX869FEu9aC31wqNeeFp7H4KOEKTNGGOahEHjc5MuwZw6AvfuqRqE\nFvFRoAJBayrHJ2MMzrf/EzN6Au7Pr8d9Z77fJYlELa2pHMPi7XDYfrkK95EZmFHfxBn/703+Lt56\n8a+oFx71wqNTRhJTTK/+ODfcgf37q7hv/hkAu30r7u8eonb+n5u81tbXYXfu8KNMkZikI4SAi9ff\nfmzVJtx7b4WefaF8OeYbYzCffoDboSPU1mK65WC/+AysizP1XkhIhJQOcTM9d7zuF0eiXni0QE6M\ni+ed3VZWYJd/hBlUiEnPpIOBnR8vho5p2K/WYXJ7YVd8gn19HuzaiRkyHPODHwZidbdIi+f94lDq\nhUeznUrMMlndMVndGx87HVMxgwpDf9c3tJgPvfpjsnuElvy8fzq8+jyM+w6s+RwSEkJ/7+jMqEhz\nKBAkqhljYHARAM6Um3Efn4X983OQlg6OA7t3htaY3leLM2IM5tQRPlcsElwKBIkZpmsOCTffHbqX\nIaNLaL3oygrsyjKMAfep2Tid0jF5X/O7VJFAUiBIzDGZWd6fDzrtZFM74/7yZ9A/H2fy1ZhOnf0q\nUSSQdHJV4oYZNAznnscxOT1xH/wpdu9uv0sSCRQFgsQVk5SEOW8ypu8A3Dt/iP3yC79LEgkMnTKS\nuGOMgf+6HN59E/eB2zH5gzGnnwUnnoxJbOd3eSK+USBIXDLGYEaMwZ76deyC13DnPQ0b1sLXhuBM\nvhq2VIZudMvqHpqRddHb8OUq6NUfZ3ix3+WLRIQCQeKaSU7BjJsI4yZi9+7GvlyK++MpYBxo1w7n\noitwX/kj1O7FFJ6OfeZRbG4fTG4vv0sXCbtABUJZWRllZWWUlJT4XYrEIZOcgpk0GVs4ErJ7YRe+\njvvb+zHnXogpHo9xEnATEnHnPUXClbf6Xa7IEZWWllJQUEBBQcFxb6upKwJOt+V7/OiFdd2mC/rU\n7cP98RU451+CGTayTWs5mPYLj3rh0WynIhF02II+7drjXHEr7tMPYz9636eqRCJDgSBynEyvfjhX\nTMX9wyO4T8zC1tf7XZJIWCgQRFrA5A3E+clD2B3bsY/NxDY0+F2SSKspEERayCQl40y5BburBvvC\nXL/LEWk1BYJIK5h27XEuvxn74T9x33oFALulUqeRJCoF6rJTkWhkOqTiXDMdd9Z0Gha9DWvLISML\nM2ocplsOZv/03CJBp0AQCQPTLQfnlruxS97FXPUj+Oxj7PKPcee/jKmswBn7bb9LFDkmBYJImJhO\n6Zgzzwk9GPp1zNCvY8dNxL1nKu7GLzHfPBe65mgFNwks7ZkiEWQyu+LcNguST8B9YDruzZdiV68M\nDUSvLMOuW+13iSKNdIQgEmGmQ0fM+ZfC+ZdiP3of94Hp4DZAbm/YuA5n+i8xGV38LlNEU1cEnW7L\n98RKL2zNdkhKxrRPwn3mUWjXHue8iw9/3dpy3HlP4Zx7IabPgNBzO7Zil/yTjoUj2N0xrY0rD6ZY\n2S/CobVTV+gIQaSNmVTvH3Iz5t9w77wRO+E8aJcEyz7EVn4FVZux7/0dUzwe98GfQkYW1NdBdSWc\nNIidf34Wc8EPoEef0PrRx7mOg62vg7WroF077FfrMd17YHr3D/M3lWijQBDxkemajRkxFvfWH0Bi\nO+iWg+nRBzKzcG74GaZHX+z478BX68FJgJ59MIntSF75Kbv+9w9QtQm2V0PvPMyZEzCnnQnly7Gf\nl0HHTpjCkbhz7gfHwZkwCdPvpND6Dr99ALu2PPSZ3XOxX6zEDBiIufTaxnCxdXW4M24NbTf0NF/7\nJG1Dp4wCTofDnljuhd1SCfv2YrJ7Nuv1B/fC1tfB8o9wn99/t/TePZjCkdgvPoMvvwiFRI/e2Fdf\ngIwuoRCo3Ytz488x7ZNC71FXh/ubuzHtkzDfuz401ffCN7DzX4Lt2zBnfRsz7juBvEIqlveL49Xa\nU0YKhIDTzu5RLzxH6oWt24dd+j5mcCEmKRnrulC+HAYMxBgTunt62YfY2r2YgUMxHVIP296d9ZPQ\nUctFV+D+9BqcSZMhuyfuY/dBYiLOJdcGbgBc+4VHgRDjtLN71AtPpHph9+7Gvf92+GodZHbDue2B\nUJi4DdhXnse+8RLOf/0AU3h62D+7pbRfeDSoLCJhY5JTcG6+C/bshvZJGGNCzzsJmG+VYAcOxX1s\nJqZ8Oabku4edQrK7d0HyCYE8tSTHFqgjhIOX0NQRQoh++/GoFx4/e2F378R96A5MWgbmv6+CxMTQ\nGtRry3Hv+3HoqqfeediaHbBnFyZ/MGbstzFp6RGpR/uFJycnR0toxjLt7B71wuN3L+y+Wuyzc7Af\nvgv79kG7RMDgXHINJJ+ArdqESe0M7ZOw772FXb0SZ+q9jYPY4eR3L4JEYwgxTju7R73wBKUXdsNa\nSM+EPXtg53ZM77zDX2Mt9rGZkJCAmXxN2E8nBaUXQaBAiHHa2T3qhSfaemH37sGdNR3S0qFyEyQl\nYUZ/C/vc49C7P86470D//BaFRbT1IpJaGwga+RGRiDPJJ+Bc+xNM50zM2edhCk7Flv4W56IpmPzB\nuE/9Gve2K7FfrvK71LimI4SA028/HvXCEwu9sNY2XsUE4L6/APuHRzDnX4IzYmyz3ycWehEuuuxU\nRKLSwWEA4Awvxub2wZ19J+4nH2DOPg9yemLatfepwvijQBCRwDC5vXB+dB/2zT/j/uZu2FoF3XJD\nd1+PGo/J6u53iTFNgSAigWKST8BMmAQTJmEbGuDLL7DvL8D9xQ3QoVNoXqXi8RgnoXEbu3YV9m//\ni5l4ESYzy8fq/zVbvgwsmAED/S7liDSGEHA6P+pRLzzx2AvrurB6Je4Lv4Md2zCnFWNXfUZC7V4a\nKtZjTvkG9qP3cK6+DdMr/FN52101kJyCSUg4/O9qtmNfmwcHphHP7Io5IQW7c0fopj1jsH98HFu2\nBBISoXsPnO9dj0npiLUW9tVikpJbXaMuO41x8fiDfzTqhSeee2GtheUfYZf+E3PSIFKyc9md0gnT\nOQP7wULcZx7DmXoPJuPwIwVbWQHbqrErPsb+4zWc/7ocM2T4UT/Lnf9yaDK/dkm4s+8MPZk/GHPa\nGZjBwzFJSdhPFuPOfQgzpAi7syY0D9S2ahgwMDS5oLWQmIg55RuYSZOhXfv94fBhKNRWfALly6BX\nf5xJkzEnntzi3igQYlw8/+AfSr3wqBeeQ3vhvjYPO/9lnIv/B/btw679HPbugcRE7Nuvh9acyOyK\nKTod96nZmMLTMaPGQ3oGrF+LrVgH9fWwayf23fnQ0AB7d+Nc9WPI7YVd+j72vb/DFysgvUtoKvFL\nr8WcNKixBluzHbv0PUzBqZCUBFu3hNa5OIj98J+hKcq752IKR8Gni3H/8Aj07Aftk2D9GpzzLsac\n+o1m90KBEOP0g+9RLzzqheeIU4F/sDC0PkRWNqZPHiQlw9aq0D0QmV29123bgn1tHnbpe6GFhnJ6\nY3J7g+NgqytxLpwCySmwZRMHljFt3Hb3LqiqCH3GCSlh+S62ZjusWo7duxfTIRX38QdCp5YGntKs\n7RUIMU4/+B71wqNeeGK5F3ZlGe7Dd+H88C5M99xjvl53KouIxChzYgFm4oW499xCw69+ga2ujOjn\nKRBERALMKT4b5+a7MV264j75q9Bkga4bmc+KyLuKiEjYmG45mPMmw7atuLdfhXvFeTT8eAr2kw+A\n/VdehYFuTBMRiQImMRHniqlQWQH982HVCtw590F2T9i5g4Tpv2z1ZygQRESihMnqDgem7xg4FOeW\ne2DjOuh/UljeX4EgIhKlmgREGGgMQUREAAWCiIjsp0AQERFAgSAiIvspEEREBFAgiIjIfgoEEREB\nFAgiIrKfAkFERIAI36m8YcMG/vKXv1BTU8PJJ5/MuHHjIvlxIiLSChE9QsjNzeWyyy7juuuuY+XK\nlZH8KBERaaVmHSHMnj2bJUuWkJaWxowZMxqfX7p0KU888QTWWkaPHs3EiRMP23bx4sW8/vrrFBcX\nh69qEREJu2YdIYwePZpp06Y1ec51XebMmcO0adOYOXMmCxcuZMOGDQAsWLCAuXPnsnXrVgoLC5k6\ndSpvv/12+KsXEZGwadYRQn5+PpWVTZduKy8vJzs7m6ysLABGjhzJokWLyM3Npbi4mOLiYpYtW8a8\nefOoq6vj1FNPDX/1IiISNi0eVK6uriYzM7PxcUZGBuXl5U1eM3DgQAYOHHjM9yorK6OsrKzxcUlJ\nSasXi44lqampfpcQGOqFR73wqBee0tLSxj8XFBRQUFDQ7G0DcdlpQUEBJSUljf8d/IXinXrhUS88\n6oVHvfCUlpY2+bf0eMIAWhEIGRkZVFVVNT6urq4mIyOjpW8nIiI+a3YgWGubLOScl5dHRUUFlZWV\n1NfXs3DhQgoLCyNSpIiIRF7C9OnTpx/rRbNmzaK0tJQtW7bwt7/9jQ4dOtCvXz+ys7N58MEH+etf\n/8oZZ5zB8OHDw1ZY165dw/Ze0U698KgXHvXCo154WtMLYw/+tV9EROJWIAaVRUTEfwoEEREBFAgi\nIrJfRGc7PV7NmRspll155ZWkpKRgjCEhIYE777yTnTt38sADD1BZWUnXrl257rrrSElJ8bvUsDvS\nfFn/6ru/+OKLvPnmmyQkJDB58mSGDBniZ/lhdaRePPfcc7zxxhukpaUBcMEFFzB06FAgtnuxZcsW\nHnroIbZv344xhrFjx3LOOefE5b5xaC/OOussJkyYEN59wwZEQ0ODveqqq+zmzZttXV2dvfHGG+36\n9ev9LqtNXXnllbampqbJc08++aSdN2+etdbaF1980T711FN+lBZxy5cvt6tXr7Y33HBD43NH++7r\n1q2zN91c6JBIAAADCUlEQVR0k62vr7ebNm2yV111lXVd15e6I+FIvSgtLbUvvfTSYa+N9V5s3brV\nrl692lpr7Z49e+zVV19t169fH5f7xtF6Ec59IzCnjA6eGykxMbFxbqR4Yg+51wNCs8WeccYZAJx5\n5pkx25P8/Hw6dOjQ5LmjfffFixczYsQIEhIS6Nq1K9nZ2YdNmxLNjtQL4LB9A2K/F507d6ZPnz4A\nJCcnk5uby5YtW+Jy3zhSL6qrq4Hw7RuBCYQjzY104MvGC2MMd9xxB1OnTuWNN94AYPv27XTu3BkI\n7RDbt2/3s8Q2dbTvXl1dTZcuXRpfFy/7yquvvspNN93Eww8/zO7du4H46sXmzZtZu3YtJ554Ytzv\nGwd6MWDAACB8+0agxhDi3c9+9jPS09PZsWMHd9xxxxEn+DPG+FBZMMTzdx8/fjyTJk3CGMMzzzzD\n7373Oy6//HK/y2oze/fu5b777mPy5MkkJycf9vfxtG8c2otw7huBOULQ3EiQnp4OQKdOnSgqKqK8\nvJzOnTuzbds2ALZt29Y4cBQPjvbdD91XtmzZEvP7SqdOnRr/0Rs7dmzjoX889KKhoYGZM2dSXFxM\nUVEREL/7xpF6Ec59IzCBEO9zI9XW1rJ3714g9BvAxx9/TK9evRg2bBhvvfUWAG+99VZM9+TQMZSj\nfffCwkLeeecd6uvr2bx5MxUVFeTl5flRcsQc2osD//gBvPfee/Ts2ROIj17Mnj2bHj16cM455zQ+\nF6/7xpF6Ec59I1BTVyxdupTHH38cay1jxoyJq8tON2/ezL333osxhoaGBkaNGsXEiRPZuXMn999/\nP1VVVWRlZXHdddcdccAx2s2aNYtly5ZRU1NDWloaJSUlFBUVHfW7v/jii8yfP5/ExMSYurQQjtyL\nsrIy1qxZgzGGrKwsvv/97zeeQ4/lXqxYsYLbb7+dXr16YYzBGMMFF1xAXl5e3O0bR+vFP/7xj7Dt\nG4EKBBER8U9gThmJiIi/FAgiIgIoEEREZD8FgoiIAAoEERHZT4EgIiKAAkFERPb7PyS/zyoxXlft\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff45cad7438>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(np.mean(performance['true_err'], axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}