How can I perform a curve fit using polyfit on an exponential decay model?

조회 수: 21(최근 30일)
Tyler
Tyler 2014년 10월 23일
편집: Matt Tearle 2014년 10월 23일
Hi, all.
I have been working with using polyfit to find curve fits for various sets of data recently. I stumbled across a problem in doing scalar addition, however. For example,
X = 1:100;
Y(1,:) = 3*exp(-0.5.*(X));
Y(2,:) = Y(1,:);
C = log(Y(2,:));
D = polyfit(X,C,1);
Z = exp(D(2)).*exp(D(1).*X);
plot(X,Z)
grid on
hold on
plot(X,Y(1,:),['g','o'])
D(1)
D(2)
will yield a perfect curve fit, yet, if a scalar of, say, 10 is added to the Y data, as in:
X = 1:100;
Y(1,:) = 3*exp(-0.5.*(X))+10;
Y(2,:) = Y(1,:)-10;
C = log(Y(2,:));
D = polyfit(X,C,1);
Z = exp(D(2)).*exp(D(1).*X)+10;
plot(X,Z)
grid on
hold on
plot(X,Y(1,:),['g','o'])
D(1)
D(2)
I cannot produce any results. Does anyone have any suggestions as to how to remedy this?
Thanks,
Tyler E.

답변(1개)

Matt Tearle
Matt Tearle 2014년 10월 23일
편집: Matt Tearle 2014년 10월 23일
The problem is that when X is big, 3*exp(-0.5.*(X)) is very small. Because of how floating-point arithmetic works, that's OK as long as it's bigger than realmin because it will be stored as 1.23e-42 or whatever, and taking the log of that is fine. But when you add 10 you get 10+(1.23e-42) which comes out as 10 because of finite precision arithmetic. So then, you subtract 10 and get 0 (the 1.23e-42 is lost as a casualty of finite-precision arithmetic). Taking the log of 0... yeah, you see the problem.
So, the remedy...
One way would be just to remove the 0s first (or, equivalently, the Inf values from the transformed data):
idx = isfinite(C);
D = polyfit(X(idx),C(idx),1);
But do you need to do this using polyfit? There are actually reasons why a linear fit to the transformed data can be a bad idea (if there is any noise in your data). Another approach would be to treat the problem as a nonlinear curve fit. You're trying to fit c1*exp(c2*x)+c3 to the data, so:
X = 1:100;
Y = 3*exp(-0.5.*(X))+10;
f = @(c,x) c(1)*exp(c(2)*x)+c(3);
D = lsqcurvefit(f,zeros(3,1),X,Y)
Z = f(D,X);
plot(X,Z)
grid on
hold on
plot(X,Y,['g','o'])
Note, this does require Optimization TB.

Community Treasure Hunt

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

Start Hunting!

Translated by