× {{alert.msg}} Never ask again
Get notified about new tutorials RECEIVE NEW TUTORIALS

MATLAB - Fit exponential curve WITHOUT toolbox

Ray Phan
Mar 31, 2015
<p>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 <a href="http://en.wikipedia.org/wiki/Least_squares" rel="nofollow">least squares</a>. 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. </p> <p>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 <code>x</code> and output <code>y</code> is:</p> <p><img src="http://latex.codecogs.com/gif.latex?y%20%3D%20Ae%5E%7Bbx%7D" alt=""></p> <p>To get this to be "linear", we can take the natural logarithm of both sides:</p> <p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D%5Cln%20y%20%3D%20%5Cln%28Ae%5E%7Bbx%7D%29%5Cend%7Balign*%7D" alt=""></p> <p>By using the fact that <code>ln(ab) = ln(a) + ln(b)</code>, we have:</p> <p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D%5Cln%20y%20%3D%20%5Cln%28A%29%20+%20%5Cln%28e%5E%7Bbx%7D%29%5Cend%7Balign*%7D" alt=""></p> <p>Also knowing that <img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D%5Cln%28e%5Ez%29%20%3D%20z%5Cln%28e%29%20%3D%20z%5Cend%7Balign*%7D" alt="">, this simplifies to:</p> <p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D%5Cln%20y%20%3D%20%5Cln%20A%20+%20bx%5Cend%7Balign*%7D" alt=""></p> <p>As you can see, the above equation is now "linear" with respect to log-space. Given a bunch of <code>x</code> and <code>y</code> values, so <code>(x_1, x_2, ..., x_n)</code> and <code>(y_1, y_2, ..., y_n)</code>, we can concatenate a bunch of equations together in a linear system:</p> <p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Balign*%7D%5Cln%20y_1%20%26%3D%20%5Cln%20A%20+%20bx_1%5C%5C%5Cln%20y_2%20%26%3D%20%5Cln%20A%20+%20bx_2%5C%5C%5Cln%20y_3%20%26%3D%20%5Cln%20A%20+%20bx_3%5C%5C%26%5Cldots%5C%5C%26%5Cldots%5C%5C%5Cln%20y_n%20%26%3D%20%5Cln%20A%20+%20bx_n%5C%5C%5Cend%7Balign*%7D" alt=""></p> <p>If we let <code>ln(A) = A'</code> for ease of notation, and rearranging this so that it's in matrix form, we get:</p> <p><img src="http://latex.codecogs.com/gif.latex?%5Cbegin%7Bbmatrix%7D%5Cln%20y_1%20%5C%5C%20%5Cln%20y_2%20%5C%5C%20%5Cln%20y_3%20%5C%5C%20%5Ccdots%20%5C%5C%20%5Ccdots%20%5C%5C%20%5Cln%20y_n%5Cend%7Bbmatrix%7D%20%3D%20%5Cbegin%7Bbmatrix%7D%201%20%26%20x_1%20%5C%5C%201%20%26%20x_2%20%5C%5C%201%20%26%20x_3%20%5C%5C%20%5Ccdots%20%26%20%5Ccdots%5C%5C%20%5Ccdots%20%26%20%5Ccdots%20%5C%5C%201%20%26%20x_n%5Cend%7Bbmatrix%7D%20%5Cbegin%7Bbmatrix%7D%20A%27%20%5C%5Cb%5Cend%7Bbmatrix%7D" alt=""></p> <p>Therefore, we simply need to solve for <code>A'</code> and <code>b</code>, and you can do that by the <a href="http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse" rel="nofollow">pseudoinverse</a>. Specifically, the above problem is of the form:</p> <p><img src="http://latex.codecogs.com/gif.latex?Y%20%3D%20MX" alt=""></p> <p>Therefore, we need to solve for <code>X</code>, and so:</p> <p><img src="http://latex.codecogs.com/gif.latex?X%20%3D%20M%5E%7B+%7DY" alt=""></p> <p><code>M^{+}</code> is the pseudoinverse of the matrix. Once you're done, simply take the <code>exp</code> operator on <code>A'</code> to get the original <code>A</code>. MATLAB has very efficient linear system solvers and least-squares solvers. Specifically, you can use the <code>\</code> or <code>ldivide</code> operator. All you have to do is create the <code>M</code> matrix from the <code>x</code> values, create a vector of <code>y</code> values and solve your system. It's really simply:</p> <pre><code>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 </code></pre> <p>Therefore, using your <code>x</code> and <code>y</code> values, this is what I get:</p> <pre><code>A = 1.3882 b = -11.508 </code></pre> <p>If we plot the above points as well as an exponential curve that fits the line, we can do:</p> <pre><code>xval = linspace(min(x), max(x)); yval = A*exp(b*xval); plot(x,y,'r.',xval,yval,'b'); </code></pre> <p>The first line of code defines a bunch of <code>x</code> values that span between the smallest and largest <code>x</code> value for our data set. For the next line, we then take the <code>x</code> 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.</p> <p>We get:</p> <p><img src="http://i.stack.imgur.com/3oIuc.png" alt="enter image description here"></p> <p>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.</p> <p>This tip was originally posted on <a href="http://stackoverflow.com/questions/29305018/MATLAB%20-%20Fit%20exponential%20curve%20WITHOUT%20toolbox/29305738">Stack Overflow</a>.</p>
comments powered by Disqus