Fixed code that explicitly calculated weight. It was buggy.
This commit is contained in:
Roger Labbe
parent
b995668db1
commit
87e595078b
136
exp/gh.py
Normal file
136
exp/gh.py
Normal file
@ -0,0 +1,136 @@
|
||||
# -*- coding: utf-8 -*-
|
||||
"""
|
||||
Created on Thu May 15 16:07:26 2014
|
||||
|
||||
@author: RL
|
||||
"""
|
||||
from __future__ import division
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy.random as random
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def g_h_filter (data, x, dx, g, h, dt=1.):
|
||||
results = []
|
||||
for z in data:
|
||||
x_est = x + (dx*dt)
|
||||
residual = z - x_est
|
||||
|
||||
dx = dx + h * (residual / float(dt))
|
||||
x = x_est + g * residual
|
||||
print('gx',dx,)
|
||||
|
||||
results.append(x)
|
||||
|
||||
|
||||
return results
|
||||
|
||||
|
||||
'''
|
||||
computation of x
|
||||
x_est = weight + gain
|
||||
residual = z - weight - gain
|
||||
x = weight + gain + g (z - weight - gain)
|
||||
|
||||
w + gain + gz -wg -ggain
|
||||
w -wg + gain - ggain + gz
|
||||
|
||||
w(1-g) + gain(1-g) +gz
|
||||
|
||||
(w+g)(1-g) +gz
|
||||
|
||||
'''
|
||||
'''
|
||||
gain computation
|
||||
|
||||
gain = gain + h/t* (z - weight - gain)
|
||||
= gain + hz/t -hweight/t - hgain/t
|
||||
|
||||
= gain(1-h/t) + h/t(z-weight)
|
||||
|
||||
'''
|
||||
'''
|
||||
gain+ h*(z-w -gain*t)/t
|
||||
|
||||
gain + hz/t -hw/t -hgain
|
||||
|
||||
gain*(1-h) + h/t(z-w)
|
||||
|
||||
|
||||
'''
|
||||
def weight2():
|
||||
w = 0
|
||||
gain = 200
|
||||
t = 10.
|
||||
weight_scale = 1./6
|
||||
gain_scale = 1./10
|
||||
|
||||
weights=[2060]
|
||||
for i in range (len(weights)):
|
||||
z = weights[i]
|
||||
w_pre = w + gain*t
|
||||
|
||||
new_w = w_pre * (1-weight_scale) + z * (weight_scale)
|
||||
|
||||
print('new_w',new_w)
|
||||
|
||||
gain = gain *(1-gain_scale) + (z - w) * gain_scale/t
|
||||
|
||||
print (z)
|
||||
print(w)
|
||||
|
||||
#gain = new_gain * (gain_scale) + gain * (1-gain_scale)
|
||||
w = new_w
|
||||
print ('w',w,)
|
||||
print ('gain=',gain)
|
||||
|
||||
|
||||
def weight3():
|
||||
w = 160.
|
||||
gain = 1.
|
||||
t = 1.
|
||||
weight_scale = 6/10.
|
||||
gain_scale = 2./3
|
||||
|
||||
weights=[158]
|
||||
for i in range (len(weights)):
|
||||
z = weights[i]
|
||||
w_pre = w + gain*t
|
||||
|
||||
new_w = w_pre * (1-weight_scale) + z * (weight_scale)
|
||||
|
||||
print('new_w',new_w)
|
||||
|
||||
gain = gain *(1-gain_scale) + (z - w) * gain_scale/t
|
||||
|
||||
print (z)
|
||||
print(w)
|
||||
|
||||
#gain = new_gain * (gain_scale) + gain * (1-gain_scale)
|
||||
w = new_w
|
||||
print ('w',w,)
|
||||
print ('gain=',gain)
|
||||
weight3()
|
||||
'''
|
||||
#zs = [i + random.randn()*50 for i in range(200)]
|
||||
zs = [158.0, 164.2, 160.3, 159.9, 162.1, 164.6, 169.6, 167.4, 166.4, 171.0]
|
||||
|
||||
#zs = [2060]
|
||||
data= g_h_filter(zs, 160, 1, .6, 0, 1.)
|
||||
|
||||
'''
|
||||
|
||||
data = g_h_filter([2060], x=0, dx=200, g=1./6, h = 1./10, dt=10)
|
||||
print data
|
||||
|
||||
|
||||
'''
|
||||
print data
|
||||
print data2
|
||||
plt.plot(data)
|
||||
plt.plot(zs, 'g')
|
||||
plt.show()
|
||||
'''
|
||||
108
g-h_filter.ipynb
108
g-h_filter.ipynb
@ -1,7 +1,7 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:221bed83659f3d65f886890673eb366441306ce17460af258c7b7626ee0dfda6"
|
||||
"signature": "sha256:79e48eadc10012625737389a3259ac9df692b99636f47612d5b5a9fdde69923f"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
@ -212,7 +212,7 @@
|
||||
"\n",
|
||||
"Should it be half way? Maybe, but in general it seems like we might know that our prediction is very accurate or very inaccurate. Probably the accuracy of our prediction differs from the accuracy of the scale. Recall what we did when A was much more accurate than B - we scaled the answer to be closer to A than B. \n",
|
||||
"\n",
|
||||
"Let's try a randomly chosen number: $\\frac{6}{10}$. Our estimate will be six tenth the prediction and the rest will be from the measurement. Let's just code that up and see the result."
|
||||
"Let's try a randomly chosen number: $\\frac{6}{10}$. Our estimate will be six tenth the measurement and the rest will be from the prediction. Let's just code that up and see the result. We have to take into account one other factor. Weight gain has units of lbs/time, so to be general we will need to add a time step $t$, which we will set to 1 (day)."
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -221,6 +221,7 @@
|
||||
"input": [
|
||||
"w = 160\n",
|
||||
"gain = 1\n",
|
||||
"t = 1\n",
|
||||
"scale = 6/10\n",
|
||||
"\n",
|
||||
"# we are storing the estimates so we can use them in the next problem.\n",
|
||||
@ -228,9 +229,11 @@
|
||||
"\n",
|
||||
"for i in range (1, len(weights)):\n",
|
||||
" z = weights[i] # most filter literature uses 'z' for measurements\n",
|
||||
" w = (w+gain) * scale + (1-scale) * z\n",
|
||||
" w = (w+gain*t) * (1-scale) + (scale) * z\n",
|
||||
" estimates.append(w)\n",
|
||||
" plt.scatter(i, w)\n",
|
||||
"\n",
|
||||
"plt.xlim([0,10])\n",
|
||||
"plt.plot(weights, c='r')\n",
|
||||
"plt.plot([0,9],[160,170],c='g')\n",
|
||||
"plt.show()"
|
||||
@ -258,17 +261,18 @@
|
||||
"import numpy.random as random\n",
|
||||
"w = 160\n",
|
||||
"gain = 1\n",
|
||||
"t=1\n",
|
||||
"weight_scale = 6/10\n",
|
||||
"gain_scale = 2/3\n",
|
||||
"\n",
|
||||
"for i in range (1, len(weights)):\n",
|
||||
"for i in range (len(weights)):\n",
|
||||
" z = weights[i]\n",
|
||||
" new_w = (w+gain) * weight_scale + z * (1-weight_scale)\n",
|
||||
" new_gain = new_w - w\n",
|
||||
" gain = new_gain * gain_scale + gain * (1-gain_scale)\n",
|
||||
" new_w = (w+gain*t) * (1-weight_scale) + z * weight_scale\n",
|
||||
" gain = gain *(1-gain_scale) + (z - w) * gain_scale/t\n",
|
||||
" w = new_w\n",
|
||||
" plt.scatter(i, w,c='b')\n",
|
||||
"\n",
|
||||
"plt.xlim([0,10])\n",
|
||||
"plt.plot(weights, c='r')\n",
|
||||
"plt.plot([0,9],[160,170],c='g')\n",
|
||||
"plt.show()"
|
||||
@ -329,26 +333,28 @@
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"def g_h_filter (data, x0, dx, g, h):\n",
|
||||
"def g_h_filter (data, x0, dx, g, h, dt=1.):\n",
|
||||
" x = x0\n",
|
||||
" results = []\n",
|
||||
" for z in data:\n",
|
||||
" x_est = (x+dx) * g + z * (1-g)\n",
|
||||
" dx_est = x_est - x\n",
|
||||
" dx = dx_est * gain_scale + dx * (1-gain_scale)\n",
|
||||
" x = x_est\n",
|
||||
" results.append(x_est)\n",
|
||||
" x_est = x + (dx*dt)\n",
|
||||
" residual = z - x_est\n",
|
||||
" dx = dx + h * (residual) / dt\n",
|
||||
" x = x_est + g * residual\n",
|
||||
" \n",
|
||||
" results.append(x)\n",
|
||||
" \n",
|
||||
" return results\n",
|
||||
"\n",
|
||||
"def plot_g_h_results (measurements, filtered_data):\n",
|
||||
" \n",
|
||||
" plt.scatter (range(len(filtered_data)), filtered_data)\n",
|
||||
" plt.plot(measurements)\n",
|
||||
"def plot_g_h_results (measurements, filtered_data): \n",
|
||||
" plt.scatter (range(len(filtered_data)), filtered_data, c='b')\n",
|
||||
" plt.plot(measurements,c='r')\n",
|
||||
" plt.show()\n",
|
||||
" \n",
|
||||
"\n",
|
||||
"data = g_h_filter (weights, 160., 1., 6/10, 2/3)\n",
|
||||
"\n",
|
||||
"plt.xlim([0,10])\n",
|
||||
"plt.plot([0,9],[160,170],c='g')\n",
|
||||
"data = g_h_filter (data=weights, x0=160, dx=1, g=6./10, h = 2./3, dt=1.)\n",
|
||||
"plot_g_h_results (weights, data)"
|
||||
],
|
||||
"language": "python",
|
||||
@ -356,6 +362,23 @@
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Note that I rewrote the equations somewhat:\n",
|
||||
"\n",
|
||||
"$$\\begin{align*}\n",
|
||||
"\\hat{x}_t &= gz + \\hat{x}_{t-1} (1-g) \\\\\n",
|
||||
"&= gz + \\hat{x}_{t-1} - g\\hat{x}_{t-1} \\\\\n",
|
||||
"&= \\hat{x}_{t-1} + g*(z-\\hat{x}_{t-1}) \\\\\n",
|
||||
"&= \\hat{x}_{t-1} + g*residual\n",
|
||||
"\\end{align*}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"I'll let you decide which form of the equation is more expressive. One form explicitly uses $g$ and $(1-g)$ to compute the point between the two values, the other finds the difference between the points, and adds a fraction of that to the first value. Both are computing the same thing. You'll see both forms in the literature, so I have used both forms in my code here."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@ -379,24 +402,6 @@
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I am using a Python idiom that might be unfamiliar to you. You can skip to the next exercise if you understand this code, otherwise here is a quick description.\n",
|
||||
"\n",
|
||||
"This code uses an idiom called *list comprehensions*. This allows you insert code right inside of your list; when encountered, the code is run and the list will contain the data generated from the code. Let's look at a simple example:\n",
|
||||
"\n",
|
||||
" [i for i in range(10)]\n",
|
||||
" \n",
|
||||
"The result of this code is the list [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]. The general structure is\n",
|
||||
" [expr for_loop]\n",
|
||||
" \n",
|
||||
"where 'for_loop' is a for loop that uses some \n",
|
||||
"\n",
|
||||
"**blah blah fill this in**"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@ -411,7 +416,7 @@
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"zs = gen_data (x0=5, dx=2, count=100, noise_factor=10)\n",
|
||||
"data = g_h_filter (data=zs, x0=100., dx=2., g=0.8, h=0.02)\n",
|
||||
"data = g_h_filter (data=zs, x0=100., dx=2., dt=1.,g=0.2, h=0.01)\n",
|
||||
"plot_g_h_results (measurements=zs, filtered_data=data)"
|
||||
],
|
||||
"language": "python",
|
||||
@ -443,8 +448,7 @@
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
@ -484,8 +488,7 @@
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
@ -519,8 +522,7 @@
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@ -528,8 +530,7 @@
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@ -537,8 +538,7 @@
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@ -546,8 +546,7 @@
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@ -555,8 +554,7 @@
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
@ -580,8 +578,7 @@
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@ -596,8 +593,7 @@
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": ""
|
||||
"outputs": []
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user