Skip to content

Instantly share code, notes, and snippets.

@pckujawa
Created December 7, 2013 20:27
Show Gist options
  • Save pckujawa/7848150 to your computer and use it in GitHub Desktop.
Save pckujawa/7848150 to your computer and use it in GitHub Desktop.
big data
Display the source blob
Display the rendered blob
Raw
{
"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": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment