MATLAB - Fit exponential curve WITHOUT toolbox

Question

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?

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 — Luis Mendo
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. — Trenera
That's a classic least-squares problem with a well-defined solution. See the link I marked as a duplicate. — rayryeng
@Vihar Ok, removing my downvote. See my suggestion in my edited first comment — Luis Mendo
Thanks for the response, but I don't really know how to adapt that to the code... — Trenera

rayryeng rayryeng · Accepted Answer · 2015-03-27T16:37:04

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:

$y = Ae^{bx}$

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

$\begin{align*}\ln y = \ln(Ae^{bx})\end{align*}$

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

$\begin{align*}\ln y = \ln(A) + \ln(e^{bx})\end{align*}$

Also knowing that $\begin{align*}\ln(e^z) = z\ln(e) = z\end{align*}$ , this simplifies to:

$\begin{align*}\ln y = \ln A + bx\end{align*}$

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:

$\begin{align*}\ln y_1 &= \ln A + bx_1\\\ln y_2 &= \ln A + bx_2\\\ln y_3 &= \ln A + bx_3\\&\ldots\\&\ldots\\\ln y_n &= \ln A + bx_n\\\end{align*}$

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

$\begin{bmatrix}\ln y_1 \\ \ln y_2 \\ \ln y_3 \\ \cdots \\ \cdots \\ \ln y_n\end{bmatrix} = \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ 1 & x_3 \\ \cdots & \cdots\\ \cdots & \cdots \\ 1 & x_n\end{bmatrix} \begin{bmatrix} A' \\b\end{bmatrix}$

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:

$Y = MX$

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

$X = M^{+}Y$

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.

MATLAB - Fit exponential curve WITHOUT toolbox

2 Answers

Caveat