Skip to content
Nov 1 / kkrizka

Non-Linear Fitting Using GNU/Octave And leasqr

I am currently taking the Optics Lab course up at SFU, which means that I have to take some data and then fit some function to it. The preferred way to do that is using MATLAB’s Curve Fitting Toolbox. However I only have limited access to that software on the slow computers in the common room. To work more efficiently, I’ve looked into alternatives that I can run on my laptop. The best choice that I’ve found is GNU/Octave. Octave is mostly compatible with MATLAB, so I can share the bulk of my code with rest of the class. The only issue is that Octave does not have it’s own “Curve Fitting Toolbox”, but instead relies on a single function, found in the optim package, called leasqr. Since there aren’t many examples of using this function, I decided the write one.

The following script creates a set of data following a line and adds some Gaussian noise. Then it fits a line to it, prints out the results including the fit parameters (m,b), covariance matrix and calculated chi square. At the end, it displays the plot of the data points, fitted line and the residuals.

clear all; clf;
 
% Function that will be fit
function [y]=line_func(x,par)
  y=par(1)*x+par(2);
end
 
% Generate a line
m=2;
b=1;
x=[0:0.1:10]';
y=m*x+b;
 
% Add some noise to the line
sigma=0.1;
weights=ones(size(x))/sigma;
y=y+randn(size(x))*sigma;
 
% Perform the fit
pin=[m,b];
[f,p,cvg,iter,corp,covp]=leasqr(x,y,pin,"line_func",.0001,20,weights);
 
% Print out the results
p
covp
sum((y-line_func(x,p)).^2.*weights.^2)
 
% Plots
subplot(2,1,1)
hold on;
errorbar(x,y,1./weights);
plot(x,line_func(x,p));
 
subplot(2,1,2);
errorbar(x,y-line_func(x,p),1./weights);

I won’t go through the code line-by-line, but I want to make a few key points.

  1. The equation that will be fit has to be defined as an Octave function that takes two arguments; an vector containing values of the independent variable (x) and a vector containing a list of parameters (p). The function should return a vector containing the dependent value for each x value input.
  2. The weights parameter contains the weighting information for each data point. In the case of leasqr function, it is vector containing 1/(uncertainity of point) for each data point. This is different from MATLAB, which uses 1/(uncertainity of point)^2.
  3. You can extract the uncertainty of the fitted parameters from the diagonals of the covariance matrix. The simplest way to do this is using the following code:
    sqrt(diag(covp))
  4. I could not find if any of the returned parameters are related to chi squared, so I just calculate them myself using the line:
    sum((y-line_func(x,p)).^2.*weights.^2)
  5. The leasqr function returns a value called the covariance matrix of the residuals. I haven’t found any other reference to this matrix through Google, so please leave a comment if you understand it. It is possibly related to the reduced chi square…

4 Comments

Leave a comment
  1. Cannes rentals / Dec 18 2010

    The program allows you to save predicted values and residuals. Saving residuals is useful for examining model assumptions.

  2. Lionel / Jul 7 2011

    Thanks a lot,that saved me a couple of hours work! Octave needs more examples such as this one.
    I was looking for a way to fit data to a custom equation, this is the way to go.

  3. Michael / Aug 23 2012

    Hi, I’m trying your code on my XP machine but get the following message:

    octave:11> curve_fit_octave
    error: Invalid call to options. Correct usage is:

    — Function File: OPT = options (“KEY1”, VALUE1, “KEY2”, VALUE2, …)

    error: called from:
    error: C:\octave\Octave3.6.1_gcc4.6.2\share\octave\3.6.1\m\help\print_usage.m at line 87, column 5
    error: C:\octave\Octave3.6.1_gcc4.6.2\share\octave\packages\control-2.3.52\options.m at line 68, column 5
    error: evaluating argument list element number 1
    error: C:\octave\Octave3.6.1_gcc4.6.2\share\octave\packages\optim-1.2.0\leasqr.m at line 574, column 5
    error: C:\code\ml\demo_mk\curve_fit_octave.m at line 21, column 25
    octave:11>

    Sorry to bother you with that but I really need curve fitting to work and before I dig into it myself I thought I ask around for advice first – looks like you’re the expert. Any suggestion is highly appreciated.
    Thanks
    Michael

    • Michael / Aug 23 2012

      I just ran your code on Octave 3.2.4 on Debian-squeeze without problems. So the problem might be related to XP, or 3.6.1, or both or something else…

Leave a comment
Cancel reply