Find slope of function given empirical data.

[Thomas]

Hello all.

Trying to find slope of function using numpy.
Getting close, but results are a bit off. Hope someone out here can
help.

import numpy as np

def deriv(y):
x = list(range(len(y)))
x.reverse() # Change from [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
x = np.array(x) # to [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
y = np.array(y) # x.reverse() is used to put point 0 at end of
list.
z = np.polyfit(x, y, 2)
print np.poly1d(z)
# Returns:
# 2
# 3.142 x - 18.85 x + 35.13
# 2
# Should be closer to 3.142 x - 6.283 +
10 ????????????????????
return [z[0] * 2, z[1]]

if __name__=='__main__':
# range(-6,5,1)
# y = 3.141592 * x ** 2 - 6.283184 * x + 10 for x in range(-6, 5,
1)
# 160.796416, 119.95572, 85.398208, 57.12388, 35.132736,
19.424776, 10.0, 6.858408, 10.0, 19.424776, 35.132736
#
# y' = 6.283184 * x - 6.283184 for x in range(-6, 5, 1)
# -43.982288, -37.699104, -31.41592, -25.132736, -18.849552,
-12.566368, -6.283184, 0.0, 6.283184, 12.566368, 18.849552
#
z = deriv([160.796416, 119.95572, 85.398208, 57.12388, 35.132736,
19.424776, 10.0, 6.858408, 10.0, 19.424776, 35.132736])
for x in range(-6,5,1):
print str(w(x)) + ',' ,
# Returns:
# -56.548656, -50.265472, -43.982288, -37.699104, -31.41592,
-25.132736, -18.849552, -12.566368, -6.283184, -1.06581410364e-14,
6.283184
# Should be:
# -43.982288, -37.699104, -31.41592, -25.132736, -18.849552,
-12.566368, -6.283184, 0.0, 6.283184, 12.566368, 18.849552
# Note that the range is offset by 2 positions

 

[duncan smith]

Hello all.

Trying to find slope of function using numpy.
Getting close, but results are a bit off. Hope someone out here can
help.

[snip]

Why are you generating y-coordinates from the x-coordinates [-6, -5, -4,
-3, -2, -1, 0, 1, 2, 3, 4]? If you're going to use the x-coordinates
[10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0] in your function, then use the same
to generate the y-coordinates. Surely, if you have empirical data
(which, for some reason, you know are well fitted by a quadratic
function?) you'd pass both the x and y coordinates to the function?
Maybe (untested),

def deriv(x, y):
z = np.polyfit(x, y, 2)
p = np.poly1d(z)
return p.deriv()

Duncan

 

[Peter Otten]

Trying to find slope of function using numpy.
Getting close, but results are a bit off. Hope someone out here can
help.

You don't make it easy to understand your post. In the future please try to
rely more on plain english than on lots of numbers and code that doesn't
run.

import numpy as np

def deriv(y):
x = list(range(len(y)))
x.reverse() # Change from [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
x = np.array(x) # to [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
y = np.array(y) # x.reverse() is used to put point 0 at end of
list.
z = np.polyfit(x, y, 2)
print np.poly1d(z)
# Returns:
# 2
# 3.142 x - 18.85 x + 35.13
# 2
# Should be closer to 3.142 x - 6.283 +
10 ????????????????????

To add one more question mark: how did you find that alternative?

Anyway, we can put both polynomials to a test:

import numpy as np
y = np.array([160.796416, 119.95572, 85.398208, 57.12388, 35.132736,19.424776, 10.0, 6.858408, 10.0, 19.424776, 35.132736])
x = np.arange(len(y), dtype=float)[::-1]
p1 = np.poly1d(np.polyfit(x, y, 2))
print p1

2
3.142 x - 18.85 x + 35.13

p2 = np.poly1d([3.142, -6.283, 10.0])
print p2

2
3.142 x - 6.283 x + 10

Now calculate the sum of the squares:
33028.342907811333

Conclusion: numpy's result is much better than what you suggest.

Peter

 

[Thomas]

Trying to find slope of function using numpy.
Getting close, but results are a bit off. Hope someone out here can
help.

You don't make it easy to understand your post. In the future please try to
rely more on plain english than on lots of numbers and code that doesn't
run.

import numpy as np

def deriv(y):
    x = list(range(len(y)))
    x.reverse()     # Change from [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    x = np.array(x) #        to   [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
    y = np.array(y) # x.reverse() is used to put point 0 at end of
list.
    z = np.polyfit(x, y, 2)
    print np.poly1d(z)
    #  Returns:
    #         2
    #  3.142 x - 18.85 x + 35.13
    #                              2
    # Should be closer to   3.142 x  - 6.283 +
10   ????????????????????

To add one more question mark: how did you find that alternative?

Anyway, we can put both polynomials to a test:

import numpy as np
y = np.array([160.796416, 119.95572, 85.398208, 57.12388,

35.132736,19.424776, 10.0, 6.858408, 10.0, 19.424776, 35.132736])>>> x = np.arange(len(y), dtype=float)[::-1]
       2
3.142 x - 18.85 x + 35.13>>> p2 = np.poly1d([3.142, -6.283, 10.0])
       2
3.142 x - 6.283 x + 10

Now calculate the sum of the squares:

5.0683524299544787e-26>>> np.sum((p2(x)-y)**2)

33028.342907811333

Conclusion: numpy's result is much better than what you suggest.

Peter

As usual, thanks to all for putting me on the right track.

Kind regards.

 

[Grant Edwards]

Trying to find slope of function using numpy. Getting close, but
results are a bit off. Hope someone out here can help.

import numpy as np

def deriv(y):
x = list(range(len(y)))
x.reverse() # Change from [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
x = np.array(x) # to [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]
y = np.array(y) # x.reverse() is used to put point 0 at end of
list.
z = np.polyfit(x, y, 2)
print np.poly1d(z)
# Returns:
# 2
# 3.142 x - 18.85 x + 35.13
# 2
# Should be closer to 3.142 x - 6.283 +
10 ????????????????????

Numpy's answer is correct.

I've no idea where you got your answer. Perhaps the commented-out
stuff was meant to convey that. If so, it went past me.

Here's the least-squares fit done by gnuplot, the resulting fuction
f(x) plotted against the data, as well as your "should be" function
plotted against the data.

As you can see, your "should be" function is way off.

For prettier results, you can change the gnuplot script to use a
different output format -- I just used the dumb terminal so I could
post it here. In any case, the values returned by numpy and gnuplot
hit the data points pretty much exactly -- far better than your
"should be" results.

NB: You might want to grab a copy of gnuplot (or a similar) tool to
use when tinkering with data analysis and visualisation. It's
very handy.


------------------------------foo.gp------------------------------
#!/usr/bin/gnuplot

f(x) = a*x**2 + b*x + c
fit f(x) "foo.dat" via a,b,c

set xra[-1:11]
set term dumb 120 40

plot f(x), "foo.dat"

plot 3.142*x**2 + 6.283*x + 10, "foo.dat"

------------------------------foo.gp------------------------------


------------------------------foo.dat------------------------------
10 160.796416
9 119.95572
8 85.398208
7 57.12388
6 35.132736
5 19.424776
4 10.0
3 6.858408
2 10.0
1 19.424776
0 35.132736
------------------------------foo.dat------------------------------

--------------------results of running foo.gp--------------------
Iteration 0
[...]
Iteration 1
[...]
Iteration 2
[...]
Iteration 3
[...]
Iteration 4
[...]
Iteration 5
[...]
Iteration 6
[...]
Iteration 7
[...]
**************************
After 8 iterations the fit converged.
final sum of squares of residuals : 7.70717e-28
rel. change during last iteration : 0

degrees of freedom (FIT_NDF) : 8
rms of residuals (FIT_STDFIT) = sqrt(WSSR/ndf) : 9.81528e-15
variance of residuals (reduced chisquare) = WSSR/ndf : 9.63396e-29

Final set of parameters Asymptotic Standard Error
======================= ==========================

a = 3.14159 +/- 3.351e-16 (1.067e-14%)
b = -18.8496 +/- 3.479e-15 (1.846e-14%)
c = 35.1327 +/- 7.478e-15 (2.128e-14%)

[...]


250 ++-------+-----------------+-----------------+------------------+-----------------+-----------------+-------++
| + + + + + f(x) ****** |
| "foo.dat" A |
| |
| |
| |
| *
200 ++ *+
| * |
| *** |
| * |
| * |
| * |
| *A |
150 ++ * ++
| * |
| ** |
| ** |
| ** |
| *A |
| ** |
| * |
100 ++ ** ++
| *** |
| **A |
| ** |
| ** |
| ** |
** **A |
50 ++**** *** ++
| ** ** |
| *A** **A* |
| *** **** |
| ***A* *A** |
| ******** ******** |
| + A********A********A + + + |
0 ++-------+-----------------+-----------------+------------------+-----------------+-----------------+-------++
0 2 4 6 8 10


500 ++-------+-----------------+-----------------+------------------+-----------------+-----------------+-------++
| + + + + 3.142*x**2 + 6.283*x + 10 ****** |
| "foo.dat" A |
| *
450 ++ **+
| * |
| *** |
400 ++ ** ++
| * |
| ** |
| ** |
350 ++ *** ++
| * |
| ** |
300 ++ ** ++
| ** |
| ** |
| *** |
250 ++ ** ++
| ** |
| ** |
| **** |
200 ++ ** ++
| ** |
| *** A |
150 ++ **** ++
| ** |
| *** A |
| **** |
100 ++ **** ++
| **** A |
| ***** |
50 ++ **** A ++
| A ******* A |
| ******* |
| ****************A* A A A + + + |
0 **-------+-----------------+--------A--------+------------------+-----------------+-----------------+-------++
0 2 4 6 8 10

--------------------results of running foo.gp--------------------