% Matlab-code used to make the illustrations in the presentation
% FYS2150-linjetilpasning-2019
%
% Alex Read 04.03.2019

% 10 fake datapoints with the same normal uncertainty at each point
x = linspace(1,10,10);
mmodel = 0.25; bmodel = 2; sigma_y = 0.5;
y0 = mmodel*x+bmodel; % The modell
y = normrnd(y0,sigma_y); % The fake data generated from the model

% Plot the data as points with error bars
figure(1)
hold off
errorbar(x,y,ones(length(y),1)*sigma_y,'k.','MarkerSize',20)
hold on
xlabel('x'), ylabel('y')

% Fit data to first order polynomial (line)
% Save the polynomial coefficients p and the error structure S
% (Matlab is using "Squires" here, so NOT using the specified
% uncertainty per measurement).
[p,S] = polyfit(x,y,1);
m=p(1);
b=p(2);

% Calculate the covariance matrix C from S (see help for polyfit)
Rinv=inv(S.R);
C = (Rinv*Rinv')*S.normr^2/S.df;

% The diagonal elements of C are the variances of the parameters.
% Non-zero off-diagonal elements indicate a covariance between the
% parameters.
dm2 = C(1,1); % Error in the slope
dc2 = C(2,2); % Error in the intercept
dmc = C(1,2);

% Evaluate and plot the model and the 
% fitted result over a much wider range,
% in order to e.g. see the x-intercept
%xwide = linspace(-15,12,100);
xwide = linspace(0,12,100);
yfit = polyval(p,xwide);
y0wide = mmodel*xwide+bmodel;
plot(xwide,y0wide,'-')
plot(xwide,yfit,'r--')

% Error in y as function of x
dy = sqrt(xwide.^2*dm2 + dc2 + 2*xwide*dmc);

% Plot error in y as function of x as a transparent green band
f=fill([xwide flip(xwide)],[yfit-dy flip(yfit+dy)],'g');
f.FaceAlpha=0.3;

% Put the grid on so we can see the y and x intercepts
grid on

% Error in x-intercept
%dx = sqrt(sigma_y^2/m^2 + dc2/m^2 + (b/m^2)^2*dm2 + 2*(b/m^2)*(-1/m)*dmc);

% Error in x given y and error in y
yobs = 2.5;
dx = sqrt(sigma_y^2/m^2 + dc2/m^2 + ((yobs-b)/m^2)^2*dm2 + 2*(-(yobs-b)/m^2)*(-1/m)*dmc);
dxp = sqrt(sigma_y^2/m^2);

% Show the x-intercept and its uncertainty
errorbar((yobs-b)/m,yobs,sigma_y,sigma_y,dx,dx,'r.','MarkerSize',20)
errorbar((yobs-b)/m,yobs,sigma_y,sigma_y,dxp,dxp,'b.','MarkerSize',10)

legend('data','modell','tilpasning','usikkerhet','y-målt','y-målt bare \sigma_y')

% Do some checks to compare the Matlab results with Squires and with the
% chi-squared method that uses a specified uncertainty per measurement
% point (here they are all equal and equal to the uncertainty specified in
% the model).

% 1) Squires
n=max(size(x));
D=sum(x.^2)-(1/n).*sum(x).^2;
E=sum(x.*y)-(1/n)*sum(x)*sum(y);
F=sum(y.^2)-(1/n)*sum(y).^2;
deltamSquires=sqrt((D*F-E^2)/(D^2*(n-2)));
deltacSquires=sqrt((n-2)^-1*(D/n+(mean(x))^2)*((D*F-E^2)/D^2));
deltamcSquires = -mean(x)*deltamSquires^2;

% 2) Matlab
deltamMatlab = sqrt(dm2);
deltacMatlab = sqrt(dc2);
deltamcMatlab = dmc;

% 3) chi-squared method (formulas in presentation)
Sigma1 = n/sigma_y^2;
deltamPres = sqrt(1/(mean(x.*x)-mean(x)^2)/Sigma1);
deltacPres = deltamPres*sqrt(mean(x.*x));
deltamcPres = -mean(x)*deltamPres^2;

fprintf('Single random sample results:\n\n')
fprintf('                   Squires  Matlab  Chi-squared\n')
fprintf('Error in slope     %7.4f %7.4f %7.4f \n',deltamSquires,deltamMatlab,deltamPres)
fprintf('Error in intercept %7.4f %7.4f %7.4f \n',deltacSquires,deltacMatlab,deltacPres)
fprintf('Covariance         %7.4f %7.4f %7.4f \n',deltamcSquires,deltamcMatlab,deltamcPres)


% Repeat the above calculations many times and make histograms
N=1000;
for I=1:N
    y = normrnd(y0,sigma_y);
    [p,S] = polyfit(x,y,1);
    m(I)=p(1);
    b(I)=p(2);

% Calculate the covariance matrix C from S (see help for polyfit)
    Rinv=inv(S.R);
    C = (Rinv*Rinv')*S.normr^2/S.df;
    
    D=sum(x.^2)-(1/n).*sum(x).^2;
    E=sum(x.*y)-(1/n)*sum(x)*sum(y);
    F=sum(y.^2)-(1/n)*sum(y).^2;
    deltamSquires2(I)=(D*F-E^2)/(D^2*(n-2));
    deltacSquires2(I)=(n-2)^-1*(D/n+(mean(x))^2)*((D*F-E^2)/D^2);
    deltamcSquires(I) = -mean(x)*deltamSquires2(I);
    
    % 2) Matlab
    deltamMatlab(I) = sqrt(dm2);
    deltacMatlab(I) = sqrt(dc2);
    deltamcMatlab(I) = dmc;
    
    % 3) Chi-squared (actually does not depend on y-data!!)
    Sigma1 = n/sigma_y^2;
    deltamPres(I) = sqrt(1/(mean(x.*x)-mean(x)^2)/Sigma1);
    deltacPres(I) = deltamPres(I)*sqrt(mean(x.*x));
    deltamcPres(I) = -mean(x)*deltamPres(I)^2;
end
figure(2), plot(m,b,'.'), xlabel('m'),ylabel('b')
fprintf('\nLarge sample results (N=%d):\n\n',N)
fprintf('                   Squires mean\n')
fprintf('Error in slope     %7.4f \n',sqrt(mean(deltamSquires2)))
fprintf('Error in intercept %7.4f \n',sqrt(mean(deltacSquires2)))
fprintf('Covariance         %7.4f \n',mean(deltamcSquires))

figure(3),histogram(deltamSquires2),xlabel('Squires (\Delta m)^2')
figure(4),histogram(deltacSquires2),xlabel('Squires (\Delta b)^2')
figure(5),histogram(deltamcSquires),xlabel('Squires (\Delta mb)')

fprintf('\nThe actual distribution of b versus m\n')
fprintf('   sqrt(variance of m) %7.4f\n',std(m))
fprintf('   sqrt(variance of b) %7.4f\n',std(b))
fprintf('   covariance of m,b   %7.4f\n',mean((m-mean(m)).*(b-mean(b))))