Segment 20 Group Activity 2014

Code:

clear all
close all
clc

%% Input

load chiSqData.mat;
xdata = A(:,1);
ydata = A(:,2);
sigmad = A(:,3);

% plotting stuff
color = 'rgbcmyk';
fign = 1;

%% Errorbars
figure(fign)
errorbar(xdata,ydata,sigmad,'o')

%% Models
nm = 4;
model = cell(nm,2);
model{1,1} = @(x,b) b(1) + b(2)*x.^2; model{1,2} = 'Parabola';
model{2,1} = @(x,b) b(1) + b(2)*x; model{2,2} = 'Linear';
model{3,1} = @(x,b) b(1)*exp(b(2)*x); model{3,2} = 'Exponential';
model{4,1} = @(x,b) b(1) + b(2)*x + b(3)*x.^2; model{4,2} = 'GenQuad';
bguess = cell(nm,1);
bguess{1} = [1 5]';
bguess{2} = [1 1]';
bguess{3} = [1 1]';
bguess{4} = [1 1 1]';

%% Initialize stuff
bfit = cell(nm,1);
chisqmin = cell(nm,1);
covar = cell(nm,1);
stderr = cell(nm,1);
lstr = cell(nm+1,1); lstr{1} = 'Data';


%% Fitting
for count = 1:length(model)
    ymodel = model{count,1};
    [bfit{count}, chisqmin{count}, covar{count}] = fitCoeff(ymodel,xdata,ydata,sigmad,bguess{count},fign,color(count));
    stderr{count} = sqrt(diag(covar{count}));
    lstr{count+1} = model{count,2};
end
figure(fign)
legend(lstr,'Location','Best')

%% Print chisqmin values
disp('Chisquarmin values:')
chisqmin

function [bfit, chisqmin, covar] = fitCoeff(ymodel,xdata,ydata,sigmad,bguess,fign,color)

    % Fitting
    chisq = @(b) sum(((ymodel(xdata,b)-ydata)./sigmad).^2);
    bfit = fminsearch(chisq,bguess);
    chisqmin = chisq(bfit);
    xfit = 0:0.01:1.1*max(xdata);
    
    % Hessian calculation
    n = length(bguess);
    unit = @(i) (1:n)' == i;
    hess = zeros(n,n);
    h = 0.1;
    for i = 1:n
        for j = 1:n
            bpp = bfit + h*(unit(i)+unit(j));
            bmm = bfit + h*(-unit(i)-unit(j));
            bpm = bfit + h*(unit(i)-unit(j));
            bmp = bfit + h*(-unit(i)+unit(j));
            hess(i,j) = (chisq(bpp)+chisq(bmm)-chisq(bpm)-chisq(bmp))./(2*h)^2;
        end 
    end
    covar = inv(0.5*hess);
    
    % Plotting
    figure(fign)
    hold on
    plot(xfit,ymodel(xfit,bfit),['-' color])
    
end