Curve Fitting using Python

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:

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. It is highly effective in the mathematical modelling of some natural processes.

Curve fit are important to know the trend followed by data and also allow additional mathematical operation which might not be possible in the initial data set. For example, the population data in country. There is no formula to find what the population is going to be in next year. However, if we fit a curve, we had a function which say how population is going to change with respect to time. The accuracy of results depends on the quality of curve fit but at least we can have a rough trend.

Goodness of a fit:

Upon finding a curve which satisfies the data points there are 4 quantities which help us measure the goodness of fit criteria or how well the equation is a representing the datapoints. They are,

1. The sum of squares due to error (SSE)
2. R-square
3. Adjusted R-square
4. Root mean squared error (RMSE)

1 . The sum of squares due to error (SSE)

The Error Sum of Squares(SSE) is the sum of squares of the deviations of the data values from the predicted values (typically predicted from a Least Squares) I.e. it's a measure of how well the independent variable fit the data.

`Err(i)=Y_i-f(x(i))`

SSE(sum of squared errrors) = `∑_(i=1)^nErr(i)^2`

When we differentiate the SSE with each and every coefficient of the proposed equation f(x) we and equate it to zero, we can find the minima and the value of that coefficient.

Upon solving the multiple system of equations we can find the values of all the coefficients.

2 . R-square

This quantity that measures how successful the fit is in explaining the variation of the data. Put another way, R-square is the square of the correlation between the response values and the predicted response values. It is also called the square of the multiple correlation coefficient and the coefficient of multiple determination.

That is the higher the value of R-square the better the fit, the lower it is the lousier the fit, R-square ranges from a value of 0 to 1

• 0 for the worst possible fit for the data set
• 1 for the best possible fit for the data set

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

Where,

SSR = Sum of squares of the regression

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

SSE = `∑_(i=1)^n(Y_i-f(x(i)))^2`

SST = SSR + SSE

3 . 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.

• If value of contributes to the output, then is a significant dependent variable
• If value of does not contribute to the output, then is an insignificant dependent variable

R-square is formulated as shown below,

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

where,

n = the number of datapoints

k = the number of independent variables

4 . Root mean squared error (RMSE)

This statistic is also known as the fit standard error and the standard error of the regression. It is an estimate of the standard deviation of the random component in the data, and is defined as,

`RMSE =√((SSE)/n)`

where,

n is the total number of datapoints available.

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.