model1-grid-approx.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model 1 - grid approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook outlines how to use model 1 as described in the paper Vincent (2015). The paper should drive one's reading, this notebooks is here to supplement one's understanding and to bridge the gap between mathematical description and algorithmic implementation.\n",
    "\n",
    "The generative capability of the model (Figure 6, left) can be used to both generate example data, but also to generate predictions. We also use the model (Figure 6, right) to conduct inference to update our beliefs about our parameter given some observed data.\n",
    "\n",
    "![Model1](img/model1.png)\n",
    "Taken from Figure 6 of Vincent (2015)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Proportion correct (Equation 2 in the paper)\n",
    "From the paper, the proportion of correct responses is given as \n",
    "\n",
    "$PC = \\Phi(\\frac{\\Delta_mu_c}{\\sqrt{2\\sigma^2}})$\n",
    "\n",
    "and we can easily define this as some Julia functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using StatsFuns\n",
    "\n",
    "# define Φ() as the cumulative normal distribution function\n",
    "Φ(x) = normcdf.(0, 1, x)\n",
    "\n",
    "# define the deterministic node PC, corresponds to Equation 2 in Vincent (2015)\n",
    "PC(Δμ, σ²::Number) = Φ(Δμ/sqrt(2σ²));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generative model (Equation 3 in the paper)\n",
    "The number of correct responses $k$ will be Binomially distributed, so the generative model for this task is\n",
    "\n",
    "$k \\sim \\mathrm{Binomial}(PC_c,T)$\n",
    "\n",
    "We can define a function `generative` which will generate random samples from $P(k|T,PC)$. Note that we define two versions of this function depending upon whether the inputs are scalar values or not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define a prior over our parameter (Equation 4 from the paper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "log_prior (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_prior(σ²) = log(1/1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "generative (generic function with 2 methods)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Distributions\n",
    "\n",
    "function generative(T::Int64, Δμ, σ²)\n",
    "    # sample from P(k|T,PC), or equivalently P(k|T,Δμ,σ²)\n",
    "    k = [rand( Binomial(T,p)) for p ∈ PC(Δμ, σ²)]\n",
    "end\n",
    "\n",
    "function generative(T::Int64, Δμ::Number, σ²::Number)\n",
    "    # sample from P(k|T,PC), or equivalently P(k|T,Δμ,σ²)\n",
    "    k = rand(Binomial(T,PC(Δμ, σ²)))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the likelihood and joint distribution (Equation 5 in the paper)\n",
    "First we define the likelihood (part of Equation 5) however we do this as the sum of log probabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "log_likelihood (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function log_likelihood(σ²::Number, Δμ::Array{Float64,1}, kvals::Array{Int64,1}, T::Integer)\n",
    "    log_likelihood = sum(binomlogpdf.(T, PC(Δμ, σ²), kvals))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then simply define the log posterior as the sum of the log prior and log likelihood. This allows us to compute $P(\\sigma^2, \\Delta\\mu, k, T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "log_posterior (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_posterior(σ², Δμ, k, T::Integer) = log_prior(σ²) + log_likelihood(σ², Δμ, k, T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For plotting the posterior, we'll need to convert back from log density, to a normalised posterior density, so below we define this utility function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "logpost2post (generic function with 1 method)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function logpost2post(log_posterior_density::Array{Float64,1})\n",
    "    posterior_density = exp.(log_posterior_density)\n",
    "    return posterior_density ./ sum(posterior_density)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STEP 1: define data\n",
    "If we wanted to generate example data from this observer we can simple sample from $P(k|T,\\Delta\\mu,\\sigma^2)$ using the `generative` function defined above. The example below shows how you would do this for $T=100$ and a set of $\\Delta\\mu$ values, with $\\sigma^2=1$.\n",
    "\n",
    "    generative(100, [0.0100, 0.0215, 0.0464, 0.1000, 0.2154, 0.4642, 1.0000, 2.1544, 4.6416, 10.0000], 1)\n",
    "    \n",
    "If you try runing this code a number of times, you'll see that we get a different set of possible experimental results $k$ each time we sample from $P(k|T,\\Delta\\mu,\\sigma^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10-element Array{Int64,1}:\n",
       "  56\n",
       "  55\n",
       "  50\n",
       "  48\n",
       "  60\n",
       "  59\n",
       "  82\n",
       "  96\n",
       " 100\n",
       " 100"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "generative(100, [0.0100, 0.0215, 0.0464, 0.1000, 0.2154, 0.4642, 1.0000, 2.1544, 4.6416, 10.0000], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But for this example, we will define our own observed data where:\n",
    "- `Δμ` = stimulus intensities\n",
    "- `k` = number of correct responses for each stimulus intensity (could range from 0 to `T`)\n",
    "- `T` = number of trials per stimulus intensity, 100 in this example\n",
    "\n",
    "Proportion correct (per stimulus level) is given by `k./T`. We'll define this data in a Dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = Dict()\n",
    "data[\"T\"] = 100\n",
    "data[\"k\"] = [50, 51, 57, 55, 63, 62, 82, 94, 99, 100]\n",
    "data[\"Δμ\"] = [0.0100, 0.0215, 0.0464, 0.1000, 0.2154, 0.4642, 1.0000, 2.1544, 4.6416, 10.0000];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can now visualise this data in a plot (see Figure 8 in the paper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7200\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7200)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7201\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7200)\" points=\"\n",
       "47.9701,360.065 580.315,360.065 580.315,11.811 47.9701,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7202\">\n",
       "    <rect x=\"47\" y=\"11\" width=\"533\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  107.034,360.065 107.034,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  179.737,360.065 179.737,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.439,360.065 252.439,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  325.142,360.065 325.142,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  397.845,360.065 397.845,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  470.547,360.065 470.547,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  543.25,360.065 543.25,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,350.209 580.315,350.209 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,284.501 580.315,284.501 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,218.792 580.315,218.792 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,153.084 580.315,153.084 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,87.3757 580.315,87.3757 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7202)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,21.6673 580.315,21.6673 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,360.065 47.9701,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  107.034,360.065 107.034,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  179.737,360.065 179.737,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.439,360.065 252.439,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  325.142,360.065 325.142,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  397.845,360.065 397.845,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  470.547,360.065 470.547,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  543.25,360.065 543.25,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,350.209 55.9553,350.209 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,284.501 55.9553,284.501 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,218.792 55.9553,218.792 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,153.084 55.9553,153.084 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,87.3757 55.9553,87.3757 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7200)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,21.6673 55.9553,21.6673 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 107.034, 373.865)\" x=\"107.034\" y=\"373.865\">-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 179.737, 373.865)\" x=\"179.737\" y=\"373.865\">-3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 252.439, 373.865)\" x=\"252.439\" y=\"373.865\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 325.142, 373.865)\" x=\"325.142\" y=\"373.865\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 397.845, 373.865)\" x=\"397.845\" y=\"373.865\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 470.547, 373.865)\" x=\"470.547\" y=\"373.865\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 543.25, 373.865)\" x=\"543.25\" y=\"373.865\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 354.709)\" x=\"41.9701\" y=\"354.709\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 289.001)\" x=\"41.9701\" y=\"289.001\">0.6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 223.292)\" x=\"41.9701\" y=\"223.292\">0.7</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 157.584)\" x=\"41.9701\" y=\"157.584\">0.8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 91.8757)\" x=\"41.9701\" y=\"91.8757\">0.9</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 26.1673)\" x=\"41.9701\" y=\"26.1673\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 314.143, 397.6)\" x=\"314.143\" y=\"397.6\">log signal instensity, ??</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7200)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">proportion correct, k/T</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"63.0365\" cy=\"350.209\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"63.0365\" cy=\"350.209\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"118.688\" cy=\"343.638\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"118.688\" cy=\"343.638\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"174.614\" cy=\"304.213\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"174.614\" cy=\"304.213\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"230.441\" cy=\"317.355\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"230.441\" cy=\"317.355\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"286.227\" cy=\"264.788\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"286.227\" cy=\"264.788\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"342.05\" cy=\"271.359\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"342.05\" cy=\"271.359\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"397.845\" cy=\"139.942\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"397.845\" cy=\"139.942\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"453.645\" cy=\"61.0923\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"453.645\" cy=\"61.0923\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"509.447\" cy=\"28.2381\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"509.447\" cy=\"28.2381\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"565.249\" cy=\"21.6673\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7202)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"565.249\" cy=\"21.6673\" r=\"3\"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "\n",
    "plot_data_space = scatter(log.(data[\"Δμ\"]), data[\"k\"]./data[\"T\"],\n",
    "    xlabel = \"log signal instensity, Δμ\", #xscale = :log,\n",
    "    ylabel = \"proportion correct, k/T\",\n",
    "    legend = false)\n",
    "\n",
    "# TODO: xaxis should be log scale\n",
    "# TODO: get latex axis labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STEP 2: parameter recovery\n",
    "Evaluate posterior over range of σ2 values. We are using grid approximation here, so we will define a fine grid of many σ2 values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5:0.0025025025025025025:3.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "σ² = linspace(0.5, 3, 10.0^3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do grid approximation\n",
    "This is really simple, all we are doing is iterating over a vector of σ2 values, evaluating the log posterior for each"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "function gridApprox(σ²_vector, data)\n",
    "    σ²_log_posterior_density = [log_posterior(σ², data[\"Δμ\"], data[\"k\"], data[\"T\"]) for σ² ∈ σ²_vector]\n",
    "end\n",
    "\n",
    "σ²_log_posterior_density = gridApprox(σ², data)\n",
    "σ²_posterior_density = logpost2post(σ²_log_posterior_density);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot our posterior over the parameter (see Figure 8 in the paper)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7400\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7400)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7401\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7400)\" points=\"\n",
       "61.1494,360.065 580.315,360.065 580.315,11.811 61.1494,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7402\">\n",
       "    <rect x=\"61\" y=\"11\" width=\"520\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  61.1494,360.065 61.1494,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  164.983,360.065 164.983,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  268.816,360.065 268.816,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  372.649,360.065 372.649,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  476.482,360.065 476.482,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  580.315,360.065 580.315,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  61.1494,290.339 580.315,290.339 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  61.1494,220.613 580.315,220.613 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  61.1494,150.886 580.315,150.886 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  61.1494,81.1588 580.315,81.1588 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,360.065 61.1494,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,360.065 61.1494,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  164.983,360.065 164.983,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  268.816,360.065 268.816,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  372.649,360.065 372.649,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  476.482,360.065 476.482,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  580.315,360.065 580.315,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,290.339 68.9369,290.339 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,220.613 68.9369,220.613 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,150.886 68.9369,150.886 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7400)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,81.1588 68.9369,81.1588 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 61.1494, 373.865)\" x=\"61.1494\" y=\"373.865\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 164.983, 373.865)\" x=\"164.983\" y=\"373.865\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 268.816, 373.865)\" x=\"268.816\" y=\"373.865\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 372.649, 373.865)\" x=\"372.649\" y=\"373.865\">2.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 476.482, 373.865)\" x=\"476.482\" y=\"373.865\">2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 580.315, 373.865)\" x=\"580.315\" y=\"373.865\">3.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 55.1494, 294.839)\" x=\"55.1494\" y=\"294.839\">0.001</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 55.1494, 225.113)\" x=\"55.1494\" y=\"225.113\">0.002</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 55.1494, 155.386)\" x=\"55.1494\" y=\"155.386\">0.003</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 55.1494, 85.6588)\" x=\"55.1494\" y=\"85.6588\">0.004</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 320.732, 397.6)\" x=\"320.732\" y=\"397.6\">sigma^2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">posterior density</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7402)\" style=\"stroke:#009af9; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  61.1494,359.994 61.6691,359.984 62.1888,359.973 62.7085,359.96 63.2281,359.946 63.7478,359.929 64.2675,359.912 64.7872,359.892 65.3069,359.869 65.8266,359.845 \n",
       "  66.3462,359.818 66.8659,359.787 67.3856,359.754 67.9053,359.717 68.425,359.677 68.9447,359.632 69.4644,359.583 69.984,359.53 70.5037,359.471 71.0234,359.407 \n",
       "  71.5431,359.337 72.0628,359.261 72.5825,359.178 73.1022,359.088 73.6218,358.99 74.1415,358.884 74.6612,358.769 75.1809,358.645 75.7006,358.511 76.2203,358.366 \n",
       "  76.74,358.211 77.2596,358.043 77.7793,357.863 78.299,357.67 78.8187,357.463 79.3384,357.241 79.8581,357.004 80.3778,356.75 80.8974,356.48 81.4171,356.192 \n",
       "  81.9368,355.885 82.4565,355.559 82.9762,355.212 83.4959,354.844 84.0155,354.453 84.5352,354.04 85.0549,353.602 85.5746,353.14 86.0943,352.651 86.614,352.136 \n",
       "  87.1337,351.592 87.6533,351.02 88.173,350.418 88.6927,349.785 89.2124,349.121 89.7321,348.423 90.2518,347.692 90.7715,346.927 91.2911,346.125 91.8108,345.288 \n",
       "  92.3305,344.412 92.8502,343.498 93.3699,342.545 93.8896,341.552 94.4093,340.518 94.9289,339.442 95.4486,338.323 95.9683,337.161 96.488,335.955 97.0077,334.703 \n",
       "  97.5274,333.406 98.047,332.063 98.5667,330.673 99.0864,329.235 99.6061,327.75 100.126,326.215 100.645,324.632 101.165,322.999 101.685,321.317 102.205,319.584 \n",
       "  102.724,317.801 103.244,315.967 103.764,314.083 104.283,312.147 104.803,310.16 105.323,308.123 105.842,306.034 106.362,303.894 106.882,301.704 107.401,299.463 \n",
       "  107.921,297.171 108.441,294.83 108.96,292.439 109.48,289.998 110,287.509 110.519,284.972 111.039,282.387 111.559,279.755 112.079,277.076 112.598,274.352 \n",
       "  113.118,271.584 113.638,268.771 114.157,265.916 114.677,263.018 115.197,260.08 115.716,257.102 116.236,254.085 116.756,251.03 117.275,247.939 117.795,244.813 \n",
       "  118.315,241.653 118.834,238.46 119.354,235.237 119.874,231.984 120.394,228.702 120.913,225.394 121.433,222.061 121.953,218.704 122.472,215.325 122.992,211.925 \n",
       "  123.512,208.507 124.031,205.072 124.551,201.621 125.071,198.156 125.59,194.68 126.11,191.193 126.63,187.698 127.149,184.196 127.669,180.689 128.189,177.179 \n",
       "  128.708,173.668 129.228,170.157 129.748,166.649 130.268,163.144 130.787,159.645 131.307,156.154 131.827,152.673 132.346,149.202 132.866,145.745 133.386,142.303 \n",
       "  133.905,138.877 134.425,135.469 134.945,132.081 135.464,128.715 135.984,125.372 136.504,122.055 137.023,118.764 137.543,115.501 138.063,112.268 138.582,109.067 \n",
       "  139.102,105.898 139.622,102.764 140.142,99.6666 140.661,96.6061 141.181,93.5844 141.701,90.603 142.22,87.6631 142.74,84.7661 143.26,81.9132 143.779,79.1057 \n",
       "  144.299,76.3447 144.819,73.6316 145.338,70.9672 145.858,68.3528 146.378,65.7894 146.897,63.278 147.417,60.8195 147.937,58.4149 148.457,56.065 148.976,53.7707 \n",
       "  149.496,51.5328 150.016,49.352 150.535,47.229 151.055,45.1645 151.575,43.1591 152.094,41.2134 152.614,39.328 153.134,37.5033 153.653,35.7397 154.173,34.0378 \n",
       "  154.693,32.3979 155.212,30.8202 155.732,29.3052 156.252,27.853 156.771,26.4638 157.291,25.1379 157.811,23.8753 158.331,22.6762 158.85,21.5406 159.37,20.4685 \n",
       "  159.89,19.4599 160.409,18.5147 160.929,17.633 161.449,16.8144 161.968,16.0589 162.488,15.3663 163.008,14.7364 163.527,14.1688 164.047,13.6634 164.567,13.2197 \n",
       "  165.086,12.8375 165.606,12.5163 166.126,12.2558 166.646,12.0554 167.165,11.9149 167.685,11.8336 168.205,11.811 168.724,11.8467 169.244,11.94 169.764,12.0904 \n",
       "  170.283,12.2973 170.803,12.5601 171.323,12.878 171.842,13.2505 172.362,13.6769 172.882,14.1564 173.401,14.6884 173.921,15.2721 174.441,15.9068 174.96,16.5917 \n",
       "  175.48,17.326 176,18.1091 176.52,18.9399 177.039,19.8179 177.559,20.742 178.079,21.7116 178.598,22.7257 179.118,23.7836 179.638,24.8843 180.157,26.0271 \n",
       "  180.677,27.211 181.197,28.4352 181.716,29.6988 182.236,31.001 182.756,32.3408 183.275,33.7173 183.795,35.1298 184.315,36.5772 184.834,38.0588 185.354,39.5736 \n",
       "  185.874,41.1207 186.394,42.6992 186.913,44.3084 187.433,45.9472 187.953,47.6148 188.472,49.3103 188.992,51.0329 189.512,52.7816 190.031,54.5557 190.551,56.3542 \n",
       "  191.071,58.1763 191.59,60.0211 192.11,61.8879 192.63,63.7757 193.149,65.6837 193.669,67.6111 194.189,69.5572 194.709,71.521 195.228,73.5017 195.748,75.4987 \n",
       "  196.268,77.511 196.787,79.538 197.307,81.5788 197.827,83.6327 198.346,85.6989 198.866,87.7767 199.386,89.8654 199.905,91.9643 200.425,94.0726 200.945,96.1896 \n",
       "  201.464,98.3147 201.984,100.447 202.504,102.586 203.023,104.731 203.543,106.882 204.063,109.037 204.583,111.197 205.102,113.36 205.622,115.526 206.142,117.694 \n",
       "  206.661,119.864 207.181,122.035 207.701,124.207 208.22,126.379 208.74,128.55 209.26,130.721 209.779,132.89 210.299,135.056 210.819,137.221 211.338,139.382 \n",
       "  211.858,141.54 212.378,143.694 212.897,145.844 213.417,147.989 213.937,150.129 214.457,152.263 214.976,154.391 215.496,156.513 216.016,158.628 216.535,160.736 \n",
       "  217.055,162.836 217.575,164.929 218.094,167.013 218.614,169.088 219.134,171.155 219.653,173.213 220.173,175.261 220.693,177.3 221.212,179.328 221.732,181.346 \n",
       "  222.252,183.354 222.772,185.351 223.291,187.337 223.811,189.311 224.331,191.274 224.85,193.225 225.37,195.165 225.89,197.092 226.409,199.007 226.929,200.909 \n",
       "  227.449,202.799 227.968,204.676 228.488,206.54 229.008,208.391 229.527,210.229 230.047,212.053 230.567,213.863 231.086,215.66 231.606,217.444 232.126,219.213 \n",
       "  232.646,220.969 233.165,222.71 233.685,224.437 234.205,226.15 234.724,227.849 235.244,229.534 235.764,231.204 236.283,232.859 236.803,234.5 237.323,236.127 \n",
       "  237.842,237.738 238.362,239.336 238.882,240.918 239.401,242.486 239.921,244.04 240.441,245.578 240.96,247.102 241.48,248.611 242,250.105 242.52,251.585 \n",
       "  243.039,253.05 243.559,254.501 244.079,255.936 244.598,257.357 245.118,258.764 245.638,260.156 246.157,261.533 246.677,262.896 247.197,264.244 247.716,265.578 \n",
       "  248.236,266.897 248.756,268.203 249.275,269.493 249.795,270.77 250.315,272.033 250.835,273.281 251.354,274.515 251.874,275.736 252.394,276.942 252.913,278.135 \n",
       "  253.433,279.314 253.953,280.479 254.472,281.631 254.992,282.769 255.512,283.893 256.031,285.004 256.551,286.102 257.071,287.187 257.59,288.259 258.11,289.317 \n",
       "  258.63,290.363 259.149,291.396 259.669,292.416 260.189,293.423 260.709,294.418 261.228,295.4 261.748,296.37 262.268,297.327 262.787,298.273 263.307,299.206 \n",
       "  263.827,300.127 264.346,301.037 264.866,301.935 265.386,302.821 265.905,303.695 266.425,304.558 266.945,305.41 267.464,306.25 267.984,307.079 268.504,307.897 \n",
       "  269.023,308.704 269.543,309.5 270.063,310.286 270.583,311.06 271.102,311.825 271.622,312.579 272.142,313.322 272.661,314.055 273.181,314.778 273.701,315.492 \n",
       "  274.22,316.195 274.74,316.888 275.26,317.572 275.779,318.246 276.299,318.91 276.819,319.566 277.338,320.211 277.858,320.848 278.378,321.476 278.898,322.094 \n",
       "  279.417,322.704 279.937,323.305 280.457,323.897 280.976,324.481 281.496,325.056 282.016,325.623 282.535,326.181 283.055,326.731 283.575,327.274 284.094,327.808 \n",
       "  284.614,328.334 285.134,328.853 285.653,329.363 286.173,329.867 286.693,330.362 287.212,330.85 287.732,331.331 288.252,331.805 288.772,332.272 289.291,332.731 \n",
       "  289.811,333.183 290.331,333.629 290.85,334.068 291.37,334.5 291.89,334.926 292.409,335.344 292.929,335.757 293.449,336.163 293.968,336.563 294.488,336.957 \n",
       "  295.008,337.344 295.527,337.726 296.047,338.102 296.567,338.471 297.086,338.835 297.606,339.193 298.126,339.546 298.646,339.893 299.165,340.235 299.685,340.571 \n",
       "  300.205,340.902 300.724,341.228 301.244,341.548 301.764,341.863 302.283,342.174 302.803,342.479 303.323,342.78 303.842,343.075 304.362,343.366 304.882,343.652 \n",
       "  305.401,343.934 305.921,344.211 306.441,344.484 306.961,344.752 307.48,345.016 308,345.275 308.52,345.53 309.039,345.781 309.559,346.028 310.079,346.271 \n",
       "  310.598,346.51 311.118,346.746 311.638,346.977 312.157,347.204 312.677,347.428 313.197,347.648 313.716,347.864 314.236,348.077 314.756,348.286 315.275,348.492 \n",
       "  315.795,348.694 316.315,348.893 316.835,349.089 317.354,349.281 317.874,349.47 318.394,349.656 318.913,349.839 319.433,350.019 319.953,350.196 320.472,350.369 \n",
       "  320.992,350.54 321.512,350.708 322.031,350.873 322.551,351.036 323.071,351.195 323.59,351.352 324.11,351.506 324.63,351.658 325.15,351.807 325.669,351.954 \n",
       "  326.189,352.098 326.709,352.239 327.228,352.378 327.748,352.515 328.268,352.649 328.787,352.781 329.307,352.911 329.827,353.039 330.346,353.164 330.866,353.288 \n",
       "  331.386,353.409 331.905,353.528 332.425,353.645 332.945,353.76 333.464,353.872 333.984,353.983 334.504,354.093 335.024,354.2 335.543,354.305 336.063,354.409 \n",
       "  336.583,354.51 337.102,354.61 337.622,354.708 338.142,354.805 338.661,354.899 339.181,354.993 339.701,355.084 340.22,355.174 340.74,355.262 341.26,355.349 \n",
       "  341.779,355.434 342.299,355.518 342.819,355.6 343.338,355.681 343.858,355.76 344.378,355.838 344.898,355.915 345.417,355.99 345.937,356.064 346.457,356.137 \n",
       "  346.976,356.208 347.496,356.278 348.016,356.347 348.535,356.415 349.055,356.481 349.575,356.546 350.094,356.61 350.614,356.673 351.134,356.735 351.653,356.796 \n",
       "  352.173,356.855 352.693,356.914 353.213,356.971 353.732,357.028 354.252,357.083 354.772,357.138 355.291,357.191 355.811,357.244 356.331,357.296 356.85,357.346 \n",
       "  357.37,357.396 357.89,357.445 358.409,357.493 358.929,357.54 359.449,357.586 359.968,357.632 360.488,357.676 361.008,357.72 361.527,357.763 362.047,357.806 \n",
       "  362.567,357.847 363.087,357.888 363.606,357.928 364.126,357.967 364.646,358.006 365.165,358.044 365.685,358.081 366.205,358.117 366.724,358.153 367.244,358.189 \n",
       "  367.764,358.223 368.283,358.257 368.803,358.29 369.323,358.323 369.842,358.355 370.362,358.387 370.882,358.418 371.401,358.448 371.921,358.478 372.441,358.507 \n",
       "  372.961,358.536 373.48,358.564 374,358.592 374.52,358.619 375.039,358.646 375.559,358.672 376.079,358.698 376.598,358.723 377.118,358.748 377.638,358.773 \n",
       "  378.157,358.797 378.677,358.82 379.197,358.843 379.716,358.866 380.236,358.888 380.756,358.91 381.276,358.931 381.795,358.952 382.315,358.973 382.835,358.993 \n",
       "  383.354,359.013 383.874,359.033 384.394,359.052 384.913,359.071 385.433,359.089 385.953,359.107 386.472,359.125 386.992,359.142 387.512,359.16 388.031,359.176 \n",
       "  388.551,359.193 389.071,359.209 389.59,359.225 390.11,359.241 390.63,359.256 391.15,359.271 391.669,359.286 392.189,359.3 392.709,359.315 393.228,359.329 \n",
       "  393.748,359.342 394.268,359.356 394.787,359.369 395.307,359.382 395.827,359.395 396.346,359.407 396.866,359.419 397.386,359.431 397.905,359.443 398.425,359.455 \n",
       "  398.945,359.466 399.464,359.477 399.984,359.488 400.504,359.499 401.024,359.51 401.543,359.52 402.063,359.53 402.583,359.54 403.102,359.55 403.622,359.56 \n",
       "  404.142,359.569 404.661,359.578 405.181,359.587 405.701,359.596 406.22,359.605 406.74,359.614 407.26,359.622 407.779,359.63 408.299,359.638 408.819,359.646 \n",
       "  409.339,359.654 409.858,359.662 410.378,359.669 410.898,359.677 411.417,359.684 411.937,359.691 412.457,359.698 412.976,359.705 413.496,359.712 414.016,359.718 \n",
       "  414.535,359.725 415.055,359.731 415.575,359.738 416.094,359.744 416.614,359.75 417.134,359.756 417.653,359.761 418.173,359.767 418.693,359.773 419.213,359.778 \n",
       "  419.732,359.783 420.252,359.789 420.772,359.794 421.291,359.799 421.811,359.804 422.331,359.809 422.85,359.814 423.37,359.818 423.89,359.823 424.409,359.828 \n",
       "  424.929,359.832 425.449,359.836 425.968,359.841 426.488,359.845 427.008,359.849 427.527,359.853 428.047,359.857 428.567,359.861 429.087,359.865 429.606,359.868 \n",
       "  430.126,359.872 430.646,359.876 431.165,359.879 431.685,359.883 432.205,359.886 432.724,359.89 433.244,359.893 433.764,359.896 434.283,359.899 434.803,359.902 \n",
       "  435.323,359.905 435.842,359.908 436.362,359.911 436.882,359.914 437.402,359.917 437.921,359.92 438.441,359.923 438.961,359.925 439.48,359.928 440,359.93 \n",
       "  440.52,359.933 441.039,359.935 441.559,359.938 442.079,359.94 442.598,359.943 443.118,359.945 443.638,359.947 444.157,359.949 444.677,359.952 445.197,359.954 \n",
       "  445.716,359.956 446.236,359.958 446.756,359.96 447.276,359.962 447.795,359.964 448.315,359.966 448.835,359.968 449.354,359.969 449.874,359.971 450.394,359.973 \n",
       "  450.913,359.975 451.433,359.976 451.953,359.978 452.472,359.98 452.992,359.981 453.512,359.983 454.031,359.984 454.551,359.986 455.071,359.987 455.591,359.989 \n",
       "  456.11,359.99 456.63,359.992 457.15,359.993 457.669,359.995 458.189,359.996 458.709,359.997 459.228,359.998 459.748,360 460.268,360.001 460.787,360.002 \n",
       "  461.307,360.003 461.827,360.005 462.346,360.006 462.866,360.007 463.386,360.008 463.905,360.009 464.425,360.01 464.945,360.011 465.465,360.012 465.984,360.013 \n",
       "  466.504,360.014 467.024,360.015 467.543,360.016 468.063,360.017 468.583,360.018 469.102,360.019 469.622,360.02 470.142,360.02 470.661,360.021 471.181,360.022 \n",
       "  471.701,360.023 472.22,360.024 472.74,360.025 473.26,360.025 473.779,360.026 474.299,360.027 474.819,360.028 475.339,360.028 475.858,360.029 476.378,360.03 \n",
       "  476.898,360.03 477.417,360.031 477.937,360.032 478.457,360.032 478.976,360.033 479.496,360.034 480.016,360.034 480.535,360.035 481.055,360.035 481.575,360.036 \n",
       "  482.094,360.036 482.614,360.037 483.134,360.038 483.654,360.038 484.173,360.039 484.693,360.039 485.213,360.04 485.732,360.04 486.252,360.041 486.772,360.041 \n",
       "  487.291,360.042 487.811,360.042 488.331,360.042 488.85,360.043 489.37,360.043 489.89,360.044 490.409,360.044 490.929,360.045 491.449,360.045 491.968,360.045 \n",
       "  492.488,360.046 493.008,360.046 493.528,360.046 494.047,360.047 494.567,360.047 495.087,360.048 495.606,360.048 496.126,360.048 496.646,360.049 497.165,360.049 \n",
       "  497.685,360.049 498.205,360.05 498.724,360.05 499.244,360.05 499.764,360.05 500.283,360.051 500.803,360.051 501.323,360.051 501.842,360.052 502.362,360.052 \n",
       "  502.882,360.052 503.402,360.052 503.921,360.053 504.441,360.053 504.961,360.053 505.48,360.053 506,360.054 506.52,360.054 507.039,360.054 507.559,360.054 \n",
       "  508.079,360.055 508.598,360.055 509.118,360.055 509.638,360.055 510.157,360.055 510.677,360.056 511.197,360.056 511.717,360.056 512.236,360.056 512.756,360.056 \n",
       "  513.276,360.057 513.795,360.057 514.315,360.057 514.835,360.057 515.354,360.057 515.874,360.057 516.394,360.058 516.913,360.058 517.433,360.058 517.953,360.058 \n",
       "  518.472,360.058 518.992,360.058 519.512,360.058 520.031,360.059 520.551,360.059 521.071,360.059 521.591,360.059 522.11,360.059 522.63,360.059 523.15,360.059 \n",
       "  523.669,360.06 524.189,360.06 524.709,360.06 525.228,360.06 525.748,360.06 526.268,360.06 526.787,360.06 527.307,360.06 527.827,360.06 528.346,360.061 \n",
       "  528.866,360.061 529.386,360.061 529.905,360.061 530.425,360.061 530.945,360.061 531.465,360.061 531.984,360.061 532.504,360.061 533.024,360.061 533.543,360.062 \n",
       "  534.063,360.062 534.583,360.062 535.102,360.062 535.622,360.062 536.142,360.062 536.661,360.062 537.181,360.062 537.701,360.062 538.22,360.062 538.74,360.062 \n",
       "  539.26,360.062 539.78,360.062 540.299,360.063 540.819,360.063 541.339,360.063 541.858,360.063 542.378,360.063 542.898,360.063 543.417,360.063 543.937,360.063 \n",
       "  544.457,360.063 544.976,360.063 545.496,360.063 546.016,360.063 546.535,360.063 547.055,360.063 547.575,360.063 548.094,360.063 548.614,360.064 549.134,360.064 \n",
       "  549.654,360.064 550.173,360.064 550.693,360.064 551.213,360.064 551.732,360.064 552.252,360.064 552.772,360.064 553.291,360.064 553.811,360.064 554.331,360.064 \n",
       "  554.85,360.064 555.37,360.064 555.89,360.064 556.409,360.064 556.929,360.064 557.449,360.064 557.968,360.064 558.488,360.064 559.008,360.064 559.528,360.064 \n",
       "  560.047,360.064 560.567,360.064 561.087,360.065 561.606,360.065 562.126,360.065 562.646,360.065 563.165,360.065 563.685,360.065 564.205,360.065 564.724,360.065 \n",
       "  565.244,360.065 565.764,360.065 566.283,360.065 566.803,360.065 567.323,360.065 567.843,360.065 568.362,360.065 568.882,360.065 569.402,360.065 569.921,360.065 \n",
       "  570.441,360.065 570.961,360.065 571.48,360.065 572,360.065 572.52,360.065 573.039,360.065 573.559,360.065 574.079,360.065 574.598,360.065 575.118,360.065 \n",
       "  575.638,360.065 576.157,360.065 576.677,360.065 577.197,360.065 577.717,360.065 578.236,360.065 578.756,360.065 579.276,360.065 579.795,360.065 580.315,360.065 \n",
       "  \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_parameter_space = plot(σ², σ²_posterior_density,\n",
    "    xlabel = \"sigma^2\",\n",
    "    ylabel = \"posterior density\",\n",
    "    legend=false)\n",
    "# TODO: get latex axis labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## STEP 3: Predictive distribution\n",
    "Now we have our posterior over parameters we can use the generative ability of the model to sample predicted behaviour $k$ given our knowledge of $\\sigma^2$. Below, we generate a set of predictions for many interpolated signal intensity levels (`Δμ_interp_vals`) for which we do not have data for. This is a useful demonstration of making and visualising the model's predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "Δμ_interp_vals = logspace(-2, 1, 20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw many samples from the posterior distribution of internal variance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100-element Array{Float64,1}:\n",
       " 0.735235\n",
       " 1.06306 \n",
       " 1.45095 \n",
       " 1.40841 \n",
       " 0.882883\n",
       " 0.917918\n",
       " 1.02052 \n",
       " 1.08058 \n",
       " 1.32332 \n",
       " 1.48849 \n",
       " 1.23323 \n",
       " 1.07808 \n",
       " 0.94044 \n",
       " ⋮       \n",
       " 1.32082 \n",
       " 1.1957  \n",
       " 1.21572 \n",
       " 1.37838 \n",
       " 0.945445\n",
       " 0.867868\n",
       " 1.06557 \n",
       " 0.79029 \n",
       " 0.755255\n",
       " 1.003   \n",
       " 1.40591 \n",
       " 1.1031  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using StatsBase # for sampling\n",
    "nSamples = 10^2\n",
    "σ²_samples = sample(σ², ProbabilityWeights(σ²_posterior_density), nSamples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will generate many predicted data points $k$ over each interpolated Δμ values, with σ2 sampled from our posterior "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20×100 Array{Int64,2}:\n",
       "  49   46   53   49   36   57   56   45  …   42   45   42   45   44   49   53\n",
       "  58   47   53   45   52   43   49   55      52   57   55   54   49   52   56\n",
       "  56   46   47   61   40   53   48   51      65   52   42   37   59   50   48\n",
       "  53   49   50   56   54   48   45   53      62   57   55   57   49   53   49\n",
       "  52   46   43   48   48   48   50   53      50   53   51   62   48   49   53\n",
       "  50   55   62   48   48   44   45   52  …   47   54   54   53   44   52   51\n",
       "  44   50   50   55   47   54   60   48      49   60   52   48   44   53   49\n",
       "  57   55   47   62   59   61   47   52      50   56   56   53   49   55   51\n",
       "  54   56   61   50   53   51   56   52      58   57   51   48   47   61   52\n",
       "  57   66   59   56   52   57   62   63      49   54   52   57   53   58   58\n",
       "  66   65   54   61   64   65   64   66  …   49   57   66   63   58   58   62\n",
       "  63   67   66   60   66   64   69   64      65   72   61   71   69   63   62\n",
       "  80   72   70   70   77   77   80   66      74   70   69   70   70   63   64\n",
       "  87   76   78   71   78   83   80   84      80   80   77   82   82   71   67\n",
       "  95   92   80   91   84   91   89   86      89   87   94   91   92   78   86\n",
       "  98   98   90   94   95   96   94   94  …   95   94   98   98  100   90   97\n",
       " 100  100   98  100   99   99  100   97      99   99  100  100   99   98   97\n",
       " 100  100  100  100  100  100  100  100     100  100  100  100  100  100  100\n",
       " 100  100  100  100  100  100  100  100     100  100  100  100  100  100  100\n",
       " 100  100  100  100  100  100  100  100     100  100  100  100  100  100  100"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# do the posterior prediction\n",
    "predk = [generative(data[\"T\"], Δμ, σ²) for Δμ ∈ Δμ_interp_vals, σ² ∈ σ²_samples]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ---- CODE BELOW IS IN PROGRESS ----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # Calculate 95% CI's for each signal level\n",
    "# #CI = prctile(predk,[5 95]) ./ data.T;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# save('output/m1MAPestimate.mat', 'vMode')\n",
    "# save(['~/Dropbox/tempModelOutputs/tempModel1run_gridApprox.mat'], '-v7.3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 600 400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7600\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"600\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7600)\" points=\"\n",
       "0,400 600,400 600,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7601\">\n",
       "    <rect x=\"120\" y=\"0\" width=\"421\" height=\"400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7600)\" points=\"\n",
       "47.9701,360.065 580.315,360.065 580.315,11.811 47.9701,11.811 \n",
       "  \" fill=\"#ffffff\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7602\">\n",
       "    <rect x=\"47\" y=\"11\" width=\"533\" height=\"349\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  107.034,360.065 107.034,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  179.737,360.065 179.737,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  252.439,360.065 252.439,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  325.142,360.065 325.142,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  397.845,360.065 397.845,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  470.547,360.065 470.547,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  543.25,360.065 543.25,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,350.209 580.315,350.209 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,284.501 580.315,284.501 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,218.792 580.315,218.792 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,153.084 580.315,153.084 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,87.3757 580.315,87.3757 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7602)\" style=\"stroke:#000000; stroke-width:0.5; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  47.9701,21.6673 580.315,21.6673 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,360.065 580.315,360.065 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,360.065 47.9701,11.811 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  107.034,360.065 107.034,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  179.737,360.065 179.737,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  252.439,360.065 252.439,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  325.142,360.065 325.142,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  397.845,360.065 397.845,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  470.547,360.065 470.547,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  543.25,360.065 543.25,354.842 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,350.209 55.9553,350.209 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,284.501 55.9553,284.501 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,218.792 55.9553,218.792 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,153.084 55.9553,153.084 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,87.3757 55.9553,87.3757 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7600)\" style=\"stroke:#000000; stroke-width:1; stroke-opacity:1; fill:none\" points=\"\n",
       "  47.9701,21.6673 55.9553,21.6673 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 107.034, 373.865)\" x=\"107.034\" y=\"373.865\">-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 179.737, 373.865)\" x=\"179.737\" y=\"373.865\">-3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 252.439, 373.865)\" x=\"252.439\" y=\"373.865\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 325.142, 373.865)\" x=\"325.142\" y=\"373.865\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 397.845, 373.865)\" x=\"397.845\" y=\"373.865\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 470.547, 373.865)\" x=\"470.547\" y=\"373.865\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:middle;\" transform=\"rotate(0, 543.25, 373.865)\" x=\"543.25\" y=\"373.865\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 354.709)\" x=\"41.9701\" y=\"354.709\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 289.001)\" x=\"41.9701\" y=\"289.001\">0.6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 223.292)\" x=\"41.9701\" y=\"223.292\">0.7</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 157.584)\" x=\"41.9701\" y=\"157.584\">0.8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 91.8757)\" x=\"41.9701\" y=\"91.8757\">0.9</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:12; text-anchor:end;\" transform=\"rotate(0, 41.9701, 26.1673)\" x=\"41.9701\" y=\"26.1673\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(0, 314.143, 397.6)\" x=\"314.143\" y=\"397.6\">log signal instensity, ??</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:16; text-anchor:middle;\" transform=\"rotate(-90, 14.4, 185.938)\" x=\"14.4\" y=\"185.938\">proportion correct, k/T</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"63.0365\" cy=\"350.209\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"63.0365\" cy=\"350.209\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"118.688\" cy=\"343.638\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"118.688\" cy=\"343.638\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"174.614\" cy=\"304.213\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"174.614\" cy=\"304.213\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"230.441\" cy=\"317.355\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"230.441\" cy=\"317.355\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"286.227\" cy=\"264.788\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"286.227\" cy=\"264.788\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"342.05\" cy=\"271.359\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"342.05\" cy=\"271.359\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"397.845\" cy=\"139.942\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"397.845\" cy=\"139.942\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"453.645\" cy=\"61.0923\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"453.645\" cy=\"61.0923\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"509.447\" cy=\"28.2381\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"509.447\" cy=\"28.2381\" r=\"3\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"565.249\" cy=\"21.6673\" r=\"4\"/>\n",
       "<circle clip-path=\"url(#clip7602)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"565.249\" cy=\"21.6673\" r=\"3\"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(plot_data_space)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Overall results should look something like this\n",
    "![Model1](img/model1results.png)\n",
    "Taken from Figure 8 of Vincent (2015)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "Vincent, B. T. (2015). A tutorial on Bayesian models of perception. Journal of Mathematical Psychology, 66, 103–114. http://doi.org/10.1016/j.jmp.2015.02.001"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}