mshvartsman/ddm.ipynb

## ddm.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Python version\n",
    "\n",
    "import numpy as np\n",
    "from numpy import sqrt, random, abs\n",
    "\n",
    "def ou_update(particle, a, s, dt, l):\n",
    "    ''' Single update for OU (special case l=0 is DDM)'''\n",
    "    return particle + dt * (l * particle + a) + random.normal(0, s)*sqrt(dt)\n",
    "\n",
    "def ddm_update(particle, a, s, dt):\n",
    "    return ou_update(particle, a, s, dt, l=0)\n",
    "\n",
    "def ddm_rt(x0, t0, a, s, z, dt):\n",
    "    samps = 0\n",
    "    particle = x0\n",
    "    while abs(particle) < z:\n",
    "        samps = samps + 1\n",
    "        particle = ddm_update(particle, a, s, dt)\n",
    "    # return -rt for errors as per convention\n",
    "    return (samps * dt + t0) if particle > 0 else -(samps * dt + t0)\n",
    "\n",
    "def ddm_distr(n, x0, t0, a, s, z, dt):\n",
    "    return np.fromiter((ddm_rt(x0, t0, a, s, z, dt) for i in range(n)), dtype='float64')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Cython version, ddm only\n",
    "import cython\n",
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%cython --cplus --compile-args=-O3 --compile-args=-ffast-math\n",
    "\n",
    "import numpy as np\n",
    "from numpy import random\n",
    "from libcpp.vector cimport vector\n",
    "from libc.math cimport exp, log, sqrt, fabs\n",
    "\n",
    "cdef vector[float] ddm_rts_c(int n, float t0, float x0, float a, float s, float z, float dt):\n",
    "    ''' Single RT for DDM'''\n",
    "    cdef int samps, i\n",
    "    cdef float particle\n",
    "    cdef float atilde = a * dt\n",
    "    cdef float stilde = s * sqrt(dt)\n",
    "    cdef vector[float] out = vector[float](n)\n",
    "    for i in range(n):\n",
    "        samps = 0\n",
    "        particle = x0\n",
    "        while fabs(particle) < z:\n",
    "            samps = samps + 1\n",
    "            particle = particle + atilde + random.normal(0, stilde)\n",
    "        # return -rt for errors as per convention\n",
    "        out[i] = (t0 + dt*samps) if particle > 0 else -(t0 + dt*samps)\n",
    "    return out\n",
    "\n",
    "def ddm_distr_cython(n, x0, t0, a, s, z, dt):\n",
    "    return np.array(ddm_rts_c(n, t0, x0, a, s, z, dt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 18.7 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "sim = ddm_distr(n=1000, x0=0, t0=0.3, a=0.5, s=1, z=3, dt = 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 1.26 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "sim = ddm_distr_cython(n=1000, x0=0, t0=0.3, a=0.5, s=1, z=3, dt = 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Analytic RT, ER: 6.221349, 0.064888; simulated RT, ER: 6.270605, 0.060700.\n",
      "Zero drift: Analytic RT, ER: 9.310000, 0.550000; simulated RT, ER: 9.432633, 0.546000.\n"
     ]
    }
   ],
   "source": [
    "from numpy import tanh, exp\n",
    "\n",
    "# Compare to analytic\n",
    "def ddm_analytic(x0, t0, a, s, z):\n",
    "    if abs(a) < 1e-8: # a close to 0 (avoid float comparison)\n",
    "        # use expression for limit a->0 from Srivastava et al. 2016\n",
    "        rt = t0 + (z**2 - x0**2)/(s**2)\n",
    "        er = (z - x0)/(2*z)\n",
    "    # expression from Bogacz et al. 2006 for nonzero drift\n",
    "    else:\n",
    "        ztilde = z/a\n",
    "        atilde = (a/s)**2\n",
    "        x0tilde = x0/a\n",
    "        rt = ztilde * tanh(ztilde * atilde) + ((2*ztilde*(1-exp(-2*x0tilde*atilde)))/(exp(2*ztilde*atilde)-exp(-2*ztilde*atilde))-x0tilde) + t0\n",
    "        er = 1/(1+exp(2*ztilde*atilde)) - ((1-exp(-2*x0tilde*atilde))/(exp(2*ztilde*atilde)-exp(-2*ztilde*atilde)))\n",
    "    return rt, er\n",
    "\n",
    "n = 10000\n",
    "x0 = -0.3\n",
    "t0 = 0.4\n",
    "a = 0.5\n",
    "s = 1\n",
    "z = 3\n",
    "dt = 0.001\n",
    "\n",
    "sim = ddm_distr_cython(n, x0, t0, a, s, z, dt)\n",
    "\n",
    "analytic_rt, analytic_er = ddm_analytic(x0, t0, a, s, z)\n",
    "\n",
    "simulated_rt = np.mean(abs(sim))\n",
    "simulated_er = np.sum(sim<0) / n\n",
    "\n",
    "print(\"Analytic RT, ER: %f, %f; simulated RT, ER: %f, %f.\" % (analytic_rt, analytic_er, simulated_rt, simulated_er))\n",
    "\n",
    "n = 1000\n",
    "x0 = -0.3\n",
    "t0 = 0.4\n",
    "a = 0\n",
    "s = 1\n",
    "z = 3\n",
    "dt = 0.0001\n",
    "\n",
    "\n",
    "sim = ddm_distr_cython(n, x0, t0, a, s, z, dt)\n",
    "\n",
    "analytic_rt, analytic_er = ddm_analytic(x0, t0, a, s, z)\n",
    "\n",
    "simulated_rt = np.mean(abs(sim))\n",
    "simulated_er = np.sum(sim<0) / n\n",
    "\n",
    "print(\"Zero drift: Analytic RT, ER: %f, %f; simulated RT, ER: %f, %f.\" % (analytic_rt, analytic_er, simulated_rt, simulated_er))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10c9d2b38>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAECCAYAAAASDQdFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmUpNV93vFv7V177z0rDGjQZYQYJIMNBoSQLCLZAhkv\nx7F9lBOTYMUcJ16UOLHig+3jHOvYkUPsxAfLETKOE1uOLAlZFhYSkiyhGSEkFsEAw5197+m9u6qr\nurq2N3+8VTM1M93T1T3VXdvzOfRh6r31Vv+6q+upW/e97309juMgIiLdwdvsAkREZOMo9EVEuohC\nX0Skiyj0RUS6iEJfRKSLKPRFRLqIf6U7GGM8wKPATUAOeNBae6Sm/T7gYaAAPG6tfayy/TeBDwAB\n4FFr7eONL19ERFajnp7+/UDIWns78BHgkWqDMcZfuf0e4G7gQ8aYIWPMO4EfruxzN7C9wXWLiMga\n1BP6dwJPAVhrnwNuqWnbBRy01qastQXgW8A7gfcCrxpjPg98AfhiQ6sWEZE1qSf0E8Bcze2iMca7\nTNt8ZdsgcDPw08BDwN9ceakiInKl6gn9FBCv3cdaW65pS9S0xYFZYAr4srW2aK09AOSMMYONKFhE\nRNZuxQO5wF7gXuAzxpjbgH01bfuBncaYXiALvAP4GLAI/Arw340xW4AI7hvBshzHcTwez+p/AhGR\n7raq4PSstOBazeyd3ZVND+AO3USttY8ZY94P/E7lG3/SWvvxyn5/ALy7sv0j1tqvrlCLMzGRXk3t\nTTE0FEd1NkY71Aiqs9FUZ2MNDcVXFfor9vSttQ7uuHytAzXtTwJPLrHfb66mEBERWX86OUtEpIso\n9EVEuohCX0Skiyj0RUS6iEJfRKSLKPRFRLqIQl9EpIso9EVEuohCX0Skiyj0RUS6iEJfRKSLKPRF\nRLqIQl9EpIso9EVEuohCX0Skiyj0RUS6iEJfRKSLKPRFRLpIPRdGF+k4juOQTqeWbIvHE3g8q7rs\nqEjbUOhLV0qnUzz93CHCkegF2xeyGe65dSeJRLJJlYmsL4W+dK1wJEokGm92GSIbSmP6IiJdRKEv\nItJFFPoiIl1EoS8i0kUU+iIiXUShLyLSRRT6IiJdRKEvItJFVjw5yxjjAR4FbgJywIPW2iM17fcB\nDwMF4HFr7WOV7S8Ac5W7HbXW/usG1y4iIqtUzxm59wMha+3txphbgUcq2zDG+Cu3bwYWgL3GmL8H\nUgDW2nevS9UiIrIm9Qzv3Ak8BWCtfQ64paZtF3DQWpuy1haAPcBduJ8KosaYLxtjvlp5sxARkSar\nJ/QTnB+mASgaY7zLtKWBJJABPmatfS/wEPDXNfuIiEiT1BPEKaB2VSqvtbZc05aoaYsDs8BB4K8B\nrLUHgSlg8xVXKyIiV6SeMf29wL3AZ4wxtwH7atr2AzuNMb1AFngH8DHgXwE3Ar9sjNmC+2YwutI3\nGhpqjxUPVWfjNKvGYLBMLDpNNNZzwXYveQYH4ySTF9bVDr9LUJ2N1i51robHcZzL3qFm9s7uyqYH\ncA/cRq21jxlj3g/8DuABPmmt/bgxJgA8DlwNlIH/ZK39zgq1OBMT6bX/JBtkaCiO6myMZtaYSs2x\nZ9/oJUsrZzNp7rxx8wXr6bfD7xJUZ6O1UZ2ruuLPij19a62DOy5f60BN+5PAkxftUwA+uJpCRERk\n/engqohIF1Hoi4h0EYW+iEgXUeiLiHQRhb6ISBdR6IuIdBGFvohIF1Hoi4h0EYW+iEgXUeiLiHQR\nhb6ISBdR6IuIdBGFvohIF1Hoi4h0EYW+iEgXUeiLiHQRhb6ISBdR6IuIdBGFvohIF1Hoi4h0EYW+\niEgX8Te7AJF25jgO6XRqybZ4PIHH49ngikQuT6EvcgXS6RRPP3eIcCR6wfaFbIZ7bt1JIpFsUmUi\nS1Poi1yhcCRKJBpvdhkiddGYvohIF1Hoi4h0EYW+iEgXUeiLiHQRhb6ISBfR7B3pWJebQ59Op8DZ\n4IJEWsCKoW+M8QCPAjcBOeBBa+2Rmvb7gIeBAvC4tfaxmrZh4HngPdbaAw2uXeSylptDDzA9OUYk\nmiAS01RL6S71DO/cD4SstbcDHwEeqTYYY/yV2+8B7gY+ZIwZqmn7OJBtcM0idavOob/4qyd86RuB\nSDeoJ/TvBJ4CsNY+B9xS07YLOGitTVlrC8Ae4K5K2x8BfwacaVy5IiJyJeoJ/QQwV3O7aIzxLtOW\nBpLGmH8JjFtrnwa0+IiISIuo50BuCqgd+PRaa8s1bYmatjgwC/wK4Bhj7gHeBvyVMeYD1trxy32j\noaH2GF9VnY2znjUGg2Vi0WmisZ4LtheKJZ49kGWkN8SP7LiwzUuewcE4yeSFdS1X53LfY7nHWW/t\n8JyD6mymekJ/L3Av8BljzG3Avpq2/cBOY0wv7tj9XcDHrLWfq97BGPNPwL9ZKfABJibSq6m9KYaG\n4qqzQda7xlQqzXxmkTK5C7a/fGiSmfkiM/NFdpycZrgvcq4tm1lkcjJNPn/+Q/Dl6lzueyz1OOut\nHZ5zUJ2Ntto3pnr+Ip8AFo0xe4H/Bvy6MebnjDEPWmuLwIeBr+C+OTxmrR29aH9NjJOWkc0Vee3o\nNH6fe/s7r41RLutPVLrHij19a60DPHTR5gM17U8CT15m/3evuTqRBnvp4ATFksONV/Uwmy1zcjLP\n/uMz3HBNf7NLE9kQOiNXusbs/CKHT6fojQXZNhBg17YwoYCPVw5PqbcvXUOhL13j9EQGgLdeO4DH\n4yHo97J9JEahWGZ2frHJ1YlsDIW+dI2J2QUAhvvC57YNJt1ZN1NzuSX3Eek0Cn3pCo7jMD6zQCTk\nJ9pz/lDWQCX0JxX60iUU+tIV5hcK5PIlhvrCF1ysvC8Wwuf1KPSlayj0pStMzLqhPtR70UlUXg/9\niRCz84sUS+WldhXpKAp96QrV8fyh3vAlbQPJHhwHZlI6mCudT6EvXWF8ZqHSq++5pG1Q4/rSRRT6\n0vEKxTKz6UUGkz34vJeu/zeQcHv/UymFvnQ+hb50vMm5BRyWHtoBSEQDBPxe9fSlKyj0peMtdxC3\nyuPxMJDoIZXJky/qYK50NoW+dLzpyrBNdU7+Uqpts/OFDalJpFkU+tLx5jJ5An4vkdDy6wv2xoIA\npLPFjSpLpCkU+tLRymWHVCZPMhq84KSsiyWrob+g0JfOptCXjjafK+I450N9OYmoQl+6g0JfOlp1\nuCYZC132fkG/j3DIp+Ed6XgKfeloqUqI90Yv39MHt7efXSxpBo90NIW+dLTzPf2VQz9ZeWOYnNV8\nfelcCn3paKlsEZ/XQzQcWPG+1XH9MYW+dDCFvnSssuOQXiiSiAbxXmbmTlUy6o77j89q4TXpXAp9\n6Vgz6TylslPX0A64yzEAjM+svaefyuR50U6wsFha82OIrKflz1YRaXNjlfDuXWHmTlU0HMDrhfE1\nDO+Uyg6vHpli3+Fpyo7D1FwP99yybdWPI7LeFPrSsaqhn6xj5g6A1+Mh3uNnfDaH4ziXPZnrYi8d\nmOD1YzOEQ35CAS+j0zlePTzGW9906X3j8cSqHlukkRT60rHOTrsXTql3eAcgHvEzN5ljdj5PX7y+\nTwhlx+HImRShgI8fv3MHc/N5vvTcCf7umdPMZEoXBPxCNsM9t+4kkUiu7ocRaRCN6UvHGpvJ4fFA\nPFJ/6MfCbj9odCpT9z4Ts3ly+RJXb4oTDPgY6gsz0utnbsFhKuMhEo2f+wpHoqv+OUQaSaEvHclx\nHM7O5Ij1+Je8cMpyEpFq6Gfr3ufkhPuJ4prN8XPbzBb3U8L+YzN1P47IRlDoS0eay7i972qI1yu+\nyp5+sVTm9NQCkZCf4b7zF2mJ9fhIRnxMzC5Q0Bm+0kIU+tKRRifd0I6vNvQjfjzAmcn6Qv+NEykK\nRYcdm+OXHJwdTPgpO+71eUVahUJfOtKZyvDManv6fp+XgWSIUxMZHMdZ8f4vHpoGYMfmxCVtg3F3\n3v9qjg+IrLcVXxHGGA/wKHATkAMetNYeqWm/D3gYKACPW2sfM8Z4gU8ABigDv2StfX0d6hdZ0pmp\nak9/5eUXLra5P8y+o7OksoXLTvcslx1eP54iEvIxkLh0pk9/zI/X4+HsdP3HB0TWWz09/fuBkLX2\nduAjwCPVBmOMv3L7PcDdwIeMMUPAfYBjrb0T9w3how2uW+SyRiczeIB42LfqfTf3u5dOPD0xf9n7\nnRhPk8uXGO4NLTnv3ufzMNTbw3RqkVxeZ+hKa6gn9O8EngKw1j4H3FLTtgs4aK1NWWsLwB7gLmvt\n3wMfqtxnB6ApDLKhzkxl6U8E8ftWP4K5qd89IHt64vLDMm8cnwVguHf5TwObByIAjKm3Ly2inldE\nApiruV2sDN8s1ZYGkgDW2rIx5i+BPwH++spLFanP/EKBVCbPSO/yF0K/nM3V0J+8fE/fnnD7MoPJ\n5U/i2jTgzsvXuL60inqOcqWAeM1tr7W2XNNWewQrDsxWb1hrf8EYMwx81xizy1p72WkMQ0PxyzW3\nDNXZOOtR48TRKQB2bE0Si/qJxi4N/4VMEK83QPyiNi95du0cxOf1MDabO1ffxXWWyg4HT88x0h9m\nuD9yyfeoPn7fQC8B/ynGZhaIx3rwkmdwME4yuT7PTTs856A6m6me0N8L3At8xhhzG7Cvpm0/sNMY\n0wtkgXcAHzPGfBDYZq39A9yDvyXcA7qXNTGRXmX5G29oKK46G2S9atx/eBKARI+X+cwiZS5dQC2T\nyeP1lgiFL2zLZhaZncmwqT/C8dEU4+MphocTl9R57GyKbK7I7mt6l/we5x9/kaHeMGcmM0xMz1Mu\nLDI5mSafb/zEuXZ4zkF1Ntpq35jq+ct7Alg0xuwF/hvw68aYnzPGPGitLQIfBr6C++bwSWvtKPA5\n4O3GmG8CXwJ+1VqrRcplQ1Tn2I/0rW14B2DrUJRcvsR0auk/2+p4/nVbV37BDSTdOqbm9BKQ5lux\np2+tdYCHLtp8oKb9SeDJi/bJAv+8EQWKrFZ1uuZIXw+n19hT2zrojsWfnpzn+p1Dl7S/URnPf9OW\nGK8evXyYV6dzTqVy9K1iHSCR9aCTs6TjjE5mSMaCREJrX0R261AMWHoGT6lc5uCpWUb6wvTWsYJn\ntac/ndJlGKX5FPrSUXL5IlOpRbYMXNlqlluH3P1PLRH6R86kWFgscf3VfXU9ViTkpyfoY2pOoS/N\np9CXjlJdHfNKQ38oGSYY8HLsbOqSthcPTADw9usG63osj8fDQLKHTK7Iok7SkiZT6EvbcxyHVGqO\nVGqOI6fcmTt9MS/pdApWXj5nSV6vhxt29DM6leXk2PnjAo7j8OKBCUJBH7uu7q/78QYS7hDPzHxh\nbQWJNIhCX9peOp3i6ecOsWffKM8fcEN/fCbDPz1/hFxu7Stc3mzcA7jf3nfm3LZTExkmZnPsvnaA\ngL/+l091XF+hL82m0JeOEI5EiUTjZCsTaUYGe+kJX9kQz9sqJ2l9+5XRc9uqQzs/8OZLZ/RcTnUG\nj0Jfmk2hLx1lbn6RYMBLT3D1C61dLNITYNeOPo6cnmN81v3E8NKBCfw+D7vfNLCqxwqH/IRDPoW+\nNJ1CXzpGqVwmnS2QjC696uVa3GKGAXjBjjM2k+XE+Dy7ru4nvMrpoB6Ph/5EDwuLJdJZBb80j0Jf\nOkYqU8CBuubO1+vt1w3i9Xr48ndP8ruPfw+AW8zqhnaqqgdzT05oxU1pHoW+dIy5eXdAP9nA0I9H\nguzeOUgqkycc9PEz79rJHbs3r+mx+ivj+memdPlEaZ61n7Io0mLmMnkAktHllzpei1/952/npdfP\n8tZr+9e0Pn9Vf7xycZZJ9fSleRT60jHm5iuh38CePsBgb5i31Xki1uVEw34CPg+nJ9XTl+bR8I50\njNn5Rfw+D9Ge1uzLeDwektEAE7M5nZkrTaPQl45QdpzKhcwbN3NnPfTGAjjAqRWuyiWyXhT60hEy\nuRLlstPwoZ1GS0YDAJwcV+hLcyj0pSOks0Wg8eP5jdYbc4eeTo4p9KU5FPrSEVKVE56S0dYO/UQk\ngNcLJ8Zb/zJ80pkU+tIR5jJuT78v3tjpmo3m83oY6e3h1HiGsrPGJUBFroBCXzrCXKaA3+chFg40\nu5QVbR2MsFgoMTGjqZuy8RT60vaKpTLphSK9sdaeuVO1dTACwAkdzJUmUOhL2xufXcRxoLfFh3aq\ntgyEATgxpnF92XgKfWl7o5W1bPpi7RH6Wwfd0Ne0TWkGhb60veoCZq1+ELcqFg7QFw8p9KUpFPrS\n9kan3dDvjbf2dM1a24djzKQXSWfzzS5FuoxCX9re6NQCPQEvPcHWXHNnKduHY4CGeGTjKfSlrWVz\nRWbm8ySirT9Vs9ZVI3EATujMXNlgCn1pa6crC5clo+3Tywe46lxPXzN4ZGO11ytF5CKnJjIAJCPr\n29N3HId0OnXJ9nQ6BWs4sXaoL0wo4NNcfdlwCn1pa6cm3NBMrHNPP51O8fRzhwhHohdsn54cIxJN\nEInFV/V4Xo+HbcNRjo2mKRRLBPy+RpYrsqwVXynGGA/wKHATkAMetNYeqWm/D3gYKACPW2sfM8b4\ngb8AdgBB4Pettf/Q+PKl250cm8frgURk/fsv4UiUSPTCcM9m1t5Tv2o4zuHTKc5MZrl60+reNETW\nqp4x/fuBkLX2duAjwCPVhkq4PwK8B7gb+JAxZgj4IDBprb0L+FHgTxtctwjlssOJ8TQjfeErunZt\ns1Rn8OjMXNlI9bxS7gSeArDWPgfcUtO2CzhorU1ZawvAHuAu4NO4vf/q9yg0rGKRitGpDPlCme3D\nkWaXsibbRzRtUzZePZ+JE8Bcze2iMcZrrS0v0ZYGktbaLIAxJg78HfBbDapX5JxjZ90e8vahCNB+\n15zdNhTD41FPXzZWPaGfAmoHHKuBX21L1LTFgVkAY8x24HPAn1pr/189xQwNtce4pupsnCupcWzu\nKAA3XjfIsTOzRGM9F7QvZIJ4vQHiF22/XJuXPIODcZLJC+saHIwTi07X/T2W2+5xFgkEygSDZYJB\nD1sGI5wYn8cfKNGbTF7xKqHt8JyD6mymekJ/L3Av8BljzG3Avpq2/cBOY0wvkMUd2vmYMWYE+DLw\ny9baf6q3mImJ1u/xDA3FVWeDXGmNbxydxuvxEPHDfGaRMrkL2jOZPF5viVA4d8m+y7VlM4tMTqbJ\n58+PfA4NxZmcTK/qeyy3fXJimidOjdLbPwCAz+OQy5f45Ode4CfeaUgkkqv/RdTU2erPOajORlvt\nG1M9of8EcI8xZm/l9gPGmJ8DopWZOh8GvgJ4gMestaPGmD8GeoGHjTG/jTuT+UettYurqk5kGaVy\nmRPjabYMRgj62+sgbk84cm4W0Eh/kRPjC2SLmj0tG2PFvzRrrQM8dNHmAzXtTwJPXrTPrwG/1ogC\nRZYyOpUlXyi3/VTHgaQ7/DOT1lwH2Rjt1UUSqTheOYi7Y1NihXu2tv5ECA8wM6/Ql42h0Je2dOxc\n6Ld3T9/v85KMBZmdL1Au60Lpsv4U+tKWjp9N4/V4zp3g1M4Gkj2Uyg5js5cecBZpNIW+tJ1iqcyJ\nsTRbBqMEA+2/Zs1Awh3XPzmebXIl0g0U+tJ2jo+lyRfLXLdt7dMbW0n1YO7JyoqhIutJoS9t5+BJ\n9yTwTgn9vngIj0c9fdkYCn1pOwdPzQJw3bbeJlfSGH6fl2TEz+nJLKVyeeUdRK6AQl/aiuM4HDw1\nx0AidG5YpBP0xYMUSg6nxjXEI+tLoS9t5ex0lvmFQsf08qsG4kEADp2eW+GeIldGoS9t5eCpynj+\n9g4L/YQb+ofPKPRlfWnBD2kbjuPw2pFxALb0+kil3IBc63VqW0ks7CMS8nHolEJf1pdCX9pGOp3i\ntaOzBPweDp2e4fAZ94DuWq9T20o8Hg87NsV4/fgcc/OLJGOhZpckHUrDO9I25jJ5svkyI30RorEE\nkWicSDROTzi68s5t4JpN7s9x6HSqyZVIJ1PoS9s4cMpdb2e4vz0vj7iSHZvcJSUO62CurCMN70jb\neOOE2wPeOtgZPftajuPQFy7h9YA9MUUqNXSuLR5PXPEVtUSqFPrSFsplhzdOzhEOeumNBZtdTsMt\nZDN8Z980iUiA42MZvvnyGXxeDwvZDPfcuvOKrqglUkvDO9IWjp1Nk8mVGOnr6dheb084wshAlLID\nuVKASDROONJ5n2qkuRT60hZePTIFwKa+zp7VMtwbBmBsWuvwyPpQ6Etb2Hd0Cq8Hhjs89EcqB6nP\nTi80uRLpVBrTl5Y3v1DgyJkU12yKrftF0B3HcU/2qhEMljfsBLBIj59kNMj4TJaSrqQl60ChLy3v\n9WPTOA5cv339r4e7kM3wzRen6e0fOLctFp3mxPETG3YC2KaBCPbELFNzC3TgMWtpMg3vSMt78cAE\nAG+5emNmsPSEI+dO/IpE40RjiQ09AWyThnhkHSn0paUtFkq8fGiK4b4wWwfDzS5nQ5wb15/SwVxp\nPIW+tLR9h6dYLJT4weuHO3aq5sV6gj764iHGZxc0ri8Np9CXlvbdN9xVNX/w+uEmV7KxNg9EKJcd\nplL5ZpciHUahLy1rMV/ilUOTjPRH2D4ca3Y5G6o6rj8+u9jkSqTTKPSlZb18eJJ8sdxVQztVI/0R\nvB44O6PQl8ZS6EvLcRyHVGqOb79yGoC3bI+QSs11xMVS6hXwexnpjzA7X2B2XkM80jgKfWk56XSK\nL+45wL5jsyQjfg6fnmHPvlH+6fkj5HLdM42xOqT16jEttSyNs+LJWcYYD/AocBOQAx601h6pab8P\neBgoAI9bax+rabsV+ANr7bsaXbh0trEUOA6Yq/uJxtyTsrKZ+SZXtbG2Dcf47v5xXj02y4/d3uxq\npFPU09O/HwhZa28HPgI8Um0wxvgrt98D3A18yBgzVGn7DeATQGcvliINV3Ycjp7N4vN6uHbL+p+F\n26pi4QDJqJ+Dp9IsLBabXY50iHpC/07gKQBr7XPALTVtu4CD1tqUtbYA7AHuqrQdAn6igbVKlzh0\nOs18rsSOTXGCAV+zy2mqLf09lMoOrx2dbnYp0iHqCf0EUDuoWDTGeJdpSwNJAGvtE4C6J7Jqz74+\nCcB123ubXEnzbR7oAeD7hyabXIl0inpCPwXUrjLltdaWa9pqP3/HgdkG1SZdaCa9yCtHZklE/Az1\n9jS7nKbriwVIRgO8fGiSQrG88g4iK6hnlc29wL3AZ4wxtwH7atr2AzuNMb1AFndo52MX7V/3BOuh\nofVfwbARVGfjXFzjF79zglLZ4a3XJEnEL1xrZyETxOsNEI/1XNH2tewTjTbme692u5c8d+zexD8+\ne5KjExnu2L3lkp+lVjs856A6m6me0H8CuMcYs7dy+wFjzM8BUWvtY8aYDwNfwQ33x6y1oxftX/fM\n6omJdL13bZqhobjqbJCLa8zkCjz57aMkIgFGeoOk53MX3D+TyeP1lgiFr2z7aveJx3oa9r1Xuz2b\nWeTGq3r5x2fhS3uO8ObNy4dQOzznoDobbbVvTCuGvrXWAR66aPOBmvYngSeX2fc4oMlmUpevv3ia\nxXyJ9968CZ+3S87CqsPmgTA7NsXZd2SauflFkjFNiJO108lZ0hIWCyW++vxJIiE/t98w1OxyWs4d\nN26m7Dg8+9pYs0uRNqfQl5bwtRdOkc4WePfNW+kJdvc0zaXc+pYR/D4Pe18dxXH0KUjWTqEvTTeX\nyfPFbx8jFg7wvh+6qtnltKRYOMDbdg5yeiLDwVNalkHWTqEvTVNdWO3TX3uDXL7Ee2/ZRDGf7aqF\n1Vbjnh/cDsCTzx5vciXSzhT60jSpVIrPfsPy7OuTxMN+HKfYlQur1eu6bb2Y7b3sOzLF8bOtP6tE\nWpNCX5qmVC7z2kl3vfgfessIsViCSDS+oRchbzf33r4DgC8+e6yZZUgbU+hL0/zDnhNMpwvs2Bxn\n61B3XRlrrd6yo48dm+K8aCc4PZlpdjnShhT60hTHzqZ44pnjhINebn3LSLPLaVmO45BOp0il5s5d\nSOY9bx/GAf72qwc0k0dWrZ4zckUaKpsr8udfeJ1S2eGWN/cR6vKVNC9nIZvhmy9O09s/cG6b4zgM\nJvy8dmyG770xzg/t0pum1E89fdlQ5bLDn3/hNcams/zYD29npE+Lqq2kJxwhEo2f+4rGEtzy5n78\nPg+f+upBsjktZiv1U+jLhvrsNw+z78gUb722n5/9kWubXU7bioX9/LObNzOXyfN33zjU7HKkjSj0\nZcM88/IZvvTcCUb6I/zSB27A6617AVa5iOM4/ODOMJv7w3zz+2f4+vNHmJtzx/01zi+XozF9aZjq\nQcelHBsv8FdPWWLhAL/607uJ9AQArQ+/VgvZDHtfnmb3tUkm5nJ86uvHSS8UKObmuefWnSQSyWaX\nKC1KoS8Nk06nePq5Q4QjF86zPzuZ4rsHs/h8Hn7lp3ezqT/SpAo7S084Qv9gH3fu9vONl87wjZcn\nuOuGvmaXJS1OwzvSUOFI9IKDjmVPD88fylIolvnQfTewc6t6oI121UicG6/tZ24+z57XpsnlS80u\nSVqYQl/WzWKhxNdeOEWuUOZHbxnius3Bc/PNU6k55ubmtMZOg7ztukGu39HHzHyBT37pMIWigl+W\npuEdWRelUplvvHiauUyea4aD5BdS7Nl34UXVctlZ8ISIxDrvknQbzePx8K4f2E4qvcDB02n++6df\n5t/91G7CIb3E5ULq6UvDOY7Dnn1nGZtZ4OpNca7fGrpkrnkkGicc0dh+I3m9Hm69vo/d1/TyxolZ\n/uhvv8/8QqHZZUmLUehLw71gJzh+Ns1wX5g7b9yEx6OpmRvF5/XwL997LXe8dRNHR1P83l9+Tyty\nygUU+tJQh85keP3YDIlokHe9fSs+n/7ENprP6+GB9+/iA3fsYHIux0f/7ws88/IZzd8XQGP60kCv\nHJnl+4fn6An6+JGbtxLSZQ+bxuvxcP87rmXH5gSf+IfX+csvvcGz+07zM3dfzUDiwgurx+MJfRrr\nIgp9aYj9x6b5P189gs/r4d03byMeCTa7pK508Qly1w4H+I8/s4tPff0w9lSa3/+bV7luSwyzLUYw\n4GUhm9EgjgMfAAALCUlEQVTJXF1GoS9X7NWjU/zPz+6jXIbbdvUxmNQias2y1KqcAFvjOaLXRLGj\ni9hT8xwdy3LjtQNcNaiD6d1GoS9X5Lv7x3jsi/sBePDH3sTUXLbJFUl1plStbGaeaMzHrp1bsCdm\n2XdkihfsBK8f89ETCvIjtyS0FlKXUOjLmhSKZT799UN87cVThAI+/u1P3sj2AR979in0W5nf5+WG\na/rZuS3Jq0em2X98mk99/TjPvDLJT73zTdy0c0Dj+x1OoS+rUnYcXrATfP5bRxidyrJ1MMpD97+V\nLYNRUqm5ZpcndQoFfNxshrhq0M/4TI4XD6f4H599hWs3x7jvtq1cszmmA7wdSqEvKyo7DmcmM7x4\nYILv7h/nzGQGrwfuuGGID9y+lVCgeO5SflpWob14Sjn6Qovc8/YhXj2e5sjoPH/yhGWk188D73sz\nb96xqdklSoMp9OWchcUi4zMLTMye/zozleXEWPrcIl4+r4ebr+snHswz2B/ge2+Mn9t/enKMSDSh\nZRXajLtaZz+bhvsZn8nygp1gbDbHH/6/13nH7hl+/M5r6YuHVn4gaQsK/S5SO51vLpPn2NkMp6cW\nODOZ5czUAtPp/CX7eICh3h5uuDrJrqsT3HB1L8V8hpePZpY8WCjtbbgvwvtuvYrDJyc5PJrlmZdH\nefa1Me7avYW73raF7cOxZpcoV2jF0DfGeIBHgZuAHPCgtfZITft9wMNAAXjcWvvYSvvIxsvli7z0\nxhmefv4E0/NlUtkLr6sa8MFwb5B42E+0x0+0x0cpN0e0x8/AYL/7GIt5Xjgwrh59h/N4PGwZ6OGn\n7trBq8ezfH7PUb724im+9uIptg/HePPWKNdtTbBlMEw87D837r/UMYClLqwTDJZJpdI6ZtAk9fT0\n7wdC1trbjTG3Ao9UtmGM8Vdu3wwsAHuNMX8P3LncPrL+HMdhanqG79sxDpxKceBUmmNj85QrF6ry\n+zxsGYyyqT9Mf6IHZ3GWcChA/+DwBY8zOV7C6/WpR9+FHMchm0lz0zUJ3nrVDbx2fI7v7J/Enpzj\n5Pg8X3tpDICAz0NP0EfA5zDSHyUUDOL1evBVvkqlAqOTaYKBAB4PBPxe+hMhSvkc77/jTWwd0Wyh\njVZP6N8JPAVgrX3OGHNLTdsu4KC1NgVgjPkW8E7ghy+zjzSI4zgsLJaYmV9kdn6RmdQiZyYzHDkz\nw5HReYql80dV+2IBEqEimweivGnHZnze82viTI4vfYlD6V5LneT1lquiDPWkWShHWCgFmZ3Pk8rk\nyeVLpBdKTM9f7u+odrVPt9Pw3MFXiIUDbBuKsm0oxrbhGNuGYmwdjGoJj3VUT+gngNq5eEVjjNda\nW16ibR5IAvHL7LOuUtk8pZKD4ziUyg5lx6Fcdig77hrvhVKZYrFMseSc+/f5bWUKJYdCscSR42dw\nPD7KZfB63d6x3+vh2qsGCYd6CAZ8BPxeQn4vgYD7B1ouO5RKZUqOQ7HkkC+UKBTL5AslFgtlCsUS\n+WKZxXyJzEKOYqlM7RpYjgOBQMAdSK9sLzsO+UKZbG6RfKFMvuh+FYplFhbdx1tKIhpg80CUzQMR\nRvojhAI+JsdH8Xp9FwS+yHKWOskrGp0n7vXRPzh0wfb59Bw3XBUlEomfe82Vyg6pdJrXT2QJhSOU\nHVjMlyg5cHpsBg8O43MF3jgxyxsnZs89lnscKczmgQixSIBwyI/X48FxwMGh8h+O4+A4sJhfBMdd\nWtrr8eD1uv8O9/Tg83rOffLwej2Uyw7ZhRylUplipcZSyf1/TyhEwO/F7/Pi93noTYYp5IuEAj6C\nAS8hv49g5d9Bv/v/gN+L1+PB4/Hg9YDH6yEeDrT0p5d6Qj+FG+JVteGdwg3+qjgws8I+6+bz3zrC\nF/YeW9fv8Z0DzVum1gP4fFQ+OkPQ77BjJMJAMkIyGiAZDTDc20MiVOTg6Dxlx316S/ks2TzkFjJ4\nvX6ymQt/hqZtz2bJ5Uob/n1Xu4+XfOv97lqsztnpCb4+eopk74XX6J2ZniQaTRCL+MEDoR6IRUM4\nmQUWFxe5ZmcfxVKQ+VyZ9EKJmfk8gUAP43MFXj48RTt6361X8TPv2tnsMpZVT+jvBe4FPmOMuQ3Y\nV9O2H9hpjOkFssA7gI9V2pbbZzmeoaErOzD4iz95E7/4kzdd0WN0ituaXUDH2d3sAuqkOuXyPCut\nsV0zE6f6LD2Ae+A2Wpmp837gd3A7op+01n58qX2stQfW4wcQEZH6rRj6IiLSOXRET0Skiyj0RUS6\niEJfRKSLKPRFRLpIUxdcM8YkgL8FYrhr9HzQWjtemeb5x7in8T1trf29JpaJMcbL+eUmQsDvWmv/\nsdXqBDDGXA98Bxi21uZbrcbKc/5/cc/vCAAfrpy13VJ1wsrrTjVTZQmUvwB2AEHg94HXgb8EysCr\n1tpfblZ9tYwxw8DzwHuAEi1YI4Ax5jeBD+D+XT4KPEML1Vp5zv837nNeBH6RNfw+m93T/wXgFWvt\nXcCngd+obP8z4Gette8AbjXGNHvy/b8A/JV67geqZ160VJ3GmDjwR7gBVdVSNQIfBr5qrb0bd/rv\no5XtrVYn1Kw7BXwE942/VXwQmKy8dt4H/Cluff/ZWvtOwGuM+fFmFgjngurjuOfxQAvWCGCMeSfw\nw5Xn+m7gKlqv1h8DfNbaO4D/AnyUNdTY7NDfx/kzehNAoRJcQWvtscr2L+P2EJrpvcAZY8wXgf8F\n/EOL1vm/cMMpC+feBFqtxkeAP6/8OwAstGidcNG6U0ArrSH1adzVbQF8uD2/H7DWfquy7Uu0xu/w\nj3Df0M/gnsvTijWC+xp/1RjzeeALwBdpvVoPAP7KJ9Ak7qfiVde4YcM7xph/Bfw67rIZ1dVl/i3w\nz4wxrwF9uGf0JnCXcahKA9c0qc6qCWDBWnuvMeYu3I9TP9+sOpep8QTwKWvtvsofBbTW77L6nD9g\nrX3BGLMJ+D/ArzS7zsu43LpTTWWtrX1j/zvgt3ADtiqNGwxNY4z5BWDcWvu0MeY/VzbXdjSbXmON\nQdze/b3AtbjB32q1zuO+Lt4ABoD7cDOzqq4aNyz0rbV/gTsGeY4x5rPAH1prP2GMuRH4HG7v6uL1\nfGbZIMvU+Sncd36stc8YY67DDYOm1LlMjQeAf22MeRDYBHwF94+ipX6XAJXn+m+Af2+t3VMJrqbV\neRlNWUOqXsaY7bivmT+11v6tMea/1jS3wu/wAaBsjLkH97jIXwG1K7W1Qo1VU8B+a20ROGCMyQHb\natpbodZfB56y1v6WMWYr8A3c4zlVddXY7OGdac73pCaAuLU2DSwaY66p9FjfC3xruQfYIHtwx9Oo\njDWfsNbO00J1WmvfbK19t7X2XcBZ4J5W/F0aY96COzTx89barwC0Yp0Vezn/vNe7htSGMMaM4A6D\n/Udr7f+ubH6p8kkU4Edp8u/QWvtOa+27Kn+T38c9NvalVqqxxh7cYyMYY7YAUeBrlbF+aI1aa/Ny\nFrfT/tJqa2z25RJ/G3jMGPPLlVoerGx/CLcn6AW+Yq39XpPqq/oE8GfGmGcrt3+p8v9Wq7OqOpwC\nbq2tVONHcWdA/Ukl4GettT9Ba/4unwDuMcbsrdx+oJnFXOQjQC/wsDHmt3Gf818F/qcxJoC7GOJn\nmljfcv4D8IlWq9Fa+6Qx5h3GmO/ivnYeAo7h5lOr1PrHwF8YY57BPR72m8ALrLJGrb0jItJFmj28\nIyIiG0ihLyLSRRT6IiJdRKEvItJFFPoiIl1EoS8i0kUU+iIiXUShLyLSRf4/+LmD5Gq7EAoAAAAA\nSUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1056068d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "sns.distplot(sim)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}