Week 5 - Curve fitting

AIM:

To write a code to perform curve fitting in Python.

OBJECTIVES:

1. Write code to fit a linear and cubic polynomial for the Cp data provided.
2. Plot the linear and cubic fit curves along with the raw data points.
3. Explain the parameters used to measure the fitness characteristics for both the curves.

Theory -

Introduction of Curve Fitting Method :-

Curve fitting is a process of determining a best possible curve for a given data. This method is useful for extrapolation or interpolation of data. In simple words, curve fitting process can be used to find a data outside the given range.

There are various types of curve fitting process and the use of these process depends on the type of data given. Some of them are as follows -

1. Linear Curve Fitting (Fitting of first order polynomial to a data)

2. Cubic Curve Fitting (Fitting of third order polynomial to a data)

3. Higher Order Curve Fitting (Fitting of higher order polynomial to a data)

4. Curve Fitting using a Splitwise Polynomial Method (Using different/same order polynomial for the different divisions (parts) of given data range)

PROGRAM:

Programming language used it Python.

# Program to measure fitness characteristics for the linear and cubic polynomial for the Cp data 
import math 
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.optimize import curve_fit 
 
# Curve fitting functions
# Linear function 
def function_1(t,a,b):  
	return a*t + b  

# Cubic function  
def function_2(t,a,b,c,d):  
	return a*pow(t,3) + b*pow(t,2) + c*t +d   
 
# 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] 

 
# Storing temperature and cp in the form of lists
temperature , cp = read_file() 
 
popt, pcov = curve_fit(function_1,temperature,cp) 
fit_cp_l = function_1(np.array(temperature), *popt) 

# For Linear curve fit 
print('For Linear Fit') 
# Printing values of parameters and covarience matrix 
print('The value of parametrs a = ' + '%.3f' %(popt[0]) + ' & b = ' + '%.3f' %(popt[1]))  
print('The Co-varience matrix is:' + 'n' + str(pcov)) 

# Plotting LInera curve fit
plt.figure(1) 
plt.plot(temperature,cp,'b--') 
plt.plot(temperature,fit_cp_l,color='red',linewidth = 1) 
plt.legend(['Actual data','Linear Curve fit']) 
plt.xlabel('Temperature (K)') 
plt.ylabel('Cp') 
plt.title('Linear Curve fitting') 
plt.show() 


popt, pcov = curve_fit(function_2,temperature,cp) 
fit_cp_c = function_2(np.array(temperature), *popt) 

# For Cubic curve fit 
print('nFor Cubic Fit:') 
# Printing values of parameters and Co-varience 
print('The value of a = ' + '%.8f' %popt[0] + ' b = ' + '%.8f' %popt[1] + ' c = ' + '%.8f' %popt[2] + ' d = ' + '%.8f' %popt[3]) 
print('The Co-varience matrix is: n' + str(pcov))

# PLotting cubic curve fit
plt.figure(2) 
plt.plot(temperature,cp,'b--') 
plt.plot(temperature,fit_cp_c,color='red',linewidth = 1) 
plt.legend(['Actual data','Cubic Curve fit']) 
plt.xlabel('Temperature (K)') 
plt.ylabel('Cp') 
plt.title('Cubic Curve fitting') 
plt.show() 
 
 
''' 
 Measuring the fitness characetristics of the curves  
 
1. Sum of squares due to error (SSE)
SSR (sum of squares of the regression)  
SST = SSE + SSR  
2. R-square 
3. Adjusted R-square
4. Root mean squared error (RMSE) 

 Finding mean of Cp data  
 mean = ( sum of all elements)/(total number of elements)

 ''' 
 
s = np.sum(cp) 
l = np.size(cp) 
m = s/l 

# For Linear curve fit 
print('nFor Linear Fit') 

linear_ssr = 0 
linear_sse = 0 

for i in range(l):     
	linear_error = abs((np.sum((cp[i] - fit_cp_l[i]))))     
	linear_sse = linear_sse+pow(linear_error,2)  						# Sum of squared errors   
	linear_ssr = linear_ssr+ np.sum(pow((fit_cp_l[i] - m),2)) 			# Sum of squares of regression


print('SSE:',linear_sse) 

linear_sst = linear_sse + linear_ssr		# Total of sum of squares

# Calculating R-square
linear_Rsquare = linear_ssr/linear_sst 
print('Rsquare :' ,linear_Rsquare)

# Calculating Adjusted R-square  
k=1;        # No. of independent variables 
AdRsquare_L = 1- (((1-linear_Rsquare)*(l-1))/(l-k-1));    # For linear fit 
print('Adjusted R-square  :' ,AdRsquare_L)

# Calculating Root mean squared error (RMSE) 
linear_rmse = pow((linear_sse/l),0.5) 
print('RMSE :',linear_rmse) 
print('') 
 
 
# For Cubic curve fit 
print('For Cubic Fit') 

cubic_sse = 0 
cubic_ssr = 0 

for j in range(l):     
	cubic_error = abs(np.sum((cp[i] - fit_cp_c[i])))     
	cubic_sse = cubic_sse+(pow(cubic_error,2))    					# Sum of squared errors  
	cubic_ssr = cubic_ssr+np.sum(pow((fit_cp_c[i] - m),2)) 			# Sum of squares of regression
 

print('SSE :',cubic_sse) 
 
cubic_sst = cubic_sse + cubic_ssr 		# Total of sum of squares

# Calculating R-square
Cubic_Rsquare= cubic_ssr/cubic_sst 
print('Rsquare :',Cubic_Rsquare) 

# Calculating Adjusted R-square  
AdRsquare_C = 1- (((1-Cubic_Rsquare)*(l-1))/(l-k-1));     # Adjusted r square for cubic fit
print('Adjusted R-square :' ,AdRsquare_C)

# Calculating Root mean squared error (RMSE) 
cubic_rmse = pow((cubic_sse/l),0.5) 
print('RMSE:',cubic_rmse)

EXPLANATION OF THE PROGRAM:

1. Program starts with importing modules. Modules needed for this task are math, matplotlib, numpy and scipy to perform curve fitting operations.

2. After importing modules, next step is to define function "function_1”for linear polynomial and ‘function_2’ for cubic polynomial. The equation for linear polynomial is:

Cp = a*t + b

The equation for cubic polynomial is:

Cp = a*t^3 + b*t^2 + c*t + d

where a , b, c & d are the constants, Cp & t are Specific Heat & temperature resp.

3. As the observed values of Cp and temperature are provided and stored in txt file locally, we need to import data in file to python. This task is done by defining another function named "read_file()". Under this function, blank array for temperature and Cp are considered. For reading file line to line, a FOR LOOP is considered as below;

for line in open('data.txt','r'):

where, data.txt is file name containing observed data for Cp & t. 'r' denotes read only command. Each line in file consists of values t & Cp separated ','. We need extract these values separately so that we can store it in respective array. For this purpose line.split(',') function is used which splits line from comma(,). Hence we get separate values of Cp & t which are appended in array where values[0] i.e. first column is values of t and values[1] i.e. second column is values of Cp.

4. After defining functions, its time now to call this function in order to calculate constants a & b and thereby defining relationship between Cp and t.

temperature, cp = read_file()

5. Then by using the “curve_fit” function we calculate values of parameters a & b, as well as co-varience matrix and using the “function_1” we calculate the predicted Cp for both linear and cubic curves.

popt, pcov = curve_fit(function_1,temperature,cp)

fit_cp_l = function_1(np.array(temperature), *popt)

6. After that we print the values of parameters a & b, as well as co-varience matrix.

7. Then graph is plotted for the original and predicted values using linear and cubic curves.

8. After that we calculate the goodness parameters for both the fits and display the parameters in form of a table. The table contains SSE, R-squre, Adjusted R_squre and RMSE data for both the curves.

RESULTS:

1. This plot shows the comparison of the linear fit with original plot.

2. This plot shows the comparison of cubic fit with original plot .

3. This is the output in the command window , which shows the values of parameters i.e. popt, the covarience matrix or error matrix i.e. pcov and the goodness factors for each curve fit.

QUESTIONS:

1. What does popt and pcov mean?

ANSWER: 'popt' represents the matrix which stores and extracts the coefficients for the fifitting functions according to the general equations defined.

popt : array, Optimal values for the parameters so that the sum of the squared residuals of f(xdata, *popt) - ydata is minimized

pcov: As the name suggests pcov is co-varience matrix i.e. error matrix for the t parameters. From that we can calculate values of standard deviation of parameters. Standard deviation is nothing but deviation of parameters from its mean value. This deviation can be on positive or negative side, hence its necessary to consider its square values, which is referred as Co-varience. It is denoted by σ². Co-varience matrix in curve fifitting shows uncertainty and interdependence of each parameters. The diagonal element of matrix are co-varience or square of standard deviation of t parameters while non-diagonal elements of matrix shows relationship between the parameters. Lesser are the values of errors, more the slope of t curve matched with curve obtained from observed data.

pcov: 2D array, The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov).

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

ANSWER: This command converts each value of temperature in the data le into an array using numpy module.

3. What does the * in *popt mean?

ANSWER: '*' , '*popt' indicates and returns each coefficients stored in popt array. The * notation will expand a list of values into the arguments of the function.

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

ANSWER: Perfect fit means, the curve should fit the original curve without showing any errors in that particular degree of polynomial. This means that every data point in the original plot should be on the curve. Perfect fit can always be a best fit but best fit cannot be a perfect fit.

We can fit a curve perfectly only by using interpolation. It is not necessary that we get a single function representing complete trend.

To improve the cubic fit we can do split wise fitting for the curve. That is splitting the curve into multiple parts and using a different cubic polynomial for different sets of temperature and specific heat values. This reduces the errors and best t is obtained.

Here i have splitted the data into 5 regions and performed a split fitting for the cubic function .

Python program for split wise fitting of cubic function:

# Program to show split curve fitting
import math 
import matplotlib.pyplot as plt 
import numpy as np 
from scipy.optimize import curve_fit 
 
# Curve fitting functions

# Cubic function  
def function_1(t,a,b,c,d):  
	return a*pow(t,3) + b*pow(t,2) + c*t +d   
 
# 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 = []  
 
# Storing temperature and cp in the form of lists
temperature , cp = read_file() 

# Spliting data into 5 parts 
temperature1 = temperature[1:640];
cp1 = cp[1:640];

temperature2=temperature[641:1280];  
cp2 = cp[641:1280]; 
 
temperature3=temperature[1281:1920];  
cp3 = cp[1281:1920];   

temperature4=temperature[1921:2560];  
cp4 = cp[1921:2560];  
 
temperature5=temperature[2561:3200]; 
cp5 = cp[2561:3200]; 

# Considering cubic polynomial and calculating predicted value of Cp for 5 parts 
# Cp=a*t^3+b*t^2+ct+d  
popt, pcov = curve_fit(function_1,temperature1,cp1) 
fit_cp_l = function_1(np.array(temperature1), *popt) 

popt, pcov = curve_fit(function_1,temperature2,cp2) 
fit_cp_2 = function_1(np.array(temperature2), *popt) 

popt, pcov = curve_fit(function_1,temperature3,cp3) 
fit_cp_3 = function_1(np.array(temperature3), *popt) 

popt, pcov = curve_fit(function_1,temperature4,cp4) 
fit_cp_4 = function_1(np.array(temperature4), *popt) 

popt, pcov = curve_fit(function_1,temperature5,cp5) 
fit_cp_5 = function_1(np.array(temperature5), *popt) 

# Plotting split curve fit
plt.figure(1) 
plt.plot(temperature,cp,'y',linewidth = 4 ) 
plt.plot(temperature1,fit_cp_l,'k--',linewidth = 1) 
plt.plot(temperature2,fit_cp_2,'k--',linewidth = 1) 
plt.plot(temperature3,fit_cp_3,'k--',linewidth = 1) 
plt.plot(temperature4,fit_cp_4,'k--',linewidth = 1) 
plt.plot(temperature5,fit_cp_5,'k--',linewidth = 1) 
plt.legend(['Actual data','Split Curve fit']) 
plt.xlabel('Temperature (K)') 
plt.ylabel('Cp') 
plt.title('Split Curve fitting') 
plt.show() 

RESULT OF THE SPLIT WISE FITTING PROGRAM:

The yellow line (thicker line) represents the original plot and the black dotted line (thin line) represents the cubic fit.

As we can see by using splitwise fitting we have improved the cubic fit substantially. We can check by comparing it with the earlier cubic plot which we got without using splitwise fitting.

CONCLUSION:

1. The linear and polynomial curves are fitted to specific heat data provided using Python.
2. Goodness parameters for each fitness curve is evaluated.
3. The Cubic curve fit is good fit than the linear curve fit, and this is verified by the goodness parameters calculated as we can see the R-square value for the Linear fit is 0.9248 and for Cubic fit is 0.9907 which is more close to 1.