Curve Fitting and Quality of Fits using Python
Introduction
Knowledge in terms of data is Power. In engineering, plethora of data is generated. However, data doesnt make sense until it could be assessed and used for the predcition of the phenomena. Curve fitting is one of the most basic operation in this approach towards utilising and making sense out of the data. The data is fitted into a form of equation, so that equation could be used for further predictions of the phenomena.
However, it is a double edged sword, underfitting and overfitting is detrimental to the prediction and hence there are several parameters to assess the quality of the data fit. Among them, R square and RMS error are the two widely used benchmarking parameter for the quality of the fit
R square varies between 0 and 1. Closer to 0 indicates a poor fit and 1 indicates the best fit. However, R square alone sometimes might not indicate the quality of spread. Hence RMS error is used which stands for root mean square error. Closer to Zero would indicate predicted curve overlays on the initial data set and there is no spread. As the value of RMS error increases, spread increases. Value closer to 1 and less than 1 would be considered a best fitAnother metrics used is adjusted R square. This is used to map the effect of adding new variables to the data set.
Sometimes R square, RMS error can lie outside the tolerance band. In this case, quality of the fit has to be improved. This could mainly be done by two means
1) Either by breaking the data into different sets and having the individual equations for each set
2) Increasing the power of fit until the error comes within tolerance.
Both the methods are demonstrated in this code
Python code
The code is written to perform the curve fitting process and analyse the quality of fit. The Temeratur vs Cp data is loaded. This data is extracted and loaded into the array. Function is created to read and import the necessary data.
Multiple functions are created to define the equation of fit. Linear, Quadratic, Cubic upto order of 6. Functions are also created to calculate the R square and Root mean sqaure error for the given fit. In the main program, Cp and temperature is passed to Curve_fit which is imported from scipy module. Based on coefficients defined, predicted Cp curve is also plotted. This data is sent to calculate R square and RMS error for the given order of fit.
In method 2, temperature range is broken into two parts Anything less than 1000K is considered as the low temperature range and separate coefficients are developed to fit and anything greater than 1000K is considered as the high temperature range and another set of coefficients are developed.
Finally, the RMS error, R square is plotted from all the case studies and best method and coefficients are downselected
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 13 15:58:04 2020
@author: sweekar
This program is coded by Sweekar and is his intellectual property, contact in case of help
[email protected]
This program loads the data abd fits the data while varying the power of the equation
of fit to reduce the RMS error
Curve fitting is an important aspect of engineering as we generate tons of data
These data only makes sense once they are able to predict future trend
To do this, Curve fitting is used to generate the equation and variable is used to predict
Two ways of achieving this feat. Either breaking curves into two datasets or increasing the exponent
Drawback of this method is that it is susceptible to overfit
Advantage: A simple code to quickly find the trend and probable power
"""
#close all
#clear all
#clc
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import math
def lin_func(t,a,b):
return a*t + b
def quad_func(t,a,b,c):
return a*pow(t,2) + b*t + c
def cubic_func(t,a,b,c,d):
return a*pow(t,3) + b*pow(t,2) + c*t +d
def order4_func(t,a,b,c,d,e):
return a*pow(t,4) + b*pow(t,3) + c*pow(t,2) + d*t+ e
def order5_func(t,a,b,c,d,e,f):
return a*pow(t,5) + b*pow(t,4) + c*pow(t,3) + d*pow(t,2) + e*t +f
def order6_func(t,a,b,c,d,e,f,g):
return a*pow(t,6) + b*pow(t,5) + c*pow(t,4) +d*pow(t,3)+e*pow(t,2)+f*t+g
def read_file():
temperature = []
cp = []
for line in open('data','r'):
values = line.split(',')
temperature.append(float(values[0]))
cp.append(float(values[1]))
return [temperature,cp]
def R_square(temperature,cp,predicted_cp):
Mean = np.mean(cp)
for i in range(0,len(temperature)):
SSR = np.sum(pow((predicted_cp-Mean),2))
SSE = np.sum(pow((cp[i]-predicted_cp[i]),2))
SST = SSR + SSE
R_sq = SSR/SST
return R_sq
def RMSE(pcov):
for i in range(0,len(pcov)):
SSE = np.sum(pcov[i])
RMSE = math.sqrt(SSE/len(pcov))
return RMSE
temperature, cp = read_file()
#Method 1
#To improve curve fit by increasing the order of polynomial fit
#Linear Fit
popt1,pcov1 = curve_fit(lin_func, temperature,cp)
predicted_cp1 = lin_func(np.array(temperature), *popt1)
R_sq1 = R_square(temperature,cp,predicted_cp1)
RMSE1 = RMSE(pcov1)
plt.figure(1)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp1,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('Linear fit')
plt.savefig('1st_order_fit')
#Quad Fit
popt2,pcov2 = curve_fit(quad_func, temperature,cp)
predicted_cp2 = quad_func(np.array(temperature), *popt2)
R_sq2 = R_square(temperature,cp,predicted_cp2)
RMSE2 = RMSE(pcov2)
plt.figure(2)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp2,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('Quad fit')
plt.savefig('2nd_order_fit')
#Cubic Fit
popt3,pcov3 = curve_fit(cubic_func, temperature,cp)
predicted_cp3 = cubic_func(np.array(temperature), *popt3)
R_sq3 = R_square(temperature,cp,predicted_cp3)
RMSE3 = RMSE(pcov3)
plt.figure(3)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp3,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('Cubic fit')
plt.savefig('3rd_order_fit')
#4th order Fit
popt4,pcov4 = curve_fit(order4_func, temperature,cp)
predicted_cp4 = order4_func(np.array(temperature), *popt4)
R_sq4 = R_square(temperature,cp,predicted_cp4)
RMSE4 = RMSE(pcov4)
plt.figure(4)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp4,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('4th order fit')
plt.savefig('4th_order_fit')
#5th Order Fit
popt5,pcov5 = curve_fit(order5_func, temperature,cp)
predicted_cp5 = order5_func(np.array(temperature), *popt5)
R_sq5 = R_square(temperature,cp,predicted_cp5)
RMSE5 = RMSE(pcov5)
plt.figure(5)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp5,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('5th order fit')
plt.savefig('5th_order_fit')
#6th Order Fit
popt6,pcov6 = curve_fit(order6_func, temperature,cp)
predicted_cp6 = order6_func(np.array(temperature), *popt6)
R_sq6 = R_square(temperature,cp,predicted_cp6)
RMSE6 = RMSE(pcov6)
plt.figure(6)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature,predicted_cp6,color='red')
plt.legend(['Actual data','Curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('6th order fit')
plt.savefig('6th_order_fit')
#Method 2
#To Divide the zone into multiple zones and solve
#1000K is taken as the mid temp, anything less than 1000 is considered as the Low temperature range
#Anything above 1000K is taken as high temp range with different polynomial fit
temperature1 = temperature[0:700]
cp1 = cp[0:700]
popt7,pcov7 = curve_fit(order4_func, temperature1,cp1)
predicted_cp = order4_func(np.array(temperature1), *popt7)
R_sq7 = R_square(temperature1,cp1,predicted_cp)
RMSE7 = RMSE(pcov7)
temperature2 = temperature[701:]
cp2 = cp[701:]
popt8,pcov8 = curve_fit(order4_func, temperature2,cp2)
predicted_cp0 = order4_func(np.array(temperature2), *popt8)
R_sq8 = R_square(temperature2,cp2,predicted_cp0)
RMSE8 = RMSE(pcov8)
R_sq = (R_sq7+R_sq8)/2
RMSE = (RMSE7+RMSE8)/2
plt.figure(7)
plt.plot(temperature,cp,color = 'blue')
plt.plot(temperature1,predicted_cp,color='red')
plt.plot(temperature2,predicted_cp0,color='yellow')
plt.legend(['Actual data','Low temp Curve fit','High temp curve fit'])
plt.xlabel('Temperature[K]')
plt.ylabel('cp')
plt.title('Split fit')
plt.savefig('CurvefitsplitMethod2')
#Plotting Data quality metrics
x = [1,2,3,4,5,6,7]
R_sqt = [R_sq1,R_sq2,R_sq3,R_sq4,R_sq5,R_sq6,R_sq]
RMSET = [RMSE1,RMSE2,RMSE3,RMSE4,RMSE5,RMSE6,RMSE]
plt.figure(8)
plt.plot(x,R_sqt,color='brown')
plt.xlabel('Cases')
plt.ylabel('R Squared')
plt.title('R Square')
plt.savefig('R_Square_dependency')
plt.figure(9)
plt.plot(x,RMSET,color='yellow')
plt.xlabel('Cases')
plt.ylabel('Root mean square error')
plt.title('RMS Error')
plt.savefig('RMSE_dependency')Results and Plots
Method 1:
Case 1: Linear Fit
Case 2: Quadratic fit
Case 3: Cubic Fit
Case 4: 4th Order fit
Case 5: 5th order fit
Case 6: 6th Order fit
Method 2: Temperature is broken into two parts and are fitted for 4th order polynomial separately in zone 1 and zone 2
Quality of the Fits
R square Plot
RMS Error Plot
Discussion and Report
1) R square is consistently greater than 0.9 for all the equation of fits. Hence R square alone doesnt represent the quality of fit
2) RMS error is close to 0.7 when Linear equation is used. This consistently reduces till order of 4 and starts increasing again
3) From Method 1, best fit could be obtained at order of 4 with data spread of 0.2
4) From Method 2, where 4th order polynomial is separately used in two regions, RMS error is equal to 1e-13 which is zero for practical considerations
5) Hence Method 2 with 4th order polynomial gives best fit corroborated by r square and RMS error values