% Utfører enkel lineær regresjon av oksidinnhold, y, mot minerali-
% seringsklasse, x (først XMET-data, deretter XRF-data)
% INPUT: .txt-fil med XRF-målinger i femte kolonne, XMET-målinger
% i sjette kolonne og mineraliseringsklasse i syvende kolonne.
% For observasjoner som ikke har XRF-målinger tilgjengelig, er
% verdien satt til 300 for å markere dette.
% OUTPUT: Én figur med regresjonslinje og spredningsplott av data for
% XMET-data til venstre og XRF-data til høyre.
% Én figur delt inn i tre rader og to kolonner, med residualplott
% (øverst), histogram (midten) og Q-Q plott (nederst) av residualene
% etter regresjon med XMET-data (venstre) og med XRF-data (høyre).
clc, clear all, close all

% Last inn data
data = load('malmdata.txt');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% XMET-data
x = data(:,7);
y = data(:,6);

% Finn grenser for plott-utsnitt
ymin = floor(min(y)); % Største heltall mindre enn minste observasjon
ymax = ceil(max(y));  % Minste heltall større enn største observasjon

% Gjennomsnitt
xbar = sum(x)/length(x);
ybar = sum(y)/length(y);

% Kvadratsummer
Sxx = sum((x-xbar).^2);
Syy = sum((y-ybar).^2);
Sxy = sum((x-xbar).*(y-ybar));

% Koeffisienter
beta1 = Sxy/Sxx;
beta0 = ybar - beta1*xbar;

% Forventet respons, SSE og R^2
yhat = beta0 + beta1.*x;
SSE = sum((y - yhat).^2);
R2 = 1 - SSE/Syy;

% Figur med regresjonslinje og spredningsplott
figure(1)
subplot(1,2,1), hold on
plot(x,y,'kx')
plot(sort(x),sort(yhat),'r-')
xlim([0 4])
box on
legend('XMET-data', 'Regresjonslinje','Location','NorthWest')
title(['{\beta_0=}' num2str(beta0) '{, \beta_1=}' num2str(beta1) '{, R^2 =}' num2str(R2)],'FontSize',12)
xlabel('Klasse')
ylabel('Oksidinnhold (XMET)')
ylim([ymin ymax])
hold off

% Residualer
epsilon = y - yhat;
epsilonmin = floor(min(epsilon));
epsilonmax = ceil(max(epsilon));

% Lag residualplott
figure(2)
subplot(3,2,1), hold on
plot(epsilon,'b.')
plot(zeros(size(epsilon)),'k--','LineWidth',2)
xlabel('Observasjon nr.')
ylabel('{\epsilon}')
ylim([epsilonmin epsilonmax])
xlim([0 length(epsilon)])
legend('Residualer','Location','SouthEast')
title('Residualplott, XMET-data')
box on
hold off

% Histogram
n = 12; % Antall delintervaller
edges = min(epsilon):range(epsilon)/(n-1):max(epsilon);
count = hist(epsilon,edges);
area = sum(count)*range(epsilon)/n;

figure(2)
subplot(3,2,3), hold on
bar(edges,count/area,1)
shading faceted
alpha(0.5)
xlabel('{\epsilon}')
ylabel('Frekvens')
title('Histogram, XMET-data')
box on
hold off

% Q-Q plott
figure(2)
subplot(3,2,5), hold on
qqplot(epsilon)
title('Q-Q plott, XMET-data')
box on
hold off

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

x = XRFdata(:,7);
y = XRFdata(:,5);

xbar = sum(x)/length(x);
ybar = sum(y)/length(y);

Sxx = sum((x-xbar).^2);
Syy = sum((y-ybar).^2);
Sxy = sum((x-xbar).*(y-ybar));

beta1 = Sxy/Sxx;
beta0 = ybar - beta1*xbar;

yhat = beta0 + beta1.*x;
SSE = sum((y - yhat).^2);
R2 = 1 - SSE/Syy;

figure(1)
subplot(1,2,2), hold on
plot(x,y,'kx')
plot(sort(x),sort(yhat),'r-')
xlim([0 4])
box on
legend('XRF-data', 'Regresjonslinje','Location','NorthWest')
title(['{\beta_0=}' num2str(beta0) '{, \beta_1=}' num2str(beta1) '{, R^2 =}' num2str(R2)],'FontSize',12)
xlabel('Klasse')
ylabel('Oksidinnhold (XRF)')
ylim([ymin ymax])
hold off

% Residualer
epsilon = y - yhat;

% Lag residualplott
figure(2)
subplot(3,2,2), hold on
plot(epsilon,'b.','LineWidth',0.5)
plot(zeros(size(epsilon)),'k--','LineWidth',2)
xlabel('Observasjon nr.')
ylabel('{\epsilon}')
ylim([epsilonmin epsilonmax])
xlim([0 length(epsilon)])
legend('Residualer','Location','SouthEast')
title('Residualplott, XRF-data')
box on
hold off

% Histogram
n = 12; % Antall delintervaller
edges = min(epsilon):range(epsilon)/(n-1):max(epsilon);
count = hist(epsilon,edges);
area = sum(count)*range(epsilon)/n;

figure(2)
subplot(3,2,4), hold on
bar(edges,count/area,1)
shading faceted
alpha(0.5)
xlabel('{\epsilon}')
ylabel('Frekvens')
title('Histogram, XRF-data')
box on
hold off

% Q-Q plott
figure(2)
subplot(3,2,6), hold on
qqplot(epsilon)
title('Q-Q plott, XRF-data')
box on
hold off

