{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multimodality and Resampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook tests the robustness of modified Liu-West (MLW) resampling to multimodality in the posterior distribution. We use as a test the familiar model\n",
    "$$\n",
    "    \\Pr(0 | \\omega, t) = \\cos^2(\\omega t),\n",
    "$$\n",
    "exploiting that the likelihood function is *even* in $\\omega$ to produce posteriors that always have exactly two modes. Because this structure is not explicitly added to the resampling, but is easy to reason about analytically, it serves as a nice test case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preamble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As usual, we start by configuring division and enabling plotting support."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "try: plt.style.use('ggplot')\n",
    "except: pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to view the incredibly verbose output from resampler debugging, it is helpful to point Python's logging functionality at a file instead of printing it inside the notebook. We use ``tempfile`` to pick a directory in a cross-platform way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Logging to C:\\Users\\Chris\\AppData\\Local\\Temp\\multimodal_testing.log.\n"
     ]
    }
   ],
   "source": [
    "import os, tempfile\n",
    "logfile = os.path.join(tempfile.gettempdir(), 'multimodal_testing.log')\n",
    "print(\"Logging to {}.\".format(logfile))\n",
    "\n",
    "import logging\n",
    "logging.basicConfig(level=logging.DEBUG, filename=logfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we import all the functionality from **QInfer** that we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Anaconda3\\lib\\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",
      "c:\\users\\chris\\dropbox\\software-projects\\python-qinfer\\src\\qinfer\\parallel.py:52: UserWarning: Could not import IPython parallel. Parallelization support will be disabled.\n",
      "  \"Could not import IPython parallel. \"\n"
     ]
    }
   ],
   "source": [
    "import qinfer as qi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model and Prior Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want the simple precession model $\\cos^2(\\omega t)$ to now extend over $\\omega \\in [-1, 1]$, so we redefine the ``are_models_valid`` method to allow for negative $\\omega$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class NotQuiteSoSimplePrecessionModel(qi.SimplePrecessionModel):\n",
    "    @staticmethod\n",
    "    def are_models_valid(modelparams):\n",
    "        return (np.abs(modelparams) <= 1)[:, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having done so, we now set a prior that explicitly includes the degeneracy in the likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = NotQuiteSoSimplePrecessionModel()\n",
    "prior = qi.UniformDistribution([-1, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def trial(a=0.98, h=None, track=True):\n",
    "    true_params = prior.sample()\n",
    "    updater = qi.SMCUpdater(model, 1000, prior, resampler=qi.LiuWestResampler(a=a, h=h, debug=True), track_resampling_divergence=track, debug_resampling=True)\n",
    "    \n",
    "    for idx_exp in range(100):\n",
    "        exp = np.array([(9/8)**idx_exp], dtype=model.expparams_dtype)\n",
    "        datum = model.simulate_experiment(true_params, exp)\n",
    "        \n",
    "        updater.update(datum, exp)\n",
    "        \n",
    "    # Since this model is always degenerate about the origin, take the mean\n",
    "    # not of the original particles, but of the absolute value of the particles.\n",
    "    updater.particle_locations = np.abs(updater.particle_locations)\n",
    "        \n",
    "    if track:\n",
    "        return [true_params, updater.est_mean()], np.mean(updater.resampling_divergences)\n",
    "    else:\n",
    "        return [true_params, updater.est_mean()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def corr(a=0.98, h=None, track=True, n_trials=100):\n",
    "    trues = np.zeros((n_trials,))\n",
    "    ests  = np.zeros((n_trials,))\n",
    "    \n",
    "    if track:\n",
    "        divs = np.zeros((n_trials, ))\n",
    "    \n",
    "    for idx_trial in range(n_trials):\n",
    "        if track:\n",
    "            print(idx_trial, end=' ')\n",
    "            (true, est), div = trial(a, h, track)\n",
    "        else:\n",
    "            true, est = trial(a, h, track)\n",
    "        trues[idx_trial] = true\n",
    "        ests[idx_trial] = est\n",
    "        if track:\n",
    "            divs[idx_trial] = div\n",
    "    \n",
    "    if track:\n",
    "        return trues, ests, divs\n",
    "    else:\n",
    "        return trues, ests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trues, ests = corr(1, 0.02, False, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00697411114983 0.0123620858504\n"
     ]
    }
   ],
   "source": [
    "bias = np.mean(np.abs(trues) - np.abs(ests))\n",
    "risk = np.mean(np.abs(np.abs(trues) - np.abs(ests)))\n",
    "print(bias, risk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9Q0/f9B/BnIPxQiUoQihKslQBxbBMK/qh2X4uyXruu\nArXLdXXeaa11Tr12XbtWqFgd6tnqzVbnj7ZYrLXduDKKd6vtWTurcu0NB2wtyo+gtgJixfiDiL+S\nfL5/uHySQJBPIL/zfNz1lg/5JHnxOpaX798yQRAEEBEROSHE2wEQEZH/YfEgIiKnsXgQEZHTWDyI\niMhpLB5EROQ0Fg8iInKa3FMftH37dtTU1GDEiBHYuHGjw3t27dqFuro6REREYOnSpRg3bpynwiMi\nIid4rOWRnZ2NwsLCPp+vra3FuXPn8Oabb+KZZ57B22+/Lfm96+vrXRFiQGAurJgLK+bCirmwGkwu\nPFY8NBoNhg0b1ufz1dXVmDFjBgAgOTkZ3d3duHTpkqT35h+DFXNhxVxYMRdWzIWVXxSP/uj1esTE\nxIjXSqUSer3eixEREVFffKZ4EBGR//DYgHl/lEolLly4IF5fuHABSqXS4b319fV2zS2tVuv2+PwF\nc2HFXFgxF1bMhZVWq0VZWZl4nZaWhrS0NEmv9WjxEAQBfe3DmJWVhc8++wzTpk1DU1MThg0bhpEj\nRzq819Ev2N7e7vJ4/ZFCoUBXV5e3w/AJzIUVc2HFXFiNGTNmwMXUY8XjjTfewPHjx9HV1YUlS5ZA\nq9XCaDRCJpMhJycH9957L2pra7F8+XJERkZiyZIlngqNiIicJAuULdnZ8riN/6qyYi6smAsr5sJq\nzJgxA34tB8yJiMhpLB5EROQ0Fg8iInIaiwcRETmNxYOIiJzG4kFERE5j8SAiIqexeBARkdNYPIiI\nyGksHkRE5DQWDyKiAGAwGHDs2DEYDAaPfJ7PbMlOREQDYzAYkJ+fj6amJqSkpKCiogJRUVFu/Uy2\nPIiI/JjBYMDHH3+MpqYmGI1GNDc3o7Gx0e2fy5YHEZEPMBgMaGhogEajkdxqsLQ4GhsbIZff/jpP\nTk5GamqqO0MFwJYHEZHXWYrAnDlzkJubi8OHD6Ojo8NuDMPRmEZDQwOamppgMplgMpmwbt06j3RZ\nAWx5EBF5naUIGI1GNDQ0YO7cuQgLC4PJZEJKSgr27NmDJ598EjqdDmq1Gvv27UNUVBQ0Gg1SUlLQ\n3NyM5ORk5ObmeqRwACweREReZykClgJiNptx48YNAEBzczP+8Y9/iOMYjY2NqK2txc9+9jMAQFFR\nEQAgIyOjz8IxkC6x/rDbiojIy6KiolBRUYG9e/ciNTUVcrkcERERkMvlSE5ORmJiot39MplM7Or6\nzW9+gzVr1vT53rZdYvn5+S6bysuWBxGRD4iKikJ6ejpeeukltLW1YdKkSfjvf/+LWbNmiV1Ulm6r\n9PR0u64uywyrzMzMXu8r9T5nsXgQEfkAg8GA2bNni91TERERMBqNSE1NRUVFBSorK9HY2IjU1FSH\n4x19zbDSaDRISkqCTqfD+PHjXTYTi8WDiMhD7jT20NDQAJ1OJ17bjnlYWgu2LQZLV5dtQemLTCaz\n+19X4JgHEZEH9Df2oNFooFarxWvbMY++WgtRUVHIzMy8Y+GwFCWTyYSWlhaXLSBky4OIyAP6G3uI\niorCvn37UFtbC5lMBrVajba2tn5bFf2R2r3lLBYPIiI3MxgM6O7uhlqtRktLS59f4lFRUeIUXACI\nj48f9Gc7073lDBYPIiIJBrpWwnbTwqSkJLz//vtIT093+B7uWI8BWLu3XInFg4ioH4PZtda2u+rk\nyZMYMmRIn4XD0zvjDgYHzImI+uFovEIqy5hDWFjYHcccBvMZ3sCWBxFRPwYz6Cx1zMFdA9vuIhME\nQfB2EK7Q3t7u7RB8gkKhQFdXl7fD8AnMhRVzYTXQXBgMBpcPOnvjM2yNGTNmwK9ltxURkQRS1lQA\ngzsOVupn+AJ2WxERuYi/DXoPBlseREQu4m+D3oPB4kFE5CJSZ1YFAnZbERG5iLtWc/siFg8iIhdy\nx2puX8RuKyIKaoOZHRXMPNryqKurQ2lpKQRBQHZ2NvLy8uye7+7uxpYtW9DZ2Qmz2YxHH30UDzzw\ngCdDJKIgEkyzo1zNYy0Ps9mMkpISFBYWYtOmTaiqqkJbW5vdPZ999hkSExPx+uuvY9WqVXjvvfdg\nMpk8FSIRBZmGhgY0NjbCaDSiqakpoGdHuZrHiodOp8Po0aMRGxsLuVyO6dOno7q62u4emUyGa9eu\nAQCuX78OhUKB0NBQT4VIREFGpVJBLr/dARMaGoqEhATxOXZn3ZnHioder0dMTIx4rVQqodfr7e55\n6KGH0NraisWLF+PFF1/E/PnzPRUeEQWh1tZWsXfDZDKJvSH9nfpHPjbbqq6uDvfccw9WrVqFjo4O\nFBcXY+PGjYiMjLS7r76+HvX19eK1VquFQqHwdLg+KTw8nLn4H+bCirmwss1FVlYWNBqNOLU2MzMT\nCoUCx48ft1vsd+bMGUyePNnLkbtHWVmZ+DgtLQ1paWmSXuex4qFUKtHZ2Sle6/V6KJVKu3sOHTok\nDqLHx8cjLi4ObW1tSEpKsrvP0S/ITd9u4wZ4VsyFFXNh1TMX5eXlYvEAbn+XjB071m6H28TExIDM\nn0KhgFarHdBrPdZtpVar0dHRgfPnz8NoNKKqqgpZWVl294waNQrffPMNAODSpUs4e/Ys7rrrLk+F\nSER+ytWbEVoW+5WXl3MGVh88uiV7XV0d3n33XQiCgJkzZyIvLw8HDhyATCZDTk4OLl68iG3btuHi\nxYsAgLy8PNx///2S3ptbst/Gf2FaMRdWgZILR8e09pxuu2fPHrS2tvY6ytXy2p7/aA1mg9mSned5\nBJhA+ZJwBebCKhBy4WhNBgB8/PHHKCgogMlkglwuh0qlQmtrq926DdvXajQalJeXszUBnudBREGg\n5461dXV1yM/PR2FhIeRyuVg4zpw5I67bqKysFFscltc2NjZyPYcLsHgQkV/ouWOt2WwWC4LJZMK6\ndetQXl6O1NRUyOVyhIaGoqCgAPn5+VCpVOJrU1NTA3q3W0/xqam6RER96bljLQC7GVG5ubniPZWV\nlSgoKBBbKW1tbeJrg2HTQk/gmEeACYS+bVdhLqwCMRcGgwE1NTUICQlBenp6r8Hx/Px8sbDYzpgK\nxFwM1GDGPNjyICK/09fguUUwnavhLRzzICK/I+W4V0frN8h1WDyIyKt6LvCTsuAvmI579VXstiIi\nr3G0wG/evHkOz9fouUCQ3VLexZYHEXlNz+6ngwcPOuyOcrTLLbulvIvFg4i8pmf306xZsxx2R0kZ\n4yDPYrcVEXmNo+4nR91RliJjmXrLMQ7v4zqPAMM57FbMhVUg5MJgMLhkjCMQcuEqXOdBRAHPMsZB\nvoFjHkTkETwTPLCweBCR29nOlsrNzcXhw4dZRPwciwcRuZ3tbKmGhgbMnTsXDz74ID777DMWET/F\n4kFEbmUwGHDt2jXcc889CA0NBQCYzWZ89913eOqppzB79mwWED/EAXMichtLd1VjYyPkcjkEQUBY\nWBhu3bol3qPT6bhVuh9iy4OI3MbSXWUymXDjxg2YzWYIgoC77rpLvEetVnPdhh9iy4OI3MayuK+p\nqQmhoaEwmUziHlY6nQ6CICAjI4NbjPghFg8ichvbFeMJCQloa2sTF/nFx8d7OzwaBBYPInIr28V9\nLBiBg2MeRETkNBYPIiJyGosHERE5jcWDiIicxuJBREROY/EgIiKnsXgQEZHTWDyIiMhpLB5EROQ0\nFg8iInIaiwcR8YhYchqLB1EQcVQkbI+Izc/PZwEhSVg8iIJEX0XC9ojY5uZmNDY2ejlS8gd33FW3\nqKgIMpms3zdZvXq1pA+rq6tDaWkpBEFAdnY28vLyet1TX1+P3bt3w2QyYfjw4Vi1apWk9yaiO3NU\nJDIzM8UzN5qbm5GcnMyDmUiSOxaPmTNnio/PnTuHf/7zn5gxYwZiY2PR2dmJL7/8EtnZ2ZI+yGw2\no6SkBEVFRYiOjsaKFSswadIkJCQkiPd0d3ejpKQEr7zyCpRKJa5cuTLAX4uIeuqrSNieuWE5a4Oo\nP3csHg888ID4uLCwEIWFhUhMTBR/dv/992P79u3QarX9fpBOp8Po0aMRGxsLAJg+fTqqq6vtisfR\no0cxZcoUKJVKAMDw4cOd+mWIqG93KhK2Z24QSSH5MKjW1la7c4cBIC4uDm1tbZJer9frERMTI14r\nlUrodDq7e9rb22EymbB69Wpcv34dDz/8MP7v//5PaohE5IDBYMDx48cxduxYh0XCYDCgoaEBGo2G\nrQ6STPKA+Y9+9CNs27YNZ8+exc2bN9He3o7t27dDo9G4LBiz2YxTp05hxYoVKCgoQHl5OTo6Olz2\n/kTBxjJI/vDDDzucScWZVjRQklseS5cuxTvvvIPnn38eZrMZoaGhmDx5Mn73u99Jer1SqURnZ6d4\nrdfrxe4p23sUCgXCw8MRHh6OCRMm4PTp072Orqyvr0d9fb14rdVqoVAopP4qAS08PJy5+B/mAjh+\n/LjdIPmZM2cwefJkyc8HIv5d2CsrKxMfp6WlIS0tTdLrJBePqKgoPPfcczCbzbhy5QqGDx+OkBDp\nM33VajU6Ojpw/vx5REdHo6qqCs8++6zdPZMmTcKuXbtgNptx69YtNDc345e//GWv93L0C3Z1dUmO\nJZApFArm4n+CPRcGgwGdnZ1Qq9VoaWlBcnIyEhMT7XIyduxYu0H0ns8HomD/u7ClUCgkjVk7Irl4\nAEBbWxu++uorXL58GQsXLkR7eztu3bqFu+++u9/XhoSEYOHChSguLoYgCJg5cyZUKhUOHDgAmUyG\nnJwcJCQkYOLEiXjhhRcQEhKCnJwcqFSqAf1iRMHM0h3V1NSEpKQklJeXO5xJxZlWNFAyQRAEKTd+\n9dVXeOeddzBlyhRUVVVh9+7daGlpwQcffICVK1e6O85+tbe3ezsEn8B/VVkFcy6OHTuGOXPmwGg0\nIiwsDPv378eECRO8HZZPCOa/i57GjBkz4NdK7ncqKyvDypUr8cwzz4jdVXfffTdOnz494A8nIvew\nrOkICwtDcnIyCwe5nORuq8uXL/fqnpLJZJJWoBORZ/XsjuK/tsnVJLc8xo8fj8OHD9v9rKqqCmq1\n2uVBEVH/+tsJ17Kmg+MY5A6SWx4LFixAcXExvvjiC9y4cQNr165Fe3s7XnnlFXfGR0QO2A6Ip6Sk\noKKigkWCPEpy8UhISMDmzZvx73//G5mZmYiJiUFmZiYiIyPdGR8ROdDXJodEniK522rXrl2IiIjA\ntGnTMHv2bEyfPh2RkZEoLS11Y3hEwcvSLdXR0dGre6rngDh3wiVPk9zy+PLLL/HUU0/1+vnhw4cx\nf/58V8ZEFLQs+0ypVCrMmzcPTU1NCA0NhdFoRGpqqtg9xfUZ5G39Fo8vvvgCAGAymcTHFj/88AOX\n+RO5iMFgwOzZs6HT6ZCQkID29nYYjUYYjUYA6NU9xZ1wyZv6LR5HjhwBABiNRvGxxYgRI7B06VL3\nREYUZGpqasRT/L7//nskJiaivb0dcrkcJpOJ3VPkU/otHpaT/P7617/iiSeecHtARMHq+vXrdteW\nzR/Gjh2L4uJipKens3uKfIbkMQ/bwiEIAmx3NXFmg0Qi6s1gMGDDhg3itUqlwtmzZ2EymXD69GkM\nGTKEhYN8iuTiodfrUVJSghMnTuDq1at2z/3tb39zeWBEwaShoUE8HE0ul6O4uBivvfYazxUnnyW5\nyfDWW29BLpejqKgIkZGR2LBhA7KysrBo0SJ3xkcUFGyn3qakpOC+++5DRUUFysvLuQCQfJLklkdT\nUxO2bduGyMhIyGQyjBs3DkuWLMErr7yCnJwcd8ZIFPD6mnrL2VTkqyS3PEJCQhAaGgoAGDZsGK5c\nuYKIiAjo9Xq3BUcUSLgXFQUSyS0PtVqN2tpaTJ48GRMnTsSf//xnhIeHIykpyZ3xEQUE7kVFgUZy\n8Vi+fLk4w2r+/PnYt28frl+/jkceecRtwREFCu5FRYFGcvEYNmyY+Dg8PByPP/64WwIi8meW7UU0\nGo1dy8IyIM7ZUxQoJBcPk8mEqqoqnDp1qtdipsWLF7s8MCJ/c6euKe5FRYFGcvHYsmULvv/+e6Sn\np2PEiBHujInIL/XXNcW9qCiQSC4edXV12L59O4YMGeLOeIj8FrumKJhILh6JiYkwGAwsHkR9YNcU\nBRPJxWPZsmXYsWMHJk6c2KvbasaMGS4PjMgfsWuKgoXk4nHo0CE0NDTg6tWrCA8PF38uk8lYPCgg\n9Zw5ZTAYUFNTA5lMhoyMDLYsKKhJLh6ffPIJNmzYAJVK5c54iHxCz5lTe/bswZNPPimet6HRaFBZ\nWckCQkFL8vYkI0eOxKhRo9wZC5HP6Dlz6uDBg+Kut8Dtvd7q6uq8GCGRd0kuHo888gi2bNmCpqYm\nnDt3zu4/okBgMBhw+PBhHDlyBCqVStzlNjk5GbNmzYJarRbvNZvNKCoq6nOfKqJAJ7nbqqSkBABw\n7NixXs/xPA/yd5bzw227pfbu3Yu2tjZx5tQHH3yAHTt2YNeuXTCZTDh58iS3GaGgJbl4sEBQILM9\njAkAdDod2traxMJgMBgwb948NDY2Qi6XQyaTcS0HBTWeH0uE2y0N224ptVptVxgsYyAmkwkmkwnr\n1q3jzrgU1O7Y8li7di0KCwsBAEVFRZDJZA7vW716tesjI/KgqKgo7Nu3D7W1tZDJZEhPT7/jxoa5\nubksHBTU7lg8bNdvzJw50+3BEHlTVFQUfvazn/X5XEVFBWpraz0cFZFvumPxuP/++8XHCQkJSE5O\n7nWPbT8xkT/q6OjA559/jpycHMTHx9/x3jVr1vBAJyI4MeZRXFzs8Odr1651WTBE7uToGNiWlhZM\nnToVL730EqZNm4aOjo4+X+9o11yiYNVv8TCbzTCbzRAEAYIgiNdmsxlnz54VzzUn8mWWFeNz5sxB\nfn4+DAYDDAYDtFotbt26BQC4ceMGDh482Od7WMY9LGs/ONOKglm/U3V//etfi4+feOIJu+dCQkKQ\nn5/v+qiIXMxRq0EQBJw/f168JywsDLNmzerzPbhrLpFVv8Vj69atEAQBr776qt2sKplMhuHDh9tt\nktifuro6lJaWQhAEZGdnIy8vz+F9Op0OK1euxHPPPYcpU6ZIfn+ivvR11kZqaioaGxsRGxuLsrKy\nfsc8uGsu0W39Fo/Y2FgAwLZt2+x+fvPmzT6n7jpiNptRUlKCoqIiREdHY8WKFZg0aRISEhJ63ffB\nBx9g4sSJkt+bqD99tRrYkiAaGMkD5u+99544s6qmpgYLFizAggULHG5X4ohOp8Po0aMRGxsLuVyO\n6dOno7q6utd9n376KaZOnYrhw4dLDY1IEkurwbZIOPoZEfVPcvE4evQoEhMTAQAfffQRli9fjj/+\n8Y/48MMPJb1er9cjJiZGvFYqldDr9b3uqa6uxoMPPig1LKJeLLOqzp49iyNHjuDw4cPcwJDIxSTv\nbXXjxg1ERESgq6sL586dw9SpUwEAnZ2dLgumtLQUc+fOFa8FQXDZe1NwsD2HQyaTiTOpUlNTsW/f\nPrYwiFxEcvEYM2YMjhw5go6ODvz0pz8FAFy5ckXygLlSqbQrNHq9Hkql0u6ekydPYvPmzRAEAV1d\nXaitrYVcLkdWVpbdffX19aivrxevtVotFAqF1F8loIWHhwdtLrq6uvDpp5+Ks6ps6XQ6nDlzBpMn\nT/ZSdN4VzH8XPTEX9srKysTHaWlpSEtLk/Q6ycVj4cKFKC0thVwux29/+1sAwH/+8x+xkPRHrVaj\no6MD58+fR3R0NKqqqvDss8/a3bN161bx8bZt25CZmdmrcACOf8Guri6pv0pAUygUQZkLS4ujsbER\nISG9e2PVajUSExODMjdA8P5dOMJcWCkUCmi12gG9VvKYh1qtxq9+9SvExcWhtLQUwO3WSHZ2trQP\nCgnBwoULUVxcjOeffx7Tp0+HSqXCgQMH8Pnnnw8oeApePVeL2+56azabER8fj9DQUIwdOxalpaXs\nsiJyMZkgcWBh//79+OSTTzBr1ixUVFRg9+7dOHPmDHbu3Nnn1iWe1N7e7u0QfEIw/Kuq5/niFRUV\nAID8/HxxHceePXtw8eJFJCYmsmggOP4upGIurMaMGTPg10rutvrkk0+wcuVKxMXFobKyEsDtzRL5\npU2e5mi1eGZmZq81G8nJyfySIHITyd1W165dw6hRo+x+ZjQaIZdLrj9ELtHXHlNcs0HkOZK/+SdM\nmICPP/4Yjz32mPiz/fv3Sx6ZJ3IV7jFF5H2SxzwuXryIDRs2oKurC3q9HnFxcRgyZAhefvlljBw5\n0t1x9ovdZ7cFYn+uwWBAQ0MDNBqNU4UiEHMxUMyFFXNh5ZExj+joaKxfvx4tLS04f/48YmJioFar\nHU6LJHIVR4PjbGkQeZ9TAxYymQxqtRpqtdpd8RABsLY2Lly4gMbGRphMJrvBcSLyLo52k8+xbW2E\nhobCZDIBAMaPH88DmIh8BIsH+RzbqbiWbUbkcjnWrFnDLisiH8EBC/I5lqm4crkcERERkMvlSElJ\nQXp6urdDI6L/YcuDvKavWVS2U3ETEhLQ1tbGKblEPobFg7yiv1lUtse99nc0LBF5HrutyCscbTFC\nRP6DxYO8oq8tRojIP7DbiryCW4wQ+Te2PMjtep69YcGNDIn8F1se5FbcXoQoMLHlQW7FgXGiwMTi\nQW7FgXGiwMRuK3IrDowTBSYWD3I72wV/RBQY2G1FA9bXLCoiCnxsedCAcBYVUXBj8SCnGAwG1NTU\niDOneEgTUXBi8aA7st35FgBmz54tTreNiIiATCbjLCqiIMTiQX3q2TW1cuVKNDc3i88bjUasX78e\nubm57LIiCjIcMKc+1dbWorGxUVzgd+PGDcjl1n9vjB8/noWDKEixeJBDBoMBq1atEs8PT0pKQmRk\nJMxmMwAgNDQUf/rTn1g4iIIUiwc51NDQgJaWFgC3C8Wrr76KjIwMcbV4amoqMjIyvBwlEXkLxzyo\n16B4Q0MDVCoVUlJS0NzcjOTkZGRkZHC1OBGJWDyCnO2guFqthiAIaGlpQUpKCvbs2dPr/HCuFici\ngMUj6NnueqvT6SAIgrh2o62tjYWCiBzimEcQMxgM6O7uhlqtRlhYGNRqtfiYazeI6E7Y8ggylvEN\nlUqFefPmoampCUlJSXj//feRnp4OABzTIKJ+sXgEuJ6D4ZbxDZVKhTNnzsBkMuHkyZMYMmSIWCzY\nVUVE/WHxCGA9B8OffPJJcT+q1tZWJCYmoq2tjV1UROQ0Fo8AZjsY3tDQgFdffRVhYWGQyWR9zqYi\nIpLCo8Wjrq4OpaWlEAQB2dnZyMvLs3v+6NGjqKysBABERkZi0aJFGDt2rCdDDCgajQZJSUniRoZm\ns7nXflTx8fFejpKI/JHHZluZzWaUlJSgsLAQmzZtQlVVFdra2uzuiYuLw+rVq/H6669jzpw52Llz\np6fCC0hRUVFYvXo1QkNDxZ8lJiZyPyoiGjSPFQ+dTofRo0cjNjYWcrkc06dPR3V1td09KSkpGDp0\nKAAgOTkZer3eU+EFrIyMDKSmpkIul2PcuHEoLy9n4SCiQfNYt5Ver0dMTIx4rVQqodPp+rz/4MGD\n4tRRGjhuKUJE7uCTA+bffvstDh06hDVr1jh8vr6+HvX19eK1VquFQqHwVHg+LTw8vFcuFAoFRo8e\n7aWIvMdRLoIVc2HFXNgrKysTH6elpSEtLU3S6zxWPJRKJTo7O8VrvV4PpVLZ677vvvsOb731FgoK\nCvr8V7KjX7Crq8u1Afs42/UbtnlSKBRBl4u+MBdWzIUVc2GlUCig1WoH9FqPjXmo1Wp0dHTg/Pnz\nMBqNqKqqQlZWlt09nZ2d2LRpE5YtW8ZZQHdgWb8xZ84c5Ofnw2AweDskIgoyHmt5hISEYOHChSgu\nLoYgCJg5cyZUKhUOHDgAmUyGnJwcfPTRRzAYDCgpKYEgCAgNDcX69es9FaLfsF2/0dzcjMbGRq4K\nJyKPkgmCIHg7CFdob2/3dggeY2l5WM7aqKioELuu2CS3Yi6smAsr5sJqzJgxA36tTw6Y051xBhUR\neRuLh5/ioUxE5E08z4OIiJzG4kFERE5j8fARBoMBx44d47RbIvILHPPwAbbnbqSkpNjNniIi8kVs\nefgAR+s2iIh8GYuHD9BoNEhJSUFYWBhP9SMiv8BuKx/AdRtE5G/Y8vACR4PjlnUbLBxE5A/Y8vAw\nDo4TUSBgy8PDODhORIGAxcPDODhORIGA3VYexsFxIgoELB5ewE0NicjfsdvKhbjFCBEFC7Y8Bsly\nlrhKpcK8efM4i4qIggKLxyDYTrtVqVRobW3l0bBEFBTYbTUIttNuW1tboVKpOIuKiIICWx6DYJl2\nazlLfM+ePWhra+MsKiIKeCweTrCMb2g0GkRFRTmcdhsfH+/tMImI3I7FQ6K+thXhtFsiCkYc85CI\n24oQEVmxeNyB7boNbitCRGTFbqse7rRug9uKEBHdxuJhQ8q6DY5vEBGx28oO120QEUnDlocNrtsg\nIpKGxcMG120QEUkTlN1Wd9r9lmeJExH1L+haHjxDnIho8IKu5cHFfkREgxd0xYOL/YiIBi/gu62k\nbGZIRETOCejiwc0MiYjcw6PFo66uDqWlpRAEAdnZ2cjLy+t1z65du1BXV4eIiAgsXboU48aNc/pz\nLK2Na9eu9RrfYNEgIho8jxUPs9mMkpISFBUVITo6GitWrMCkSZOQkJAg3lNbW4tz587hzTffRHNz\nM95++22sXbtW0vsbDAbU1taiu7sbr732GnQ6HZKSkqBWq9HS0sLxDSIiF/JY8dDpdBg9ejRiY2MB\nANOnT0d1dbVd8aiursaMGTMAAMnJyeju7salS5cwcuTIft8/NzcXDQ0Ndj87efIk9uzZg6FDh3J8\ng4jIhTxWPPR6PWJiYsRrpVIJnU7X7z16vV5S8Whubra7lsvlSE5ORkZGBosGEZGLBcxU3eTkZLvH\ne/fu5QLrFzbHAAAHBklEQVRAIiI38VjLQ6lUorOzU7zW6/VQKpW97rlw4YJ4feHChV73AEB9fT3q\n6+vFa61WixMnTrghav+kUCi8HYLPYC6smAsr5sKqrKxMfJyWloa0tDRJr/NYy0OtVqOjowPnz5+H\n0WhEVVUVsrKy7O7JysrCl19+CQBoamrCsGHDHHZZpaWlQavViv/Z/vLBjrmwYi6smAsr5sKqrKzM\n7rtUauEAPNjyCAkJwcKFC1FcXAxBEDBz5kyoVCocOHAAMpkMOTk5uPfee1FbW4vly5cjMjISS5Ys\n8VR4RETkBI+u80hPT8cbb7xh97Of//zndtcLFy70ZEhERDQAATFg7kxTK9AxF1bMhRVzYcVcWA0m\nFzJBEAQXxkJEREEgIFoeRETkWSweRETkNL/aVddTGyv6g/5ycfToUVRWVgIAIiMjsWjRIowdO9Yb\nobqdlL8L4PYWOStXrsRzzz2HKVOmeDhKz5CSi/r6euzevRsmkwnDhw/HqlWrvBCp+/WXi+7ubmzZ\nsgWdnZ0wm8149NFH8cADD3gnWDfavn07ampqMGLECGzcuNHhPQP63hT8hMlkEpYtWyb88MMPwq1b\nt4QXXnhBaG1ttbunpqZGWLdunSAIgtDU1CQUFBR4I1S3k5KLxsZG4erVq4IgCEJtbW1Q58Jy3+rV\nq4X169cLX3/9tRcidT8pubh69arw+9//Xrhw4YIgCIJw+fJlb4TqdlJy8fe//13Yu3evIAi387Bg\nwQLBaDR6I1y3OnHihHDq1CnhD3/4g8PnB/q96TfdVrYbK8rlcnFjRVt9bawYaKTkIiUlBUOHDgVw\nOxd6vd4bobqdlFwAwKeffoqpU6di+PDhXojSM6Tk4ujRo5gyZYq4c0Og5kNKLmQyGa5duwYAuH79\nOhQKBUJDQ70RrltpNBoMGzasz+cH+r3pN8Wjr00Tnb0nEDj7ex48eBDp6emeCM3jpP5dVFdX48EH\nH/R0eB4lJRft7e0wGAxYvXo1VqxYgcOHD3s6TI+QkouHHnoIra2tWLx4MV588UXMnz/fw1H6hoF+\nb/pN8aCB+fbbb3Ho0CHMnTvX26F4TWlpqd3vLwTx7HSz2YxTp05hxYoVKCgoQHl5OTo6OrwdllfU\n1dXhnnvuwc6dO7FhwwaUlJTg+vXr3g7Lb/jNgLkrN1b0d1JyAQDfffcd3nrrLRQUFATs7sJScnHy\n5Els3rwZgiCgq6sLtbW1kMvlvfZW83dS/z+iUCgQHh6O8PBwTJgwAadPn0Z8fLynw3UrKbk4dOiQ\nOIgeHx+PuLg4tLW1ISkpyaOxettAvzf9puXhyo0V/Z2UXHR2dmLTpk1YtmxZwH0x2JKSi61bt2Lr\n1q34y1/+gqlTp+Lpp58OuMIBSMvFpEmT0NDQALPZjBs3bqC5uRkqlcpLEbuPlFyMGjUK33zzDQDg\n0qVLOHv2LO666y5vhOt2giD02eIe6PemX60wr6urw7vvviturJiXl2e3sSIAlJSUoK6uTtxYcfz4\n8V6O2j36y8WOHTvwr3/9C7GxsRAEAaGhoVi/fr23w3YLKX8XFtu2bUNmZmZAT9XtLxf79u3DoUOH\nEBISglmzZuHhhx/2ctTu0V8uLl68iG3btuHixYsAgLy8PNx///1ejtr13njjDRw/fhxdXV0YMWIE\ntFotjEbjoL83/ap4EBGRb/CbbisiIvIdLB5EROQ0Fg8iInIaiwcRETmNxYOIiJzG4kFERE5j8SAi\nIqexeBANwNKlS/Htt996Owwir2HxIHIxs9ns7RCI3I4rzImctHXrVhw5cgRhYWEIDQ3FnDlzsHfv\nXixevBgfffQR4uLioNVqsWXLFmzfvl183dKlS7FkyRL8+Mc/hiAIqKysxMGDB9Hd3Y2f/OQnWLRo\n0R3PXSDyJWx5EDlp2bJlGDVqFF5++WXs3r0b9913HwDgxIkT2Lx5MwoLC/t9j/379+PYsWNYs2YN\ndu7ciWHDhuGdd95xd+hELsPiQeQiWq0W4eHhCAsL6/feAwcO4IknnkB0dDTkcjkef/xxfP311+zy\nIr/hN+d5EPk629PY+tPZ2YmNGzdCJpOJP5PL5bh8+TKio6PdER6RS7F4EA2A7Ze+IxEREbh586Z4\nbTabceXKFfF61KhRWLJkCVJSUtwWI5E7sduKaABGjhyJc+fOAXB8rO3o0aNx8+ZN1NbWwmQyoby8\nHEajUXw+JycHH374oXja3ZUrV3Ds2DHPBE/kAmx5EA1AXl4edu3ahffffx+PPfZYr+eHDh2Kp59+\nGjt27IDZbEZubq5dt9YvfvELAEBxcTEuXryIESNGYNq0aQF5wiEFJk7VJSIip7HbioiInMbiQURE\nTmPxICIip7F4EBGR01g8iIjIaSweRETkNBYPIiJyGosHERE5jcWDiIic9v8OoNIDoElgtQAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x17a9d5bfba8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.abs(trues), ests, 'k.')\n",
    "plt.xlabel('true')\n",
    "plt.ylabel('estimated')\n",
    "plt.xlim((0, 1))\n",
    "plt.ylim((0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.19281919156\n"
     ]
    }
   ],
   "source": [
    "trues, ests = corr(0.98, None, False, 100)\n",
    "bias = np.mean(np.abs(trues) - np.abs(ests))\n",
    "print(bias)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X9Qk/cdB/B3IAJVghJEEYK1JSCWrWLFH1M7q9KuXdeK\npcu5s73VUWetdnVbuxWtWDys15u91eq0v7A4+2Pjylq8XVuPs7OtzK5YyNaiAQFtRYwVU5WIoiTP\n/rAkBELyPPnx5Nf7deddQp4kn3yE55Pvz0chCIIAIiIiCaICHQAREYUeFg8iIpKMxYOIiCRj8SAi\nIslYPIiISDIWDyIikkwp1xvt2LED9fX1GDlyJDZv3uz0mJ07d0Kv1yM2NhYrV67EhAkT5AqPiIgk\nkK3lMW/ePKxdu3bIxxsaGnD69Gm8+OKL+PWvf41XX31V9Gs3Njb6IsSwwFzYMRd2zIUdc2HnTS5k\nKx7Z2dkYMWLEkI/X1dVh7ty5AIDMzEx0d3fj3Llzol6bvwx2zIUdc2HHXNgxF3YhUTzcMZlMSEpK\nst1Xq9UwmUwBjIiIiIYSNMWDiIhCh2wD5u6o1WqcPXvWdv/s2bNQq9VOj21sbHRobul0Or/HFyqY\nCzvmwo65sGMu7HQ6HSorK233c3JykJOTI+q5shYPQRAw1D6MeXl52Lt3L2bNmoXm5maMGDECo0aN\ncnqssw/Y0dHh83hDkUqlQldXV6DDCArMhR1zYcdc2KWmpnpcTGUrHlu2bMHhw4fR1dWFFStWQKfT\nobe3FwqFAvn5+bjlllvQ0NCAxx57DHFxcVixYoVcoRERkUSKcNmSnS2Pa/ityo65sGMu7JgLu9TU\nVI+fywFzIiKSjMWDiIgkY/EgIiLJWDyIiEgyFg8iIpKMxYOIiCRj8SAiIslYPIiISDIWDyIikozF\ng4iIJGPxICIiyVg8iIhIMhYPIiKSjMWDiIgkY/EgIiLJWDyIiEgyFg8iIpKMxYOIiCRj8SAiIslY\nPIiISDIWDyIikozFg4iIJGPxICIiyVg8iIhIMhYPIiKSjMWDiIgkY/EgIiLJWDyIiEgyFg8iIpKM\nxYOIiCRj8SAiIslYPIiISDIWDyIikozFg4iIJGPxICIiyVg8iIhIMqWcb6bX61FRUQFBEDBv3jwU\nFBQ4PN7d3Y2tW7eis7MTVqsV99xzD2677TY5QyQiIhFkKx5WqxXl5eUoKSlBYmIiiouLMW3aNKSl\npdmO2bt3L9LT0/HHP/4RFy5cwOrVq3HrrbciOjparjCJiEgE2bqtWlpaMG7cOCQnJ0OpVGL27Nmo\nq6tzOEahUODSpUsAgMuXL0OlUrFwEBEFIdmKh8lkQlJSku2+Wq2GyWRyOObOO+9Ee3s7li9fjief\nfBIPPfSQXOEREZEEso55uKPX63HDDTdg/fr1MBqNKCsrw+bNmxEXF+dwXGNjIxobG233dTodVCqV\n3OEGpZiYGObie8yFHXNhx1w4qqystN3OyclBTk6OqOfJVjzUajU6Oztt900mE9RqtcMx+/fvtw2i\np6SkYMyYMTh58iQyMjIcjnP2Abu6uvwUeWhRqVTMxfeYCzvmwo65sFOpVNDpdB49V7ZuK61WC6PR\niDNnzqC3txe1tbXIy8tzOGb06NH48ssvAQDnzp3DqVOnMHbsWLlCJCIikWRreURFRaGoqAhlZWUQ\nBAHz58+HRqNBTU0NFAoF8vPzUVhYiO3bt+OJJ54AACxZsgTx8fFyhUhERCIpBEEQAh2EL3R0dAQ6\nhKDAJrkdc2HHXNgxF3apqakeP5crzImISDIWDyIikozFg4iIJGPxICIiyVg8iIhIMhYPIiKSjMWD\niIgkY/EgIiLJWDyIiEgyFg8iIpKMxYOIiCRj8SAiIslYPIiISDIWDyIikozFg4iIJGPxICIiyVg8\niIhIMhYPIiKSjMWDiIgkY/EgIiLJWDyIiEgyFg8iIpKMxYOIiCRj8SAiIslYPIiISDIWDyIikozF\ng4iIJGPxICIiyVg8iIhIMhYPIiKSjMWDiIgkY/EgIiLJWDyIiEgyFg8iIpKMxYOIiCRTunqwpKQE\nCoXC7YuUlpaKejO9Xo+KigoIgoB58+ahoKBg0DGNjY3YtWsXLBYLEhISsH79elGvTURE8nFZPObP\nn2+7ffr0afzrX//C3LlzkZycjM7OTnz88ceYN2+eqDeyWq0oLy9HSUkJEhMTUVxcjGnTpiEtLc12\nTHd3N8rLy/H0009DrVbjwoULHn4sIiLyJ5fF47bbbrPdXrt2LdauXYv09HTbz+bMmYMdO3ZAp9O5\nfaOWlhaMGzcOycnJAIDZs2ejrq7OoXgcOHAAM2bMgFqtBgAkJCRI+jBERCQPl8Wjv/b2dowdO9bh\nZ2PGjMHJkydFPd9kMiEpKcl2X61Wo6WlxeGYjo4OWCwWlJaW4vLly7jrrrvw4x//WGyIREQkE9ED\n5jfddBO2b9+OU6dO4cqVK+jo6MCOHTuQnZ3ts2CsViuOHTuG4uJirFmzBlVVVTAajT57fSIi8g3R\nLY+VK1fitddew+9+9ztYrVZER0dj+vTpePTRR0U9X61Wo7Oz03bfZDLZuqf6H6NSqRATE4OYmBhM\nmjQJx48fR0pKisNxjY2NaGxstN3X6XRQqVRiP0pYi4mJYS6+x1zYMRd2zIWjyspK2+2cnBzk5OSI\nep7o4hEfH4/Vq1fDarXiwoULSEhIQFSU+Jm+Wq0WRqMRZ86cQWJiImpra/H44487HDNt2jTs3LkT\nVqsVV69exdGjR/Gzn/1s0Gs5+4BdXV2iYwlnKpWKufgec2HHXNgxF3YqlUrUmLUzoosHAJw8eRIH\nDx7E+fPnUVRUhI6ODly9ehXXX3+92+dGRUWhqKgIZWVlEAQB8+fPh0ajQU1NDRQKBfLz85GWlobJ\nkyfjiSeeQFRUFPLz86HRaDz6YERE5D8KQRAEMQcePHgQr732GmbMmIHa2lrs2rULra2teOutt7Bu\n3Tp/x+lWR0dHoEMICvxWZcdc2DEXdsyFXWpqqsfPFd3yqKysxLp16zBhwgQcPHgQAHD99dfj+PHj\nHr85ERGFJtGDFufPnx/UPaVQKEStQCciovAiunjceOON+OSTTxx+VltbC61W6/OgiIgouInutlq6\ndCnKysrw0UcfoaenBxs3bkRHRweefvppf8ZHRERBSPSAOQD09PTgiy++QGdnJ5KSkjB16lTExcX5\nMz7ROGB+DQcD7ZgLO+bCjrmw82bAXHS31c6dOxEbG4tZs2bh3nvvxezZsxEXF4eKigqP35yIiEKT\n6OLx8ccfO/35wHEQIiIKf27HPD766CMAgMVisd3u8+2333KZPxFRBHJbPD799FMAQG9vr+12n5Ej\nR2LlypX+iYyIiIKW2+LRdyW/v/3tb1i8eLHfAyIiouAneqpu/8IhCAL6T9KSskEiERGFPtHFw2Qy\noby8HEeOHMHFixcdHvv73//u88CIiCh4iW4yvPLKK1AqlSgpKUFcXByee+455OXlYdmyZf6Mj4iI\ngpDo4tHc3IwVK1ZgwoQJUCgUmDBhAlasWIF//vOf/oyPiIiCkOjiERUVhejoaADAiBEjcOHCBcTG\nxsJkMvktOCIiCk6ixzy0Wi0aGhowffp0TJ48GX/+858RExODjIwMf8ZHRERBSPTeVhcvXoQgCIiP\nj8eVK1ewZ88eXL58GXfffTcSExP9Hadb3NvqGu7bY8dc2DEXdsyFnSwXgxoxYoTtdkxMDO6//36P\n35SIiEKb6OJhsVhQW1uLY8eO4fLlyw6PLV++3OeBERFR8BJdPLZu3YpvvvkGubm5GDlypD9jIiKi\nICe6eOj1euzYsQPXXXedP+MhChpmsxkGgwHZ2dmIj48PdDhEQUX0VN309HSYzWZ/xkIUNMxmMxYt\nWoTCwkIsWrSIv/tEA4hueaxatQovvfQSJk+ePKjbau7cuT4PjCiQDAYDmpub0dvbi6NHj6KpqQlT\np04NdFhEQUN08di/fz8MBgMuXryImJgY288VCgWLB4Wd7OxsZGVl4ejRo8jMzMTEiRMDHRJRUBFd\nPN5//30899xz0Gg0/oyHKCjEx8fj3XffRVNTEyZOnMgxD6IBRI95jBo1CqNHj/ZnLERBJT4+HlOn\nTmXhILfMZjMOHToUUWNjolsed999N7Zu3YqFCxcOGvMYO3aszwMjIgoFfZMrmpubkZWVhXfffTci\nvnCILh7l5eUAgEOHDg16jNfzIKJIFamTK0QXDxYIIooUUtb4ROrkCtHFg4goEkjthorUyRUui8fG\njRuxdu1aAEBJSQkUCoXT40pLS30fGRFRANTX16OpqQkWi0V0N1Tf5IpI4rJ49F+/MX/+fL8HQ0QU\nSGazGaWlpbBYLACAG2+80WU3VCRvYeOyeMyZM8d2Oy0tDZmZmYOOaWlp8X1UREQBYDAYbOc0pVKJ\nDRs2DFkUInWWVR/R6zzKysqc/nzjxo0+C4aIKJD6Br+HDRuGrKws5ObmDnls/1lWzc3NqK6ujqh1\nHm6Lh9VqhdVqhSAIEATBdt9qteLUqVO265oTUWBE4gI1f+kb/K6qqnLbkugrNEqlEtHR0SguLo6o\nTTTdzrb6xS9+Ybu9ePFih8eioqKwaNEi30dFRKJEeteJP4gd/O4rNNXV1SguLpY0wB4OYyVui8e2\nbdsgCAKeeeYZh1lVCoUCCQkJDpskuqPX61FRUQFBEDBv3jwUFBQ4Pa6lpQXr1q3D6tWrMWPGDNGv\nTxRpInWBWrCIj4/HwoULUVFRIXqdR7gUfLfFIzk5GQCwfft2h59fuXJlyKm7zlitVpSXl6OkpASJ\niYkoLi7GtGnTkJaWNui4t956C5MnTxb92kSRKhQXqIXDt+7+pK7zCJeCL3rA/K9//attFkJ9fT2W\nLl2KpUuXOt2uxJmWlhaMGzcOycnJUCqVmD17Nurq6gYd9+GHH2LmzJlISEgQGxpRxJLSR++OHGMn\n4XqRLSmbaPYflHdV8IN9LEt08Thw4ADS09MBAO+88w4ee+wx/OEPf8Dbb78t6vkmkwlJSUm2+2q1\nGiaTadAxdXV1uOOOO8SGRRTxxJ64XJ2M5DqpO/vWLZU/TqpynqjFFPxQKLKii0dPTw9iY2PR1dWF\n06dPY+bMmbj55pvR2dnps2AqKiqwZMkS231BEHz22kSRrKury+XJyBcndTHEfuseij9OqoE4Ubsr\n+HL9f3hD9N5Wqamp+PTTT2E0GnHzzTcDAC5cuCB6wFytVjsUGpPJBLVa7XBMW1sbXnjhBQiCgK6u\nLjQ0NECpVCIvL8/huMbGRjQ2Ntru63Q6qFQqsR8lrMXExDAX3/NVLrq6unD48GHcdNNNQZdbsbHV\n19c7nIxOnDiB6dOn2x7Py8tDdna2rd9+6tSpfvmsKpUKNTU1OHLkCCZNmiT5PQ4fPjzk5xCbi4G/\nF65eM1Dk+v8AgMrKStvtnJwc5OTkiHqe6OJRVFSEiooKKJVKPPLIIwCA//73v7ZC4o5Wq4XRaMSZ\nM2eQmJiI2tpaPP744w7HbNu2zXZ7+/btmDp16qDCATj/gF1dXWI/SlhTqVTMxfd8kYtgnhkjJbbM\nzEyHgfX09PRBuamqqrKdrAD735Q/BrgnTZrk8B5ijR8/3unnkJKLgb8XQ71moA31/+Gt/v+fKpUK\nOp3Oo9cRXTy0Wi1+/vOfo7a2FhUVFXjqqaeQmpqKxMREUc+PiopCUVERysrKIAgC5s+fD41Gg5qa\nGigUCuTn53v0AYj8Se6ZMQNP1K5O3FJiU6lUbmcEOVvfEGzFc6iZTd78PwXrrrj+2Gxx4P9n/x4c\nqUQXjw8++ADvv/8+FixYgP/85z8ArjX/Xn/99SG3LhkoNzcXW7ZscfjZ7bff7vTYRx99VGxoRH4j\n51TYgX/Yu3fvxoMPPjjkiVtqbJ6cjIJxWqmzzzFULsS2miJlV9yB/5/eEF083n//faxbtw5jxoxB\ndXU1gGubJXZ0dHgVAJE7cqwLGOo95PhW2vfe3d3dDn/Y+/btc3niliO2UFlH4iwXcreaQmH9ysD/\nT2+ILh6XLl3C6NGjHX7W29sLpZLXkyL/keME4O49/PmttP97a7VaZGRkoK2tDZmZmViwYIHbE7e/\nvzEHa5dOn4En7P65kLPVZDQaUVhYiPb29qDo3hvKwP9Pb4ieqjtp0iS89957Dj/74IMPRI/ME3lC\njimLDQ0NaGpq8tnuqFLWDPT/fK2trdiwYYNt/n9KSorPFgB6Q8oCODm5m2Lr7bRgKXEUFhbi+PHj\ntt+hYJxa28dX/5+ii8evfvUrfP7551i5ciUuX76Mxx9/HAcPHsQvf/lLrwIgcsXfJwCz2Yz169fb\nLv7j7e6oUtcMDPx8ubm5Dn/YwXriDgbuvlj4avW9uy8DBoMBJ06csN3XaDRB273nS6L7nBITE7Fp\n0ya0trbizJkzSEpKglarRVSU6PpDJJm/u00MBgNaW1sBXCscV69ehdVqHbKbw12/ttSukmDvFvIX\nV3kUO3YgZjzG2249Md2m2dnZmDhxIpqbm6HRaFBVVRUR/4+SBiwUCgW0Wi20Wq2/4iEaxJ/9+v1P\nQBkZGRAEwTbmMPBkJPZEInWAOVJm+vRxlUcpY1zeFl4xRUrMlwE5J1W4K6hyDtpztJsihrM/LGcD\niEOdBILlRBLqXOXRk5abJ4VXbJES+2VArkkVrmKVe3YZ+5woIrgai+g/ruBqjEHs+EskjFN4s5Gg\nqzzKNcgtdiKGL3ct9pTYWOXeD4stD4oIvpi2Gemtir6Wm0ajcbl40Z2+PDY0NAz5mL9zLKV7MdDd\nimJjlXtNjkIIk61ruVjxGu5tZdc/F30tj74/rGCdh+8v3v5e9O8S0Wg0aG9vR29vL4YNG4aqqqoh\nT65D9cEHctuTvlyYzeaQ+SIgNlapnyk1NdXjmNhtRRHBF90PwX5xHmf6Yvb2C0X/llt7ezs0Go2o\nixkN1VUYDFuOh1L3othYxRznq99jdltRxJDS/eBsg8Jg2iBQjP4xZ2dnezWFdGCXyO7du3Hy5EmX\n33BddRWGyrYn4SYgGyMShQuz2Yz6+nooFApMmTLF6TqDgYWi/4mwubkZer0ec+bM8SoGf0+p7B9z\nU1OTV9tzOBuLSElJcfkcVwUi0sePAiUgGyMShQOz2Yx7773X1k2SnZ2N6upqh5OXs2/M2dnZ0Gq1\nMBgM6O3tRUlJCfbs2eNx95ccrZj+J++JEyd6/e1e6sCxuwIR6IHoSBSQjRGJwoHBYEBLS4vtfktL\ny6Bv5M6+McfHx2P9+vV44IEHYLFY0NbWJvmbfF9r49KlS7Js2Nf/5B2ok3SoF4hQ2ClXir7fCb1e\nD6vV6tVrsXhQ0PHnH2xfC6Kv5aHVagd9Ix/qG3NWVhbS09PR3t4uuZ++f2sjIyMDWq0Wra2tg17H\n15+97+TNWXjSheI4l1ilpaVobm7G4sWLPX4NFg8KKv7e2jo+Ph579uxBQ0MDFAoFcnNznb7+wG/M\nZrMZDz74oG2m0e7duyXF1b8rrK2tDbt378bw4cMdilM4n6xCUTBeCMsX+n8ub3CqLgUNKVtbezPd\nMD4+HrfeeivmzJkj+uTc/w/u5MmTOHnypKT3HLhyesqUKYOmVAbD9FWyk2u1u9z6fy5vsOVBQUPs\n1taB+Ibu7dRSMbOLOH01uITrjLD+n8sbXGEeZkK5b3vgKuaqqiqn00E/+eQT28C1qxXOvs6FHCuS\n/fUeofx74WvMhZ03K8xZPMJMqP9huDt5ms1mLFy4EAaDAQAwceLEIafMhnoufIm5sGMu7LwpHuy2\noqDibmpn/6m20dHR2LBhQ9h0JxCFEg6YU0jRaDSIjo4GACiVSl6YjChAWDxIVt5uyta3mysAWCwW\nHD16NOQ2KyQKByweNIi/do91tcuqWH3Xix42bBgyMjLwzDPPePV6ROQZjnmQA39Og/X1BZm6u7vx\nwAMPhN0iLqJQwJYHOfDnQjVfLbrqG1SfMmVKWC7iIgoFbHmQA38uVPP1oqtwXcRFFAq4ziPMDDWH\nXcqGe6F0eU5XOJ/fjrmwYy7suM6DXJI6jhHq22gTkf9xzCMCcMM9ovDlr9mR7rB4RIBw3R2UKNL5\nYvq7p9htFQE4sEwUngJ5zRG2PCJE3ziGp9fc5ipuouATyF4FtjzIpb5mcVNTE9LT04fcJt2T1w2n\na0MTBUIgexXY8iCXDAYDmpqaYLFYcPz4cRQWFnrdAglkPy1RuPGmV8EbsrY89Ho9KioqIAgC5s2b\nh4KCAofHDxw4gOrqagBAXFwcli1bhvHjx8sZIg2QnZ2N9PR0HD9+HMC1jQm97VcN12tDE0US2Voe\nVqsV5eXlWLt2LZ5//nnU1tYOug70mDFjUFpaij/96U8oLCzEyy+/LFd4NIT4+HhUVVVhwoQJUCqV\nyMrK8rpflbO/iEKfbC2PlpYWjBs3DsnJyQCA2bNno66uDmlpabZjsrKybLczMzNhMpnkCo9cSElJ\nwd69e4fsV5U6fsHZX0ShT7aWh8lkQlJSku2+Wq12WRz27duH3NxcOUKjIRiNRrzxxhswGo1D9qt6\nOn4RqH5ab3DWGZFdUM62+uqrr7B//35s2LDB6eONjY1obGy03dfpdFCpVHKFF9RiYmJ8kotTp05h\n1qxZ6OnpQWxsLP73v/9h3Lhxg447fPiww/jFiRMnMH36dK/f3xd8lQsA6OrqQmFhoa2FtXfv3pD6\nnfNlLkIdc+GosrLSdjsnJwc5OTminidb8VCr1ejs7LTdN5lMUKvVg477+uuv8corr2DNmjVDfit1\n9gG50dk1UjZ9c9Xd9N5776GnpwcA0NPTg+rqaixZsmTQa4wfP95hF9709PQh31/u6bm+3ADv0KFD\nMBgM6O3tRVNTE7744ouQGuTnZoB2zIWdSqWCTqfz6LmydVtptVoYjUacOXMGvb29qK2tRV5ensMx\nnZ2deP7557Fq1SqfrCWgobnrbsrPz0dsbCwAIDY2FgsWLHD6On3jF1VVVS43XAz16bkc5CdyJFvL\nIyoqCkVFRSgrK4MgCJg/fz40Gg1qamqgUCiQn5+Pd955B2azGeXl5RAEAdHR0di0aZNcIUYUd9Nl\nU1JS8O9//xv79u3DggULXBZzMbvwhvr0XA7yEzni9TyClKddPGKb5P23addoND5bOe7u/fq6t3x5\neduhsHvCjrmwYy7svLmeB4tHEPL0OuJmsxnffPMNxo8fL+p4o9GIwsJCnDhxAhMnTvT7CV3ui0zx\nJGHHXNgxF3beFA9uTxKEnHXxuJsm2ldw7rrrLtFjCu3t7Whvb4fFYpHlOh+hOD2XiJxj8QhCAwdn\n09LS3A42e3LBJw4CE5GngnKdR6QbODgrZrC5rxD0jSmIKQQcBCYiT3HMIwSIHWw2m804ceIE0tPT\nWQjAvu3+mAs75sKOA+YI7+IBiB9s7v+HEenXzOBJwo65sGMu7DhgHmI82SNJ6mBzqC/KI6LgxuIh\nM7lO6p4MoBMRicXiISOz2Yz33nsPTU1NTk/qvty1lTOpiMifONtKJv0X/imVSigUCoeTuqcLA4fC\nmVRE5E8sHjLp340EAM8++ywWLlwI4NqOrd3d3T7f+0nMnlNERJ5g8ZDJwHUYfYWjr7Wh1WqRkZGB\ntrY2djMRUdBj8ZCJs26kQ4cO2Vobra2teOONN3Ddddexm4mIgh6Lh4wGdiMNbI3k5uayaBBRSGDx\n8DExC/P6H8NBbSIKRSwePiRmxpSzYzioTUShhus8fEjMwjwu3iOicMDi4SNmsxmXLl1CRkaGy4V5\nXLxHROGA3VYe6j9uAThOud29ezemTJnidAyDi/eIKByweHhg4LhFSUmJw5Tb4cOHuywKXLxHRKGO\n3VYeGDhuAUB0V5Qv968iIgoUtjw8MHB9xpQpU0R1Rfl6/yoiokBh8fDAUOMW7rqixFxOlogoFLDb\nykNSL84EcKYVEYUPtjzckHopV1fHc6YVEYULFg8XpI5RiDmeM62IKByw28oFqavBuXqciCIFi4cL\n2dnZ0Gq1UCqVyMjIcDtGwTENIooU7LaC63EKQRBs/9zhmAYRRYqIb3kYjUb85Cc/wX333YdFixY5\nLN4zGAxobW2FxWJBW1ubqG4oT2ZhERGFmoguHmazGYWFhTh+/DgsFguam5sdCgS7oYiInIvobiuD\nwYD29nbbfY1G41Ag2A1FRORcRLc8+loWSqUSEyZMQFVV1ZBTa1k4iIjsIqLlMdSAOFsWRESeCdvi\n0VcwNBoNlixZYtvEsLq6elAB4aI9IiJpZC0eer0eFRUVEAQB8+bNQ0FBwaBjdu7cCb1ej9jYWKxc\nuRITJkyQ9B5msxkNDQ1Yv349WltbkZqaim+++QbAtTEOvV6POXPm+OLjEBFFLNmKh9VqRXl5OUpK\nSpCYmIji4mJMmzYNaWlptmMaGhpw+vRpvPjiizh69CheffVVbNy4UfR79G0P0tTUBIvFAgAOA+IA\nRK3XICIi12QbMG9pacG4ceOQnJwMpVKJ2bNno66uzuGYuro6zJ07FwCQmZmJ7u5unDt3TtTrG41G\nbNmyBQaDwVY4oqOjkZmZiczMTERHR2PixImYMmWKbz8YEVEEkq3lYTKZkJSUZLuvVqvR0tLi9hiT\nyYRRo0a5ff1Zs2ahp6cHCoXCVjSeeeYZW7HgoDgRke+EzYB5T08PgGvdUo888gh+85vfOBQKDooT\nEfmObMVDrVajs7PTdt9kMkGtVg865uzZs7b7Z8+eHXQMADQ2NqKxsdF2X6fTcSyjH5VKFegQggZz\nYcdc2DEXdpWVlbbbOTk5yMnJEfU82cY8tFotjEYjzpw5g97eXtTW1iIvL8/hmLy8PHz88ccAgObm\nZowYMcJpl1VOTg50Op3tX/8PH+mYCzvmwo65sGMu7CorKx3OpWILByBjyyMqKgpFRUUoKyuDIAiY\nP38+NBoNampqoFAokJ+fj1tuuQUNDQ147LHHEBcXhxUrVsgVHhERSSDrmEdubi62bNni8LPbb7/d\n4X5RUZGcIRERkQfCYm8rKU2tcMdc2DEXdsyFHXNh500uFAJHmomISKKwaHkQEZG8WDyIiEiykFok\nKMfGiqHVEv0bAAAFz0lEQVTCXS4OHDiA6upqAEBcXByWLVuG8ePHByJUvxPzewFc2yJn3bp1WL16\nNWbMmCFzlPIQk4vGxkbs2rULFosFCQkJWL9+fQAi9T93ueju7sbWrVvR2dkJq9WKe+65B7fddltg\ngvWjHTt2oL6+HiNHjsTmzZudHuPReVMIERaLRVi1apXw7bffClevXhWeeOIJob293eGY+vp64dln\nnxUEQRCam5uFNWvWBCJUvxOTi6amJuHixYuCIAhCQ0NDROei77jS0lJh06ZNwmeffRaASP1PTC4u\nXrwo/Pa3vxXOnj0rCIIgnD9/PhCh+p2YXPzjH/8Q3nzzTUEQruVh6dKlQm9vbyDC9asjR44Ix44d\nE37/+987fdzT82bIdFv5e2PFUCImF1lZWRg+fDiAa7kwmUyBCNXvxOQCAD788EPMnDkTCQkJAYhS\nHmJyceDAAcyYMcO2c0O45kNMLhQKBS5dugQAuHz5MlQqFaKjowMRrl9lZ2djxIgRQz7u6XkzZIrH\nUJsmSj0mHEj9nPv27UNubq4coclO7O9FXV0d7rjjDrnDk5WYXHR0dMBsNqO0tBTFxcX45JNP5A5T\nFmJyceedd6K9vR3Lly/Hk08+iYceekjmKIODp+fNkCke5JmvvvoK+/fvx5IlSwIdSsBUVFQ4fH4h\ngmenW61WHDt2DMXFxVizZg2qqqpgNBoDHVZA6PV63HDDDXj55Zfx3HPPoby8HJcvXw50WCEjZAbM\nfbmxYqgTkwsA+Prrr/HKK69gzZo1YbsVvZhctLW14YUXXoAgCOjq6kJDQwOUSuWgvdVCndi/EZVK\nhZiYGMTExGDSpEk4fvw4UlJS5A7Xr8TkYv/+/bZB9JSUFIwZMwYnT55ERkaGrLEGmqfnzZBpefhy\nY8VQJyYXnZ2deP7557Fq1aqwOzH0JyYX27Ztw7Zt2/CXv/wFM2fOxMMPPxx2hQMQl4tp06bBYDDA\narWip6cHR48ehUajCVDE/iMmF6NHj8aXX34JADh37hxOnTqFsWPHBiJcvxMEYcgWt6fnzZBaYa7X\n6/H666/bNlYsKChw2FgRAMrLy6HX620bK954440Bjto/3OXipZdewueff47k5GQIgoDo6Ghs2rQp\n0GH7hZjfiz7bt2/H1KlTw3qqrrtc7NmzB/v370dUVBQWLFiAu+66K8BR+4e7XHz33XfYvn07vvvu\nOwBAQUEB5syZE+CofW/Lli04fPgwurq6MHLkSOh0OvT29np93gyp4kFERMEhZLqtiIgoeLB4EBGR\nZCweREQkGYsHERFJxuJBRESSsXgQEZFkLB5ERCQZiweRB1auXImvvvoq0GEQBQyLB5GPWa3WQIdA\n5HdcYU4k0bZt2/Dpp59i2LBhiI6ORmFhId58800sX74c77zzDsaMGQOdToetW7dix44dtuetXLkS\nK1aswA9+8AMIgoDq6mrs27cP3d3d+OEPf4hly5a5vO4CUTBhy4NIolWrVmH06NF46qmnsGvXLvzo\nRz8CABw5cgQvvPAC1q5d6/Y1PvjgAxw6dAgbNmzAyy+/jBEjRuC1117zd+hEPsPiQeQjOp0OMTEx\nGDZsmNtja2pqsHjxYiQmJkKpVOL+++/HZ599xi4vChkhcz0PomDX/2ps7nR2dmLz5s1QKBS2nymV\nSpw/fx6JiYn+CI/Ip1g8iDzQ/6TvTGxsLK5cuWK7b7VaceHCBdv90aNHY8WKFcjKyvJbjET+xG4r\nIg+MGjUKp0+fBuD8srbjxo3DlStX0NDQAIvFgqqqKvT29toez8/Px9tvv2272t2FCxdw6NAheYIn\n8gG2PIg8UFBQgJ07d+KNN97AfffdN+jx4cOH4+GHH8ZLL70Eq9WKhQsXOnRr/fSnPwUAlJWV4bvv\nvsPIkSMxa9assLzCIYUnTtUlIiLJ2G1FRESSsXgQEZFkLB5ERCQZiwcREUnG4kFERJKxeBARkWQs\nHkREJBmLBxERScbiQUREkv0fxCWI4v/e4lAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x17a9d665470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.abs(trues), ests, 'k.')\n",
    "plt.xlabel('true')\n",
    "plt.ylabel('estimated')\n",
    "plt.xlim((0, 1))\n",
    "plt.ylim((0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WIP: Mutliparameter Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also want to try with a multicos model, since that model admits a degeneracy between parameters that is broken by experimental variety. To arrive at a good estimate, then, a resampler must preserve unusual posterior structures that arise in the approach to a unimodal final posterior.\n",
    "\n",
    "It's worth noting here that the traditional LW parameters work quite well for a better experimental protocol (choosing lots of varied experiments), but that this protocol may not always be available. Thus, MLW is a resource here to compensate for experimental restrictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class MultiCosModel(qi.FiniteOutcomeModel):\n",
    "    \n",
    "    def __init__(self, n_terms=2):\n",
    "        self._n_terms = n_terms\n",
    "        super(MultiCosModel, self).__init__()\n",
    "    \n",
    "    @property\n",
    "    def n_modelparams(self):\n",
    "        return self._n_terms\n",
    "    \n",
    "    @property\n",
    "    def is_n_outcomes_constant(self):\n",
    "        return True\n",
    "    def n_outcomes(self, expparams):\n",
    "        return 2\n",
    "    \n",
    "    def are_models_valid(self, modelparams):\n",
    "        return np.all(np.logical_and(modelparams > 0, modelparams <= 1), axis=1)\n",
    "    \n",
    "    @property\n",
    "    def expparams_dtype(self):\n",
    "        return [('ts', '{}float'.format(self._n_terms))]\n",
    "    \n",
    "    def likelihood(self, outcomes, modelparams, expparams):\n",
    "        # We first call the superclass method, which basically\n",
    "        # just makes sure that call count diagnostics are properly\n",
    "        # logged.\n",
    "        super(MultiCosModel, self).likelihood(outcomes, modelparams, expparams)\n",
    "        \n",
    "        # Next, since we have a two-outcome model, everything is defined by\n",
    "        # Pr(0 | modelparams; expparams), so we find the probability of 0\n",
    "        # for each model and each experiment.\n",
    "        #\n",
    "        # We do so by taking a product along the modelparam index (len 2,\n",
    "        # indicating omega_1 or omega_2), then squaring the result.\n",
    "        pr0 = np.prod(\n",
    "            np.cos(\n",
    "                # shape (n_models, 1, 2)\n",
    "                modelparams[:, np.newaxis, :] *\n",
    "                # shape (n_experiments, 2)\n",
    "                expparams['ts']\n",
    "            ), # <- broadcasts to shape (n_models, n_experiments, 2).\n",
    "            axis=2 # <- product over the final index (len 2)\n",
    "        ) ** 2 # square each element\n",
    "        \n",
    "        # Now we use pr0_to_likelihood_array to turn this two index array\n",
    "        # above into the form expected by SMCUpdater and other consumers\n",
    "        # of likelihood().\n",
    "        return qi.FiniteOutcomeModel.pr0_to_likelihood_array(outcomes, pr0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mp_model = qi.BinomialModel(MultiCosModel(3))\n",
    "mp_prior = qi.UniformDistribution([[0, 1]]*3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def multicos_trial(a, h, n_exps=200):\n",
    "    true = mp_prior.sample()\n",
    "    mp_updater = qi.SMCUpdater(mp_model, 2000, mp_prior,\n",
    "        resampler=qi.LiuWestResampler(a=a, h=h, debug=True),\n",
    "        debug_resampling=True\n",
    "    )\n",
    "    \n",
    "    while True:\n",
    "        exp = np.array([(20 * np.pi * np.random.random(3), 40)], dtype=mp_model.expparams_dtype)\n",
    "        # To test the resampler, we take data according to one experiment design until it has to resample, then switch.\n",
    "        # Since the MLW resampler works quite well, we need to set an upper limit.\n",
    "        for idx in range(500):\n",
    "            datum = mp_model.simulate_experiment(true, exp)\n",
    "            mp_updater.update(datum, exp)\n",
    "            \n",
    "            if len(mp_updater.data_record) >= n_exps:\n",
    "                return true, mp_updater.est_mean()\n",
    "            \n",
    "            if mp_updater.just_resampled and idx >= 100:\n",
    "                break\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def nanmean(arr):\n",
    "    return np.nansum(arr) / np.sum(np.isfinite(arr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def mp_risk(a, h, n_trials=10, n_exps=200):\n",
    "    mp_errors = np.empty((n_trials,))\n",
    "    for idx_trial in range(n_trials):\n",
    "        true, est = multicos_trial(a, h, n_exps)\n",
    "        mp_errors[idx_trial] = np.sum((true - est)**2)\n",
    "        \n",
    "    n_nan = np.sum(np.isnan(mp_errors))\n",
    "    if n_nan > 0:\n",
    "        print(\"{} NaNs observed.\".format(n_nan))\n",
    "        \n",
    "    return nanmean(mp_errors)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.2"
  },
  "widgets": {
   "state": {},
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}