% Regner ut sannsynlighetsmaksimeringsestimat for koeffisientene b0 og b1,
% samt standardavvikene s1, s2 og s3 ved å maksimere likelihoodfunksjonen
% numerisk. Bruker kun XMET-data.
% INPUT: .txt-fil med XMET-målinger i sjette kolonne og mineraliserings-
% klasse i syvende kolonne.
% OUTPUT: Vektoren b = [b0 b1 s1 s2 s3] som inneholder estimatene, og
% et plott av y = b0 + b1*x for x = 1, 2, 3, med errorbars som indikerer
% estimert standardavvik for hver verdi av x.
clc, clear all, close all

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

find1 = find(data(:,7) == 1); % Finn radene med klasse 1-observasjoner
find2 = find(data(:,7) == 2); % Tilsvarende for klasse 2...
find3 = find(data(:,7) == 3); % ...og klasse 3.

n1 = length(find1); % Antall klasse 1-observasjoner
n2 = length(find2);
n3 = length(find3);

x1 = ones(size(find1)); % Forklaringsvariabelen, mineraliseringsklasse
x2 = 2.*ones(size(find2));
x3 = 3.*ones(size(find3));

y1 = data(find1,6); % Responsen, oksidinnhold målt med XMET
y2 = data(find2,6);
y3 = data(find3,6);

% Definerer -ln(L) hvor L er likelihoodfunksjonen til residualene epsilon_i = y_i - (b0 + b1*x_i), i = 1,...,n1,...,n1+n2,...,n1+n2+n3
negloglike = @(b) -(-(1/2).*(n1.*log(2.*pi.*b(3).^2)+n2.*log(2.*pi.*b(4).^2)+n3.*log(2.*pi.*b(5).^2)) -(1/2).*(sum(((y1-b(1)-b(2).*x1)/b(3)).^2) + sum(((y2-b(1)-b(2).*x2)/b(4)).^2) + sum(((y3-b(1)-b(2).*x3)/b(5)).^2)));

x0 = [-1, 1, 1, 2, 3]; % Oppgir et startpunkt for optimeringsfunksjonen, får samme resultat med et annet startpunkt

[b,fval] = fminunc(negloglike,x0) % Kaller fminunc (function minimization, unconstrained) fra Optimization Toolbox
% b er vektoren med sannsynlighetsmaksimeringsestimatene, mens fval er den minimale funksjonsverdien av -ln(L)

% Plott med errorbars
figure
hold on
% handles til plott-elementene lagres i variablene til venstre for likhetstegnene
p1 = plot(x1,y1,'g+','LineWidth',0.6);
p2 = plot(x2,y2,'g+','LineWidth',0.6);
p3 = plot(x3,y3,'g+','LineWidth',0.6);
p4 = errorbar([1 2 3],[b(1)+b(2) b(1)+2*b(2) b(1)+3*b(2)],[b(3) b(4) b(5)],'ks--');

set(get(get(p2,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); % Utelat disse (p2, p3) fra tegnforklaringen
set(get(get(p3,'Annotation'),'LegendInformation'),'IconDisplayStyle','off'); % siden det er tre serier med like symboler
legend('XMET-data','{y = \beta_0 + \beta_1 x}','Location','NorthWest') % Lag tegnforklaringen med kun to linjer

xlabel('Klasse, x')
ylabel('Oksidinnhold, y')
% Bruker [] for å sette sammen flere strenger til én streng. Bruker num2str() for å få med estimatene i tittelen.
title(['{\beta_0 = }' num2str(b(1)) '{, \beta_1 = }' num2str(b(2)) '{, \sigma_1 = }' num2str(b(3)) '{, \sigma_2 = }' num2str(b(4)) '{, \sigma_3 = }' num2str(b(5))],'FontSize',12)
box on
hold off
