% Estimerer koeffisientene beta0 og beta1 i den lineære modellen
% for enkel lineær regresjon av oksidinnhold (først XMET-, deretter
% XRF-målinger) mot mineraliseringsklasse.
% Kontruerer 100(1-alpha)%-konfidensintervall for koeffisientene.
% INPUT: .txt-fil med XRF-målinger i femte kolonne, XMET-målinger i
% sjette kolonne og mineraliseringsklasse i syvende kolonne.
% Hvis XRF-måling ikke er tilgjengelig for en observasjon, er
% verdien satt til 300
% OUTPUT: Estimat, nedre grense og øvre grense for hver parameter
% og for hvert datagrunnlag (XMET/XRF) skrives til kommandovinduet
clc, clear all, close all

alpha = 0.05;

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

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

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

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

% Beregn estimatene for parametrene
beta1 = Sxy/Sxx;
beta0 = ybar - beta1*xbar;

% Finn estimert respons
yhat = beta0 + beta1.*x;

% Antall observasjoner
n = length(y);

% Estimat for standardavviket til residualene
s2 = sum((y-yhat).^2)/(n-2);

% Finn nedre (l) og øvre (u) grense for konfidensintervallene
l1 = beta1 - tinv(1-alpha/2,n-2).*sqrt(s2/Sxx);
u1 = beta1 + tinv(1-alpha/2,n-2).*sqrt(s2/Sxx);

l0 = beta0 - tinv(1-alpha/2,n-2).*sqrt(s2*sum(x.^2)/(n.*Sxx));
u0 = beta0 + tinv(1-alpha/2,n-2).*sqrt(s2*sum(x.^2)/(n.*Sxx));

% Skriv resultatet til skjerm
fprintf('XMET-data:\n')
fprintf('Koeffisient\tEstimat\t\t%i-prosent konfidensintervall\n',100*(1-alpha))
fprintf('beta0\t\t%.4f\t\t(%.4f, %.4f)\n',beta0,l0,u0)
fprintf('beta1\t\t%.4f\t\t(%.4f, %.4f)\n\n',beta1,l1,u1)

% 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);

l1 = beta1 - tinv(1-alpha/2,n-2).*sqrt(s2/Sxx);
u1 = beta1 + tinv(1-alpha/2,n-2).*sqrt(s2/Sxx);

l0 = beta0 - tinv(1-alpha/2,n-2).*sqrt(s2*sum(x.^2)/(n.*Sxx));
u0 = beta0 + tinv(1-alpha/2,n-2).*sqrt(s2*sum(x.^2)/(n.*Sxx));

fprintf('XRF-data:\n')
fprintf('Koeffisient\tEstimat\t\t%i-prosent konfidensintervall\n',100*(1-alpha))
fprintf('beta0\t\t%.4f\t\t(%.4f, %.4f)\n',beta0,l0,u0)
fprintf('beta1\t\t%.4f\t\t(%.4f, %.4f)\n\n',beta1,l1,u1)
