How to use polyfit to get the difference out of two functions?

Hi

I have two datasets (X,Y). For every value of Y that corresponds to a value of X.

I am thinking of applying a polyfit command to get the best curve fitting for the data, and as a result I will come up with a function (slope+intercept).

How do I find the difference between those two functions? What command do I apply to subtract Function A from Function B?

[Walter Roberson]

One of the functions is the one implied by polyfit(x,y,1) . But what is the other function?

A collection of points is not considered a function except in mathematical theory (in which case, the collection of ordered values is considered to define the function, and formulae are just shorthand for the collection of values.)

[Stelios Fanourakis]

@Walter. The datasets are curves. If I plot them in excel and apply a polynomial curve fitting of high degree (e.g. 6) then I get the best line fit and a function. I need to get the same function in Matlab. And then subtract one function from another since the datasets are very close each other, to find their difference. How do I do that?

Answer 1

[Walter Roberson]

P1 = polyfit((1:length(X)).', X(:), 6);

Likewise for Y producing P2

Pd = polyval(P1, 1:length(X)) - polyval(P2, 1:length(X))

[Stelios Fanourakis]

@Walter.

I guess this is the correct answer but I need a few clarifications.

The two datasets are comprised a series of X data that corresponds to a series of Y data.

I guess the polyfit command would be like

P1 = polyfit((1:length(X1)).', Y1(:), 6);
P2 = polyfit((1:length(X2).',Y2(:),6);

Pd = polyval(P1, 1:length(X)) - polyval(P2, 1:length(X2))

Am I doing something wrong? I need to estimate numerically the difference between two curves.

Can you also remind me what the ' means in programming?

[Stelios Fanourakis]

Is it correct as stated?

[Walter Roberson]

Supposing that X and Y have two columns, with X(K,1) corresponding to Y(K,1) and X(K,2) corresponding to Y(K,2) then

P1 = polyfit(X(:,1), Y(:,1), 6);
P2 = polyfit(X(:,2), Y(:,2), 6);
allX = union(X(:,1), X(:,2));
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd);

In MATLAB, ' is complex conjugate transpose, formal name ctranspose() . But I did not use it. I used .' which is regular transpose without complex conjugate, formal name transpose() .

polyfit() can accept X and Y that are both rows, or X and Y that are both columns, but it does not like a mix such as a row X = 1:20 and a column Y = rand(20,1) . I do not know the orientation of your dataset so I cannot assume that your Y is a row, so it is not safe for me to program polyfit(1:length(Y), Y, 6) because Y might be a column. Therefore I need to force both X and Y to be rows, or force both to be columns. Using (:) forces a column so Y(:) to force a column for Y is easy, but the two ways to force Y to be a row are Y(:).' or else reshape(Y,1,[]) which are a bit ugly. So it is easiest to force X to be a column by using the .' operator on what is known to be a row, 1:length(X), and use Y(:) . I could also have written

P1 = polyfit(1:length(X), Y(:).', 6)

Now, if someone happened to be writing this as part of a sequence of statements and had reason to know that Y was a row already, then they could just use

P1 = polyfit(1:length(X), Y, 6)

but since I do not know because the person asking the problem did not specify the shape, then I have to put in the protections to make the code work either way.

[Stelios Fanourakis]

You can say Y are rows and X columns. For every Y is a X and the same applies at both data sets

[Walter Roberson]

%"Y are rows and X columns"
P1 = polyfit(X(:,1), Y(1,:).', 6);
P2 = polyfit(X(:,2), Y(2,:).', 6);
allX = unique(X);
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd)

[Stelios Fanourakis]

@Walter

I use this code

clc;
clear all;
 
x1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','M2:M384');
y1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','A2:A376');
 
x2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','N2:N867');
y2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','B2:B868');
 
P1 = polyfit(x1(:,1), y1(1,:).', 6);
P2 = polyfit(x2(:,2), y2(2,:).', 6);
allX = unique(X);
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd)

And I get the error

Error using  .'  (line 191)
Undefined function 'transpose' for input arguments of type 'table'.

Error in curvefitdif (line 10)
P1 = polyfit(x1(:,1), y1(1,:).', 6);

[Stelios Fanourakis]

I also get the error

Error using polyfit (line 44)
X and Y vectors must be the same size.

Error in curvefitdif (line 10)
P1 = polyfit(x1(:,1), y1(1,:), 6);

Why?

[Walter Roberson]

You seem to be reading 384-2+1 = 383 values for x1

You seem to be reading 376-2+1 = 375 values for y1

You seem to be reading 867-2+1 = 866 values for x2

You seem to be reading 868-2+1 = 867 values for y2.

I do not know what it means to fit 383 x values against 375 y values.

When you use readtable(), you get back a table() object as the result, not a numeric array. You will need to use {} indexing to get to the numeric results, such as

x1 = x1{:,};  %replace the table object with numeric object it contains

[Walter Roberson]

x1r = readtable('ValidationTest.xls', 'Sheet',5, 'Range','M2:M384');
y1r = readtable('ValidationTest.xls', 'Sheet',5, 'Range','A2:A376');
 
x2r = readtable('ValidationTest.xls', 'Sheet',5, 'Range','N2:N867');
y2r = readtable('ValidationTest.xls', 'Sheet',5, 'Range','B2:B868');

h1 = min([size(x1r,1), size(y1r,1)]);
h2 = min([size(x2r,1), size(y2r,1)]);

x1 = x1r{1:h1, 1};
y1 = y1r{1:h1, 1};
x2 = x2r{1:h2, 1};
y2 = y2r{1:h2, 1};

P1 = polyfit(x1(:,1), y1(1,:), 6);
P2 = polyfit(x2(:,2), y2(2,:), 6);
allX = unique(X);
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd)

[Stelios Fanourakis]

@walter

I adjusted a bit the code

x1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','M2:M384');
y1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','A2:A384');
 
x2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','N2:N868');
y2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','B2:B868');
 
 
P1 = polyfit(x1{:,1}, y1{1,:}, 6);
P2 = polyfit(x2{:,2}, y2{2,:}, 6);
allX = unique(X);
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd)

But still I get the error

Error using polyfit (line 44)
X and Y vectors must be the same size.

Error in curvefitdif (line 11)
P1 = polyfit(x1{:,1}, y1{1,:}, 6);

Why are they not the same size now??

Since we are talking about x1, y1.

x2,y2 are different dataset

[Stelios Fanourakis]

@Walter

I don't know what to do.

Even if I assign the same number of cells to all variables x1,y1,x2,y2 still I get the error

Error using polyfit (line 44)
X and Y vectors must be the same size.

Error in curvefitdif (line 19)
P1 = polyfit(x1(:,1), y1(1,:), 6);

And it especially stucks at the specific line of code, where it shouldn't since I have made both x1 and y1 same size.

Is it something wrong with my Matlab?

[Walter Roberson]

x1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','M2:M384');
y1 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','A2:A384');
 
x2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','N2:N868');
y2 = readtable('ValidationTest.xls', 'Sheet',5, 'Range','B2:B868');

P1 = polyfit(x1{:,:}, y1{:,:}, 6);
P2 = polyfit(x2{:,:}, y2{:,:}, 6);

allX = unique([x1{:,}; x2{:,:}]);
Pd = polyval(P1, allX) - polyval(P2, allX);
plot(allX, Pd)

[Stelios Fanourakis]

I get the error

Error using unique (line 117)
Invalid input. Valid flags are 'rows', 'first', 'last', 'stable', 'sorted', 'legacy'.

Error in curvefitdif (line 14)
allX = unique(x1{:,:}, x2 {:,:});

[Stelios Fanourakis]

@Walter

I get NaN for P2. Why?

P1 =

   -0.0000    0.0000   -0.0019    0.0344   -0.2979    1.1208    4.9503


P2 =

   NaN   NaN   NaN   NaN   NaN   NaN   NaN

[Stelios Fanourakis]

I found the error of Nan but now I get the new error

Error using horzcat
Dimensions of arrays being concatenated are not consistent.

Error in curvefitdif (line 16)
allX = unique([x1{:,:}, x2{:,:}]);

[Walter Roberson]

You did not copy the assignment to allX from my code. My code clearly had

allX = unique([x1{:,:}; x2{:,:}]);

with [] in place and using semi-colon between the two arrays. You used

allX = unique(x1{:,:}, x2{:,:});

with no [] and using comma instead of semi-colon. Then you put back the [] but you used

allX = unique([x1{:,:}, x2{:,:}]);

with comma instead of semi-colon.

My guess about the P2 and NaN is that when you said to use N2:N868 and B2:B868 that you should have said N2:N867 and B2:B867 . The code I offered you with h1 and h2 and min() of size() would have taken care of that problem for you, but you decided that your version was more correct for your needs and you got burned because you did not understand what it was you were discarding.

[Stelios Fanourakis]

Yes thank you so much Walter. Now I see a result. I used to come up with various errors of brackets and columns that's why I experimented a bit but now it is sorted.

I am trying to plot everything now in one graph. How to plot P1, P2 and their difference together in one plot?

Actually, when I try to use those lines

Pd = polyval(P1, allX) - polyval(P2, allX)
 
P11 = polyval(P1, allX)
P22 = polyval(P2, allX)
 
hold on
plot(allX, Pd)
plot(allX, P11)
plot(allX, P22)

I get three plots. But one of the datasets is not correct

[Stelios Fanourakis]

Finally, I think it's great. Thank you very much Walter. I'll come back if I notice something not right.

[Stelios Fanourakis]

@Walter.

I just changed excel sheet to read a different set of datasets with the same row number and now I get again NaN as result for both P1 and P2. Why is this?

If I apply the same code for Sheet 5 and different columns but same number of rows I get numerical values. Why am I not getting them for sheet 6?

x1 = readtable('ValidationTest.xls', 'Sheet',6, 'Range','N2:N384 '); 
y1 = readtable('ValidationTest.xls', 'Sheet',6, 'Range','A2:A384 '); 
 
x2 = readtable('ValidationTest.xls', 'Sheet',6, 'Range','O2:O867'); 
y2 = readtable('ValidationTest.xls', 'Sheet',6, 'Range','B2:B867');
 
 
 
P1 = polyfit(x1{:,:}, y1{:,:}, 6 ) 
P2 = polyfit(x2{:,:}, y2{:,:}, 6)
 
allX = unique([x1{:,:}; x2{:,:}]);
 
Pd = polyval(P1, allX) - polyval(P2, allX)
 
P11 = polyval(P1, allX)
P22 = polyval(P2, allX)
 
hold on
plot(allX, Pd, 'b')
plot(allX, P11, 'g')
plot(allX, P22, 'r')

[Walter Roberson]

Dang, I posted complete code for this, but somehow it did not get saved!

[Walter Roberson]

%after the readtable
mask = isnan(x1{:,:}) | isnan(y1{:,:});
x1(mask,:) = [];
y1(mask,:) = [];
mask = isnan(x2{:,:}) | isnan(y2{:,:});
x2(mask,:) = [];
y2(mask,:) = [];

[Stelios Fanourakis]

I get the below error

P1 =

     0     0     0     0     0     0     0

Warning: Polynomial is not unique; degree >= number of data points. 
> In polyfit (line 74)
  In curvefitdif (line 19) 

P2 =

     0     0     0     0     0     0     0


Pd =

  0×1 empty double column vector


P11 =

  0×1 empty double column vector


P22 =

  0×1 empty double column vector

[Stelios Fanourakis]

I also get the error

Error using isnan
Too many input arguments.

Error in see (line 38)
mask = isnan(x1{:,:}) | isnan(y1{:,:});

Why?

Answer 2

[John D'Errico]

Polynomials are linear in the coefficients. So the difference of two polynomials is obtained by just subtracting the coefficients. If they are different order polynomials, just pad the low order polynomial with zeros as necessary to correspond to the higher order coefficients that are effectively zero.

[John D'Errico]

You claim to have two polynomials, each of degree 6. You want the difference between them. SUBTRACT THE COEFFICIENTS.

x = rand(1,50);
y1 = rand(1,50);
y2 = rand(1,50);

Fit.

p1 = polyfit(x,y1,6);
p2 = polyfit(x,y2,6);

Subtract.

pdiff = p1 - p2;

pdiff is a polynomial of order 6, that represents the difference between the two. You can evaluate it.

fplot(@(x) polyval(p1,x),[0 1])
hold on
fplot(@(x) polyval(p2,x),[0 1],'g')
fplot(@(x) polyval(pdiff,x),[0 1],'r')
axis([0 1 -1,1])

Is there a reason why you refuse to believe me? You can evaluate that difference function at any point.

[Stelios Fanourakis]

@John

pdiff gives me this

 =

   NaN   NaN   NaN   NaN   NaN   NaN   NaN

[John D'Errico]

So effectively, your question has nothing to do with your real problem?

Problem 1: Read the data in. xlsread should suffice.

Problem 2: Extract the data from a table. There is no need to create a table. See problem 1.

Problem 3: Apply polyfit. Um, use polyfit? Read the help for polyfit. Or read my comment.

Problem 4: Subtracting the polynomials. I told you how to do that.

You are making this wildly overly complex.

[Stelios Fanourakis]

See my most recent comment. I get NaN as I use

pdiff = P1 - P2

[Stelios Fanourakis]

It has to do with P2. I don;t know why it gives NaN. I am trying to find out at the moment

[Stelios Fanourakis]

Yours John was also good but I couldn't make the fplots work correctly. I was getting some results, but they were almost linearly, which in my case the data are curves.

[John D'Errico]

If pdiff gives a result of NaN, then one or both of the polynomials was already garbage. LOOK at the cofficients that you got from those two fits.

NaN can only ever result as a subtraction by things like inf - inf, or anything plus a NaN. So the problem lies in your original fits.

[Stelios Fanourakis]

Actually, the problem lies where I use readtable and import the columns from excel.

They are imported as NaNs

Answer 3

[Stelios Fanourakis]

Actually, both datasets are high degree polynomials e.g 6. They are almost identical and need to find their difference. How do I do that in Matlab?

[John D'Errico]

Did I not say just how to do that in my answer?