classes/09d-HandsOn-Scipy-Convolution.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:eb98df4455e903a248b7b27a98649f9f0a743294233d0b561c0e240ad420754a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Hands-on: Numerical Python -- Some SciPy -- Convolutions Excercise"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Disclaimer/Appraisal**: This excercise is a based on [Practical NeuroImaging](http://github.com/practical-neuroimaging/pna2015) course tought at Berkeley by Matthew Brett (nibabel).  His copyright and distributed under  Creative Common Attribution (CC-BY 2.0) license.  The BOLD data ``files/ds114_sub009_t2r1.nii.gz`` comes from the OpenFMRI project dataset ``ds114``: https://openfmri.org/dataset/ds000114 . It is a recompressed and renamed copy of ``sub009/BOLD/task002_run001/bold.nii.gz``.\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Scipy"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "*Scipy* is a large library of scientific routines that builds on top of numpy. \n",
      "\n",
      "You can think of numpy as being a subset of MATLAB, and numpy + scipy as being as being roughly equivalent to MATLAB plus the MATLAB toolboxes."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Scipy has many sub-packages, for doing things like reading MATLAB `.mat` files (`scipy.io`) or working with sparse matrices (`scipy.sparse`). We are going to be using the functions and objects for working with statistical distributions in `scipy.stats`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`scipy.stats` contains objects for working with many different distributions.  We are going to be working with `scipy.stats.gamma`, which implements the [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In particular we are interested in the [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) (PDF) of the gamma distribution.\n",
      "\n",
      "Because this is a function, we need to pass it an array of values at which it will evaluate.\n",
      "\n",
      "We can also pass various parameters which change the shape, location and width of the gamma PDF.  The most important is the first parameter (after the input array) known as the [shape parameter](https://en.wikipedia.org/wiki/Shape_parameter) ($k$ in the [wikipedia page on gamma distributions](https://en.wikipedia.org/wiki/Gamma_distribution))."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we chose some x values at which to sample from the gamma PDF:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.arange(0, 25, 0.1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we plot the gamma PDF for shape values of 2, 4, 8, 12."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x, gamma.pdf(x, 2), label='k=2')\n",
      "plt.plot(x, gamma.pdf(x, 4), label='k=4')\n",
      "plt.plot(x, gamma.pdf(x, 8), label='k=8')\n",
      "plt.plot(x, gamma.pdf(x, 12), label='k=12')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will see more from the statistics distributions next week."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Correlation coefficient"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "By \"correlation coefficient\", we mean the [Pearson product moment correlation coefficient](https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient).\n",
      "\n",
      "Let's say I have two arrays, `x` and `y`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = np.array([ 20.57,  18.33,  20.8 ,  18.77,  18.46,  21.09,  19.96,  20.61,\n",
      "               18.73,  19.6 ,  18.15,  19.7 ,  20.36,  20.39,  19.67,  20.73])\n",
      "y = np.array([ 28.98,  30.13,  29.64,  29.87,  28.78,  33.48,  31.36,  31.04,\n",
      "               30.43,  31.69,  27.08,  29.25,  29.5 ,  30.04,  29.74,  30.06])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(x, y, 'o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The numpy routine `np.corrcoef` gives me the correlation *matrix* between `x` and `y`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "corr_matrix = np.corrcoef(x, y)\n",
      "corr_matrix"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "`corr_matrix[0, 0]` is the correlation of `x` with itself, `corr_matrix[1, 0]` (and `corr_matrix[0, 1]`) is the correlation of `x` with `y`."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Be careful - `np.correlate` is something else entirely.\n",
      "\n",
      "**Question: ** what is it?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## dot products"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "If I have two vectors $\\mathbf{a}$ with elements $a_0, a_1, ... a_{n-1}$, and $\\mathbf{b}$ with elements $b_0, b_1, ... b_{n-1}$ then the [dot product](https://en.wikipedia.org/wiki/Dot_product) is defined as:\n",
      "\n",
      "$$\n",
      "\\mathbf{a}\\cdot \\mathbf{b} = \\sum_{i=0}^{n-1} a_ib_i = a_0b_0 + a_1b_1 + \\cdots + a_{n-1}b_{n-1}\n",
      "$$\n",
      "\n",
      "In code:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = np.arange(5)\n",
      "b = np.arange(10, 15)\n",
      "np.dot(a, b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The same thing as\n",
      "np.sum(a * b)  # Elementwise multiplication"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Matrix multiplication operates by taking dot products of the rows of the first array (matrix) with the columns of the second.\n",
      "\n",
      "Let's say I have a matrix $\\mathbf{X}$, and $X_{i,:}$ is row $i$ in $\\mathbf{X}$.  I have a matrix $\\mathbf{Y}$, and $Y_{:,j}$ is column $j$ in $\\mathbf{Y}$.  The output matrix $\\mathbf{Z} = \\mathbf{X} \\mathbf{Y}$ has entry $Z_{i,j} = X_{i,:} \\cdot Y_{:, j}$.   We will see this often over the next few weeks."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X = np.array([[0, 1, 2], [3, 4, 5]])\n",
      "X"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Y = np.array([[7, 8], [9, 10], [11, 12]])\n",
      "Y"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X.dot(Y)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[0, :].dot(Y[:, 0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "X[1, :].dot(Y[:, 0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Convolution and the tail values"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's say I have a neural vector with a couple of spikes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "times = np.arange(0, 60, 0.5)  # samples every 0.5 seconds\n",
      "neural_vector = np.zeros(times.shape)\n",
      "neural_vector[10] = 1  # At 5 seconds\n",
      "neural_vector[20] = 1  # At 10 seconds\n",
      "plt.plot(times, neural_vector)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Then I have my HRF function, sampled every half second, to match:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def hrf(t):\n",
      "    \"A hemodynamic response function\"\n",
      "    values = t ** 8.6 * np.exp(-t / 0.547)\n",
      "    # Scale max to 1\n",
      "    return values / np.max(values)\n",
      "\n",
      "hrf_times = np.arange(0, 20, 0.5)\n",
      "hrf_samples = hrf(hrf_times)\n",
      "plt.plot(hrf_times, hrf_samples)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The input neural vector is length 120, and the HRF vector is length 40:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(len(neural_vector))\n",
      "print(len(hrf_samples))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For reasons that may be familiar to you now, when we convolve the neural vector with the hrf signal, we get an output that is length 120 + 40 - 1 = 159:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hemodynamic_prediction = np.convolve(neural_vector, hrf_samples)\n",
      "len(hemodynamic_prediction)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is because of the HRF falling off the end of the input vector.  The value at index 120 in the new vector refers to time 60, and value 121 refers to time 60.5 seconds.  To retain only the values in the new hemodynamic vector that refer to times up to (not including) 60s, we can just drop the last `len(hrf_signal) - 1` values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hemodynamic_for_60s = hemodynamic_prediction[:len(neural_vector)]\n",
      "plt.plot(times, neural_vector, label='neural vector')\n",
      "plt.plot(times, hemodynamic_for_60s, label='hemodynamic prediction')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Making our own hemodynamic response function"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The file `files/mt_hrf_estimates.txt` contains the estimated FMRI signal from voxels in the MT motion area at 0, 2, 4, ..., 28 seconds after a brief moving visual stimulus (see: http://nipy.org/nitime/examples/event_related_fmri.html).\n",
      "\n",
      "Here are the first four rows.  The numbers are in exponential floating point format:\n",
      "\n",
      "```\n",
      "1.394086932967517900e-01\n",
      "3.938585701431884800e-01\n",
      "5.012927230566770476e-01\n",
      "5.676763716149294536e-01\n",
      "```\n",
      "\n",
      "Read the values from the file into an array `mt_hrf_estimates`.  Make a new array `mt_hrf_times` with the times of acquisition (0, 2, ...).  Plot them together to see the HRF estimate at these times:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the estimated values from the text file into an array\n",
      "# Make an array of corresponding times\n",
      "# Plot signal by time"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We want to make a hemodynamic response function that matches this shape.\n",
      "\n",
      "Our function will accept an array that gives the times we want to calculate the HRF for, and returns the values of the HRF for those times.  We will assume that the true HRF starts at zero, and gets to zero sometime before 35 seconds.\n",
      "\n",
      "Like SPM, I'm going to try using the sum of two [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution) probability density functions."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import the gamma density function\n",
      "from scipy.stats import gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here's my first shot:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - my attempt - you can do better than this\n",
      "def not_great_hrf(times):\n",
      "    \"\"\" Return values for HRF at given times \"\"\"\n",
      "    # Gamma pdf for the peak\n",
      "    peak_values = gamma.pdf(times, 6)\n",
      "    # Gamma pdf for the undershoot\n",
      "    undershoot_values = gamma.pdf(times, 12)\n",
      "    # Combine them\n",
      "    values = peak_values - 0.35 * undershoot_values\n",
      "    # Scale max to 0.6\n",
      "    return values / np.max(values) * 0.6"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - plot the data against our estimate\n",
      "plt.plot(mt_hrf_times, not_great_hrf(mt_hrf_times), label='not_great_hrf')\n",
      "plt.plot(mt_hrf_times, mt_hrf_estimates, label='mt_hrf_estimates')\n",
      "plt.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now see if you can make a better function by playing with the Gamma distribution pdf parameter, and the mix of the two gamma distribution functions.  Call your function `mt_hrf`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Your \"mt_hrf\" function here"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot your function against the mt_hrf_estimates data to test"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For extra points - other than looking at these plots, how would you convince me your function is better than mine?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Getting a predicted neural time course"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to be analyzing the data for the 4D image `files/ds114_sub009_t2r1.nii.gz`.\n",
      "\n",
      "Load the image into an image object with nibabel, and get the image data array.  Print the shape."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Load the image with nibabel, get the image data array, print the data shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Last week we read in the stimulus file for this run to make an on - off timeseries.\n",
      "\n",
      "The stimulus file is `files/ds114_sub009_t2r1_cond.txt`.\n",
      "\n",
      "Here's a version of the same thing we did last week, as a function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - read in stimulus file, return neural prediction\n",
      "def events2neural(task_fname, tr, n_trs):\n",
      "    \"\"\" Return predicted neural time course from event file\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    task_fname : str\n",
      "        Filename of event file\n",
      "    tr : float\n",
      "        TR in seconds\n",
      "    n_trs : int\n",
      "        Number of TRs in functional run\n",
      "        \n",
      "    Returns\n",
      "    -------\n",
      "    time_course : array shape (n_trs,)\n",
      "        Predicted neural time course, one value per TR\n",
      "    \"\"\"\n",
      "    task = np.loadtxt(task_fname)\n",
      "    if task.ndim != 2 or task.shape[1] != 3:\n",
      "        raise ValueError(\"Is {0} really a task file?\", task_fname)\n",
      "    task[:, :2] = task[:, :2] / tr\n",
      "    time_course = np.zeros(n_trs)\n",
      "    for onset, duration, amplitude in task:\n",
      "        time_course[onset:onset + duration] = amplitude\n",
      "    return time_course"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use this function to read `files/ds114_sub009_t2r1_cond.txt` and return a predicted neural time course.\n",
      "\n",
      "The TR for this run is 2.5.  You know the number of TRs from the image data shape above."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Read the stimulus data file and return a predicted neural time course.\n",
      "# Plot the predicted neural time course."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We discovered last week that the first volume in this run was bad.  Use your slicing skills to make a new array called `data_no_0` that is the data array without the first volume:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make new array excluding the first volume\n",
      "# data_no_0 = ?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Our neural prediction time series currently has one value per volume, including the first volume.  To match the data,\n",
      "make a new neural prediction variable that does not include the first value of the time series.  Call this new variable `neural_prediction_no_0`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Knock the first element off the neural prediction time series\n",
      "# neural_prediction_no_0 = ?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For now, we're going to play with data for a single voxel.\n",
      "\n",
      "Last week, we subtracted the rest from the task scans, something like this:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - subtracting rest from task scans\n",
      "task_scans = data_no_0[..., neural_prediction_no_0 == 1]\n",
      "rest_scans = data_no_0[..., neural_prediction_no_0 == 0]\n",
      "difference = np.mean(task_scans, axis=-1) - np.mean(rest_scans, axis=-1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - showing slice 14 from the difference image\n",
      "plt.imshow(difference[:, :, 14], cmap='gray')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It looks like there's a voxel that is greater for activation than rest at about (i, j, k) == (45, 43, 14).\n",
      "\n",
      "Get and plot the values for this voxel position, for every volume in the 4D data (not including the first volume).  You can do it with a loop, but slicing is much nicer."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get the values for (i, j, k) = (45, 43, 14) and every volume\n",
      "# Plot the values (voxel time course)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Correlate the predicted neural time series with the voxel time course:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Correlation of predicted neural time course with voxel signal time course"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we will do a predicted hemodynamic time course using convolution.\n",
      "\n",
      "Next we need to get the HRF vector to convolve with.\n",
      "\n",
      "Remember we have defined the HRF as a function of time, not TRs.\n",
      "\n",
      "For our convolution, we need to *sample* the HRF at times corresponding the start of the TRs.\n",
      "\n",
      "So, we need to sample at times (0, 2.5, ...)\n",
      "\n",
      "Make a vector of times at which to sample the HRF.  We want to sample every TR up until (but not including) somewhere near 35 seconds (where the HRF should have got close to zero again)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make a vector of times at which to sample the HRF"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Sample your HRF function at these times and plot:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Sample HRF at given times\n",
      "# Plot HRF samples against times"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Convolve the predicted neural time course with the HRF samples:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Convolve predicted neural time course with HRF samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The default output of convolve is longer than the input neural prediction vector, by the length of the convolving vector (the HRF samples) minus 1.  Knock these last values off the result of the convolution to get a vector the same length as the neural prediction:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Remove extra tail of values put there by the convolution"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot the convolved neural prediction, and then, on the same plot, plot the unconvolved neural prediction."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot convolved neural prediction and unconvolved neural prediction"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Does the new convolved time course correlate better with the voxel time course?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Correlation of the convolved time course with voxel time course"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plot the hemodynamic prediction against the actual signal (voxel values).  Remember to use a marker such as '+' to give you a scatter plot.  How does it look?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}