
% Example commands to fit a line to data.
% Once beyond Phys1001, your data points errors vary from point to point
% So how do you fit data then?
% The example uses two artificial data sets designed ot have proper statistical properties
% likely want to read data in from a file instead.

xdata = [1.0, 2.0, 3.0, 5.0, 7.0, 10.0]
ydata = [1.64, 3.22, 3.87, 6.85, 10.96, 14.37]
%sigma = [0.20, 0.40, 0.50, 0.60, 0.70, 0.80]

% Alternate, less precise data gives larger uncertainties on slope
%ydata = [1.43, 2.10, 3.20, 6.43, 12.33, 12.75]
sigma = [0.40, 0.80, 1.00, 1.20, 1.40, 1.60]


% simple polynomial fitting with no weights.
% gives the right central value if the error bars are the same size.
% but how often does that happen ?  never in physics.
simplefit = polyfit(xdata',ydata',1)

% what we really want requires the curve fitting toolbox
% which sometimes costs extra money or has license restrictions
% octave might have a variant that is free and suitable.
%myoptions = fitoptions('poly1')
% the weights are 1/(sigma^2)
%myoptions.weights=sigma.^(-2)
%[betterfit, gof, output] = fit(xdata',ydata','poly1',myoptions)

myoptions = fitoptions('Method', 'NonlinearLeastSquares','StartPoint',[1,1])
myoptions.weights = sigma.^(-2)
ft = fittype('a*x + b','options',myoptions)
[betterfit, gof, output] = fit(xdata',ydata',ft)


myErrorbarplot = errorbar(xdata,ydata,sigma,'LineStyle','none', ...
                        'LineWidth', 2.0, 'Marker', 'o' )
% make some of the pot minimally handsome
set(myErrorbarplot, 'LineStyle' , 'none', 'Color', [0. 0. 0.], ...
                  'MarkerSize',6, 'MarkerFaceColor', [0. 0. 0.])
set(gcf,'color','white')
set(gca, 'Box', 'off', 'LineWidth',2, 'FontSize', 16)

hold on
myBetterLine = plot(betterfit,xdata,ydata)
hold off
% Make things even prettier.
myXLabel = xlabel('time (milliseconds)','FontSize',16);
myYLabel = ylabel('distance (micrometers)','FontSize',16);
myTitle = title('illustration of weighted fit','FontSize',16);
set(myBetterLine, 'LineWidth', 2.0, 'Color', [0. 0. 1.])
legend('off')


% But the above is not what a physicist wants.  We have to do one more step.
% The above best fit is right, but the uncertainties are not right.
% you can figure that out by artificially making sigma smaller
% and contrary to what you expect, the parameter uncertainties dont shrink.
% curve_fit is thinking of sigma as weights (relative size is important)
% and not as error bars (absolute size is also important).
% this will scale everything to the right absolute size.

perr_scale_factor=1.0/sqrt((gof.sse/gof.dfe))

% One sigma is 0.6827 contour is what we want.
% Two sigma is 0.9545 contour is the matlab default
ci = confint(betterfit,0.6827)
% and now scale it.
p1err = (0.5 * (ci(2) - ci(1))) * perr_scale_factor
p2err = (0.5 * (ci(4) - ci(3))) * perr_scale_factor

disp( ['contours ', num2str(ci(1)), ' ', num2str(ci(2)),' ', num2str(ci(3)),' ',num2str(ci(4))] )

% another thing we want is hiding in the output structure
% this is the covariance matrix, from which we can obtain
% a simple but very common expression for the approximate error
% by using just the diagonal elements.
% many fitting packages give this result by default.
% btw, I dont know why CovB is scaled right, but the above are not.
J=output.Jacobian
CovB=inv(J'*J)
ApproxErr=sqrt(diag(CovB))

betterfit

disp( ['CarefulErr 1sigma ', num2str(p1err), ' ', num2str(p2err)] )
disp( ['ApproxErr 1sigma ', num2str(ApproxErr(1)), ' ', num2str(ApproxErr(2)) ])
disp( ['Maybe use Approx error for your work in this class.'])

% print chisquare and dof.
disp( ['weighted chisquare ', num2str(gof.sse), ' for dof ', num2str(gof.dfe)])



% An attempt to do the same thing using a different fit routine
% this fit routine is nlinfit from the 
% Statistics and machine Learning Toolbox, probably also not free.
% This function may be available in the forge optim package of octave
% it will fit a nonlinear function, but here I'll pass a linear one.

%modelfun = @(b,x)(b(2)*x+b(1));
%beta0 = [1.0,1.0]


%[beta,myR,myJ,myCovB,MSE,ErrorModelInfo] = nlinfit(xdata',ydata',modelfun,beta0','Weights',sigma')
%% a related version is fitnlm
%% you could try to hack the standard polyfit
%% www.walkingrandomly.com/?p=2942


%unscalederr = sqrt(diag(myCovB))

%inv(myJ'*myJ)*MSE

%chisquare = myR'*myR

%[m,n] = size(J)

%scalederr = unscalederr/sqrt((chisquare/(m-n)))
