1

I'm trying to make a 3 dimensional 3rd degree polynomial fit using scipy curve_fit, but since I don't have much experience with python curve fitting tools I'm having some trouble doing that.

I have two arguments for a function which are vectors, and a matrix representing the value at each point. I'm trying to do a 3rd degree polynomial fit to each argument but I'm having some troubles. I'm trying to use the curve_fit function because thats what I found online, if there are any better suggestions like polyfit I'll gladly hear them. Here is my code:

import numpy as np
from scipy.optimize import curve_fit

def polyfit(data,a1,a2,a3,a4,a5,a6,a7):
    return a1*np.power(data[0],3) + a2*np.power(data[0],2) + a3*data[0] \
        + a4*np.power(data[1],3) + a5*np.power(data[1],2) + a6*data[1] + a7

xtemp = np.array([10,25,35,45,55])
yexpT = np.array([1,5,10,100,500,1000,5000,7500])
dataMat = np.array([[1475,1475,1478,1528,1776,2086,4588,6146] ,
            [1575,1575,1578,1728,1976,2086,6588,7146],
            [1675,1675,1778,1928,2176,2586,7588,9146],
            [2575,2575,2578,2728,2976,3086,9588,11146],
            [3575,3575,3578,3728,4976,8086,12588,15146]]) 
ymesh,xmesh = np.meshgrid(yexpT,xtemp)

popt,pcov = curve_fit(polyfit,[xmesh,ymesh],dataMat)

As you can see, the xtemp is the number of lines and yexpT is the number of columns.

I tried to find help for example in the following link: Python surface fitting of variables of different dimensionto get unknown parameters? Yet still to no avail.

Thanks in advance to anyone who is able to help.

mashtock
  • 400
  • 4
  • 21
  • Id rather use `scipy.optimize.least_squares()`...with using `numpy.meshgrid()` Checking how to implement this, likely gives you the hint to solve your problem. – mikuszefski Feb 03 '22 at 14:09
  • @mikuszefski I don't mind as long as you manage to get a fitted function. I think I found someone who did a really complexed kind of solution. If it's good then i'll point it out. It seems like my question is a really complex one to do in python compare to matlab. – mashtock Feb 03 '22 at 14:11
  • The only thing is to see how to define the function for residuals. like def residuals( params, xlist, ylist, zmatrix):... – mikuszefski Feb 03 '22 at 14:14
  • @mikuszefski you can take for each line or each column the std or the mad if you'd like. – mashtock Feb 03 '22 at 14:17
  • I remember an online example how to directly use curve_fit on 3d data, but cannot find it right now. Don't have a version on my disk, as I use `least_squares`. I'll try to give the main lines of thought as an answer – mikuszefski Feb 03 '22 at 14:24
  • ...no mixed terms? – mikuszefski Feb 03 '22 at 14:30
  • @mikuszefski no mixed terms, thanks for the help btw! I try to find something that will work here: https://stackoverflow.com/questions/7997152/python-3d-polynomial-surface-fit-order-dependent/32297563#32297563 – mashtock Feb 03 '22 at 14:34

2 Answers2

1

This is how a simple solution with least_squares would look like:

import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
import numpy as np
from scipy.optimize import least_squares


def poly( x, y,
    c0,
    a1, a2, a3,
    b1, b2, b3
):
    "non-mixing test plynomial"
    out = c0
    out += a1 * x**1 + a2 * x**2 + a3 * x**3
    out += b1 * y**1 + b2 * y**2 + b3 * y**3
    return out



### test grid
xL = np.linspace( -2, 5, 17 )
yL = np.linspace( -0, 6, 16 )

XX, YY = np.meshgrid(  xL, yL )

### test data
ZZ = poly(
    XX, YY,
    1.1,
    2.2, -0.55, 0.09,
    -3.1, +0.75, -0.071
)


### test noise
znoise = np.random.normal( size=( len(xL) * len(yL) ), scale=0.5 )
znoise.resize( len(yL), len(xL) )
ZZ += znoise


### residual finction, i.e. the important stuff
def res( params, X, Y, Z ):
    th = poly( X, Y, *params )
    diff = Z - th
    return np.concatenate( diff )

### fit
sol = least_squares( res, x0=7*[0], args=( XX, YY, ZZ ) )
print( sol.x )

### plot it
fig = plt.figure( )
ax = fig.add_subplot( 1, 1, 1, projection="3d" )
ax.scatter( XX, YY, ZZ, color='r' )
ax.plot_surface( XX, YY, poly( XX, YY, *(sol.x) ) )

plt.show()
mikuszefski
  • 3,943
  • 1
  • 25
  • 38
  • Wow maybe because I'm stupid or tired, the least square function implementation in python looks really complicated. I'll try to understand better your solution but it's the end of the day for me and it's a weekend so it might take some time :P. Really thank you though! – mashtock Feb 03 '22 at 15:05
  • @mashtock Don't worry, I am sure, the moment you look at it after a good rest, you'll notice that it is rather easy. In caste of question, do not hesitate. Cheers – mikuszefski Feb 03 '22 at 15:16
  • Weekend is over for me and I tried to run your solution and thanks again for your solution it seems like it's working. Unfortunately for my data set the fit was far from good :(. So I'll have to try something else. By something else I mean that I'll try to find a fit without giving the poly model myself. – mashtock Feb 06 '22 at 07:13
  • @mashtock Sorry to hear that. Hve to confess that now is the first time I really lokk at your data in detail. Looks more like a jump happening between 1000 and 5000, which is likely to move as a function of x. Unfortunately, no data in that region. I guess there are plenty of functions fitting this, but a justified model would really help. Cheers – mikuszefski Feb 07 '22 at 06:43
0

BTW, this is a pure linear fit of course. So with the definitions from my other post, one could just do:

### fitting:
X = np.concatenate( XX )
Y = np.concatenate( YY )
Z = np.concatenate( ZZ )

VT = np.array( [
    np.ones( len(xL) * len(yL) ),
    X, X**2, X**3,
    Y, Y**2, Y**3
] )
V = np.transpose( VT )
eta = np.dot( VT, Z )
A = np.dot( VT, V )
sol = np.linalg.solve( A, eta )
print( sol )

### plot it
fig = plt.figure( )
ax = fig.add_subplot( 1, 1, 1, projection="3d" )
ax.scatter( XX, YY, ZZ, color='r' )
ax.plot_surface( XX, YY, poly( XX, YY, *(sol) ) )
plt.show()

I put a bit more info in here

mikuszefski
  • 3,943
  • 1
  • 25
  • 38