MATLAB: Function fitting on a set of data points

curve fittingMATLAB

Hello,

I am trying to fit a function on a set of data point using some parameters to be optimized and the function lsqcurvefit.

The function is essentially:

y = a+b*cos(x)+c*cos(2*x)+d*sin(x)-e*sin(2*x);

So I first defined the parameters in terms on one variable, as following:

x(1) = a

x(2) = b

x(3) = c

x(4) = d

x(5) = e

And here goes the code:

DATA = load('data_f.txt');
t = DATA(:,1);
y = DATA(:,2);
plot(t,y,'ro')
F = @(x,xdata) +x(1) - x(2)*cos(xdata*(pi/180)) + x(3)*cos(2*xdata*(pi/180)) - x(4)*sin(xdata*(pi/180)) - x(5)*sin(2*xdata*(pi/180));
x0 = [1 1 1 1 1];
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)
hold on 
plot(t,F(t,y))
hold off

Problem is that the curve looks horrible but the parameters look fine compared to the fitting done in other programs. Is there something I am not seeing?

Thank you very much in advance!

Stay safe,

Alex

Best Answer

What did you do wrong? It looks like Jon may be correct, in that you did not use the newly estimated set of parameters from the fit.

xy = [0.0	6.08
15	5.3
30	4.28
45	2.89
60	1.31
75	-0.16
90	-0.47
105	0.23
120	1.68
135	3.55
150	5.32
165	6.51
180	6.86
195	6.29
210	4.92
225	3.04
240	1.19
255	0
270	-0.09
285	0.75
300	2.21
315	3.89
330	5.29
345	6.09
360	6.08];
x = xy(:,1);
y = xy(:,2);

FIRST, PLOT YOUR DATA. LOOK AT IT. Does this simple trig model make sense?

plot(x,y,'ro')

And, yes, the result does seem to be a very simple sinusoidally varying curve.

However, there is no need to use lsqcurvefit to fit this, since it is a linear model in the parameters. Using a nonlinear estimation tool to fit a linear model is like using a Mack truck to carry a pea to Boston. A bit of over kill. But, yes, you could use lsqcurvefit. I won't bother to do so here.

A = [ones(size(x)),sind(x),sind(2*x),cosd(x),cosd(2*x)];
coefs = A\y
coefs = 5×1
    3.1957
   -0.1826
   -0.4175
   -0.2026
    3.3489
plot(x,y,'ro')
hold on
fun = @(u,coefs) coefs(1) + coefs(2)*sind(u) + coefs(3)*sind(2*u) + coefs(4)*cosd(u) + coefs(5)*cosd(2*u);
plot(x,fun(x,coefs),'b-')
hold off

While there is a little lack of fit, it is not terrible. Look at the residual errors.

plot(x,y - fun(x,coefs),'o-')

Lack of fit is usually seen as patterning in the residuals. Here we see exactly that. So your model is probably not the correct model for this data, merely a good approximation to that data. By way of comparison, how much of the information in your original signal has been explained by the model?

[var(y),var(y - fun(x,coefs))]
ans = 1×2
    6.1724    0.0445

As a percentage of the total variance...

100*(1 - var(y - fun(x,coefs))/var(y))
ans = 99.2797

So around 99% of the total variance in y has been explained by your model. So pretty good, but not perfect.

Related Solutions

MATLAB: Curve fitting and optimisation code problems for f(x)= a.*exp(-b./(x.^p)) to dataset

A serious problem with this is what I would suggest looks like almost two functions merged into one. Worse, the two "curves" at the bottom end are not consistent in either case with what we see at the top end, so schizophrenic data. And worse yet, you have a least two large outliers. Large residual nonlinear least square problems are difficult to solve, because one data point can really swings things around. It has an effectively huge residual. What else? Oh yeah, your data varies by an order of magnitude. The result of that is the points at one end of the curve impart effectively more influence on the result than they should.

First what happens to your model if we log y? This actually solves one of the issues I mention above, since we no longer have data that varies by powers of 10 when the fit is performed. Effectively, I am suggesting we think of your data in a semi-logy form.

semilogy(x,y,'o')
grid on
xlabel X
ylabel Y

The strange behavior of your data should be clear here. Thus we see the outliers (be careful not to miss the outlier on the very bottom.) We see data that on the right half seems to be nearly a perfectly straight line. The left half we see two pieces that clearly represent two completely different and inconsistent processes. They are also inconsistent with the data on the right. Two different things are happening there, neither of which is consistent with the right half of your data.

What will happen when we try to estimate a model, is your curve will not represent all three of those behaviors.

Before I get into a fit, first, log your model. That is, when we form the semilogy plot, we are changing your model from this:

y = a*exp(-b/(x^p))

into this:

log(y) = log(a) - b*x^(-p)

Dividing by x^p is just a multiply by x^(-p). And log(a) is still just a constant that we don't know. Effectively, we can think of the model in a logged form as almost a polynomial. Most seem to call that a power-law model. Again though, these curves really do seem to be purely linear when we look at the semi-logy plot. That implies p==-1 is about as close as my eye wants to fit this. Remember, you have that -p in there, so two negatives do make a positive in this case. As such, the proper model for your curve prbably has p as very close to -1.

Time to do a first cut at a model. I'll just use the basic fitting tools from the plot figure itself to create that model.

plot(x,log(y),'o')
grid on

On the right half, again, that data really is nearly perfectly linear, implying p==-1. But on the left half, which of those curves should it follow? A curve fit cannot follow both curves. And if it chases one, then it misses the other by much more. And worse, if you do go chase one of those lower personalities, then the curve will miss the upper part by more yet.

The funny thing is, those two nasty looking outlier points almost seem to be the least of your problems. Even so, I would still delete them from the analysis before proceeding further. They are still influencing your estimates.

What would I do?

Delete the two outliers.
Determine why it is your data is so schizophrenic. What happens that drives the two behaviors? Why is there a split at roughly x~150? Choose the process you wish to accept as truth. Only when you excise the split personalities in your data can you build a model that works.
Finally, either split the problem into two pieces, performing a broken-stick regression with the break between the two halves of your curve, or try to build once global model that won't fit either half very well.

So first, in your data set (ignoring the NaNs) we see

Right half of the curve: points # 10:47

The two outliers: # 1, 49

Left curve (lower branch: # 2:9

Left curve, upper bransh: # 48:60

So if I had to guess, you are doing something strange to generate the data, that creates different behaviors that are inconsistent with each other. A model is not going to bring it all together in some magic way.

I might split the data like this:

ind1 = 10:47;
ind2 = 2:9;
ind3 = 48:60;
ind4 = [1 49];

If we just consider the piece at the right, we would see:

p1 = polyfit(x(ind1),log(y(ind1)),1)
p1 =
        0.0185627813572888         -27.2729761320311

Nothing more sophisticated than polyfit was needed. So a very good model for that curve segment will be simply

y = exp(-27.2729761320311) * exp(0.0185627813572888 * x)

Note the fact that x in that model did not even require a negative exponent.

Plot the data now, with this model. As you will see, the fit really is very good. I'll even view it in an unlogged form, to show I was sucessful in posing a reasonable model for the entire process, even though the model I would pose is far simpler than what you have intended.

plot(x(ind1),y(ind1),'bo',x(ind2),y(ind2),'go',x(ind3),y(ind3),'ko',x(ind4),y(ind4),'r*')
grid on
hold on
xhat = 40:230;
yhat = exp(p1(2))*exp(p1(1)*xhat);
plot(xhat,yhat,'b-')

In fact, the model fits splendidly well for the upper half. There simply is no good way to resolve the split personality in the lower half though. At least not so until you understand why the split exists, and how to properly deal with it. My suggestion at this point would be to get better data.

Finally, can I reconcile your claim that p should be on the order of 0.25 to 0.5? That seems to be wildly inconsistent with the data you have provided.

MATLAB: Problems using the curve fitting app

Without even looking at your data, just your question, you need to understand how any such optimization tool works. Pretend you are a blind person, tasked with finding the location of lowest altitude on the surface of the earth. Yes, I know, it happens to be deep under water in the Pacific ocean. I'll give you scuba gear.

ll you are allowed to do is to use a cane, to infer the local clope of the surface. Then you will walk downhill. I told you to use the scuba gear! Things will get wet.

But now, suppose you were started out in some unfortunate place. Perhaps the vicinity of the Dead Sea? Thus a local depression where as you walk downhill, while you may get wet, you will also get stuck in the wrong place?

The point is, you need to use intelligent starting values for any such solver. Start it in a bad (wildly wrong) place, and you should expect to get complete crapola as a result.

DISCLAIMER: No blind individuals werre harmed in this thought experiment.

Having said all that, now we can look at your data.

plot(x,y)

Sigh. Can I suggest you learn to work in better units? Numbers with hugely varying scales will always cause you numerical problems on a computer. So instead of working with numbers that are on the order of 10^-11 seconds, use picoseconds or at least nanoseconds. Since a picosecond is 10^-12 seconds, that seems perfect here.

xp = x/1e-12;

If nanoseconds would make you happier, then use them instead.

Next, consider shifting/translating the model. Since your data starts at roughly 18, then translate the model to reflect that. Otherwise, you will be using exponentials that start at 18. That causes you to have strange values appear in the coefficient a.

Sorry though. I won't use the app here, as I cannot as easily show and explain the results.

ft = fittype('a*cos(2*pi*f*(x-x0)).*exp(-(x-x0)/tp)','indep','x','problem','x0')
ft
ft = 
     General model:
     ft(a,f,tp,x0,x) = a*cos(2*pi*f*(x-x0)).*exp(-(x-x0)/tp)

Here, when I plot the data, I can infer intelligent start points for the parameters.

plot(xp,y)

a is the value at x == x0. That should be on the order of 1.

f and tp are rate parameters. tp is an exponential decay. 10 should be reasonable.

The most sensitive one is f, which as you have written it, should be roughly the inverse of the period, so a frequency. We can think of f as the number of oscillations per picosecond. Since the period of oscillation seems to be roughly 2, 0.5 is an intelligent start point. To get a decent estimate of the period of oscillation, I'll use the intersections tool written by Doug Schwarz, as found on the file exchange.

xpint = intersections(xp,y,[0 40],[0 0])
xpint =
          18.6739787975366
          19.7766155889863
          20.8497137648191
          21.9085954944223
          22.9711077650663
          24.0312152815146
          25.0899375680738
          26.1465450240289
          27.2047855420838
          28.2608584039546
          29.3180979049269
          30.3759822488665
          31.4383534867009
          32.4978593885987
          33.5605930233276
           34.622453287097
          35.6820730748939
          36.7440919564783
          37.8046527675794
          
diff(xpint(1:2:end))
ans =
          2.17573496728251
          2.12139400024718
          2.11882980300753
          2.11484797400998
          2.11331236284309
          2.12025558177397
          2.12223953662676
          2.12148005156622
          2.12257969268551

The period seems to be vaguely constant, though perhaps the first period is off from the remainder.

1/median(diff(xpint(1:2:end)))
ans =
         0.471388153206563

Time to fit a model.

mdl = fit(xp,y,ft,'start',[1,.5,10],'problem',18)
mdl = 
     General model:
     mdl(x) = a*cos(2*pi*f*(x-x0))*exp(-(x-x0)/tp)
     Coefficients (with 95% confidence bounds):

       a =      0.9925  (0.9156, 1.07)
       f =      0.4445  (0.4405, 0.4485)
       tp =       3.367  (2.971, 3.764)
     Problem parameters:
       x0 =          18

That seems to fit. At least, when I plot it, things seem to be somewhat reaonable.

plot(mdl)
hold on
plot(xp,y)

Now getting reasonable. This model suffers from a flaw, in that it assumes a fixed period of oscillation in that decay. Instead, your oscillations seem to be also slowing down. You could repair that, using a time varying coefficient f. Perhaps f(xp) might take the form f0 - f1*xp.

ft = fittype('a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)','indep','x','problem','x0')
ft = 
     General model:
     ft(a,f0,f1,tp,x0,x) = a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)
mdl = fit(xp,y,ft,'start',[1,.5,.01,10],'problem',18)
mdl = 
     General model:
     mdl(x) = a*cos(2*pi*(f0 - f1*(x - x0)).*(x-x0)).*exp(-(x-x0)/tp)
     Coefficients (with 95% confidence bounds):
       a =      0.9983  (0.9463, 1.05)
       f0 =      0.4179  (0.4129, 0.4229)
       f1 =   -0.005473  (-0.006309, -0.004638)
       tp =       3.494  (3.218, 3.771)
     Problem parameters:
       x0 =          18
       
figure
plot(mdl)
hold on
plot(xp,y)

That was a bit of a hack in the model. It does better, but we could surely do better yet with some effort. The problem is your fit will be jerked around by that first period. Since the absolute numbers are the largest there, the fit will be most impacted by errors in that first period. That would suggest a simple weighted fit might be as good, where the weight on the latter part of the series is increased.

Best Answer

Related Solutions

MATLAB: Curve fitting and optimisation code problems for f(x)= a.*exp(-b./(x.^p)) to dataset

MATLAB: Problems using the curve fitting app

Related Question