Setting a polynomial using np.polyfit in 3 dimensions

Question

I have a data array with dimensions (N, 3) for some integer N that defines the trajectory of the particle in three-dimensional space, i.e. each row entry is the coordinates (x, y, z) of the particle. This trajectory is smooth and uncomplicated, and I want to be able to match the polynomial to this data.

I can only do this with coordinates (x, y) using np.polyfit :

import numpy as np #Load the data some_file = 'import_file.txt' data = np.loadtxt(some_file) x = data[:,0] y = data[:,1] #Fit a 4th order polynomial fit = np.polyfit(x,y,4) 

It gives me polynomial coefficients, no problem.

How would I extend this to my case when I need a polynomial that describes the coordinates x, y, z?

Answer 1

You have several options. First, let me expand your 2D case fit = np.polyfit(x,y,4) . This means that you describe the position of the particle y as a function of x. This is normal since it will not move back to x. (That is, it can have a single y value for each x). Since the motion in space is decomposed into three independent coordinates, we can place the coordinates independently to get a 3D model:

 fitxy = np.polyfit(x, y, 4) fitxz = np.polyfit(x, z, 4) 

Now both y and z are a polynomial function of x. As mentioned earlier, this has the disadvantage that the particle can move monotonously along x.

A real physical particle will not behave like this. They usually bounce in all three dimensions, walking the way they like. However, there is a 4th dimension in which they never turn: time.

So add the time:

 t = np.arange(data.shape[0]) # simple assumption that data was sampled in regular steps fitx = np.polyfit(t, x, 4) fity = np.polyfit(t, y, 4) fitz = np.polyfit(t, z, 4) 

Now the particle is modeled for free movement in space, as a function of time.

Answer 2

You can do this with sklearn using PolynomialFeatures .

For example, consider the following code:

 from sklearn.preprocessing import PolynomialFeatures from sklearn.linear_model import LinearRegression X = np.random.rand(100,2) y=X[:,0]**2-3*X[:,1]**3+2*X[:,0]*X[:,1]-5 poly = PolynomialFeatures(degree=3) X_t = poly.fit_transform(X) clf = LinearRegression() clf.fit(X_t, y) print(clf.coef_) print(clf.intercept_) 

he prints

 [ 0.00000000e+00 1.08482012e-15 9.65543103e-16 1.00000000e+00 2.00000000e+00 -1.18336405e-15 2.06115185e-15 1.82058329e-15 2.33420247e-15 -3.00000000e+00] -5.0 

which are exactly the coefficients.