% Plotter residualplott, histogram og Q-Q plott for differansene y_2i - y_1i,
% i = 1,...,103, hvor y_2i er XMET-observasjonene, og y_1i er XRF-observasjonene
% INPUT: .txt-fil med XRF-målinger i femte kolonne og XMET-målinger i sjette
% kolonne. For observasjoner hvor det ikke finnes XRF-målinger, er verdien satt
% til 300 for å markere dette.
% OUTPUT: Tredelt figur med residualplott (øverst), histogram med relative
% frekvenser (nederst, venstre) og Q-Q plott (nederst, høyre)
clc, clear all, close all

data = load('malmdata.txt');

% Finn XRF-observasjonene
findXRF = find(data(:,5) < 300);

y1 = data(findXRF,5); % Lagre XRF-målingene
y2 = data(findXRF,6); % Lagre XMET-målingene

epsilon = y2 - y1;

% Estimater for forventningsverdien og variansen til epsilon
epsilon_bar = mean(epsilon);
tau_hat = std(epsilon);

% Histogram med n delintervaller
n = 12;
edges = linspace(min(epsilon),max(epsilon),n+1);
count = histc(epsilon,edges);
area = sum(count).*range(epsilon)/n; % Beregn arealet av søylene for standardisering.

% Lag figur
figure(1)

subplot(2,2,1:2), hold on
plot(epsilon,'b*','LineWidth',1.2)
plot(zeros(length(epsilon),1),'k--','LineWidth',1)
title('{Støy}','FontSize',14)
xlim([0 length(epsilon)+1])
xlabel('Observasjon nr.')
ylabel('{\epsilon}')
box on
hold off

subplot(2,2,3), hold on
bar(edges,count/area,1); % Søylediagram med relative frekvenser
% Bredden til søylene settes til 1 slik at de kommer helt inntil hverandre.
shading faceted % markerer kanten på søylene
alpha(0.5) & gjør søylene delvis gjennomsiktige
title('Histogram for {\epsilon = y_2 - y_1}')
xlabel('{\epsilon}')
ylabel('Relativ frekvens')
box on
hold off

subplot(2,2,4)
qqplot(epsilon)
box on
