On 3/17/2011 1:42 AM, Astan Chee wrote:
>
> I have 2 points in 3D space and a bunch of points in-between them. I'm
> trying to fit a polynomial curve on it. Currently I'm looking through
> numpy but I don't think the function exists to fit a function like this:
> y = ax**4 + bx**3 + cx**2 + dx + e
You can use np.polyfit, which uses np.linalg.lstsq.
For a degree M-1 polynomial and K samples such that K > M, this is an
over-determined least squares problem:
t = [t1, t2, ..., tk] # K sample times
Y.shape == (K, N) # data samples Y in R_N
# design matrix X
# you don't have to calculate this since polyfit does it
# for you internally with vander(t, M)
X = [[x**m for m in range(M-1,-1,-1)] for x in t]
X.shape == (K, M)
# unknown coefficient matrix B
B.shape == (M, N)
Y =~ dot(X, B)
Use polyfit to form the least squares solution. For example, a 4th order fit:
B = np.polyfit(t, Y, 4)
# or for more information
B, r, rankX, sX, rcond = np.polyfit(t, Y, 4, full=True)
where r is the vector of residuals; rankX and sX are the rank and singular
values of the Van der Monde matrix; and rcond was used to threshold the
singular values.
Interpolation:
t2 = np.linspace(t[0], t[-1], len(t)*1000)
Y2 = np.dot(np.vander(t2, M), B)
or
Y2 = np.polyval(B, t2[:, newaxis])
polyval uses Horner's method, which calculates the following:
t = t2[:, newaxis]
y = zeros_like(t)
for i in range(len(B)):
y = y * t + B[i]
return y
y is initially a length len(t) column vector of zeros and B[0] is a length N
row vector, so the first iteration just tiles B[0] in a len(t) by N matrix.
Otherwise it's a normal Horner evaluation.
Other than minimizing the residuals, you can examine the fit visually with a
simple matplotlib plot:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(*Y2.T)
ax.plot(*Y.T, marker='.', markerfacecolor='r')
Or you could create three 2D plots of (t2, Y2[:,i]) and (t, Y[:,i]).
--
http://mail.python.org/mailman/listinfo/python-list
