Fitting 3d data

Question

I would like to fit a function to a 3d data. I read the data with pandas:

df = pd.read_csv('data.csv')
Ca = df.Ca
q = df.q
L = df.L0

Then, I define my 3d function (z=f(x,y)) as:

def func(q, Ca, l0, v0, beta):
    return l0 + q*v0*(1+beta/(q*Ca))

then I use curve_fit to find the best fit parameters:

from scipy.optimize import curve_fit
guess = (1,1,1)
popt, pcov = curve_fit(func, q,Ca,L, guess)

And it gives me the following errors:

ValueError: `sigma` has incorrect shape.

Do you know what is the mistake and how to solve it? Thanks a lot for your help

Answer 1

I'm not sure in the correctness of your Model-func - if proposed objective function can converge at all - & am sure in uncorrectness of the way you're passing the params. & having empty plot from James Phillips's code, just plot his Model-func to fit some data

from scipy.optimize import curve_fit
import numpy as np

x = np.random.random(100)
y = np.random.random(100)
z = x**2 + y**2 + np.random.normal(0, 0.1, size=100)

def func(xy, a, b):
    x, y = xy
    return a * x * y  + b

# Perform curve fitting
popt, pcov = curve_fit(func, (x, y), z)

# As a quick check of whether the model may be overparameterized (resulting depedencies among dims, or Multicollinearity), calculate the condition number of the covariance matrix:
print(np.linalg.cond(pcov))

# Print optimized parameters
print(popt)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Create 3D plot of the data points and the fitted curve/surface
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue')

# make grid
x_range = np.linspace(x.min(), x.max(), 50)
y_range = np.linspace(y.min(), y.max(), 50)
X, Y = np.meshgrid(x_range, y_range)
Z = func((X, Y), *popt);
# plot fitted surface
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue')
ax.plot_surface(X, Y, Z, linewidth=2,  alpha = 0.3)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

P.S. but always pay attention to the Condition Number of your pcov from fit - np.linalg.cond(pcov) - see in docs

large value is cause for concern. The diagonal elements of the covariance matrix, which is related to uncertainty of the fit, gives more information: np.diag(pcov) array([1.48814742e+29, 3.78596560e-02, 5.39253738e-03, 2.76417220e+28]) # may vary - in e.g. Note that the first and last terms are much larger than the other elements, suggesting that the optimal values of these parameters are ambiguous and that only one of these parameters is needed in the model.

OR use/change scaling

P.P.S. ref. curve_fit

Answer 2

Here is a graphical 3D fitter with 3D scatter plot, 3D surface plot, and 3D contour plot.

import numpy, scipy, scipy.optimize
import matplotlib
from mpl_toolkits.mplot3d import  Axes3D
from matplotlib import cm # to colormap 3D surfaces from blue to red
import matplotlib.pyplot as plt

graphWidth = 800 # units are pixels
graphHeight = 600 # units are pixels

# 3D contour plot lines
numberOfContourLines = 16


def SurfacePlot(func, data, fittedParameters):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

    matplotlib.pyplot.grid(True)
    axes = Axes3D(f)

    x_data = data[0]
    y_data = data[1]
    z_data = data[2]

    xModel = numpy.linspace(min(x_data), max(x_data), 20)
    yModel = numpy.linspace(min(y_data), max(y_data), 20)
    X, Y = numpy.meshgrid(xModel, yModel)

    Z = func(numpy.array([X, Y]), *fittedParameters)

    axes.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=1, antialiased=True)

    axes.scatter(x_data, y_data, z_data) # show data along with plotted surface

    axes.set_title('Surface Plot (click-drag with mouse)') # add a title for surface plot
    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label
    axes.set_zlabel('Z Data') # Z axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems


def ContourPlot(func, data, fittedParameters):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    x_data = data[0]
    y_data = data[1]
    z_data = data[2]

    xModel = numpy.linspace(min(x_data), max(x_data), 20)
    yModel = numpy.linspace(min(y_data), max(y_data), 20)
    X, Y = numpy.meshgrid(xModel, yModel)

    Z = func(numpy.array([X, Y]), *fittedParameters)

    axes.plot(x_data, y_data, 'o')

    axes.set_title('Contour Plot') # add a title for contour plot
    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    CS = matplotlib.pyplot.contour(X, Y, Z, numberOfContourLines, colors='k')
    matplotlib.pyplot.clabel(CS, inline=1, fontsize=10) # labels for contours

    plt.show()
    plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems


def ScatterPlot(data):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)

    matplotlib.pyplot.grid(True)
    axes = Axes3D(f)
    x_data = data[0]
    y_data = data[1]
    z_data = data[2]

    axes.scatter(x_data, y_data, z_data)

    axes.set_title('Scatter Plot (click-drag with mouse)')
    axes.set_xlabel('X Data')
    axes.set_ylabel('Y Data')
    axes.set_zlabel('Z Data')

    plt.show()
    plt.close('all') # clean up after using pyplot or else thaere can be memory and process problems


def func(data, a, alpha, beta):
    x = data[0]
    y = data[1]
    return a * (x**alpha) * (y**beta)


if __name__ == "__main__":
    xData = numpy.array([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0])
    yData = numpy.array([11.0, 12.1, 13.0, 14.1, 15.0, 16.1, 17.0, 18.1, 90.0])
    zData = numpy.array([1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.0, 9.9])

    data = [xData, yData, zData]

    initialParameters = [1.0, 1.0, 1.0] # these are the same as scipy default values in this example

    # here a non-linear surface fit is made with scipy's curve_fit()
    fittedParameters, pcov = scipy.optimize.curve_fit(func, [xData, yData], zData, p0 = initialParameters)

    ScatterPlot(data)
    SurfacePlot(func, data, fittedParameters)
    ContourPlot(func, data, fittedParameters)

    print('fitted prameters', fittedParameters)

    modelPredictions = func(data, *fittedParameters) 

    absError = modelPredictions - zData

    SE = numpy.square(absError) # squared errors
    MSE = numpy.mean(SE) # mean squared errors
    RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
    Rsquared = 1.0 - (numpy.var(absError) / numpy.var(zData))
    print('RMSE:', RMSE)
    print('R-squared:', Rsquared)