% Regner ut estimat og konfidensintervall for forventet oksidinnhold for
% hver mineraliseringsklasse
% INPUT: .txt-fil med XRF-målinger i femte kolonne, XMET-målinger i sjette
% kolonne og mineraliseringsklasse i syvende kolonne.
% XRF-måleverdien er satt til 300 for å markere de observasjonene hvor XRF-
% målinger ikke er tilgjengelige.
% OUTPUT: Matrisen CI som inneholder estimatene i første kolonne, nedre
% grenser for konfidensintervall i andre kolonne, øvre grenser for konfidens-
% intervall i tredje kolonne og antall frihetsgrader i fjerde kolonne.
% Rekkefølgen på radene er klasse 1, 2, 3 (XMET), klasse 1, 2, 3 (XRF).
clc, clear all, close all

% Last inn data
data = load('malmdata.txt');
findXRF = find(data(:,5) < 300);
XRF = data(findXRF,:);

% Del inn i klasser
% Finn først indeksene til klasse 1-, 2- og 3-observasjonene
find_XMET_1 = find(data(:,7) == 1);
find_XMET_2 = find(data(:,7) == 2);
find_XMET_3 = find(data(:,7) == 3);

XMET_1 = data(find_XMET_1,6);
XMET_2 = data(find_XMET_2,6);
XMET_3 = data(find_XMET_3,6);

find_XRF_1 = find(XRF(:,7) == 1);
find_XRF_2 = find(XRF(:,7) == 2);
find_XRF_3 = find(XRF(:,7) == 3);

XRF_1 = XRF(find_XRF_1,5);
XRF_2 = XRF(find_XRF_2,5);
XRF_3 = XRF(find_XRF_3,5);

% Regn ut estimat for forventningsverdier
ybar_XMET_1 = mean(XMET_1);
ybar_XMET_2 = mean(XMET_2);
ybar_XMET_3 = mean(XMET_3);

ybar_XRF_1 = mean(XRF_1);
ybar_XRF_2 = mean(XRF_2);
ybar_XRF_3 = mean(XRF_3);

% Regn ut estimat for varians
s2_XMET_1 = var(XMET_1);
s2_XMET_2 = var(XMET_2);
s2_XMET_3 = var(XMET_3);

s2_XRF_1 = var(XRF_1);
s2_XRF_2 = var(XRF_2);
s2_XRF_3 = var(XRF_3);

% Sett konfidensnivå, og lagre antall observasjoner i hver klasse
alpha = 0.05;
n_XMET_1 = length(XMET_1);
n_XMET_2 = length(XMET_2);
n_XMET_3 = length(XMET_3);

n_XRF_1 = length(XRF_1);
n_XRF_2 = length(XRF_2);
n_XRF_3 = length(XRF_3);

% Regn ut halve lengden av konfidensintervallene ved å evaluere
% den inverse kumulative fordelingsfunksjonen til t-fordelingen
% med gitt antall frihetsgrader
d_XMET_1 = tinv(1-alpha/2,n_XMET_1-1).*sqrt(s2_XMET_1/n_XMET_1);
d_XMET_2 = tinv(1-alpha/2,n_XMET_2-1).*sqrt(s2_XMET_2/n_XMET_2);
d_XMET_3 = tinv(1-alpha/2,n_XMET_3-1).*sqrt(s2_XMET_3/n_XMET_3);

d_XRF_1 = tinv(1-alpha/2,n_XRF_1-1).*sqrt(s2_XRF_1/n_XRF_1);
d_XRF_2 = tinv(1-alpha/2,n_XRF_2-1).*sqrt(s2_XRF_2/n_XRF_2);
d_XRF_3 = tinv(1-alpha/2,n_XRF_3-1).*sqrt(s2_XRF_3/n_XRF_3);

% Konstruer konfidensintervallene
format short g;

CI = [ybar_XMET_1 ybar_XMET_1-d_XMET_1 ybar_XMET_1+d_XMET_1 n_XMET_1-1;
      ybar_XMET_2 ybar_XMET_2-d_XMET_2 ybar_XMET_2+d_XMET_2 n_XMET_2-1;
      ybar_XMET_3 ybar_XMET_3-d_XMET_3 ybar_XMET_3+d_XMET_3 n_XMET_3-1;
      ybar_XRF_1 ybar_XRF_1-d_XRF_1 ybar_XRF_1+d_XRF_1 n_XRF_1-1;
      ybar_XRF_2 ybar_XRF_2-d_XRF_2 ybar_XRF_2+d_XRF_2 n_XRF_2-1;
      ybar_XRF_3 ybar_XRF_3-d_XRF_3 ybar_XRF_3+d_XRF_3 n_XRF_3-1]


