{
 "metadata": {
  "name": "",
  "signature": "sha256:284f3b63ff744d2d35699612393519fbd58f3db058038c92e29a8f723aa8c61a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Writing a custom function to do Pearson correlation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The numpy `np.corrcoef` function is frustrating, because it returns the correlation value of every input vector with every other, when often we just want the single correlation coefficient between two vectors.\n",
      "\n",
      "Returning all the coefficients, as numpy does, may also be inefficient.\n",
      "\n",
      "Here are some timings for a custom correlation routine, to check whether it is faster, when avoiding calculating all the coefficients."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import numpy.linalg as npl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "vec1 = np.random.normal(size=(100000,))\n",
      "vec2 = np.random.normal(size=(100000,))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use the IPython ``timeit`` magic command to time execution of `np.corrcoef` on these vectors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit np.corrcoef(vec1, vec2)[0, 1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 1.1 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The formula for Pearson's product-moment correlation coefficient from [wikipedia](https://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient) is:\n",
      "\n",
      "$$\n",
      "r = r_{xy} =\\frac{\\sum ^n _{i=1}(x_i - \\bar{x})(y_i - \\bar{y})}{\\sqrt{\\sum ^n _{i=1}(x_i - \\bar{x})^2} \\sqrt{\\sum ^n _{i=1}(y_i - \\bar{y})^2}}\n",
      "$$\n",
      "\n",
      "where $\\bar{x}$ is the mean:\n",
      "\n",
      "$$\n",
      "\\bar{x} = \\frac{1}{n} \\sum ^n _{i=1} x_i\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mycorr(x, y):\n",
      "    x_mean0 = x.ravel() - x.mean()\n",
      "    y_mean0 = y.ravel() - y.mean()\n",
      "    corr_top = x_mean0.dot(y_mean0)\n",
      "    # Use dot product on self to do sum of squares efficiently\n",
      "    corr_bottom = (np.sqrt(x_mean0.dot(x_mean0)) *\n",
      "                   np.sqrt(y_mean0.dot(y_mean0)))\n",
      "    return corr_top / corr_bottom"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mycorr(vec1, vec2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 563 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The custom routine is twice as fast as `np.corrcoef` for this simple case of just two vectors."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}