
% This version is for Octave, tested with Octave 4.0, needs optim package
% Is not exactly compatible with running on Matlab, but it is close.

pkg load optim
% run runtestfitO.m

% Example commands to fit a line to data,
% Like a spreadsheet will do, but better.
% Once beyond Phys1001, your data points' errors probably 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
% You can put your own data in, or 
% 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
% The ydata here was synthetically genreated with that sigma
%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]

% You can put the alternate data and sigma in to see.
% Or you can put the wrong sigma with the more precise data above
% Just to sense what would happen.

% simple polynomial fitting with no weights.
% will give an answer that is rough in many cases, but not right.
% 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)

disp ("The above is the incorrect fit and physics errors from an unweighted fit to data")
disp (" ")

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% An attempt to do the same thing using a different fit routine
% this fit routine is nlinfit from the 
% Statistics and machine Learning Toolbox in Matlab
% and a variant of it is in octave in the optim package
% for these two options next need to install then load pkg optim 
% pkg install -forge optim   (if you haven't installed it yet)
% pkg load optim
% it will fit an arbitrary function, but here I'll use y=(b1)x + b2.

disp (" ")
disp ("Now start a correct weighted fit, printing out stuff as we go")
disp (" ")

modelfun = @(b,x)(b(1)*x+b(2));  % this is just a line y=mx+b
beta0 = [1.0,1.0];  % the fit wants starting values
weights = sigma.^(-2) 

[beta,myR,myJ,myCovB,MSE] = nlinfit(xdata',ydata',modelfun, beta0', [], 'weights',weights')
% a related version is fitnlm, which may appear in the Octave statistics package
% which could/would be equivalent to the Curve Fitting Toolbox Matlab version


% the covariance matrix doesn't know about our sigma error bars yet
% you can tell if you rerun this putting in the "wrong" sigma above.
disp (" ")
disp ("The correct fit but incorrect error directly from covariance matrix")
disp ("Doesn't yet know we had assigned uncertainties to the data points")
disp ("The errors above are just sqrt(diagonal elements)"

unscalederr = sqrt(diag(myCovB))

disp ("So fix it up to use the sigma on the data points")

% the fit saves and returns a set of "residuals"
% which is simply the difference (data - best fit) for each point.
% you could calculate those by hand pretty easily if you needed.
% today we want the chisquare is just the sum of the squares of those.
chisquare = myR'*myR
[m,n] = size(myJ);
dof = m - n

disp (" ")

disp ("The correct fit and approximate physics errors from an unweighted fit to data")
disp ("Scaled the covariance matrix by this chisquare/dof based on our sigmas")
disp ("First one is b(1) the slope, second one is b(2) intercept")

bestfit = beta
correctscalederrors = unscalederr/sqrt((chisquare/dof))


% how to print out the value of my best fit function
%disp ("Example to write out function values")
%modelfun(beta',xdata)


% Octave did use Gnuplot for graphics.
% But might use a different interface now.

myErrorbarplot = errorbar(xdata,ydata,sigma,"~");
% make some of the pot minimally handsome
set(myErrorbarplot, 'LineStyle' , 'none', 'Color', [0. 0. 0.], ...
                  'MarkerSize',6, 'MarkerFaceColor', [0. 0. 0.], ...
                  'LineWidth', 2.0, 'Marker', 'o' )

% also set some generic graphics options for the axes and titles
set(gcf,'color','white')
set(gca, 'Box', 'off', 'LineWidth',2, 'FontSize', 16)

% now plot the best fit line.
hold on
myLineplot = plot(xdata,modelfun(beta,xdata));
set(myLineplot,'Color', [0. 0. 0.7], 'LineWidth', 2.0 )
hold off

% Make things even prettier and communicative
myXLabel = xlabel('time (milliseconds)','FontSize',16);
myYLabel = ylabel('distance (micrometers)','FontSize',16);
myTitle = title('illustration of weighted fit in Octave','FontSize',16);
legend('off')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Another method available in octave base, only for polynomials.
% Saves you from having to code in a function, still allows
% weighted fits.
% you could try to hack the standard polyfit
% www.walkingrandomly.com/?p=2942
% This whole thing is commented out using #{ }#
#{
[poly1,S1] = wpolyfit(xdata', ydata', sigma',1)

% Errors unscaled
dp = sqrt(sumsq(inv(S1.R'))'/S1.df)*S1.normr

inv(S1.R')
% this is the diagonal elements of the covariance matrix.
sqrt(sumsq(inv(S1.R'))')

%  But thats not what I want actually, so undo it.

chisquare2 = S1.normr * S1.normr

newdp = dp/(sqrt(chisquare2/S1.df))

disp ("The above is the correct physics errors from a fit to data")
#}