Nonlinear data-fitting

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

[Walter Roberson]

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]

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]

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.

[Walter Roberson]

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]

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]

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]

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]

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]

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]

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.

[Walter Roberson]

The first section of it looks like it should just be polyfit with degree 2.

[gjashta]

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]

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]

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.