1

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?

10
  • 2
    Asking people to "write the complete code for that" is not how this site works, as you probably know. You should try yourself and then ask something more specific if you need. My suggestion (I assume you want to fit the curve in a least-squares sense): Write down the equation of the sum-square-error as a function of the two parametes a and b of the exponential, y = a * exp(b * x). To minimize error, differentiate with respect to a and b and equate fo 0 Commented Mar 27, 2015 at 16:03
  • As you might be able to see from previous posts of mine, this isn't a request that I usually do. But, when the code is (probably) highly influenced by math (that I suck at) and note coding, I don't believe I would be able to adjust the code to my needs. If you have some suggestions I would be more than happy to hear them. Commented Mar 27, 2015 at 16:06
  • That's a classic least-squares problem with a well-defined solution. See the link I marked as a duplicate. Commented Mar 27, 2015 at 16:07
  • @Vihar Ok, removing my downvote. See my suggestion in my edited first comment Commented Mar 27, 2015 at 16:07
  • Thanks for the response, but I don't really know how to adapt that to the code... Commented Mar 27, 2015 at 16:08

2 Answers 2

8

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.

Sign up to request clarification or add additional context in comments.

5 Comments

May be worth mentioning that the solution to the modified problem does not, in general, coincide with that of the original problem... even though it's often good enough, in practice.
@Jubobs - Can you explain why it doesn't coincide with the original problem? No accusations are being made here. I'm just unsure as to why it doesn't. Some insight would be awesome and I'll certainly modify my post.
@rayreng In other words, the values of parameters A and b that minimise the residual of the modified problem (i.e. the square root of the sum of the (\ln y_i - \ln A - b x_i)^2, in LaTeX notation) do not, in general minimise the residual of the original problem (i.e. the square root of the sum of the (y_i - A e^{b x_i))^2). The linearization step is just a trick that allows us to solve a related (but not equivalent) problem with linear least squares.
@Jubobs - Ah. Now that makes perfect sense. Thanks. I'll make a modification.
@Jubobs - Merci beaucoup mon ami :)
1

And for the more stupid like me, the entire code (thanks to rayryeng) is:

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

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

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

4 Comments

Oh c'mon the math isn't that bad :)... but I'm glad you posted a tl;dr answer!
Hah...I can argue about the "fun" in math :)
Haha fair enough. Did the math make sense though? What I wrote that is.
@rayryeng: How can I modify the code above, in order to have the function as: a+b* exp(c*xval) ? I am missing the first constant.Thanks!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.