lsftest.m

% lsftest - Program for testing least squares fit routines
clear; help lsftest; % Clear memory and print header
N = 50;          % Number of data points
rand('seed',0);  % Initialize random number generator
rand('normal');  % Gaussian generator with unit variance
fprintf('Enter the coefficients of the quadratic \n');
fprintf('  Y(x) = c(1) + c(2)*x + c(3)*x^2 \n');
c = input('as [c(1) c(2) c(3)] - ');
alpha = input('Enter estimated error bar - ');
x = 1:N;               % x(i)=i, i=1,...,N
r = alpha*rand(1,N);   % Gaussian distributed random vector
y = c(1) + c(2)*x + c(3)*x.^2 + r;
sigma = alpha*ones(1,N);    % Constant error bar
%%%%%% Linear regression (Straight line) fit %%%%%
M = 2;  % Number of fit parameters (M=2 is straight line)
[a_fit sig_a yy chisqr] = linreg(x,y,sigma);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Polynomial fit %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%M = 3;  % Number of fit parameters (M=3 is quadratic)
%[a_fit sig_a yy chisqr] = pollsf(x,y,sigma,M);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Fit parameters:\n');
for i=1:M
  fprintf(' a(%g) = %g +/- %g \n',i,a_fit(i),sig_a(i));
end
fprintf('Chi square = %g, N-M = %g \n',chisqr,N-M);
plot(x,y,'o',x,yy,'-');    % Plot data and curve fit
errorbar(x,y,sigma)        % Add error bars to plot
xlabel('x'); ylabel('y');
title('Data (circle); curve fit (solid line)');