% Utfører simpel lineær regresjon på snødata
% Dette scriptet utfører simpel lineær regresjon ved hjelp av matriser,
% Se W.M.M.Y kap. 12.3 for detaljer.
% INPUT: .txt-filer med observasjoner fra Rybekken og Sylsjødammen med snødybde,
% snøtetthet, SWE og DOY i henholdsvis andre, tredje, fjerde og femte kolonne
% OUTPUT: Tre figurer som hver er delt inn i fire underfigurer:
% 1) SWE ~ dybde
% 2) SWE ~ tetthet
% 3) tetthet ~ dybde
% 4) tetthet ~ DOY
% Den første figuren viser regresjonslinje og spredningsplott av datapunkter,
% den andre viser Q-Q plott av residualene etter regresjonen, og den tredje
% plotter residualene direkte.

clc, clear all, close all

% Last inn data
sy = load('sylsjodammen.txt');
ry = load('rybekken.txt');
kombinert = [ry; sy];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1) Regresjon av SWE mot dybde
x = kombinert(:,2);
y = kombinert(:,4);

% Lagre antall observasjoner
n = length(x);

% Oprett design-matrise
X = [ones(n,1) x];

% Løs likningssystem X'X(beta) = X'y
beta = (X'*X)\X'*y;
beta0 = beta(1);
beta1 = beta(2);

% Regner ut feilen og R^2
y_hat = beta0 + beta1.*x;
e = y - y_hat; % Residualer
y_bar = mean(y);
R2 = 1 - sum(e.^2)/sum((y-y_bar).^2);

% Lag figur med plott av regresjonslinje og datapunkter
figure(1)

subplot(2,2,1), hold on
plot(x,y,'bx','LineWidth',1.2);
plot(sort(x),beta0 + beta1.*sort(x),'r-','LineWidth',2);
title(['{SWE \sim dybde, R^2 = }' num2str(R2)],'FontWeight','bold')
xlabel('{Dybde (cm)}')
ylabel('{SWE (cm)}')
legend('Data','Regresjonslinje','Location','NorthEast')
set(gcf,'Color','White')
box on
hold off

% Lager figur med Q-Q plott av residualer
figure(2)

subplot(2,2,1)
qqplot(e)
title('{Residualer, SWE \sim dybde}','FontWeight','Bold')
box on

% Lag figur med residualplott
figure(3)

subplot(2,2,1), hold on
plot(e,'b+')
plot(zeros(length(e)),'k--')
xlim([0 length(e)])
xlabel('{Obs. nr.}')
ylabel('{SWE (cm)}')
legend('Residualer','Location','NorthEast')
title('{Residualer, SWE \sim dybde}','FontWeight','Bold')
box on
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2) Regresjon av SWE mot tetthet
x = kombinert(:,3);

% Ny design-matrise
X = [ones(n,1) x];

% Løs likningssystem X'X(beta) = X'y
beta = (X'*X)\X'*y;
beta0 = beta(1);
beta1 = beta(2);

% Regner ut feilen og R^2
y_hat = beta0 + beta1.*x;
e = y - y_hat;
y_bar = mean(y);
R2 = 1 - sum(e.^2)/sum((y-y_bar).^2);

% Lag resten av figuren
figure(1)

subplot(2,2,2), hold on
plot(x,y,'bx','LineWidth',1.2);
plot(sort(x),beta0 + beta1.*sort(x),'r-','LineWidth',2);
title(['{SWE \sim tetthet, R^2 = }' num2str(R2)],'FontWeight','bold')
xlabel('{Tetthet (g/cm^3)}')
ylabel('{SWE (cm)}')
legend('Data','Regresjonslinje','Location','NorthEast')
set(gcf,'Color','White')
box on
hold off

% Lager figur med Q-Q plott av residualer
figure(2)

subplot(2,2,2)
qqplot(e)
title('{Residualer, SWE \sim tetthet}','FontWeight','Bold')
box on

% Lag figur med residualplott
figure(3)

subplot(2,2,2), hold on
plot(e,'b+')
plot(zeros(length(e)),'k--')
xlim([0 length(e)])
xlabel('{Obs. nr.}')
ylabel('{SWE (cm)}')
legend('Residualer','Location','NorthEast')
title('{Residualer, SWE \sim tetthet}','FontWeight','Bold')
box on
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3) Regresjon av tetthet mot dybde
x = kombinert(:,2);
y = kombinert(:,3);

% Ny design-matrise
X = [ones(n,1) x];

% Løs likningssystem X'X(beta) = X'y
beta = (X'*X)\X'*y;
beta0 = beta(1);
beta1 = beta(2);

% Regner ut feilen og R^2
y_hat = beta0 + beta1.*x;
e = y - y_hat;
y_bar = mean(y);
R2 = 1 - sum(e.^2)/sum((y-y_bar).^2);

% Lag resten av figuren
figure(1)

subplot(2,2,3), hold on
plot(x,y,'bx','LineWidth',1.2);
plot(sort(x),beta0 + beta1.*sort(x),'r-','LineWidth',2);
title(['{tetthet \sim dybde, R^2 = }' num2str(R2)],'FontWeight','bold')
ylabel('{Tetthet (g/cm^3)}')
xlabel('{Dybde (cm)}')
legend('Data','Regresjonslinje','Location','NorthEast')
set(gcf,'Color','White')
box on
hold off

% Lager figur med Q-Q plott av residualer
figure(2)

subplot(2,2,3)
qqplot(e)
title('{Residualer, tetthet \sim dybde}','FontWeight','Bold')
box on

% Lag figur med residualplott
figure(3)

subplot(2,2,3), hold on
plot(e,'b+')
plot(zeros(length(e)),'k--')
xlim([0 length(e)])
ylabel('{Tetthet (g/cm^3)}')
xlabel('{Obs. nr.}')
legend('Residualer','Location','NorthEast')
title('{Residualer, tetthet \sim dybde}','FontWeight','Bold')
box on
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  4) Regresjon av tetthet mot DOY
x = kombinert(:,5);

% Ny design-matrise
X = [ones(n,1) x];

% Løs likningssystem X'X(beta) = X'y
beta = (X'*X)\X'*y;
beta0 = beta(1);
beta1 = beta(2);

% Regner ut feilen og R^2
y_hat = beta0 + beta1.*x;
e = y - y_hat;
y_bar = mean(y);
R2 = 1 - sum(e.^2)/sum((y-y_bar).^2);

% Lag resten av figuren
figure(1)

subplot(2,2,4), hold on
plot(x,y,'bx','LineWidth',1.2);
plot(sort(x),beta0 + beta1.*sort(x),'r-','LineWidth',2);
title(['{tetthet \sim DOY, R^2 = }' num2str(R2)],'FontWeight','bold')
ylabel('{Tetthet (g/cm^3)}')
xlabel('{DOY}')
legend('Data','Regresjonslinje','Location','NorthEast')
set(gcf,'Color','White')
box on
hold off

% Lager figur med Q-Q plott av residualer
figure(2)

subplot(2,2,4)
qqplot(e)
title('{Residualer, tetthet \sim DOY}','FontWeight','Bold')
box on

% Lag figur med residualplott
figure(3)

subplot(2,2,4), hold on
plot(e,'b+')
plot(zeros(length(e)),'k--')
xlim([0 length(e)])
ylabel('{Tetthet (g/cm^3)}')
xlabel('{Obs. nr.}')
legend('Residualer','Location','NorthEast')
title('{Residualer, tetthet \sim DOY}','FontWeight','Bold')
box on
hold off
