%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Eksperimentel Fysik, Niels Bohr Institutet, Københavns Universitet
% Opgaver tirsdag d. 22. maj og torsdag d. 24. maj (uge 5), 2012
%
% Dette dokument indeholder løsninger til opgaverne stillet i kovarians og
% andet mere avanceret statistik, lavet af Morten Dam Jørgensen (http://mdj.dk). 
% 
% Opgaverne er uploadet i et .pdf-dokument på siden
%   http://mdj.dk/teaching/experimental-physics/
% hvor de forskellige datafiler også ligger pakket ned i .zip filer.
% 
% Pia Jensen, bozack@bozack.dk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all                      %tømmer alle gemte variable
close all                      %luk alle eksterne vinduer (figurer o.lign.)
clc                            %rens kommando-vinduet
datadir = 'data/kovarians/';   %mit data-directory


%% Opgave 1 ---------------------------------------------------------------------
% Disse øvelser viser hvordan lineære korrelationer kan beregnes og visualiseres,
% metoder som er anvendelige på det meste data med flere parametre. 
% 
% Datasættet er hentet fra: http://www.infoplease.com/ipa/A0004379.html

% 1.1
% Start med at indlæse datafilen
lande  = importdata([datadir,'lande.txt'],';');  %indlæs data til cell array
sorted = sortrows(lande.data,1);                 %sorterer efter landeareal
A      = sorted(:,1);                            %lands areal
P      = sorted(:,2);                            %antal indbyggere

% Plot populationen i forhold til landenes størrelse
figure
subplot(1,2,1)
plot(A,P,'.')
xlabel('lands størrelse [miles^2]')
ylabel('antal indbyggere')
title('Lineær sammenhæng')

% 1.2
% Beregn den linære korrelation, manuelt. Er de to variable korrelerede? 
X1    = A - mean(A);
X2    = P - mean(P);
C_man = [ mean( X1 .* X1 )  mean( X1 .* X2 ) ; ...   %manuel kovarians matrice
          mean( X2 .* X1 )  mean( X2 .* X2 ) ];
r_man = mean( X1 .* X2 ) / ( std(A) * std(P) );      %manuel lineær korrelation

% 1.3
% Beregn kovariansmatricen og korrelationsmatricen med MATLAB's indbyggede
% funktioner, og sammenlign resultatet med dine beregninger.
C_aut = cov(A,P);   %MATLAB udregner kovarians matricen
r_aut = corr(A,P);  %MATLAB udregner den lineære korrelation

% 1.4
% Tag logaritmen af begge værdier, og plot resultatet igen. Diskuter forskellen, 
% er korrelationen tydligere nu? 
subplot(1,2,2)
plot(log(A),log(P),'.')
xlabel('log( lands størrelse [miles^2] )')
ylabel('log( antal indbyggere )')
title('Dobbeltlogaritmisk sammenhæng')

% Beregn korrelationen for de logaritmiske værdier. Hvorfor er korrelationen 
% anderledes end i opgave 1.2?
C_log = cov(log(A),log(P));
r_log = corr(log(A),log(P));

% 1.5
% Lav et lineært fit af resultatet i opgave 1.4.
[cfun good] = fit(log(A),log(P),'poly1');
hold on
plot(log(A),cfun(log(A)),'-g')

% Hvad er Chi^2 for fittet, beskriver en lineær relation forholdet?
good.rsquare;   %denne er forholdsvist tæt på 1, så lineær relation er ok


%% Opgave 2 ---------------------------------------------------------------------
% Disse øvelser viser mere med korrelationer.

% 2.1
% Indlæs datasættet "opg2data.txt"
tal = importdata([datadir,'opg2data.txt'],',');  %indlæs data til cell array
x   = tal(:,1);                                  %første kolonne gemt
y   = tal(:,2);                                  %anden kolonne gemt

% Beregn korrelationerne mellem de to variable. Er variablerne korrelerede?
C = cov(x,y);    %kovarians matrice
r = corr(x,y);   %lineær korrelation
% Variablene er ikke korrelerede! (ca. 1.1% korrelation, hvilket ikke er noget)

% 2.2
% Lav et scatter plot af de to variable. Ved visuel inspektion, forklar hvorfor 
% der ikke var en korrelation mellem de to akser i beregningen.
figure
plot(x,y,'.')
xlabel('data fra første kolonne')
ylabel('data fra anden kolonne')
title('Scatter plot af data')

% 2.3
% Datasættet er genereret med følgende udtryk
%   n = 2000; 
%   x = linspace(-1, 1, n);
%   y = - 5 * (x.^2 - 1/2).^2 + unifrnd(-1, 1, [1 n])/3;
%
% Der er tydeligvis en afhængighed mellem x og y aksen. I hvilke tilfælde vil en 
% lineær korrelation være et acceptabelt mål for korrelationen mellem to værdier?

% I tilfælde af overordnet linearitet, uden symmetrier og højere-ordens effekter.


%% Opgave 3 ---------------------------------------------------------------------
% Disse øvelser handler om radial acceleration (med Principal Component Analysis)
%
% Datasæt kommer fra gruppen Allan Hansen og co.

% 3.1
% Indlæs datasættet "ipod_4.txt". De fire variable er hhv. acceleration langs x, 
% y og z akserne samt tid i sekunder.
data = importdata([datadir,'ipod_4.txt']);

% Plot de tre accelerationskomponenter som funktion af tiden.
figure(1)
subplot(3,1,1)
plot(data(:,4), data(:,1))
xlabel('acceleration x-akse (g)')

subplot(3,1,2)
plot(data(:,4), data(:,2))
xlabel('acceleration y-akse (g)')

subplot(3,1,3)
plot(data(:,4), data(:,3))
xlabel('acceleration z-akse (g)')

% PCA Algoritme start
% 3.2
% Find de dominerede komponenter, ved at benytte PCA metoden beskrevet i
% øvelsesvejledningen (afsnit 2).
%
% 3.2.1
% Transponér input matricen sådan at den er M x N, hvor M svarer til antallet af
% parametre og N antallet af målinger.
dat   = data';
[M N] = size(dat);

% Centrér dataen ved først at beregne middelværdien for de fire parametre,
% og fratræk denne de enkelte komponenter. 
% Hint: data_center = data - repmat(middelværdierne, 1, antallet_af_målinger)
m   = mean(dat,2);
dat = dat - repmat(m, 1, N);

% 3.2.2
% Beregn kovariansmatricen af den centrerede data
cor = cov(dat');

% Beregn egenværdierne og egenvektorerne af kovariansmatricen, med den
% indbyggede funktion i MATLAB.
[V E] = eig(cor);

% Gem diagonalen af egenværdimatricen i en vektor (M x 1)
eval = diag(E);

% 3.2.3
% Sorter egenværdierne efter faldende orden, og gem de indexer som 
% sorteringsmetoden retunerer i en vektor.
% Hint: [smidvæk indices] = sort(egenværdier, 'descend')
[drop indices] = sort(eval, 'descend');

% 3.2.4
% Sorter egenvektorerne og egenværdierne med index-vektoren
% Hint: egenvec = egenvec(:,indices); egenværdier = egenværdier(indices)
egvec = V(:, indices);
eval  = eval(indices);

% 3.2.85
% Projekter den oprindelige data langs den nye basis (egenvektor matricen)
% Hint: data_pca = egenvec' * data
data_pca = egvec' * dat;

% Plotning
% 3.3
% Plot hvert komponent af det transformerede datasæt som funktion af tiden. Ved 
% at sammenligne værdierne fra egenvektorne med de fire plots, forklar hvad de 
% enkelte komponenter beskriver.
figure(2)
subplot(4,1,1)
plot(data(:,4), data_pca(1, :))
xlabel('1. Komponent')

subplot(4,1,2)
plot(data(:,4), data_pca(2, :))
xlabel('2. Komponent')

subplot(4,1,3)
plot(data(:,4), data_pca(3, :))
xlabel('3. Komponent')

subplot(4,1,4)
plot(data(:,4), data_pca(4, :))
xlabel('4. Komponent')

% 3.4 (frivillig!)
% Benyt MATLABs indbyggede pca routiner
% 
%   [C, latent, explained] = pcacov(cov(data'))
%   biplot(egenvec(:,1:2), 'scores', signals(1:2,:)', 'varlabels',datalabels)
[C, latent, explained] = pcacov(cov(data));

figure(3)
datalabels = char('x','y','z','t');
biplot(egvec(:,1:2), 'scores', dat(1:2,:)', 'varlabels',datalabels)

% 3.5 (frivillig!)
%
% PCA teknikken afhænger af lineære relationer mellem de forskellige
% variable, der findes en anden metode til ikke-lineære data, der benytte
% det såkaldte "kernel-kneb".
% I zip-filen medfølger en funktion kaldet "kernelpca_tutorial.m", gentag
% opgave 3.5 for resultatet fra funktionen. 
%
% Mere info: http://en.wikipedia.org/wiki/Kernel_principal_component_analysis
data_kpca = kernelpca_tutorial(data', 4);

figure(4)
subplot(4,1,1)
plot(data(:,4), data_kpca(1, :))
xlabel('1. Komponent')

subplot(4,1,2)
plot(data(:,4), data_kpca(2, :))
xlabel('2. Komponent')

subplot(4,1,3)
plot(data(:,4), data_kpca(3, :))
xlabel('3. Komponent')

subplot(4,1,4)
plot(data(:,4), data_kpca(4, :))
xlabel('4. Komponent')