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.
- 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.
- 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.
- 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:
- I could not find if any of the returned parameters are related to chi squared, so I just calculate them myself using the line:
- 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…