CODEWITHSUNDEEP
pythonpython-3.xscipycurve-fittingdata-fitting
mashtock
mashtock

Reputation: 400

Multivariant different dimensional fit

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.

Upvotes: 1

Views: 430

Answers (2)

mikuszefski
mikuszefski

Reputation: 4043

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

Upvotes: 0

mikuszefski
mikuszefski

Reputation: 4043

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()

Upvotes: 1

Related Questions