Regression - Orientation tuning.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "$\\newcommand{\\normalN}{\\mathcal{N}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Introduction\n",
    "\n",
    "The working example we will use in this notebook is describing the orientation tuning of a neuron. To do this, we might perform an experiment where a monochrome grating is shown to the animal while the recording from the neuron of interest. Repeating the experiment with different grating orientations, we would obtain the firing rate $f$ of this neuron as a function of the orientation $θ$.\n",
    "\n",
    "<img src='./orientation-tuning_plot.png' style='width: 50%'/>\n",
    "\n",
    "Clearly this neuron responds preferentially to grating orientations around 90°, and with about twice the firing rate as the least preferential direction. The latter is an estimation on the *contrast ratio*, which is the specific quantity we want to evaluate. But to do so with any scientific credibility, we need to do a proper data analysis that takes into account the observation noise – “eyeballing” the number just won't cut it.\n",
    "\n",
    "*Remark* – Contrast is a measure of the amount of *signal* in an output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Loading and plotting the data\n",
    "\n",
    "This first step to any data analysis is to grab the data. We've provided the data for the plot above in the file `orientation-tuning_data.npz`. It comprises two Numpy arrays:\n",
    "\n",
    "* `θarr`: the array of orientations\n",
    "* `farr`: the array of corresponding measured firing rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f2841bc39e8>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGvNJREFUeJzt3X+QXWV9x/HPl7DgRh1XydbCSgydoWhL2kR2KJrR0tAW\n/DGQRtQw/SFUm6m19ccfaGhtaTudEodO7VgdNa2MWh2MFaQUVECDMsMIdQMoUMDBUpDVmhUJlGbV\nJH77x71LNjf33HvuPec5z3nOeb9mMtncPbnn2b33nu/zfJ/v8xxzdwEA2uuo2A0AAMRFIACAliMQ\nAEDLEQgAoOUIBADQcgQCAGg5AgEAtByBAABajkAAAC13dOwG5LFq1Spfs2ZN7GYAQFJ27979A3ef\nHnZcEoFgzZo1mpubi90MAEiKmT2c5zhSQwDQcgQCAGg5AgEAtByBAABajkAAAC1HIACAlkuifBQo\n6po753X5DQ/ou3sXdcLUpC4++xRtWj8Tu1lALRAI0HjX3DmvS66+W4v7D0qS5vcu6pKr75YkggEg\nUkNogctveODpILBkcf9BXX7DA5FaBNQLgQCN9929iyM9DrQNqSE03glTk5rvc9E/YWoyQmsOYd4C\ndcGIAI138dmnaHJixWGPTU6s0MVnnxKpRYfmLeb3Lsp1aN7imjvno7UJ7cWIAMkYtwe9dEydet+D\n5i2KtouRBkZFIEASBlX+SMMv8pvWz9TqYhhq3oIKKYyDQIAkZPWg//Lae/XjAz9N7sIXat4i5EgD\nzcUcAZKQ1VPeu7i/lqWh19w5rw3bd+mkbddrw/ZdR+T+Q81bUCGFcRAIkIRRe8oxL3x5JoI3rZ/R\nZZvXamZqUiZpZmpSl21eW7jXnvV7il0hhXoLFgjM7Aoz22Nm9yx77HIzu9/MvmlmnzOzqVDnR7Nk\n9aCfu3Ki7/ExL3x5F7BtWj+jW7dt1EPbX61bt20sJXVTxwop1F/IOYKPSfqApE8se+wmSZe4+wEz\ne6+kSyS9O2Ab0BBZlT+SDpsclaq98PWr0ImZnqljhRTqL1ggcPdbzGxNz2M3LvvnbZLOD3V+NM+g\nyp8YF76sCp2plRN6fN/+I46vapRStwop1F/MqqHfl7Qz65tmtlXSVklavXp1VW1CgmJd+LJSQMce\nfZQmJ1ZEG6UAo4oyWWxmfybpgKRPZR3j7jvcfdbdZ6enp6trHJI2rFqnTFmpnicW9weZCAZCqXxE\nYGYXSnqNpLPc3as+P5qr6sVUg9YCkJ5BSiodEZjZOZLeJelcd99X5bnRfFVvN02FDpoi2IjAzK6U\ndKakVWb2qKRL1akSOlbSTWYmSbe5+x+GagPapepqHSp00BQhq4Yu6PPwR0OdD4ix3TQpIDQBK4vR\nGCmnaqqc5AZ6sekcGiPVVA07hiI2AgEaJcVUTd4dQ7nPAEIhEACR5ZnkZtSAkJgjACLLs2No1aWx\naBcCASrFpOiR8kxyc58BhERqCJVJLb1RVU4+zyR3jNJYtAeBAJVJ6TaKVQetYZPcF599StTtttFs\nBAJUJqX0Rt2C1rBRAxVFKIJAgMqklN6oY9DKGjWklnJD/TBZjMqktPI3pXv/UlGEoggEqEyoG7aH\nkFLQquPoBWkhNYRK1WXl7/Kc+nMmJ2Qm7d23/4j8egp595RSbqgnAgFapzenvnfx0P2Fe/Prdbzw\n96KiCEURCNA6/XLqy9W1pDVL3tELlUXIQiBA6+TJnaeWXx82eqGyCIMwWYwkFdmqIk/uvGn5dSqL\nMAgjAiSnaO+2X059uWH59awUS51TL1QWYRACAQoJcfEb9pxFV/325tQHVQ31a1u/IDT38A911e75\n2qZeqCzCIAQCjC1E3jnPc5bRux23IigrCF15+3d00P2Ix8eZdA4RXKkswiDMEWBsIfLOeZ4z5Krf\nYXMPWcGmNwgMO37Q+S+5+m7N712U61AgLLpdd0qL+VA9RgQYW4i8c57nDNW7zTMayUqxrDDrGwxG\nDU4hN7tLZV0EqseIAGML0TPP85yherd5RiNZW09c8CsnZm5JMUqFE5O6iIERAcYWomee9zlD9G7z\nXIQHLd6afeHzjnhc0kjzKEzqIgYCAcYWYj+emHv85L0IZwWhfo9v2L5rpFQPk7qIgUCAQkL0zGPl\nskNchEdN9aS02R2ag0AAdIW4CI+T6mFSF1UjEADLlH0RJtWDFBAIgIBI9SAFBAIgMFI9qLtg6wjM\n7Aoz22Nm9yx77HVmdq+Z/dTMZkOdGwCQX8gFZR+TdE7PY/dI2izploDnBQCMIFhqyN1vMbM1PY/d\nJ0lmFuq0aJg6b+2cMn6vWK62cwRmtlXSVklavXp15NYghGEXI+6qFQa/V/Sq7V5D7r7D3WfdfXZ6\nejp2c1CyPLtscletMPi9oldtRwRotjy7bDZxA7Y6pGSa+HtFMbUdEaDZ8lyMQt53IIZQ9xoYVdN+\nryguZPnolZK+JukUM3vUzN5kZr9lZo9Keqmk683shlDnR73luRhlbfmc6qrcMlMyo2xt3atpv1cU\nF7Jq6IKMb30u1DlRTJVpizxbLzRtVW5ZKZmik71N+72iOOYIIKn6SpK8F6Mmrcot614DZdzFrEm/\nVxRHIICksLdIzJLSxaiM0VJZG9A1cbK3DpPobUYggKRmXlzKUtZoqayUTNPuYsa6hvgIBJDUvItL\nmcocLZUxCmra1tYxRqM4HOWjkEQlySB1Gy1tWj+jyzav1czUpEzSzNSkLtu8NtmLZt1+v23EiACS\nqCQZpI6jpZTmV4ap4++3bQgEeFqTLi5laloqpkx1mkTH+AgEwBApj5ZCVuPUbRId4zN3j92GoWZn\nZ31ubi52M1otzwWFEsB66b1QS52edlnzCRu27+qb0pmZmtSt2zYWfn4UZ2a73X3oTcCYLMZQefbI\nKbqPTpEtE9Bf6F1GsyZz5/cu8homhkCAofJcUIpcdOqyGVvThK7GGTSZy2uYFgIBhsrT8yty0WF/\n/DBC7zLar+R4OV7DdBAIMFSent/UyomR/+8S6sjDCL02ZPl6hiy8hmkgEGCoPD0/d4190WF//DCq\nWHi2af2Mbt22MTMY8Bqmgaoh5LJUEdSvSkSSTNL73rDu6aqh50xOyEzau2//0Aqi0NUtKUulEovX\nsJ6oGkKp8vT8lo553xvW6ccHfqrH9+3PNfnbtC0TypLSJDqvYdoYESQqVk8xT8+P+vJy8HtEUYwI\nGixmTzFPz4/68nIwiY6qsMVEgmJv2ztsT6KsTcQk9pofRYzN2EKMNFOZ52gzRgQJqntPkfryclS9\nNXiIkWYZz8mq8/AIBAmqe7kl9eXlqHoCNsTCvqLPmdKEecpIDSUohW17l9JHWROedQladVfl1uAh\nRppFnzN2GrQtGBEkKKVSvZh3PiOlMJoQI82iz1n3NGhTMCJIVCo3kYm11zw3RB9diJFm0efk7mXV\nYB0BGoka/PEsr/AZZXV43ucc9XlYsVxM3nUEjAjQSKQUxrM00ixzRFVk9Mrdy6pBIEAjkVIopk6T\ntKmkQVPGZDEaKeYkdRMwomoXRgRoJFIK+fXL4TOiapdggcDMrpD0Gkl73P3U7mPPk7RT0hpJ/y3p\n9e7+eKg2oN1IKQyXNRfw2tNmdNXu+bGrfdhWIi0hU0Mfk3ROz2PbJH3Z3U+W9OXuv5ETdfEoW9Zc\nwM33L4y9VoXVwOkJNiJw91vMbE3Pw+dJOrP79cclfUXSu0O1oUmoi0cIg+YCxh1R1WmiGflUPVn8\nfHf/Xvfr/5H0/IrPnyxu8I4QQqwmZqI5PdGqhryzki1zNZuZbTWzOTObW1hYqLBl9cSHCyGEqK6q\n+6aIOFLVgeD7Zna8JHX/3pN1oLvvcPdZd5+dnp6urIF1xYcLIYTYt4rS3fRUXT56raQ3Stre/fvf\nKj5/slLYcRRpKru6itLd9OQKBGY2K+nlkk6QtCjpHkk3DSr9NLMr1ZkYXmVmj0q6VJ0A8Bkze5Ok\nhyW9vlDrWyTUh4syP4RA6W5aBm46Z2YXSfoTSQ9J2q1OKucZkn5e0gZ1AsKfu/sjIRvJpnNhsKEX\n0GxlbTq3UtIGd+87I2lm6ySdLCloIGibqnrplPkhJkaj9TEwELj7B4d8/65ym4NQ6wX6feioREIs\nrIupl4FVQ2a20swuNbNvmtkPzeyx7tfvMTPKVQIIsV4ga6Xn1MqJvsdTiYTQWBdTL8PKRz+tzsKv\nl0l62N2Pk/TK7vd2hmxYW4XopWd96Nw1sMyPLS0QCqPRehkWCE5294+4+1PqLv5y93l3/xtJLw7e\nuhaqcqXnE4v7M2vI2S8GIbEupl6GBYIHzexNZnbM0gNmNm1m75bEGC6Aqld6blo/o1u3bdRD21+t\nW7dtPKxMlaE7QmHRWb0MCwRbJL1A0tclnWRmj0naJekYSW8I3LZWqstKT4buCGnQ+5yUZPW4eX1L\njFqqx83fEQNrW8pVyjoCM3uPpA9mrSA2s42SVrr7deM1E1UZdaUnW1ogBta2xDFsQdndkq4zsx9J\nukPSgjori0+WtE7SlyT9bdAWIgr2i0EIw0ampCTjGBYIznf3DWb2LnW2lzhe0pOSPilpa9aKYzQD\n+8VgSRmrgPMsIht2r2RWI4cxbLL4NDM7QdJvq7Nz6EckfUKdyWPqvIAWKKuUOE8l2qDCBkqawxkW\nCD6szr2FXyRpbtmf3d2/ATRcWaXEedI+g6qJKGkOZ9heQ++X9H4z+5C7v6WiNgGokbLy9sPSPkuy\nUpLMH4ST6w5lBAGgvcpaBVx0ERmrkcOJds9iAGkoaxVw0cWSrEYOp+pbVQJITJmlxEUq0ShpDoeV\nxQDQUHlXFpMaAoCWIxAAQMs1eo6gyCrEJqxgbMLPgPrhfdU8jQ0ERe6J2oT7qTbhZ0D98L5qpsYG\ngiK7GDZhB8RhqzDp0WEcTfhshJLySKmxgaDIKsQmrGDMautSD44eHcbRhM+GVP5FO/WRUmMni4us\nQmzCCsastq4wY78WjK0Jn40Qm9elvg9SYwNBkVWITVjBmPUzHMxYN5Jajw5xNOGzEeKinfpIqbGB\noMhy9hD3Da5a1s8w04AeHeJpwmcjxEU79ZFSY+cIpOLL2VN6c/eT9TNwC0oUkfpnI+8uqKNI/dau\njR0RoL8m9OiAIkKkt1L/XLHXEIDWSbnUcxR59xqKkhoys7dL+gNJJumf3P0fYrQDQDOMemFPPb1V\ntsoDgZmdqk4QOF3STyR90cyuc/cHq24LgPSlXsNfBzHmCF4s6XZ33+fuByR9VdLmCO0A0ACp1/DX\nQYxAcI+kl5vZcWa2UtKrJJ0YoR0AGiD1Gv46qDwQuPt9kt4r6UZJX5R0l6SDvceZ2VYzmzOzuYWF\nhYpbCSAVqdfw10GU8lF3/6i7n+bur5D0uKRv9Tlmh7vPuvvs9PR09Y0EkIQmrHaOLVbV0M+4+x4z\nW63O/MAZMdoBIH3cy7i4WCuLrzKz4yTtl/RWd98bqR0AGoBy0GKiBAJ3f3mM8wIAjsQWEwDQcgQC\nAGi5Ru8+CqB92rKPUJlaGQh4owDNxHYT42ldaijEbeoA1APbTYyndYGANwrQXGw3MZ7WpYbGeaOU\nmUoiLQWEE+LuY23QuhHBqPuSlJlKIi0FhFXFdhPX3DmvDdt36aRt12vD9l2N+Py2bkQw6r1Fh6WS\nRundZz3XO3bepctveIDRAVBQ6O0mqpqMrjpz0LpAMOobJStltPQGGOUNMSj9RHUDUI6Q200M6hhm\nnXPUi3qMyqfWBQJptDdKVs5xhdnIb4is58r7/wHENeoc47CLer8gMU6wKap1cwSjyso5HnTve/yg\nXn+/5xrl/wOIK2su0aW+8wWDLupZc4ZZncWQ1wYCwRCb1s/oss1rNTM1KZM0MzX59L/7OcoscxJp\n+XNlGae6oYmTV0AdDerM9Sv+GDSCyAoSK8z6/p+QlU+tTA2NKiuV1DvpLOnpkUJWXm/puXqHjNJ4\n1Q2spASqs3yOsV/PvTeFM6icNStIHHTX5MSKwteGUTAiGFPvSKFfFB+0UC1rpNF78R7W22eBHFCt\nTetndOu2jerfbz98FDConDWrh7886zDo2lAmRgQFLB8pnLTt+r7HDMrrDZu0ztPbZyUlkC1kGWae\nxWvDqhSzsgJV32in9YGgrDfKoDfFuOfIUz3ASkqgv9Bp07xrkrIu6nW6xWarA0GZb5SsN8WvvWh6\n7HPk6e2PukAOaIvQZZhlXMjrcovNVgeCMt8oWW+KvOfoN2ooY+gJNNWwkXYVadO6XMiLanUgKPuN\n0u9N8c6ddw09R9bI5LWnzeiq3fNjDz2Bpsozmidtml+rq4ZG3YAu1DmyRg03379QefUAkII81XJF\nN6Br0/qcVo8Iqsiv5znHoJEJvX3gSHlG80XSpm1bn9PqQFBFfj3PORjCAqPJ+5kZtyMVY7+fmFod\nCKRq8uvDzkHlDzCaMj8z/Sad27Y+p/WBoA6o/AFGU9ZnJisFNLVyQo/v23/E8U0dpRMIaoK5AGA0\nZXxmslJAxx59VOX7/cTU6qohAO2Wlep5YnF/oYq91CqOGBEAaK1Bk87jjjhSrDhiRACgtULc7D7F\nHYFbMyKo+mbQAOovRKFGihVHUQKBmb1T0pvVucPb3ZIucvcfhTpfikM1ANUou1AjxXVBlaeGzGxG\n0tskzbr7qZJWSNoS8pwpDtUApClEuim0WKmhoyVNmtl+SSslfTfkyVIcqo2LFBgQV4rrgioPBO4+\nb2Z/J+kRSYuSbnT3G0OeM8Wh2jhIgQH1kNq6oBipoedKOk/SSZJOkPRMM/udPsdtNbM5M5tbWFgo\ndM4Uh2rjIAUGYBwxUkO/Lukhd1+QJDO7WtLLJH1y+UHuvkPSDkmanZ31IidMcag2jjalwICmiZnW\njREIHpF0hpmtVCc1dJakudAnTW2oNo62pMCApomd1q08NeTut0v6rKQ71CkdPUrdnj+KaUsKDAgl\n1tYQsdO6UaqG3P1SSZfGOHeTtSUFBoQQs1ceO63bmpXFbdGGFBgQQsyb0cRO67LXEIBWGJb2idkr\nj53WZUQAoPHypH1i9spjp3UJBAAaL0/aJ/YtY2OmdQkEy7A9A9BMedI+sXvlMREIumLX8QIIJ2/a\np63FFkwWd8Wu4wUQTuzJ2LpjRNAVu44XQDhtTvvkQSDoil3HCyCstqZ98iA11MXQEUBbMSLoYugI\noK0IBMswdATQRqSGAKDlCAQA0HKkhkbE6mMATUMgGAGrjwE0EYFgBKPsV87IAUAqCAQjyLv6mJED\ngJQwWTyCrFXGvY+zbxGAfmLdE3kYAsEI8q4+Zt8iAL2WMgXzexflOpQpqEMwIBCMYNP6GV22ea1m\npiZlkmamJnXZ5rVHpHvyjhwAtEedMwXMEYwoz+rj2Hc6AlCuMoo/6pwpIBAEwL5FQHOUVfxR5x2O\nCQSBsG8R0AyjlI0PUudMAYEAAAYoK6VT50wBgQAABigzpVPXTAFVQwAwQBtuWsWIAAAGqHNKpywE\nAgAYoq4pnbKQGgKAlqs8EJjZKWZ217I/T5rZO6puBwCgo/LUkLs/IGmdJJnZCknzkj5XdTsAAB2x\nU0NnSfq2uz8cuR0A0FqxA8EWSVf2+4aZbTWzOTObW1hYqLhZANAe0QKBmR0j6VxJ/9rv++6+w91n\n3X12enq62sYBQIuYu8c5sdl5kt7q7r+Z49gFSUXSR6sk/aDA/4+N9seX+s9A++OK1f4XuvvQnnTM\ndQQXKCMt1CvPDzKImc25+2yR54iJ9seX+s9A++Oqe/ujpIbM7JmSfkPS1THODwA4JMqIwN3/T9Jx\nMc4NADhc7KqhquyI3YCCaH98qf8MtD+uWrc/2mQxAKAe2jIiAABkaFQgMLNzzOwBM3vQzLb1+f6x\nZraz+/3bzWxN9a3MlqP9F5rZwrJ9mt4co51ZzOwKM9tjZvdkfN/M7P3dn++bZvaSqts4SI72n2lm\nTyz7/f9F1W0cxMxONLObzew/zexeM3t7n2Nq+xrkbH9tXwMze4aZ/YeZfaPb/r/qc0w9r0Hu3og/\nklZI+rakn5N0jKRvSPqFnmP+SNKHu19vkbQzdrtHbP+Fkj4Qu60DfoZXSHqJpHsyvv8qSV+QZJLO\nkHR77DaP2P4zJV0Xu50D2n+8pJd0v362pG/1eQ/V9jXI2f7avgbd3+mzul9PSLpd0hk9x9TyGtSk\nEcHpkh509/9y959I+rSk83qOOU/Sx7tff1bSWWZmFbZxkDztrzV3v0XSDwcccp6kT3jHbZKmzOz4\nalo3XI7215q7f8/d7+h+/b+S7pPUu4l+bV+DnO2vre7v9KnuPye6f3onYWt5DWpSIJiR9J1l/35U\nR76Jnj7G3Q9IekL1KWPN035Jem13SP9ZMzuxmqaVJu/PWGcv7Q79v2Bmvxi7MVm6KYf16vRKl0vi\nNRjQfqnGr4GZrTCzuyTtkXSTu2f+/ut0DWpSIGiDf5e0xt1/SdJNOtSzQDXuUGfJ/i9L+kdJ10Ru\nT19m9ixJV0l6h7s/Gbs9oxrS/lq/Bu5+0N3XSXqBpNPN7NTYbcqjSYFgXtLyHvILuo/1PcbMjpb0\nHEmPVdK64Ya2390fc/cfd//5z5JOq6htZcnzGtWWuz+5NPR3989LmjCzVZGbdRgzm1DnIvopd++3\ncr/Wr8Gw9qfwGkiSu++VdLOkc3q+VctrUJMCwdclnWxmJ3V3Nt0i6dqeY66V9Mbu1+dL2uXdWZsa\nGNr+nlzuuerkUFNyraTf61aunCHpCXf/XuxG5WVmP7uUzzWz09X5/ET/EC/ptu2jku5z97/POKy2\nr0Ge9tf5NTCzaTOb6n49qc42Ovf3HFbLa1Bjbl7v7gfM7I8l3aBOBc4V7n6vmf21pDl3v1adN9m/\nmNmD6kwKbonX4sPlbP/bzOxcSQfUaf+F0Rrch5ldqU5Vxyoze1TSpepMmMndPyzp8+pUrTwoaZ+k\ni+K0tL8c7T9f0lvM7ICkRUlb6vAhXmaDpN+VdHc3Ty1JfypptZTEa5Cn/XV+DY6X9HHr3HnxKEmf\ncffrUrgGsbIYAFquSakhAMAYCAQA0HIEAgBoOQIBALQcgQAAWo5AAIzJzF5vZvd1/7wudnuAcVE+\nCozBzJ4t6SuSflWdXSe/KukVyzYdA5LBiAAYzzmSvuDuT3V3yvySpFdGbhMwlsasLAYqtlrSRWb2\nmu6/pyV9P2J7gLExIgDG9yF3X9fdbfIjsRsDjItAAIznER2+j/+Jkh6O1BagECaLgTF098yfU2ej\nNJP0NUnrmSxGipgjAMbg7k+Z2Xsk3dZ96BKCAFLFiAAAWo45AgBoOQIBALQcgQAAWo5AAAAtRyAA\ngJYjEABAyxEIAKDlCAQA0HL/D5ilBalDlPJlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f2843e71710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data = np.load(\"orientation-tuning_data.npz\")\n",
    "θarr = data['θarr']\n",
    "farr = data['farr']\n",
    "\n",
    "plt.scatter(θarr, farr)\n",
    "plt.xlabel(\"θ\")\n",
    "plt.ylabel(\"f(θ)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Defining a model\n",
    "\n",
    "## On the ubiquity of models\n",
    "Any interpretation of data implies an underlying model of how the different variables in that data are related. For example, a typically used measure to show that two variables are related is the *correlation coefficient*, which takes values between -1 and 1 (with 0 corresponding fully uncorrelated variables). However, applying this measure to data implicitly assumes a *linear* relationship between the variables—effectively, computing the correlation is akin to imposing a linear model to data. That is why, despite the plot above clearly showing a relationship between $θ$ and $f(θ)$, the correlation between those two quantities is a minuscule $r = -0.043$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.       , -0.0429521],\n",
       "       [-0.0429521,  1.       ]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.corrcoef(θarr, farr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Thankfully, a model doesn't have a be perfect to be useful. In data analysis, the game is not to find the “perfect“ model, but rather a “good enough” model, that explains the data sufficiently well. We can then compute statistics that tell us how well *that* model fits the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Our toy model\n",
    "When choosing a model, a good principle is to apply Occam's razor: make it as simple as possible, adding only as much complexity as required to explain the feature that interests us.\n",
    "In this case, we need a model that minimally:\n",
    "\n",
    "- captures the orientation corresponding to maximum neural response $θ_0$;\n",
    "- is periodic (rotating the grating by 180° should give the same response);\n",
    "- ensures $f(θ)$ is always positive.\n",
    "\n",
    "A possible model with these characteristics would be\n",
    "$$f(θ) = α \\cos( 2(θ - θ_0) ) + β \\,,$$\n",
    "\n",
    "It is important to stress that this only *a* model: plenty of other models would also appropriately explain our data. They may involve different parameters, which would take on different values, but they should all agree on measurable quantities: **fitting any sufficiently expressive model to the same data will produce similar estimates for the contrast and preferred orientation.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "To translate this model into Python, we define a function `f` which returns the firing rate. You should be able to recognize the components of the equation above in the definition below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def f(θ, α, β, θ0):\n",
    "    return α * np.cos((θ - θ0)*2) + β"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Determining the contrast ratio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The contrast ratio is given by \n",
    "$$c = \\frac{f_{max}}{f_{min}} \\,.$$\n",
    "Since $f_{max}$ and $f_{min}$ respectively occur when $θ = θ_0$ and $θ = θ_0 \\pm \\tfrac{\\pi}{4}$, we find that\n",
    "\\begin{equation}\n",
    "c = \\frac{f(θ_0)}{f(θ_0 + \\tfrac{\\pi}{4})} = \\frac{α + β}{-α + β} \\,. \\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "Thus, to obtain the contrast ratio $c$, we only need to determine $α$ and $β$. This is good news, as our model is linear in both of these parameters, which makes regression easy. For this reason, we will assume that $θ_0$ is known and set it to 90° for the purposes of this exercise. Of course, in a realistic setting $θ_0$ would also need to be determined from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Preferred orientation obtained through a magical oracle\n",
    "θ0 = 0.5*np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Finding $α$ and $β$ using least-squares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "To find the values of $α$ and $β$, we need a way of quantifying “how wrong” our model is. There are many ways of doing this, and perhaps the most common is the squared distance:\n",
    "\\begin{equation}\n",
    "D = \\sum_i (\\hat{f}(θ_i) - f(θ_i))^2 \\,, \\tag{2}\n",
    "\\end{equation}\n",
    "where the $f(θ_i)$ are the measured data points and $\\hat{f}(θ_i)$ are the corresponding estimates made by our model. A model that perfectly predicts the data would have $D = 0$, and any deviation from the data would cause $D$ to become greater than 0. Thus the parameter values which minimize $D$ would be “optimal” in that they produce the model which is closest to perfect – this is called a *least-squares estimate* for those parameters.\n",
    "\n",
    "To find the minimum of a function, we set its derivative to 0. (c.f. your favourite introductory calculus resource.) Here we want to minimize $D$ for both $α$ and $β$, so we set both those derivatives to zero:\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "\\frac{\\partial\\,D(α,β,θ_0)}{\\partial α} &= \\sum_i \\biggl[ 2 \\,\\cos( 2(θ_i - θ_0))^2 \\, α + 2 \\, \\cos( 2(θ_i - θ_0) ) \\, β  - 2 \\, \\cos( 2(θ_i - θ_0) ) f(θ_i) \\biggr] \\\\\n",
    "\\frac{\\partial\\,D(α,β,θ_0)}{\\partial β} &= \\sum_i \\biggl[ 2 \\,\\cos( 2(θ_i - θ_0)) \\, α + 2 \\, β  - 2 f(θ_i) \\biggr] \\,.\n",
    "\\end{aligned} \\tag{3}\n",
    "\\end{equation}\n",
    "We can rewrite this as a system of two equations with two unknowns:\n",
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "0 &= a_{11} α + a_{12} β + a_{13} \\\\\n",
    "0 &= a_{21} α + a_{22} β + a_{23} \\,,\n",
    "\\end{aligned} \\tag{4}\n",
    "\\end{equation}\n",
    "where\n",
    "\\begin{align*}\n",
    "a_{11} &= \\sum_i 2 \\,\\cos( 2(θ_i - θ_0))^2 & a_{12} &= \\sum_i 2 \\, \\cos( 2(θ_i - θ_0) ) & a_{13} &= - \\sum_i 2 \\, \\cos( 2(θ_i - θ_0) ) f(θ_i) \\,; \\\\\n",
    "a_{21} &= \\sum_i 2 \\,\\cos( 2(θ_i - θ_0)) & a_{12} &= \\sum_i 2 & a_{13} &= - \\sum_i 2 f(θ_i) \\,, \\\\\n",
    "\\end{align*}\n",
    "\n",
    "Solving this system for $α$ and $β$ will give us the desired parameters.\n",
    "\n",
    "*Remark* – We should note that choosing to minimize the squared error is equivalent to assuming that the discrepancies between the data and our estimate are due to *additive Gaussian noise*. Thus the type of noise we expect should dictate our choice of error measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 1. Determine the contrast ratio\n",
    "\n",
    "Solve (on paper) the system of equations (4) for $α$ and $β$ in terms of the $a_{ij}$. Then use Python, along with the data loaded at the beginning of this notebook, to compute the least-squares estimates for $α$ and $β$.\n",
    "\n",
    "Finally, use Eq.&nbsp;(1) to compute the contrast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Type your solution in here. Add more code cells with the '+' button in the toolbar.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2a. Perform least-squares using an optimizer\n",
    "\n",
    "Use *SciPy*'s [`minimize` function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize) to directly minimize Eq.&nbsp;(2); it should find the same values for $α$ and $β$ as in Exercise&nbsp;1. (You can leave the `method` parameter to its default value, which will automatically select a suitable optimizer.)\n",
    "\n",
    "You can start from the code suggested below, or go your own way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Find the optimal parameters\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "def D(params):\n",
    "    return sum( _____________________ for θ, fobs in zip(θarr, farr))\n",
    "\n",
    "res = minimize(_____________)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Print them\n",
    "print(\"α = {}\".format(res._____))\n",
    "print(\"β = {}\".format(_________))\n",
    "print(\"contrast: {}\".____________)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 2b. Optimize using the gradient\n",
    "\n",
    "Repeat Exercise&nbsp;2, this time specifying the 'CG' method. Note that this method makes use of the derivative (which we calculed in Eq.&nbsp;(3)), so you will need to define an function for the gradient. (In the *SciPy* documentation, the derivative is referred to as the “gradient” or the “jacobian”.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Solution to Exercise 2b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 3a. Make your own data\n",
    "\n",
    "Produce your own dataset, using the same model but with different parameters. Plot the frequency response curve ($f(θ)$ vs $θ$), and see confirm that the code you wrote for exercises 1 and 2 is able to recover the correct parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Solution to Exercise 3a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## 3b. Quantifying error\n",
    "\n",
    "Using your answer from exercise&nbsp;3a, generate multiple sets of 100 data points for the same parameters. What is the variance in the estimate of $α$ and $β$ ? Plot a histogram ([`matplotlib.pyplot.hist`](https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist) of the estimated parameters – are they centered around the true values of $α$ and $β$ ?\n",
    "\n",
    "Compute the variance on the estimate of $α$ for sample sizes ranging between 5 and 5000; polt the result as a function of sample size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Solution to Exercise 3b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true,
    "heading_collapsed": true
   },
   "source": [
    "# Resources\n",
    "\n",
    "The following functions were used to generate data for this notebook and are not part of the exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "class NeuronModel:\n",
    "    \n",
    "    def __init__(self, θ0, α, β, σ):\n",
    "        self.α = α\n",
    "        self.β = β\n",
    "        self.σ = σ\n",
    "        self.θ0 = θ0\n",
    "        \n",
    "    def f(self, θ):\n",
    "        return np.clip(self.α * np.cos((θ - self.θ0)*2) + self.β + self.σ * np.random.normal(0, 1),\n",
    "                       0, np.inf)\n",
    "        return retval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "neuron = NeuronModel(np.pi/2, 1.8, 10, 0.6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true,
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEPCAYAAABP1MOPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHEtJREFUeJzt3X2wXVV5x/HfYwh4Yy0XydWRKzTY0uBLpgncWjXiaHwB\nXyoxWoGxjlVqRmurMlM0jFb8wylx6NTq1GmNI2MdHMQWiFRaXhQqIxXqDQQTysuAgHJRc5UEq7nK\nJTz945wTTs49++x9ztlrv63vZ8bh5mTffZYr5+xnr2c9a21zdwEA4vWUshsAACgXgQAAIkcgAIDI\nEQgAIHIEAgCIHIEAACJHIACAyBEIACByBAIAiNxhZTcgi5UrV/qqVavKbgYA1MqOHTt+5u5TacfV\nIhCsWrVKs7OzZTcDAGrFzB7MchypIQCIHIEAACJHIACAyBEIACByBAIAiByBAAAiV4vyUSAP22+b\n04XX3K2H9y3omMkJnXvqam1cN112s4DSEQgQhe23zem8y3dpYfGAJGlu34LOu3yXJBEMED1SQ4jC\nhdfcfTAIdCwsHtCF19xdUouA6iAQIAoP71sY6nUgJqSGEIVjJic01+eif8zkRAmtaWHOAlXBiABR\nOPfU1ZpYvuyQ1yaWL9O5p64upT2dOYu5fQtyPTlnsf22uVLag7gxIkCtjHoX3TmmKnfgg+YsxmkT\nowyMgkCA2hhU+SOlX+Q3rpuuzEUxxJwFlVEYFYEAtZF0F/2JK+/Qbx5/olYXwBBzFqFGGWg+5ghQ\nG0l3y/sWFitXGrr9tjmt33q9jt9yldZvvX5J7j/EnAWVURgVgQC1MezdclkXwCwTwRvXTeuCTWs0\nPTkhkzQ9OaELNq0Z6849qX/KrIxCPQQLBGZ2kZntMbPdXa9daGZ3mdn3zewKM5sM9f5onqS76KNW\nLO97fFkXwKyL1zaum9ZNWzbo/q1v0E1bNoydvqlaZRTqI+QcwZck/aOkL3e9dp2k89z9cTP7lKTz\nJH0kYBvQIEmVP5IOmSSVirsA9qvSKStFU7XKKNRHsEDg7jea2aqe167t+uPNkt4a6v3RTIMqf4q+\nACZV6UyuWK69+xeXHF/ECKVKlVGojzKrht4t6dKkvzSzzZI2S9Jxxx1XVJtQU2VcAJNSQEcc9hRN\nLF9WyggFGEUpk8Vm9lFJj0v6StIx7r7N3WfcfWZqaqq4xqH20ip28pKU6nl0YTH3iWAgpMJHBGb2\nZ5LeKOlV7u5Fvz+archFVYPWApCiQZ0UOiIws9MkfVjSm9x9f5HvjTgUud00VTpoimAjAjO7RNIr\nJK00s4ckna9WldARkq4zM0m62d3fG6oNiE+RFTtU6aApQlYNndXn5S+Gej9AKn67aVJAaAJWFqNR\n6pquKWqCG+iHTefQKHVM17BrKMpGIEDj1C1dk3XXUJ41gFAIBEDJskxwM2pASMwRACXLsmtokWWx\niA+BAIVjYvRQWSa4edYAQiI1hELVKcVRVE4+ywR30WWxiAuBAIWqy+MUiw5YaRPc5566urStttF8\nBAIUqi4pjqoFrEGjBqqJMC4CAQpVlxRHFQNWv1FDnVJtqC4mi1Gouqz8rcvzf6kmQh4IBChUiIe2\nh1CXgFXFkQvqh9QQCleFlb/defUjJ5bLTNq3f3FJjr3qufe6pNpQbQQCRKc3r75v4cnnC/fm2Kt2\n4e9FNRHyQCBAdPrl1btVsZw1SdaRC5VFGIRAgOhkyZ/XKceeNnKhsghpmCxGbY26VUWW/HmTcuxU\nFiENIwLU0jh3uf3y6t3ScuxJaZaqpl+oLEIaAgHGEuril3becVb+9ubVB1UN9WtXvwA0++AjumzH\nXCXTL1QWIQ2BACMLlXvOct5x73JHrQhKCkCX3PIjHXBf8vook855B1cqi5CGOQKMLFTuOct5Q638\nTZt3SAo0vUEg7fhB73/e5bs0t29BrieD4DhbdddlER/Kw4gAIwuVe85y3hB3uVlGIklplmVmfYPB\nsIEp1GZ3dVgTgfIwIsDIQt2VZzlviLvcLCORpK0nzvqjYxO3pBimuomJXZSBEQFGFir3nPW8ed/l\nZrkID1rANfM7z1jyuqSh5lGY2EUZCAQYWaj9eMra5yfrRTgpAPV7ff3W64dK9TCxizIQCDCWULnn\nMnLaIS7Cw6Z66rLZHZqFQAC0hbgIj5LqYWIXRSMQAF3yvgiT6kEdEAiAgEj1oA4IBEBgpHpQdcHW\nEZjZRWa2x8x2d732J2Z2h5k9YWYzod4bAJBdyAVlX5J0Ws9ruyVtknRjwPcFAAwhWGrI3W80s1U9\nr90pSWYW6m3RQFXd3rmu6E/0quwcgZltlrRZko477riSW4NQ0i5KPF0rX/Qn+qnsXkPuvs3dZ9x9\nZmpqquzmIIAsO23ydK180Z/op7IjAjRflp02m7YJW9lpmab1J/JR2REBmi/LRSnUDqdlCPGsgWE1\nqT+Rn5Dlo5dI+q6k1Wb2kJmdbWZvNrOHJL1E0lVmdk2o90f1ZbkoJW37XMeVuXmmZYbZ2rpbk/oT\n+QlZNXRWwl9dEeo9Mb4iUxdZtl9o0srcvNIy40z4Nqk/kR/mCHBQ0RUlWS9KTVmZm9ezBsZ9illT\n+hP5IRDgoFCPSRykLhelPEZKeW1A17QJ37In0EEgQJemXWDyktdIKa+0TJOeYsa6hmogEOCgJl1g\n8pTnSCmPEVCTtrYuYxSKpSgfxUFUlPRXtZHSxnXTumDTGk1PTsgkTU9O6IJNa2p54axa38aKEQEO\noqKkvyqOlOoyt5Kmin0bIwIBDtGUC0yempSKyVOVJtAxHgIBkKKuI6WQ1ThVm0DHeMzdy25DqpmZ\nGZ+dnS27GVHLclGhDLA6ei/UUutOO6+5hPVbr++b0pmenNBNWzaMfX7kw8x2uHvqQ8CYLEaqLHvk\n5LGPzqjbJmCp0LuMJk3mzu1b4N+uhggESJXlojLuhacKG7I1SehqnEGTufzb1Q+BAKmy3P2Ne+Fh\nn/x8hd5ltF+pcTf+7eqFQIBUWe7+JlcsH/p3u1FPnq/Qa0K61zIk4d+uPggESJXl7s9dY1142Cc/\nX0UsOtu4blo3bdmQGAz4t6sPqoaQSaciqF+liCSZpE+fsfZg1dCRE8tlJu3bv5ipgih0lUtd1aES\ni3+76qJqCLnKcvfXOebTZ6zVbx5/Qnv3L2ae+G3Stgl5qcsEOv929ceIoMbKuFvMcvdHjXk+6EeM\nixFBw5V1t5jl7o8a83wwgY6isMVETZW5fW/afkRJG4lJ7Dc/jKI3ZAsxwqzDHAcYEdRWle8WqTHP\nR5HbgocYYY57TlaaF4dAUFNVLrekxjwfRU7ChljQN8456zJR3hSkhmqq6tv3dtJHSROeVQhYdVDU\ntuAhRpjjnJMnlxWLEUFN1aVkr+ynnpFeyCbECHOcc1Y59dlEjAhqrA4PkSlzv3kejJ5diBHmOOfk\nyWXFYh0BGos6/OF0V/gMuzI8yzmHOQ+rlfORdR0BIwI0FumF4XRGmHmOpEYdtfLksmIRCNBYpBdG\nU5WJ2jqkPpuCyWI0VtkT1XXFSCo+jAjQWKQX0vXL4TOSik+wQGBmF0l6o6Q97v7C9mvPkHSppFWS\nHpD0NnffG6oNAOmFZElzAW85eVqX7ZgbuYKIbSXqJ2Rq6EuSTut5bYukb7n7CZK+1f4zhkBdPPKS\nNBdww13zI69RYUVwPQUbEbj7jWa2qufl0yW9ov3zv0j6L0kfCdWGpqEuHnkaNBcw6kiqKhPNGE7R\nk8XPcvcft3/+iaRnFfz+tcYD3pGnEKuJmWiup9Kqhry1ki1xNZuZbTazWTObnZ+fL7Bl1cWXDHkK\nUVVV5c0QkazoQPBTM3u2JLX/uyfpQHff5u4z7j4zNTVVWAOrjC8Z8hRivypKduup6PLRKyW9U9LW\n9n+/XvD711rVdxxF/eRdVUXJbj1lCgRmNiPpFEnHSFqQtFvSdYNKP83sErUmhlea2UOSzlcrAHzN\nzM6W9KCkt43V+siE+pJR7oc8UbJbPwM3nTOzd0n6K0n3S9qhVirnqZJ+X9J6tQLC37j7D0M2kk3n\nwmFzL6C58tp0boWk9e7edzbSzNZKOkFS0EAQo6Lu0in3QxkYhVbLwEDg7p9L+fud+TYHUrj1Av2+\nfFQioWish6megVVDZrbCzM43s++b2SNm9vP2zx8zM0pVAgmxXiBpxefkiuV9j6cSCaGwHqZ60spH\nv6rWwq+XSnrQ3Y+W9Lr2310asmExC3GXnvTlc1diuR/bWSAERqHVkxYITnD3z7v7L9Ve/OXuc+7+\nSUnPC966SBW54vPRhcW+teSS2DMGQbAepnrSAsG9Zna2mR3eecHMpszsI5IYxwVS9IrPjeumddOW\nDbp/6xt005YN2rhumuE7gmHRWfWkBYIzJT1H0vckHW9mP5d0vaTDJZ0RuG3RqsKKT4bvCGXQ55t0\nZDl4eH1EhinZ48HvKBprWvKXyzoCM/uYpM8lrSA2sw2SVrj7N0ZrJoo0zIpPtrNA0VjTUp60BWW7\nJH3DzH4t6VZJ82qtLD5B0lpJ35T0t0FbiFKwZwzyljYiJR1ZnrRA8FZ3X29mH1Zre4lnS/qFpIsl\nbU5acYxmYM8YSPmsAs6yiGzQs5JZiRxW2mTxyWZ2jKS3q7Vz6OclfVmtyWNqvYCGy+vRk1mq0JIK\nGl554hSlzIGlBYJ/VuvZwidKmu363472fwE0WF5lxFnSPknVRDfcNU8pc2Bpew19VtJnzeyf3P19\nBbUJQEXklbcflPbp1i8dec6l/bc0Y+4gP5meUEYQAOKU1yrgcRaRsRI5vNKeWQyg+vJaBTzOIklW\nIodX9KMqAdRInmXEo1ahUcocHiuLAaChsq4sJjUEAJEjEABA5KKYIxhnVWLdVzTWvf2oFj5PzdT4\nQDDO81Hr/mzVurcf1cLnqbkaHwjG2dGw7rshpq0K5c4Ow6j79yGkuo+UGh8IxlkZWffdEJPa2bmT\n484Ow6j796Ej74t2E0ZKjZ8sHmdVYt1XNCa1c5kZe7dgaHX/Pkj5baLXrQmPdW18IBhnVWLdVzQm\ntf9AwtqRut3ZoVh1/z5IYS7aTRgpNT4QjLO0PcSzg4uU1P7pBtzZoXh1/z5IYS7aTRgpNX6OQBrv\nASt1fzhLUvt5DCVGUffvQ9ZdUIfRhMe6Nn5EgKWacGcHjCJEeqsJ3yf2GgIQlbqXeg4j615DpaSG\nzOyDkt4jySR9wd3/oYx2AKi/YS/sdU9vhVB4IDCzF6oVBF4k6TFJV5vZN9z93qLbAqDemlDDXwVl\nzBE8T9It7r7f3R+X9G1Jm0poB4Caa0INfxWUEQh2SzrFzI42sxWSXi/p2BLaAaDmmlDDXwWFBwJ3\nv1PSpyRdK+lqSTslHeg9zsw2m9msmc3Oz88X3EoAddCEGv4qKKV81N2/6O4nu/vLJe2VdE+fY7a5\n+4y7z0xNTRXfSACV14TVzlVQVtXQM919j5kdp9b8wIvLaAeAeuN5xvkoa2XxZWZ2tKRFSe93930l\ntQNAzVEOOr5SAoG7n1LG+wIAlmKLCQCIHIEAACIXxe6jAOIQ0z5CeYo6EPChAZqD7SZGF21qKMQj\n6wCUh+0mRhdtIOBDAzQL202MLtrU0CgfmrxSSaSkgPyFePpYLKIdEQy7R0leqSRSUkAYRWw3sf22\nOa3fer2O33KV1m+9vjHf22hHBMM+ZzQtlZT1Dj/pPB+6dKcuvOZuRgfAiEJvN1HEZHRZ2YJoA8Gw\nH5qklFHnw5D1wzEo9USVAzCekNtNDLoZTHrPYS7sZVY9RRsIpOE+NEn5x2VmQ304ks6T5XcBlGfY\necVBF3Zp6U3oKIEmL9HOEQwrKf94wL3v8Ukfjn7nyfq7AMqTNH/oUt/5gqQL+yeuvKPvPGHSDWIR\n1wMCQUYb103rgk1rND05IZM0PTlx8M/9JH04us+TZJQqh6ZOYgFVMegmrl/RR9IFfN/CYt8Ascys\n7/FFVD1FnRoaVlIqqXfSuSMpx9c5T+/QURqtyoEVlUB43fOK/e7ee9M4aWngXgfcNbF82djXg1Ew\nIhhT2h3+oEVqSaOM3ot32t0+i+OAYmxcN62btmxQ/3v3Q0cBSenko1Ys7/u73VmGQdeDEBgR5KBz\nh3/8lqvUb8ZgUI4vbcI6y90+KyqBpUKWYmZZvJZUmSgtzSJ07vzLesgOgaAtjw/NoA/HqOfPUknA\nikrgUKHTpVnXIQ26sFdpdwECgfL70CR9OF554tTI589ytz/s4jig6UKXYo67eK1qj9ckECi/D03S\nhyPr+fuNGsYZglbpgwbkJcvouoh0adUu5uMgECjfD02/D8c5l+5MPX/SqOQtJ0/rsh1zYw1BgabI\nOnonXTocqoY0/AZ0Ic6fNGq44a750ioJgKrJWiE37gZ0sa3LYUSg8Dn2LOcfNCrhbh9oyTp6Hydd\nGuO6HAKBwufYs5yfoSyQbpjvyag3UGXu+VMWAkFb6LvutPNT+QOky/t70m/iOcZ1OQSCiqDyB0iX\n5/ckKQU0uWK59u5fXHJ8k0fnBIIKYS4ASJfX9yQpBXTEYU8pbc+fslA1BCBKSameRxcWx6rUq2PF\nESMCAFEaNPE86qijrhVHjAgARCnEw+7ruhNwdCOCsh4ODaBaQhRo1LXiqJRAYGbnSPpztR7ktUvS\nu9z916Hft67DNgBh5F2gUdf1QIWnhsxsWtIHJM24+wslLZN0ZhHvXddhG4B6CJFuKkJZqaHDJE2Y\n2aKkFZIeLuJN6zpsGxbpL6AcdV0PVHggcPc5M/s7ST+UtCDpWne/toj3ruuwbRikv4By1XE9UBmp\noaMknS7peEnHSHqamf1pn+M2m9msmc3Oz8/n8t51HbYNg/QXgGGVkRp6taT73X1ekszsckkvlXRx\n90Huvk3SNkmamZnp9yjgodV12DaMWNJfQFNUIZVbRiD4oaQXm9kKtVJDr5I0W9Sb13HYNowY0l9A\nU1QllVt4asjdb5H0b5JuVat09Clq3/ljfDGkv4AQytgaoiqp3FKqhtz9fEnnl/HeTRdD+gvIW1l3\n5lVJ5Ua3sjgGTU9/AXkr62E0VUnlstcQgMZLS/uUdWdelVQuIwIAjZYl7VPWnXlVUrkEAgCNliXt\nU+ajYquQyiUQ9FGFul4A+ciS9qnKnXlZCAQ9qlLXCyAfWdM+VbgzLwuTxT2qUtcLIB9VmZCtMkYE\nPapS1wsgH7GnfbIgEPSoSl0vgPzEnPbJgtRQD4aRAGLDiKAHw0gAsSEQ9MEwEkBMSA0BQOQIBAAQ\nOVJDI2L1MYCmIBCMgNXHAJqEQDCCrHuXM2oAUAcEghFkWX3MqAFAXTBZPIKkVcbdr7NnEYBuZTwT\nOSsCwQiyrD5mzyIAHZ0Mwdy+BbmezBBUJRgQCEawcd20Lti0RtOTEzJJ05MTumDTmkNSPllGDQDi\nUPUMAXMEI0pbfVzmE48A5CePoo+qZwgIBIGwZxFQf3kVfVR9V2MCQUDsWQTUW9ZS8TRVzxAQCAAg\nQV4pnapnCAgEAJAgz5ROlTMEVA0BQIJYHlTFiAAAElQ9pZMXAgEADFDllE5eSA0BQOQKDwRmttrM\ndnb97xdm9qGi2wEAaCk8NeTud0taK0lmtkzSnKQrim4HAKCl7NTQqyTd5+4PltwOAIhW2YHgTEmX\n9PsLM9tsZrNmNjs/P19wswAgHqUFAjM7XNKbJP1rv793923uPuPuM1NTU8U2DgAiYu5ezhubnS7p\n/e7+2gzHzkvKI320UtLPcjhPDOirbOin7Oir7PLqq99x99Q76TLXEZylhLRQryz/R7Iws1l3n8nj\nXE1HX2VDP2VHX2VXdF+Vkhoys6dJeo2ky8t4fwDAk0oZEbj7ryQdXcZ7AwAOVXbVUNG2ld2AGqGv\nsqGfsqOvsiu0r0qbLAYAVENsIwIAQA8CAQBEjkAAAJFrXCAws9PM7G4zu9fMtpTdnjKl9YWZHWtm\nN5jZ/5rZHWb2wa6/e8DMdrV3iJ0ttuXlyvIZirV/MnymEncXjrXPkpjZRWa2x8x2l90WuXtj/idp\nmaT7JD1X0uGSbpf0/LLbVdW+kPRsSSe1f366pHs6x0h6QNLKsv9/VLHfYu2fYb9f7eN/otbq1ij7\nLKU/Xy7pJEm7y25L00YEL5J0r7v/wN0fk/RVSad3H2Bmv2tm8+27k51m9oiZ3Wdmv11Ki8NJ7Qt3\n/7G739r++f8k3Skp86OYGtqXqf2WVQP7Z9i+GXp34Qb2WSJ3v1HSI0l/X2RfNC0QTEv6UdefH1LP\nhc3d75P0HUnvcPe1kr4vaaO7/6KwVhYjtS+6mdkqSesk3dJ+ySVda2Y7zGxzv99paF9m7bcY+2eo\nz5SW7i4cY5+NrMi+iPWZxS+Q1MnLPU/S3SW2pXRm9luSLpP0oa4P2cvcfc7MninpOjO7q30H0yvW\nvqR/BujaXfi8rpfps+EV0hdNGxHMSTq268/Pab92kJlNSHqqu+81s2Ml/aw9zG2a1L6QJDNbrlYQ\n+Iq7H9z7yd3n2v/do9YT5F7U53eb2JeZ+i3S/snUN22vk3Sru/+080KkfTayIvuiaYHge5JOMLPj\n23ckZ0q6sueY56uVC5daEfZONVNqX5iZSfqipDvd/e+7Xn+amT2987Ok1+rJu5JuTezLLP0Wa/9k\n+X51HLK7cMR9No7C+qJRgcDdH5f0l5KuUavTvubud/Qc1j3UWpB0kpmdWFwrizGoL8zsP8zsGEnr\nJb1D0oaucr/XS3qWpO+Y2e2S/kfSVe5+dZ+3aVxfZuy3KPsnY98k7S4cZZ8NYmaXSPqupNVm9pCZ\nnd1zSGF9wV5DABC5Ro0IAADDIxAAQOQIBAAQOQIBAESOQAAAkSMQAEDkCAQAEDkCARrLzJaZ2Wes\n9ayFXWb23D7HTJjZt81s2Rjv8wkz++v2z/89Tpv7nPtwM7vRzGLdFwwFIBCgyc6T9AN3f4Gkz0r6\niz7HvFvS5e5+oPtFaxn6++HuLx2ppcnne0zStySdked5gW4EAjRSe5uDN7v7Z9ov3S/p9/oc+nZJ\nX2//ziprPX3ry2ot7T/WzLa3t02+o3vrZDP7qJndY2bfkbS66/Vfdv285Hfb73GnmX2h/fq17c3F\nOvvxXGVmt5vZbjPrXPy3t9sJBMFwE031arUu5Dvbf36GpG92H9DeOO257v5A18snSHqnu9/cPubd\n7v5I+2L9PTO7TNIqtTZcW6vWd+hWSTv6tKHf73be4yx3f4+ZfU3SWyRdLOk0SQ+7+xva731k+/jd\nkv5wxH4AUjEiQFOtlfRxd1/bfqjHtZJ29hyzUtK+ntce7ASBtg+0N0q7Wa0tmE+QdIqkK9x9f/v5\nDUk7cPb7XUm63907bdmhVmCRpF2SXmNmnzKzU9z9UUlqp60e6+zeCeSNQICmOkrSfklqT7S+VtK/\n9xyzIOmpPa/9qvODmb1CrZHFS9z9DyTd1uf4vlJ+9zddhx5Qe2Tu7veo9QzbXZI+aWYf7zruCEm/\nzvLewLAIBGiqeyS9uP3zOWpte3x/9wHuvlfSMjNLurgfKWmvu+9vb//bOd+Nkja2K46eLumPh/jd\nRO1tnPe7+8WSLlQrKMjMjlbroSSLaecARsEcAZrqEkn/aWb3qrXne99n5KqVMnqZeuYP2q6W9F4z\nu1OtRwTeLEnufquZXSrpdkl71HpgS6bfTbFG0oVm9oSkRUnva7/+SklXZfh9YCQ8jwBRM7OTJJ3j\n7u8ouy1JzOxySVvaqSMgd6SGEDV3v1XSDeMsKAupXdm0nSCAkBgRAEDkGBEAQOQIBAAQOQIBAESO\nQAAAkSMQAEDkCAQAEDkCAQBE7v8BOpD2y/yL2vwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe54ae50d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#θarr = np.random.uniform(0, np.pi, 15)\n",
    "θarr = np.linspace(0, np.pi, 100)\n",
    "farr = np.array([neuron.f(θ) for θ in θarr])\n",
    "plt.scatter(θarr/np.pi, farr)\n",
    "plt.xlabel(\"$θ$ (radians)\")\n",
    "plt.ylabel(\"f(θ)\")\n",
    "# Magic incantation for radians from http://stackoverflow.com/a/40642200\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_formatter(plt.FormatStrFormatter('%g $\\pi$'))\n",
    "ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(base=0.25))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "hidden": true
   },
   "outputs": [],
   "source": [
    "np.savez(\"regression_orientation-tuning_data.npz\", θarr=θarr, farr=farr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true,
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (mackelab)",
   "language": "python",
   "name": "mackelab"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}