Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
278 views
in Technique[技术] by (71.8m points)

numerical methods - MATLAB - Fit exponential curve WITHOUT toolbox

I want to fit a decaying exponential to the plotted data. I do NOT have the Curve Fitting or Optimization Toolboxes.

x = [0    0.0036    0.0071    0.0107    0.0143    0.0178    0.0214    0.0250    0.0285    0.0321    0.0357    0.0392    0.0428    0.0464    0.0464];
y = [1.3985    1.3310    1.2741    1.2175    1.1694    1.1213    1.0804    1.0395    1.0043    0.9691    0.9385    0.9080    0.8809    0.7856    0.7856];

figure()
plot(x,y,'*')

How could I achieve this in MATLAB?

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)

Assuming that you have Gaussian distributed errors between the input and output points, and assuming that the errors are additive, you can solve this by classic least squares. It boils down to having an overdetermined linear system of equations where each constraint defines one input-output observation. Solving this overdetermined linear system with the least amount of residual error is the solution you're looking for.

Caveat

Jubobs made a very interesting point in his comment to me below. The parameters that minimize this residual error don't, in general, minimize the residual error of the original problem. This linearization step allows us to solve the parameters in an easier way, but it isn't the equivalent problem. However, it is usually accepted in practice as the solution is good enough.


To get this into a linear system, we need to do some clever rearranging. Because you want to fit a series of points using an exponential model, the relationship between the input x and output y is:

To get this to be "linear", we can take the natural logarithm of both sides:

By using the fact that ln(ab) = ln(a) + ln(b), we have:

Also knowing that , this simplifies to:

As you can see, the above equation is now "linear" with respect to log-space. Given a bunch of x and y values, so (x_1, x_2, ..., x_n) and (y_1, y_2, ..., y_n), we can concatenate a bunch of equations together in a linear system:

If we let ln(A) = A' for ease of notation, and rearranging this so that it's in matrix form, we get:

Therefore, we simply need to solve for A' and b, and you can do that by the pseudoinverse. Specifically, the above problem is of the form:

Therefore, we need to solve for X, and so:

M^{+} is the pseudoinverse of the matrix. Once you're done, simply take the exp operator on A' to get the original A. MATLAB has very efficient linear system solvers and least-squares solvers. Specifically, you can use the or ldivide operator. All you have to do is create the M matrix from the x values, create a vector of y values and solve your system. It's really simply:

x = ...; %// Define numbers here - either row or column vectors
y = ...;
M = [ones(numel(x),1), x(:)]; %// Ensure x is a column vector
lny = log(y(:)); %// Ensure y is a column vector and take ln

X = M  lny; %// Solve for parameters
A = exp(X(1)); %// Solve for A
b = X(2); %// Get b

Therefore, using your x and y values, this is what I get:

A =

    1.3882

b = 

   -11.508

If we plot the above points as well as an exponential curve that fits the line, we can do:

xval = linspace(min(x), max(x));
yval = A*exp(b*xval);
plot(x,y,'r.',xval,yval,'b');

The first line of code defines a bunch of x values that span between the smallest and largest x value for our data set. For the next line, we then take the x values and run them through our exponential model. Finally, we plot both the original data points, as well as the exponential curve with the parameters found through the above procedure together. The points are in red while the line is in blue.

We get:

enter image description here

I think that looks pretty good! For those people who have noticed, the above plot looks a bit different than a normal plot and figure window that is produced by MATLAB. That plot was generated in Octave as I don't have MATLAB on the computer that I'm currently working on. However, the above code should still work in MATLAB.


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...