curvefitting.ipynb

{
  "metadata": {
    "name": "Curve-Fitting"
  }, 
  "nbformat": 3, 
  "nbformat_minor": 0, 
  "worksheets": [
    {
      "cells": [
        {
          "cell_type": "heading", 
          "level": 1, 
          "metadata": {}, 
          "source": [
            "Curve-Fitting"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Here we are going to take a look at a few examples of fitting a function to\ndata. The first example shows how to fit an HRF model to noisy peristimulus\ntime-series data.\n\n", 
            "First, importing the necessary pieces:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "import numpy as np\nfrom scipy.stats import norm\n\nfrom mvpa2.support.pylab import pl\nfrom mvpa2.misc.plot import plot_err_line, plot_bars\nfrom mvpa2.misc.fx import *\nfrom mvpa2 import cfg"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "BOLD-Response parameters"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Let's generate some noisy \"trial time courses\" from a simple gamma function\n(40 samples, 6s time-to-peak, 7s FWHM and no additional scaling:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "a = np.asarray([single_gamma_hrf(np.arange(20), A=6, W=7, K=1)] * 40)\n# get closer to reality with noise\na += np.random.normal(size=a.shape)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Fitting a gamma function to this data is easy (using resonable seeds for the\nparameter search (5s time-to-peak, 5s FWHM, and no scaling):"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "fpar, succ = least_sq_fit(single_gamma_hrf, [5,5,1], a)"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "With these parameters we can compute high-resultion curves for the estimated\ntime course, and plot it together with the \"true\" time course, and the data:"
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "x = np.linspace(0,20)\ncurves = [(x, single_gamma_hrf(x, 6, 7, 1)),\n          (x, single_gamma_hrf(x, *fpar))]\n\n# plot data (with error bars) and both curves\nplot_err_line(a, curves=curves, linestyle='-')\n\n# add legend to plot\npl.legend(('original', 'fit'))\npl.title('True and estimated BOLD response')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n\\[Visit [http://pymvpa.org/examples/curvefitting.html](http://pymvpa.org/examples/curvefitting.html) to view this figure\\]"
          ]
        }, 
        {
          "cell_type": "heading", 
          "level": 2, 
          "metadata": {}, 
          "source": [
            "Searchlight accuracy distributions"
          ]
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "When doing a searchlight analysis one might have the idea that the\nresulting accuracies are actually sampled from two distributions: one\ncauses by an actual signal source and the chance distribution.  Let's\nassume the these two distributions can be approximated by a Gaussian,\nand take a look at a toy example, how we could explore the data.\n\n", 
            "First, we generate us a few searchlight accuracy maps that might have\nbeen computed in the folds of a cross-validation procedure. We\ngenerate the data from two normal distributions. The majority of\ndatapoints comes from the chance distribution that is centered at 0.5.\nA fraction of the data is samples from the \"signal\" distribution\nlocated around 0.75."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "nfolds = 10\nraw_data = np.vstack([np.concatenate((np.random.normal(0.5, 0.08, 10000),\n                                      np.random.normal(0.75, 0.05, 500)))\n                        for i in range(nfolds)])"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "Now we bin the data into one histogram per fold and fit a dual Gaussian\n(the sum of two Gaussians) to the total of 10 histograms."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "histfit = fit2histogram(raw_data,\n                        dual_gaussian, (1000, 0.5, 0.1, 1000, 0.8, 0.05),\n                        nbins=20)\nH, bin_left, bin_width, fit = histfit"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "All that is left to do is composing a figure -- showing the accuracy\nhistogram and its variation across folds, as well as the two estimated\nGaussians."
          ]
        }, 
        {
          "cell_type": "code", 
          "collapsed": false, 
          "input": [
            "# new figure\npl.figure()\n\n# Gaussian parameters\nparams = fit[0]\n\n# plot the histogram\nplot_bars(H.T, xloc=bin_left, width=bin_width, yerr='std')\n\n# show the Gaussians\nx = np.linspace(0, 1, 100)\n# first gaussian\npl.plot(x, params[0] * norm.pdf(x, params[1], params[2]), \"r-\", zorder=2)\npl.axvline(params[1], color='r', linestyle='--', alpha=0.6)\n# second gaussian\npl.plot(x, params[3] * norm.pdf(x, params[4], params[5]), \"b-\", zorder=3)\npl.axvline(params[4], color='b', linestyle='--', alpha=0.6)\n# dual gaussian\npl.plot(x, dual_gaussian(x, *params), \"k--\", alpha=0.5, zorder=1)\npl.xlim(0, 1)\npl.ylim(ymin=0)\n\npl.title('Dual Gaussian fit of searchlight accuracies')"
          ], 
          "language": "python", 
          "metadata": {}, 
          "outputs": []
        }, 
        {
          "cell_type": "markdown", 
          "metadata": {}, 
          "source": [
            "\n\n", 
            "And this is how it looks like.\n\n\\[Visit [http://pymvpa.org/examples/curvefitting.html](http://pymvpa.org/examples/curvefitting.html) to view this figure\\]"
          ]
        }
      ], 
      "metadata": {}
    }
  ]
}