%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Beispiel: Lineare Regression mit Matlab
% 
% This Matlab-program was developed using "Matlab_R2014b" on a MacBook Pro
% running OS X 10.9.5
%
% written by T. Ihn, D-PHYS ETH Zurich, 5 Oct 2014
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Daten aus einem csv-File einlesen
%
data = csvread('daten.csv');

%
% Mit Datenlisten arbeiten
%

% Länge der Datenliste (Zahl der Datenpunkte) ermitteln:
n = length(data);
disp(['n = ' num2str(n)]);

% Auf einzelne Elemente (Datenpunkte) der Datenliste zugreifen:
disp(['3. Datenpunkt:']);
disp(data(3,:))

% Auf Bereiche in der Datenliste zugreifen:
disp('2.-4. Datenpunkt:');
disp(data(2:4,:));

% Nur die x-Werte der Datenpunkte 2 bis 4:
disp('x-Werte der Datenpunkte 2-4:');
disp(data(2:4,1));
% alle x-Werte:
x = data(:,1);

% Nur die y-Werte der Datenpunkte 2 bis 4:
disp('y-Werte der Datenpunkte 2-4:');
disp(data(2:4,2));
% alle y-Werte:
y = data(:,2);

%
% Datenliste graphisch darstellen:
%
figure(1)
set(gca,...
    'FontName','Times New Roman',...
    'FontSize',24,...
    'XTick',0:0.25:1,...
    'Ytick',0:3,...
    'Color',[0.8 0.8 0.8],...
    'LineWidth',2,...
    'TickLength',[0.02 0.0250],...
    'XLim',[0 1],...
    'YLim',[0 3],...
    'Units','inches',...
    'Position',[1 0.8 5.77 3.57],...
    'Layer','Top')
plot(x,y,'bo','MarkerFaceColor','b');
axis([0 9 0 25])
xlabel('x');
ylabel('y');

%
% Statistik der Datenliste:
%

% Mittelwert der x-Werte:
Meanx = mean(x);
disp(['x-Mittelwert = ' num2str(Meanx)]);

% Mittelwert der y-Werte:
Meany = mean(y);
disp(['y-Mittelwert = ' num2str(Meany)]);

% Varianz der x-Werte:
Varx = var(x,1);
disp(['x-Varianz = ' num2str(Varx)]);

% Varianz der y-Werte:
Vary = var(y,1);
disp(['y-Varianz = ' num2str(Vary)]);

% Empirischer Korrelationskoeffizient
rho = corrcoef(x,y);
rho = rho(1,2);
disp(['Linearer Korrelationskoeffizient = ' num2str(rho)]);
disp(' ');
%
% Lineare Regression, Variante I: Schätzwerte und ihre Unsicherheiten (gemäss
% Vorlesung)
%
disp('Lineare Regression, Variante I: gemäss Vorlesung')

A0 = sqrt(Vary/Varx)*rho;
sigmaA = sqrt(C0/Varx/(n-2));
disp(['A0 = ' num2str(A0) ' +/- ' num2str(sigmaA)]);

B0 = Meany;
sigmaB = sqrt(C0/(n-2));
disp(['B0 = ' num2str(B0) ' +/- ' num2str(sigmaB)]);

C0 = Vary*(1-rho^2);
disp(['C0 = ' num2str(C0)]);
disp(' ');

%
% Fitkurve graphisch darstellen
% 
xFit = 0:0.1:9;
yFit = A0*(xFit-Meanx) + B0;
figure(2)
set(gca,...
    'FontName','Times New Roman',...
    'FontSize',24,...
    'XTick',0:0.25:1,...
    'Ytick',0:3,...
    'Color',[0.8 0.8 0.8],...
    'LineWidth',2,...
    'TickLength',[0.02 0.0250],...
    'XLim',[0 1],...
    'YLim',[0 3],...
    'Units','inches',...
    'Position',[1 0.8 5.77 3.57],...
    'Layer','Top')
plot(xFit,yFit,'r-','LineWidth',2);
axis([0 9 0 25])
xlabel('x');
ylabel('y');

%
% Daten mit Fitkurve graphisch darstellen
%
figure(3)
set(gca,...
    'FontName','Times New Roman',...
    'FontSize',24,...
    'XTick',0:0.25:1,...
    'Ytick',0:3,...
    'Color',[0.8 0.8 0.8],...
    'LineWidth',2,...
    'TickLength',[0.02 0.0250],...
    'XLim',[0 1],...
    'YLim',[0 3],...
    'Units','inches',...
    'Position',[1 0.8 5.77 3.57],...
    'Layer','Top')
plot(xFit,yFit,'r-','LineWidth',2);
axis([0 9 0 25])
xlabel('x');
ylabel('y');
hold on
plot(x,y,'bo','MarkerFaceColor','b');
hold off
saveas(gcf,'LinearRegressionMatlabForStudents','pdf');

%
% Variablen löschen
%
clear n

%
% Lineare Regression, Variante II: mit polyfit-Funktion
%
disp('Lineare Regression, Variante II: mit polyfit-Funktion')
[P,S]=polyfit(x-Meanx,y,1);
Rinv = inv(S.R);
cov=(Rinv*Rinv')*(S.normr)^2/S.df;
disp(['A0 = ' num2str(P(1)) ' +/- ' num2str(sqrt(cov(1,1)))]);
disp(['B0 = ' num2str(P(2)) ' +/- ' num2str(sqrt(cov(2,2)))]);
disp(' ');

%
% Lineare Regression, Variante III: mit fitlm-Funktion
%
disp('Lineare Regression, Variante III: mit fitlm-Funktion')
LM = fitlm(x-Meanx,y,'linear')
figure(4)
LM.plot;
title('Output of fitlm');
disp(['A0 = ' num2str(LM.Coefficients.Estimate(2)) ' +/- ' num2str(LM.Coefficients.SE(2))]);
disp(['B0 = ' num2str(LM.Coefficients.Estimate(1)) ' +/- ' num2str(LM.Coefficients.SE(1))]);
