Week 5 - Curve fitting
Objective: To perform Curve Fitting for Linear and Cubic Polynomial for the Cp data.
Theory:
Curve fitting is the process of constructing a curve or mathematical function, that has the best fit to a series of data points, possibly subject to constraints. Curve fitting can involve either interpolation, where an exact fit to the data is required, or smoothing, in which a "smooth" function is constructed that approximately fits the data.
Parameters used to measure the fitness of a curve:
- The sum of squares due to error (SSE)
- R-square
- Adjusted R-square
- Root mean squared error (RMSE)
The syntax for the curve_fit is
curve_fit(f, xdata, ydata)
ydata = f(xdata, *params)
where,
f : i.e the model function, f(x, ...) must take the independent variable as the first argument and the parameters to fit as separate remaining arguments.
xdata : An M-length sequence or an (k,M)-shaped array for functions with k predictors. The independent variable where the data is measured.
ydata : M-length sequence, The dependent data — nominally f(xdata, ...)
4. Write code to fit a linear and cubic polynomial for the Cp data. Explain if your results are good or bad.
Program:
Code for the Linear Fitting:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Curve fit function
def linear(t, a, b):
return a*t + b
# 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(linear, temperature, cp)
fit_cp = linear(np.array(temperature),*popt)
plt.plot(temperature, cp,color="green",linewidth=3)
plt.plot(temperature,fit_cp,color="magenta",linewidth=3)
plt.title('Comparing Actual data and Linear fit')
plt.legend(['Actual data','Linear fit'])
plt.xlabel('Temperature [K]')
plt.ylabel('Cp')
plt.show()
# Goodness of the fit
error = np.array(cp) - np.array(fit_cp)
error_square = np.square(error)
regression = np.array(cp) - np.mean(cp)
ssr = np.sum(pow(regression,2))
sse = np.sum(error_square)
mse = np.mean(error_square)
sst = sse + ssr
r_square = ssr/sst
print('Goodness of fit for Linear polynomial is')
print(r_square)
Code for the Quadratic Fitting:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Curve fit function
def quad(t, a, b, c):
return a*pow(t,2) + b*t + c
# 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(quad, temperature, cp)
fit_cp = quad(np.array(temperature),*popt)
plt.plot(temperature, cp,color="blue",linewidth=3)
plt.plot(temperature,fit_cp,color="red",linewidth=3)
plt.title('Comparing Actual data and Quadratic fit')
plt.legend(['Actual data','Quadratic fit'])
plt.xlabel('Temperature [K]')
plt.ylabel('Cp')
plt.show()
# Goodness of the fit
error = np.array(cp) - np.array(fit_cp)
error_square = np.square(error)
regression = np.array(cp) - np.mean(cp)
ssr = np.sum(pow(regression,2))
sse = np.sum(error_square)
mse = np.mean(error_square)
sst = sse + ssr
r_square = ssr/sst
print('Goodness of fit for Quadratic polynomial is')
print(r_square)
Code for the Cubic Fitting:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Curve fit function
def cubic(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]
# Main Program
temperature,cp = read_file()
popt, pcov = curve_fit(cubic, temperature, cp)
fit_cp = cubic(np.array(temperature),*popt)
plt.plot(temperature, cp,color="magenta",linewidth=3)
plt.plot(temperature,fit_cp,color="blue",linewidth=3)
plt.title('Comparing Actual data and Cubic fit')
plt.legend(['Actual data','Cubic fit'])
plt.xlabel('Temperature [K]')
plt.ylabel('Cp')
plt.show()
# Goodness of the fit
error = np.array(cp) - np.array(fit_cp)
error_square = np.square(error)
regression = np.array(cp) - np.mean(cp)
ssr = np.sum(pow(regression,2))
sse = np.sum(error_square)
mse = np.mean(error_square)
sst = sse + ssr
r_square = ssr/sst
print('Goodness of fit for Cubic polynomial is')
print(r_square)
The Procedure followed while writing the above code is cited below:
- The code is written to load the file 'data' in python which contains the raw data of multiple cp and temperatures.
- Two variables 'temperature' and 'cp' are created and the values are assigned to them by reading the file using the open function and then extracting the information from each line which is run in the function 'read_file'.
- In this step, we create a function 'linear' with the parameters t,a,b where t is the variable and a,b are the coefficients of the linear equation.
- Then we use the 'curve_fit' function to get the curve fit by passing the linear function, temperature and the cp later storing them in popt and pcov.
- The 'fit_cp' is then calculated by using the linear function and adding the temperature and the popt as it's parameters. The '*' ahead of popt convey that all the data of the popt to be included.
- Then, code is written to generate the plots as mentioned with titles, axes with labels and legend which differentiates the curves from each other. The plots are shown at the end of the result.
The results are good but not perfect because the plots as we can observe below are close enough but not exactly the same. So we can conclude that the fit is not perfect.
Result:
The Plot for Actual Data and Linear Fit is shown above.
The Plot for Actual Data and Quadratic Fit is shown above.
The Plot for Actual Data and Cubic Fit is shown above.
1. What does popt and pcov mean?
popt - It is an array containing Optimal values for the parameters so that the sum of the squared error of f(xdata, *popt) - ydata(i.e fit_cp) is minimized. Refer to the above-mentioned syntax.
pcov - It is a two-dimensional array that helps in calculating the estimated covariance of popt. The diagonals provide the variance of the parameter estimate.
2. What does np.array(temperature) do?
The 'np.array' helps in converting the temperature data into arrays which can be utilized to perform array operations. NumPy is used to work with arrays.
3. What does the * in *popt mean?
When the '*' is utilized all the values contained in the popt are extracted and used in the corresponding statement.
5. What needs to be done in order to make the curve fit perfect?
The curve fit can be made close to perfect increasing the order of the polynomials to higher degree. Another way to accomplish this task is by splitting the data into multiple groups which would bring the closer to being perfect. Unlike the MatLab where there is toolbox provided for curve fitting it's difficult to achieve the perfect fit. But we can achieve the best fit in this case. The value of R^2 must be equal to 1 to obtain the perfect fit.
References:
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.curve_fit.html
Project Link: https://skill-lync.com/projects/curve-fitting-using-python-49