big data

hw5-restart.ipynb

{
 "metadata": {
  "name": "hw5-restart"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "# HW5\n\n## Thoughts\nIf we are truly streaming data, we would have a fixed window on y (the results of our observation). For example, we might keep only the latest three values - the present value and two past ones. Or we might keep three values but think of them as one each of past, present, and future values. How we interpret the values depends on our field of view, known as $\\Delta$, the length of which is our window on y."
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "%pylab inline",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "Populating the interactive namespace from numpy and matplotlib\n"
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "import numpy as np\nimport scipy\nimport scipy.linalg\nfrom scipy.linalg import (toeplitz, svd, inv)\nfrom numpy import convolve as conv\nimport random",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "class functionInZ(object):\n    def __init__(self, vector, ixAtZero=0):\n        assert ixAtZero >= 0 and ixAtZero < len(vector)\n        self.vector = vector\n        self.ixAtZero = ixAtZero\n        \n    def star(self):\n        mirroredZeroIx = len(self.vector) - 1 - self.ixAtZero\n        return functionInZ(list(reversed(self.vector)), mirroredZeroIx)\n                           \n    def conv(self, other, mode='full'):\n        result = np.convolve(self.vector, other.vector, mode=mode)\n        newZeroIx = self.ixAtZero + other.ixAtZero  # ???\n        return functionInZ(result, newZeroIx)\n    \n    def highestIndex(self): return len(self) - 1 - self.ixAtZero\n    \n    def atIndex(self, key): \n        return self.vector[key + self.ixAtZero]  # right?\n    \n    def __getitem__(self, key): \n        return self.vector[key]\n    \n    def __len__(self): return len(self.vector)\n        \n    def __str__(self): return str(self.__dict__)\n        \n    def __repr__(self): print str(self)\n        \nv = functionInZ  # alias",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": "Given some a priori information (f, the covariance function of x) about the signal measured (x), the multiplicative noise of the measurement system (a), and the additive noise of the system (s, the covariance function of $\\nu$ in $y = A x + \\nu$), we can construct an estimate ($\\hat{x}$) for the original signal.\n\nThe equation for that estimate is\n\n$\\hat{x} = r \\ast y$\n\nwhere $r = P^{-1} q_\\Delta$, $P_{row,col} = p_{row-col}$, $p = a \\ast q + s$, $q = f \\ast a^*$ ($a^*$ is $a$ reflected across $i=0$), and $q_\\Delta$ is the values of $q$ within the field of view $\\Delta$.\n\n(What is eqn for variance of estimate?)"
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "# Instead of making up f and hoping it is of the correct form, generate it\nb = v([1, 1], 0)  # centered at i=0\nbStar = b.star()\nf = b.conv(bStar)\nprint f",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "{'ixAtZero': 1, 'vector': array([1, 2, 1])}\n"
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "a = v([1, 1], 0)  # asymm\nvariance_nu = 1  # s = \\sigma^2 = 1 for now\n\nlenMuStream = 100  # TODO replace with lazy iterator\nmu = np.random.randn(lenMuStream)\nmuIter = iter(mu)\ndef getNextMu(): return muIter.next()\n\nnu = np.random.rand(lenMuStream) * variance_nu  # could be less\nnuIter = iter(nu)\ndef getNextNu(): return nuIter.next()\n\nmus = []; xs = []; ys = []\n\n# prep arrays\nmus += [getNextMu() for i in range(len(b))]  # is len(b) right?",
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "def getNextX():\n    global mus\n    mus = mus[1:] + [getNextMu()]\n    print 'mus', mus\n    x = b.conv(v(mus, 0))\n    print 'x', x\n    return x[0]  # ok to only take this one element?\n\n[getNextX() for i in xrange(3)]",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "mus [2.1318788191415918, -0.14259514002645449]\nx {'ixAtZero': 0, 'vector': array([ 2.13187882,  1.98928368, -0.14259514])}\nmus [-0.14259514002645449, -1.6236415185022821]\nx {'ixAtZero': 0, 'vector': array([-0.14259514, -1.76623666, -1.62364152])}\nmus [-1.6236415185022821, 0.075522405441503765]\nx {'ixAtZero': 0, 'vector': array([-1.62364152, -1.54811911,  0.07552241])}\n"
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 63,
       "text": "[2.1318788191415918, -0.14259514002645449, -1.6236415185022821]"
      }
     ],
     "prompt_number": 63
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "\ndef getNextMeasurement(a, x, nu0):\n    y = a.conv(v(x, 0))  # should x be at i=0?\n    return y[0] + nu0\n\ndef getEstimate(ys):\n    return 0\n\n\ntry:\n    while(True):\n        xs += getNextX()\n        nu0 = getNextNu()\n        ys += getNextMeasurement(a, xs, nu0)\n        xHats += getEstimate(ys)\nexcept StopIteration as exc:\n    pass\n    \nplot(ys); title('y');\nplot(xHats); title(r'\\hat{x}')",
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "ename": "AssertionError",
       "evalue": "",
       "output_type": "pyerr",
       "traceback": [
        "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
        "\u001b[1;32m<ipython-input-69-808135105636>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     12\u001b[0m         \u001b[0mxs\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mgetNextX\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m         \u001b[0mnu0\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetNextNu\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m         \u001b[0mys\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mgetNextMeasurement\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mxs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnu0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     15\u001b[0m         \u001b[0mxHats\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mgetEstimate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mys\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m \u001b[1;32mexcept\u001b[0m \u001b[0mStopIteration\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mexc\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m<ipython-input-69-808135105636>\u001b[0m in \u001b[0;36mgetNextMeasurement\u001b[1;34m(a, x, nu0)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mgetNextMeasurement\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnu0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m     \u001b[0my\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0ma\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m  \u001b[1;31m# should x be at i=0?\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mnu0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;32m<ipython-input-24-3092ff4d3c79>\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, vector, ixAtZero)\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;32mclass\u001b[0m \u001b[0mfunctionInZ\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mobject\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0m__init__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvector\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mixAtZero\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m         \u001b[1;32massert\u001b[0m \u001b[0mixAtZero\u001b[0m \u001b[1;33m>=\u001b[0m \u001b[1;36m0\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mixAtZero\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvector\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvector\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mvector\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mixAtZero\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mixAtZero\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
        "\u001b[1;31mAssertionError\u001b[0m: "
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": "mus [2.3128279236018972, -0.19689004934565829]\nx {'ixAtZero': 0, 'vector': array([ 2.31282792,  2.11593787, -0.19689005])}\n"
      }
     ],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": "",
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}