Nonlinear data-fitting

Question

gjashta 2020-1-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/502381-nonlinear-data-fitting

评论： gjashta 2020-2-4

采纳的回答： Walter Roberson

Variables.mat

在 MATLAB Online 中打开

I have the data below, where y=variables(:,2), x=variables(:,1) and t=length(x).

Do you have any idea how to get e better fitting in mdl2?

myfun=@(b,x) b(1)+b(2)*t+b(3)*t.^2;   
InitGuess=[8.2075e+05 1 1];
mdl1=fitnlm(t,y,myfun,InitGuess);
YY=feval(mdl1,t);
figure,plot(t,1-YY./y'-b',t,0*t,'-r'), title('Relativ error ')
figure, plot(t,y,'-g',t,YY,'-b'), title('Data and model'), legend('Data','Model')
%%
Z=y-YY;
figure,plot(Z)
X=[t x];
f1=@(a,X) sin(a*X);
f2=@(a,X) cos(a*X);  
myfun=@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1); 
InitGuess=[-7.6342e-08 1 1 1 1 1 1];
mdl2=fitnlm(X,Z,myfun,InitGuess)
ZZ=feval(mdl2,X);

2 个评论
显示无隐藏无

Walter Roberson 2020-1-29

在 MATLAB Online 中打开

myfun=@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/31.5,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1);

Has subexpression

b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/31.5,X(:,1))

The f2 parts are the same so that is (b(2)+b(3)) times the f2 part. You would then combine b(2)+b(3) in to a single parameter.

Or is there a mistake in the formula?

gjashta 2020-1-29

编辑：gjashta 2020-1-29

Yeah, Sorry Walter Roberson. There is a mistake. I just corrected it. In b(3) and b(4) I am trying to use the same “a” for the trigonometric expression as in b(1) and b(2). Do you have any idea how to have a good model with less parameters that fits better the data Z ( the residuals from the first section)??

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2020-1-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/502381-nonlinear-data-fitting#answer_412751

The second system of equations has no constant terms, so the model fits perfectly if you set all of the parameters to 0.

This is not just an accident: you can do a calculus analysis of the sum of squares of the entries and show that all zeroes is the only critical point.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2020-1-30

在 MATLAB Online 中打开

The posted code has a lot of mistakes. Code with fewer mistakes would be

filestruct = load('Variables.mat');
variables = filestruct.Data1Rm;
y = variables(:,2);
x = variables(:,1);
t = length(x);
[sx, xind] = sort(x);
sy = y(xind);
mdl1 = polyfit(sx, sy, 2);
YY = polyval(mdl1, sx);
figure
plot(sx, sy, 'b*', sx, YY, '-r')
legend({'actual', 'predicted by degree 2 polynomial'});
Z=sy-YY;
figure,plot(sx, Z)
title('actual minus predicted')
X=[(1:t).', x];
f1=@(a,X) sin(a*X);
f2=@(a,X) cos(a*X);  
myfun=@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1); 
InitGuess=[-7.6342e-08 1 1 1 1 1 1];
mdl2=fitnlm(X,Z,myfun,InitGuess)
ZZ=feval(mdl2,X);

However as I indicated above, your myfun has no constant addative term that is independent of the b values, so it is fit perfectly by b being all 0.

Your original code tried to fit sample number as the independent variable in the quadratic function, instead of using x. That makes no sense to me, and using sample number as the first column of X in your second fit makes no sense to me either.

The fit reached by polyfit() will, by the way, be quite different than the fit reached by your original quadratic fit, because your starting point was not near to the actual solution. The actual solution, here found with polyfit, does not require any searching: it is completely deterministic solved by constructing a vandermode matrix.

The first fit you do, to the quadratic equation, plays no part in the second fit, and it is not obvious why you bother doing that first fit.

Walter Roberson 2020-2-1

在 MATLAB Online 中打开

You construct a whole series of equations of the form

M = [
b1 * t11 + b2 * t12 + b3 * t13 + b4 * t14 + b5 * t15 + b6 * t16 + b7 * t17
b1 * t21 + b2 * t22 + b3 * t23 + b4 * t24 + b5 * t25 + b6 * t26 + b7 * t27
b1 * t31 + b2 * t32 + b3 * t33 + b4 * t34 + b5 * t35 + b6 * t36 + b7 * t37
...]

and so on, where the t* variables are effectively constant for this purpose (you construct them from values that are read in from the file.)

Each row of that is to equal 0, so

M == zeros(size(M,1), 1)

will ideally be true.

This is a linear system of equations that could be written

A = [
  t11, t12, t13, t14, t15, t16, t17
  t21, t22, t23, t24, t25, t26, t27
  t31, t32, t33, t34, t35, t36, t37
  ...]
A * [b1; b2; b3; b4; b5; b6; b7] == zeros(size(M,1), 1)

which then could potentially be solved as

B = A \ zeros(size(M,1), 1)

However, the special case of A*x = 0 for known matrix A and unknown x, is finding the null space of A. It is obvious that if you set x all 0 then A*x = 0 will be true, so it is not even worth doing the A \ zeros(size(M,1), 1) if that is the only solution you are looking for: b(1) = b(2) = b(3) = b(4) = b(5) = b(6) = b(7) = 0 will always be a solution of that system of equations, no matter what the input values are.

The question then becomes whether there are other b solutions:

b = sym('b',[1 7]);
MF = myfun(b,X);
[A,B] = equationsToMatrix(MF,b);
null(A)

ans =
 
Empty sym: 7-by-0

No, there are no other solutions.

rank(A) is 7, same as the number of columns, rank(A) would have to be less than the number of columns in order for the null space to be non-empty.

Therefore, there is only a single solution for your myfun function coming out to all 0, and that solution is b(1) = b(2) = b(3) = b(4) = b(5) = b(6) = b(7) = 0 . So, there is no point in running the fitting: the unique answer will be 0 (unless the input file has fewer than 7 rows, or has a lot of rows that are linear combinations of each other.)

There would be different results if instead of

myfun=@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1); 

you had

myfun=@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1) + SOMETHING_DERIVED_FROM_X; 

where SOMETHING_DERIVED_FROM_X is not being multiplied by one of the b entries. In such a case, instead of

A * [b1; b2; b3; b4; b5; b6; b7] == zeros(size(M,1), 1)

you would be working with

A * [b1; b2; b3; b4; b5; b6; b7] == column vector of SOMETHING_DERIVED_FROM_X

which you would then work with as

B = A \ column vector of SOMETHING_DERIVED_FROM_X

The result would not be a vector of coefficients that made the system exactly equal to 0: it will be a vector of coefficients that minimizes the least-squared error of the system -- the best fit. Skipping the fitting call, it would go like

b = sym('b',[1 7]);
MF = myfun(b,X);
[A,B] = equationsToMatrix(MF,b);
fit_coefficients = double(A)\double(B);

Walter Roberson 2020-2-3

在 MATLAB Online 中打开

In the function handle you have + 2*X(:,2) . In the line

A* [b1; b2; b3; b4; b5; b6; b7] ==2*X(:,2);

you have 2*X(:,2) on the right hand side: if you were to bring it to the left hand side to have an equation equal to 0, then it would correspond to - 2*X(:,2) rather than to + 2*X(:,2)

B = A \ 2*X(:,2);

That corresponds to - 2*X(:,2) not to + 2*X(:,2)

filestruct = load('Variables.mat');
variables = filestruct.Data1Rm;
y = variables(:,2);
x = variables(:,1);
t = length(x);
X=[(1:t).', x];
myfun2 =@(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1) + 2*X(:,2)
b = sym('b',[1 7]);
MF = myfun2(b,X);
[A,B] = equationsToMatrix(MF,b);
fit_coefficients = double(A)\double(B);
sum(myfun2(fit_coefficients,X).^2)

ans =
      1.24829349642173e-26

That is a pretty good fit.

Note that your fitting function makes no reference to y. Your myfun functions are predicting x not y.

Predicting y might be something like

filestruct = load('Variables.mat');
variables = filestruct.Data1Rm;
y = variables(:,2);
x = variables(:,1);
t = length(x);
X = [(1:t).', x, y];
myfun2 = @(b,X)  b(1)*f1(2*pi/30.5,X(:,1)).*X(:,2) + b(2)*f2(2*pi/31.5,X(:,1)) + ...
+ b(3)*f2(pi/2.3,X(:,1))+ b(4)*f1(pi/2.3,X(:,2)).*X(:,1)+b(5)*X(:,2).*X(:,1) +b(6)*X(:,2)+b(7)*X(:,1) - X(:,3)
b = sym('b',[1 7]);
MF = myfun2(b, X);
[A,B] = equationsToMatrix(MF, b);
fit_coefficients = double(A)\double(B);

fit_coefficients =
         -1231.92429160698
          20370.1016972963
         -16.9180535390798
          873.045987898082
           -465.9789817208
          54259.3096584326
          6265.66604056945

  ypred = myfun2(fit_coefficients,X) + X(:,3);
  scatter(x, y, 'b');
  hold on
  scatter(x, ypred, 'r');
  hold off
  legend({'original', 'predicted'})

gjashta 2020-2-4

residuals.fig

I have attached the plot of (y-ypred) and it's quite big.I also tried to use the residuals from myfun1 as y but still I got e big variance between the y and predicted y.

Walter Roberson when you are Predicting y you are not taking into account the myfun1, right?

What is your R-squared of y explained by x and t ??

请先登录，再进行评论。

Answer 2

Hiro Yoshino 2020-1-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/502381-nonlinear-data-fitting#answer_412569

I just wonder if this is a linear model?

(1) Is b a coefficient vector?

(2) Do you need to estimate "a" too?

if (1) yes, (2) no, then this is a linear model and you can solve this analytically.

4 个评论
显示 2更早的评论隐藏 2更早的评论

gjashta 2020-1-29

Thanks Walter Roberson

I can do that for the first section, but in the second section I want to use both variables (t,x) to explain y.

Hiro Yoshino 2020-1-30

sorry for late.

Well, these are linear models, i.e., you don't need to run optimization to obtain parameters.

The solutions are analytically calculated.

As long as the parameters $$\mathbf{b}$$ are linear with respect to the given data $$X$$, the problem is called "linear problem". The solution can be given by

.

in matlab, fitlm is the one you should apply to this problem.

请先登录，再进行评论。