Curve fitting Using Python

Aim: Performing Curve Fitting in PYTHON

Objective:

1. Getting best fit for given set of data
2. Understanding the parameters to check the goodness of fit

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. After fitting data with one or more models, we should evaluate the goodness of fit.

Goodness of Fit Statistics

After using graphical methods to evaluate the goodness of fit, we should examine the goodness of fit statistics. Some of them are as follows:

• Residuals: The residuals from a fitted model are defined as the differences between the response data and the fit to the response data at each predictor value. `residual = data – fit`.Residuals are also known as error which are calculated as `"Error"=|y(i)-"value of" f(x)"at" x_i |=|y(i)-f(x_i)|`.
• The sum of squares due to error (SSE): A value closer to 0 indicates that the model has a smaller random error component, and that the fit will be more useful for prediction. It is calculated as follows: `SSE= sum_(i=0)^N(("Error"(i)))^2`.
• SSR (Sum of Square regression) `SSR= sum_(i=0)^N(f(x_i )-Mean)^2`.
• SST (Sum of squared total): `SST=SSE+SSR`
• R-square: R-square can take on any value between 0 and 1, with a value closer to 1 indicating that a greater proportion of variance is accounted for by the model. For example, an R-square value of 0.8234 means that the fit explains 82.34% of the total variation in the data about the average. Mathematically `R^2=(SSR)/(SST)`.
• 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=s=sqrt((SSE)/N)`.

  where N is the total number of data points available.

Problem Statement:

1. Perform curve fitting on specific heat and temperature data provided in the file.
2. What does popt and pcov mean? (Syntax of curve fit in Python)
3. What does np.array(temperature) do?
4. What does the * in *popt mean?
5. Write code to fit a linear and cubic polynomial for the Cp data. Explain if your results are good or bad.
6. What needs to be done in order to make the curve fit perfect?

Program Explanation

1. Program started with importing the necessary modules which were pyplot, numpy, scipy and tabulate.
2. In order to perform curve fit, we need a function which relates closely to the generated experimental curve, so far now as we don’t have the experimental curve yet, it was assumed that the variation would be linear or quadratic or cubic or fourth order. So, four different fitness functions with their corresponding equations were defined. All these functions took input arguments (which were temperature and coefficient values) and returned the function value at that input values.
3. The file which was provided contained two column with values cp and temperature, so a function named read_file() was defined that took the filename as input argument and returned two output arguments named temperature and cp, which contained the experimental values. Inside the function, for loop was implemented in order to extract these values line by line from the file and the values were stored in a variable of cp and temperature.
4. Then the curve_fit function was called which took the fitting function (linear, quadratic, cubic or fourth in our case), x values (temperature in our case) and y values (cp in our case) and returned popt and pcov as output arguments. popt provided the optimum values of coefficients of the fitness function, while pcov provided the covariance in the curve fit data.
5. In order to plot the curve which was calculated by the curve_fit function, the function values were calculated again by calling the fitness function with input arguments as array(temperature)(this meant that for all the elements contained in the array) and *opt(this meant to consider all the elements available in the list as input arguments).
6. The plots were made with both the curves (actual curve, and the curve developed by curve_fit function). This gave us a visual representation of how close the curve is to the actual curve. These steps were followed for all the fitness function and corresponding graphs were plotted.
7. In order to get to know the fitness characteristics of curve, SSE (Sum of square error), SSR (Sum of square Regression), R Square and RMSE (Root mean square error) were calculated by using the formulas from the theory and were displayed on the screen. Theses values were displayed in the form of table which was developed by using tabulate function.

Program

"""
	Curvefit by Saurav Sankhe
"""
# Importing modules
import numpy as np 
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit
from tabulate import tabulate

# Curve Fit function
# Linear function
def linear(t,a,b):
	return a*t + b

# Quadrtic function
def quadratic(t,a,b,c):
	return a*pow(t,2) + b*t + c

# Cubic function
def cubic(t,a,b,c,d):
	return a*pow(t,3) + b*pow(t,2) + c*t + d

# Fourth order Function
def fourth(t,a,b,c,d,e):
	return a*pow(t,4) + b*pow(t,3) + c*pow(t,2) + d*t + e

# 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
# =============================================================================
# For Linear Fit
temperature, cp = read_file()
popt, pcov_l = curve_fit(linear,temperature,cp)

linear_cp = linear(np.array(temperature), *popt)


plt.figure()
plt.plot(temperature, cp,color='blue',linewidth=2)
plt.plot(temperature, linear_cp,color='red',linewidth=2)
plt.legend(['Actual data','Quadratic fit'])
plt.xlabel('temperature [K]')
plt.ylabel('Cp')
plt.savefig('1_Linear fit')


# Main program
# For Quadratic Fit
popt, pcov_q = curve_fit(quadratic,temperature,cp)


quadratic_cp = quadratic(np.array(temperature), *popt)


plt.figure()
plt.plot(temperature, cp,color='blue',linewidth=2)
plt.plot(temperature, quadratic_cp,color='red',linewidth=2)
plt.legend(['Actual data','Quadratic fit'])
plt.xlabel('temperature [K]')
plt.ylabel('Cp')
plt.savefig('2_Quadratic fit')


# Main program
# For Cubic Fit
popt, pcov_c = curve_fit(cubic,temperature,cp)


cubic_cp = cubic(np.array(temperature), *popt)


plt.figure()
plt.plot(temperature, cp,color='blue',linewidth=2)
plt.plot(temperature, cubic_cp,color='red',linewidth=2)
plt.legend(['Actual data','Cubic fit'])
plt.xlabel('temperature [K]')
plt.ylabel('Cp')
plt.savefig('3_Cubic fit')

# Main program
# For fourth order polynomial
popt, pcov_f = curve_fit(fourth,temperature,cp)


fourth_cp = fourth(np.array(temperature), *popt)


plt.figure()
plt.plot(temperature, cp,color='blue',linewidth=2)
plt.plot(temperature, fourth_cp,color='red',linewidth=2)
plt.legend(['Actual data','Fourth order fit'])
plt.xlabel('temperature [K]')
plt.ylabel('Cp')
plt.savefig('4_Fourth Order fit')


# =============================================================================
# # Comparison of Fits
# =============================================================================
plt.figure()
plt.plot(temperature, cp,color='blue',linewidth=4)
plt.plot(temperature, quadratic_cp,color='red',linewidth=2,linestyle='-')
plt.plot(temperature, cubic_cp,color='green',linewidth=2,linestyle=':')
plt.plot(temperature, fourth_cp,color='yellow',linewidth=2,linestyle='-.')
plt.legend(['Actual data','Linear fit','Cubic fit','Fourth order fit'])
plt.xlabel('temperature [K]')
plt.ylabel('Cp')
plt.savefig('5_Comparison')


# =============================================================================
# # Error Characteristics
# =============================================================================
# Mean value
cp_mean = np.mean(cp)
print('Mean Value of cp is',cp_mean,'\n')
length = np.size(cp)


# For Linear Fit
residual_linear = np.array(cp) - np.array(linear_cp)
SSE_linear = np.sum(pow(np.array(residual_linear),2))

SSR_linear = np.sum(pow((np.array(linear_cp) - cp_mean),2))

SST_linear = SSE_linear + SSR_linear
RSquare_linear = SSR_linear/SST_linear

RMSE_linear = pow((SSE_linear/length),0.5)


# For Quadratic Fit
residual_quadratic = np.array(cp) - np.array(quadratic_cp)
SSE_quadratic = np.sum(pow(np.array(residual_quadratic),2))

SSR_quadratic = np.sum(pow((np.array(quadratic_cp) - cp_mean),2))

SST_quadratic = SSE_quadratic + SSR_quadratic
RSquare_quadratic = SSR_quadratic/SST_quadratic

RMSE_quadratic = pow((SSE_quadratic/length),0.5)


# For Cubic Fit
residual_cubic = np.array(cp) - np.array(cubic_cp)
SSE_cubic = np.sum(pow(np.array(residual_cubic),2))

SSR_cubic = np.sum(pow((np.array(cubic_cp) - cp_mean),2))

SST_cubic = SSE_cubic + SSR_cubic
RSquare_cubic = SSR_cubic/SST_cubic

RMSE_cubic = pow((SSE_cubic/length),0.5)



# For Fourth Order Fit
residual_fourth = np.array(cp) - np.array(fourth_cp)
SSE_fourth = np.sum(pow(np.array(residual_fourth),2))

SSR_fourth = np.sum(pow((np.array(fourth_cp) - cp_mean),2))

SST_cubic = SSE_fourth + SSR_fourth
RSquare_fourth = SSR_fourth/SST_cubic

RMSE_fourth = pow((SSE_fourth/length),0.5)


# =============================================================================
# # Printing the output in the form of table using tabulate module
# =============================================================================
table = [["Linear",SSE_linear,SSR_linear, RSquare_linear, RMSE_linear],\
["Quadratic",SSE_quadratic,SSR_quadratic, RSquare_quadratic, RMSE_quadratic],\
["Cubic",SSE_cubic,SSR_cubic, RSquare_cubic, RMSE_cubic],\
["Fourth Order",SSE_fourth,SSR_fourth, RSquare_fourth, RMSE_fourth]]

headers = ['Type','SSE','SSR','R Square','RMSE']

print(tabulate(table, headers, tablefmt="github"))

Output

Plots Generated

As it can be seen from the above table that as the degree of the fitness function increases, the curve gets closer to the actual curve. This can be depicted from the above graphs as well as the values of R Square (value gets close to 1) from the table given in Command prompt.

Conclusion:

1. What does popt and pcov mean?
   • These are the return values from the curve_fit function available in scipy module.
   • popt it is an array of optimal values for the parameters to the fitness function in which the sum of the squared residuals of model is minimized. These are nothing but the values of the coeeficients which are calculated with optimum fits. The number of elements in thses array depends upon the degree of the function.
   • pcov is a 2-dimensional array of the covariance of the popt. The diagonals provide the variance of the parameter estimate. Covariance indicates the uncertainties and correlations between parameters. This is most useful when the data has uncertainties.
2. What does np.array(temperature) do?
   • np is the short form provided while importing the numpy module in python. Also, while extracting the values of temperature in python, by default, these temperature values are stored in a datatype called as list. In order to perform the necessary operations on these values, it is necessary to convert it into a data type of array. So, array(temperature) converts this data type. Also, it helps in elements wise operations in any calculations.
3. What does the * in *popt mean?
   • It means that the star (*) is used as a reference or as an argument to a function. The star (*) notation will expand a list of values into the arguments of the function. This is useful if the function has more than one or two parameters. In general, the star (*) notation is used so that the variable which is placed after the star (*) is being referenced.
4. What needs to be done in order to make the curve fit perfect?
   • The curve fit can never be perfect. The main reason for that is the actual experimental value and the ideal curve-fit value always vary due to random errors. However, the value of actual data may come close to the ideal curve fit value. This can be seen by the code and graph shown below. As the order of the equations increases, the curve gets closer to the experimental data. There are other methods like dividing the data and curve fitting them separately or extrapolation.

References

1. Curve_fit function from scipy module
2. Tabulate library