Curve fitting using Python

Aim – To write a program to fit a linear and cubic polynomial for C_p.

Objective –

1. To write a code to fit a linear and cubic polynomial for the given C_p data.
2. Plot the linear and cubic fit curves.
3. Explain the parameters used to measure the fitness characteristics

Theory –

Curve fitting is the process of constructing a curve, or mathematical functions, which possess the closest proximity to the real series of data. By curve fitting, we can mathematically construct the functional relationship between the observed dataset and parameter values, etc.

After fitting data with one or more models, we should evaluate the goodness of fit. We can evaluate the goodness of fit in two ways – Graphical and Numerical.

Generally speaking, graphical measures are more beneficial than numerical measures because they allow you to view the entire data set at once, and they can easily display a wide range of relationships between the model and the data. In practice, depending on your data and analysis requirements, you might need to use both types to determine the best fit.

After using graphical methods to evaluate the goodness of fit, you should examine the goodness of fit statistics. The Curve Fitting Toolbox supports these goodness of fit statistics for parametric models:

• The sum of squares due to error (SSE)
• R-square
• Adjusted R-square
• Root mean squared error (RMSE)

For the current fit, these statistics are displayed in the Results pane in the Curve Fitting app. For all fits in the current curve-fitting session, you can compare the goodness-of-fit statistics in the Table of fits.

Assume data points (x1,y1) , (x2,y2) , (x3,y3) , ... (xn,yn) and assume the best fit for these set of points will be y = f(x), the criterion under which the equation of the curve will be the best fit possible is as follows,

Suppose each of these n points has an individual error (difference in y co-ordinate values)

The error for a general "i" point can be given by the function,

Error(i) = |y(i) − value of f(x) at x(i)| = |y(i) − f(x(i))|

then,

SSE (Sum of squared errors) - This statistic measures the total deviation of the response values from the fit to the response values.

`SSE = ∑_(i=1)^n (Err0r (i))^2`

R-square - It measures how successful the fit is in explaining the variation of the data. The higher the value of R-square the better the fit. R-square ranges from a value of 0 to 1.

The value of R-square can be calculated using the following formula,

`R^2= (SSR)/(SST)`

where,

SSR = sum of squares of the regression, which is basically the square of the difference between the value of the fitting function at a given point and the arithmetic mean of the data set.

`SSR= ∑_(i=1)^n(f(x(i))-Mean)^2`

and, `SST=SSR+SSE`

Adjusted R-square –

Adjusted R-square is a normalizing variation of R-square which neglects the variation of R-square due to insignificant independent variables.

`R_(adj)^2 = 1-[ (1- R^2 )(n-1)/(n-k-1) ]`

where,

n - the number of data points

k - the number of independent variables

RMSE (Root Mean Squared Error) –

This statistic is also known as the fit standard error and the standard error of the regression.

`RMSE= √(SSE)/n`

Using these four quantities, we can effectively deduce the goodness of a fit. Typically R-square is greater than 0.95 then the fit is considered good.

Q 1. What does ‘popt’ and ‘pcov’ mean?

Answer – ‘popt’ – ‘popt’ is a one dimensional array of the best estimates for the parameter values, each entry matches the order in the function definition.

‘pcov’ – The estimated covariance of ‘popt’. ‘pcov’ is the covariance matrix showing the uncertainty and interdependence of each parameter in ‘popt’.

Q2. What does np.array(temperature) do?

Answer – It creates a list of temperature values as an array so that we can use it in our program.

Q3. What does the * in *popt mean?

Answer – It refers to the coefficient values stored in ‘popt’. *in *popt uses the values of coefficients which were stored in ‘popt’ and then use them for making the plot.

Q4. What needs to be done in order to make the curve fit perfect?

Answer – To get the best fit of the curve we should go first for the geometric fitting and then after analyzing the data we should opt for the numerical values of the polynomial equations. In order to make best curve fit, we can either increase the degree of polynomial or we can use the split data-set method.

Code –

import math
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

# Linear curve fit function
def func_l(t1, a1, b1):
	return a1*t1 + b1

# Cubic curve fit function
def func_c(t2, a2, b2, c2, d2):
	return a2*pow(t2,3) + b2*pow(t2,2) + c2*t2 + d2

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

temperature, cp = read_file()
popt, pcov = curve_fit(func_l,temperature,cp)
fit_cpl = func_l(np.array(temperature), *popt)

popt, pcov = curve_fit(func_c,temperature,cp)
fit_cpc = func_c(np.array(temperature), *popt)

# Fitness parameters - Linear Curvefit
SSE_linear = [];
SSR_linear = [];
for i in range(len(cp)):
	SSE_l = 0
	SSE_l = pow(cp[i]- fit_cpl[i],2)

	SSR_l = 0
	SSR_l = pow(fit_cpl[i]- (sum(cp)/len(cp)),2)

SSE_linear.append(SSE_l)
SSR_linear.append(SSR_l)
SST_linear = np.array(SSE_linear) + np.array(SSR_linear)

R_square_linear = (np.array(SSR_linear)/np.array(SST_linear))
print("R_square_linear =",R_square_linear)

# Fitness parameters - Linear Curvefit
SSE_cubic = [];
SSR_cubic = [];
for i in range(len(cp)):
	SSE_c = 0
	SSE_c = pow(cp[i]- fit_cpc[i],2)

	SSR_c = 0
	SSR_c = pow(fit_cpc[i]- (sum(cp)/len(cp)),2)

SSE_cubic.append(SSE_c)
SSR_cubic.append(SSR_c)
SST_cubic = np.array(SSE_cubic) + np.array(SSR_cubic)

R_square_cubic = (np.array(SSR_cubic)/np.array(SST_cubic))
print("R_square_cubic =",R_square_cubic)


#Plotting
plt.plot(temperature,cp,'r-',linewidth = 5)
plt.plot(temperature,fit_cpl,'g-',linewidth = 3)
plt.plot(temperature,fit_cpc,'b-',linewidth = 3)
plt.legend(['Original dataset','Linear curvefit','Cubic curvefit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp (KJ/kmol-K)')
plt.show()

Explanation –

First of all, we import the modules which are going to be used in the program. Then we define the linear and cubic curve-fir functions with names as ‘func_l’ and ‘func_c’ respectively.

Then we define another function with name ‘read_file’ ehich gives us the list of cp and temperature values as array.

After that, we use inbuilt ‘curve_fit’ command to generate the coefficients of the polynomials with the use of ‘popt’ and ‘pcov’ commands and use it to create out fitness curve.

To check the degree of fitness, we use the fitness parameter ‘R-square’ which can be obtained by calculating SSE, SSR and SST respectively.

After that we use the ‘plot’ command, with legend, labels to create the desired curve_fit plot.

Result –

1. Linear and Cubic curve-fit plot –

2. R-square values –

For linear fit = 0.9591

For cubic fit = 0.9907

Discussion –

From the results, we can say that higher the degree of polynomial, better it fits the given curve.

R-square values for both linear and cubic fit are greater than 0.95, hence our fit is good.

Reference –

http://phy224.ca/19-curvefit/index.html#:~:text=popt%20is%20a%20one%20dimensional,of%20each%20parameter%20in%20popt

https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html?highlight=popt