opgaver:Uge5

Antal indbyggere som en funktion af landeareal.

Datafilen kan indlæses med koden

lande = importdata('lande.txt',';');

For nemmere at kunne arbejde med data sorteres de nu i rækkefølge efter landeareal, og population og areal skrives ind i hvert sit array:

sorted = sortrows(lande.data,1);

A = sorted(:,1);

P = sorted(:,2);

Plottet laves med

figure

plot(A,P,'.')

xlabel('lands størrelse [miles^2]')

ylabel('antal indbyggere')

title('Lineær sammenhæng')

Dette plot kan ses her til højre.

Beregn først kovarians matricen (ligning 1 i multivariatstatistik.pdf dokumentet).

Covariansen mellem to vektorer er givet ved $\Sigma_{i,j} = \text{cov}(X_i,X_j) = E[(X_i-\mu_i)(X_j-\mu_j)]$, hvor $E[\cdot]$ er forventningsværdien - simpelthen gennemsnittet. Dermed kan man manuelt udregne covariansmatricen ved først at udregne

X1 = A - mean(A);

X2 = P - mean(P);

og kan man simpelthen manuelt indsætte i kovarians matricen,

C_man = [ mean( X1 .* X1 ) mean( X1 .* X2 ) ; ...

mean( X2 .* X1 ) mean( X2 .* X2 ) ];

Tilsvarende kan man manuelt finde den lineære korrelation med

r_man = mean( X1 .* X2 ) / ( std(A) * std(P) );

Værdien af r_man er ca. 0.45, så de variable er korrelerede, men ikke super-godt.

Læs om funktionerne cov og corr i MATLAB hjælpen.

Man bruger simpelthen de to funktioner cov og corr på A og P vektorerne,

C_aut = cov(A,P);

r_aut = corr(A,P);

Antal indbyggere som en funktion af landeareal, logaritmisk.

Man plotter simpelthen

figure

plot(log(A),log(P),'.')

xlabel('log( lands størrelse [miles^2] )')

ylabel('log( antal indbyggere )')

title('Dobbeltlogaritmisk sammenhæng')

Nu er det meget nemmere at se på data, da den tydeligvis er mere eksponentielt fordelt. Ud fra figuren forventer man også at korrelationen mellem de logaritmiske vektorer er mere lineær, da den faktisk ser lineær ud nu!

C_log = cov(log(A),log(P));

r_log = corr(log(A),log(P));

Det viser sig at r_log er omkring 0.86, hvilket præcis er som forventet - korrelationen er bedre for de logaritmiske.

Antal indbyggere som en funktion af landeareal, logaritmisk, sammen med et lineært fit.

Fittet laves bare med

[cfun good] = fit(log(A),log(P),'poly1');

og kan plottes oven i den forrige figur med

hold on

plot(log(A),cfun(log(A)),'-g')

For at finde ud af hvor godt fittet er, kan man udregne $\chi^2$ for det. Man kan også bruge den $R^2$ værdi som fit funktionen giver. Denne er givet i

good.rsquare;

Denne viser sig at være meget tæt på 1, hvilket fortæller at fittet er ret godt.

Altså er en lineær relation god til at beskrive forholdet mellem logaritmen af antallet af indbyggere i et lang og logaritmen af arealet af landet.

Først indlæses datasættet med

tal = importdata([datadir,'opg2data.txt'],',');

x = tal(:,1);

y = tal(:,2);

Korrelationerne mellem de to variable findes med de indbyggede MATLAB funktioner,

C = cov(x,y);

r = corr(x,y);

Det viser sig at der er cirka 1.1 % korrelation mellem variablene, og derfor er de slet ikke korrelerede!

Disse variable er ikke lineært korrelerede.

Plottet laves simpelthen med koden

figure

plot(x,y,'.')

xlabel('data fra første kolonne')

ylabel('data fra anden kolonne')

title('Scatter plot af data')

Dette plot er vist her til højre.

Det er tydeligt at se at der er en eller anden sammenhæng mellem de to akser, men den er absolut ikke lineær. Den lineære korrelation tager kun (logisk nok) linearitet med i beregningen, så alle højere ordens korrelationer bliver ikke fundet. For at kunne se højere ordens korrelationer skal man bruge mere avancerede statistiske værktøjer, som der bliver taget hul på i næste opgave.

Lineær korrelation er et acceptabelt mål for korrelation i tilfælde hvor der er overordnet linearitet, uden symmetrier og højere-ordens effekter.

De tre accelerations-retninger.

Data indlæses med

data = importdata([datadir,'ipod_4.txt']);

og man kan plotte de tre retnings-akser i samme figur med subplot

figure

subplot(3,1,1)

plot(data(:,4), data(:,1))

title('acceleration x-akse (g)')

subplot(3,1,2)

plot(data(:,4), data(:,2))

title('acceleration y-akse (g)')

subplot(3,1,3)

plot(data(:,4), data(:,3))

title('acceleration z-akse (g)')

Disse tre plots kan ses i figuren til højre.

For at centrere dataen kan det være praktisk at bruge en kode a la

data_center = data - repmat(middelværdierne, 1, antallet_af_målinger)

Prøv at læse lidt om repmat i MATLAB hjælpen.

Transponeringen gøres med den simple ' operator, og size kan bruges til at få M og N,

dat = data';

[M N] = size(dat);

Centreringen af dataen kan nu gøres som vist i hintet ovenfor, ved først at finde gennemsnittene i et nyt array,

m = mean(dat,2);

og derefter bruge repmat til at lave et array der kan trækkes fra de oprindelige data (leg lidt med den, hvis du ikke forstår hvordan den virker).

dat = dat - repmat(m, 1, N);

Man kan finde egenværdier og egenvektorer for en matrice i MATLAB ved at bruge eig funktionen. Man kan trække diagonalen af en matrice ud til en vektor med diag funktionen.

Kovariansmatricen udregnes som sædvanligt, med

cor = cov(dat');

Egenværdierne og egenvektorerne for denne matrice findes ved hjælp af den indbyggede funktion i MATLAB, der giver

[V E] = eig(cor);

Her er V en matrice der indeholder egenvektorerne, mens E er en matrice der indeholder egenværdierne i diagonalen. Denne diagonal gemmes med

eval = diag(E);

Læs lidt om sort funktionen i MATLAB hjælpen. Man kan få den til at spytte indekser ud ved hjælp kode på formen

[smidvæk indices] = sort(egenværdier, 'descend')

Dette gøres simpelthen med koden

[drop indices] = sort(eval, 'descend');

hvor drop ikke skal bruges til noget videre.

Brug de indices som du gemte i del 3 til at bytte rundt på dine egenværdier og egenvektorer. Dette kan gøres med kode som

egenvec = egenvec(:,indices); egenværdier = egenværdier(indices)

Egenvektorerne i den ønskede rækkefølge findes med

egvec = V(:, indices);

mens egenværdierne i samme rækkefølge findes med

eval = eval(indices);

For at projicere skal man bare gange egenvektor matricen (transponeret, for at dimensionerne passer) sammen med dataen, som

data_pca = egvec' * dat;

De fire komponenter efter PCA analysen.

Igen kan man plotte flere plots i samme figur ved at bruge subplot.

figure

subplot(4,1,1)

plot(data(:,4), data_pca(1, :))

title('1. Komponent')

subplot(4,1,2)

plot(data(:,4), data_pca(2, :))

title('2. Komponent')

subplot(4,1,3)

plot(data(:,4), data_pca(3, :))

title('3. Komponent')

subplot(4,1,4)

plot(data(:,4), data_pca(4, :))

title('4. Komponent')

For at kunne vurdere hvad de forskellige komponenter er, er det nødvendigt at se på egenvektorerne fundet med PCA metoden. I MATLAB findes det at

egvec =

0.0019 -0.7804 -0.6252 0.0092

-0.0098 -0.6253 0.7803 -0.0096

0.0001 -0.0012 -0.0132 -0.9999

1.0000 -0.0046 0.0088 0.0000

I første omgang kan man altså se at første komponent efter analysen næsten kun indeholder 4.-komponenten af dataen (tid), hvilket også giver fin mening med den monotont stigende kurve set i figuren. På samme måde kan man se at den fjerde komponent efter analysen næsten kun indeholder (den negative) 3.-komponent af dataen ($z$-accelerationen). Denne kurve ligner meget godt den støj der også blev set på den oprindelige $z$-acceleration komponent i dataen.

Derimod er anden og tredje komponenterne efter analysen givet som et miks imellem $x$- og $y$-accelerationerne fra den oprindelige data.