Curve fitting using Python

The project aims at finding the best fit curve for a given data and finding the error between the original and said curve.

The data being used is that of Specific heat (Cp) vs Temperature. The plot will be made to show the variation of Cp with Temp first. 

Code:

import numpy as nm
from numpy import array
from scipy.integrate import odeint
from scipy.optimize import curve_fit
from statistics import mean
import matplotlib.pyplot as plt
import math
import time

# Curve fit- Linear Poly
def curflin(t, a, b):
	return a*t + b

# Curve fit- quadratic Poly
def curfquad(t, a, b, c):
	return a*pow(t, 2) + b*t + c

# Curve fit- cubic Poly
def curfcub(t, a, b, c, d):
	return a*pow(t, 3) + b*pow(t, 2) + c*t + d

# Curve fit- 4th degree Poly
def curf4(t, a, b, c, d, e):
	return a*pow(t, 4) + b*pow(t, 3) + c*pow(t, 2) + d*t + e

# Read file
def readfile():
	temp = []
	cp = []
	for line in open('data', 'r'):
		values = line.split(',')
		temp.append(float(values[0]))
		cp.append(float(values[1]))

	return [temp, cp]

# Main code
temp, cp = readfile()

# Linear
popt, pcov = curve_fit(curflin, temp, cp)
fitcp1 = curflin(nm.array(temp), *popt)

# Quadratic
popt, pcov = curve_fit(curfquad, temp, cp)
fitcp2 = curfquad(nm.array(temp), *popt)

# Cubic
popt, pcov = curve_fit(curfcub, temp, cp)
fitcp3 = curfcub(nm.array(temp), *popt)

# 4th Degree
popt, pcov = curve_fit(curf4, temp, cp)
fitcp4 = curf4(nm.array(temp), *popt)

print(fitcp1)


# Errors Calculation begins here

# Mean for SSR (sum of squares of regression)
meanavg = mean(cp)

sse = array([0, 0, 0]) # Sum of squared errors
ssr = array([0, 0, 0]) # Sum of squared regression

sserror1 = []
sserror2 = []
sserror3 = []
ssregress1 = []
ssregress2 = []
ssregress3 = []

i = 0
l = len(cp)
for i in range(l): 
	sserror1.append(cp[i] - fitcp1[i]) # sum of suard errors = Difference between actual data points and fitting function
	sserror2.append(cp[i] - fitcp2[i])
	sserror3.append(cp[i] - fitcp3[i])

	#sse1 = sse1 + pow(sserror1[i], 2)
	#sse1 = sse2 + pow(sserror2[i], 2)
	#sse1 = sse3 + pow(sserror3[i], 2)
	sse[0] = sse[0] + pow(sserror1[i], 2)
	sse[1] = sse[1] + pow(sserror2[i], 2)
	sse[2] = sse[2] + pow(sserror3[i], 2)

	ssregress1.append(fitcp1[i] - meanavg) # sum of suared regression = difference between points of predicted function and mean
	ssregress2.append(fitcp2[i] - meanavg)
	ssregress3.append(fitcp3[i] - meanavg)
	
	ssr[0] = ssr[0] + pow(ssregress1[i], 2)
	ssr[1] = ssr[1] + pow(ssregress2[i], 2)
	ssr[2] = ssr[2] + pow(ssregress3[i], 2)
  	
print('Sum of squared errors (SSE): ', sse)
print('Sum of squares of regression (SSR): ', ssr)

#array(sse)
#array(ssr)

sst = array([0, 0, 0])
sst = array(sse) + array(ssr) # using "nm.array()" is also possible in place of "array()"

print('SST = SSR + SSE = ', sst)

rsquare = ssr/sst # measure of successful fitting
print('R Square:', rsquare)

# Adjusted r square
rsqadj = 1 - (((1-rsquare)*(3200-1))/(3200 - 1 - 1))
print('R Square Adjusted:', rsqadj)

# Root mean square errors = square root of sse divided by nuumber of data points available

rmse = pow((sse/3200), 0.5)

plt.plot(temp, cp, color = 'black', linewidth = 3)
plt.plot(temp, fitcp1, color = 'green', linewidth = 2)
plt.plot(temp, fitcp2, color = 'red', linewidth = 2)
plt.plot(temp, fitcp3, color = 'blue', linewidth = 2)
plt.plot(temp, fitcp4, color = 'yellow', linewidth = 2)
plt.legend(['Actual Data', 'Linear fit', 'Quad fit', 'Cubic Fit', '4th Degree fit'])
plt.xlabel('Temperature (K)')
plt.ylabel('Cp')
plt.show()

Best Fit: Refers to the curve that fits the data and has the least error with respect to the original curve.

Perfect Fit: Perfect fit is the fit with zero error, or as we will calculate later, an R^2 value equal to 1. A perfect fit isn't possible to achieve as there will always be approximations, and some unrealistic points as mentioned in the linked article: Curve Fitting; best fit and perfect fit. The higher the degree of the polynomial, the higher the R^2 value but greater is the approximation, which is oftentimes shown as a warning in the command window. Basically, the polynomial is sensitive to change.

The goal here is to check the goodness of the fit and assess the error.

1. First, the given Cp vs temp data is read and extracted in 2 lists Cp and temp.
2. Then, the function is defined for each degree polynomial, and the same is passed to curve_fit command to find the coefficients, predicted fit and optimal solution.
3. np.array(temp) converts the temp list to an array so that it can be passed to the functions.
4. The * in popt will expand a list of values into the arguments of the function. This is useful in our case where the function has more than one or two parameters. 

What is popt:

Optimal values for the parameters so that the sum of the squared residuals of "f(xdata, *popt) - ydata" is minimized. Basically, popt stores the best fit values of the coefficients of the equation.

What is pcov:

This is the covariance matrix, which indicates the uncertainties and correlations between parameters. This is most useful when the data has uncertainties.

Error calculation:

1. First, we take 2 arrays One is sse (sum of squared errors) which is the square of the difference between original data and fit curve, and other is ssr (sum of squared regression) which is the difference of fit curve from the mean of original data.
2. 3 values in the array are to store errors for linear, quadratic and cubic vectors.
3. In a for loop from 1 to length of the cp data (Which is 3200), we calculate sse and ssr in every iteration and store it in the respective arrays.
4. sse, ssr and sst, which is the sum of sse and ssr, are returned.
5. R^2 is the measure of goodness of the fit, which is calculated as "ssr/sst". This ranges from 0 to 1, with 1 indicating perfect fit.
6. Additionally, Adjusted R^2 and RMSE (root mean square error) are calculated, which do the adjustment to account for independent variables.

Output:

Sum of squared errors (SSE):  [2161516  538534   92718]

Sum of squares of regression (SSR):  [26629000 28251972 28697788]

SST = SSR + SSE =  [28790516 28790506 28790506]

R Square: [0.92492264 0.98129474 0.99677956]

R Square Adjusted: [0.92489916 0.98128889 0.99677856]

RMSE: [25.98987784 12.97273583  5.38278506]

R^2 values for various fits:

1. Linear: 0.9249
2. Quadratic: 0.98129
3. Cubic: 0.99677

Ideally, increasing the degree of the polynomial increases the goodness of the fit and gets closer and closer to being the best fit. This is demonstrated above, where cubic shows the highest R^2 value amongst the 3, and we know that R^2 value is the means to judge the goodness of the fit.

However, we also know that we get warnings when we increase the degree to above 2. This is because of the approximations that the solver makes to fit the curve.

As it stands, the R^2 values obtained show that our results are good since they are very close to 1.

Plot:

Errors:

1- Calling error:

"fitc1" changed to "fitcp1"

2- pow syntax error:

Added comma.

Other References:

1. Polynomials
2. Introduction to Scipy