%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eksperimentel Fysik, Niels Bohr Institutet, Københavns Universitet
% Tirsdag d. 3. maj og torsdag d. 5. maj, 2011
%
% Dette dokument indeholder løsninger til opgaven med fitning. Disse opgaver
% blev lavet af Kim Lefmann og lagt op som et PDF dokument på Absalon.
% 
% Pia Jensen, bozack@bozack.dk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all     %tømmer alle gemte variable
close all     %luk alle eksterne vinduer (figurer o.lign.)
clc           %rens kommando-vinduet


%% Opgave 1 og 2 ----------------------------------------------------------------
%  I opgave 1 skal der bare fittes en enkelt gang, mens der i opgave 2 skal gøres
%  det samme mange gange i træk, og bagefter skal der ses på statistikken

% Sæt konstanterne
A_model = 0;
B_model = 1;
sigma   = 0.5;

% Konstanter for gentagne fits (opgave 2 del)
N_iter     = 1000;
A_liste    = zeros(1,N_iter);
B_liste    = zeros(1,N_iter);
chi2_liste = zeros(1,N_iter);

for i = 1:N_iter   % Loop over forskelige fits, opgave 2 del
    
    % Bestem funktionspunkter
    x       = -2:0.1:2;
    y_model = A_model + B_model*x;
    sigma_y = sigma*ones(size(y_model));
    y       = y_model + sigma_y.*randn(size(y_model));
    
    % Plot data og model
    errorbar(x,y,sigma_y,'bo','MarkerFaceColor', [0.4 0.4 0.8],...
          'MarkerEdgeColor',[0.4 0.4 0.8])
    hold on
    plot(x,y_model,'r-','LineWidth',2)
    
    % Beregn regressionsparametre
    N     = length(x);
    Sy    = sum(y);
    Sx    = sum(x);
    Sxx   = sum(x.^2);
    Sxy   = sum(x.*y);
    Delta = N*Sxx-Sx^2;
    
    % Beregn modelparametre
    A_fit = ( Sxx*Sy - Sx*Sxy )/Delta;
    B_fit = ( N*Sxy - Sx*Sy )/Delta;
    
    % Beregn modellinie
    y_fit = A_fit + B_fit*x;
    plot(x,y_fit,'g-','LineWidth',2)
    xlabel('Datapunkter x')
    ylabel('Datapunkter y')
    hold off
    
    % Formatér figur
    figure(1)
    ah1 = gca;
    set(ah1,'Box','off')
    set(ah1,'XGrid','on')
    set(ah1,'YGrid','on')
    set(ah1,'GridLineStyle',':')
    set(ah1,'FontSize',8)
    
    % Beregn chi-square
    chi2    = sum((y-y_fit).^2 ./ sigma_y.^2);
    chi2red = chi2/(N-2);
    
    % Gem resultater
    A_liste(i)    = A_fit;
    B_liste(i)    = B_fit;
    chi2_liste(i) = chi2red;
    
end   % ... loop over forskellige fits

% Statistik over forskellige fitteresultater
figure(1)
hist(A_liste,[-0.4:0.02:0.4])
xlabel('A-parameter')
ylabel('Antal i histogram-bin')
figure(2)
hist(B_liste,[0.6:0.02:1.4])
xlabel('B-parameter')
ylabel('Antal i histogram-bin')
figure(3)
hist(chi2_liste,[0:0.05:2])
xlabel('Chi^2')
ylabel('Antal i histogram-bin')

% Korrelation mellem A og B fra fittene
figure(4)
plot(A_liste,B_liste,'.','MarkerFaceColor', [0.4 0.4 0.8],...
          'MarkerEdgeColor',[0.4 0.4 0.8])
xlabel('A-parameter')
ylabel('B-parameter')

% Formatér figurer
figure(1)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',8)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(2)
ah1 = gca;
%set(ah1,'FaceColor',[0.4 0.4 0.8],'EdgeColor',[0.4 0.4 0.8])
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',8)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(3)
ah1 = gca;
%set(ah1,'FaceColor',[0.4 0.4 0.8],'EdgeColor',[0.4 0.4 0.8])
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',8)
h = findobj(gca,'Type','patch');
set(h,'FaceColor',[0.4 0.4 0.8],'EdgeColor','k')

figure(4)
ah1 = gca;
set(ah1,'Box','off')
set(ah1,'XGrid','on')
set(ah1,'YGrid','on')
set(ah1,'GridLineStyle',':')
set(ah1,'FontSize',8)


%% Opgave 3 ---------------------------------------------------------------------
%  I opgave 3 skulle man prøve at bruge MATLAB's indbyggede fit() funktion til
%  at lave fittet

% Brug MATLABs indbyggede fit 
[res good] = fit(x',y','poly1')

% Plot resultatet
errorbar(x,y,sigma_y,'ro')
hold on
plot(x,y_model,'b-')
plot(x,res(x),'g-')