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.

But it's not a function, ie it doesn't have a unique x for each y, so you can't fit a polynomial function to it. Maybe you can you translate this to polar coordinates? — Mark Mikofski
– Mark Mikofski, Commented Sep 20, 2019 at 18:22
@MarkMikofski sorry I noticed that the question was lacking further explanation. I updated the question. — amaatouq
– amaatouq, Commented Sep 20, 2019 at 18:29
Thanks for the explanation, so IIUC, point 1 (0.126, 12) is at the top of the plot, point 2 (0.12, 0.11) is the next going counter clockwise around to the last point 12 (0.186, 0.09) at the bottom right. So this means that a line passing through these points curves so that at a given x, say 0.12, there would be two values of y, picture a line vertically from 0.12 crossing the fit line twice. That is not a function. — Mark Mikofski
– Mark Mikofski, Commented Sep 20, 2019 at 18:35
Since the range of x and y are about the same, I would recommend switching to a polar system of coordinates. (1) find the average x and average y, and shift all of the points by those averages so the center is now at (0, 0) then calculate the radii to each point, r=np.sqrt(x*x + y*y), and the angle to each point theta=np.arctan2(x, y) then it will be a function you can fit with a polynomial — Mark Mikofski
– Mark Mikofski, Commented Sep 20, 2019 at 18:38

Mark Mikofski · Accepted Answer · 2019-09-20 19:26:58Z

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'])

James Phillips · Accepted Answer · 2019-09-20 22:02:56Z

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)

Collectives™ on Stack Overflow

curve fitting and matplotlib

2 Answers 2

No Centering

With Centering

Final Output Fit Line

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

No Centering

With Centering

Final Output Fit Line

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related