How to use polyfit to get the difference out of two functions?
15 次查看(过去 30 天)
显示 更早的评论
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?
2 个评论
Walter Roberson
2019-4-21
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
2019-4-21
@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?
采纳的回答
Walter Roberson
2019-4-21
P1 = polyfit((1:length(X)).', X(:), 6);
Likewise for Y producing P2
Pd = polyval(P1, 1:length(X)) - polyval(P2, 1:length(X))
23 个评论
Stelios Fanourakis
2019-4-21
@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?
Walter Roberson
2019-4-21
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
2019-4-21
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
2019-4-21
%"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
2019-4-21
编辑:Stelios Fanourakis
2019-4-21
@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
2019-4-21
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
2019-4-21
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
2019-4-21
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
2019-4-21
编辑:Stelios Fanourakis
2019-4-21
@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
2019-4-21
@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
2019-4-22
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
2019-4-22
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
2019-4-22
编辑:Stelios Fanourakis
2019-4-22
@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
2019-4-22
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
2019-4-22
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
2019-4-22
编辑:Stelios Fanourakis
2019-4-22
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
2019-4-22
Finally, I think it's great. Thank you very much Walter. I'll come back if I notice something not right.
Stelios Fanourakis
2019-4-22
@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
2019-4-22
%after the readtable
mask = isnan(x1{:,:}) | isnan(y1{:,:});
x1(mask,:) = [];
y1(mask,:) = [];
mask = isnan(x2{:,:}) | isnan(y2{:,:});
x2(mask,:) = [];
y2(mask,:) = [];
Stelios Fanourakis
2019-4-23
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
2019-4-24
I also get the error
Error using isnan
Too many input arguments.
Error in see (line 38)
mask = isnan(x1{:,:}) | isnan(y1{:,:});
Why?
更多回答(2 个)
John D'Errico
2019-4-21
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.
8 个评论
John D'Errico
2019-4-21
编辑:John D'Errico
2019-4-21
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
2019-4-21
编辑:Stelios Fanourakis
2019-4-22
@John
pdiff gives me this
=
NaN NaN NaN NaN NaN NaN NaN
John D'Errico
2019-4-22
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
2019-4-22
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
2019-4-22
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
2019-4-23
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
2019-4-23
Actually, the problem lies where I use readtable and import the columns from excel.
They are imported as NaNs
Stelios Fanourakis
2019-4-21
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?
1 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)