Modeling a polynomial to relate specific heat with temperature

Aim: Fitting a linear polynomial and a cubic polynomial to the given data set of Temperature vs Specific Heat.

Objective:

1) Curve Fitting

2) Determining the Goodness of the curve

3) Improving the cubic polynomial fit

Procedure:

Curve Fitting:

1) The data file is read using "load('filename')" command (line 6 - Main).

2) Temperature and Specific Heat Arrays were generated using the data (lines 7 & 8 - Main).

3) Matlab's polyfit function was used to fit the data and obtain coefficients for the curve(line 11 - Main).

syntax: polyfit(x,y,n) where x and y are the data arrays and n is the degree of polynimal required.

4) polyval fucntion was used to obtain prediction data set for the fit(line 12 - Main)

syntax: polyval(x,y)

5) The arrays thus obtained were plotted along with the original data set(lines 19 -26 - Main)

Warning:

Solution:

A three variable output for polyfit and a four variable input for polyval were used to center and scale the data to the standard deviation of the fit [1] (lines 15,16 - Main).

Main Program:

clear all
close all
clc

%Input Data
Data = load('data');
Temp = Data(:,1);
Cp(:,1) = Data(:,2);

%Linear Fit
coeff1 = polyfit(Temp,Cp(:,1),1);
Cp(:,2) = polyval(coeff1,Temp);

%Cubic Fit
[coeff2,S,M] = polyfit(Temp,Cp(:,1),3);
Cp(:,3) = polyval(coeff2,Temp,[],M);

%Plot
plot(Temp,Cp(:,1),'linewidth',2)
hold on
plot(Temp,Cp(:,2),'linewidth',2)
plot(Temp,Cp(:,3),'linewidth',2)
xlabel('temperature (K)')
ylabel('Specific Heat (J/mol-K)')
legend('Data points','Linear Fit','Cubic Fit','location','SouthEast')
title('Specific Heat vs Temperature')

figure
%%Mean square error and residual plotting
E(1) = deviation(Temp,Cp(:,2),Cp(:,1));
hold on
E(2) = deviation(Temp,Cp(:,3),Cp(:,1));
title('Residuals vs Temperature')

%Results
disp('Mean square error:')
disp(E)

Fitness Characteristics:

Characteristics of the curve can be determined by using the difference between the predicted values and the experimental values/original data set.

1) Residuals: The error at x(i) is the difference between predicted Y(i) and actual y(i)[2].

2) Mean Square Error: Is the mean of squares of the residuals. This parameter defines how perfectly the curve fits to the given data[2].

Deviation:

Deviation function was used to determine the mean squared error(E) and plot the residuals(R). (lines 30,32 - main)

function E = deviation(X,A,B)
%Residual
R = A-B;
%Residual Plot
plot(X,R)
xlabel('Temperature (K)')
ylabel('Residual')

%Linear Square Error
E = 0;
for i=1:length(A)
E = E + R(i)^2;
end
E = E/length(A);%mean squared error
end

Temperature, Cp(calculated), Cp(given) where given as inputs to the function.

Improving the fit:

Zero mean squared error would imply a perfect fit where in each and every point in the data set is exactly on the curve. This can be achieved by increasing order of the polynomial.[3]

Although that might seem like a meticulous way of fitting the curve one must understand that curve fitting is basically done to hypothesise a relation between the variables, and a perfect fit will lead to a higher error when using the relation on new data (even when the data is within the range of the training set) and will therefore make a bad model.

Best fit is one which can predict the outcome for a data point outside the training set with least amount of error. Cross validation is widely used for this purpose, where a part of the data is set aside for testing purpose. K-fold cross validation is a widely used method in statistics to determine the best fit. In this process data is divided into K folds(sets) and k iterations of cross validation are carried out by taking a different fold as the test set each time. Average mean square error for all iternations is taken and comapred to determine the optimum degree for the polynomial. Same procedure can be used to refine/improve the fit for a given degree of polynomial.[4]

Main 2:

The following program reads the data file(line 6) to deliver the Cp vs Temperature plots for the original data set along with cubic polynomials fit over the 5-fold data sets.

clear all
close all
clc

%Input Data
Data1 = load('data');
Data = rearrange(Data1);%Generating random sequence to obtain folds with random data sets

%Data manipulation generating 5 data folds
%set 1
Set1a = Data(1:4*length(Data(:,1))/5,:);
Set1b = Data(1+4*length(Data(:,1))/5:length(Data(:,1)),:); %Data spliting
Set1a = foldsort(Set1a);%Soritng
Set1b = foldsort(Set1b);

%set 2
Set2a = Data(1:3*length(Data(:,1))/5,:);
Set2a(1+3*length(Data(:,1))/5:4*length(Data(:,1))/5,:) = Data(1+4*length(Data(:,1))/5:length(Data(:,1)),:);
Set2b = Data(1+3*length(Data(:,1))/5:4*length(Data(:,1))/5,:); %Data splitting
Set2a = foldsort(Set2a);%sorting
Set2b = foldsort(Set2b);

%set 3
Set3a = Data(1:2*length(Data(:,1))/5,:);
Set3a(1+2*length(Data(:,1))/5:4*length(Data(:,1))/5,:) = Data(1+3*length(Data(:,1))/5:length(Data(:,1)),:);
Set3b = Data(1+2*length(Data(:,1))/5:3*length(Data(:,1))/5,:); %Data splitting
Set3a = foldsort(Set3a);%sorting
Set3b = foldsort(Set3b);

%set 4
Set4a = Data(1:length(Data(:,1))/5,:);
Set4a(1+length(Data(:,1))/5:4*length(Data(:,1))/5,:) = Data(1+2*length(Data(:,1))/5:length(Data(:,1)),:);
Set4b = Data(1+length(Data(:,1))/5:2*length(Data(:,1))/5,:); %Data splitting
Set4a = foldsort(Set4a);%sorting
Set4b = foldsort(Set4b);

%set 5
Set5a = Data(1+length(Data(:,1))/5:length(Data(:,1)),:);
Set5b = Data(1:length(Data(:,1))/5,:); %Data splitting
Set5a = foldsort(Set5a);%sorting
Set5b = foldsort(Set5b);

%Fitting with centering and scaling
[fit(:,6),~,M] = polyfit(Data1(:,1),Data1(:,2),3);
cp6 = polyval(fit(:,6),Data1(:,1),[],M);

[fit(:,1),~,m1] = polyfit(Set1a(:,1),Set1a(:,2),3);
cp1 = polyval(fit(:,1),Data1(:,1),[],m1);

[fit(:,2),~,m2] = polyfit(Set2a(:,1),Set2a(:,2),3);
cp2 = polyval(fit(:,2),Data1(:,1),[],m2);

[fit(:,3),~,m3] = polyfit(Set3a(:,1),Set3a(:,2),3);
cp3 = polyval(fit(:,3),Data1(:,1),[],m3);

[fit(:,4),~,m4] = polyfit(Set4a(:,1),Set4a(:,2),3);
cp4 = polyval(fit(:,4),Data1(:,1),[],m4);

[fit(:,5),~,m5] = polyfit(Set5a(:,1),Set5a(:,2),3);
cp5 = polyval(fit(:,5),Data1(:,1),[],m5);


%Curve Plot
plot(Data1(:,1),Data1(:,2),'linewidth',0.1)
hold on
plot(Data1(:,1),cp6,'linewidth',3.5)
plot(Data1(:,1),cp1,'linewidth',3)
plot(Data1(:,1),cp2,'linewidth',2.5)
plot(Data1(:,1),cp3,'linewidth',2)
plot(Data1(:,1),cp4,'linewidth',1.5)
plot(Data1(:,1),cp5,'linewidth',1)
legend('Data','FIT','Fit-1','Fit-2','Fit-3','Fit-4','Fit-5','location','southeast')
xlabel('Temperature')
ylabel('Specific Heat')
title('Fit Comparision Plots')

%Mean square error in Test Set and residual plotting
figure
hold on
E(6) = deviation(Data1(:,1),cp6,Data1(:,2),'Single Fold Data');
E(1) = deviation(Set1b(:,1),polyval(fit(:,1),Set1b(:,1),[],m1),Set1b(:,2),'Fold 5 as Validation Set');
E(2) = deviation(Set2b(:,1),polyval(fit(:,2),Set2b(:,1),[],m2),Set2b(:,2),'Fold 4 as Validation Set');
E(3) = deviation(Set3b(:,1),polyval(fit(:,3),Set3b(:,1),[],m3),Set3b(:,2),'Fold 3 as Validation Set');
E(4) = deviation(Set4b(:,1),polyval(fit(:,4),Set4b(:,1),[],m4),Set4b(:,2),'Fold 2 as Validation Set');
E(5) = deviation(Set5b(:,1),polyval(fit(:,5),Set5b(:,1),[],m5),Set5b(:,2),'Fold 1 as Validation Set');
title('Residual Comparision Plots')


[G,i] = greater(E);%Comparision

%Results
disp('Best fit coefficients:')
disp(fit(:,i))
disp('Mean square error for best fit:')
disp(G)

1) The data is first randomly arranged using the rearrange function(line 7-Main2):

function Data = rearrange(Data1)
a = randperm(length(Data1));
for i = 1:length(Data1)
    Data(a(i),:) = Data1(i,:);
end

Note: Array a had random integers from 1 to length of input array, by assigning the 'i' th value from inital data set to 'a(i)'th value of output data set, the array is shuffled row wise.

2) The rearranged data is then split into 5 folds where Set#a is the training set for '#'th iteration and Set#b is the testing/validation set.(lines 11,12,17,18,19,24,25,26,31,32,33,38 and 39 - Main2)

3)The folds are then sorted to achieve a continuity in the plots using foldsort function (lines 13,14,20,21,27,28,34,35,40 and 41 - Main2):

Note: Two loops were used to swap numbers in an ascending order. variable b was used to facilitate the swap.

function A = foldsort(A)
i = 1;
j = 1;
for i = 1:length(A(:,1))
    for j = 1:length(A(:,1))
        if A(i,1) < A(j,1);
            b = A(i,:);
            A(i,:) = A(j,:);
            A(j,:) = b;
        end
    end
end

4) The training arrays thus obtianed were used to fit cubic equations using polyfit function (lines 44,47,50,53,56 and 59 -Main2)

5) Prediction arrays were generated using polyval function and then plotted.

6) Devation function was called to plot the residuals and return the mean squared errors for all the fits over the testing sets generated in step 2. (lines 80-85-Main 2)

7) Returned mean square errors were then compared to determine the best fit using greater function to get the fit with least mean square error (line 89 - Main 2):

function [B,j] = greater(A)
B = 1000000;
for i = 1:length(A)
    if B>A(i)
        B=A(i);
    end
end
for i = 1:length(A)
    if B == A(i)
        j = i;
    end
end
end

Note: Least mean square error (B) was found by storing the lesser of B , A(i) in B for all values of i. to determine the fit that gives B, the number of the fit(i) that has the least value was saved in j. B, j were returned to main2 and printed on screen (lines 92 to 95 - Main 2)

Results:

Main:

Main 2:

Discussion:

1) The plots show that increasing the degree of polynomial gives a better fitting curve

2) The Mean Squar Error comparison clearly shows that the fit characteristics can be improved by using K-fold cross validation over the same data set.

3) Higher numeber of folds can give better results.

Scope:

The data used in generating folds is randomly arranged when running the program and hence we get a different least mean square and its corresponding fit coefficients every time the program is run. Instead the program can be turned into a function with a few tweaks and called inside a loop to get a finer fit.

References:

1) https://in.mathworks.com/help/matlab/ref/polyfit.html

2) Emperical laws and Curve fitting - Higher Engineering Mathematics by B.S. Grewal, 44th Edition.

3)https://www.theanalysisfactor.com/assessing-the-fit-of-regression-models/

4)https://towardsdatascience.com/overfitting-vs-underfitting-a-complete-example-d05dd7e19765