MATLAB: Does Curve Fitting Toolbox produce a good exponential fit graph with a wrong equation

Curve Fitting Toolboxexponential fitextracting equationwrong confidence boundswrong equation

Hello everyone,

I have created a graph for my data points using second order exponent (a * exp(bx) + c * exp(dx)). You can see it on the first picture. It looks good to me, there is no problem with the graph. However, I need to get the equation of this graph and that is where the problem comes out.

On the second picture you can see the constants for my equation with unusually large 95% confidence bounds. The graph shows these boundaries to be very small, almost insignificant. Also if I plug in some value of x in the equation and the graph, for example x=25, the results are very different:

Equation: -1.297 * exp(-2.486 * 25) + 1.325 * exp(-2.48 * 25) = 2.47 * 10^-20

Graph: at x = 25, y = 0.183 (shown on the first picture)

You can find my variables by clicking on this link and downloading the excel file.

https://drive.google.com/file/d/1YVOYX0Mno68Ab5pr53ytZEC5_GuEKC6H/view?usp=sharing

The value on the graph is correct. My question is what is the problem and how can I solve it?

Best Answer

NEVER use 4 or 5 digit approximations to coefficients. In fact, even if you display the coefficients in the command window using format long, STILL don't copy them out, because you won't get the exact coefficient!

I don't have your data, so I can't show results for your data. I'll just make up some garbage data.

x = rand(10,1);
y = rand(10,1);
F = fit(x,y,'exp1')
F = 
     General model Exp1:
     poly(x) = a*exp(b*x)
     Coefficients (with 95% confidence bounds):
       a =      0.6954  (0.1274, 1.263)
       b =     -0.7068  (-2.599, 1.186)

Now, NEVER just copy those coefficients.

Instead, use coeffvalues to extract them.

C = coeffvalues(F)
ans =
      0.69543     -0.70679

That is not their true values though. Instead, the numbers are closer to this:

format long g
coeffvalues(F)
ans =
         0.695433813207908        -0.706794760999295

So C has the final set of coefficients used for the model.

Again however, try not to use 0.695433813207908 as the number. Instead, the model is:

y = C(1)*exp(C(2)*x)

Use it that way.

Ok, now to what I think is the more important part of your question. Why are those bounds so wide?

LOOK AT YOUR CURVE!!!!!!

Your curve basically looks like a negative exponential, with a minor tweak at the beginning. So the fit with one exponential term in there was probably not terribly bad. But sums of exponentials are moderately difficult things to estimate. Why do I say that?

Think of it as arising from the fact that exp(x) looks a lot like exp(2*x), or exp(k*x) in general. So the nonlinear solver can have some problems.

Now, look at the coefficients that it found! Don't just build a model and then assume it is correct! In fact, the model you estimated was probably complete crap! (Sorry, but true.)

LOOK AT THE COEFFICIENTS.

b = -2.486, with bounds [-56,51]
d = -2.48 with bounds of [-55,50]

Both exponentials are almost identical in shape, with the exception that the coefficients out front had opposite signs. This is a sign that the fit gave you lower error, but it is probably meaningless. Those coefficients have no value at all. Were the noise in the data to change by a small amount, you would get different numbers.

Again, those hugely wide bounds on the coefficients suggests that the model is terribly poorly estimated. (I say poorly because you really have no clue as to the actual values of those coefficients, if that was indeed the true model for this process.) Your data is simply insufficient to estimate a two term exponential fit. The noise is way too high, and you have way too little data for this model.

This is not a problem of the curve fitting toolbox, or of MATLAB at all, but in your choice of model for this data.

Is there a better choice for the model? Perhaps. But I don't know why you chose the model you did. Perhaps you have some reason for thinking it should be a sum of two exponentials. Honestly, I doubt that is true. Most of the time, people use an exponential model, then decide that it the curve is basically a negative exponential, then a sum of two of them MUST be better. Sorry, but that does not follow.

Related Solutions

MATLAB: Curve fitting error:nlinfit rank deficient

Whats wrong? Mainly really poor starting values. Not enough data is another good reason, and a model that is not easily fit on data that varies over multiple orders of magnitude.

Other problems? A better algorithm is often a good choice. You might also consider using a weighted nonlinear regression.

Here is my fit, using a better algorithm for this class of problem.

funlist = {@(c,xdata) 1 - exp(-xdata/c), 1};
[INLP,ILP] = fminspleas(funlist,1719,tr,si1)
INLP =
       427.15
ILP =
       1528.4
        175.6
plot(tr,si1,'o')
trev = 32:2050;
hold on
plot(trev,ILP(1)*(1-exp(-trev/INLP)) + ILP(2),'r-')
grid on

The fit here seems reasonable. Still not nearly sufficient data to do the fit well.

I used my own fminspleas (it lives on the file exchange, for user download.) It is better for this kind of problem because it is more robust to poor starting values.

So what is wrong here? Why did nlinfit fail? Look at the initial parameters you chose,

[0, 1719, 50] % start point

[1528.4, 427.8, 175.6] % estimated values

So when I used fminspleas, it has the ability to think about this as a ONE variable nonlinear optimization. p(1) and p(3) are intrinsically linear parameters, that can be estimated conditionally on the value of p(2). So fminspleas is able to succeed even when nlinfit gave up the ghost. But see how far away it had to move for the nonlinear parameter from your starting values. nlinfit has it much harder, since you gave it crap to start with in the estimation.

Could you have found more intelligent starting values? Probably so. Given this model...

fun1=@(p,x) p(1).*(1-exp(-x./p(2)))+p(3);

we see that at x == 0, fun1 will return exactly p(3).

As a lmit, as x--> inf, the function will approach p(1) + p(3).

If we look at the plot of the data, I might have guessed p(3) should have been roughly 100, or maybe even a little higher.. You tried 50. Regardless, that was not a terrible choice for p(3).

But if I had to guess a value for the asymptote as x approaches inf. I might have guessed roughly 1800, just from the plot.

So if we could have predicted

p(3) ~ 100

p(1) + p(3) ~ 1800

then your starting values were clearly way off. p(1) should be something like 1700, just looking at the data.

How about the starting value for p(2)? Lets see. An exponential rate parameter is easy to estimate. What is the point where a negative exponential has died off by say 90%? That is, can we find the value where

exp(-u) = 0.1

Easy.

-log(.1)
ans =
       2.3026

So, if I had to guess a start point for p(2), since the exponential has died by about 90% at x==1500, then we can find an very rough initial value for p(2) as

p2init = -1500/log(0.1)
p2init =
       651.44

Lets try nlinfit now.

fun1=@(p,x) p(1).*(1-exp(-x./p(2)))+p(3);
p0= [1700,-1500/log(0.1),100];
options = statset('Funvalcheck','off');
beta0 = p0;
[beta1,R,J,CovB,MSE] = nlinfit(tr,si1,fun1,beta0);%nonlinear fitting
beta1
beta1 =
       1528.4       427.15        175.6

Amazing. It worked, and did not even complain. ;-)

Better (more) data would be useful. And the model does not seem to fit terribly well, with not really the correct shape. But it is way better than the sharp stick in the eye that you got when you tried it. The biggest lessons you need to take from this are to use a better algorithm to fit exponential models, or if you are going to use an algorithm that cannot handle poor starting vlaues, then choose more intelligent starting values. Otherwise, expect to see crap out the end far too often.

What could you do to improve this mess? For example, if I were your own consultant, I would tell you to get better data first. The process has not yet gone far enough to know where it will level off. For example, let me use a constrained regression spline model to fit your data really well. I'll let it tell me where the data seeem to want to level off,

slm = slmengine(tr,si1,'knots',[0 50 100 200 400 800 1600 3200 5000],'increasing','on','concavedown','on','leftminvalue',0,'plot','on');

I used my SLM toolbox for the fit, providing only information that seemed logical from what I can see. That is, we have a curve that is a monotone increasing function, always non-negative, but always a negatively curved function.

As you can see, the spline model fits quite well compared to the negaive exponential. But it does not seem to know where to level the curve off on the top end. To me, that tells me you need a bit more data on the top end. A little more data down near zero would help too, in order to better resolve the behavior at both ends, as well as to pin down the rate parameter more accurately. It is the data near each end of the curve that provides the most information about p(1) and p(3) in your model, and the data in the middle that are telling the estimation the value of the rate parameter.

Is there a better choice of model? Perhaps. But I don't really know if you had a good reason to choose that model, or if you just thought it might work well. It is ok, but has too much lack of fit for my eyes. One option might be...

ft = fittype('a + b/(1+ exp(-x/c))')
ft = 
     General model:
     ft(a,b,c,x) = a + b/(1+ exp(-x/c))
mdl = fit(tr',si1',ft,'startpoint',[-1000 2900 600])
mdl = 
     General model:
     mdl(x) = a + b/(1+ exp(-x/c))
     Coefficients (with 95% confidence bounds):

       a =       -1204  (-1513, -893.7)
       b =        2871  (2477, 3264)
       c =       284.1  (187.6, 380.6)
plot(tr,si1,'o'),hold on
plot(mdl)

Sigh. Not really any better in terms of fitting the roll over. Well, that was only one of several classic forms you might have tried. Here is another:

ft = fittype('a + b*erf(x/c)')
ft = 
     General model:
     ft(a,b,c,x) = a + b*erf(x/c)
mdl = fit(tr',si1',ft,'startpoint',[0 1800 1000])
mdl = 
     General model:
     mdl(x) = a + b*erf(x/c)
     Coefficients (with 95% confidence bounds):
       a =       244.2  (73.16, 415.3)
       b =        1405  (1177, 1633)
       c =       664.1  (415.2, 913)
plot(mdl)
hold on
plot(tr,si1,'o')

Even more so, just too much lack of fit. What it tells me is your data is not rolling over sufficiently well to pin down the shape of that curve.

Still, I'd just use a spline. AND get better/more data where it counts.

MATLAB: How to plot the confidence bounds for an exponential fit

You mean something like this?

a = [1.842; 1.755; 1.929];
b = [-0.002508; -0.002784; -0.002232]; 
x = linspace(-10,50,100);
y_fit = a(1)*exp(b(1)*x);
y_cbl = a(2)*exp(b(2)*x);
y_cbh = a(3)*exp(b(3)*x);
plot(x,y_fit,'r',x,y_cbl,'b:',x,y_cbh,'b:')

Best Answer

Related Solutions

MATLAB: Curve fitting error:nlinfit rank deficient

MATLAB: How to plot the confidence bounds for an exponential fit

Related Question