Nonlinear data-fitting

2020 1 28
2 个回答
回答已采纳
更新时间：2020 2 4
2 次查看（30 天）

0 个投票

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
Walter Roberson 2020-1-29
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 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)??

请先登录，再进行评论。

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

 采纳的回答

Walter Roberson
Walter Roberson 2020-1-30

0 个投票

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
Walter Roberson 2020-1-30
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.
gjashta
gjashta 2020-1-30
编辑：gjashta 2020-1-31
Walter Roberson, I am using the first fit to remove some deterministic patterns and then the data after the transformation in the first fit I am fitting in the second fit.
Are you suggesting to use InitGuess=[0 0 0 0 0 0 0] ???
Walter Roberson
Walter Roberson 2020-2-1
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);
gjashta
gjashta 2020-2-3
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 2020-2-3
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 2020-2-4
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 ??

请先登录，再进行评论。

更多回答（1 个）

Hiro Yoshino
Hiro Yoshino 2020-1-29

0 个投票

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更早的评论

Walter Roberson
Walter Roberson 2020-1-29
The first section of it looks like it should just be polyfit with degree 2.
gjashta
gjashta 2020-1-29
编辑：gjashta 2020-1-29
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 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
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by