Started work on EKF. Mostly code to generate nonlinear transfer functions plots.

This commit is contained in:
Roger Labbe committed 2014-05-18 23:03:38 -07:00
1 parent 3b8654f70b
commit 0c4f93d3c5
4 files changed
+617

No files matched your search

+382
View File
@@ -0,0 +1,382 @@
{
"metadata": {
"name": "",
"signature": "sha256:a543c963e2375474d406b99de024c80765b4fd311e4b9ef7f8df686b92a8e0ac"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<center><h1>Kalman and Bayesian Filters in Python</h1></center>\n",
"<center><a href =\"http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/toc.ipynb\">Table of Contents</a></center>"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#format the book\n",
"%matplotlib inline\n",
"import book_format\n",
"book_format.load_style()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"html": [
"<style>\n",
" div.cell{\n",
" width: 850px;\n",
" margin-left: 0% !important;\n",
" margin-right: auto;\n",
" }\n",
" h1 {\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
"\t}\n",
"\t\n",
" div.input_area {\n",
" background: #F6F6F9;\n",
" border: 1px solid #586e75;\n",
" }\n",
"\n",
" .text_cell_render h1 {\n",
" font-weight: 200;\n",
" font-size: 36pt;\n",
" line-height: 100%;\n",
" color:#c76c0c;\n",
" margin-bottom: 0.5em;\n",
" margin-top: 1em;\n",
" display: block;\n",
" white-space: wrap;\n",
" } \n",
" h2 {\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" text-indent:1em;\n",
" }\n",
" .text_cell_render h2 {\n",
" font-weight: 200;\n",
" font-size: 20pt;\n",
" font-style: italic;\n",
" line-height: 100%;\n",
" color:#c76c0c;\n",
" margin-bottom: 1.5em;\n",
" margin-top: 0.5em;\n",
" display: block;\n",
" white-space: nowrap;\n",
" } \n",
" h3 {\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" }\n",
" .text_cell_render h3 {\n",
" font-weight: 300;\n",
" font-size: 18pt;\n",
" line-height: 100%;\n",
" color:#d77c0c;\n",
" margin-bottom: 0.5em;\n",
" margin-top: 2em;\n",
" display: block;\n",
" white-space: nowrap;\n",
" }\n",
" h4 {\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" }\n",
" .text_cell_render h4 {\n",
" font-weight: 300;\n",
" font-size: 16pt;\n",
" color:#d77c0c;\n",
" margin-bottom: 0.5em;\n",
" margin-top: 0.5em;\n",
" display: block;\n",
" white-space: nowrap;\n",
" }\n",
" h5 {\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" }\n",
" .text_cell_render h5 {\n",
" font-weight: 300;\n",
" font-style: normal;\n",
" color: #1d3b84;\n",
" font-size: 16pt;\n",
" margin-bottom: 0em;\n",
" margin-top: 1.5em;\n",
" display: block;\n",
" white-space: nowrap;\n",
" }\n",
" div.text_cell_render{\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" line-height: 135%;\n",
" font-size: 125%;\n",
" width:750px;\n",
" margin-left:auto;\n",
" margin-right:auto;\n",
" text-align:justify;\n",
" text-justify:inter-word;\n",
" }\n",
" div.output_subarea.output_text.output_pyout {\n",
" overflow-x: auto;\n",
" overflow-y: scroll;\n",
" max-height: 300px;\n",
" }\n",
" div.output_subarea.output_stream.output_stdout.output_text {\n",
" overflow-x: auto;\n",
" overflow-y: scroll;\n",
" max-height: 300px;\n",
" }\n",
" code{\n",
" font-size: 78%;\n",
" }\n",
" .rendered_html code{\n",
" background-color: transparent;\n",
" }\n",
" ul{\n",
" margin: 2em;\n",
" }\n",
" ul li{\n",
" padding-left: 0.5em; \n",
" margin-bottom: 0.5em; \n",
" margin-top: 0.5em; \n",
" }\n",
" ul li li{\n",
" padding-left: 0.2em; \n",
" margin-bottom: 0.2em; \n",
" margin-top: 0.2em; \n",
" }\n",
" ol{\n",
" margin: 2em;\n",
" }\n",
" ol li{\n",
" padding-left: 0.5em; \n",
" margin-bottom: 0.5em; \n",
" margin-top: 0.5em; \n",
" }\n",
" ul li{\n",
" padding-left: 0.5em; \n",
" margin-bottom: 0.5em; \n",
" margin-top: 0.2em; \n",
" }\n",
" a:link{\n",
" font-weight: bold;\n",
" color:#447adb;\n",
" }\n",
" a:visited{\n",
" font-weight: bold;\n",
" color: #1d3b84;\n",
" }\n",
" a:hover{\n",
" font-weight: bold;\n",
" color: #1d3b84;\n",
" }\n",
" a:focus{\n",
" font-weight: bold;\n",
" color:#447adb;\n",
" }\n",
" a:active{\n",
" font-weight: bold;\n",
" color:#447adb;\n",
" }\n",
" .rendered_html :link {\n",
" text-decoration: underline; \n",
" }\n",
" .rendered_html :hover {\n",
" text-decoration: none; \n",
" }\n",
" .rendered_html :visited {\n",
" text-decoration: none;\n",
" }\n",
" .rendered_html :focus {\n",
" text-decoration: none;\n",
" }\n",
" .rendered_html :active {\n",
" text-decoration: none;\n",
" }\n",
" .warning{\n",
" color: rgb( 240, 20, 20 )\n",
" } \n",
" hr {\n",
" color: #f3f3f3;\n",
" background-color: #f3f3f3;\n",
" height: 1px;\n",
" }\n",
" blockquote{\n",
" display:block;\n",
" background: #fcfcfc;\n",
" border-left: 5px solid #c76c0c;\n",
" font-family: 'Open sans',verdana,arial,sans-serif;\n",
" width:680px;\n",
" padding: 10px 10px 10px 10px;\n",
" text-align:justify;\n",
" text-justify:inter-word;\n",
" }\n",
" blockquote p {\n",
" margin-bottom: 0;\n",
" line-height: 125%;\n",
" font-size: 100%;\n",
" }\n",
"</style>\n",
"<script>\n",
" MathJax.Hub.Config({\n",
" TeX: {\n",
" extensions: [\"AMSmath.js\"]\n",
" },\n",
" tex2jax: {\n",
" inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
" displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
" },\n",
" displayAlign: 'center', // Change this to 'center' to center equations.\n",
" \"HTML-CSS\": {\n",
" styles: {'.MathJax_Display': {\"margin\": 4}}\n",
" }\n",
" });\n",
"</script>\n"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 1,
"text": [
"<IPython.core.display.HTML at 0x1bdd710>"
]
}
],
"prompt_number": 1
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Kalman filter that we have developed to this point is extremely good, but it is also limited. Its derivation is in the linear space, and hence it only works for linear problems. Let's be a bit more rigorous here. You can, and we have in this book, apply the Kalman filter to nonlinear problems. For example, in the g-h filter chapter we explored using a g-h filter in a problem with constant acceleration. It 'worked', in that it remained numerically stable and the filtered output did track the input, but there was always a lag. It is easy to prove that there will always be a lag when $\\mathbf{\\ddot{x}}>0$. The filter no longer produces an optimal result. If we make our time step arbitrarily small we can still handle many problems, but typically we are using Kalman filters with physical sensors and solving real-time problems. Either fast enough sensors do not exist, are prohibitively expensive, or the computation time required is excessive. It is not a workable solution.\n",
"\n",
"The early adopters of Kalman filters were the radar people, and this fact was not lost on them. Radar is inherently nonlinear. Radars measure the slant range to an object, and we are typically interested in the aircraft's position over the ground. We invoke Pythagoras and get the nonlinear equation:\n",
"$$x=\\sqrt{slant^2 - altitude^2}$$\n",
"\n",
"So shortly after the Kalman filter was enthusiastically taken up by the radar industry people began working on how to extend the Kalman filter into nonlinear problems. It is still an area of ongoing research, and in the Unscented Kalman filter chapter we will implement a powerful, recent result of that research. But in this chapter we will cover the most common form, the Extended Kalman filter, or EKF. Today, most real world \"Kalman filters\" are actually EKFs. The Kalman filter in your car's and phone's GPS is an EKF, for example. \n",
"\n",
"###The Problem with Nonlinearity\n",
"\n",
"You may not realize it, but the only math you really know how to do is linear math. Equations of the form \n",
"$$ A\\mathbf{x}=\\mathbf{b}$$.\n",
"\n",
"That may strike you as hyperbole. After all, in this book we have integrated a polynomial to get distance from velocity and time:\n",
" We know how to integrate a polynomial, for example, and so we are able to find the closed form equation for distance given velocity and time:\n",
"$$\\int{(vt+v_0)}\\,dt = \\frac{a}{2}t^2+v_0t+d_0$$\n",
"\n",
"That's nonlinear. But it is also a very special form. You spent a lot of time, probably at least a year, learning how to integrate various terms, and you still can not integrate some arbitrary equation - no one can. We don't know how. If you took freshman Physics you perhaps remember homework involving sliding frictionless blocks on a plane and other toy problems. At the end of the course you were almost entirely unequipped to solve real world problems because the real world is nonlinear, and you were taught linear, closed forms of equations. It made the math tractable, but mostly useless. \n",
"\n",
"The mathematics of the Kalman filter is beautiful in part due to the Gaussian equation being so special. It is nonlinear, but when we add and multipy it using linear algebra we get another Gaussian equation as a result. That is very rare. $\\sin{x}*\\sin{y}$ does not yield a $\\sin(\\cdot)$ as an output.\n",
"\n",
"### The Effect of Nonlinear Transfer Functions on Gaussians\n",
"\n",
"Unfortunately Gaussians are not closed under an arbitrary nonlinear function. Recall the equations of the Kalman filter - at each step of its evolution we do things like pass the covariances through our process function to get the new covariance at time $k$. Our process function was always linear, so the output was always another Gaussian. Let's look at that on a graph. I will take an arbitrary Gaussian and pass it through the function $f(x) = 2x + 1$ and plot the result. We know how to do this analytically, but lets do this with sampling. I will generate 500,000 points on the Gaussian curve, pass it through the function, and then plot the results. I will do it this way because the next example will be nonlinear, and we will have no way to compute this analytically."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from numpy.random import normal\n",
"\n",
"normals = normal(loc=0.0, scale=1, size=500000)\n",
"ys = 2*normals + 1\n",
"\n",
"plt.hist(ys,1000)\n",
"plt.show()"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAssAAAFyCAYAAAADJZf7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3VGMVOd99/HfmbM7w8Iua2PYl111afymAhpCCOBWdcCV\nwG1vbF9EEUTxRRurqKjOjS8sldveNKiqJVrJbomSm0qtZKo4UoMvotq9CW61aoqFUCTULQrmBc8y\ngcUzOzucOTNn5r0YzmF2mLNnZ3d25plnvh8pss9hWA7Os//57zO/53mcer1eFwAAAICnpPr9AAAA\nAICpaJYBAACAGDTLAAAAQAyaZQAAACAGzTIAAAAQg2YZAAAAiLFqs3z37l0dP35cX/3qV3X06FF9\n9NFHkiTXdXX48GEdPnxYb731VvT6S5cuae/evdq3b58uX76ceB8AAAAwmbPaPsu5XE737t3TwYMH\ndfv2bX3jG9/QnTt3NDExoaWlpRWv9X1f+/fv19zcnDzP04kTJ/S///u/sfcBAAAA042s9otTU1Oa\nmpqSJO3Zs0e+78v3/bavnZub04EDB7Rr1y5J0uzsrK5du6ZCodD2/qFDh7r59wAAAAC6btVmudnP\nfvYzHT16VOl0Wp7n6ejRoxobG9P3v/99vfTSS7p3756mp6d18eJF7dixQ7t371Y2m1WxWGx7n2YZ\nAAAApltTs7ywsKC3335b//qv/yqpkWWemprSL37xC33zm9/U/Py8wjTH2bNnJUkffPDBiq/RfN9x\nnK79BQAAAIDNktgse56nU6dO6Z133tHzzz8vSVE044UXXtDMzIw+++wzzczMKJvNRr9vYWFBMzMz\nWlpaeur+9PT0U3/O+++/r507d274LwQAAADE8TxPr7zyyppfv2qzXK/X9cYbb+j111/XH/3RH0mS\nHj58qC1btmhsbEy3bt3S3bt3tWfPHj3//PP65S9/qV//+tfyPE937tzR1772Nfm+3/Z+q507d+rI\nkSMd/nWBzffTn/5Ur732Wr8fA2iL8QlTMTZhqqtXr3b0+lWb5U8++UQ//vGPdePGDf3gBz+Q4zh6\n99139cYbbyiTych1Xf3oRz/S2NiYJOn8+fM6duyYJOnChQuSpHQ63fY+MChu377d70cAYjE+YSrG\nJmyxarN8/Pjxtrtf3Lhxo+3rT58+rdOnT6/5PgAAAGAyTvADEvzBH/xBvx8BiMX4hKkYm7AFzTKQ\nYN++ff1+BCAW4xOmYmzCFjTLQIIrV670+xGAWIxPmIqxCVvQLAMAAAAxaJaBBMePH+/3IwCxGJ8w\nFWMTtqBZBgAAAGLQLAMJyN3BZIxPmIqxCVvQLAMAAAAxaJaBBOTuYDLGJ0zF2IQtaJYBAACAGDTL\nQAJydzAZ4xOmYmzCFjTLAAAAQAyaZSABuTuYjPEJUzE2YQuaZQAAACAGzTKQgNwdTMb4hKkYm7AF\nzTIAAAAQg2YZSEDuDiZjfMJUjE3YgmYZAAAAiEGzDCQgdweTMT5hKsYmbEGzDAAAAMSgWQYSkLuD\nyRifMBVjE7agWQYAAABi0CwDCcjdwWSMT5iKsQlb0CwDAAAAMWiWgQTk7mAyxidMxdiELWiWAQAA\ngBg0y0ACcncwGeMTpmJswhY0ywAAAEAMmmUgAbk7mIzxCVMxNmELmmUAAAAgBs0ykIDcHUzG+ISp\nGJuwBc0yAAAAEINmGUhA7g4mY3zCVIxN2IJmGQAAAIhBswwkIHcHkzE+YSrGJmxBswwAAADEoFkG\nEpC7g8kYnzAVYxO2oFkGAAAAYtAsAwnI3cFkjE+YirEJW9AsAwAAADFoloEE5O5gMsYnTMXYhC1o\nlgEAAIAYNMtAAnJ3MBnjE6ZibMIWNMsAAABADJplIAG5O5iM8QlTMTZhC5plAAAAIAbNMpCA3B1M\nxviEqRibsAXNMgAAABCDZhlIQO4OJmN8wlSMTdiCZhkAAACIQbMMJCB3B5MxPmEqxiZsQbMMAAAA\nxKBZBhKQu4PJGJ8wFWMTtqBZBgAAAGLQLAMJyN3BZIxPmIqxCVvQLAMAAAAxaJaBBOTuYDLGJ0zF\n2IQtaJYBAACAGKs2y3fv3tXx48f11a9+VUePHtVHH30kSbp06ZL27t2rffv26fLly9HrO70PDAJy\ndzAZ4xOmYmzCFiOr/eLo6Kj+/u//XgcPHtTt27f1jW98Q7/61a907tw5zc3NyfM8nThxQq+++qp8\n3+/oPgDALtlCWZI0vT3T5ycBgO5ZtVmemprS1NSUJGnPnj3yfV//+Z//qQMHDmjXrl2SpNnZWV27\ndk2FQqGj+4cOHdrMvxfQNeTuYDKTxmeu6EuiWUaDSWMT2IhVm+VmP/vZz3T06FHlcjlNT0/r4sWL\n2rFjh3bv3q1sNqtisdjRfZplAAAAmG5NC/wWFhb09ttv67333ovunT17VqdOnXrqtWu57zjOep8X\n6DlydzAZ4xOmYmzCFokzy57n6dSpU3rnnXf0/PPP6/PPP1c2m41+fWFhQTMzM1paWlrz/enp6bZ/\n1ptvvqk9e/ZIkiYnJ3Xw4MHoY5zwm45rrrnmmmszryf+7yGjnofr/l6HTHkerof3+vr168rn85Kk\n27dv68yZM+qEU6/X63G/WK/X9frrr+v3f//39ed//ueSJN/3tX///mjB3smTJzU/P9/x/VYff/yx\njhw50tHDAwDMce3zJUnSoZmJFfdZ+AfAJFevXtXLL7+85tePrPaLn3zyiX784x/rxo0b+sEPfiDH\ncfThhx/q/PnzOnbsmCTpwoULkqR0Ot3RfQDAcGDhH4BBturMci8xswxTXblyJfo4BzCNSeMzbmY5\n7j7sZtLYBJp1OrPMCX4AAABADJplIAEzIzAZ4xOmYmzCFjTLAIBItlCOFuQBAGiWgUSt2yABJun2\n+MwV/WhBHrAR1E7YgmYZABCLmWYAw27VreMAkLuD2TZ7fOaKvral3cTXZQtl+UFNaTd+DiZbKKvk\nB9qadtlGbghQO2ELZpYBAKsazyQ3y7miLz9YfSfSXNHXAjEPAAOGZhlIQO4OJuv1+Ow0lpEtlBWY\nsZ0/eozaCVvQLAMAJD2JUripJ02xH9RWvCZpAWD4e5tfH9RolgEMLjLLQAJydzBZN8dnGKUIajUF\ntUZD3BytyBbKqtTqGk05sV8j7wUKan5iJjlbKGs842oiMxJdk2e2C7UTtqBZBoAhFc4AtzancWv0\nckVf9XpdkrPi97aLWoTNb+vMdPPXcpx01Cznir6KfqBxmmUAhiGGASQgdweTbWR8Nkcq2kUu4rgp\nKbtUjn5va9TCTUl38p4WYhb9Ne+uwdZ09qJ2whbMLAMAntrNwk09iWAs+8GKhjbvBUq7jly3fRwj\n7wWr/lmOo6gxb5d/bjfjHTcLDgCbjZllIAG5O5ism+PTdZ40v43scaNZXixVurrd22KpIj+oa9kP\n2u6U0W4RIScLDh5qJ2xBswwAQ6Rd7KF5F4z1fL3W+EZz072axVJl1Z0ysoWyHizTIAPoL5plIAG5\nO5is0/EZN2ubdKDIal9vRXzDcdbVdMd97Qrbzg0saidsQWYZAIbARhfRbUu7WvZXzyJLK3fScB2n\nqweSrOXYbQDoNmaWgQTk7mCy1vEZt7vERjO/aznyutWKxjmljhrnzEhKflBbsbhwPc+A/qF2whY0\nywBgkXZN8WpHTnfaxK739zYvGFzL11kqV+UH9a4vLgSATtEsAwnI3cFkaxmfqx05vVoT22y5zQEj\na/29Sbr1dWAWaidsQbMMAEgUbvcGAMOGBX5AAnJ3MFm78emmtGrO1001Yg4mWutCQpiP2glb0CwD\ngGUasYZGzrdUcdvGJ0xtSMczNMsAzEIMA0hA7g4mSxqftsUn4nb7gHmonbAFzTIADKG4o6ZNx7HX\nAHqNZhlIQO4OJlvv+Ew6atokbkrRnsutkRKYi9oJW9AsA4CFbDrtrpHBfhIpsenvBsB8NMtAAnJ3\nMFnc+LTttLvm0wBt+7vZitoJW9AsA8AA2zY1a9WCt7VELZqPwAaAzUazDCQgdweTpcZ3RAveVjvW\nWtrY0da9spbdOzgCezBQO2ELmmUAGHDhISSrHWst2XWstMkHqwCwC80ykIDcHUyWz+eV94Khm2kN\n/87EMcxF7YQtaJYBAAOpNY7BgSUANgPNMpCA3B1MNjk5Kddx+v0YfRNGUCQOLDENtRO2GOn3AwAA\nNsYd4mmPvBdIokEGsHmGuMQCa0PuDqbKFsoqPvL6/Rh915zZbp5pRn9RO2ELmmUAGFC5oq9qfXgj\nGO0M42JHAJuLZhlIQO4Opmi3gG1khDQdzETthC1olgFgQHSygG0tJ+HZYpgXOALYfDTLQAJyd+i3\ndjPK4Wl91Wr7gznWchKeLYZ5gaPJqJ2wBZ/fAYDh2s0mJ53WBwDoDn4eBxKQu4PJyCw3uCkpqNej\nf28+CpvDSvqD2glb0CwDAAZe3guimfa8F2jZD6Jf47ASABtBswwkIHcHk8VlloF+o3bCFjTLAAAA\nQAyaZSABuTuYwE2p7VZwZJZhKmonbEGzDACGal6YlveCaCu4bKH8VOPcvMBt2LmOo2U/YFEfgK6g\nWQYSkLtDv8QtTMsV/ahxDjPLzQvchp2bauwzzaK+/qJ2whY0ywAwIFzHYQa5A2whB6AbaJaBBOTu\n0E+t+we3ziCTWY7XvIXctrTLFnI9Ru2ELaiyAGCY5tnPvBco7Tp9fBo7jGfcqHEO//tOb8/085EA\nDAhmloEE5O7QK9lCWUvlakczoOyz3DlmmHuD2glbMLMMAIbIFX05Trrfj2GVZT9ou+UeAKwVM8tA\nAnJ3MBmZ5dUtlirRziHoLWonbEGzDACwHjPMANaLZhlIQO4OvdRpU0dmeW3CGeZtabffjzLine truncated
"text": [
"<matplotlib.figure.Figure at 0x2c48990>"
]
}
],
"prompt_number": 2
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is an unsuprising result. The result of passing the Gaussian through $f(x)=2x+1$ is another Gaussian centered around 1. Let's look at the input, transfer function, and output at once."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from nonlinear_plots import plot_transfer_func\n",
"\n",
"def g(x):\n",
" return 2*x+1\n",
"\n",
"plot_transfer_func (normals, g, lims=(-10,10), num_bins=100)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAsUAAAF9CAYAAADsuhWdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X18FfWZ///3OSEnt+SOEJJAwl0IchMQKCBKq0B1raJb\n3P5ge7NqK9UK335Lv9sb29r+pK3Kaumi/aF1qy2L1l3cXW1XXG/BqvEG1CBGESGBJBCSkJBbQk5O\nbub3h19SuT8TZs7MOfN6Ph48Hsw5c2auXHzmw3Um18z4DMMwBAAAAHiY3+kAAAAAAKdRFAMAAMDz\nKIoBAADgeRTFAAAA8DyKYgAAAHgeRTEAAAA8j6IYAAAAnkdRjKg1ZswYrV692vb9+P1+bdy40fb9\nAAAA51AUI2r5fD7b93H82TY84wYAgNhGUYyIu/POO1VQUKDExERNmzZNzzzzzAnv+/1+vfrqqwPL\nGzZskN//16E6ZswY+f1+VVdXa/Xq1fL7/fL7/fr5z38+sM5ll12mZcuW6e///u+VnJyscePG6Ykn\nnhh4v6qqSn6/XzU1NQOv3XHHHRo7duwJccTFxUmSvv71rw/sh7PGAADEHopiRNQf/vAH3XnnnVqz\nZo0++OADXXXVVbruuuu0b9++sLfx7rvvqq6uTqNGjdL3vvc91dfXq76+Xv/4j/94wnpPPvmkJk+e\nrJ07d+ob3/iGvvrVr2rv3r1n3fanzz7X19errq5OknTfffcN7Gfp0qUmfmIAABANKIoRUQ888ID+\n4R/+QV/96ldVVFSkNWvWqLCwUA899FDY2xg2bJhGjBihuLg4paamKicnRzk5OUpJSTlhvQkTJuhn\nP/uZJkyYoNtvv11FRUX63e9+d9Ztf7pNIicnRyNGjJAkpaenD+wnMTHRxE8MAACiAUUxIqqiokLT\np08/4bXp06ef8wzuYEydOvWU5YqKCsv3AwAAoh9FMRxnGMZZL5rr7++3bD/HnW5/Vu0HAABEH4pi\nRFRRUZHee++9gWXDMLRz504VFRUNvJaRkaH29vaB5erq6tNuKxAIqKen54z7Ki8vP2H5gw8+GNhP\nRkaGJKmjo+OE/ZyuWI6Pj1dvb+/ZfiwAABDlKIoRUbfeeqsee+wxPfbYY9qzZ49uu+02HThwQLfc\ncsvAOrNnz9bjjz8uwzBUXV2tf/3Xfz3ttoqLi/Xiiy+qtrZWwWBQfX19J7xfUVGhX/ziF9qzZ49+\n+ctfqrKyUsuXL5f0SY/whAkT9Oijj0qSdu7cqT//+c9n3M9///d/q6mpScFgkDPKAADEIIpiRNQ3\nvvEN/ehHP9KPfvQjlZSU6Nlnn9V//ud/aty4cQPr/PrXv9ZHH32k4cOH68Ybb9T1119/2jO4d911\nl/x+vyZOnKjk5GTdeeedA+/5fD4tWbJEO3fu1IUXXqjf//73evTRR1VcXDywzr/8y7/oqaeeUk5O\njn7605/qa1/72mlj/s1vfqOqqioVFBQoOTlZjz32mIUZAQAAbuAzBvlUgu9973t67LHHNHz48IFf\nUz/xxBO6/fbb5fP5tHbtWi1evNjSYIFwLViwQGPHjtXvf/97p0MBXIE5GwDObtBF8ZtvvqlAIKAb\nb7xR5eXlCoVCuuCCC7Rt2zYFg0EtWLCAK/3hmMsuu0xjx47VH/7wB6dDAVyBORsAzm7Q7RPz5s3T\nsGHDBpa3bdumKVOmaPjw4SooKFBBQYF27txpSZCAWT6fLyKPgQaiBXM2AJzdEKs2VF9fr7y8PD30\n0EPKyspSbm6u6urqTrknLRAJL7/8stMhAK7GnA0AJ7KsKD7u+F0EnnzySc7UAYDLMWcDwCcsK4rz\n8/NVV1c3sHz8LMTJ/vgfTyk3K82q3QJAxASDQV199dVOh2GJcOfsTZs2KTs7O5KhAYBlzMzblhXF\ns2fP1ocffqjGxkYFg0EdPHhQ06ZNO2W93Kw0zZw506rdxrz169dr5cqVTocRNciXOeTLnLKyMqdD\nsEy4c3Z2djZztgkcU+aQr08eYvVoWb22VDTrzivHa1R64lnXJ2fmmJm3B10Ur1y5Uk899ZSamppU\nUFCgBx54QGvWrNEll1wiSVq3bt1gNw0AsBhzNuA+vf2G1r1Wo+rWoNZdW6zMpHinQ/K0QRfF69ev\n1/r16095fenSpecVEE5UWFjodAhRhXyZQ768gzk7MjimzPFyvjpDffrFlv0KxPl0z1VFSoqPC+tz\nXs6Z3Sy/0A7W+vQT2HBu5Msc8gVYi2PKHK/mq6kzpNuf36fJI1K0ct4oxfnDv8jVqzmLBB7z7HKN\njY1OhxBVyJc55AuwFseUOV7MV1VLl1Y9vUeXjc/Qty82VxBL3sxZpHCmGAAAIALeO9Shu7ZW6VsX\njdTCoiynw8FJKIpdbv78+U6HEFXIlznkC7AWx5Q5XsrX1opm/fatWv144RhdmD900NvxUs4ijaIY\nAADAJoZhaNP7Ddr8UZP+6aoijc1KcjoknAE9xS5XWlrqdAhRhXyZQ74Aa3FMmRPr+errN/SbNw7q\nL5UtWndNsSUFcaznzEmcKQYAALBYV0+f7n65St29htYuLlZKILxbrsE5nCl2OXqHzCFf5pAvwFoc\nU+bEar5aunr0g/+pUGrCEP3yb8ZZWhDHas7cgDPFAAAAFqltC+onz1dqwfgsXT8zVz6fuVuuwTmc\nKXY5eofMIV/mkC/AWhxT5sRavnY1dOr/bN6rZdNG6IZZebYUxLGWMzdx5ExxY2dIw1MCTuwaAADA\ncqVVrbqv9IC+f2mh5hSkOx0OBsGRoriunaI4XPQOmUO+zCFfgLU4psyJlXz96cNGbdrZoDuvHK/i\n7GRb9xUrOXMjR9on9jYdk2EYTuwaAADAEv2GoX/ZVqundzXqn6+ZYHtBDHs5UhQ3HA3p8NEeJ3Yd\ndegdMod8mUO+AGtxTJkTzfkK9fbr7q1V2t3YqX++pli5QxMist9ozpnbOVIUf35Cll7Z3+LErgEA\nAM5Le7BXtz1XIUlac2WR0hK5mVcscKQoLs5OVuWRLj28vVbtwV4nQoga9A6ZQ77MIV+AtTimzInG\nfNV3dOu7T+/RBcNT9KOFYxQYEtlSKhpzFi0cuyXbjxaMUdGwZJVWtdJfDAAAXG9v0zF99+m9Wjwp\nWzfPHSk/9yCOKY7ep3je6HRVHOlSTz9F8ZnQO2QO+TKHfAHW4pgyJ5rytf1Am378XKVWzhulJVNz\nHIsjmnIWbRwtihOG+PU3xVn69as12t/c5WQoAAAAp/Xsx0e09tUa3XH5WM0fm+F0OLCJ453hE4en\n6LqpPj2354gS4vy6MD9VM0emOR2Wa9A7ZA75Mod8AdbimDLH7fkyDEOPltVrS0Wz1i6eoFHpiU6H\n5PqcRTPHi2JJKh6erOLhyWrt6tH6Nw6quiWoxCF+feGCbKdDAwAAHtTbb2jdazWqbg1q3bXFykyK\ndzok2MzR9omTZSTF6zvzC7SoKEsH27pV196tho6Q02E5it4hc8iXOeQLsBbHlDluzVdnqE+3P1+p\n9u5e3XNVkasKYrfmLBa4qiiWpNSEIUpLHKIrirP0QcNRbXj3ELdtAwAAEdHUGdI/bt6r/LQE/b+f\nH6ek+DinQ0KEuKJ94nRGZyZpdGaSclIC+q8PDivO51O/YWhkeoIunzDM6fAiht4hc8iXOeQLsBbH\nlDluy1dVS5duf75Siydla9m0EfK58JZrbstZLHFtUXzc9Pyhmp4/VJLU09evX26pUl37X1sqEob4\ntWz6CKfCAwAAMeC9Qx26a2uVvnXRSC0synI6HDjAlvaJuLg4zZgxQzNmzNCqVass2258nF+rrxin\n62flDfw51tOn9w51nPKn8sgxy/brJHqHzCFf5pAvSPbN2V7EMWWOW/K1taJZd22t0o8XjnF9QeyW\nnMUiW84UJycna8eOHXZs+hTXTMrWofbuU15/5qMj+t/zkyMSAwBEs0jO2YCbGIahTe836OldTfqn\nq4o0NivJ6ZDgINe3T5xLdkpA2SmBU15v7+7TxnfrBpaPHOvR5ROyNDU3NZLhnTd6h8whX+aQL8Ba\nHFPmOJmvvn5D6988qF0NR3XftcWnrSXciDFmH1uK4mAwqFmzZikpKUl33323PvvZz9qxm7OaPyZD\n88f89akzNa1BPfXBYZXVdpzxM8Hefv3t5OEaMTQ6DgwAsIIb5mwgkrp6+nT3y1Xq7jW0dnGxUgLc\nYQI2FcW1tbXKycnRO++8oyVLlqiiokIJCQkD769YsUKFhYWSpPT0dJWUlAx88zneK2PH8nfmF571\n/XcPtuuhl3ZoWMDQhAkTJEl79+6VJC2cXaJJOSm2xne65QcffDBi+YmFZfJFvqzOT3l5+cB8dcUV\nVygWuXXOjsZljin356uzV9rcNlwFGYlalHRAO7bXuyYf4SyXl5fr1ltvdU08blsuLy9XW1ubJKmm\npkbLly9XuHyGYRhhrz0Ic+fO1caNGzVx4kRJ0pYtWzRz5kw7dzloff2GDnee/mEhT+9q0t+f4y4X\ngSF+JQ6x9trF0tLSgX9snBv5Mod8mVNWVqZFixY5HYatomnOdiOOKXMina/atqB+8nylFozP0vUz\nc115y7VzYYyZY2beHmL1zltaWpSYmKikpCRVVVWptrZ24AyD28X5fcobmnDa9y7ISdbWypazfv5Q\ne7dWzBtlaUwMfHPIlznkC9E8Z7sRx5Q5kczXroZOrX5pn26clacvXJAdsf1ajTFmH8uL4t27d+vr\nX/+6EhISFBcXp0ceeURJSdF/Nefnxmaec50/7qg/4eK+MzlyrEc3zx1JDxMAx8XqnA18WmlVq+4r\nPaDvX1qoOQXpTocDl7K8KJ43b552795t9Wajwldn5Ia13gt7jujxHfVKCqMorqmuVuHo0RLine truncated
"text": [
"<matplotlib.figure.Figure at 0x2ccc410>"
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The plot labelled 'input' is the histogram of the original data. This is passed through the transfer function $f(x)=2x+1$ which is displayed in the chart to the upper right. The red lines shows how one value, $x=0$ is passed through the function. Each value from input is passed through in the same way to the output function on the left. The output looks like a Gaussian, and is in fact a Gaussian. We can see that it is altered -the variance in the output is larger than the variance in the input, and the mean has been shifted from 0 to 1, which is what we would expect given the transfer function $f(x)=2x+1$ The $2x$ affects the variance, and the $+1$ shifts the mean.\n",
"\n",
"Now let's look at a nonlinear function and see how it affects the probability distribution."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from nonlinear_plots import plot_transfer_func\n",
"\n",
"def g(x):\n",
" return (np.cos(4*(x/2+0.7)))*np.sin(0.3*x)-0.9*x\n",
"\n",
"plot_transfer_func (normals, g, lims=(-4,4), num_bins=300)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAr0AAAF9CAYAAAAJJNDxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VOXdPvD7zJpMlsm+kslCEtYQMOyOKKi0KvQVa3Fr\n+aHNiwVti5W22ooVq77Y1ortq1hbrUVfW1cqSqnFIGhUEIOEIFsCJCEb2fdMZv39gUQSssxMZuY8\nM7k/18V1cTInc+48OXnyzTPfOUdyOBwOEBEREREFMIXcAYiIiIiIvI1FLxEREREFPBa9RERERBTw\nWPQSERERUcBj0UtEREREAY9FLxEREREFPBa9RERERBTwWPSS8NLS0rBhwwavH0ehUGDLli1ePw4R\nERH5HoteEp4kSV4/xvl7tPBeLURERIGJRS95zaOPPoqUlBQEBQVh2rRp2L59e7/HFQoFPvzww77t\nF198EQrF16dkWloaFAoFKioqsGHDBigUCigUCjz88MN9+1xxxRW46aabcPPNN0On0yEjIwOvvfZa\n3+Pl5eVQKBSorKzs+9hDDz2E9PT0fjmUSiUA4Pbbb+87Dld9iYiIAgeLXvKKv/71r3j00UexceNG\nHD58GNdeey1uuOEGnDp1yunnKCoqQm1tLcaNG4d169ahrq4OdXV1uPfee/vt99Zbb2Hy5MkoLi7G\nHXfcgdtuuw2lpaXDPveFq8d1dXWora0FADz11FN9x1m+fLkLXzERERGJjEUvecUzzzyD733ve7jt\nttuQmZmJjRs3wmAw4E9/+pPTzxEdHY34+HgolUqEhoYiLi4OcXFxCAkJ6bdfVlYWHnzwQWRlZeGB\nBx5AZmYm/vznPw/73Be2McTFxSE+Ph4AoNfr+44TFBTkwldMREREImPRS15RVlaG3Nzcfh/Lzc0d\ncQXWHVOnTr1ou6yszOPHISIiIv/Fopd8xuFwDPumNLvd7rHjnDfY8Tx1HCIiIvIfLHrJKzIzM3Hw\n4MG+bYfDgeLiYmRmZvZ9LCIiAu3t7X3bFRUVgz6XRqOBxWIZ8lglJSX9tg8fPtx3nIiICABAR0dH\nv+MMVgyr1WpYrdbhviwiIiLyUyx6yStWr16Nl19+GS+//DJOnDiB++67D2fOnMGdd97Zt8+sWbPw\nyiuvwOFwoKKiAn/7298Gfa7s7Gzs3LkT1dXVMJlMsNls/R4vKyvDr3/9a5w4cQKPPPIITp48ifz8\nfADnenSzsrLw0ksvAQCKi4vx9ttvD3mcbdu2obGxESaTiSvCREREAYRFL3nFHXfcgfvvvx/3338/\ncnJysGPHDrzxxhvIyMjo2+f3v/89jh49itjYWKxcuRIrVqwYdAX2scceg0KhwIQJE6DT6fDoo4/2\nPSZJEpYtW4bi4mJMnz4dL7zwAl566SVkZ2f37fPcc89h69atiIuLw/r16/Hd73530Mx//OMfUV5e\njpSUFOh0Orz88sseHBEiIiKSk+QYxdX4Ozo6MGHCBNx7770XXUaKyBcWLlyI9PR0vPDCC3JHIRIe\n52wiGstGtdL76KOPYubMmT65YxbRYBwOB++iRuQkztlENJa5XfQeP34cDQ0NyMvLY9FBspEkib/A\niZzAOZuIxjq3i977778fDz30kAejELnugw8+YGsDkRM4ZxPRWOdW0fvOO+8gOzsbKSkpXDEgIhIc\n52wiIkDlzid99tlnePPNN/H222+jsbERCoUCSUlJuOWWW/r2eemNt5EUGeqxoEREvmQymXDdddfJ\nHcMjnJmzX331VcTExMiYkojIfc7M2aO6egMAbNiwAWFhYfjJT37S7+MFBQW45JJLRvPUHvH000/j\nrrvukjsGcwiWgTmYYyQHDhzAlVdeKXcMjxN9zh5IlPNhIFFzAeJmYy7XMJdrnJmzeZ1eIiIiIgp4\nbrU3XOhXv/qVJ3J4jcFgkDsCAOYQLQPAHAMxx9gg+pw9kKjng6i5AHGzMZdrmMvzAn6l98I7c8mJ\nOcTKADDHQMxBIhL1fBA1FyBuNuZyDXN5nleL3m1HGlDfafbmIUbU0NAg6/HPYw6xMgDMMRBzkDfY\n7A40dJlxurkHh+s68WlFG/ZWtuHI2S5UtZnQZrIOe0UJUc8HUXMB4mZjLtcwl+eNur1hOFE6NW5/\n/Qj+fstUhAd59VBERCQzh8OBqrZefFHTgZNNPTjV3IOKFhN0agXCtCqEapUI1SjhANBusqK914Y2\nkxUqhYSJsTpMigvBpPgQ5CSEQqXgTWeIyLNGffWGoRQUFODd5kj8/IpUaFUB30VBRAEmUK/eMBR3\nr95gttlRVNWB/Wfasb+qHTaHA3nJYciK0SEjKhjpUcEI0SiHfY6GLjOO1nfhWH03Suo6UddhxhUZ\nkbg6KwpZMcG86yIRjciZOdury69BagULXiKiAGN3OHC4rhMFZS0oLG9FWmQw5hnC8etvZCA1Isjl\nIjU2RIPYdA0WpEcCAGrae1FQ1oxHd52GRqXArdPjsSA9Ekqu/hLRKHi1IhXhzj+FhYVyRwDAHKJl\nAJhjIOagkTR0mfHygVr8v1eP4JlPq5AcrsXmZRPxxJIs3DgtHmmRnlmVTQrX4nuXJOLF5ZNhDG3F\n21824s63jmHPqRbYBfi9Aoh9noqajblcw1ye59WV3hid2ptPT0REXmazO/DZmXb861gjjtR34fL0\nSDx4VTqyYnReP7YkSRgfYsOKxVn4vKoDWw7U4tXis1hrNCA71vvHJ6LA4tWe3i1Vodj0Lf+9tAUR\njV1jvae3scuMHcebsON4E2JD1Lh2YgwWpEcgWD18f643ORwOFJS14Ll91bgyMxIr8hJlzUNE4vDq\nHdmampowa9YsTJ8+Hbm5uXjttdcu2ketZP8VEZEInJmzzTY7PjrdivXvncSdbx1Da48Vv16cgae+\nNQHfyI6WvcCUJAlXZUXhuW9PRKvJilVvHsMXNR2yZiIi/+F20avX67Fnzx4cPHgQu3btwt133w27\n3d5vn4kCvPwkSu8Jc4iVAWCOgZgjsDkzZ9/29y+x7UgDFmRE4OWbp+CHl6ZgfLS88/hg50NEsBo/\nvyINP7x0HB7fXY4tRbWw2X3b6yvyeSpqNuZyDXN5nts9vSqVCirVuU9vaWmBVqu9aJ8eq/2ijxER\nke85M2f/4b+ykRh28cdFNTtFj6evn4iNH5Tj/n+X4b4r0hDF95IQ0RBG1dPb2dmJefPm4eTJk3jl\nlVdw/fXX9z1WUFCAzaeC8ecbJ3kkKBGRLwViT+9Ic7Y71+kVgc3uwP99UYcdx5vwi0VpyEkIlTsS\nEfmYV3t6ASA0NBQlJSU4cOAA1q1bh66urn6Pq5SSz19yIiKiwY00Z/srpULCirxE/OQyAx5+/zT+\nc6JJ7khEJCCPXLJs4sSJSE1NxdGjRzFz5sy+j1dXnMbG32yDEnbo9Xrk5OTAaDQC+LonxNvb5z/m\nq+MNtb1582ZZvn4Rx2NgFrnGo6SkBKtXr5bt+BwPscajpKQEbW1tAIDKykrk5+cjUA01Z69ZswYG\ngwEAZJuzRztn/e66PKz/zyl8+mUZFsVYcNllgT2ni/QzNNK2KHMMxyswxsudOdvt9oaamhpotVpE\nR0ejrq4OM2fORHFxMaKjowGce6lsvzkOq+Yky3oLycLCwr5BkhNziJWBOZhjJIHW3uDMnC1ie4M7\n50NrjwUb3j+NyGA1fnZFKoK8cGdQUc7TwYiajblcw1yucWbOdrvo3bt3L1atWgXg3LUTH3jgAdx0\n0019jxcUFGBnexR+MHcc9EEqdw5BRCSbQCt6nZmzRSx63WW22fHEh5Vo6DRjw+IMhGn5e4gokDkz\nZ7s9C8ydOxeHDh0adp9glRJtPVYWvUREMnNmzg4kGqUCP78iFX/aV41175bisW9mIjqEV3YgGss8\n/5rPBean6VFQ1uzNQ4zowt4TOTGHWBkA5hiIOUhEozkfFJKEH8xJxsLMSNzz7glUt5mEyOVtomZj\nLtcwl+d5tejVKCWkRwV78xBERERDkiQJN+cm4JbpCbh3eylONnXLHYmIZDKq6/QOp6CgAJNycvHM\np1W4d0GqNw5BROQ1gdbTO5JA6+kdzIenW/C/H1dhw+IMTIoLkTsOEXmQ16/TO5LyFhM6e23ePAQR\nEZFTFqRH4t4FBjz4n1MorumQOw4R+ZhXi96MqGC0mazoschX+IrSe8IcYmUAmGMg5iARefp8mGPQ\n4xeL0vDIrnJ8dqbN7ecR+TwVNRtzuYa5PM+rRa9WpUCUTg2T1e7NwxARETltRlIYHro6Hb/dU4mP\ny1vljkNEPuLVnt5LLrkEP3jrGO6ePw5TeS90IvIj7OkNfCcauvHAeyexZt44XDE+Uu44RDQKsvf0\nAkBWTDDae63ePgwREZFLsmN12HhNJp7dV4WdpU1yxyEiL3O76K2urobRaMTUqVORl5eH999/f9D9\nbs5NwIen5Hv5SJTeE+YQKwPAHAMxR2Bzds4WjbfPh4zoYPzmmiz8dX8tth9rdPrzRD5PRc3GXK5h\nLs9z+1ZparUamzdvRk5ODiorKzF//nxUVVVdtF9rjwXROt4Fh4hITs7O2WORITIIv70uC/ftKEOv\n1Y4bpsbJHYmIvMBjPb1xcXGorq6GWn2uwL2wP+y3eyrw08t5rV4i8h+B3tM73Jw9VtV3mvHzf5Xh\nqqwo3Do9HpIkyR2JiJzks57e9957D3l5eX2T50UHkYCmLosnDkVERKM00pw9VsWFavDEkizsOdWC\nF/bXwEvv8yYimbjd3nBeXV0d1q1bh23btl302Jo1a2AwGNCm0mP9UQNunaKH0WgE8HVPiLe3z3/M\nV8cbavv8y4pyHV+k8RiYRa7xKCkpwerVq2U7PsdDrPEoKSlBW9u567ZWVlYiPz8fgciZORLine truncated
"text": [
"<matplotlib.figure.Figure at 0x49bc6d0>"
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This result may be somewhat suprising to you. The transfer function looks \"fairly\" linear - it is pretty close to a straight line, but the probability distribution of the output is completely different from a "
]
}
],
"metadata": {}
}
]
}
+28
View File
@@ -0,0 +1,28 @@
{
"metadata": {
"name": "",
"signature": "sha256:6f62f259e5289b9332939632d63f121b9911a52f22ec7c03038d30998a233885"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With small exceptions the only mathematics that we truly know how to do is linear algebra, solving\n",
"$$ A\\mathbf{x}=\\mathbf{b}$$\n",
"\n",
"All the nonlinear forms that you have dealt with have been very special. We know how to integrate a polynomial, for example, and so we are able to find the closed form equation for distance given velocity and time:\n",
"$$\\int{(vt+v_0)}\\,dt = \\frac{a}{2}t^2+v_0t+d_0$$\n",
"\n",
"But recall when you leaned that. You learned rules for integrating sums, a different rule for integrating variables raised to a constant power, and so on. And eventually you were handed large tables of integrals. Over the centuries different mathematicians discovered how to integrate differernt equations, and our knowledge of how to integrate was slowly built up. But if I was to give you an arbitrary equation - truly arbitrary, not a polynomial "
]
}
],
"metadata": {}
}
]
}
+104
View File
@@ -0,0 +1,104 @@
# -*- coding: utf-8 -*-
"""
Created on Sun May 18 11:09:23 2014
@author: rlabbe
"""
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import normal
def plot_transfer_func(data, f, lims,num_bins=1000):
ys = f(data)
#plot output
plt.subplot(2,2,1)
plt.hist(ys, num_bins, orientation='horizontal',histtype='step')
plt.ylim(lims)
plt.gca().xaxis.set_ticklabels([])
# plot transfer function
plt.subplot(2,2,2)
x = np.arange(lims[0], lims[1],0.1)
y = f(x)
plt.plot (x,y)
isct = f(0)
plt.plot([0,0,lims[0]],[lims[0],isct,isct],c='r')
plt.xlim(lims)
# plot input
plt.subplot(2,2,4)
plt.hist(data, num_bins, histtype='step')
plt.xlim(lims)
plt.gca().yaxis.set_ticklabels([])
plt.show()
normals = normal(loc=0.0, scale=1, size=5000000)
#rint h(normals).sort()
def f(x):
return 2*x + 1
def g(x):
return (cos(4*(x/2+0.7)))*sin(0.3*x)-0.9*x
return (cos(4*(x/3+0.7)))*sin(0.3*x)-0.9*x
#return -x+1.2*np.sin(0.7*x)+3
return sin(5-.2*x)
def h(x): return cos(.4*x)*x
plot_transfer_func (normals, g, lims=(-4,4),num_bins=500)
del(normals)
#plt.plot(g(np.arange(-10,10,0.1)))
'''
ys = f(normals)
r = np.linspace (min(normals), max(normals), num_bins)
h= np.histogram(ys, num_bins,density=True)
print h
print len(h[0]), len(h[1][0:-1])
#plot output
plt.subplot(2,2,1)
h = np.histogram(ys, num_bins,normed=True)
p, = plt.plot(h[0],h[1][1:])
plt.ylim((-10,10))
plt.xlim((max(h[0]),0))
# plot transfer function
plt.subplot(2,2,2)
x = np.arange(-10,10)
y = 1.2*x + 1
plt.plot (x,y)
plt.plot([0,0],[-10,f(0)],c='r')
plt.ylim((-10,10))
# plot input
plt.subplot(2,2,4)
h = np.histogram(normals, num_bins,density=True)
plt.plot(h[1][1:],h[0])
plt.xlim((-10,10))
plt.show()
'''
+103
View File
@@ -0,0 +1,103 @@
# -*- coding: utf-8 -*-
"""
Created on Sun May 18 11:09:23 2014
@author: rlabbe
"""
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from numpy.random import normal
def plot_transfer_func(data, f, lims,num_bins=1000):
ys = f(data)
#plot output
plt.subplot(2,2,1)
plt.hist(ys, num_bins, orientation='horizontal',histtype='step')
plt.ylim(lims)
plt.gca().xaxis.set_ticklabels([])
plt.title('output')
# plot transfer function
plt.subplot(2,2,2)
x = np.arange(lims[0], lims[1],0.1)
y = f(x)
plt.plot (x,y)
isct = f(0)
plt.plot([0,0,lims[0]],[lims[0],isct,isct],c='r')
plt.xlim(lims)
plt.ylim(lims)
# plot input
plt.subplot(2,2,4)
plt.hist(data, num_bins, histtype='step')
plt.xlim(lims)
plt.gca().yaxis.set_ticklabels([])
plt.title('input')
plt.show()
'''
normals = normal(loc=0.0, scale=1, size=5000000)
#rint h(normals).sort()
def f(x):
return 2*x + 1
def g(x):
return (cos(4*(x/2+0.7)))*sin(0.3*x)-0.9*x
return (cos(4*(x/3+0.7)))*sin(0.3*x)-0.9*x
#return -x+1.2*np.sin(0.7*x)+3
return sin(5-.2*x)
def h(x): return cos(.4*x)*x
plot_transfer_func (normals, g, lims=(-4,4),num_bins=500)
del(normals)
#plt.plot(g(np.arange(-10,10,0.1)))
'''
'''
ys = f(normals)
r = np.linspace (min(normals), max(normals), num_bins)
h= np.histogram(ys, num_bins,density=True)
print h
print len(h[0]), len(h[1][0:-1])
#plot output
plt.subplot(2,2,1)
h = np.histogram(ys, num_bins,normed=True)
p, = plt.plot(h[0],h[1][1:])
plt.ylim((-10,10))
plt.xlim((max(h[0]),0))
# plot transfer function
plt.subplot(2,2,2)
x = np.arange(-10,10)
y = 1.2*x + 1
plt.plot (x,y)
plt.plot([0,0],[-10,f(0)],c='r')
plt.ylim((-10,10))
# plot input
plt.subplot(2,2,4)
h = np.histogram(normals, num_bins,density=True)
plt.plot(h[1][1:],h[0])
plt.xlim((-10,10))
plt.show()
'''