curve fitting and matplotlib

Question

I have these these three variables and I want to plot the relationship between two of them two variables:

x = [0.125735  , 0.11753342, 0.11572967, 0.11963533, 0.1255283 ,
       0.13183589, 0.13904629, 0.14754317, 0.15548172, 0.16429631,
       0.17474308, 0.18641375]
y = [0.11917991, 0.10663986, 0.09897077, 0.09291739, 0.08743263,
       0.08346636, 0.08161819, 0.08132199, 0.08216186, 0.0834759 ,
       0.08551088, 0.08770163]
z = [1, 2, 3, 4, 5,
       6, 7, 8, 9, 10 ,
       11, 12]

the image shows x,y. I want to fit a line that goes through all of these points and label each point with z.

Answer 1

Here is a graphical 3D surface fitter with 3D scatterplot, 3D surface plot, and contour plot. The contour plot shows that the surface is clearly curved, which is why this example's equation fit better than a flat surface plane. You should be able to click-drag on the 3D plots and rotate them in 3-space for visual inspection. Using your data and a simple power equation "z = a * pow(x, b) + c * pow(y, d)" with fitted parameters a = 2.14091547e+02, b = 1.56841786e+00, c = -2.24366942e+03, and d = 2.69437535e+00 yields RMSE = 0.1122 and R-sqiared = 0.9989

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


xData = numpy.array([0.125735  , 0.11753342, 0.11572967, 0.11963533, 0.1255283 , 0.13183589, 0.13904629, 0.14754317, 0.15548172, 0.16429631, 0.17474308, 0.18641375], dtype=float)
yData = numpy.array([0.11917991, 0.10663986, 0.09897077, 0.09291739, 0.08743263, 0.08346636, 0.08161819, 0.08132199, 0.08216186, 0.0834759 , 0.08551088, 0.08770163], dtype=float)
zData = numpy.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 , 11, 12], dtype=float)


def func(data, a, b, c, d):
    x = data[0]
    y = data[1]
    return a * numpy.power(x, b) + c * numpy.power(y, d)


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 there 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 there 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 there can be memory and process problems


if __name__ == "__main__":

    data = [xData, yData, zData]

    initialParameters = [100.0, 1.0, 1000.0, 1.0]

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

Answer 2

No Centering

Here's a quick example of fitting the points in polar coordinates:

x = [0.125735  , 0.11753342, 0.11572967, 0.11963533, 0.1255283 ,
     0.13183589, 0.13904629, 0.14754317, 0.15548172, 0.16429631,
     0.17474308, 0.18641375]
y = [0.11917991, 0.10663986, 0.09897077, 0.09291739, 0.08743263,
     0.08346636, 0.08161819, 0.08132199, 0.08216186, 0.0834759 ,
     0.08551088, 0.08770163]
z = [1, 2, 3, 4, 5,
     6, 7, 8, 9, 10 ,
     11, 12]

# you need numpy
import numpy as np
import matplotlib.pyplot as plt

x = np.array(x)
y = np.array(y)

r = np.sqrt(x*x + y*y)
theta = np.arctan2(x, y)

plt.scatter(theta, r, z, z)
plt.colorbar()
plt.grid()
plt.title('polar coords')
plt.xlabel('$\\theta$ [rad]')
plt.ylabel('r')

p = np.polyfit(theta, r, 2)
xfit = np.linspace(0.8, 1.15, 15)
yfit = np.polyval(p, x)

plt.plot(xfit, yfit, '--')
plt.legend(['original data', 'fit'])

With Centering

We might do better if we center the points first:

# find the averages to find the centroid of the data
x_avg = x.mean()
y_avg = y.mean()

# center the data
x_star = x - x_avg
y_star = y - y_avg

# now find the radii
r = np.sqrt(x_star*x_star + y_star*y_star)

# make sure points are between 0 adn 360[degrees]
theta = np.degrees(np.arctan2(y_star, x_star)) % 360

plt.scatter(theta, r, z, z)
plt.colorbar()
plt.grid()
plt.title('polar coords')
plt.ylabel('r')
plt.xlabel('$\\theta$ [degrees]')

# fit with 3rd order polynomial
# because there are 2 inflection points
p = np.polyfit(theta, r, 3)

# plot fit
x_fit = np.linspace(90, 360, 270)
y_fit = np.polyval(p, x_fit)
plt.plot(x_fit, y_fit, '--')
plt.legend(['original data', 'fit'])

Final Output Fit Line

Here is the final output fit line through the original data:

x_out = y_fit * np.cos(np.radians(x_fit)) + x_avg
y_out = y_fit * np.sin(np.radians(x_fit)) + y_avg

plt.scatter(x, y, z, z)
plt.plot(x_out, y_out, '--')
plt.colorbar()
plt.grid()
plt.title('output fit line')
plt.legend(['original data', 'fit'])