MATLAB Answers

Nonlinear data-fitting

48 views (last 30 days)
gjashta
gjashta on 28 Jan 2020
Commented: gjashta on 4 Feb 2020
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 Comments

Walter Roberson
Walter Roberson on 29 Jan 2020
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
gjashta on 29 Jan 2020

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)??

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 30 Jan 2020
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 Comments

Show 3 older comments
gjashta
gjashta on 3 Feb 2020
Walter Roberson I have constructed the series of equations and I followed all your steps. Then instead of the second 'myfun' above I used this AS YOU SUGGESTED:
y = variables(:,2);
x = variables(:,1);
t = length(x);
X=[t 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) + 2*X(:,2)
A* [b1; b2; b3; b4; b5; b6; b7] ==2*X(:,2);
B = A \ 2*X(:,2);
b = sym('b',[1 7]);
MF = myfun(b,X);
[A,B] = equationsToMatrix(MF,b);
fit_coefficients = double(A)\double(B);
But, the coefficients are very small in the end and the estimated y is very small compare with real y. I think I am doing something wrong.
Can you show me a plot of your model and the real data, how they look???
Walter Roberson
Walter Roberson on 3 Feb 2020
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
gjashta on 4 Feb 2020
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 ??

Sign in to comment.

More Answers (1)

Hiro
Hiro on 29 Jan 2020
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 Comments

Show 1 older comment
gjashta
gjashta on 29 Jan 2020
Hiro Yoshino , (1) yes.
(2)No, I am guessing 'a' in the code above. Maybe that's why I am not getting a good fit.
gjashta
gjashta on 29 Jan 2020
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
Hiro on 30 Jan 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!