Fitting polynomial curve

[eryksun ()]

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]).