% Estimerer koeffisientene beta0 og beta1 i den lineære regresjons-
% modellen for regresjon av oksidinnhold, y, på mineraliserings-
% klasse, x (først for XMET-data, deretter for XRF-data).
% Konstruerer 100(1-alpha)%-konfidensintervaller for estimatene.
% INPUT: .txt-fil med XRF-data i femte kolonne, XMET-data i sjette
% kolonne og mineraliseringsklasse i syvende kolonne. For de obs-
% ervasjonene hvor XRF-målinger ikke foreligger, er XRF-måleverdien
% satt til 300 for å indikere dette.
% OUTPUT: Estimatene samt nedre og øvre grense for konfidensinter-
% vallene skrives til kommandovinduet.
% En todelt figur med forventet respons, konfidensgrenser og data-
% punkter for XMET-data (venstre) og XRF-data (høyre).
% En figur med forventet respons, estimert standardavvik
% og datapunkter for XMET-data.
clc, clear all, close all

alpha = 0.05;

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

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

ymin = floor(min(y));
ymax = ceil(max(y));

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;

n = length(y);
s2 = sum((y-yhat).^2)/(n-2);

yhat1 = beta0 + beta1;
yhat2 = beta0 + 2*beta1;
yhat3 = beta0 + 3*beta1;

conflim = @(x) tinv(1-alpha/2,n-2).*sqrt(s2.*(1/n + (x-xbar).^2/Sxx));

l1 = yhat1 - conflim(1);
u1 = yhat1 + conflim(1);

l2 = yhat1 - conflim(2);
u2 = yhat1 + conflim(2);

l3 = yhat1 - conflim(3);
u3 = yhat1 + conflim(3);

fprintf('XMET-data:\n')
fprintf('Forventningsverdi\tEstimat\t\t%i-prosent konfidensintervall\n',100*(1-alpha))
fprintf('mu1\t\t\t%.4f\t\t(%.4f, %.4f)\n',yhat1,l1,u1)
fprintf('mu2\t\t\t%.4f\t\t(%.4f, %.4f)\n',yhat2,l2,u2)
fprintf('mu3\t\t\t%.4f\t\t(%.4f, %.4f)\n\n',yhat3,l3,u3)

% Lag figur som viser forventet respons og konfidensgrenser
figure(1)
subplot(1,2,1), hold on
plot(x,y,'g+','LineWidth',0.6)
errorbar([1 2 3],[yhat1 yhat2 yhat3],conflim(1:3),'ks-')
title([num2str(100*(1-alpha)) '-prosent konfidensgrenser, XMET-data'],'FontWeight','bold')
box on
legend('XMET-data','Forventet respons','Location','NorthWest')
xlabel('Klasse')
ylabel('Oksidinnhold (XMET)')
ylim([ymin ymax])
hold off

% Lag figur som viser estimert standardavvik (for sammenlikning med
% alternativ modell senere)
figure(2), hold on
plot(x,y,'g+','LineWidth',0.6)
errorbar([1 2 3],[yhat1 yhat2 yhat3],sqrt(s2)*ones(1,3),'ks--')
legend('XMET-data','Forventet respons','Location','NorthWest')
xlabel('Klasse')
ylabel('Oksidinnhold (XMET)')
title(['{\beta_0 = }' num2str(beta0) '{, \beta_1 = }' num2str(beta1) '{, s = }' num2str(sqrt(s2))])
ylim([ymin ymax])
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;

n = length(y);
s2 = sum((y-yhat).^2)/(n-2);

yhat1 = beta0 + beta1;
yhat2 = beta0 + 2*beta1;
yhat3 = beta0 + 3*beta1;

conflim = @(x) tinv(1-alpha/2,n-2).*sqrt(s2.*(1/n + (x-xbar).^2/Sxx));

l1 = yhat1 - conflim(1);
u1 = yhat1 + conflim(1);

l2 = yhat1 - conflim(2);
u2 = yhat1 + conflim(2);

l3 = yhat1 - conflim(3);
u3 = yhat1 + conflim(3);

fprintf('XRF-data:\n')
fprintf('Forventningsverdi\tEstimat\t\t%i-prosent konfidensintervall\n',100*(1-alpha))
fprintf('mu1\t\t\t%.4f\t\t(%.4f, %.4f)\n',yhat1,l1,u1)
fprintf('mu2\t\t\t%.4f\t\t(%.4f, %.4f)\n',yhat2,l2,u2)
fprintf('mu3\t\t\t%.4f\t\t(%.4f, %.4f)\n\n',yhat3,l3,u3)

figure(1)
subplot(1,2,2), hold on
plot(x,y,'g+','LineWidth',0.6)
errorbar([1 2 3],[yhat1 yhat2 yhat3],conflim(1:3),'ks-')
title([num2str(100*(1-alpha)) '-prosent konfidensgrenser, XRF-data'],'FontWeight','bold')
box on
legend('XRF-data','Forventet respons','Location','NorthWest')
xlabel('Klasse')
ylabel('Oksidinnhold (XRF)')
ylim([ymin ymax])
hold off
