playground.PCA.ipynb

{
 "metadata": {
  "name": "playground.PCA"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "from 'A Tutorial on PCA' by Lindsay I Smith"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "import numpy\n",
      "from numpy.fft import *\n",
      "from pprint import pprint\n",
      "import cmath"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 164
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "1. Background"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "1.1 Statistics"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "1.1.1. Standard Deviation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#sample\n",
      "X = [1, 2, 4, 6, 12, 15, 25, 45, 68, 67, 65, 98]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mean(data):\n",
      "    return float(sum(data)) / len(data)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 145
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sigma(sample):\n",
      "    mu = mean(sample)\n",
      "    return math.sqrt(sum([(x - mu)**2 for x in sample]) * 1./(len(sample) - 1.))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#tests\n",
      "set_1 = [0, 8, 12, 20]\n",
      "print 'set 1:\\n', set_1, '\\nstandard deviation: ', sigma(set_1) \n",
      "set_2 = [8, 9, 11, 12]\n",
      "print 'set 2:\\n', set_2, '\\nstandard deviation: ', sigma(set_2) "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "set 1:\n",
        "[0, 8, 12, 20] \n",
        "standard deviation:  8.32666399786\n",
        "set 2:\n",
        "[8, 9, 11, 12] \n",
        "standard deviation:  1.82574185835\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "1.1.2. Variance"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def variance(sample):\n",
      "    mu = mean(sample)\n",
      "    return (sum([(x - mu)**2 for x in sample]) * 1.) / (len(sample) - 1.)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 138
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#exercises:\n",
      "a = [12, 23, 34, 44, 59, 70, 98]\n",
      "b = [12, 15, 25, 27, 32, 88, 99]\n",
      "c = [15, 35, 78, 82, 90, 95, 97]\n",
      "for sample in [a, b, c]:\n",
      "    print sample, '\\nmean: ', mean(sample), '\\nstandard deviation: ', sigma(sample), '\\nvariance: ', variance(sample)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[12, 23, 34, 44, 59, 70, 98] \n",
        "mean:  48.5714285714 \n",
        "standard deviation:  29.5401260596 \n",
        "variance:  872.619047619\n",
        "[12, 15, 25, 27, 32, 88, 99] \n",
        "mean:  42.5714285714 \n",
        "standard deviation:  35.6036381234 \n",
        "variance:  1267.61904762\n",
        "[15, 35, 78, 82, 90, 95, 97] \n",
        "mean:  70.2857142857 \n",
        "standard deviation:  32.1751160874 \n",
        "variance:  1035.23809524\n"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "1.1.3. Covariance"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "One of the bottlenecks is calculation of the covariance matrix.  A straightforward approach is n^2d^2. where d is the number of dimensions.\n",
      "A method for efficient calculation uses fft to calculate the covariance between 2 matrices resulting in nlog(n)d^2."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "nlog(n) version:  \n",
      "fft code: http://jeremykun.com/2012/07/18/the-fast-fourier-transform/  \n",
      "fast covariance: http://research.google.com/pubs/pub36416.html  \n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy_pad as npp\n",
      "\n",
      "epsilon = 1.e-14\n",
      "\n",
      "def pad(n):\n",
      "    p = 1\n",
      "    while (p < n):\n",
      "        p <<= 1\n",
      "    p <<= 1\n",
      "    npp.pad\n",
      "    return  [0.] * (p - n)\n",
      "\n",
      "def getmax(A):\n",
      "    R = numpy.real(A)\n",
      "    maxi = numpy.amax(R)\n",
      "    mini = numpy.amin(R)\n",
      "    maxcorr = if abs(maxi) > abs(mini) else mini\n",
      "    index = np.argmax(R)\n",
      "    #print 'maxcorr: ', maxcorr\n",
      "    #pprint(A)\n",
      "    \n",
      "    return maxcorr, index\n",
      "    \n",
      "def fft_covariance(sample):\n",
      "    '''\n",
      "        returns the covariance matrix and index of each item for a dimension-list of sample-lists\n",
      "        [d0, d1, ...,dd] where d0 is the first dimension containing n samples.\n",
      "        n.b. numpy 1.7rc1 has numpy.pad to pad array.  currently running 1.6.2 (released version)\n",
      "        so imported rc1 version.\n",
      "    '''\n",
      "        \n",
      "    A = numpy.array(sample)\n",
      "    d, n = A.shape\n",
      "    Q = numpy.ones(1) if(d==2) else numpy.ones((d, d))\n",
      "    idx = [[1]] if(d==2) else [list([1.] * d) for i in range(d)]\n",
      "    nm1 = float(n - 1)\n",
      "    #prep for fft\n",
      "    A =  (A.transpose() - numpy.mean(A, axis=1)).transpose()\n",
      "    A = npp.pad(A,(0,len(pad(n))), 'constant', constant_values=(0,0))[:d, : ]\n",
      "    #2x2return\n",
      "    if (d == 2):\n",
      "        corr = ifft(fft(A[0, :]) * np.conjugate(fft(A[1, :])))\n",
      "        maxcorr = getmax(corr)\n",
      "        Q[0] =  maxcorr[0] / nm1\n",
      "        idx[0] = maxcorr[1]\n",
      "    else:\n",
      "        for i in range(d):\n",
      "            for j in range(i, d):\n",
      "                corr = ifft(fft(A[i, :]) * np.conjugate(fft(A[j, : ])))\n",
      "                maxcorr = getmax(corr)\n",
      "                Q[i][j] = maxcorr[0] / nm1 \n",
      "                idx[i][j] = maxcorr[1]\n",
      "                if i != j:\n",
      "                    Q[j][i] = Q[i][j]\n",
      "    return Q, idx\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 617
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "H = [9., 15., 25., 14., 10., 18., 0., 16., 5., 19., 16., 20.]\n",
      "M = [39., 56., 93., 61., 50., 75., 32., 85., 42., 70., 66., 80.]\n",
      "W = [120., 260., 490., 255., 150., 250., 0., 266., 88., 270., 270., 500. ]\n",
      "\n",
      "\n",
      "p = fft_covariance([H, M])\n",
      "print p\n",
      "\n",
      "one = [1,2,1]\n",
      "two = [-1, 1, 3]\n",
      "three = [1, 3, -1]\n",
      "c3 = covariance([one, two, three])\n",
      "print 'C3'\n",
      "pprint(c3)\n",
      "\n",
      "c5 = fft_covariance([one, two, three])\n",
      "print 'c5:'\n",
      "pprint(c5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(array([ 122.9469697]), [0])\n",
        "C3\n",
        "[[0.3333333333333333, 0.0, 1.0], [0.0, 4.0, -2.0], [1.0, -2.0, 4.0]]\n",
        "c5:\n",
        "(array([[ 0.33333333,  0.66666667,  1.        ],\n",
        "       [ 0.66666667,  4.        ,  2.        ],\n",
        "       [ 1.        ,  2.        ,  4.        ]]),\n",
        " [[0, 7, 0], [1.0, 0, 6], [1.0, 1.0, 0]])\n"
       ]
      }
     ],
     "prompt_number": 618
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def omega(p, q, inverse=False):\n",
      "    two = 2. if inverse else -2.\n",
      "    return cmath.exp((two * cmath.pi * 1j * q) / p)\n",
      "\n",
      "def my_fft(sample, inverse=False):\n",
      "    n = len(sample)\n",
      "    if n == 1:\n",
      "        return sample\n",
      "    else:\n",
      "        Feven = fft([sample[i] for i in xrange(0, n, 2)])\n",
      "        Fodd = fft([sample[i] for i in xrange(1, n, 2)])\n",
      "\n",
      "        combined = [0] * n\n",
      "        for m in xrange(n/2):\n",
      "            combined[m] = Feven[m] + omega(n, -m, inverse) * Fodd[m]\n",
      "            combined[m + n/2] = Feven[m] - omega(n, -m, inverse) * Fodd[m]\n",
      "\n",
      "        return combined"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 582
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "\n",
      "def pad(V):\n",
      "    n = len(V)\n",
      "    p = 1\n",
      "    while (p < n):\n",
      "        p <<= 1\n",
      "    p <<= 1\n",
      "    return  [0.] * (p - n)\n",
      "    \n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 473
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "d^2 version:\n",
      "To reduce the time the matrix is symmetrical complete for unordered pairs."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def sansmu(V):\n",
      "    vbar = mean(V)\n",
      "    return [v - vbar for v in V]\n",
      "\n",
      "\n",
      "def covariance(sample):\n",
      "    '''\n",
      "    calculate the covariance matrix for a given sample\n",
      "    '''\n",
      "    d = len(sample)\n",
      "    n = len(sample[0])\n",
      "    nm1 = float(n - 1)\n",
      "    Q =  [list([1.] * d) for i in range(d)]\n",
      "    mu0 = [sansmu(s) for s in sample]\n",
      "    for i in range(d):\n",
      "        for j in range(i, d):\n",
      "            Q[i][j] = sum([x * y  for x,y in zip(mu0[i], mu0[j])]) / nm1\n",
      "            if i != j:\n",
      "                Q[j][i] = Q[i][j]\n",
      "    return Q\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 292
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "        \n",
      "H = [9., 15., 25., 14., 10., 18., 0., 16., 5., 19., 16., 20.]\n",
      "M = [39., 56., 93., 61., 50., 75., 32., 85., 42., 70., 66., 80.]\n",
      "\n",
      "c1= covariance([H, M])\n",
      "pprint(c1)\n",
      "\n",
      "X = [10, 39, 19, 23, 28]\n",
      "Y = [43, 13, 32, 21, 20]\n",
      "c2 = covariance([X, Y])\n",
      "pprint(c2)\n",
      "\n",
      "one = [1,2,1]\n",
      "two = [-1, 1, 3]\n",
      "three = [1, 3, -1]\n",
      "c3 = covariance([one, two, three])\n",
      "pprint(c3)\n",
      "\n",
      "c5 = fft_covariance([one, two, three])\n",
      "print c5\n",
      "\n",
      "sone = sigma(one)\n",
      "stow = sigma(two)\n",
      "sthree = sigma(three)\n",
      "\n",
      "p = pad(one)\n",
      "c41 = [c - mean(one) for c in one]\n",
      "c42 = [c - mean(two) for c in one]\n",
      "c43 = [c - mean(three) for c in one]\n",
      "c41.extend(p)\n",
      "c42.extend(p)\n",
      "c43.extend(p)\n",
      "\n",
      "c4a = np.array([c41, c42, c43])\n",
      "print c4a.shape\n",
      "print c4a.size\n",
      "\n",
      "c4m = fftn(c4a, axes=[1])#, axes=1)\n",
      "corrm =ifftn(c4m *(np.conjugate(c4m)), axes=[1])\n",
      "\n",
      "realcorr = np.real(corrm)\n",
      "pprint(realcorr)\n",
      "\n",
      "pprint(realcorr / 2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "ValueError",
       "evalue": "setting an array element with a sequence.",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-585-1d7e1ab821a5>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     17\u001b[0m \u001b[0mpprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mc3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     18\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 19\u001b[1;33m \u001b[0mc5\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfft_covariance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtwo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthree\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     20\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mc5\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     21\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m<ipython-input-584-e8fd7093775e>\u001b[0m in \u001b[0;36mfft_covariance\u001b[1;34m(sample)\u001b[0m\n\u001b[0;32m     50\u001b[0m                 \u001b[0mmaxcorr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcorr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     51\u001b[0m                 \u001b[0mQ\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmaxcorr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m/\u001b[0m \u001b[0mnm1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 52\u001b[1;33m                 \u001b[0midx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmaxcorr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     53\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0mj\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     54\u001b[0m                     \u001b[0mQ\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mQ\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mValueError\u001b[0m: setting an array element with a sequence."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[47.71969696969697, 122.9469696969697],\n",
        " [122.9469696969697, 370.08333333333337]]\n",
        "[[115.7, -120.55], [-120.55, 138.7]]\n",
        "[[0.3333333333333333, 0.0, 1.0], [0.0, 4.0, -2.0], [1.0, -2.0, 4.0]]\n",
        "(0,)\n",
        "[]\n"
       ]
      }
     ],
     "prompt_number": 585
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Matrix Algebra"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}