Higher Order Polynomial vs Piecewise Curve Fitting using PYTHON
Objective
The aim of the project is to write a PYTHON code to perform curve fitting for the provided data points from the thermodynamic data file, i.e. ‘data’.
- Curve Fitting with higher order polynomials
- Curve Fitting with two defined functions using if and conditional loops
Also, get some understanding about curve fitting, measures for ‘Goodness-of-fit’ and regression model.
Curve Fitting
Curve fitting is a process of constructing a curve or mathematical function that best fits to series of data points. The mathematical function established from a given set of data points provides relationship between dependent and independent variables. Fitted curves can be used as an aid for data visualization and also to predict the data beyond the range using extrapolation, while being subjected to degree of uncertainty. Further, regression analysis can be performed on fitted data for the uncertainty presented in fitted model due to some random errors.
The data points can be fitted to fitting functions such as:
- Polynomial Function –
- Linear function (first order polynomial): `y=ax+b`
- Quadratic function (second order polynomial): `y=ax^2+bx+c`
- and so on …
- Trigonometric Functions – terms such as sine and cosine
- Exponential and Logarithmic Functions – terms such as `e^x`, ln(x), log(x) etc.
- Spectroscopy related Functions – such as Gaussian (normal distribution), Lorentzian (Cauchy distribution), Voigt (involving both Gaussian and Lorentzian), etc.
Measures for ‘Goodness-of-fit’
- Standard Error (S):
Also called as standard error of estimate (SEE) or root mean square error (RMSE). It provides absolute measure of average distance of data points from fitted curve.
The standard error has the unit of dependent variable and its value is preferred to be lower i.e. ‘0’ being perfect. It is valid for both linear and non-linear models.
- R-squared (R2):
Most well-known and widely used amongst the goodness-of-fit. R-squared provides the relative measure to higher percentage of the dependent variable variance with respect to fitted curve.
Higher R-squared values i.e. greater than 95% for fitted curve is considered to be good. On the contrary to standard error, R-squared is valid for only linear models. Yes, linear models can use polynomials to model curvature. Then, what does linear means? Explained in next section.
- Biased or Unbiased Model:
The regression model is said to be unbiased, if the mean of sampling distribution is equal to the mean of estimated values obtained through regression model, else regression model is biased.
Linear vs Non-Linear Regression Models
Based on the terms, its intuitive that linear model produces straight lines and non-linear curvatures. But it’s not true, rather there is ‘non’ such difference. Here, we will be discussing with respect to single regression model irrespective of being multi-variate regression model.
Linear Regression Models follows a particular form:
- Suppose, ‘y’ being dependent variable (predicted or response variable)
- Always contain constant, βo also called as intercept
- Weight or slope or parameter (β1, β2, β3 and so on…, also called as coefficients) multiplied by independent variable (suppose, ‘x’ also called as predictor variable)
It may be in polynomial form, reciprocal form, logarithmic form, etc.
Examples:
In simpler words, Non-Linear Regression Models are those models which are not linear models and follows a particular form, i.e. `y=f(x,β_i)`. Does it sound great? Nope, then best way to say about any model is to identify whether it is linear or not. Since, non-linear models are quite complex and their enormous variants exist.
One interesting thing about some of the non-linear regression model is that it can be transformed into linear regression model.
Examples:
CASE-1
Code Structure: Curve fit using higher order polynomials
PYTHON Code: Curve fit using higher order polynomials
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# curve fit function
def func1(t,a,b):
return a*t+b
def func2(t,a,b,c):
return a*pow(t,2)+b*t+c
def func3(t,a,b,c,d):
return a*pow(t,3)+b*pow(t,2)+c*t+d
def func4(t,a,b,c,d,e):
return a*pow(t,4)+b*pow(t,3)+c*pow(t,2)+d*t+e
def func5(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 func6(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 func7(t,a,b,c,d,e,f,g,h):
return a*pow(t,7)+b*pow(t,6)+c*pow(t,5)+d*pow(t,4)+e*pow(t,3)+f*pow(t,2)+g*t+h
def func8(t,a,b,c,d,e,f,g,h,i):
return a*pow(t,8)+b*pow(t,7)+c*pow(t,6)+d*pow(t,5)+e*pow(t,4)+f*pow(t,3)+g*pow(t,2)+h*t+i
def func8(t,a,b,c,d,e,f,g,h,i,j):
return a*pow(t,9)+b*pow(t,8)+c*pow(t,7)+d*pow(t,6)+e*pow(t,5)+f*pow(t,4)+g*pow(t,3)+h*pow(t,2)+i*t+j
# reading thermodynamic data file
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]
# main program
temperature,cp = read_file()
popt,pcov = curve_fit(func1,temperature,cp)
fit_cp1 = func1(np.array(temperature),*popt)
popt,pcov = curve_fit(func2,temperature,cp)
fit_cp2 = func2(np.array(temperature),*popt)
popt,pcov = curve_fit(func3,temperature,cp)
fit_cp3 = func3(np.array(temperature),*popt)
popt,pcov = curve_fit(func4,temperature,cp)
fit_cp4 = func4(np.array(temperature),*popt)
popt,pcov = curve_fit(func5,temperature,cp)
fit_cp5 = func5(np.array(temperature),*popt)
popt,pcov = curve_fit(func6,temperature,cp)
fit_cp6 = func6(np.array(temperature),*popt)
popt,pcov = curve_fit(func7,temperature,cp)
fit_cp7 = func7(np.array(temperature),*popt)
popt,pcov = curve_fit(func8,temperature,cp)
fit_cp8 = func8(np.array(temperature),*popt)
# nth order polynomial names
poly_name = [\'Linear Curve Fit\',\'Quadratic Curve Fit\',\'Cubic Curve Fit\',\'Quartic Curve Fit\',
\'Quintic Curve Fit\',\'Hexic Curve Fit\',\'Septic Curve Fit\',\'Octic Curve Fit\']
mean_sample = []
mean_model = []
error = []
sse = []
ssr = []
sst = []
r2 = []
rmse = []
for m in range(8):
# Goodness-of-fit
mean_sample.append(sum(cp)/len(cp))
mean_model.append(sum(eval(\'fit_cp\'+str(m+1)))/len(eval(\'fit_cp\'+str(m+1))))
error.append(abs(cp-eval(\'fit_cp\'+str(m+1))))
sse.append(sum(pow(cp-eval(\'fit_cp\'+str(m+1)),2)))
ssr.append(sum(pow(eval(\'fit_cp\'+str(m+1))-np.mean(cp),2)))
sst.append(sse[m]+ssr[m])
r2.append(ssr[m]/sst[m])
rmse.append(math.sqrt(sse[m]/len(cp)))
# Plot actual data set and fitted polynomial
plt.figure()
plt.plot(temperature,cp,color=\"blue\")
plt.plot(temperature,eval(\'fit_cp\'+str(m+1)),color=\"red\")
plt.xlabel(\'Temperature [K]\')
plt.ylabel(\'Cp\')
plt.legend([\'Actual Data\', poly_name[m]])
plt.xlim([200,3650])
plt.ylim([900,1350])
if round(r2[m],4) < 0.95:
plt.title(\'R-Square is \'+str(round(r2[m],4))+\', implies curve fit is not good\')
elif round(r2[m],4) == 1:
plt.title(\'R-Square is \'+str(round(r2[m],4))+\', implies curve fit is perfect\')
else:
plt.title(\'R-Square is \'+str(round(r2[m],4))+\', implies curve fit is good\')
# plt.savefig(str(m+1)+\'.png\')
print(\'Mean of sample data set:\\t\'+ str(round(mean_sample[m],4)))
print(\'Mean of fitted data:\\t\\t\'+ str(round(mean_model[m],4)))
print(\'Standard Error, S:\\t\\t\\t\'+ str(round(rmse[m],4)))
print(\'R-squared:\\t\\t\\t\\t\\t\'+ str(round(r2[m],4)))
plt.show()Discussion: Curve fit using higher order polynomials
Number of polynomial functions were considered for curve fitting in order to achieve best fit i.e. first order (linear) to eighth order (octic) polynomial. In order to get the essence of polynomial, highest order polynomial considered here is represented below:
`c_p=at^8+bt^7+ct^6+dt^5+et^4+ft^3+g t^2+ht+i, t∈[300,3500]`
Here,
a, b, c, …, and i are unknown parameters or coefficients evaluated after curve fitting
cp is estimated response
t is the independent variable or estimator
After defining above functions, the curve fitting is performed using ‘curve_fit’ tool present in the optimization and root finding module of SciPy library represented by ‘scipy.optimize’. The syntax for curve fitting is given below:
popt,pcov = curve_fit(function,x,y)
The returns from above syntax are popt and pcov. The ‘popt’ array gives optimal values for the unknown parameters or coefficients so that the squared residuals of f(xdata,*popt) i.e. ydata is minimized. The ‘pcov’ matrix represents the estimated covariance of popt.
So based on above step, the evaluated functions that fits the given data have been plotted and their coefficients have also been tabulated below.
The evaluated functions goodness-of-fit verified by computing standard error, R-squared and biasness of estimated data whose outcomes are tabulated below.
So, from the above tabulated data, it can be observed that the higher order polynomial best fits the given data i.e. here, in this case 8th order polynomial best fits. But coefficients related to higher order polynomial are almost zero i.e. of order e-10 or higher. Thus, it doesn’t make sense to fit given data with higher order polynomial.
Instead, we can split the given data into pieces and fit them with relevant function such as linear, quadratic, cubic, exponential, etc.
Case-2
Code Structure: Piecewise fit using if and conditional loop
PYTHON Code: Piecewise fit using if and conditional loop
import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# curve fit function
def func1(t,a,b,c):
return a*pow(t,2)+b*t+c
def func2(t,m,n,o):
return m*pow(t,2)+n*t+o
# reading thermodynamic data file
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]
temp1 = []
temp2 = []
cp1 = []
cp2 = []
# main program
temperature,cp = read_file()
for i in range(0,len(temperature)):
if temperature[i] < 1323:
temp1.append(temperature[i])
cp1.append(cp[i])
else:
temp2.append(temperature[i])
cp2.append(cp[i])
popt1,pcov1 = curve_fit(func1,temp1,cp1)
fit_cp1 = func1(np.array(temp1),*popt1)
popt2,pcov2 = curve_fit(func2,temp2,cp2)
fit_cp2 = func2(np.array(temp2),*popt2)
fit_cp = []
for i in range(0,2):
if i == 0:
for j in range(0,len(temp1)):
fit_cp.append(fit_cp1[j])
else:
for j in range(0,len(temp2)):
fit_cp.append(fit_cp2[j])
# Goodness-of-fit
mean_sample = sum(cp)/len(cp)
mean_model = sum(fit_cp)/len(fit_cp)
error = abs(np.array(cp)-np.array(fit_cp))
sse = sum(pow(np.array(cp)-np.array(fit_cp),2))
ssr = sum(pow(fit_cp-np.mean(cp),2))
sst = sse+ssr
r2 = ssr/sst
rmse = math.sqrt(sse/len(cp))
# Plot actual data set and fitted curve
plt.figure()
plt.plot(temperature,cp,color=\"blue\",linewidth=1)
plt.plot(temperature,fit_cp,color=\"red\",linewidth=1)
plt.legend([\'Actual Dataset\',\'Fitted Curve\'])
plt.xlabel(\'Temperature, K\')
plt.ylabel(\'Cp\')
# Plot title based on R-squared goodness-of-
if round(r2,4) < 0.95:
plt.title(\'R-Square is \'+str(round(r2,4))+\', implies curve fit is not good\')
elif round(r2,4) == 1:
plt.title(\'R-Square is \'+str(round(r2,4))+\', implies curve fit is perfect\')
else:
plt.title(\'R-Square is \'+str(round(r2,4))+\', implies curve fit is good\')
print(\'Mean of sample data set:\\t\'+ str(round(mean_sample,4)))
print(\'Mean of fitted data:\\t\\t\'+ str(round(mean_model,4)))
print(\'Standard Error, S:\\t\\t\\t\'+ str(round(rmse,4)))
print(\'R-squared:\\t\\t\\t\\t\\t\'+ str(round(r2,4)))
plt.savefig(\'Piecewise fit.png\')
plt.show()Discussion: Piecewise fit using if and conditional loop
The above figure shows the plot for specific heat capacity at constant pressure with respect to temperature for given dataset. As discussed previously, the higher order polynomials are not preferred irrespective of them being estimating perfectly. It is better to fit data set with multiple function by splitting the dataset at appropriate point.
The below figure shows the fitted curve using two polynomial functions. The baseline of 1323 K temperature was considered to split the data set of cp and temperature into two pieces based on visual inspection and iteration to achieve the best fit possible.
The two quadratic functions considered for curve fitting are given by:
`c_(p1)=at^2+bt+c,t∈[300,1323]`
`c_(p2)=mt^2+nt+o,t∈[1323,3500]`
Here,
a, b, c, m, n, and o are unknown parameters or coefficients evaluated after curve fitting
cp1 and cp2 are estimated response
t is the independent variable or estimator
After defining above functions, the curving fitting is performed as discussed earlier.
Thus, the evaluated functions that fits the given data are:
`c_(p1)=-(1.495×10^(-4))t^2+(4.953×10^(-1))t+(8.060×10^(+2)),t∈[300,1323]`
`c_(p2)=-(5.516×10^(-6))t^2+(9.750×10^(-2))t+(1.036×10^(+4)),t∈[1323,3500]`
The above functions goodness-of-fit verified by computing standard error, R-squared and biasness of estimated data whose outcomes are shown below in figure.
As the mean of sample data and fitted data is equal, it can be inferred that fitted data is unbiased. The standard error is 1.36 which states that estimated ydata may get deviated by actual ydata by the distance of ±1.36 which is considerable. Also, the R-squared error is almost equal to 1.0.
Overall, we can say that the functions evaluated are best and unbiased fit to the actual dataset.