PCA für zweidimensionale Daten

Vorbetrachtung

Bei der Hauptkomponentenanalyse (engl.: PCA - Principal Component Analysis) wird ein mehrdimensionaler Datensatz auf 2 (oder 3) Dimensionen projiziert, um einen visuellen Eindruck von der Varianz der Daten zu ermöglichen. Liegen zweidimensionale Datensätze vor, ist eine Hauptkomponentenanalyse eigentlich nicht sinnvoll, da der Datensatz ja ohnehin komplett visualisiert werden kann. Trotzdem soll hier anhand eines 2D-Datensatzes das Grundprinzip der Hauptkomponentenanalyse demonstriert werden. Wir werden hier also eine zweidimensionale Punktwolke auf eine Linie projizieren. Im Allgemeinen besteht die PCA darin, einen mehrdimensionalen Datensatz so zu im Raum zu rotieren, dass die größtmögliche Varianz nach Projektion auf 2 (oder 3) Dimensionen sichtbar wird.

Parameter for 2D-normalverteilte Daten

Wir simulieren zweidimensionale normalverteilte Daten mit den Mittelwerten \(\mu_x = 1\) und \(\mu_y = 3\). Die Kovarianzmatrix (\(\Sigma\)) muss symmetrisch und positiv definit sein. Letzteres wird mit Hilfe der Funktion is.positive.definite aus dem Package matrixcalc überprüft.

mux = 1 
muy = 3
mu = c(mux, muy)    
Sigma <- matrix(c(4.0,  1.1, 1.1,  2.0), byrow = TRUE, nrow = 2)  
isSymmetric(Sigma)

## [1] TRUE

library(matrixcalc) 
is.positive.definite(Sigma, tol = 1e-8)  

## [1] TRUE

Simulation der Daten

Wir setzen einen seed um die Datenerzeugung wiederholbar zu machen. Zur Erzeugung der Daten wird mvrnorm aus dem Package MASS benutzt.

set.seed(77)
library(MASS) 
data2 = mvrnorm(n = 100, mu = mu, Sigma = Sigma)
data2 = data.frame(x = data2[,1], y = data2[,2])
head(data2, 3)

##            x        y
## 1  1.7432282 4.197688
## 2 -1.9204242 3.889230
## 3 -0.1751774 2.306802

Zentrieren und Skalieren der Daten

Die Hauptkomponentenanalyse beruht im Wesentlichen auf einer Analyse der Varianz von Datensätzen entlang verschiedener Raumrichtungen (siehe unten). Die Daten müssen daher zuvor standardisiert werden, d.h. die Daten müssen so transformiert werden, dass jede Variable den Mittelwert 0 und die Varianz 1 hat. Dies ermöglicht, vereinfacht ausgedrückt, eine “Gleichbehandlung” aller Variablen. Wir definieren zunächst eine Funktion, welche die Standardisierung durchführen kann (dies kann auch mit der in R eingebauten Funktion scale gemacht werden):

zscore <- function(vector) new_vector = (vector - mean(vector))/sd(vector)  

Wir wenden die Transformation nun auf die Daten an und überprüfen das Ergebnis:

points = apply(data2, 2, zscore)
colMeans(points)  

##            x            y 
## 1.932482e-17 5.141720e-17

apply(points, 2, sd) 

## x y 
## 1 1

Die Mittelwerte sind nah genug bei Null, Rundungsfehler können nicht vermieden werden.

Daten plotten

Wir werfen eine Blick auf unsere Daten:

plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")

Projektion von Punkten auf eine Gerade

Bei der Hauptkompontenanalyse wird eine Gerade gesucht, die derart ausgerichtet ist, dass die Projektionen aller Punkte auf diese Gerade eine maximale Varianz haben. Dies ist dann die erste Hauptkomponente. Als “Handwerkszeug” brauchen wir daher zunächst eine Funktion zur Projektion eines Punktes \((x_0,y_0)\) auf eine Gerade \(A \cdot x + B\).
Dazu suchen wir eine zweite Gerade \(a \cdot x + b\) ( blau im Bild ), die senkrecht auf unserer vorgebenen Gerade ( rot im Bild ) steht und durch den Punkt \((x_0, y_0)\) geht. Der Schnittpunkt \((x_s, y_s)\) dieser zweiten Geraden mit der vorgegebenen Gerade ist dann der projizierte Punkt.
Wir können hierbei ausnutzen, dass der Anstieg der zweiten Geraden durch \(a = -1/A\) mit dem Anstieg der vorgebenen Gerade zusammenhängt, da beide senkrecht zueinander stehen. Zum Verständnis dieses Zusammenhanges bedenke man, dass der Anstieg einer Geraden gleich dem Tangens des Steigungswinkels ist, und dass allgemein \(\tan{(\phi+\pi/2)}=-1/tan{(\phi)}\) gilt.

projectedPoint <- function(A, B, x0, y0) {       # line y = A*x + B
  if(A == 0) {            #  1/A not defined, the line is identical with the x-axis
     coords = c(x0, 0)    # ... and the projections are simply the x-values 
  } else {
    a = -1/A              # slope of the projection line (rectangular) 
    b = y0 - a*x0         # intercept of the projection line 
    xs = (b - B)/(A - a)  # x coordinate of projected point               
    ys = A*xs + B         # y coordinate of projected point
    coords = c(xs, ys)
  }
  names(coords) = c("xs", "ys")
  return(coords)
}

Die Funktion projectedPoint soll für den Punkt points[25,] und die Gerade \(y = 0.7 \cdot x\) probiert werden (wir können den Achsenabschnitt \(B = 0\) wählen, da die Daten standardisiert sind).

x0 = points[25,1]
y0 = points[25,2]
projected = projectedPoint(A = 0.7, B = 0, x0 = x0, y0 = y0)
xs = projected[1]
ys = projected[2]
cat(paste0("xs = ", xs, "  ys = ", ys))

## xs = 1.71762212204329  ys = 1.2023354854303

Wir sollten uns mit einem Bild überzeugen, ob wir richtig gerechnet haben:

plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
abline(a = 0, b = 0.7, lwd = 1.8, col = "red")  # we project onto this line 
x0 = points[25,1]  # we project this point 
y0 = points[25,2]
points(x0, y0, col = "green", pch = 18, cex = 1.5)     # the point to project 
points(xs, ys, col = "green", pch = 18, cex = 1.5)     # the projection of this point
segments(x0, y0, xs, ys, col = "darkgreen", lwd = 1.3) # line from (x0,y0) to (xs, ys)

Bemerkung: Beide Achsen des Plots müssen unbedingt gleiche Skalen haben, um das Bild nicht zu verfälschen. Wir wählen daher asp = 1 im plot-Kommando.

Die Projektion scheint geklappt zu haben.
Wir können nun eine Funktion schreiben, die die Projektionen aller Punkte des Datensatzes findet.
Dies kann ganz einfach mit der R-Funktion Vectorize realisiert werden. Wir können damit unsere Funktion projectedPoint so modifizieren, dass sie statt skalarer Werte für \(x_0\) und \(y_0\) Vektoren als Input akzeptiert:

projectedPoints <- Vectorize(projectedPoint, vectorize.args = c("x0", "y0"))

(Die neue Funktion heißt projectedPoints, mit einem zusätzlichen s !)

Lassen wir es auf einen Versuch ankommen:

projected_points = projectedPoints(A = 0.7, B = 0, x0 = points[,1], y0 = points[,2]) # Vektoren als Input
projected_points = t(projected_points)
head(projected_points, 3)

##              xs         ys
## [1,]  0.5945086  0.4161560
## [2,] -0.6297718 -0.4408403
## [3,] -0.6184749 -0.4329324

Wir können nun die Projektionen aller Punkte zeichnen:

plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
abline(a = 0, b = 0.7, lwd = 1.8, col = "red") 
segments(points[,1], points[,2], projected_points[,1], projected_points[,2], col = "darkgreen", lwd = 1.3)

Auch das scheint zu funktionieren.

Berechnung der Varianz der Projektionen

Wie bereits oben erwähnt, wird bei der Hauptkompontenanalyse eine Gerade gesucht, die derart ausgerichtet ist, dass die Projektionen aller Punkte auf diese Gerade eine maximale Varianz haben. Wir brauchen also eine Funktion für die Berechnung der Varianz der projizierten Punkte.
Die empirische Varianz für \(n\) Datenpunkte \(u_i\) mit dem Mittelwert \(\mu\) wird folgendermaßen berechnet:

\[ S^2 = \frac{1}{n-1} \cdot \sum_{i=1}^n \left( u_i - \mu \right)^2 \]

Da die Daten bereits zentriert wurden (\(\mu = 0\)), vereinfacht sich die Formel zu

\[ S^2 = \frac{1}{n-1} \cdot \sum_{i=1}^n u_i^2 \]

Die projizierten Punkte \(u_i\) liegen auf der Geraden \(y = A \cdot x + B\). Die Skizze zeigt, wie ein \(u_i\) mit seinen Koordinaten auf den Achsen zusammenhängt. Wir bezeichen die Varianz mit \(S^2\), da sie das Quadrat der Standardabweichung \(S\) ist.

Die Funktion zur Berechnung der Varianz der projizierten Punkte kann also so aussehen:

getVariance <- function(projected_points) {
  n = nrow(projected_points)
  sum(projected_points[,1]^2 + projected_points[,2]^2)/(n - 1)
}

Wir sollten das gleich einmal mit den oben berechneten projizierten Punkten probieren:

S_squared = getVariance(projected_points)
S_squared

## [1] 1.472078

Berechnung des “Reconstruction Error”

Bevor wir damit beginnen, die optimale Gerade \(y = A \cdot x + B\) zu suchen, wollen wir noch den sogenannten Reconstruction Error (RE) berechnen. Das ist der Fehler, den wir machen, wenn wir einen Punkt vom 2-dimensionalen Raum auf die Gerade projizieren. Der RE ist definiert als die mittlere quadratische Länge der Verbindungslinien zwischen den Originalpunkten \((x_{0,i}, y_{0,i})\) und den projizierten Punkten \((x_{s,i}, y_{s,i})\), denn das ist ja die Strecke, um die wir den Punkt verschieben.
In der folgenden Skizze ist diese Verbindungslinie \(w_i\) rot markiert.

Mit dem ïn der Skizze gezeigten Ausdruck für \(w_i\) wird der RE:

\[ RE = \frac{1}{n} \cdot \sum_{i=1}^n w_i^2 = \frac{1}{n} \cdot \sum_{i=1}^n \left( x_{0,i}^2 + y_{0,i}^2 - x_{s,i}^2 - y_{s,i}^2 \right) \]

Eine Funktion zur Berechnung des RE könnte also so aussehen:

getRE <- function(points, projected_points) { # reconstruction error
  n = nrow(points)
  RE = sum(points[,1]^2 + points[,2]^2 - projected_points[,1]^2 - projected_points[,2]^2)/n 
}

Eine Probe:

RE = getRE(points, projected_points)
RE

## [1] 0.5226427

Summe von Varianz und Reconstruction Error

Ein interessanter Fakt, ist dass die (leicht modifizierte) Summe von RE und Varianz unabhängig von der Ausrichtung der Geraden ist, also insbesondere unabhängig vom Anstieg \(A\). Es gilt für beliebige \(A\):

\[ RE + \left( 1 - \frac{1}{n}\right) \cdot S^2 = C \]

wobei \(C\) eine Konstante ist.

Bemerkung: Wir haben die Varianz mit dem Vorfaktor \(1/(n-1)\) definiert. Mit dem manchmal benutzten, doch nicht korrekten Vorfaktor \(1/n\) würden wir die Beziehung \(RE + S^2 = C\) erhalten.

Wir probieren die oben postulierte Konstanz aus, indem wir den Anstieg der Geraden variieren und jeweils die Summe \(RE + \left( 1 - 1/n\right) \cdot S^2\) berechnen:

df = data.frame(slope = numeric(), RE = numeric(), S2 = numeric(), Sum = numeric())
n = nrow(points)
for (slope in seq(from = 0, to = 3, by = 3/100)) {
  projected_points = projectedPoints(A = slope, B = 0, x0 = points[,1], y0 = points[,2])
  projected_points = t(projected_points)
  S_squared = getVariance(projected_points)  # variance of projected points                
  RE = getRE(points, projected_points)     # reconstruction error      
  C = RE + (1-1/n)*S_squared  
  new_row = data.frame(slope = slope, RE = RE, S2 = S_squared, Sum = C)
  df = rbind(df, new_row)
}
# print(df)

Schauen wir uns das Ergebnis an:

ylim = c( min(df$Sum, df$S2, df$RE ),  max(df$Sum, df$S2, df$RE))
plot(df$slope, df$Sum, col = "red", pch = 20, ylim = ylim,
     xlab = "slope of line", ylab = "RE, S2, and Sum of both", type = "b", lwd = 2, cex = 0.7)
points(df$slope, df$S2, col = "green", pch = 20, type = "b", lwd = 2, cex = 0.7)
points(df$slope, df$RE, col = "blue", pch = 20, type = "b", lwd = 2, cex = 0.7)
legend(x = 0.79, y = 1.2, c("RE", "S2", "Sum"), col = c("blue", "green", "red"), lwd = 2, cex = 1.2)

Wir können ganz leicht beweisen, dass der Ausdruck \(RE + \left( 1 - 1/n\right) \cdot S^2\) bei Drehung der Geraden \(y = A \cdot x + B\) gleich bleiben muss. Wir brauchen nur den Satz des Pythagoras. Schauen wir noch einmal auf die Skizze:

Es ist \(M^2 = u_i^2 + w_i^2\) (oberes Dreieck), andererseits ist \(M^2 = x_{0,i}^2 + y_{0,i}^2\) (unteres Dreieck). Wie können also schreiben

\[ u_i^2 + w_i^2 = x_{0,i}^2 + y_{0,i}^2 \]

Wenn wir über beide Seiten summieren, erhalten wir:

\[ \sum_i u_i^2 + \sum_i w_i^2 = \sum_i \left( x_{0,i}^2 + y_{0,i}^2 \right) \]

Die rechte Seite der Gleichung ist konstant, denn unsere Originalpunkte ändern sich ja nicht bei Drehung der Geraden \(y = A \cdot x + B\)

\[ \sum_i u_i^2 + \sum_i w_i^2 = C \]

Wir benutzen nun den Ausdruck für die Varianz \(S^2 = 1/(n-1) \cdot \sum u_i^2\) und die Definition für den Reconstruction Error \(RE = 1/n \cdot \sum w_i^2\) und erhalten:

\[ (n-1) \cdot S^2 + n \cdot RE = C \]

Wir teilen nur noch durch \(n\) und erhalten den gesuchten Ausdruck (\(C/n\) ist wieder eine Konstante):

\[ RE + \left( 1 - \frac{1}{n}\right) \cdot S^2 = C \]

Gerade mit der größten Varianz der Projektionen

Nun wollen wir die Gerade (d.h. den Anstieg \(A\)) suchen, für welche die Varianz am größten wird: Wir drehen einfach eine durch den Nullpunkt gehende Gerade und registrieren für jede Neigung die empirische Varianz. Wir benutzen dazu die oben definierte Funktion getVariance. Die Variable speed bestimmt, wie schnell die Drehung vonstatten geht. Die Idee für eine solche Simulation stammt von dieser Webseite.

speed = 3
varTable <- data.frame(phi = numeric(0), slope = numeric(0), variance = numeric(0))
for (phi in seq(0, 180, 3)) {   # the step determines the precision of the calculation of the optimum slope
  slope = tan(phi*pi/180)       # tangens needs angle in radians
  projected_points = projectedPoints(A = slope, B = 0, x0 = points[,1], y0 = points[,2]) 
  projected_points = t(projected_points)
  S2 = getVariance(projected_points)
  txt = paste("Angle = ", phi, " Slope = ", round(slope,3), " Variance of projection = ", signif(S2,4))
  
  varRow = data.frame(phi = phi, slope = slope, variance = S2) 
  varTable <- rbind(varTable, varRow) 
  
  plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1)
  abline(h=0, lty = 2, col ="darkgrey")
  abline(v=0, lty = 2, col ="darkgrey")
  abline(a = 0, b = slope, lwd = 1.8, col = "red")
  segments(points[,1], points[,2], projected_points[,1], projected_points[,2], col = "darkgreen", lwd = 1.3)
  mtext(side = 3, txt)
  Sys.sleep(1/speed)
} 

Im Dataframe varTable haben wir die Varianz der Projektionen der Punkte in Abhängigkeit vom Neigungswinkel und von Anstieg regisriert. Wir können jetzt ermitteln, welche Gerade die größte Varianz liefert, indem wir die Varianz über dem Anstieg auftragen:

# plot(varTable$phi, varTable$variance, col = "red", pch = 18, xlab = "angle", ylab = "Variance of projection")  
plot(varTable$slope, varTable$variance, col = "red", pch = 18, xlab = "slope", ylab = "Variance of projection", xlim = c(-3,3)) 
ind = which(varTable$variance == max(varTable$variance))
opt_slope = varTable$slope[ind] 
opt_slope

## [1] 1

abline(v = opt_slope, lty = 3, col = "darkgrey")   

Konstruktion der Hauptkomponenten

Der Anstieg 1 (Winkel 45 Grad) führt zur höchsten Varianz der Projektionen, wie aus der Grafik ersichtlich ist. (Das ist kein Wunder, da ganz am Anfang die Kovarianz-Matrix \(\Sigma\) dementsprechend gewählt wurde¹!) Damit ist die gesuchte optimale Gerade einfach \(y = x\). Wir definieren nun einen Vektor, der in die Richtung dieser Geraden zeigt, indem wir einfach einen Punkt auf der Geraden wählen (die ja durch den Nullpunkt geht). Anschließend normieren wir den Vektor:

normvec <- function(vec) sqrt(sum(vec^2))  # Norm of a vector. Alternative: norm(vec, type = "E")
u = c(1,1)
u = u/normvec(u)
normvec(u)  # Norm 

## [1] 1

u

## [1] 0.7071068 0.7071068

Dieser Vektor ist die 1. Hauptkomponte. Er repräsentiert die Richtung, in der die Punktmenge am meisten variiert. Wollte man wirklich eine eindimensionale Repräsentation der Daten haben, würde man diese Richtung wählen und die Projektionen aller Punkte darauf einzeichnen. (Wie bereits erwähnt, ist das nicht sinnvoll, da die zweidimensionale Darstellung genügend aussagt - wir werden jedoch bei höherdimensionalen Daten das gleiche Prinzip anwenden.)

Eine zweidimensionale Punktmenge hat zwei Haupkomponenten. Die zweite Hauptkomponenete muss senkrecht auf der ersten stehen, d.h. deren Anstieg ist \(a = -1\) (das Produkt beider Anstiege muss \(-1\) sein, da die Geraden senkrecht aufeinander stehen, siehe dazu die Bemerkung oben). Da auch die 2. Hauptkomponente durch den Ursprung geht, ist es wieder leicht, einen Vektor zu konstruieren, der in ihre Richtung zeigt und die Norm \(1\) besitzt, indem man zunächst eine beliebigen Punkt auf dieser Geraden wählt:

v = c(1,-1)
v = v/normvec(v)
normvec(v)  # Norm 

## [1] 1

v

## [1]  0.7071068 -0.7071068

Wir können nun noch die Daten zusammen mit beiden Hauptkomponenten zeichnen:

title = "Dataset and Principal Components"
plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "PC1", ylab = "PC2", asp = 1, main = title, font.main = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
abline(a = 0, b = 1, col = "red", lty = 2, lwd = 1.5)
abline(a = 0, b = -1, col = "red", lty = 2, lwd = 1.5)

Vergleich mit linearer Regression

Wenn man sich die oben gezeigten Plots ansieht, könnte man annehmen, dass das Auffinden der 1. Hauptkomponente im zweidimensionalen Fall identisch mit der Berechnung der Regressionsgeraden ist. Jedoch gibt es einen wesentlichen Unterschied. Wie hatten oben festgestellt, dass das Maximum der Varianz mit dem Minimum des Reconstruction Errors (RE) einhergeht, wobei der RE die mittlere quadratische Entfernung der Originalpunkte von den projizierten Punkten ist. PCA minimimiert also die orthogonalen (senkrechten) Abstände. Bei der linearen Regression werden jedoch die vertikalen Abstände der Originalpunkte von den Regressionspunkten minimiert:

Wir können zum Vergleich die Regressionsgerade für unsere Daten zusammen mit der 1. Hauptkomponente zeichnen. Zunächst führen wir die Regression durch:

lm.out = lm(y ~ x, data = as.data.frame(points))
# str(lm.out)
lm.out$coefficients

## (Intercept)           x 
## 4.16199e-17 5.02426e-01

title = "PC1 and Regression line"
plot(points[,1], points[,2], col = "blue", pch = 20, xlab = "x", ylab = "y", asp = 1, main = title, font.main = 1)
abline(h=0, lty = 2, col ="darkgrey")
abline(v=0, lty = 2, col ="darkgrey")
abline(a = 0, b = 1, col = "red", lty = 1, lwd = 1.5)                                              # 1. Hauptkomponente
abline(a = lm.out$coefficients[1], b = lm.out$coefficients[2], col = "green", lty = 1, lwd = 1.5)  # Regressionsgerade
legend("bottomright", c("1. Hauptkomponente","Regressionsgerade"), col = c("red", "green"), lty = 1, lwd = 2)

Im Übrigen wird bei der Regression vorausgesetzt, dass die \(x\)-Werte der Daten exakt bekannt sind, während bei der PCA angenommen wird, dass sowohl die \(x\)-Werte als auch die \(y\)-Werte verrauscht sind.

Vergleich mit einer in R eingebauten Funktion

Die Funktion prcomp ist eine in R etablierte Funktion zur Hauptkomponentenanalyse. Mit den Optionen center = TRUE und scale = TRUE wird Standardisierung vereinbart, wir können also die ursprünglichen (unskalierten) Daten data2 in die Funktion schicken:

pca <- prcomp(data2, center = TRUE, scale = TRUE) 
summary(pca)

## Importance of components:
##                           PC1    PC2
## Standard deviation     1.2257 0.7054
## Proportion of Variance 0.7512 0.2488
## Cumulative Proportion  0.7512 1.0000

pca$rotation

##          PC1        PC2
## x -0.7071068  0.7071068
## y -0.7071068 -0.7071068

Fällt etwas auf?
Die in pca$rotation gespeicherten Vektoren entsprechen den von uns gefundenen Hauptkomponenten. Die Koordinaten des Vektors PC1 haben lediglich beide ein negatives Vorzeichen, der Vektor repräsentiert jedoch ebenfalls die Gerade \(y = x\), genau wie der von uns gefundene.

Vergleich mit den Eigenvektoren der Kovarianz-Matrix

Zum Schluss sollen die Eigenvektoren der Kovarianz-Matrix der standardisierten Daten berechnet werden. Die Kovarianz kann mit Hilfe der R-Funktion cov berechnet werden, während eigen zur Berechnung der Eigenwerte und Eigenvektoren dient:

covp = cov(points)      # covariance matrix
ev = eigen(covp)  
evalues = ev$values     # eigenvalues
evectors = ev$vectors   # eigenvectors
evectors

##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

Die Eigenvektoren der Kovarianz-Matrix entsprechen ebenfalls den zuvor berechneten Vektoren.

Zusammenfassung

Wir haben:

1. die erste Hauptkomponente eines Datensatzes experimentell gefunden, indem wir einen (normierten) Vektor ermittelt haben, der so positioniert ist, dass die Projektionen aller Datenpunkte auf diesen Vektor eine maximale Varianz haben. Die zweite Hauptkomponente wurde so gelegt, dass sie senkrecht auf der ersten steht (und ebenfalls normiert ist),
   
2. die Funktion prcomp benutzt, um die Haupkomponenten zu bestimmen,
   
3. die Eigenvektoren der Kovarianz-Matrix berechnet.

Alle 3 Vektorenpaare sind identisch (bis auf die Vorzeichen der mit prcomp berechneten 1. Hauptkomponente). Im Theorie-Teil werden wir sehen, wie das alles zusammenhängt.

Anhang

Länge der Projektion

Wir wollen einen Ausdruck für die Länge der Projektion eines Punktes \(\vec{x}\) auf einen Vektor \(\vec{v}\) finden, d.h. wir suchen die Größe \(a\) in der folgenden Grafik:

Das Skalarprodukt von \(\vec{x}\) und \(\vec{v}\) ist

\[\begin{align*} \vec{x} \cdot \vec{v} &= |\vec{x}| \cdot |\vec{v}| \cdot \cos{\varphi} \\ &= |\vec{x}| \cdot |\vec{v}| \cdot \frac{a}{|\vec{x}|} \\ &= |\vec{v}| \cdot a \end{align*}\]

Wir bekommen also:

\[ a = \frac{\vec{x} \cdot \vec{v}}{|\vec{v}| } \] Wenn der Vektor \(\vec{v}\) die Norm 1 hat, vereinfacht sich dies zu

\[ a = \vec{x} \cdot \vec{v} \]

Der Skalar \(a\) ist die Länge der Projektion des Vektors \(\vec{x}\) auf die Gerade \(\vec{v}\). Diese Formel wird in der Funktion proVar benutzt.

Die Norm eines Vektors

Die Norm eines \(n\)-dimensionalen Vektors ist seine Länge im \(n\)-dimensionalen Raum. Sei \(\vec{u}\) ein Vektor mit den Komponenten \(u_1, u_2, \dots u_n\):

\[ \begin{align} \vec{u} &= \begin{bmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{bmatrix} \end{align} \]

Dann ist seine Norm

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2 + \ldots + u_n^2} \]

Im Dreidimensionalen wird das

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2 + u_3^2} \]

… und im Zweidimensionalen:

\[ \lVert \vec{u} \rVert = \sqrt{u_1^2 + u_2^2} \]

Dies ist der Satz des Pythagoras.

Die Norm kann auch mit der R-internen Funktion norm berechnet werden.

Die Funktion prcomp in R

Die R-Funktion prcomp führt eine Hauptkomponentenanalyse von mehrdimensionalen Daten durch. Die offizielle Dokumentation ist hier.

Eigenwerte und Eigenvektoren

Die Eigenvektoren einer Matrix \(\boldsymbol{A}\) sind die nicht-trivialen (\(\vec{x} \ne 0\)) Lösungen der Gleichung

\[ \boldsymbol{A} \cdot \vec{x} = \lambda \cdot \vec{x} \]

Diese Gleichung ist im Allgemeinen nur für bestimmte Werte von \(\lambda\) (die Eigenwerte) lösbar. Die Eigenwerte ergeben sich aus der Gleichung

\[ \det(\boldsymbol{A} - \lambda \cdot \boldsymbol{I}) = 0 \] wobei \(\boldsymbol{I}\) die Einheitsmatrix ist. Mehr dazu bei Wikipedia

Note that all code snippets are designed for comprehensibility, not for efficiency!

uwe.menzel@matstat.org

---

1. Das kann gemacht werden, indem man ausgehend von den Eigenvektoren die Kovarianzmatrix konstruiert, also den umgekehrten Weg beschreitet.↩︎