Week 5 - Curve fitting

Aim:

To write code in python for curve fitting.

Objective:

1. To write code to fit a linear and cubic polynomial for the Cp data.
2. To plot the linear and cubic polynomial fit curves along with the initial data points.
3. To measure the fitness characteristics for both the curves.

Theory:

Curve fitting is the way we model or represent a data spread by assigning a ‘best fit‘ function (curve) along with the entire range. Ideally, it will capture the trend in the data and allow us to make predictions of how the data series will behave in the future.

Types of curve fitting include:

• Interpolation, where you discover a function that is an exact fit to the data points. Since this assumes no measurement error, it has limited applicability to real-life scenarios.
• Smoothing is when we find a function that is an approximate fit to the data points, but we give room for error and we allow our actual points to be near, but not necessarily on the line; given the error is minimized overall.

Linear and Polynomial Curve fitting:

Linear curve fitting, or linear regression, is when the data is fit to a straight line. Although there might be some curve to your data, a straight line provides a reasonable enough fit to make predictions. Since the equation of a generic straight line is always given by `f(x) = ax + b`, the question becomes: what a and b will give us the best fit line for our data?

Considering the vertical distance from each point to a prospective line as an error, and summing them up over our range, gives us a concrete number that expresses how far from ‘best’ the prospective line is.

A line that provides a minimum error can be considered the best straight line.

Since it’s the distance from our points to the line we’re interested in—whether it is positive or negative distance is not relevant—we square the distance in our error calculations. This also allows us to weight greater errors more heavily. So this method is called the least square approach.

Polynomial curve fitting is when we fit our data to the graph of a polynomial function. The same least square method can be used to find the polynomial, of a given degree, that has a minimum total error.

To measure the goodness of my given fit curve criteria we usually need four parameters to choose the best fit for the curve and also tell us how well the equation representing the given data points:

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

Governing Equation Used:

1. Error`(i)` = `|y(i) - f(x(i))|`,
2. `SSE` (Sum of square error) = `sum_(i = 1)^n`(Error`(i))^2`,
3. `SSR` (Sum of a square of regression) = `sum_(i=1)^n (f(x(i)) - Mean)^2`,
4. `SST` (Sum of the squared total) = `SSE + SSR`,
5. `R^2` (R - squared) = `(SSR)/(SST)`,
6. `RMSE` (Root mean squared error) = `sqrt((SSE)/n)`,

where `n` is the total number of datapoints available.

The initial data of Temperature (K) and specific heat (kJ/kcalK) is given by: https://drive.google.com/file/d/1vWi9Vj3AdQf0gqUK1u9jwS8zg92BHGR0/view

Solution Procedure:

1. At first, the required modules are imported viz., math, numpy, matplotlib, scipy to perform their respective functions.
2. After that two new functions are defined which return linear and cubic equations.
3. The data of Cp vs Temp. the file is loaded in this Python program from the right directory and a function read_file is created.
4. For loop is used to run each line by using an "open" command with "r" to read it.
5. From the data, we derived two individuals column matrices of temperature and Cp.
6. For curve fitting, 'popt' and 'pcov' are used as output for storing parameters obtained by 'curve_fit' which takes input as the required function (linear or cubic polynomial), temperature, and Cp.
7. Now another variable named 'fit_cp' stores the new values of Cp calculated by using temperature array.
8. Plots of both initial dataset and predicted dataset are made and compared to find fitness characteristics for both the curves viz., linear and cubic polynomial.
9. Various parameters such as `SSE, SSR, SST, R^2, RMSE` are calculated to measure the fitness characteristics for both the curves.
10. The plot which has the highest `R^2` or closet to 1 is considered to be a good fit.

Python Code:

# A program to 	measure fitness characteristics for the linear and cubical polynomial for the Cp data

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

# Curve fit function
# Linear function
def function_1(t,a,b):
	return a*t + b
# Cubic function
def function_3(t,a,b,c,d):
	return a*t**3 + b*t**2 + c*t + d

# Reading thermodynamics 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(function_1,temperature,cp)
fit_cp_l = function_1(np.array(temperature), *popt)
plt.figure(1)
plt.plot(temperature, cp, 'm--')
plt.plot(temperature, fit_cp_l, color = 'green', linewidth = 1)
plt.legend(['Actual data', 'Curve fit'])
plt.xlabel('Temperature (K)', color = 'maroon')
plt.ylabel('Cp (kJ/kcalK)', color = 'maroon')
plt.title('Linear Curve fitting', color = 'blue')
plt.grid()
plt.show()


popt,pcov = curve_fit(function_3,temperature,cp)
fit_cp_c = function_3(np.array(temperature), *popt)
plt.figure(3)
plt.plot(temperature, cp, 'm--')
plt.plot(temperature, fit_cp_c, color = 'green', linewidth = 1)
plt.legend(['Actual data', 'Curve fit'])
plt.xlabel('Temperature (K)', color = 'maroon')
plt.ylabel('Cp (kJ/kcalK)', color = 'maroon')
plt.title('Cubic Curve fitting', color = 'blue')
plt.grid()
plt.show()


'''
Measuring the fitness characteristics of all curves

SSE (sum of error squared)
SSR (sum of squares of the regression)
SST = SSE + SSR
R^2 = SSR/SST
RMSE = (SSE/l)^0.5

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

s = np.sum(cp)
l = np.size(cp)
m = s/l
print('Mean of all cp :', m)
print('')
# For Linear curve fit

linear_sse = 0
linear_ssr = 0
for i in range(l):
	linear_error = abs((np.sum((cp[i] - fit_cp_l[i]))))
	linear_sse = linear_sse + (linear_error**2)
	linear_ssr = linear_ssr + np.sum(pow((fit_cp_l[i] - m),2))

linear_sst = linear_sse + linear_ssr
print('linear_sst :', linear_sst)
linear_R2 = linear_ssr/linear_sst
print('linear_R2 :', linear_R2)
linear_rmse = (linear_sse/l)**0.5
print('linear_RMSE :', linear_rmse)
print('')

# For Cubic curve 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 + (cubic_error)**2
	cubic_ssr = cubic_ssr + np.sum(pow((fit_cp_c[i] - m), 2))

cubic_sst = cubic_sse + cubic_ssr
print('cubic_sst :', cubic_sst)
cubic_R2 = cubic_ssr/cubic_sst
print('cubic_R2 :', cubic_R2)
cubic_rmse = (cubic_sse/l)**0.5
print('cubic_RMSE :', cubic_rmse)

Results:

The following outputs are the results of above program...

The curve obtained from the program are:

The Linear Curve fitting plot is.

The Cubic Curve fitting plot is.

Now the answers to the following questions:

Q1. What does 'popt' and 'pcov' mean?

Ans. 'popt' represents the matrix that stores and extracts the coefficients for the fitting functions according to the equation generally defined.

popt: array

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

pcov: 2-D array

The estimated covariance of popt. The diagonals provide the variance of the parameter estimate. To compute one standard deviation errors on the parameters used perr = np.sqrt(np.diag(pcov)).How the sigma parameter affects the estimated covariance depends on absolute_sigma argument. As if true, sigma is used in an absolute sense and the estimated parameters covariance pcov reflects these absolute values. If false, only the relative magnitudes of the sigma values matter. The returned parameter covariance matrix pcov is based on scaling sigma by a constant factor.

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

Ans. Using numpy module, this command converts each value of temperature in the data file into an array.

3. What does the * in *popt mean?

Ans. * in *popt indicates the return of each coefficient stored in popt array.

4. Write code to fit a linear and cubic polynomial for the Cp data. Explain if your results are good or bad.

Ans. I already mention the code above and we could see that on being compared with linear and cubic polynomial we could see from the resulting plot that cubic fit better than a linear curve.

5. What needs to be done in order to make the curve fit perfectly?

Ans. There may be many ways in which curve can be fit perfectly and among them, two ways are:-

• Firstly we could find the curve fit more when we increase the order of the polynomial, the error produced will be less and the curve fit perfectly. As the value of `R^2` tends to 1 and we assumed to be a perfect fit.
• The previous method sometimes becomes cumbersome as we need to define functions every time for another polynomial. So another method is that one of the data says temperature can be split into multiple domains (let say) temperature column has 5000 data in it, we can split it into 10 domains of 500 values and if we plot them it will fit perfectly on a curve.

Conclusion:

Hence we conclude that to make the curve fit good `R^2` should tend to 1. And from the result, we could see that the value of `R^2` is approximately equal to 0.924877 for linear and 0.99071 for cubic and it is nearly 0.99. Therefore we can say that the cubic polynomial fit better than linear polynomial. If the order of polynomial increases than the value of `R^2` will inch towards 1 and hence could consider a more fit curve than lower-order polynomial.

Reference:

1. https://skill-lync.com/projects/Curve-Fitting-and-criterion-to-choose-best-fit-41860
2. https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html#scipy.optimize.curve_fit